MagickCore 7.1.1-43
Convert, Edit, Or Compose Bitmap Images
Loading...
Searching...
No Matches
distort.c
1/*
2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3% %
4% %
5% %
6% DDDD IIIII SSSSS TTTTT OOO RRRR TTTTT %
7% D D I SS T O O R R T %
8% D D I SSS T O O RRRR T %
9% D D I SS T O O R R T %
10% DDDD IIIII SSSSS T OOO R R T %
11% %
12% %
13% MagickCore Image Distortion Methods %
14% %
15% Software Design %
16% Cristy %
17% Anthony Thyssen %
18% June 2007 %
19% %
20% %
21% Copyright @ 1999 ImageMagick Studio LLC, a non-profit organization %
22% dedicated to making software imaging solutions freely available. %
23% %
24% You may not use this file except in compliance with the License. You may %
25% obtain a copy of the License at %
26% %
27% https://imagemagick.org/script/license.php %
28% %
29% Unless required by applicable law or agreed to in writing, software %
30% distributed under the License is distributed on an "AS IS" BASIS, %
31% WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
32% See the License for the specific language governing permissions and %
33% limitations under the License. %
34% %
35%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36%
37%
38*/
39␌
40/*
41 Include declarations.
42*/
43#include "MagickCore/studio.h"
44#include "MagickCore/artifact.h"
45#include "MagickCore/cache.h"
46#include "MagickCore/cache-view.h"
47#include "MagickCore/channel.h"
48#include "MagickCore/colorspace-private.h"
49#include "MagickCore/composite-private.h"
50#include "MagickCore/distort.h"
51#include "MagickCore/exception.h"
52#include "MagickCore/exception-private.h"
53#include "MagickCore/gem.h"
54#include "MagickCore/image.h"
55#include "MagickCore/linked-list.h"
56#include "MagickCore/list.h"
57#include "MagickCore/matrix.h"
58#include "MagickCore/matrix-private.h"
59#include "MagickCore/memory_.h"
60#include "MagickCore/monitor-private.h"
61#include "MagickCore/option.h"
62#include "MagickCore/pixel.h"
63#include "MagickCore/pixel-accessor.h"
64#include "MagickCore/resample.h"
65#include "MagickCore/resample-private.h"
66#include "MagickCore/registry.h"
67#include "MagickCore/resource_.h"
68#include "MagickCore/semaphore.h"
69#include "MagickCore/shear.h"
70#include "MagickCore/string_.h"
71#include "MagickCore/string-private.h"
72#include "MagickCore/thread-private.h"
73#include "MagickCore/token.h"
74#include "MagickCore/transform.h"
75␌
76/*
77 Numerous internal routines for image distortions.
78*/
79static inline void AffineArgsToCoefficients(double *affine)
80{
81 /* map external sx,ry,rx,sy,tx,ty to internal c0,c2,c4,c1,c3,c5 */
82 double tmp[4]; /* note indexes 0 and 5 remain unchanged */
83 tmp[0]=affine[1]; tmp[1]=affine[2]; tmp[2]=affine[3]; tmp[3]=affine[4];
84 affine[3]=tmp[0]; affine[1]=tmp[1]; affine[4]=tmp[2]; affine[2]=tmp[3];
85}
86
87static inline void CoefficientsToAffineArgs(double *coeff)
88{
89 /* map internal c0,c1,c2,c3,c4,c5 to external sx,ry,rx,sy,tx,ty */
90 double tmp[4]; /* note indexes 0 and 5 remain unchanged */
91 tmp[0]=coeff[3]; tmp[1]=coeff[1]; tmp[2]=coeff[4]; tmp[3]=coeff[2];
92 coeff[1]=tmp[0]; coeff[2]=tmp[1]; coeff[3]=tmp[2]; coeff[4]=tmp[3];
93}
94static void InvertAffineCoefficients(const double *coeff,double *inverse)
95{
96 /* From "Digital Image Warping" by George Wolberg, page 50 */
97 double determinant;
98
99 determinant=PerceptibleReciprocal(coeff[0]*coeff[4]-coeff[1]*coeff[3]);
100 inverse[0]=determinant*coeff[4];
101 inverse[1]=determinant*(-coeff[1]);
102 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[2]*coeff[4]);
103 inverse[3]=determinant*(-coeff[3]);
104 inverse[4]=determinant*coeff[0];
105 inverse[5]=determinant*(coeff[2]*coeff[3]-coeff[0]*coeff[5]);
106}
107
108static void InvertPerspectiveCoefficients(const double *coeff,
109 double *inverse)
110{
111 /* From "Digital Image Warping" by George Wolberg, page 53 */
112 double determinant;
113
114 determinant=PerceptibleReciprocal(coeff[0]*coeff[4]-coeff[3]*coeff[1]);
115 inverse[0]=determinant*(coeff[4]-coeff[7]*coeff[5]);
116 inverse[1]=determinant*(coeff[7]*coeff[2]-coeff[1]);
117 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[4]*coeff[2]);
118 inverse[3]=determinant*(coeff[6]*coeff[5]-coeff[3]);
119 inverse[4]=determinant*(coeff[0]-coeff[6]*coeff[2]);
120 inverse[5]=determinant*(coeff[3]*coeff[2]-coeff[0]*coeff[5]);
121 inverse[6]=determinant*(coeff[3]*coeff[7]-coeff[6]*coeff[4]);
122 inverse[7]=determinant*(coeff[6]*coeff[1]-coeff[0]*coeff[7]);
123}
124
125/*
126 * Polynomial Term Defining Functions
127 *
128 * Order must either be an integer, or 1.5 to produce
129 * the 2 number_valuesal polynomial function...
130 * affine 1 (3) u = c0 + c1*x + c2*y
131 * bilinear 1.5 (4) u = '' + c3*x*y
132 * quadratic 2 (6) u = '' + c4*x*x + c5*y*y
133 * cubic 3 (10) u = '' + c6*x^3 + c7*x*x*y + c8*x*y*y + c9*y^3
134 * quartic 4 (15) u = '' + c10*x^4 + ... + c14*y^4
135 * quintic 5 (21) u = '' + c15*x^5 + ... + c20*y^5
136 * number in parenthesis minimum number of points needed.
137 * Anything beyond quintic, has not been implemented until
138 * a more automated way of determining terms is found.
139
140 * Note the slight re-ordering of the terms for a quadratic polynomial
141 * which is to allow the use of a bi-linear (order=1.5) polynomial.
142 * All the later polynomials are ordered simply from x^N to y^N
143 */
144static size_t poly_number_terms(double order)
145{
146 /* Return the number of terms for a 2d polynomial */
147 if ( order < 1 || order > 5 ||
148 ( order != floor(order) && (order-1.5) > MagickEpsilon) )
149 return 0; /* invalid polynomial order */
150 return((size_t) floor((order+1)*(order+2)/2));
151}
152
153static double poly_basis_fn(ssize_t n, double x, double y)
154{
155 /* Return the result for this polynomial term */
156 switch(n) {
157 case 0: return( 1.0 ); /* constant */
158 case 1: return( x );
159 case 2: return( y ); /* affine order = 1 terms = 3 */
160 case 3: return( x*y ); /* bilinear order = 1.5 terms = 4 */
161 case 4: return( x*x );
162 case 5: return( y*y ); /* quadratic order = 2 terms = 6 */
163 case 6: return( x*x*x );
164 case 7: return( x*x*y );
165 case 8: return( x*y*y );
166 case 9: return( y*y*y ); /* cubic order = 3 terms = 10 */
167 case 10: return( x*x*x*x );
168 case 11: return( x*x*x*y );
169 case 12: return( x*x*y*y );
170 case 13: return( x*y*y*y );
171 case 14: return( y*y*y*y ); /* quartic order = 4 terms = 15 */
172 case 15: return( x*x*x*x*x );
173 case 16: return( x*x*x*x*y );
174 case 17: return( x*x*x*y*y );
175 case 18: return( x*x*y*y*y );
176 case 19: return( x*y*y*y*y );
177 case 20: return( y*y*y*y*y ); /* quintic order = 5 terms = 21 */
178 }
179 return( 0 ); /* should never happen */
180}
181static const char *poly_basis_str(ssize_t n)
182{
183 /* return the result for this polynomial term */
184 switch(n) {
185 case 0: return(""); /* constant */
186 case 1: return("*ii");
187 case 2: return("*jj"); /* affine order = 1 terms = 3 */
188 case 3: return("*ii*jj"); /* bilinear order = 1.5 terms = 4 */
189 case 4: return("*ii*ii");
190 case 5: return("*jj*jj"); /* quadratic order = 2 terms = 6 */
191 case 6: return("*ii*ii*ii");
192 case 7: return("*ii*ii*jj");
193 case 8: return("*ii*jj*jj");
194 case 9: return("*jj*jj*jj"); /* cubic order = 3 terms = 10 */
195 case 10: return("*ii*ii*ii*ii");
196 case 11: return("*ii*ii*ii*jj");
197 case 12: return("*ii*ii*jj*jj");
198 case 13: return("*ii*jj*jj*jj");
199 case 14: return("*jj*jj*jj*jj"); /* quartic order = 4 terms = 15 */
200 case 15: return("*ii*ii*ii*ii*ii");
201 case 16: return("*ii*ii*ii*ii*jj");
202 case 17: return("*ii*ii*ii*jj*jj");
203 case 18: return("*ii*ii*jj*jj*jj");
204 case 19: return("*ii*jj*jj*jj*jj");
205 case 20: return("*jj*jj*jj*jj*jj"); /* quintic order = 5 terms = 21 */
206 }
207 return( "UNKNOWN" ); /* should never happen */
208}
209static double poly_basis_dx(ssize_t n, double x, double y)
210{
211 /* polynomial term for x derivative */
212 switch(n) {
213 case 0: return( 0.0 ); /* constant */
214 case 1: return( 1.0 );
215 case 2: return( 0.0 ); /* affine order = 1 terms = 3 */
216 case 3: return( y ); /* bilinear order = 1.5 terms = 4 */
217 case 4: return( x );
218 case 5: return( 0.0 ); /* quadratic order = 2 terms = 6 */
219 case 6: return( x*x );
220 case 7: return( x*y );
221 case 8: return( y*y );
222 case 9: return( 0.0 ); /* cubic order = 3 terms = 10 */
223 case 10: return( x*x*x );
224 case 11: return( x*x*y );
225 case 12: return( x*y*y );
226 case 13: return( y*y*y );
227 case 14: return( 0.0 ); /* quartic order = 4 terms = 15 */
228 case 15: return( x*x*x*x );
229 case 16: return( x*x*x*y );
230 case 17: return( x*x*y*y );
231 case 18: return( x*y*y*y );
232 case 19: return( y*y*y*y );
233 case 20: return( 0.0 ); /* quintic order = 5 terms = 21 */
234 }
235 return( 0.0 ); /* should never happen */
236}
237static double poly_basis_dy(ssize_t n, double x, double y)
238{
239 /* polynomial term for y derivative */
240 switch(n) {
241 case 0: return( 0.0 ); /* constant */
242 case 1: return( 0.0 );
243 case 2: return( 1.0 ); /* affine order = 1 terms = 3 */
244 case 3: return( x ); /* bilinear order = 1.5 terms = 4 */
245 case 4: return( 0.0 );
246 case 5: return( y ); /* quadratic order = 2 terms = 6 */
247 default: return( poly_basis_dx(n-1,x,y) ); /* weird but true */
248 }
249 /* NOTE: the only reason that last is not true for 'quadratic'
250 is due to the re-arrangement of terms to allow for 'bilinear'
251 */
252}
253␌
254/*
255%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
256% %
257% %
258% %
259% A f f i n e T r a n s f o r m I m a g e %
260% %
261% %
262% %
263%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
264%
265% AffineTransformImage() transforms an image as dictated by the affine matrix.
266% It allocates the memory necessary for the new Image structure and returns
267% a pointer to the new image.
268%
269% The format of the AffineTransformImage method is:
270%
271% Image *AffineTransformImage(const Image *image,
272% AffineMatrix *affine_matrix,ExceptionInfo *exception)
273%
274% A description of each parameter follows:
275%
276% o image: the image.
277%
278% o affine_matrix: the affine matrix.
279%
280% o exception: return any errors or warnings in this structure.
281%
282*/
283MagickExport Image *AffineTransformImage(const Image *image,
284 const AffineMatrix *affine_matrix,ExceptionInfo *exception)
285{
286 double
287 distort[6];
288
289 Image
290 *deskew_image;
291
292 /*
293 Affine transform image.
294 */
295 assert(image->signature == MagickCoreSignature);
296 assert(affine_matrix != (AffineMatrix *) NULL);
297 assert(exception != (ExceptionInfo *) NULL);
298 assert(exception->signature == MagickCoreSignature);
299 if (IsEventLogging() != MagickFalse)
300 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
301 distort[0]=affine_matrix->sx;
302 distort[1]=affine_matrix->rx;
303 distort[2]=affine_matrix->ry;
304 distort[3]=affine_matrix->sy;
305 distort[4]=affine_matrix->tx;
306 distort[5]=affine_matrix->ty;
307 deskew_image=DistortImage(image,AffineProjectionDistortion,6,distort,
308 MagickTrue,exception);
309 return(deskew_image);
310}
311␌
312/*
313%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
314% %
315% %
316% %
317+ G e n e r a t e C o e f f i c i e n t s %
318% %
319% %
320% %
321%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
322%
323% GenerateCoefficients() takes user provided input arguments and generates
324% the coefficients, needed to apply the specific distortion for either
325% distorting images (generally using control points) or generating a color
326% gradient from sparsely separated color points.
327%
328% The format of the GenerateCoefficients() method is:
329%
330% Image *GenerateCoefficients(const Image *image,DistortMethod method,
331% const size_t number_arguments,const double *arguments,
332% size_t number_values, ExceptionInfo *exception)
333%
334% A description of each parameter follows:
335%
336% o image: the image to be distorted.
337%
338% o method: the method of image distortion/ sparse gradient
339%
340% o number_arguments: the number of arguments given.
341%
342% o arguments: the arguments for this distortion method.
343%
344% o number_values: the style and format of given control points, (caller type)
345% 0: 2 dimensional mapping of control points (Distort)
346% Format: u,v,x,y where u,v is the 'source' of the
347% the color to be plotted, for DistortImage()
348% N: Interpolation of control points with N values (usually r,g,b)
349% Format: x,y,r,g,b mapping x,y to color values r,g,b
350% IN future, variable number of values may be given (1 to N)
351%
352% o exception: return any errors or warnings in this structure
353%
354% Note that the returned array of double values must be freed by the
355% calling method using RelinquishMagickMemory(). This however may change in
356% the future to require a more 'method' specific method.
357%
358% Because of this, this method should not be classed as stable or used
359% outside other MagickCore library methods.
360*/
361
362static inline double MagickRound(double x)
363{
364 /*
365 Round the fraction to nearest integer.
366 */
367 if ((x-floor(x)) < (ceil(x)-x))
368 return(floor(x));
369 return(ceil(x));
370}
371
372static double *GenerateCoefficients(const Image *image,
373 DistortMethod *method,const size_t number_arguments,const double *arguments,
374 size_t number_values,ExceptionInfo *exception)
375{
376 double
377 *coeff;
378
379 size_t
380 i;
381
382 size_t
383 number_coefficients, /* number of coefficients to return (array size) */
384 cp_size, /* number floating point numbers per control point */
385 cp_x,cp_y, /* the x,y indexes for control point */
386 cp_values; /* index of values for this control point */
387 /* number_values Number of values given per control point */
388
389 if ( number_values == 0 ) {
390 /* Image distortion using control points (or other distortion)
391 That is generate a mapping so that x,y->u,v given u,v,x,y
392 */
393 number_values = 2; /* special case: two values of u,v */
394 cp_values = 0; /* the values i,j are BEFORE the destination CP x,y */
395 cp_x = 2; /* location of x,y in input control values */
396 cp_y = 3;
397 /* NOTE: cp_values, also used for later 'reverse map distort' tests */
398 }
399 else {
400 cp_x = 0; /* location of x,y in input control values */
401 cp_y = 1;
402 cp_values = 2; /* and the other values are after x,y */
403 /* Typically in this case the values are R,G,B color values */
404 }
405 cp_size = number_values+2; /* each CP definition involves this many numbers */
406
407 /* If not enough control point pairs are found for specific distortions
408 fall back to Affine distortion (allowing 0 to 3 point pairs)
409 */
410 if ( number_arguments < 4*cp_size &&
411 ( *method == BilinearForwardDistortion
412 || *method == BilinearReverseDistortion
413 || *method == PerspectiveDistortion
414 ) )
415 *method = AffineDistortion;
416
417 number_coefficients=0;
418 switch (*method) {
419 case AffineDistortion:
420 case RigidAffineDistortion:
421 /* also BarycentricColorInterpolate: */
422 number_coefficients=3*number_values;
423 break;
424 case PolynomialDistortion:
425 /* number of coefficients depend on the given polynomial 'order' */
426 i = poly_number_terms(arguments[0]);
427 number_coefficients = 2 + i*number_values;
428 if ( i == 0 ) {
429 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
430 "InvalidArgument","%s : '%s'","Polynomial",
431 "Invalid order, should be integer 1 to 5, or 1.5");
432 return((double *) NULL);
433 }
434 if ((number_arguments < (1+i*cp_size)) ||
435 (((number_arguments-1) % cp_size) != 0)) {
436 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
437 "InvalidArgument", "%s : 'require at least %.20g CPs'",
438 "Polynomial", (double) i);
439 return((double *) NULL);
440 }
441 break;
442 case BilinearReverseDistortion:
443 number_coefficients=4*number_values;
444 break;
445 /*
446 The rest are constants as they are only used for image distorts
447 */
448 case BilinearForwardDistortion:
449 number_coefficients=10; /* 2*4 coeff plus 2 constants */
450 cp_x = 0; /* Reverse src/dest coords for forward mapping */
451 cp_y = 1;
452 cp_values = 2;
453 break;
454#if 0
455 case QuadrilateralDistortion:
456 number_coefficients=19; /* BilinearForward + BilinearReverse */
457#endif
458 break;
459 case ShepardsDistortion:
460 number_coefficients=1; /* The power factor to use */
461 break;
462 case ArcDistortion:
463 number_coefficients=5;
464 break;
465 case ScaleRotateTranslateDistortion:
466 case AffineProjectionDistortion:
467 case Plane2CylinderDistortion:
468 case Cylinder2PlaneDistortion:
469 number_coefficients=6;
470 break;
471 case PolarDistortion:
472 case DePolarDistortion:
473 number_coefficients=8;
474 break;
475 case PerspectiveDistortion:
476 case PerspectiveProjectionDistortion:
477 number_coefficients=9;
478 break;
479 case BarrelDistortion:
480 case BarrelInverseDistortion:
481 number_coefficients=10;
482 break;
483 default:
484 perror("unknown method given"); /* just fail assertion */
485 }
486
487 /* allocate the array of coefficients needed */
488 coeff=(double *) AcquireQuantumMemory(number_coefficients,sizeof(*coeff));
489 if (coeff == (double *) NULL)
490 {
491 (void) ThrowMagickException(exception,GetMagickModule(),
492 ResourceLimitError,"MemoryAllocationFailed","%s",
493 "GenerateCoefficients");
494 return((double *) NULL);
495 }
496
497 /* zero out coefficients array */
498 for (i=0; i < number_coefficients; i++)
499 coeff[i] = 0.0;
500
501 switch (*method)
502 {
503 case AffineDistortion:
504 {
505 /* Affine Distortion
506 v = c0*x + c1*y + c2
507 for each 'value' given
508
509 Input Arguments are sets of control points...
510 For Distort Images u,v, x,y ...
511 For Sparse Gradients x,y, r,g,b ...
512 */
513 if ( number_arguments%cp_size != 0 ||
514 number_arguments < cp_size ) {
515 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
516 "InvalidArgument", "%s : 'require at least %.20g CPs'",
517 "Affine", 1.0);
518 coeff=(double *) RelinquishMagickMemory(coeff);
519 return((double *) NULL);
520 }
521 /* handle special cases of not enough arguments */
522 if ( number_arguments == cp_size ) {
523 /* Only 1 CP Set Given */
524 if ( cp_values == 0 ) {
525 /* image distortion - translate the image */
526 coeff[0] = 1.0;
527 coeff[2] = arguments[0] - arguments[2];
528 coeff[4] = 1.0;
529 coeff[5] = arguments[1] - arguments[3];
530 }
531 else {
532 /* sparse gradient - use the values directly */
533 for (i=0; i<number_values; i++)
534 coeff[i*3+2] = arguments[cp_values+i];
535 }
536 }
537 else {
538 /* 2 or more points (usually 3) given.
539 Solve a least squares simultaneous equation for coefficients.
540 */
541 double
542 **matrix,
543 **vectors,
544 terms[3];
545
546 MagickBooleanType
547 status;
548
549 /* create matrix, and a fake vectors matrix */
550 matrix=AcquireMagickMatrix(3UL,3UL);
551 vectors=(double **) AcquireQuantumMemory(number_values,
552 sizeof(*vectors));
553 if (matrix == (double **) NULL || vectors == (double **) NULL)
554 {
555 matrix = RelinquishMagickMatrix(matrix, 3UL);
556 vectors = (double **) RelinquishMagickMemory(vectors);
557 coeff = (double *) RelinquishMagickMemory(coeff);
558 (void) ThrowMagickException(exception,GetMagickModule(),
559 ResourceLimitError,"MemoryAllocationFailed",
560 "%s", "DistortCoefficients");
561 return((double *) NULL);
562 }
563 /* fake a number_values x3 vectors matrix from coefficients array */
564 for (i=0; i < number_values; i++)
565 vectors[i] = &(coeff[i*3]);
566 /* Add given control point pairs for least squares solving */
567 for (i=0; i < number_arguments; i+=cp_size) {
568 terms[0] = arguments[i+cp_x]; /* x */
569 terms[1] = arguments[i+cp_y]; /* y */
570 terms[2] = 1; /* 1 */
571 LeastSquaresAddTerms(matrix,vectors,terms,
572 &(arguments[i+cp_values]),3UL,number_values);
573 }
574 if ( number_arguments == 2*cp_size ) {
575 /* Only two pairs were given, but we need 3 to solve the affine.
576 Fake extra coordinates by rotating p1 around p0 by 90 degrees.
577 x2 = x0 - (y1-y0) y2 = y0 + (x1-x0)
578 */
579 terms[0] = arguments[cp_x]
580 - ( arguments[cp_size+cp_y] - arguments[cp_y] ); /* x2 */
581 terms[1] = arguments[cp_y] +
582 + ( arguments[cp_size+cp_x] - arguments[cp_x] ); /* y2 */
583 terms[2] = 1; /* 1 */
584 if ( cp_values == 0 ) {
585 /* Image Distortion - rotate the u,v coordinates too */
586 double
587 uv2[2];
588 uv2[0] = arguments[0] - arguments[5] + arguments[1]; /* u2 */
589 uv2[1] = arguments[1] + arguments[4] - arguments[0]; /* v2 */
590 LeastSquaresAddTerms(matrix,vectors,terms,uv2,3UL,2UL);
591 }
592 else {
593 /* Sparse Gradient - use values of p0 for linear gradient */
594 LeastSquaresAddTerms(matrix,vectors,terms,
595 &(arguments[cp_values]),3UL,number_values);
596 }
597 }
598 /* Solve for LeastSquares Coefficients */
599 status=GaussJordanElimination(matrix,vectors,3UL,number_values);
600 matrix = RelinquishMagickMatrix(matrix, 3UL);
601 vectors = (double **) RelinquishMagickMemory(vectors);
602 if ( status == MagickFalse ) {
603 coeff = (double *) RelinquishMagickMemory(coeff);
604 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
605 "InvalidArgument","%s : 'Unsolvable Matrix'",
606 CommandOptionToMnemonic(MagickDistortOptions, *method) );
607 return((double *) NULL);
608 }
609 }
610 return(coeff);
611 }
612 case RigidAffineDistortion:
613 {
614 double
615 inverse[6],
616 **matrix,
617 terms[5],
618 *vectors[1];
619
620 MagickBooleanType
621 status;
622
623 /*
624 Rigid affine (also known as a Euclidean transform), restricts affine
625 coefficients to 4 (S, R, Tx, Ty) with Sy=Sx and Ry = -Rx so that one has
626 only scale, rotation and translation. No skew.
627 */
628 if (((number_arguments % cp_size) != 0) || (number_arguments < cp_size))
629 {
630 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
631 "InvalidArgument", "%s : 'require at least %.20g CPs'",
632 CommandOptionToMnemonic(MagickDistortOptions,*method),2.0);
633 coeff=(double *) RelinquishMagickMemory(coeff);
634 return((double *) NULL);
635 }
636 /*
637 Rigid affine requires a 4x4 least-squares matrix (zeroed).
638 */
639 matrix=AcquireMagickMatrix(4UL,4UL);
640 if (matrix == (double **) NULL)
641 {
642 coeff=(double *) RelinquishMagickMemory(coeff);
643 (void) ThrowMagickException(exception,GetMagickModule(),
644 ResourceLimitError,"MemoryAllocationFailed","%s",
645 CommandOptionToMnemonic(MagickDistortOptions,*method));
646 return((double *) NULL);
647 }
648 /*
649 Add control points for least squares solving.
650 */
651 vectors[0]=(&(coeff[0]));
652 for (i=0; i < number_arguments; i+=4)
653 {
654 terms[0]=arguments[i+0];
655 terms[1]=(-arguments[i+1]);
656 terms[2]=1.0;
657 terms[3]=0.0;
658 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+2]),4UL,1UL);
659 terms[0]=arguments[i+1];
660 terms[1]=arguments[i+0];
661 terms[2]=0.0;
662 terms[3]=1.0;
663 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+3]),4UL,1UL);
664 }
665 /*
666 Solve for least-squares coefficients.
667 */
668 status=GaussJordanElimination(matrix,vectors,4UL,1UL);
669 matrix=RelinquishMagickMatrix(matrix,4UL);
670 if (status == MagickFalse)
671 {
672 coeff=(double *) RelinquishMagickMemory(coeff);
673 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
674 "InvalidArgument","%s : 'Unsolvable Matrix'",
675 CommandOptionToMnemonic(MagickDistortOptions,*method));
676 return((double *) NULL);
677 }
678 /*
679 Convert (S, R, Tx, Ty) to an affine projection.
680 */
681 inverse[0]=coeff[0];
682 inverse[1]=coeff[1];
683 inverse[2]=(-coeff[1]);
684 inverse[3]=coeff[0];
685 inverse[4]=coeff[2];
686 inverse[5]=coeff[3];
687 AffineArgsToCoefficients(inverse);
688 InvertAffineCoefficients(inverse,coeff);
689 *method=AffineDistortion;
690 return(coeff);
691 }
692 case AffineProjectionDistortion:
693 {
694 /*
695 Arguments: Affine Matrix (forward mapping)
696 Arguments sx, rx, ry, sy, tx, ty
697 Where u = sx*x + ry*y + tx
698 v = rx*x + sy*y + ty
699
700 Returns coefficients (in there inverse form) ordered as...
701 sx ry tx rx sy ty
702
703 AffineProjection Distortion Notes...
704 + Will only work with a 2 number_values for Image Distortion
705 + Can not be used for generating a sparse gradient (interpolation)
706 */
707 double inverse[8];
708 if (number_arguments != 6) {
709 coeff = (double *) RelinquishMagickMemory(coeff);
710 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
711 "InvalidArgument","%s : 'Needs 6 coeff values'",
712 CommandOptionToMnemonic(MagickDistortOptions, *method) );
713 return((double *) NULL);
714 }
715 /* FUTURE: trap test for sx*sy-rx*ry == 0 (determinant = 0, no inverse) */
716 for(i=0; i<6UL; i++ )
717 inverse[i] = arguments[i];
718 AffineArgsToCoefficients(inverse); /* map into coefficients */
719 InvertAffineCoefficients(inverse, coeff); /* invert */
720 *method = AffineDistortion;
721
722 return(coeff);
723 }
724 case ScaleRotateTranslateDistortion:
725 {
726 /* Scale, Rotate and Translate Distortion
727 An alternative Affine Distortion
728 Argument options, by number of arguments given:
729 7: x,y, sx,sy, a, nx,ny
730 6: x,y, s, a, nx,ny
731 5: x,y, sx,sy, a
732 4: x,y, s, a
733 3: x,y, a
734 2: s, a
735 1: a
736 Where actions are (in order of application)
737 x,y 'center' of transforms (default = image center)
738 sx,sy scale image by this amount (default = 1)
739 a angle of rotation (argument required)
740 nx,ny move 'center' here (default = x,y or no movement)
741 And convert to affine mapping coefficients
742
743 ScaleRotateTranslate Distortion Notes...
744 + Does not use a set of CPs in any normal way
745 + Will only work with a 2 number_valuesal Image Distortion
746 + Cannot be used for generating a sparse gradient (interpolation)
747 */
748 double
749 cosine, sine,
750 x,y,sx,sy,a,nx,ny;
751
752 /* set default center, and default scale */
753 x = nx = (double)(image->columns)/2.0 + (double)image->page.x;
754 y = ny = (double)(image->rows)/2.0 + (double)image->page.y;
755 sx = sy = 1.0;
756 switch ( number_arguments ) {
757 case 0:
758 coeff = (double *) RelinquishMagickMemory(coeff);
759 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
760 "InvalidArgument","%s : 'Needs at least 1 argument'",
761 CommandOptionToMnemonic(MagickDistortOptions, *method) );
762 return((double *) NULL);
763 case 1:
764 a = arguments[0];
765 break;
766 case 2:
767 sx = sy = arguments[0];
768 a = arguments[1];
769 break;
770 default:
771 x = nx = arguments[0];
772 y = ny = arguments[1];
773 switch ( number_arguments ) {
774 case 3:
775 a = arguments[2];
776 break;
777 case 4:
778 sx = sy = arguments[2];
779 a = arguments[3];
780 break;
781 case 5:
782 sx = arguments[2];
783 sy = arguments[3];
784 a = arguments[4];
785 break;
786 case 6:
787 sx = sy = arguments[2];
788 a = arguments[3];
789 nx = arguments[4];
790 ny = arguments[5];
791 break;
792 case 7:
793 sx = arguments[2];
794 sy = arguments[3];
795 a = arguments[4];
796 nx = arguments[5];
797 ny = arguments[6];
798 break;
799 default:
800 coeff = (double *) RelinquishMagickMemory(coeff);
801 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
802 "InvalidArgument","%s : 'Too Many Arguments (7 or less)'",
803 CommandOptionToMnemonic(MagickDistortOptions, *method) );
804 return((double *) NULL);
805 }
806 break;
807 }
808 /* Trap if sx or sy == 0 -- image is scaled out of existence! */
809 if ( fabs(sx) < MagickEpsilon || fabs(sy) < MagickEpsilon ) {
810 coeff = (double *) RelinquishMagickMemory(coeff);
811 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
812 "InvalidArgument","%s : 'Zero Scale Given'",
813 CommandOptionToMnemonic(MagickDistortOptions, *method) );
814 return((double *) NULL);
815 }
816 /* Save the given arguments as an affine distortion */
817 a=DegreesToRadians(a); cosine=cos(a); sine=sin(a);
818
819 *method = AffineDistortion;
820 coeff[0]=cosine/sx;
821 coeff[1]=sine/sx;
822 coeff[2]=x-nx*coeff[0]-ny*coeff[1];
823 coeff[3]=(-sine)/sy;
824 coeff[4]=cosine/sy;
825 coeff[5]=y-nx*coeff[3]-ny*coeff[4];
826 return(coeff);
827 }
828 case PerspectiveDistortion:
829 { /*
830 Perspective Distortion (a ratio of affine distortions)
831
832 p(x,y) c0*x + c1*y + c2
833 u = ------ = ------------------
834 r(x,y) c6*x + c7*y + 1
835
836 q(x,y) c3*x + c4*y + c5
837 v = ------ = ------------------
838 r(x,y) c6*x + c7*y + 1
839
840 c8 = Sign of 'r', or the denominator affine, for the actual image.
841 This determines what part of the distorted image is 'ground'
842 side of the horizon, the other part is 'sky' or invalid.
843 Valid values are +1.0 or -1.0 only.
844
845 Input Arguments are sets of control points...
846 For Distort Images u,v, x,y ...
847 For Sparse Gradients x,y, r,g,b ...
848
849 Perspective Distortion Notes...
850 + Can be thought of as ratio of 3 affine transformations
851 + Not separable: r() or c6 and c7 are used by both equations
852 + All 8 coefficients must be determined simultaneously
853 + Will only work with a 2 number_valuesal Image Distortion
854 + Can not be used for generating a sparse gradient (interpolation)
855 + It is not linear, but is simple to generate an inverse
856 + All lines within an image remain lines.
857 + but distances between points may vary.
858 */
859 double
860 **matrix,
861 *vectors[1],
862 terms[8];
863
864 size_t
865 cp_u = cp_values,
866 cp_v = cp_values+1;
867
868 MagickBooleanType
869 status;
870
871 if ( number_arguments%cp_size != 0 ||
872 number_arguments < cp_size*4 ) {
873 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
874 "InvalidArgument", "%s : 'require at least %.20g CPs'",
875 CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0);
876 coeff=(double *) RelinquishMagickMemory(coeff);
877 return((double *) NULL);
878 }
879 /* fake 1x8 vectors matrix directly using the coefficients array */
880 vectors[0] = &(coeff[0]);
881 /* 8x8 least-squares matrix (zeroed) */
882 matrix = AcquireMagickMatrix(8UL,8UL);
883 if (matrix == (double **) NULL) {
884 coeff=(double *) RelinquishMagickMemory(coeff);
885 (void) ThrowMagickException(exception,GetMagickModule(),
886 ResourceLimitError,"MemoryAllocationFailed",
887 "%s", "DistortCoefficients");
888 return((double *) NULL);
889 }
890 /* Add control points for least squares solving */
891 for (i=0; i < number_arguments; i+=4) {
892 terms[0]=arguments[i+cp_x]; /* c0*x */
893 terms[1]=arguments[i+cp_y]; /* c1*y */
894 terms[2]=1.0; /* c2*1 */
895 terms[3]=0.0;
896 terms[4]=0.0;
897 terms[5]=0.0;
898 terms[6]=-terms[0]*arguments[i+cp_u]; /* 1/(c6*x) */
899 terms[7]=-terms[1]*arguments[i+cp_u]; /* 1/(c7*y) */
900 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_u]),
901 8UL,1UL);
902
903 terms[0]=0.0;
904 terms[1]=0.0;
905 terms[2]=0.0;
906 terms[3]=arguments[i+cp_x]; /* c3*x */
907 terms[4]=arguments[i+cp_y]; /* c4*y */
908 terms[5]=1.0; /* c5*1 */
909 terms[6]=-terms[3]*arguments[i+cp_v]; /* 1/(c6*x) */
910 terms[7]=-terms[4]*arguments[i+cp_v]; /* 1/(c7*y) */
911 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_v]),
912 8UL,1UL);
913 }
914 /* Solve for LeastSquares Coefficients */
915 status=GaussJordanElimination(matrix,vectors,8UL,1UL);
916 matrix = RelinquishMagickMatrix(matrix, 8UL);
917 if ( status == MagickFalse ) {
918 coeff = (double *) RelinquishMagickMemory(coeff);
919 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
920 "InvalidArgument","%s : 'Unsolvable Matrix'",
921 CommandOptionToMnemonic(MagickDistortOptions, *method) );
922 return((double *) NULL);
923 }
924 /*
925 Calculate 9'th coefficient! The ground-sky determination.
926 What is sign of the 'ground' in r() denominator affine function?
927 Just use any valid image coordinate (first control point) in
928 destination for determination of what part of view is 'ground'.
929 */
930 coeff[8] = coeff[6]*arguments[cp_x]
931 + coeff[7]*arguments[cp_y] + 1.0;
932 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
933
934 return(coeff);
935 }
936 case PerspectiveProjectionDistortion:
937 {
938 /*
939 Arguments: Perspective Coefficients (forward mapping)
940 */
941 if (number_arguments != 8) {
942 coeff = (double *) RelinquishMagickMemory(coeff);
943 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
944 "InvalidArgument", "%s : 'Needs 8 coefficient values'",
945 CommandOptionToMnemonic(MagickDistortOptions, *method));
946 return((double *) NULL);
947 }
948 /* FUTURE: trap test c0*c4-c3*c1 == 0 (determinate = 0, no inverse) */
949 InvertPerspectiveCoefficients(arguments, coeff);
950 /*
951 Calculate 9'th coefficient! The ground-sky determination.
952 What is sign of the 'ground' in r() denominator affine function?
953 Just use any valid image coordinate in destination for determination.
954 For a forward mapped perspective the images 0,0 coord will map to
955 c2,c5 in the distorted image, so set the sign of denominator of that.
956 */
957 coeff[8] = coeff[6]*arguments[2]
958 + coeff[7]*arguments[5] + 1.0;
959 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
960 *method = PerspectiveDistortion;
961
962 return(coeff);
963 }
964 case BilinearForwardDistortion:
965 case BilinearReverseDistortion:
966 {
967 /* Bilinear Distortion (Forward mapping)
968 v = c0*x + c1*y + c2*x*y + c3;
969 for each 'value' given
970
971 This is actually a simple polynomial Distortion! The difference
972 however is when we need to reverse the above equation to generate a
973 BilinearForwardDistortion (see below).
974
975 Input Arguments are sets of control points...
976 For Distort Images u,v, x,y ...
977 For Sparse Gradients x,y, r,g,b ...
978
979 */
980 double
981 **matrix,
982 **vectors,
983 terms[4];
984
985 MagickBooleanType
986 status;
987
988 /* check the number of arguments */
989 if ( number_arguments%cp_size != 0 ||
990 number_arguments < cp_size*4 ) {
991 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
992 "InvalidArgument", "%s : 'require at least %.20g CPs'",
993 CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0);
994 coeff=(double *) RelinquishMagickMemory(coeff);
995 return((double *) NULL);
996 }
997 /* create matrix, and a fake vectors matrix */
998 matrix=AcquireMagickMatrix(4UL,4UL);
999 vectors=(double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
1000 if (matrix == (double **) NULL || vectors == (double **) NULL)
1001 {
1002 matrix = RelinquishMagickMatrix(matrix, 4UL);
1003 vectors = (double **) RelinquishMagickMemory(vectors);
1004 coeff = (double *) RelinquishMagickMemory(coeff);
1005 (void) ThrowMagickException(exception,GetMagickModule(),
1006 ResourceLimitError,"MemoryAllocationFailed",
1007 "%s", "DistortCoefficients");
1008 return((double *) NULL);
1009 }
1010 /* fake a number_values x4 vectors matrix from coefficients array */
1011 for (i=0; i < number_values; i++)
1012 vectors[i] = &(coeff[i*4]);
1013 /* Add given control point pairs for least squares solving */
1014 for (i=0; i < number_arguments; i+=cp_size) {
1015 terms[0] = arguments[i+cp_x]; /* x */
1016 terms[1] = arguments[i+cp_y]; /* y */
1017 terms[2] = terms[0]*terms[1]; /* x*y */
1018 terms[3] = 1; /* 1 */
1019 LeastSquaresAddTerms(matrix,vectors,terms,
1020 &(arguments[i+cp_values]),4UL,number_values);
1021 }
1022 /* Solve for LeastSquares Coefficients */
1023 status=GaussJordanElimination(matrix,vectors,4UL,number_values);
1024 matrix = RelinquishMagickMatrix(matrix, 4UL);
1025 vectors = (double **) RelinquishMagickMemory(vectors);
1026 if ( status == MagickFalse ) {
1027 coeff = (double *) RelinquishMagickMemory(coeff);
1028 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1029 "InvalidArgument","%s : 'Unsolvable Matrix'",
1030 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1031 return((double *) NULL);
1032 }
1033 if ( *method == BilinearForwardDistortion ) {
1034 /* Bilinear Forward Mapped Distortion
1035
1036 The above least-squares solved for coefficients but in the forward
1037 direction, due to changes to indexing constants.
1038
1039 i = c0*x + c1*y + c2*x*y + c3;
1040 j = c4*x + c5*y + c6*x*y + c7;
1041
1042 where i,j are in the destination image, NOT the source.
1043
1044 Reverse Pixel mapping however needs to use reverse of these
1045 functions. It required a full page of algebra to work out the
1046 reversed mapping formula, but resolves down to the following...
1047
1048 c8 = c0*c5-c1*c4;
1049 c9 = 2*(c2*c5-c1*c6); // '2*a' in the quadratic formula
1050
1051 i = i - c3; j = j - c7;
1052 b = c6*i - c2*j + c8; // So that a*y^2 + b*y + c == 0
1053 c = c4*i - c0*j; // y = ( -b +- sqrt(bb - 4ac) ) / (2*a)
1054
1055 r = b*b - c9*(c+c);
1056 if ( c9 != 0 )
1057 y = ( -b + sqrt(r) ) / c9;
1058 else
1059 y = -c/b;
1060
1061 x = ( i - c1*y) / ( c1 - c2*y );
1062
1063 NB: if 'r' is negative there is no solution!
1064 NB: the sign of the sqrt() should be negative if image becomes
1065 flipped or flopped, or crosses over itself.
1066 NB: technically coefficient c5 is not needed, anymore,
1067 but kept for completeness.
1068
1069 See Anthony Thyssen <A.Thyssen@griffith.edu.au>
1070 or Fred Weinhaus <fmw@alink.net> for more details.
1071
1072 */
1073 coeff[8] = coeff[0]*coeff[5] - coeff[1]*coeff[4];
1074 coeff[9] = 2*(coeff[2]*coeff[5] - coeff[1]*coeff[6]);
1075 }
1076 return(coeff);
1077 }
1078#if 0
1079 case QuadrilateralDistortion:
1080 {
1081 /* Map a Quadrilateral to a unit square using BilinearReverse
1082 Then map that unit square back to the final Quadrilateral
1083 using BilinearForward.
1084
1085 Input Arguments are sets of control points...
1086 For Distort Images u,v, x,y ...
1087 For Sparse Gradients x,y, r,g,b ...
1088
1089 */
1090 /* UNDER CONSTRUCTION */
1091 return(coeff);
1092 }
1093#endif
1094
1095 case PolynomialDistortion:
1096 {
1097 /* Polynomial Distortion
1098
1099 First two coefficients are used to hole global polynomial information
1100 c0 = Order of the polynomial being created
1101 c1 = number_of_terms in one polynomial equation
1102
1103 Rest of the coefficients map to the equations....
1104 v = c0 + c1*x + c2*y + c3*x*y + c4*x^2 + c5*y^2 + c6*x^3 + ...
1105 for each 'value' (number_values of them) given.
1106 As such total coefficients = 2 + number_terms * number_values
1107
1108 Input Arguments are sets of control points...
1109 For Distort Images order [u,v, x,y] ...
1110 For Sparse Gradients order [x,y, r,g,b] ...
1111
1112 Polynomial Distortion Notes...
1113 + UNDER DEVELOPMENT -- Do not expect this to remain as is.
1114 + Currently polynomial is a reversed mapped distortion.
1115 + Order 1.5 is fudged to map into a bilinear distortion.
1116 though it is not the same order as that distortion.
1117 */
1118 double
1119 **matrix,
1120 **vectors,
1121 *terms;
1122
1123 size_t
1124 nterms; /* number of polynomial terms per number_values */
1125
1126 ssize_t
1127 j;
1128
1129 MagickBooleanType
1130 status;
1131
1132 /* first two coefficients hold polynomial order information */
1133 coeff[0] = arguments[0];
1134 coeff[1] = (double) poly_number_terms(arguments[0]);
1135 nterms = (size_t) coeff[1];
1136
1137 /* create matrix, a fake vectors matrix, and least sqs terms */
1138 matrix=AcquireMagickMatrix(nterms,nterms);
1139 vectors=(double **) AcquireQuantumMemory(number_values,
1140 sizeof(*vectors));
1141 terms=(double *) AcquireQuantumMemory(nterms,sizeof(*terms));
1142 if ((matrix == (double **) NULL) || (vectors == (double **) NULL) ||
1143 (terms == (double *) NULL))
1144 {
1145 matrix = RelinquishMagickMatrix(matrix, nterms);
1146 vectors = (double **) RelinquishMagickMemory(vectors);
1147 terms = (double *) RelinquishMagickMemory(terms);
1148 coeff = (double *) RelinquishMagickMemory(coeff);
1149 (void) ThrowMagickException(exception,GetMagickModule(),
1150 ResourceLimitError,"MemoryAllocationFailed",
1151 "%s", "DistortCoefficients");
1152 return((double *) NULL);
1153 }
1154 /* fake a number_values x3 vectors matrix from coefficients array */
1155 for (i=0; i < number_values; i++)
1156 vectors[i] = &(coeff[2+i*nterms]);
1157 /* Add given control point pairs for least squares solving */
1158 for (i=1; i < number_arguments; i+=cp_size) { /* NB: start = 1 not 0 */
1159 for (j=0; j < (ssize_t) nterms; j++)
1160 terms[j] = poly_basis_fn(j,arguments[i+cp_x],arguments[i+cp_y]);
1161 LeastSquaresAddTerms(matrix,vectors,terms,
1162 &(arguments[i+cp_values]),nterms,number_values);
1163 }
1164 terms = (double *) RelinquishMagickMemory(terms);
1165 /* Solve for LeastSquares Coefficients */
1166 status=GaussJordanElimination(matrix,vectors,nterms,number_values);
1167 matrix = RelinquishMagickMatrix(matrix, nterms);
1168 vectors = (double **) RelinquishMagickMemory(vectors);
1169 if ( status == MagickFalse ) {
1170 coeff = (double *) RelinquishMagickMemory(coeff);
1171 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1172 "InvalidArgument","%s : 'Unsolvable Matrix'",
1173 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1174 return((double *) NULL);
1175 }
1176 return(coeff);
1177 }
1178 case ArcDistortion:
1179 {
1180 /* Arc Distortion
1181 Args: arc_width rotate top_edge_radius bottom_edge_radius
1182 All but first argument are optional
1183 arc_width The angle over which to arc the image side-to-side
1184 rotate Angle to rotate image from vertical center
1185 top_radius Set top edge of source image at this radius
1186 bottom_radius Set bottom edge to this radius (radial scaling)
1187
1188 By default, if the radii arguments are nor provided the image radius
1189 is calculated so the horizontal center-line is fits the given arc
1190 without scaling.
1191
1192 The output image size is ALWAYS adjusted to contain the whole image,
1193 and an offset is given to position image relative to the 0,0 point of
1194 the origin, allowing users to use relative positioning onto larger
1195 background (via -flatten).
1196
1197 The arguments are converted to these coefficients
1198 c0: angle for center of source image
1199 c1: angle scale for mapping to source image
1200 c2: radius for top of source image
1201 c3: radius scale for mapping source image
1202 c4: centerline of arc within source image
1203
1204 Note the coefficients use a center angle, so asymptotic join is
1205 furthest from both sides of the source image. This also means that
1206 for arc angles greater than 360 the sides of the image will be
1207 trimmed equally.
1208
1209 Arc Distortion Notes...
1210 + Does not use a set of CPs
1211 + Will only work with Image Distortion
1212 + Can not be used for generating a sparse gradient (interpolation)
1213 */
1214 if ( number_arguments >= 1 && arguments[0] < MagickEpsilon ) {
1215 coeff = (double *) RelinquishMagickMemory(coeff);
1216 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1217 "InvalidArgument","%s : 'Arc Angle Too Small'",
1218 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1219 return((double *) NULL);
1220 }
1221 if ( number_arguments >= 3 && arguments[2] < MagickEpsilon ) {
1222 coeff = (double *) RelinquishMagickMemory(coeff);
1223 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1224 "InvalidArgument","%s : 'Outer Radius Too Small'",
1225 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1226 return((double *) NULL);
1227 }
1228 coeff[0] = -MagickPI2; /* -90, place at top! */
1229 if ( number_arguments >= 1 )
1230 coeff[1] = DegreesToRadians(arguments[0]);
1231 else
1232 coeff[1] = MagickPI2; /* zero arguments - center is at top */
1233 if ( number_arguments >= 2 )
1234 coeff[0] += DegreesToRadians(arguments[1]);
1235 coeff[0] /= Magick2PI; /* normalize radians */
1236 coeff[0] -= MagickRound(coeff[0]);
1237 coeff[0] *= Magick2PI; /* de-normalize back to radians */
1238 coeff[3] = (double)image->rows-1;
1239 coeff[2] = (double)image->columns/coeff[1] + coeff[3]/2.0;
1240 if ( number_arguments >= 3 ) {
1241 if ( number_arguments >= 4 )
1242 coeff[3] = arguments[2] - arguments[3];
1243 else
1244 coeff[3] *= arguments[2]/coeff[2];
1245 coeff[2] = arguments[2];
1246 }
1247 coeff[4] = ((double)image->columns-1.0)/2.0;
1248
1249 return(coeff);
1250 }
1251 case PolarDistortion:
1252 case DePolarDistortion:
1253 {
1254 /* (De)Polar Distortion (same set of arguments)
1255 Args: Rmax, Rmin, Xcenter,Ycenter, Afrom,Ato
1256 DePolar can also have the extra arguments of Width, Height
1257
1258 Coefficients 0 to 5 is the sanitized version first 6 input args
1259 Coefficient 6 is the angle to coord ratio and visa-versa
1260 Coefficient 7 is the radius to coord ratio and visa-versa
1261
1262 WARNING: It is possible for Radius max<min and/or Angle from>to
1263 */
1264 if ( number_arguments == 3
1265 || ( number_arguments > 6 && *method == PolarDistortion )
1266 || number_arguments > 8 ) {
1267 (void) ThrowMagickException(exception,GetMagickModule(),
1268 OptionError,"InvalidArgument", "%s : number of arguments",
1269 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1270 coeff=(double *) RelinquishMagickMemory(coeff);
1271 return((double *) NULL);
1272 }
1273 /* Rmax - if 0 calculate appropriate value */
1274 if ( number_arguments >= 1 )
1275 coeff[0] = arguments[0];
1276 else
1277 coeff[0] = 0.0;
1278 /* Rmin - usually 0 */
1279 coeff[1] = number_arguments >= 2 ? arguments[1] : 0.0;
1280 /* Center X,Y */
1281 if ( number_arguments >= 4 ) {
1282 coeff[2] = arguments[2];
1283 coeff[3] = arguments[3];
1284 }
1285 else { /* center of actual image */
1286 coeff[2] = (double)(image->columns)/2.0+image->page.x;
1287 coeff[3] = (double)(image->rows)/2.0+image->page.y;
1288 }
1289 /* Angle from,to - about polar center 0 is downward */
1290 coeff[4] = -MagickPI;
1291 if ( number_arguments >= 5 )
1292 coeff[4] = DegreesToRadians(arguments[4]);
1293 coeff[5] = coeff[4];
1294 if ( number_arguments >= 6 )
1295 coeff[5] = DegreesToRadians(arguments[5]);
1296 if ( fabs(coeff[4]-coeff[5]) < MagickEpsilon )
1297 coeff[5] += Magick2PI; /* same angle is a full circle */
1298 /* if radius 0 or negative, its a special value... */
1299 if ( coeff[0] < MagickEpsilon ) {
1300 /* Use closest edge if radius == 0 */
1301 if ( fabs(coeff[0]) < MagickEpsilon ) {
1302 coeff[0]=MagickMin(fabs(coeff[2]-image->page.x),
1303 fabs(coeff[3]-image->page.y));
1304 coeff[0]=MagickMin(coeff[0],
1305 fabs(coeff[2]-image->page.x-image->columns));
1306 coeff[0]=MagickMin(coeff[0],
1307 fabs(coeff[3]-image->page.y-image->rows));
1308 }
1309 /* furthest diagonal if radius == -1 */
1310 if ( fabs(-1.0-coeff[0]) < MagickEpsilon ) {
1311 double rx,ry;
1312 rx = coeff[2]-image->page.x;
1313 ry = coeff[3]-image->page.y;
1314 coeff[0] = rx*rx+ry*ry;
1315 ry = coeff[3]-image->page.y-image->rows;
1316 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1317 rx = coeff[2]-image->page.x-image->columns;
1318 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1319 ry = coeff[3]-image->page.y;
1320 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1321 coeff[0] = sqrt(coeff[0]);
1322 }
1323 }
1324 /* IF Rmax <= 0 or Rmin < 0 OR Rmax < Rmin, THEN error */
1325 if ( coeff[0] < MagickEpsilon || coeff[1] < -MagickEpsilon
1326 || (coeff[0]-coeff[1]) < MagickEpsilon ) {
1327 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1328 "InvalidArgument", "%s : Invalid Radius",
1329 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1330 coeff=(double *) RelinquishMagickMemory(coeff);
1331 return((double *) NULL);
1332 }
1333 /* conversion ratios */
1334 if ( *method == PolarDistortion ) {
1335 coeff[6]=(double) image->columns/(coeff[5]-coeff[4]);
1336 coeff[7]=(double) image->rows/(coeff[0]-coeff[1]);
1337 }
1338 else { /* *method == DePolarDistortion */
1339 coeff[6]=(coeff[5]-coeff[4])/image->columns;
1340 coeff[7]=(coeff[0]-coeff[1])/image->rows;
1341 }
1342 return(coeff);
1343 }
1344 case Cylinder2PlaneDistortion:
1345 case Plane2CylinderDistortion:
1346 {
1347 /* 3D Cylinder to/from a Tangential Plane
1348
1349 Projection between a cylinder and flat plain from a point on the
1350 center line of the cylinder.
1351
1352 The two surfaces coincide in 3D space at the given centers of
1353 distortion (perpendicular to projection point) on both images.
1354
1355 Args: FOV_arc_width
1356 Coefficients: FOV(radians), Radius, center_x,y, dest_center_x,y
1357
1358 FOV (Field Of View) the angular field of view of the distortion,
1359 across the width of the image, in degrees. The centers are the
1360 points of least distortion in the input and resulting images.
1361
1362 These centers are however determined later.
1363
1364 Coeff 0 is the FOV angle of view of image width in radians
1365 Coeff 1 is calculated radius of cylinder.
1366 Coeff 2,3 center of distortion of input image
1367 Coefficients 4,5 Center of Distortion of dest (determined later)
1368 */
1369 if (number_arguments < 1) {
1370 coeff = (double *) RelinquishMagickMemory(coeff);
1371 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1372 "InvalidArgument", "%s : 'Needs at least 1 argument'",
1373 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1374 return((double *) NULL);
1375 }
1376 if ( arguments[0] < MagickEpsilon || arguments[0] > 160.0 ) {
1377 coeff=(double *) RelinquishMagickMemory(coeff);
1378 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1379 "InvalidArgument", "%s : Invalid FOV Angle",
1380 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1381 return((double *) NULL);
1382 }
1383 coeff[0] = DegreesToRadians(arguments[0]);
1384 if ( *method == Cylinder2PlaneDistortion )
1385 /* image is curved around cylinder, so FOV angle (in radians)
1386 * scales directly to image X coordinate, according to its radius.
1387 */
1388 coeff[1] = (double) image->columns/coeff[0];
1389 else
1390 /* radius is distance away from an image with this angular FOV */
1391 coeff[1] = (double) image->columns / ( 2 * tan(coeff[0]/2) );
1392
1393 coeff[2] = (double)(image->columns)/2.0+image->page.x;
1394 coeff[3] = (double)(image->rows)/2.0+image->page.y;
1395 coeff[4] = coeff[2];
1396 coeff[5] = coeff[3]; /* assuming image size is the same */
1397 return(coeff);
1398 }
1399 case BarrelDistortion:
1400 case BarrelInverseDistortion:
1401 {
1402 /* Barrel Distortion
1403 Rs=(A*Rd^3 + B*Rd^2 + C*Rd + D)*Rd
1404 BarrelInv Distortion
1405 Rs=Rd/(A*Rd^3 + B*Rd^2 + C*Rd + D)
1406
1407 Where Rd is the normalized radius from corner to middle of image
1408 Input Arguments are one of the following forms (number of arguments)...
1409 3: A,B,C
1410 4: A,B,C,D
1411 5: A,B,C X,Y
1412 6: A,B,C,D X,Y
1413 8: Ax,Bx,Cx,Dx Ay,By,Cy,Dy
1414 10: Ax,Bx,Cx,Dx Ay,By,Cy,Dy X,Y
1415
1416 Returns 10 coefficient values, which are de-normalized (pixel scale)
1417 Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Xc, Yc
1418 */
1419 /* Radius de-normalization scaling factor */
1420 double
1421 rscale = 2.0/MagickMin((double) image->columns,(double) image->rows);
1422
1423 /* sanity check number of args must = 3,4,5,6,8,10 or error */
1424 if ( (number_arguments < 3) || (number_arguments == 7) ||
1425 (number_arguments == 9) || (number_arguments > 10) )
1426 {
1427 coeff=(double *) RelinquishMagickMemory(coeff);
1428 (void) ThrowMagickException(exception,GetMagickModule(),
1429 OptionError,"InvalidArgument", "%s : number of arguments",
1430 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1431 return((double *) NULL);
1432 }
1433 /* A,B,C,D coefficients */
1434 coeff[0] = arguments[0];
1435 coeff[1] = arguments[1];
1436 coeff[2] = arguments[2];
1437 if ((number_arguments == 3) || (number_arguments == 5) )
1438 coeff[3] = 1.0 - coeff[0] - coeff[1] - coeff[2];
1439 else
1440 coeff[3] = arguments[3];
1441 /* de-normalize the coefficients */
1442 coeff[0] *= pow(rscale,3.0);
1443 coeff[1] *= rscale*rscale;
1444 coeff[2] *= rscale;
1445 /* Y coefficients: as given OR same as X coefficients */
1446 if ( number_arguments >= 8 ) {
1447 coeff[4] = arguments[4] * pow(rscale,3.0);
1448 coeff[5] = arguments[5] * rscale*rscale;
1449 coeff[6] = arguments[6] * rscale;
1450 coeff[7] = arguments[7];
1451 }
1452 else {
1453 coeff[4] = coeff[0];
1454 coeff[5] = coeff[1];
1455 coeff[6] = coeff[2];
1456 coeff[7] = coeff[3];
1457 }
1458 /* X,Y Center of Distortion (image coordinates) */
1459 if ( number_arguments == 5 ) {
1460 coeff[8] = arguments[3];
1461 coeff[9] = arguments[4];
1462 }
1463 else if ( number_arguments == 6 ) {
1464 coeff[8] = arguments[4];
1465 coeff[9] = arguments[5];
1466 }
1467 else if ( number_arguments == 10 ) {
1468 coeff[8] = arguments[8];
1469 coeff[9] = arguments[9];
1470 }
1471 else {
1472 /* center of the image provided (image coordinates) */
1473 coeff[8] = (double)image->columns/2.0 + image->page.x;
1474 coeff[9] = (double)image->rows/2.0 + image->page.y;
1475 }
1476 return(coeff);
1477 }
1478 case ShepardsDistortion:
1479 {
1480 /* Shepards Distortion input arguments are the coefficients!
1481 Just check the number of arguments is valid!
1482 Args: u1,v1, x1,y1, ...
1483 OR : u1,v1, r1,g1,c1, ...
1484 */
1485 if ( number_arguments%cp_size != 0 ||
1486 number_arguments < cp_size ) {
1487 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1488 "InvalidArgument", "%s : 'requires CP's (4 numbers each)'",
1489 CommandOptionToMnemonic(MagickDistortOptions, *method));
1490 coeff=(double *) RelinquishMagickMemory(coeff);
1491 return((double *) NULL);
1492 }
1493 /* User defined weighting power for Shepard's Method */
1494 { const char *artifact=GetImageArtifact(image,"shepards:power");
1495 if ( artifact != (const char *) NULL ) {
1496 coeff[0]=StringToDouble(artifact,(char **) NULL) / 2.0;
1497 if ( coeff[0] < MagickEpsilon ) {
1498 (void) ThrowMagickException(exception,GetMagickModule(),
1499 OptionError,"InvalidArgument","%s", "-define shepards:power" );
1500 coeff=(double *) RelinquishMagickMemory(coeff);
1501 return((double *) NULL);
1502 }
1503 }
1504 else
1505 coeff[0]=1.0; /* Default power of 2 (Inverse Squared) */
1506 }
1507 return(coeff);
1508 }
1509 default:
1510 break;
1511 }
1512 /* you should never reach this point */
1513 perror("no method handler"); /* just fail assertion */
1514 return((double *) NULL);
1515}
1516␌
1517/*
1518%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1519% %
1520% %
1521% %
1522+ D i s t o r t R e s i z e I m a g e %
1523% %
1524% %
1525% %
1526%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1527%
1528% DistortResizeImage() resize image using the equivalent but slower image
1529% distortion operator. The filter is applied using a EWA cylindrical
1530% resampling. But like resize the final image size is limited to whole pixels
1531% with no effects by virtual-pixels on the result.
1532%
1533% Note that images containing a transparency channel will be twice as slow to
1534% resize as images one without transparency.
1535%
1536% The format of the DistortResizeImage method is:
1537%
1538% Image *DistortResizeImage(const Image *image,const size_t columns,
1539% const size_t rows,ExceptionInfo *exception)
1540%
1541% A description of each parameter follows:
1542%
1543% o image: the image.
1544%
1545% o columns: the number of columns in the resized image.
1546%
1547% o rows: the number of rows in the resized image.
1548%
1549% o exception: return any errors or warnings in this structure.
1550%
1551*/
1552MagickExport Image *DistortResizeImage(const Image *image,const size_t columns,
1553 const size_t rows,ExceptionInfo *exception)
1554{
1555#define DistortResizeImageTag "Distort/Image"
1556
1557 Image
1558 *resize_image,
1559 *tmp_image;
1560
1561 RectangleInfo
1562 crop_area;
1563
1564 double
1565 distort_args[12];
1566
1567 VirtualPixelMethod
1568 vp_save;
1569
1570 /*
1571 Distort resize image.
1572 */
1573 assert(image != (const Image *) NULL);
1574 assert(image->signature == MagickCoreSignature);
1575 assert(exception != (ExceptionInfo *) NULL);
1576 assert(exception->signature == MagickCoreSignature);
1577 if (IsEventLogging() != MagickFalse)
1578 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1579 if ((columns == 0) || (rows == 0))
1580 return((Image *) NULL);
1581 /* Do not short-circuit this resize if final image size is unchanged */
1582
1583 (void) memset(distort_args,0,sizeof(distort_args));
1584 distort_args[4]=(double) image->columns;
1585 distort_args[6]=(double) columns;
1586 distort_args[9]=(double) image->rows;
1587 distort_args[11]=(double) rows;
1588
1589 vp_save=GetImageVirtualPixelMethod(image);
1590
1591 tmp_image=CloneImage(image,0,0,MagickTrue,exception);
1592 if (tmp_image == (Image *) NULL)
1593 return((Image *) NULL);
1594 (void) SetImageVirtualPixelMethod(tmp_image,TransparentVirtualPixelMethod,
1595 exception);
1596
1597 if ((image->alpha_trait & BlendPixelTrait) == 0)
1598 {
1599 /*
1600 Image has no alpha channel, so we are free to use it.
1601 */
1602 (void) SetImageAlphaChannel(tmp_image,SetAlphaChannel,exception);
1603 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1604 MagickTrue,exception),
1605 tmp_image=DestroyImage(tmp_image);
1606 if (resize_image == (Image *) NULL)
1607 return((Image *) NULL);
1608 (void) SetImageAlphaChannel(resize_image,OffAlphaChannel,exception);
1609 }
1610 else
1611 {
1612 /*
1613 Image has transparency so handle colors and alpha separately.
1614 Basically we need to separate Virtual-Pixel alpha in the resized
1615 image, so only the actual original images alpha channel is used.
1616
1617 distort alpha channel separately
1618 */
1619 Image
1620 *resize_alpha;
1621
1622 (void) SetImageAlphaChannel(tmp_image,ExtractAlphaChannel,exception);
1623 (void) SetImageAlphaChannel(tmp_image,OpaqueAlphaChannel,exception);
1624 resize_alpha=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1625 MagickTrue,exception),
1626 tmp_image=DestroyImage(tmp_image);
1627 if (resize_alpha == (Image *) NULL)
1628 return((Image *) NULL);
1629
1630 /* distort the actual image containing alpha + VP alpha */
1631 tmp_image=CloneImage(image,0,0,MagickTrue,exception);
1632 if (tmp_image == (Image *) NULL)
1633 return((Image *) NULL);
1634 (void) SetImageVirtualPixelMethod(tmp_image,
1635 TransparentVirtualPixelMethod,exception);
1636 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1637 MagickTrue,exception),
1638 tmp_image=DestroyImage(tmp_image);
1639 if (resize_image == (Image *) NULL)
1640 {
1641 resize_alpha=DestroyImage(resize_alpha);
1642 return((Image *) NULL);
1643 }
1644 /* replace resize images alpha with the separately distorted alpha */
1645 (void) SetImageAlphaChannel(resize_image,OffAlphaChannel,exception);
1646 (void) SetImageAlphaChannel(resize_alpha,OffAlphaChannel,exception);
1647 (void) CompositeImage(resize_image,resize_alpha,CopyAlphaCompositeOp,
1648 MagickTrue,0,0,exception);
1649 resize_alpha=DestroyImage(resize_alpha);
1650 resize_image->alpha_trait=image->alpha_trait;
1651 resize_image->compose=image->compose;
1652 }
1653 (void) SetImageVirtualPixelMethod(resize_image,vp_save,exception);
1654
1655 /*
1656 Clean up the results of the Distortion
1657 */
1658 crop_area.width=columns;
1659 crop_area.height=rows;
1660 crop_area.x=0;
1661 crop_area.y=0;
1662
1663 tmp_image=resize_image;
1664 resize_image=CropImage(tmp_image,&crop_area,exception);
1665 tmp_image=DestroyImage(tmp_image);
1666 if (resize_image != (Image *) NULL)
1667 {
1668 resize_image->page.width=0;
1669 resize_image->page.height=0;
1670 }
1671 return(resize_image);
1672}
1673␌
1674/*
1675%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1676% %
1677% %
1678% %
1679% D i s t o r t I m a g e %
1680% %
1681% %
1682% %
1683%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1684%
1685% DistortImage() distorts an image using various distortion methods, by
1686% mapping color lookups of the source image to a new destination image
1687% usually of the same size as the source image, unless 'bestfit' is set to
1688% true.
1689%
1690% If 'bestfit' is enabled, and distortion allows it, the destination image is
1691% adjusted to ensure the whole source 'image' will just fit within the final
1692% destination image, which will be sized and offset accordingly. Also in
1693% many cases the virtual offset of the source image will be taken into
1694% account in the mapping.
1695%
1696% If the '-verbose' control option has been set print to standard error the
1697% equivalent '-fx' formula with coefficients for the function, if practical.
1698%
1699% The format of the DistortImage() method is:
1700%
1701% Image *DistortImage(const Image *image,const DistortMethod method,
1702% const size_t number_arguments,const double *arguments,
1703% MagickBooleanType bestfit, ExceptionInfo *exception)
1704%
1705% A description of each parameter follows:
1706%
1707% o image: the image to be distorted.
1708%
1709% o method: the method of image distortion.
1710%
1711% ArcDistortion always ignores source image offset, and always
1712% 'bestfit' the destination image with the top left corner offset
1713% relative to the polar mapping center.
1714%
1715% Affine, Perspective, and Bilinear, do least squares fitting of the
1716% distortion when more than the minimum number of control point pairs
1717% are provided.
1718%
1719% Perspective, and Bilinear, fall back to a Affine distortion when less
1720% than 4 control point pairs are provided. While Affine distortions
1721% let you use any number of control point pairs, that is Zero pairs is
1722% a No-Op (viewport only) distortion, one pair is a translation and
1723% two pairs of control points do a scale-rotate-translate, without any
1724% shearing.
1725%
1726% o number_arguments: the number of arguments given.
1727%
1728% o arguments: an array of floating point arguments for this method.
1729%
1730% o bestfit: Attempt to 'bestfit' the size of the resulting image.
1731% This also forces the resulting image to be a 'layered' virtual
1732% canvas image. Can be overridden using 'distort:viewport' setting.
1733%
1734% o exception: return any errors or warnings in this structure
1735%
1736% Extra Controls from Image meta-data (artifacts)...
1737%
1738% o "verbose"
1739% Output to stderr alternatives, internal coefficients, and FX
1740% equivalents for the distortion operation (if feasible).
1741% This forms an extra check of the distortion method, and allows users
1742% access to the internal constants IM calculates for the distortion.
1743%
1744% o "distort:viewport"
1745% Directly set the output image canvas area and offset to use for the
1746% resulting image, rather than use the original images canvas, or a
1747% calculated 'bestfit' canvas.
1748%
1749% o "distort:scale"
1750% Scale the size of the output canvas by this amount to provide a
1751% method of Zooming, and for super-sampling the results.
1752%
1753% Other settings that can effect results include
1754%
1755% o 'interpolate' For source image lookups (scale enlargements)
1756%
1757% o 'filter' Set filter to use for area-resampling (scale shrinking).
1758% Set to 'point' to turn off and use 'interpolate' lookup
1759% instead
1760%
1761*/
1762MagickExport Image *DistortImage(const Image *image, DistortMethod method,
1763 const size_t number_arguments,const double *arguments,
1764 MagickBooleanType bestfit,ExceptionInfo *exception)
1765{
1766#define DistortImageTag "Distort/Image"
1767
1768 double
1769 *coeff,
1770 output_scaling;
1771
1772 Image
1773 *distort_image;
1774
1775 RectangleInfo
1776 geometry; /* geometry of the distorted space viewport */
1777
1778 MagickBooleanType
1779 viewport_given;
1780
1781 PixelInfo
1782 invalid; /* the color to assign when distort result is invalid */
1783
1784 assert(image != (Image *) NULL);
1785 assert(image->signature == MagickCoreSignature);
1786 assert(exception != (ExceptionInfo *) NULL);
1787 assert(exception->signature == MagickCoreSignature);
1788 if (IsEventLogging() != MagickFalse)
1789 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1790 /*
1791 Handle Special Compound Distortions
1792 */
1793 if ( method == ResizeDistortion )
1794 {
1795 if ( number_arguments != 2 )
1796 {
1797 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1798 "InvalidArgument","%s : '%s'","Resize",
1799 "Invalid number of args: 2 only");
1800 return((Image *) NULL);
1801 }
1802 distort_image=DistortResizeImage(image,(size_t)arguments[0],
1803 (size_t)arguments[1], exception);
1804 return(distort_image);
1805 }
1806
1807 /*
1808 Convert input arguments (usually as control points for reverse mapping)
1809 into mapping coefficients to apply the distortion.
1810
1811 Note that some distortions are mapped to other distortions,
1812 and as such do not require specific code after this point.
1813 */
1814 coeff = GenerateCoefficients(image, &method, number_arguments,
1815 arguments, 0, exception);
1816 if ( coeff == (double *) NULL )
1817 return((Image *) NULL);
1818
1819 /*
1820 Determine the size and offset for a 'bestfit' destination.
1821 Usually the four corners of the source image is enough.
1822 */
1823
1824 /* default output image bounds, when no 'bestfit' is requested */
1825 geometry.width=image->columns;
1826 geometry.height=image->rows;
1827 geometry.x=0;
1828 geometry.y=0;
1829
1830 if ( method == ArcDistortion ) {
1831 bestfit = MagickTrue; /* always calculate a 'best fit' viewport */
1832 }
1833
1834 /* Work out the 'best fit', (required for ArcDistortion) */
1835 if ( bestfit ) {
1836 PointInfo
1837 s,d,min,max; /* source, dest coords --mapping--> min, max coords */
1838
1839 MagickBooleanType
1840 fix_bounds = MagickTrue; /* enlarge bounds for VP handling */
1841
1842 s.x=s.y=min.x=max.x=min.y=max.y=0.0; /* keep compiler happy */
1843
1844/* defines to figure out the bounds of the distorted image */
1845#define InitalBounds(p) \
1846{ \
1847 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1848 min.x = max.x = p.x; \
1849 min.y = max.y = p.y; \
1850}
1851#define ExpandBounds(p) \
1852{ \
1853 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1854 min.x = MagickMin(min.x,p.x); \
1855 max.x = MagickMax(max.x,p.x); \
1856 min.y = MagickMin(min.y,p.y); \
1857 max.y = MagickMax(max.y,p.y); \
1858}
1859
1860 switch (method)
1861 {
1862 case AffineDistortion:
1863 case RigidAffineDistortion:
1864 { double inverse[6];
1865 InvertAffineCoefficients(coeff, inverse);
1866 s.x = (double) image->page.x;
1867 s.y = (double) image->page.y;
1868 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1869 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1870 InitalBounds(d);
1871 s.x = (double) image->page.x+image->columns;
1872 s.y = (double) image->page.y;
1873 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1874 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1875 ExpandBounds(d);
1876 s.x = (double) image->page.x;
1877 s.y = (double) image->page.y+image->rows;
1878 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1879 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1880 ExpandBounds(d);
1881 s.x = (double) image->page.x+image->columns;
1882 s.y = (double) image->page.y+image->rows;
1883 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1884 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1885 ExpandBounds(d);
1886 break;
1887 }
1888 case PerspectiveDistortion:
1889 { double inverse[8], scale;
1890 InvertPerspectiveCoefficients(coeff, inverse);
1891 s.x = (double) image->page.x;
1892 s.y = (double) image->page.y;
1893 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1894 scale=PerceptibleReciprocal(scale);
1895 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1896 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1897 InitalBounds(d);
1898 s.x = (double) image->page.x+image->columns;
1899 s.y = (double) image->page.y;
1900 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1901 scale=PerceptibleReciprocal(scale);
1902 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1903 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1904 ExpandBounds(d);
1905 s.x = (double) image->page.x;
1906 s.y = (double) image->page.y+image->rows;
1907 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1908 scale=PerceptibleReciprocal(scale);
1909 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1910 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1911 ExpandBounds(d);
1912 s.x = (double) image->page.x+image->columns;
1913 s.y = (double) image->page.y+image->rows;
1914 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1915 scale=PerceptibleReciprocal(scale);
1916 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1917 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1918 ExpandBounds(d);
1919 break;
1920 }
1921 case ArcDistortion:
1922 { double a, ca, sa;
1923 /* Forward Map Corners */
1924 a = coeff[0]-coeff[1]/2; ca = cos(a); sa = sin(a);
1925 d.x = coeff[2]*ca;
1926 d.y = coeff[2]*sa;
1927 InitalBounds(d);
1928 d.x = (coeff[2]-coeff[3])*ca;
1929 d.y = (coeff[2]-coeff[3])*sa;
1930 ExpandBounds(d);
1931 a = coeff[0]+coeff[1]/2; ca = cos(a); sa = sin(a);
1932 d.x = coeff[2]*ca;
1933 d.y = coeff[2]*sa;
1934 ExpandBounds(d);
1935 d.x = (coeff[2]-coeff[3])*ca;
1936 d.y = (coeff[2]-coeff[3])*sa;
1937 ExpandBounds(d);
1938 /* Orthogonal points along top of arc */
1939 for( a=(double) (ceil((double) ((coeff[0]-coeff[1]/2.0)/MagickPI2))*MagickPI2);
1940 a<(coeff[0]+coeff[1]/2.0); a+=MagickPI2 ) {
1941 ca = cos(a); sa = sin(a);
1942 d.x = coeff[2]*ca;
1943 d.y = coeff[2]*sa;
1944 ExpandBounds(d);
1945 }
1946 /*
1947 Convert the angle_to_width and radius_to_height
1948 to appropriate scaling factors, to allow faster processing
1949 in the mapping function.
1950 */
1951 coeff[1] = (double) (Magick2PI*image->columns/coeff[1]);
1952 coeff[3] = (double)image->rows/coeff[3];
1953 break;
1954 }
1955 case PolarDistortion:
1956 {
1957 if (number_arguments < 2)
1958 coeff[2] = coeff[3] = 0.0;
1959 min.x = coeff[2]-coeff[0];
1960 max.x = coeff[2]+coeff[0];
1961 min.y = coeff[3]-coeff[0];
1962 max.y = coeff[3]+coeff[0];
1963 /* should be about 1.0 if Rmin = 0 */
1964 coeff[7]=(double) geometry.height/(coeff[0]-coeff[1]);
1965 break;
1966 }
1967 case DePolarDistortion:
1968 {
1969 /* direct calculation as it needs to tile correctly
1970 * for reversibility in a DePolar-Polar cycle */
1971 fix_bounds = MagickFalse;
1972 geometry.x = geometry.y = 0;
1973 geometry.height = (size_t) ceil(coeff[0]-coeff[1]);
1974 geometry.width = (size_t) ceil((coeff[0]-coeff[1])*
1975 (coeff[5]-coeff[4])*0.5);
1976 /* correct scaling factors relative to new size */
1977 coeff[6]=(coeff[5]-coeff[4])*PerceptibleReciprocal(geometry.width); /* changed width */
1978 coeff[7]=(coeff[0]-coeff[1])*PerceptibleReciprocal(geometry.height); /* should be about 1.0 */
1979 break;
1980 }
1981 case Cylinder2PlaneDistortion:
1982 {
1983 /* direct calculation so center of distortion is either a pixel
1984 * center, or pixel edge. This allows for reversibility of the
1985 * distortion */
1986 geometry.x = geometry.y = 0;
1987 geometry.width = (size_t) ceil( 2.0*coeff[1]*tan(coeff[0]/2.0) );
1988 geometry.height = (size_t) ceil( 2.0*coeff[3]/cos(coeff[0]/2.0) );
1989 /* correct center of distortion relative to new size */
1990 coeff[4] = (double) geometry.width/2.0;
1991 coeff[5] = (double) geometry.height/2.0;
1992 fix_bounds = MagickFalse;
1993 break;
1994 }
1995 case Plane2CylinderDistortion:
1996 {
1997 /* direct calculation center is either pixel center, or pixel edge
1998 * so as to allow reversibility of the image distortion */
1999 geometry.x = geometry.y = 0;
2000 geometry.width = (size_t) ceil(coeff[0]*coeff[1]); /* FOV * radius */
2001 geometry.height = (size_t) (2*coeff[3]); /* input image height */
2002 /* correct center of distortion relative to new size */
2003 coeff[4] = (double) geometry.width/2.0;
2004 coeff[5] = (double) geometry.height/2.0;
2005 fix_bounds = MagickFalse;
2006 break;
2007 }
2008 case ShepardsDistortion:
2009 case BilinearForwardDistortion:
2010 case BilinearReverseDistortion:
2011#if 0
2012 case QuadrilateralDistortion:
2013#endif
2014 case PolynomialDistortion:
2015 case BarrelDistortion:
2016 case BarrelInverseDistortion:
2017 default:
2018 /* no calculated bestfit available for these distortions */
2019 bestfit = MagickFalse;
2020 fix_bounds = MagickFalse;
2021 break;
2022 }
2023
2024 /* Set the output image geometry to calculated 'bestfit'.
2025 Yes this tends to 'over do' the file image size, ON PURPOSE!
2026 Do not do this for DePolar which needs to be exact for virtual tiling.
2027 */
2028 if ( fix_bounds ) {
2029 geometry.x = (ssize_t) floor(min.x-0.5);
2030 geometry.y = (ssize_t) floor(min.y-0.5);
2031 geometry.width=(size_t) ceil(max.x-geometry.x+0.5);
2032 geometry.height=(size_t) ceil(max.y-geometry.y+0.5);
2033 }
2034
2035 } /* end bestfit destination image calculations */
2036
2037 /* The user provided a 'viewport' expert option which may
2038 overrides some parts of the current output image geometry.
2039 This also overrides its default 'bestfit' setting.
2040 */
2041 { const char *artifact=GetImageArtifact(image,"distort:viewport");
2042 viewport_given = MagickFalse;
2043 if ( artifact != (const char *) NULL ) {
2044 MagickStatusType flags=ParseAbsoluteGeometry(artifact,&geometry);
2045 if (flags==NoValue)
2046 (void) ThrowMagickException(exception,GetMagickModule(),
2047 OptionWarning,"InvalidSetting","'%s' '%s'",
2048 "distort:viewport",artifact);
2049 else
2050 viewport_given = MagickTrue;
2051 }
2052 }
2053
2054 /* Verbose output */
2055 if (IsStringTrue(GetImageArtifact(image,"verbose")) != MagickFalse) {
2056 ssize_t
2057 i;
2058 char image_gen[MagickPathExtent];
2059 const char *lookup;
2060
2061 /* Set destination image size and virtual offset */
2062 if ( bestfit || viewport_given ) {
2063 (void) FormatLocaleString(image_gen,MagickPathExtent,
2064 " -size %.20gx%.20g -page %+.20g%+.20g xc: +insert \\\n",
2065 (double) geometry.width,(double) geometry.height,(double) geometry.x,
2066 (double) geometry.y);
2067 lookup="v.p{xx-v.page.x-0.5,yy-v.page.y-0.5}";
2068 }
2069 else {
2070 image_gen[0] = '\0'; /* no destination to generate */
2071 lookup = "p{xx-page.x-0.5,yy-page.y-0.5}"; /* simplify lookup */
2072 }
2073
2074 switch (method)
2075 {
2076 case AffineDistortion:
2077 case RigidAffineDistortion:
2078 {
2079 double
2080 *inverse;
2081
2082 inverse=(double *) AcquireQuantumMemory(6,sizeof(*inverse));
2083 if (inverse == (double *) NULL)
2084 {
2085 coeff=(double *) RelinquishMagickMemory(coeff);
2086 (void) ThrowMagickException(exception,GetMagickModule(),
2087 ResourceLimitError,"MemoryAllocationFailed","%s","DistortImages");
2088 return((Image *) NULL);
2089 }
2090 InvertAffineCoefficients(coeff, inverse);
2091 CoefficientsToAffineArgs(inverse);
2092 (void) FormatLocaleFile(stderr, "Affine projection:\n");
2093 (void) FormatLocaleFile(stderr,
2094 " -distort AffineProjection \\\n '");
2095 for (i=0; i < 5; i++)
2096 (void) FormatLocaleFile(stderr, "%.*g,",GetMagickPrecision(),
2097 inverse[i]);
2098 (void) FormatLocaleFile(stderr, "%.*g'\n",GetMagickPrecision(),
2099 inverse[5]);
2100 (void) FormatLocaleFile(stderr,
2101 "Equivalent scale, rotation(deg), translation:\n");
2102 (void) FormatLocaleFile(stderr," %.*g,%.*g,%.*g,%.*g\n",
2103 GetMagickPrecision(),sqrt(inverse[0]*inverse[0]+
2104 inverse[1]*inverse[1]),GetMagickPrecision(),
2105 RadiansToDegrees(atan2(inverse[1],inverse[0])),
2106 GetMagickPrecision(),inverse[4],GetMagickPrecision(),inverse[5]);
2107 inverse=(double *) RelinquishMagickMemory(inverse);
2108 (void) FormatLocaleFile(stderr,"Affine distort, FX equivalent:\n");
2109 (void) FormatLocaleFile(stderr, "%s", image_gen);
2110 (void) FormatLocaleFile(stderr,
2111 " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2112 (void) FormatLocaleFile(stderr," xx=%+.*g*ii %+.*g*jj %+.*g;\n",
2113 GetMagickPrecision(),coeff[0],GetMagickPrecision(),coeff[1],
2114 GetMagickPrecision(),coeff[2]);
2115 (void) FormatLocaleFile(stderr," yy=%+.*g*ii %+.*g*jj %+.*g;\n",
2116 GetMagickPrecision(),coeff[3],GetMagickPrecision(),coeff[4],
2117 GetMagickPrecision(),coeff[5]);
2118 (void) FormatLocaleFile(stderr," %s' \\\n",lookup);
2119 break;
2120 }
2121 case PerspectiveDistortion:
2122 {
2123 double
2124 *inverse;
2125
2126 inverse=(double *) AcquireQuantumMemory(8,sizeof(*inverse));
2127 if (inverse == (double *) NULL)
2128 {
2129 coeff=(double *) RelinquishMagickMemory(coeff);
2130 (void) ThrowMagickException(exception,GetMagickModule(),
2131 ResourceLimitError,"MemoryAllocationFailed","%s",
2132 "DistortCoefficients");
2133 return((Image *) NULL);
2134 }
2135 InvertPerspectiveCoefficients(coeff, inverse);
2136 (void) FormatLocaleFile(stderr,"Perspective Projection:\n");
2137 (void) FormatLocaleFile(stderr,
2138 " -distort PerspectiveProjection \\\n '");
2139 for (i=0; i < 4; i++)
2140 (void) FormatLocaleFile(stderr, "%.*g, ",GetMagickPrecision(),
2141 inverse[i]);
2142 (void) FormatLocaleFile(stderr, "\n ");
2143 for ( ; i < 7; i++)
2144 (void) FormatLocaleFile(stderr, "%.*g, ",GetMagickPrecision(),
2145 inverse[i]);
2146 (void) FormatLocaleFile(stderr, "%.*g'\n",GetMagickPrecision(),
2147 inverse[7]);
2148 inverse=(double *) RelinquishMagickMemory(inverse);
2149 (void) FormatLocaleFile(stderr,"Perspective Distort, FX Equivalent:\n");
2150 (void) FormatLocaleFile(stderr,"%.1024s",image_gen);
2151 (void) FormatLocaleFile(stderr,
2152 " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2153 (void) FormatLocaleFile(stderr," rr=%+.*g*ii %+.*g*jj + 1;\n",
2154 GetMagickPrecision(),coeff[6],GetMagickPrecision(),coeff[7]);
2155 (void) FormatLocaleFile(stderr,
2156 " xx=(%+.*g*ii %+.*g*jj %+.*g)/rr;\n",
2157 GetMagickPrecision(),coeff[0],GetMagickPrecision(),coeff[1],
2158 GetMagickPrecision(),coeff[2]);
2159 (void) FormatLocaleFile(stderr,
2160 " yy=(%+.*g*ii %+.*g*jj %+.*g)/rr;\n",
2161 GetMagickPrecision(),coeff[3],GetMagickPrecision(),coeff[4],
2162 GetMagickPrecision(),coeff[5]);
2163 (void) FormatLocaleFile(stderr," rr%s0 ? %s : blue' \\\n",
2164 coeff[8] < 0.0 ? "<" : ">", lookup);
2165 break;
2166 }
2167 case BilinearForwardDistortion:
2168 {
2169 (void) FormatLocaleFile(stderr,"BilinearForward Mapping Equations:\n");
2170 (void) FormatLocaleFile(stderr,"%s", image_gen);
2171 (void) FormatLocaleFile(stderr," i = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
2172 coeff[0],coeff[1],coeff[2],coeff[3]);
2173 (void) FormatLocaleFile(stderr," j = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
2174 coeff[4],coeff[5],coeff[6],coeff[7]);
2175#if 0
2176 /* for debugging */
2177 (void) FormatLocaleFile(stderr, " c8 = %+lf c9 = 2*a = %+lf;\n",
2178 coeff[8], coeff[9]);
2179#endif
2180 (void) FormatLocaleFile(stderr,
2181 "BilinearForward Distort, FX Equivalent:\n");
2182 (void) FormatLocaleFile(stderr,"%s", image_gen);
2183 (void) FormatLocaleFile(stderr,
2184 " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",0.5-coeff[3],0.5-
2185 coeff[7]);
2186 (void) FormatLocaleFile(stderr," bb=%lf*ii %+lf*jj %+lf;\n",
2187 coeff[6], -coeff[2], coeff[8]);
2188 /* Handle Special degenerate (non-quadratic) or trapezoidal case */
2189 if (coeff[9] != 0)
2190 {
2191 (void) FormatLocaleFile(stderr,
2192 " rt=bb*bb %+lf*(%lf*ii%+lf*jj);\n",-2*coeff[9],coeff[4],
2193 -coeff[0]);
2194 (void) FormatLocaleFile(stderr,
2195 " yy=( -bb + sqrt(rt) ) / %lf;\n",coeff[9]);
2196 }
2197 else
2198 (void) FormatLocaleFile(stderr," yy=(%lf*ii%+lf*jj)/bb;\n",
2199 -coeff[4],coeff[0]);
2200 (void) FormatLocaleFile(stderr,
2201 " xx=(ii %+lf*yy)/(%lf %+lf*yy);\n",-coeff[1],coeff[0],
2202 coeff[2]);
2203 if ( coeff[9] != 0 )
2204 (void) FormatLocaleFile(stderr," (rt < 0 ) ? red : %s'\n",
2205 lookup);
2206 else
2207 (void) FormatLocaleFile(stderr," %s' \\\n", lookup);
2208 break;
2209 }
2210 case BilinearReverseDistortion:
2211 {
2212#if 0
2213 (void) FormatLocaleFile(stderr, "Polynomial Projection Distort:\n");
2214 (void) FormatLocaleFile(stderr, " -distort PolynomialProjection \\\n");
2215 (void) FormatLocaleFile(stderr, " '1.5, %lf, %lf, %lf, %lf,\n",
2216 coeff[3], coeff[0], coeff[1], coeff[2]);
2217 (void) FormatLocaleFile(stderr, " %lf, %lf, %lf, %lf'\n",
2218 coeff[7], coeff[4], coeff[5], coeff[6]);
2219#endif
2220 (void) FormatLocaleFile(stderr,
2221 "BilinearReverse Distort, FX Equivalent:\n");
2222 (void) FormatLocaleFile(stderr,"%s", image_gen);
2223 (void) FormatLocaleFile(stderr,
2224 " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2225 (void) FormatLocaleFile(stderr,
2226 " xx=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",coeff[0],coeff[1],
2227 coeff[2], coeff[3]);
2228 (void) FormatLocaleFile(stderr,
2229 " yy=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",coeff[4],coeff[5],
2230 coeff[6], coeff[7]);
2231 (void) FormatLocaleFile(stderr," %s' \\\n", lookup);
2232 break;
2233 }
2234 case PolynomialDistortion:
2235 {
2236 size_t nterms = (size_t) coeff[1];
2237 (void) FormatLocaleFile(stderr,
2238 "Polynomial (order %lg, terms %lu), FX Equivalent\n",coeff[0],
2239 (unsigned long) nterms);
2240 (void) FormatLocaleFile(stderr,"%s", image_gen);
2241 (void) FormatLocaleFile(stderr,
2242 " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2243 (void) FormatLocaleFile(stderr, " xx =");
2244 for (i=0; i < (ssize_t) nterms; i++)
2245 {
2246 if ((i != 0) && (i%4 == 0))
2247 (void) FormatLocaleFile(stderr, "\n ");
2248 (void) FormatLocaleFile(stderr," %+lf%s",coeff[2+i],
2249 poly_basis_str(i));
2250 }
2251 (void) FormatLocaleFile(stderr,";\n yy =");
2252 for (i=0; i < (ssize_t) nterms; i++)
2253 {
2254 if ((i != 0) && (i%4 == 0))
2255 (void) FormatLocaleFile(stderr,"\n ");
2256 (void) FormatLocaleFile(stderr," %+lf%s",coeff[2+i+(int) nterms],
2257 poly_basis_str(i));
2258 }
2259 (void) FormatLocaleFile(stderr,";\n %s' \\\n", lookup);
2260 break;
2261 }
2262 case ArcDistortion:
2263 {
2264 (void) FormatLocaleFile(stderr,"Arc Distort, Internal Coefficients:\n");
2265 for (i=0; i < 5; i++)
2266 (void) FormatLocaleFile(stderr,
2267 " c%.20g = %+lf\n",(double) i,coeff[i]);
2268 (void) FormatLocaleFile(stderr,"Arc Distort, FX Equivalent:\n");
2269 (void) FormatLocaleFile(stderr,"%s", image_gen);
2270 (void) FormatLocaleFile(stderr," -fx 'ii=i+page.x; jj=j+page.y;\n");
2271 (void) FormatLocaleFile(stderr," xx=(atan2(jj,ii)%+lf)/(2*pi);\n",
2272 -coeff[0]);
2273 (void) FormatLocaleFile(stderr," xx=xx-round(xx);\n");
2274 (void) FormatLocaleFile(stderr," xx=xx*%lf %+lf;\n",coeff[1],
2275 coeff[4]);
2276 (void) FormatLocaleFile(stderr,
2277 " yy=(%lf - hypot(ii,jj)) * %lf;\n",coeff[2],coeff[3]);
2278 (void) FormatLocaleFile(stderr," v.p{xx-.5,yy-.5}' \\\n");
2279 break;
2280 }
2281 case PolarDistortion:
2282 {
2283 (void) FormatLocaleFile(stderr,"Polar Distort, Internal Coefficients\n");
2284 for (i=0; i < 8; i++)
2285 (void) FormatLocaleFile(stderr," c%.20g = %+lf\n",(double) i,
2286 coeff[i]);
2287 (void) FormatLocaleFile(stderr,"Polar Distort, FX Equivalent:\n");
2288 (void) FormatLocaleFile(stderr,"%s", image_gen);
2289 (void) FormatLocaleFile(stderr,
2290 " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",-coeff[2],-coeff[3]);
2291 (void) FormatLocaleFile(stderr," xx=(atan2(ii,jj)%+lf)/(2*pi);\n",
2292 -(coeff[4]+coeff[5])/2 );
2293 (void) FormatLocaleFile(stderr," xx=xx-round(xx);\n");
2294 (void) FormatLocaleFile(stderr," xx=xx*2*pi*%lf + v.w/2;\n",
2295 coeff[6] );
2296 (void) FormatLocaleFile(stderr," yy=(hypot(ii,jj)%+lf)*%lf;\n",
2297 -coeff[1],coeff[7] );
2298 (void) FormatLocaleFile(stderr," v.p{xx-.5,yy-.5}' \\\n");
2299 break;
2300 }
2301 case DePolarDistortion:
2302 {
2303 (void) FormatLocaleFile(stderr,
2304 "DePolar Distort, Internal Coefficients\n");
2305 for (i=0; i < 8; i++)
2306 (void) FormatLocaleFile(stderr," c%.20g = %+lf\n",(double) i,
2307 coeff[i]);
2308 (void) FormatLocaleFile(stderr,"DePolar Distort, FX Equivalent:\n");
2309 (void) FormatLocaleFile(stderr,"%s", image_gen);
2310 (void) FormatLocaleFile(stderr," -fx 'aa=(i+.5)*%lf %+lf;\n",
2311 coeff[6],+coeff[4]);
2312 (void) FormatLocaleFile(stderr," rr=(j+.5)*%lf %+lf;\n",
2313 coeff[7],+coeff[1]);
2314 (void) FormatLocaleFile(stderr," xx=rr*sin(aa) %+lf;\n",
2315 coeff[2]);
2316 (void) FormatLocaleFile(stderr," yy=rr*cos(aa) %+lf;\n",
2317 coeff[3]);
2318 (void) FormatLocaleFile(stderr," v.p{xx-.5,yy-.5}' \\\n");
2319 break;
2320 }
2321 case Cylinder2PlaneDistortion:
2322 {
2323 (void) FormatLocaleFile(stderr,
2324 "Cylinder to Plane Distort, Internal Coefficients\n");
2325 (void) FormatLocaleFile(stderr," cylinder_radius = %+lf\n",coeff[1]);
2326 (void) FormatLocaleFile(stderr,
2327 "Cylinder to Plane Distort, FX Equivalent:\n");
2328 (void) FormatLocaleFile(stderr, "%s", image_gen);
2329 (void) FormatLocaleFile(stderr,
2330 " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",-coeff[4],
2331 -coeff[5]);
2332 (void) FormatLocaleFile(stderr," aa=atan(ii/%+lf);\n",coeff[1]);
2333 (void) FormatLocaleFile(stderr," xx=%lf*aa%+lf;\n",
2334 coeff[1],coeff[2]);
2335 (void) FormatLocaleFile(stderr," yy=jj*cos(aa)%+lf;\n",coeff[3]);
2336 (void) FormatLocaleFile(stderr," %s' \\\n", lookup);
2337 break;
2338 }
2339 case Plane2CylinderDistortion:
2340 {
2341 (void) FormatLocaleFile(stderr,
2342 "Plane to Cylinder Distort, Internal Coefficients\n");
2343 (void) FormatLocaleFile(stderr," cylinder_radius = %+lf\n",coeff[1]);
2344 (void) FormatLocaleFile(stderr,
2345 "Plane to Cylinder Distort, FX Equivalent:\n");
2346 (void) FormatLocaleFile(stderr,"%s", image_gen);
2347 (void) FormatLocaleFile(stderr,
2348 " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",-coeff[4],
2349 -coeff[5]);
2350 (void) FormatLocaleFile(stderr," ii=ii/%+lf;\n",coeff[1]);
2351 (void) FormatLocaleFile(stderr," xx=%lf*tan(ii)%+lf;\n",coeff[1],
2352 coeff[2] );
2353 (void) FormatLocaleFile(stderr," yy=jj/cos(ii)%+lf;\n",coeff[3]);
2354 (void) FormatLocaleFile(stderr," %s' \\\n", lookup);
2355 break;
2356 }
2357 case BarrelDistortion:
2358 case BarrelInverseDistortion:
2359 {
2360 double
2361 xc,
2362 yc;
2363
2364 /*
2365 NOTE: This does the barrel roll in pixel coords not image coords
2366 The internal distortion must do it in image coordinates,
2367 so that is what the center coeff (8,9) is given in.
2368 */
2369 xc=((double)image->columns-1.0)/2.0+image->page.x;
2370 yc=((double)image->rows-1.0)/2.0+image->page.y;
2371 (void) FormatLocaleFile(stderr, "Barrel%s Distort, FX Equivalent:\n",
2372 method == BarrelDistortion ? "" : "Inv");
2373 (void) FormatLocaleFile(stderr, "%s", image_gen);
2374 if ( fabs(coeff[8]-xc-0.5) < 0.1 && fabs(coeff[9]-yc-0.5) < 0.1 )
2375 (void) FormatLocaleFile(stderr," -fx 'xc=(w-1)/2; yc=(h-1)/2;\n");
2376 else
2377 (void) FormatLocaleFile(stderr," -fx 'xc=%lf; yc=%lf;\n",coeff[8]-
2378 0.5,coeff[9]-0.5);
2379 (void) FormatLocaleFile(stderr,
2380 " ii=i-xc; jj=j-yc; rr=hypot(ii,jj);\n");
2381 (void) FormatLocaleFile(stderr,
2382 " ii=ii%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2383 method == BarrelDistortion ? "*" : "/",coeff[0],coeff[1],coeff[2],
2384 coeff[3]);
2385 (void) FormatLocaleFile(stderr,
2386 " jj=jj%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2387 method == BarrelDistortion ? "*" : "/",coeff[4],coeff[5],coeff[6],
2388 coeff[7]);
2389 (void) FormatLocaleFile(stderr," p{ii+xc,jj+yc}' \\\n");
2390 }
2391 default:
2392 break;
2393 }
2394 }
2395 /*
2396 The user provided a 'scale' expert option will scale the output image size,
2397 by the factor given allowing for super-sampling of the distorted image
2398 space. Any scaling factors must naturally be halved as a result.
2399 */
2400 { const char *artifact;
2401 artifact=GetImageArtifact(image,"distort:scale");
2402 output_scaling = 1.0;
2403 if (artifact != (const char *) NULL) {
2404 output_scaling = fabs(StringToDouble(artifact,(char **) NULL));
2405 geometry.width=(size_t) (output_scaling*geometry.width+0.5);
2406 geometry.height=(size_t) (output_scaling*geometry.height+0.5);
2407 geometry.x=(ssize_t) (output_scaling*geometry.x+0.5);
2408 geometry.y=(ssize_t) (output_scaling*geometry.y+0.5);
2409 if ( output_scaling < 0.1 ) {
2410 coeff = (double *) RelinquishMagickMemory(coeff);
2411 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
2412 "InvalidArgument","%s", "-set option:distort:scale" );
2413 return((Image *) NULL);
2414 }
2415 output_scaling = 1/output_scaling;
2416 }
2417 }
2418#define ScaleFilter(F,A,B,C,D) \
2419 ScaleResampleFilter( (F), \
2420 output_scaling*(A), output_scaling*(B), \
2421 output_scaling*(C), output_scaling*(D) )
2422
2423 /*
2424 Initialize the distort image attributes.
2425 */
2426 distort_image=CloneImage(image,geometry.width,geometry.height,MagickTrue,
2427 exception);
2428 if (distort_image == (Image *) NULL)
2429 {
2430 coeff=(double *) RelinquishMagickMemory(coeff);
2431 return((Image *) NULL);
2432 }
2433 /* if image is ColorMapped - change it to DirectClass */
2434 if (SetImageStorageClass(distort_image,DirectClass,exception) == MagickFalse)
2435 {
2436 coeff=(double *) RelinquishMagickMemory(coeff);
2437 distort_image=DestroyImage(distort_image);
2438 return((Image *) NULL);
2439 }
2440 if ((IsPixelInfoGray(&distort_image->background_color) == MagickFalse) &&
2441 (IsGrayColorspace(distort_image->colorspace) != MagickFalse))
2442 (void) SetImageColorspace(distort_image,sRGBColorspace,exception);
2443 if (distort_image->background_color.alpha_trait != UndefinedPixelTrait)
2444 distort_image->alpha_trait=BlendPixelTrait;
2445 distort_image->page.x=geometry.x;
2446 distort_image->page.y=geometry.y;
2447 ConformPixelInfo(distort_image,&distort_image->matte_color,&invalid,
2448 exception);
2449
2450 { /* ----- MAIN CODE -----
2451 Sample the source image to each pixel in the distort image.
2452 */
2453 CacheView
2454 *distort_view;
2455
2456 MagickBooleanType
2457 status;
2458
2459 MagickOffsetType
2460 progress;
2461
2462 PixelInfo
2463 zero;
2464
2465 ResampleFilter
2466 **magick_restrict resample_filter;
2467
2468 ssize_t
2469 j;
2470
2471 status=MagickTrue;
2472 progress=0;
2473 GetPixelInfo(distort_image,&zero);
2474 resample_filter=AcquireResampleFilterTLS(image,UndefinedVirtualPixelMethod,
2475 MagickFalse,exception);
2476 distort_view=AcquireAuthenticCacheView(distort_image,exception);
2477#if defined(MAGICKCORE_OPENMP_SUPPORT)
2478 #pragma omp parallel for schedule(static) shared(progress,status) \
2479 magick_number_threads(image,distort_image,distort_image->rows,1)
2480#endif
2481 for (j=0; j < (ssize_t) distort_image->rows; j++)
2482 {
2483 const int
2484 id = GetOpenMPThreadId();
2485
2486 double
2487 validity; /* how mathematically valid is this the mapping */
2488
2489 MagickBooleanType
2490 sync;
2491
2492 PixelInfo
2493 pixel; /* pixel color to assign to distorted image */
2494
2495 PointInfo
2496 d,
2497 s; /* transform destination image x,y to source image x,y */
2498
2499 ssize_t
2500 i;
2501
2502 Quantum
2503 *magick_restrict q;
2504
2505 q=QueueCacheViewAuthenticPixels(distort_view,0,j,distort_image->columns,1,
2506 exception);
2507 if (q == (Quantum *) NULL)
2508 {
2509 status=MagickFalse;
2510 continue;
2511 }
2512 pixel=zero;
2513
2514 /* Define constant scaling vectors for Affine Distortions
2515 Other methods are either variable, or use interpolated lookup
2516 */
2517 switch (method)
2518 {
2519 case AffineDistortion:
2520 case RigidAffineDistortion:
2521 ScaleFilter( resample_filter[id],
2522 coeff[0], coeff[1],
2523 coeff[3], coeff[4] );
2524 break;
2525 default:
2526 break;
2527 }
2528
2529 /* Initialize default pixel validity
2530 * negative: pixel is invalid output 'matte_color'
2531 * 0.0 to 1.0: antialiased, mix with resample output
2532 * 1.0 or greater: use resampled output.
2533 */
2534 validity = 1.0;
2535
2536 for (i=0; i < (ssize_t) distort_image->columns; i++)
2537 {
2538 /* map pixel coordinate to distortion space coordinate */
2539 d.x = (double) (geometry.x+i+0.5)*output_scaling;
2540 d.y = (double) (geometry.y+j+0.5)*output_scaling;
2541 s = d; /* default is a no-op mapping */
2542 switch (method)
2543 {
2544 case AffineDistortion:
2545 case RigidAffineDistortion:
2546 {
2547 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2548 s.y=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2549 /* Affine partial derivatives are constant -- set above */
2550 break;
2551 }
2552 case PerspectiveDistortion:
2553 {
2554 double
2555 p,n,r,abs_r,abs_c6,abs_c7,scale;
2556 /* perspective is a ratio of affines */
2557 p=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2558 n=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2559 r=coeff[6]*d.x+coeff[7]*d.y+1.0;
2560 /* Pixel Validity -- is it a 'sky' or 'ground' pixel */
2561 validity = (r*coeff[8] < 0.0) ? 0.0 : 1.0;
2562 /* Determine horizon anti-alias blending */
2563 abs_r = fabs(r)*2;
2564 abs_c6 = fabs(coeff[6]);
2565 abs_c7 = fabs(coeff[7]);
2566 if ( abs_c6 > abs_c7 ) {
2567 if ( abs_r < abs_c6*output_scaling )
2568 validity = 0.5 - coeff[8]*r/(coeff[6]*output_scaling);
2569 }
2570 else if ( abs_r < abs_c7*output_scaling )
2571 validity = 0.5 - coeff[8]*r/(coeff[7]*output_scaling);
2572 /* Perspective Sampling Point (if valid) */
2573 if ( validity > 0.0 ) {
2574 /* divide by r affine, for perspective scaling */
2575 scale = 1.0/r;
2576 s.x = p*scale;
2577 s.y = n*scale;
2578 /* Perspective Partial Derivatives or Scaling Vectors */
2579 scale *= scale;
2580 ScaleFilter( resample_filter[id],
2581 (r*coeff[0] - p*coeff[6])*scale,
2582 (r*coeff[1] - p*coeff[7])*scale,
2583 (r*coeff[3] - n*coeff[6])*scale,
2584 (r*coeff[4] - n*coeff[7])*scale );
2585 }
2586 break;
2587 }
2588 case BilinearReverseDistortion:
2589 {
2590 /* Reversed Mapped is just a simple polynomial */
2591 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2]*d.x*d.y+coeff[3];
2592 s.y=coeff[4]*d.x+coeff[5]*d.y
2593 +coeff[6]*d.x*d.y+coeff[7];
2594 /* Bilinear partial derivatives of scaling vectors */
2595 ScaleFilter( resample_filter[id],
2596 coeff[0] + coeff[2]*d.y,
2597 coeff[1] + coeff[2]*d.x,
2598 coeff[4] + coeff[6]*d.y,
2599 coeff[5] + coeff[6]*d.x );
2600 break;
2601 }
2602 case BilinearForwardDistortion:
2603 {
2604 /* Forward mapped needs reversed polynomial equations
2605 * which unfortunately requires a square root! */
2606 double b,c;
2607 d.x -= coeff[3]; d.y -= coeff[7];
2608 b = coeff[6]*d.x - coeff[2]*d.y + coeff[8];
2609 c = coeff[4]*d.x - coeff[0]*d.y;
2610
2611 validity = 1.0;
2612 /* Handle Special degenerate (non-quadratic) case
2613 * Currently without horizon anti-aliasing */
2614 if ( fabs(coeff[9]) < MagickEpsilon )
2615 s.y = -c/b;
2616 else {
2617 c = b*b - 2*coeff[9]*c;
2618 if ( c < 0.0 )
2619 validity = 0.0;
2620 else
2621 s.y = ( -b + sqrt(c) )/coeff[9];
2622 }
2623 if ( validity > 0.0 )
2624 s.x = ( d.x - coeff[1]*s.y) / ( coeff[0] + coeff[2]*s.y );
2625
2626 /* NOTE: the sign of the square root should be -ve for parts
2627 where the source image becomes 'flipped' or 'mirrored'.
2628 FUTURE: Horizon handling
2629 FUTURE: Scaling factors or Derivatives (how?)
2630 */
2631 break;
2632 }
2633#if 0
2634 case BilinearDistortion:
2635 /* Bilinear mapping of any Quadrilateral to any Quadrilateral */
2636 /* UNDER DEVELOPMENT */
2637 break;
2638#endif
2639 case PolynomialDistortion:
2640 {
2641 /* multi-ordered polynomial */
2642 ssize_t
2643 k;
2644
2645 ssize_t
2646 nterms=(ssize_t)coeff[1];
2647
2648 PointInfo
2649 du,dv; /* the du,dv vectors from unit dx,dy -- derivatives */
2650
2651 s.x=s.y=du.x=du.y=dv.x=dv.y=0.0;
2652 for(k=0; k < nterms; k++) {
2653 s.x += poly_basis_fn(k,d.x,d.y)*coeff[2+k];
2654 du.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k];
2655 du.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k];
2656 s.y += poly_basis_fn(k,d.x,d.y)*coeff[2+k+nterms];
2657 dv.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k+nterms];
2658 dv.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k+nterms];
2659 }
2660 ScaleFilter( resample_filter[id], du.x,du.y,dv.x,dv.y );
2661 break;
2662 }
2663 case ArcDistortion:
2664 {
2665 /* what is the angle and radius in the destination image */
2666 s.x = (double) ((atan2(d.y,d.x) - coeff[0])/Magick2PI);
2667 s.x -= MagickRound(s.x); /* angle */
2668 s.y = hypot(d.x,d.y); /* radius */
2669
2670 /* Arc Distortion Partial Scaling Vectors
2671 Are derived by mapping the perpendicular unit vectors
2672 dR and dA*R*2PI rather than trying to map dx and dy
2673 The results is a very simple orthogonal aligned ellipse.
2674 */
2675 if ( s.y > MagickEpsilon )
2676 ScaleFilter( resample_filter[id],
2677 (double) (coeff[1]/(Magick2PI*s.y)), 0, 0, coeff[3] );
2678 else
2679 ScaleFilter( resample_filter[id],
2680 distort_image->columns*2, 0, 0, coeff[3] );
2681
2682 /* now scale the angle and radius for source image lookup point */
2683 s.x = s.x*coeff[1] + coeff[4] + image->page.x +0.5;
2684 s.y = (coeff[2] - s.y) * coeff[3] + image->page.y;
2685 break;
2686 }
2687 case PolarDistortion:
2688 { /* 2D Cartesian to Polar View */
2689 d.x -= coeff[2];
2690 d.y -= coeff[3];
2691 s.x = atan2(d.x,d.y) - (coeff[4]+coeff[5])/2;
2692 s.x /= Magick2PI;
2693 s.x -= MagickRound(s.x);
2694 s.x *= Magick2PI; /* angle - relative to centerline */
2695 s.y = hypot(d.x,d.y); /* radius */
2696
2697 /* Polar Scaling vectors are based on mapping dR and dA vectors
2698 This results in very simple orthogonal scaling vectors
2699 */
2700 if ( s.y > MagickEpsilon )
2701 ScaleFilter( resample_filter[id],
2702 (double) (coeff[6]/(Magick2PI*s.y)), 0, 0, coeff[7] );
2703 else
2704 ScaleFilter( resample_filter[id],
2705 distort_image->columns*2, 0, 0, coeff[7] );
2706
2707 /* now finish mapping radius/angle to source x,y coords */
2708 s.x = s.x*coeff[6] + (double)image->columns/2.0 + image->page.x;
2709 s.y = (s.y-coeff[1])*coeff[7] + image->page.y;
2710 break;
2711 }
2712 case DePolarDistortion:
2713 { /* @D Polar to Cartesian */
2714 /* ignore all destination virtual offsets */
2715 d.x = ((double)i+0.5)*output_scaling*coeff[6]+coeff[4];
2716 d.y = ((double)j+0.5)*output_scaling*coeff[7]+coeff[1];
2717 s.x = d.y*sin(d.x) + coeff[2];
2718 s.y = d.y*cos(d.x) + coeff[3];
2719 /* derivatives are useless - better to use SuperSampling */
2720 break;
2721 }
2722 case Cylinder2PlaneDistortion:
2723 { /* 3D Cylinder to Tangential Plane */
2724 double ax, cx;
2725 /* relative to center of distortion */
2726 d.x -= coeff[4]; d.y -= coeff[5];
2727 d.x /= coeff[1]; /* x' = x/r */
2728 ax=atan(d.x); /* aa = atan(x/r) = u/r */
2729 cx=cos(ax); /* cx = cos(atan(x/r)) = 1/sqrt(x^2+u^2) */
2730 s.x = coeff[1]*ax; /* u = r*atan(x/r) */
2731 s.y = d.y*cx; /* v = y*cos(u/r) */
2732 /* derivatives... (see personal notes) */
2733 ScaleFilter( resample_filter[id],
2734 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2735#if 0
2736if ( i == 0 && j == 0 ) {
2737 fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2738 fprintf(stderr, "phi = %lf\n", (double)(ax * 180.0/MagickPI) );
2739 fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n",
2740 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2741 fflush(stderr); }
2742#endif
2743 /* add center of distortion in source */
2744 s.x += coeff[2]; s.y += coeff[3];
2745 break;
2746 }
2747 case Plane2CylinderDistortion:
2748 { /* 3D Cylinder to Tangential Plane */
2749 /* relative to center of distortion */
2750 d.x -= coeff[4]; d.y -= coeff[5];
2751
2752 /* is pixel valid - horizon of a infinite Virtual-Pixel Plane
2753 * (see Anthony Thyssen's personal note) */
2754 validity = (double) (coeff[1]*MagickPI2 - fabs(d.x))/output_scaling + 0.5;
2755
2756 if ( validity > 0.0 ) {
2757 double cx,tx;
2758 d.x /= coeff[1]; /* x'= x/r */
2759 cx = 1/cos(d.x); /* cx = 1/cos(x/r) */
2760 tx = tan(d.x); /* tx = tan(x/r) */
2761 s.x = coeff[1]*tx; /* u = r * tan(x/r) */
2762 s.y = d.y*cx; /* v = y / cos(x/r) */
2763 /* derivatives... (see Anthony Thyssen's personal notes) */
2764 ScaleFilter( resample_filter[id],
2765 cx*cx, 0.0, s.y*cx/coeff[1], cx );
2766#if 0
2767/*if ( i == 0 && j == 0 )*/
2768if ( d.x == 0.5 && d.y == 0.5 ) {
2769 fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2770 fprintf(stderr, "radius = %lf phi = %lf validity = %lf\n",
2771 coeff[1], (double)(d.x * 180.0/MagickPI), validity );
2772 fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n",
2773 cx*cx, 0.0, s.y*cx/coeff[1], cx);
2774 fflush(stderr); }
2775#endif
2776 }
2777 /* add center of distortion in source */
2778 s.x += coeff[2]; s.y += coeff[3];
2779 break;
2780 }
2781 case BarrelDistortion:
2782 case BarrelInverseDistortion:
2783 { /* Lens Barrel Distortion Correction */
2784 double r,fx,fy,gx,gy;
2785 /* Radial Polynomial Distortion (de-normalized) */
2786 d.x -= coeff[8];
2787 d.y -= coeff[9];
2788 r = sqrt(d.x*d.x+d.y*d.y);
2789 if ( r > MagickEpsilon ) {
2790 fx = ((coeff[0]*r + coeff[1])*r + coeff[2])*r + coeff[3];
2791 fy = ((coeff[4]*r + coeff[5])*r + coeff[6])*r + coeff[7];
2792 gx = ((3*coeff[0]*r + 2*coeff[1])*r + coeff[2])/r;
2793 gy = ((3*coeff[4]*r + 2*coeff[5])*r + coeff[6])/r;
2794 /* adjust functions and scaling for 'inverse' form */
2795 if ( method == BarrelInverseDistortion ) {
2796 fx = 1/fx; fy = 1/fy;
2797 gx *= -fx*fx; gy *= -fy*fy;
2798 }
2799 /* Set the source pixel to lookup and EWA derivative vectors */
2800 s.x = d.x*fx + coeff[8];
2801 s.y = d.y*fy + coeff[9];
2802 ScaleFilter( resample_filter[id],
2803 gx*d.x*d.x + fx, gx*d.x*d.y,
2804 gy*d.x*d.y, gy*d.y*d.y + fy );
2805 }
2806 else {
2807 /* Special handling to avoid divide by zero when r==0
2808 **
2809 ** The source and destination pixels match in this case
2810 ** which was set at the top of the loop using s = d;
2811 ** otherwise... s.x=coeff[8]; s.y=coeff[9];
2812 */
2813 if ( method == BarrelDistortion )
2814 ScaleFilter( resample_filter[id],
2815 coeff[3], 0, 0, coeff[7] );
2816 else /* method == BarrelInverseDistortion */
2817 /* FUTURE, trap for D==0 causing division by zero */
2818 ScaleFilter( resample_filter[id],
2819 1.0/coeff[3], 0, 0, 1.0/coeff[7] );
2820 }
2821 break;
2822 }
2823 case ShepardsDistortion:
2824 { /* Shepards Method, or Inverse Weighted Distance for
2825 displacement around the destination image control points
2826 The input arguments are the coefficients to the function.
2827 This is more of a 'displacement' function rather than an
2828 absolute distortion function.
2829
2830 Note: We can not determine derivatives using shepards method
2831 so only a point sample interpolation can be used.
2832 */
2833 double
2834 denominator;
2835
2836 size_t
2837 k;
2838
2839 denominator = s.x = s.y = 0;
2840 for(k=0; k<number_arguments; k+=4) {
2841 double weight =
2842 ((double)d.x-arguments[k+2])*((double)d.x-arguments[k+2])
2843 + ((double)d.y-arguments[k+3])*((double)d.y-arguments[k+3]);
2844 weight = pow(weight,coeff[0]); /* shepards power factor */
2845 weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
2846
2847 s.x += (arguments[ k ]-arguments[k+2])*weight;
2848 s.y += (arguments[k+1]-arguments[k+3])*weight;
2849 denominator += weight;
2850 }
2851 s.x /= denominator;
2852 s.y /= denominator;
2853 s.x += d.x; /* make it as relative displacement */
2854 s.y += d.y;
2855 break;
2856 }
2857 default:
2858 break; /* use the default no-op given above */
2859 }
2860 /* map virtual canvas location back to real image coordinate */
2861 if ( bestfit && method != ArcDistortion ) {
2862 s.x -= image->page.x;
2863 s.y -= image->page.y;
2864 }
2865 s.x -= 0.5;
2866 s.y -= 0.5;
2867
2868 if ( validity <= 0.0 ) {
2869 /* result of distortion is an invalid pixel - don't resample */
2870 SetPixelViaPixelInfo(distort_image,&invalid,q);
2871 }
2872 else {
2873 /* resample the source image to find its correct color */
2874 status=ResamplePixelColor(resample_filter[id],s.x,s.y,&pixel,
2875 exception);
2876 if (status == MagickFalse)
2877 SetPixelViaPixelInfo(distort_image,&invalid,q);
2878 else
2879 {
2880 /* if validity between 0.0 & 1.0 mix result with invalid pixel */
2881 if ( validity < 1.0 ) {
2882 /* Do a blend of sample color and invalid pixel */
2883 /* should this be a 'Blend', or an 'Over' compose */
2884 CompositePixelInfoBlend(&pixel,validity,&invalid,(1.0-validity),
2885 &pixel);
2886 }
2887 SetPixelViaPixelInfo(distort_image,&pixel,q);
2888 }
2889 }
2890 q+=(ptrdiff_t) GetPixelChannels(distort_image);
2891 }
2892 sync=SyncCacheViewAuthenticPixels(distort_view,exception);
2893 if (sync == MagickFalse)
2894 status=MagickFalse;
2895 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2896 {
2897 MagickBooleanType
2898 proceed;
2899
2900#if defined(MAGICKCORE_OPENMP_SUPPORT)
2901 #pragma omp atomic
2902#endif
2903 progress++;
2904 proceed=SetImageProgress(image,DistortImageTag,progress,image->rows);
2905 if (proceed == MagickFalse)
2906 status=MagickFalse;
2907 }
2908 }
2909 distort_view=DestroyCacheView(distort_view);
2910 resample_filter=DestroyResampleFilterTLS(resample_filter);
2911
2912 if (status == MagickFalse)
2913 distort_image=DestroyImage(distort_image);
2914 }
2915
2916 /* Arc does not return an offset unless 'bestfit' is in effect
2917 And the user has not provided an overriding 'viewport'.
2918 */
2919 if ( method == ArcDistortion && !bestfit && !viewport_given ) {
2920 distort_image->page.x = 0;
2921 distort_image->page.y = 0;
2922 }
2923 coeff=(double *) RelinquishMagickMemory(coeff);
2924 return(distort_image);
2925}
2926␌
2927/*
2928%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2929% %
2930% %
2931% %
2932% R o t a t e I m a g e %
2933% %
2934% %
2935% %
2936%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2937%
2938% RotateImage() creates a new image that is a rotated copy of an existing
2939% one. Positive angles rotate counter-clockwise (right-hand rule), while
2940% negative angles rotate clockwise. Rotated images are usually larger than
2941% the originals and have 'empty' triangular corners. X axis. Empty
2942% triangles left over from shearing the image are filled with the background
2943% color defined by member 'background_color' of the image. RotateImage
2944% allocates the memory necessary for the new Image structure and returns a
2945% pointer to the new image.
2946%
2947% The format of the RotateImage method is:
2948%
2949% Image *RotateImage(const Image *image,const double degrees,
2950% ExceptionInfo *exception)
2951%
2952% A description of each parameter follows.
2953%
2954% o image: the image.
2955%
2956% o degrees: Specifies the number of degrees to rotate the image.
2957%
2958% o exception: return any errors or warnings in this structure.
2959%
2960*/
2961MagickExport Image *RotateImage(const Image *image,const double degrees,
2962 ExceptionInfo *exception)
2963{
2964 Image
2965 *distort_image,
2966 *rotate_image;
2967
2968 double
2969 angle;
2970
2971 PointInfo
2972 shear;
2973
2974 size_t
2975 rotations;
2976
2977 /*
2978 Adjust rotation angle.
2979 */
2980 assert(image != (Image *) NULL);
2981 assert(image->signature == MagickCoreSignature);
2982 if (IsEventLogging() != MagickFalse)
2983 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2984 assert(exception != (ExceptionInfo *) NULL);
2985 assert(exception->signature == MagickCoreSignature);
2986 angle=fmod(degrees,360.0);
2987 while (angle < -45.0)
2988 angle+=360.0;
2989 for (rotations=0; angle > 45.0; rotations++)
2990 angle-=90.0;
2991 rotations%=4;
2992 shear.x=(-tan((double) DegreesToRadians(angle)/2.0));
2993 shear.y=sin((double) DegreesToRadians(angle));
2994 if ((fabs(shear.x) < MagickEpsilon) && (fabs(shear.y) < MagickEpsilon))
2995 return(IntegralRotateImage(image,rotations,exception));
2996 distort_image=CloneImage(image,0,0,MagickTrue,exception);
2997 if (distort_image == (Image *) NULL)
2998 return((Image *) NULL);
2999 (void) SetImageVirtualPixelMethod(distort_image,BackgroundVirtualPixelMethod,
3000 exception);
3001 rotate_image=DistortImage(distort_image,ScaleRotateTranslateDistortion,1,
3002 &degrees,MagickTrue,exception);
3003 distort_image=DestroyImage(distort_image);
3004 return(rotate_image);
3005}
3006␌
3007/*
3008%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3009% %
3010% %
3011% %
3012% S p a r s e C o l o r I m a g e %
3013% %
3014% %
3015% %
3016%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3017%
3018% SparseColorImage(), given a set of coordinates, interpolates the colors
3019% found at those coordinates, across the whole image, using various methods.
3020%
3021% The format of the SparseColorImage() method is:
3022%
3023% Image *SparseColorImage(const Image *image,
3024% const SparseColorMethod method,const size_t number_arguments,
3025% const double *arguments,ExceptionInfo *exception)
3026%
3027% A description of each parameter follows:
3028%
3029% o image: the image to be filled in.
3030%
3031% o method: the method to fill in the gradient between the control points.
3032%
3033% The methods used for SparseColor() are often simular to methods
3034% used for DistortImage(), and even share the same code for determination
3035% of the function coefficients, though with more dimensions (or resulting
3036% values).
3037%
3038% o number_arguments: the number of arguments given.
3039%
3040% o arguments: array of floating point arguments for this method--
3041% x,y,color_values-- with color_values given as normalized values.
3042%
3043% o exception: return any errors or warnings in this structure
3044%
3045*/
3046MagickExport Image *SparseColorImage(const Image *image,
3047 const SparseColorMethod method,const size_t number_arguments,
3048 const double *arguments,ExceptionInfo *exception)
3049{
3050#define SparseColorTag "Distort/SparseColor"
3051
3052 double
3053 *coeff;
3054
3055 Image
3056 *sparse_image;
3057
3058 size_t
3059 number_colors;
3060
3061 SparseColorMethod
3062 sparse_method;
3063
3064 assert(image != (Image *) NULL);
3065 assert(image->signature == MagickCoreSignature);
3066 assert(exception != (ExceptionInfo *) NULL);
3067 assert(exception->signature == MagickCoreSignature);
3068 if (IsEventLogging() != MagickFalse)
3069 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
3070
3071 /* Determine number of color values needed per control point */
3072 number_colors=0;
3073 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3074 number_colors++;
3075 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3076 number_colors++;
3077 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3078 number_colors++;
3079 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3080 (image->colorspace == CMYKColorspace))
3081 number_colors++;
3082 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3083 (image->alpha_trait != UndefinedPixelTrait))
3084 number_colors++;
3085
3086 /*
3087 Convert input arguments into mapping coefficients, in this case
3088 we are mapping (distorting) colors, rather than coordinates.
3089 */
3090 { DistortMethod
3091 distort_method;
3092
3093 distort_method=(DistortMethod) method;
3094 if ( distort_method >= SentinelDistortion )
3095 distort_method = ShepardsDistortion; /* Pretend to be Shepards */
3096 coeff = GenerateCoefficients(image, &distort_method, number_arguments,
3097 arguments, number_colors, exception);
3098 if ( coeff == (double *) NULL )
3099 return((Image *) NULL);
3100 /*
3101 Note some Distort Methods may fall back to other simpler methods,
3102 Currently the only fallback of concern is Bilinear to Affine
3103 (Barycentric), which is also sparse_colr method. This also ensures
3104 correct two and one color Barycentric handling.
3105 */
3106 sparse_method = (SparseColorMethod) distort_method;
3107 if ( distort_method == ShepardsDistortion )
3108 sparse_method = method; /* return non-distort methods to normal */
3109 if ( sparse_method == InverseColorInterpolate )
3110 coeff[0]=0.5; /* sqrt() the squared distance for inverse */
3111 }
3112
3113 /* Verbose output */
3114 if (IsStringTrue(GetImageArtifact(image,"verbose")) != MagickFalse) {
3115
3116 switch (sparse_method) {
3117 case BarycentricColorInterpolate:
3118 {
3119 ssize_t x=0;
3120 (void) FormatLocaleFile(stderr, "Barycentric Sparse Color:\n");
3121 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3122 (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf' \\\n",
3123 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
3124 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3125 (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf' \\\n",
3126 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
3127 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3128 (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf' \\\n",
3129 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
3130 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3131 (image->colorspace == CMYKColorspace))
3132 (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf' \\\n",
3133 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
3134 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3135 (image->alpha_trait != UndefinedPixelTrait))
3136 (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf' \\\n",
3137 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
3138 break;
3139 }
3140 case BilinearColorInterpolate:
3141 {
3142 ssize_t x=0;
3143 (void) FormatLocaleFile(stderr, "Bilinear Sparse Color\n");
3144 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3145 (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3146 coeff[ x ], coeff[x+1],
3147 coeff[x+2], coeff[x+3]),x+=4;
3148 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3149 (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3150 coeff[ x ], coeff[x+1],
3151 coeff[x+2], coeff[x+3]),x+=4;
3152 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3153 (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3154 coeff[ x ], coeff[x+1],
3155 coeff[x+2], coeff[x+3]),x+=4;
3156 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3157 (image->colorspace == CMYKColorspace))
3158 (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3159 coeff[ x ], coeff[x+1],
3160 coeff[x+2], coeff[x+3]),x+=4;
3161 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3162 (image->alpha_trait != UndefinedPixelTrait))
3163 (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3164 coeff[ x ], coeff[x+1],
3165 coeff[x+2], coeff[x+3]),x+=4;
3166 break;
3167 }
3168 default:
3169 /* sparse color method is too complex for FX emulation */
3170 break;
3171 }
3172 }
3173
3174 /* Generate new image for generated interpolated gradient.
3175 * ASIDE: Actually we could have just replaced the colors of the original
3176 * image, but IM Core policy, is if storage class could change then clone
3177 * the image.
3178 */
3179
3180 sparse_image=CloneImage(image,0,0,MagickTrue,exception);
3181 if (sparse_image == (Image *) NULL)
3182 return((Image *) NULL);
3183 if (SetImageStorageClass(sparse_image,DirectClass,exception) == MagickFalse)
3184 { /* if image is ColorMapped - change it to DirectClass */
3185 sparse_image=DestroyImage(sparse_image);
3186 return((Image *) NULL);
3187 }
3188 if (IsGrayColorspace(sparse_image->colorspace) != MagickFalse)
3189 (void) SetImageColorspace(sparse_image,sRGBColorspace,exception);
3190 { /* ----- MAIN CODE ----- */
3191 CacheView
3192 *sparse_view;
3193
3194 MagickBooleanType
3195 status;
3196
3197 MagickOffsetType
3198 progress;
3199
3200 ssize_t
3201 j;
3202
3203 status=MagickTrue;
3204 progress=0;
3205 sparse_view=AcquireAuthenticCacheView(sparse_image,exception);
3206#if defined(MAGICKCORE_OPENMP_SUPPORT)
3207 #pragma omp parallel for schedule(static) shared(progress,status) \
3208 magick_number_threads(image,sparse_image,sparse_image->rows,1)
3209#endif
3210 for (j=0; j < (ssize_t) sparse_image->rows; j++)
3211 {
3212 MagickBooleanType
3213 sync;
3214
3215 PixelInfo
3216 pixel; /* pixel to assign to distorted image */
3217
3218 Quantum
3219 *magick_restrict q;
3220
3221 ssize_t
3222 i;
3223
3224 q=GetCacheViewAuthenticPixels(sparse_view,0,j,sparse_image->columns,1,
3225 exception);
3226 if (q == (Quantum *) NULL)
3227 {
3228 status=MagickFalse;
3229 continue;
3230 }
3231 GetPixelInfo(sparse_image,&pixel);
3232 for (i=0; i < (ssize_t) sparse_image->columns; i++)
3233 {
3234 GetPixelInfoPixel(sparse_image,q,&pixel);
3235 switch (sparse_method)
3236 {
3237 case BarycentricColorInterpolate:
3238 {
3239 ssize_t x=0;
3240 if ((GetPixelRedTraits(sparse_image) & UpdatePixelTrait) != 0)
3241 pixel.red = coeff[x]*i +coeff[x+1]*j
3242 +coeff[x+2], x+=3;
3243 if ((GetPixelGreenTraits(sparse_image) & UpdatePixelTrait) != 0)
3244 pixel.green = coeff[x]*i +coeff[x+1]*j
3245 +coeff[x+2], x+=3;
3246 if ((GetPixelBlueTraits(sparse_image) & UpdatePixelTrait) != 0)
3247 pixel.blue = coeff[x]*i +coeff[x+1]*j
3248 +coeff[x+2], x+=3;
3249 if (((GetPixelBlackTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3250 (sparse_image->colorspace == CMYKColorspace))
3251 pixel.black = coeff[x]*i +coeff[x+1]*j
3252 +coeff[x+2], x+=3;
3253 if (((GetPixelAlphaTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3254 (sparse_image->alpha_trait != UndefinedPixelTrait))
3255 pixel.alpha = coeff[x]*i +coeff[x+1]*j
3256 +coeff[x+2], x+=3;
3257 break;
3258 }
3259 case BilinearColorInterpolate:
3260 {
3261 ssize_t x=0;
3262 if ((GetPixelRedTraits(sparse_image) & UpdatePixelTrait) != 0)
3263 pixel.red = coeff[x]*i + coeff[x+1]*j +
3264 coeff[x+2]*i*j + coeff[x+3], x+=4;
3265 if ((GetPixelGreenTraits(sparse_image) & UpdatePixelTrait) != 0)
3266 pixel.green = coeff[x]*i + coeff[x+1]*j +
3267 coeff[x+2]*i*j + coeff[x+3], x+=4;
3268 if ((GetPixelBlueTraits(sparse_image) & UpdatePixelTrait) != 0)
3269 pixel.blue = coeff[x]*i + coeff[x+1]*j +
3270 coeff[x+2]*i*j + coeff[x+3], x+=4;
3271 if (((GetPixelBlackTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3272 (image->colorspace == CMYKColorspace))
3273 pixel.black = coeff[x]*i + coeff[x+1]*j +
3274 coeff[x+2]*i*j + coeff[x+3], x+=4;
3275 if (((GetPixelAlphaTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3276 (sparse_image->alpha_trait != UndefinedPixelTrait))
3277 pixel.alpha = coeff[x]*i + coeff[x+1]*j +
3278 coeff[x+2]*i*j + coeff[x+3], x+=4;
3279 break;
3280 }
3281 case InverseColorInterpolate:
3282 case ShepardsColorInterpolate:
3283 { /* Inverse (Squared) Distance weights average (IDW) */
3284 double
3285 denominator;
3286
3287 size_t
3288 k;
3289
3290 if ((GetPixelRedTraits(sparse_image) & UpdatePixelTrait) != 0)
3291 pixel.red=0.0;
3292 if ((GetPixelGreenTraits(sparse_image) & UpdatePixelTrait) != 0)
3293 pixel.green=0.0;
3294 if ((GetPixelBlueTraits(sparse_image) & UpdatePixelTrait) != 0)
3295 pixel.blue=0.0;
3296 if (((GetPixelBlackTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3297 (image->colorspace == CMYKColorspace))
3298 pixel.black=0.0;
3299 if (((GetPixelAlphaTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3300 (sparse_image->alpha_trait != UndefinedPixelTrait))
3301 pixel.alpha=0.0;
3302 denominator = 0.0;
3303 for (k=0; k<number_arguments; k+=2+number_colors)
3304 {
3305 double weight =
3306 ((double) i-arguments[ k ])*((double) i-arguments[ k ])
3307 + ((double) j-arguments[k+1])*((double) j-arguments[k+1]);
3308 ssize_t x = (ssize_t) k+2;
3309
3310 weight = pow(weight,coeff[0]); /* inverse of power factor */
3311 weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
3312 if ((GetPixelRedTraits(sparse_image) & UpdatePixelTrait) != 0)
3313 pixel.red += arguments[x++]*weight;
3314 if ((GetPixelGreenTraits(sparse_image) & UpdatePixelTrait) != 0)
3315 pixel.green += arguments[x++]*weight;
3316 if ((GetPixelBlueTraits(sparse_image) & UpdatePixelTrait) != 0)
3317 pixel.blue += arguments[x++]*weight;
3318 if (((GetPixelBlackTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3319 (image->colorspace == CMYKColorspace))
3320 pixel.black += arguments[x++]*weight;
3321 if (((GetPixelAlphaTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3322 (sparse_image->alpha_trait != UndefinedPixelTrait))
3323 pixel.alpha += arguments[x++]*weight;
3324 denominator += weight;
3325 }
3326 if ((GetPixelRedTraits(sparse_image) & UpdatePixelTrait) != 0)
3327 pixel.red/=denominator;
3328 if ((GetPixelGreenTraits(sparse_image) & UpdatePixelTrait) != 0)
3329 pixel.green/=denominator;
3330 if ((GetPixelBlueTraits(sparse_image) & UpdatePixelTrait) != 0)
3331 pixel.blue/=denominator;
3332 if (((GetPixelBlackTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3333 (image->colorspace == CMYKColorspace))
3334 pixel.black/=denominator;
3335 if (((GetPixelAlphaTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3336 (sparse_image->alpha_trait != UndefinedPixelTrait))
3337 pixel.alpha/=denominator;
3338 break;
3339 }
3340 case ManhattanColorInterpolate:
3341 {
3342 double
3343 minimum = MagickMaximumValue;
3344
3345 size_t
3346 k;
3347
3348 /*
3349 Just use the closest control point you can find!
3350 */
3351 for (k=0; k<number_arguments; k+=2+number_colors)
3352 {
3353 double distance = fabs((double)i-arguments[ k ])+
3354 fabs((double)j-arguments[k+1]);
3355 if ( distance < minimum ) {
3356 ssize_t x=(ssize_t) k+2;
3357 if ((GetPixelRedTraits(sparse_image) & UpdatePixelTrait) != 0)
3358 pixel.red=arguments[x++];
3359 if ((GetPixelGreenTraits(sparse_image) & UpdatePixelTrait) != 0)
3360 pixel.green=arguments[x++];
3361 if ((GetPixelBlueTraits(sparse_image) & UpdatePixelTrait) != 0)
3362 pixel.blue=arguments[x++];
3363 if (((GetPixelBlackTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3364 (image->colorspace == CMYKColorspace))
3365 pixel.black=arguments[x++];
3366 if (((GetPixelAlphaTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3367 (sparse_image->alpha_trait != UndefinedPixelTrait))
3368 pixel.alpha=arguments[x++];
3369 minimum = distance;
3370 }
3371 }
3372 break;
3373 }
3374 case VoronoiColorInterpolate:
3375 default:
3376 {
3377 double
3378 minimum = MagickMaximumValue;
3379
3380 size_t
3381 k;
3382
3383 /*
3384 Just use the closest control point you can find!
3385 */
3386 for (k=0; k<number_arguments; k+=2+number_colors) {
3387 double distance =
3388 ((double) i-arguments[ k ])*((double) i-arguments[ k ])
3389 + ((double) j-arguments[k+1])*((double) j-arguments[k+1]);
3390 if ( distance < minimum ) {
3391 ssize_t x = (ssize_t) k+2;
3392 if ((GetPixelRedTraits(sparse_image) & UpdatePixelTrait) != 0)
3393 pixel.red=arguments[x++];
3394 if ((GetPixelGreenTraits(sparse_image) & UpdatePixelTrait) != 0)
3395 pixel.green=arguments[x++];
3396 if ((GetPixelBlueTraits(sparse_image) & UpdatePixelTrait) != 0)
3397 pixel.blue=arguments[x++];
3398 if (((GetPixelBlackTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3399 (image->colorspace == CMYKColorspace))
3400 pixel.black=arguments[x++];
3401 if (((GetPixelAlphaTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3402 (sparse_image->alpha_trait != UndefinedPixelTrait))
3403 pixel.alpha=arguments[x++];
3404 minimum = distance;
3405 }
3406 }
3407 break;
3408 }
3409 }
3410 /* set the color directly back into the source image */
3411 if ((GetPixelRedTraits(sparse_image) & UpdatePixelTrait) != 0)
3412 pixel.red=(MagickRealType) ClampPixel((double) QuantumRange*
3413 pixel.red);
3414 if ((GetPixelGreenTraits(sparse_image) & UpdatePixelTrait) != 0)
3415 pixel.green=(MagickRealType) ClampPixel((double) QuantumRange*
3416 pixel.green);
3417 if ((GetPixelBlueTraits(sparse_image) & UpdatePixelTrait) != 0)
3418 pixel.blue=(MagickRealType) ClampPixel((double) QuantumRange*
3419 pixel.blue);
3420 if (((GetPixelBlackTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3421 (image->colorspace == CMYKColorspace))
3422 pixel.black=(MagickRealType) ClampPixel((double) QuantumRange*
3423 pixel.black);
3424 if (((GetPixelAlphaTraits(sparse_image) & UpdatePixelTrait) != 0) &&
3425 (image->alpha_trait != UndefinedPixelTrait))
3426 pixel.alpha=(MagickRealType) ClampPixel((double) QuantumRange*
3427 pixel.alpha);
3428 SetPixelViaPixelInfo(sparse_image,&pixel,q);
3429 q+=(ptrdiff_t) GetPixelChannels(sparse_image);
3430 }
3431 sync=SyncCacheViewAuthenticPixels(sparse_view,exception);
3432 if (sync == MagickFalse)
3433 status=MagickFalse;
3434 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3435 {
3436 MagickBooleanType
3437 proceed;
3438
3439#if defined(MAGICKCORE_OPENMP_SUPPORT)
3440 #pragma omp atomic
3441#endif
3442 progress++;
3443 proceed=SetImageProgress(image,SparseColorTag,progress,image->rows);
3444 if (proceed == MagickFalse)
3445 status=MagickFalse;
3446 }
3447 }
3448 sparse_view=DestroyCacheView(sparse_view);
3449 if (status == MagickFalse)
3450 sparse_image=DestroyImage(sparse_image);
3451 }
3452 coeff = (double *) RelinquishMagickMemory(coeff);
3453 return(sparse_image);
3454}
_AffineMatrix
Definition geometry.h:95
_CacheView
Definition cache-view.c:66
_ExceptionInfo
Definition exception.h:102
_Image
Definition image.h:132
_PixelInfo
Definition pixel.h:182
_PointInfo
Definition geometry.h:123
_RectangleInfo
Definition geometry.h:130
_ResampleFilter
Definition resample.c:94