2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 %
13 % MagickCore Image Distortion Methods %
21 % Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization %
22 % dedicated to making software imaging solutions freely available. %
24 % You may not use this file except in compliance with the License. You may %
25 % obtain a copy of the License at %
27 % http://www.imagemagick.org/script/license.php %
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. %
35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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/hashmap.h"
55 #include "MagickCore/image.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/pixel-private.h"
65 #include "MagickCore/resample.h"
66 #include "MagickCore/resample-private.h"
67 #include "MagickCore/registry.h"
68 #include "MagickCore/resource_.h"
69 #include "MagickCore/semaphore.h"
70 #include "MagickCore/shear.h"
71 #include "MagickCore/string_.h"
72 #include "MagickCore/string-private.h"
73 #include "MagickCore/thread-private.h"
74 #include "MagickCore/token.h"
75 #include "MagickCore/transform.h"
78 Numerous internal routines for image distortions.
80 static inline double MagickMin(const double x,const double y)
82 return( x < y ? x : y);
84 static inline double MagickMax(const double x,const double y)
86 return( x > y ? x : y);
89 static inline void AffineArgsToCoefficients(double *affine)
91 /* map external sx,ry,rx,sy,tx,ty to internal c0,c2,c4,c1,c3,c5 */
92 double tmp[4]; /* note indexes 0 and 5 remain unchanged */
93 tmp[0]=affine[1]; tmp[1]=affine[2]; tmp[2]=affine[3]; tmp[3]=affine[4];
94 affine[3]=tmp[0]; affine[1]=tmp[1]; affine[4]=tmp[2]; affine[2]=tmp[3];
97 static inline void CoefficientsToAffineArgs(double *coeff)
99 /* map internal c0,c1,c2,c3,c4,c5 to external sx,ry,rx,sy,tx,ty */
100 double tmp[4]; /* note indexes 0 and 5 remain unchanged */
101 tmp[0]=coeff[3]; tmp[1]=coeff[1]; tmp[2]=coeff[4]; tmp[3]=coeff[2];
102 coeff[1]=tmp[0]; coeff[2]=tmp[1]; coeff[3]=tmp[2]; coeff[4]=tmp[3];
104 static void InvertAffineCoefficients(const double *coeff,double *inverse)
106 /* From "Digital Image Warping" by George Wolberg, page 50 */
109 determinant=PerceptibleReciprocal(coeff[0]*coeff[4]-coeff[1]*coeff[3]);
110 inverse[0]=determinant*coeff[4];
111 inverse[1]=determinant*(-coeff[1]);
112 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[2]*coeff[4]);
113 inverse[3]=determinant*(-coeff[3]);
114 inverse[4]=determinant*coeff[0];
115 inverse[5]=determinant*(coeff[2]*coeff[3]-coeff[0]*coeff[5]);
118 static void InvertPerspectiveCoefficients(const double *coeff,
121 /* From "Digital Image Warping" by George Wolberg, page 53 */
124 determinant=PerceptibleReciprocal(coeff[0]*coeff[4]-coeff[3]*coeff[1]);
125 inverse[0]=determinant*(coeff[4]-coeff[7]*coeff[5]);
126 inverse[1]=determinant*(coeff[7]*coeff[2]-coeff[1]);
127 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[4]*coeff[2]);
128 inverse[3]=determinant*(coeff[6]*coeff[5]-coeff[3]);
129 inverse[4]=determinant*(coeff[0]-coeff[6]*coeff[2]);
130 inverse[5]=determinant*(coeff[3]*coeff[2]-coeff[0]*coeff[5]);
131 inverse[6]=determinant*(coeff[3]*coeff[7]-coeff[6]*coeff[4]);
132 inverse[7]=determinant*(coeff[6]*coeff[1]-coeff[0]*coeff[7]);
136 * Polynomial Term Defining Functions
138 * Order must either be an integer, or 1.5 to produce
139 * the 2 number_valuesal polynomial function...
140 * affine 1 (3) u = c0 + c1*x + c2*y
141 * bilinear 1.5 (4) u = '' + c3*x*y
142 * quadratic 2 (6) u = '' + c4*x*x + c5*y*y
143 * cubic 3 (10) u = '' + c6*x^3 + c7*x*x*y + c8*x*y*y + c9*y^3
144 * quartic 4 (15) u = '' + c10*x^4 + ... + c14*y^4
145 * quintic 5 (21) u = '' + c15*x^5 + ... + c20*y^5
146 * number in parenthesis minimum number of points needed.
147 * Anything beyond quintic, has not been implemented until
148 * a more automated way of determining terms is found.
150 * Note the slight re-ordering of the terms for a quadratic polynomial
151 * which is to allow the use of a bi-linear (order=1.5) polynomial.
152 * All the later polynomials are ordered simply from x^N to y^N
154 static size_t poly_number_terms(double order)
156 /* Return the number of terms for a 2d polynomial */
157 if ( order < 1 || order > 5 ||
158 ( order != floor(order) && (order-1.5) > MagickEpsilon) )
159 return 0; /* invalid polynomial order */
160 return((size_t) floor((order+1)*(order+2)/2));
163 static double poly_basis_fn(ssize_t n, double x, double y)
165 /* Return the result for this polynomial term */
167 case 0: return( 1.0 ); /* constant */
169 case 2: return( y ); /* affine order = 1 terms = 3 */
170 case 3: return( x*y ); /* bilinear order = 1.5 terms = 4 */
171 case 4: return( x*x );
172 case 5: return( y*y ); /* quadratic order = 2 terms = 6 */
173 case 6: return( x*x*x );
174 case 7: return( x*x*y );
175 case 8: return( x*y*y );
176 case 9: return( y*y*y ); /* cubic order = 3 terms = 10 */
177 case 10: return( x*x*x*x );
178 case 11: return( x*x*x*y );
179 case 12: return( x*x*y*y );
180 case 13: return( x*y*y*y );
181 case 14: return( y*y*y*y ); /* quartic order = 4 terms = 15 */
182 case 15: return( x*x*x*x*x );
183 case 16: return( x*x*x*x*y );
184 case 17: return( x*x*x*y*y );
185 case 18: return( x*x*y*y*y );
186 case 19: return( x*y*y*y*y );
187 case 20: return( y*y*y*y*y ); /* quintic order = 5 terms = 21 */
189 return( 0 ); /* should never happen */
191 static const char *poly_basis_str(ssize_t n)
193 /* return the result for this polynomial term */
195 case 0: return(""); /* constant */
196 case 1: return("*ii");
197 case 2: return("*jj"); /* affine order = 1 terms = 3 */
198 case 3: return("*ii*jj"); /* bilinear order = 1.5 terms = 4 */
199 case 4: return("*ii*ii");
200 case 5: return("*jj*jj"); /* quadratic order = 2 terms = 6 */
201 case 6: return("*ii*ii*ii");
202 case 7: return("*ii*ii*jj");
203 case 8: return("*ii*jj*jj");
204 case 9: return("*jj*jj*jj"); /* cubic order = 3 terms = 10 */
205 case 10: return("*ii*ii*ii*ii");
206 case 11: return("*ii*ii*ii*jj");
207 case 12: return("*ii*ii*jj*jj");
208 case 13: return("*ii*jj*jj*jj");
209 case 14: return("*jj*jj*jj*jj"); /* quartic order = 4 terms = 15 */
210 case 15: return("*ii*ii*ii*ii*ii");
211 case 16: return("*ii*ii*ii*ii*jj");
212 case 17: return("*ii*ii*ii*jj*jj");
213 case 18: return("*ii*ii*jj*jj*jj");
214 case 19: return("*ii*jj*jj*jj*jj");
215 case 20: return("*jj*jj*jj*jj*jj"); /* quintic order = 5 terms = 21 */
217 return( "UNKNOWN" ); /* should never happen */
219 static double poly_basis_dx(ssize_t n, double x, double y)
221 /* polynomial term for x derivative */
223 case 0: return( 0.0 ); /* constant */
224 case 1: return( 1.0 );
225 case 2: return( 0.0 ); /* affine order = 1 terms = 3 */
226 case 3: return( y ); /* bilinear order = 1.5 terms = 4 */
228 case 5: return( 0.0 ); /* quadratic order = 2 terms = 6 */
229 case 6: return( x*x );
230 case 7: return( x*y );
231 case 8: return( y*y );
232 case 9: return( 0.0 ); /* cubic order = 3 terms = 10 */
233 case 10: return( x*x*x );
234 case 11: return( x*x*y );
235 case 12: return( x*y*y );
236 case 13: return( y*y*y );
237 case 14: return( 0.0 ); /* quartic order = 4 terms = 15 */
238 case 15: return( x*x*x*x );
239 case 16: return( x*x*x*y );
240 case 17: return( x*x*y*y );
241 case 18: return( x*y*y*y );
242 case 19: return( y*y*y*y );
243 case 20: return( 0.0 ); /* quintic order = 5 terms = 21 */
245 return( 0.0 ); /* should never happen */
247 static double poly_basis_dy(ssize_t n, double x, double y)
249 /* polynomial term for y derivative */
251 case 0: return( 0.0 ); /* constant */
252 case 1: return( 0.0 );
253 case 2: return( 1.0 ); /* affine order = 1 terms = 3 */
254 case 3: return( x ); /* bilinear order = 1.5 terms = 4 */
255 case 4: return( 0.0 );
256 case 5: return( y ); /* quadratic order = 2 terms = 6 */
257 default: return( poly_basis_dx(n-1,x,y) ); /* weird but true */
259 /* NOTE: the only reason that last is not true for 'quadratic'
260 is due to the re-arrangement of terms to allow for 'bilinear'
265 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
269 % A f f i n e T r a n s f o r m I m a g e %
273 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
275 % AffineTransformImage() transforms an image as dictated by the affine matrix.
276 % It allocates the memory necessary for the new Image structure and returns
277 % a pointer to the new image.
279 % The format of the AffineTransformImage method is:
281 % Image *AffineTransformImage(const Image *image,
282 % AffineMatrix *affine_matrix,ExceptionInfo *exception)
284 % A description of each parameter follows:
286 % o image: the image.
288 % o affine_matrix: the affine matrix.
290 % o exception: return any errors or warnings in this structure.
293 MagickExport Image *AffineTransformImage(const Image *image,
294 const AffineMatrix *affine_matrix,ExceptionInfo *exception)
303 Affine transform image.
305 assert(image->signature == MagickSignature);
306 if (image->debug != MagickFalse)
307 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
308 assert(affine_matrix != (AffineMatrix *) NULL);
309 assert(exception != (ExceptionInfo *) NULL);
310 assert(exception->signature == MagickSignature);
311 distort[0]=affine_matrix->sx;
312 distort[1]=affine_matrix->rx;
313 distort[2]=affine_matrix->ry;
314 distort[3]=affine_matrix->sy;
315 distort[4]=affine_matrix->tx;
316 distort[5]=affine_matrix->ty;
317 deskew_image=DistortImage(image,AffineProjectionDistortion,6,distort,
318 MagickTrue,exception);
319 return(deskew_image);
323 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
327 + G e n e r a t e C o e f f i c i e n t s %
331 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
333 % GenerateCoefficients() takes user provided input arguments and generates
334 % the coefficients, needed to apply the specific distortion for either
335 % distorting images (generally using control points) or generating a color
336 % gradient from sparsely separated color points.
338 % The format of the GenerateCoefficients() method is:
340 % Image *GenerateCoefficients(const Image *image,DistortImageMethod method,
341 % const size_t number_arguments,const double *arguments,
342 % size_t number_values, ExceptionInfo *exception)
344 % A description of each parameter follows:
346 % o image: the image to be distorted.
348 % o method: the method of image distortion/ sparse gradient
350 % o number_arguments: the number of arguments given.
352 % o arguments: the arguments for this distortion method.
354 % o number_values: the style and format of given control points, (caller type)
355 % 0: 2 dimensional mapping of control points (Distort)
356 % Format: u,v,x,y where u,v is the 'source' of the
357 % the color to be plotted, for DistortImage()
358 % N: Interpolation of control points with N values (usally r,g,b)
359 % Format: x,y,r,g,b mapping x,y to color values r,g,b
360 % IN future, variable number of values may be given (1 to N)
362 % o exception: return any errors or warnings in this structure
364 % Note that the returned array of double values must be freed by the
365 % calling method using RelinquishMagickMemory(). This however may change in
366 % the future to require a more 'method' specific method.
368 % Because of this this method should not be classed as stable or used
369 % outside other MagickCore library methods.
372 static inline double MagickRound(double x)
375 Round the fraction to nearest integer.
377 if ((x-floor(x)) < (ceil(x)-x))
382 static double *GenerateCoefficients(const Image *image,
383 DistortImageMethod *method,const size_t number_arguments,
384 const double *arguments,size_t number_values,ExceptionInfo *exception)
393 number_coeff, /* number of coefficients to return (array size) */
394 cp_size, /* number floating point numbers per control point */
395 cp_x,cp_y, /* the x,y indexes for control point */
396 cp_values; /* index of values for this control point */
397 /* number_values Number of values given per control point */
399 if ( number_values == 0 ) {
400 /* Image distortion using control points (or other distortion)
401 That is generate a mapping so that x,y->u,v given u,v,x,y
403 number_values = 2; /* special case: two values of u,v */
404 cp_values = 0; /* the values i,j are BEFORE the destination CP x,y */
405 cp_x = 2; /* location of x,y in input control values */
407 /* NOTE: cp_values, also used for later 'reverse map distort' tests */
410 cp_x = 0; /* location of x,y in input control values */
412 cp_values = 2; /* and the other values are after x,y */
413 /* Typically in this case the values are R,G,B color values */
415 cp_size = number_values+2; /* each CP defintion involves this many numbers */
417 /* If not enough control point pairs are found for specific distortions
418 fall back to Affine distortion (allowing 0 to 3 point pairs)
420 if ( number_arguments < 4*cp_size &&
421 ( *method == BilinearForwardDistortion
422 || *method == BilinearReverseDistortion
423 || *method == PerspectiveDistortion
425 *method = AffineDistortion;
429 case AffineDistortion:
430 /* also BarycentricColorInterpolate: */
431 number_coeff=3*number_values;
433 case PolynomialDistortion:
434 /* number of coefficents depend on the given polynomal 'order' */
435 if ( number_arguments <= 1 && (number_arguments-1)%cp_size != 0)
437 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
438 "InvalidArgument","%s : '%s'","Polynomial",
439 "Invalid number of args: order [CPs]...");
440 return((double *) NULL);
442 i = poly_number_terms(arguments[0]);
443 number_coeff = 2 + i*number_values;
445 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
446 "InvalidArgument","%s : '%s'","Polynomial",
447 "Invalid order, should be interger 1 to 5, or 1.5");
448 return((double *) NULL);
450 if ( number_arguments < 1+i*cp_size ) {
451 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
452 "InvalidArgument", "%s : 'require at least %.20g CPs'",
453 "Polynomial", (double) i);
454 return((double *) NULL);
457 case BilinearReverseDistortion:
458 number_coeff=4*number_values;
461 The rest are constants as they are only used for image distorts
463 case BilinearForwardDistortion:
464 number_coeff=10; /* 2*4 coeff plus 2 constants */
465 cp_x = 0; /* Reverse src/dest coords for forward mapping */
470 case QuadraterialDistortion:
471 number_coeff=19; /* BilinearForward + BilinearReverse */
474 case ShepardsDistortion:
475 number_coeff=1; /* The power factor to use */
480 case ScaleRotateTranslateDistortion:
481 case AffineProjectionDistortion:
482 case Plane2CylinderDistortion:
483 case Cylinder2PlaneDistortion:
486 case PolarDistortion:
487 case DePolarDistortion:
490 case PerspectiveDistortion:
491 case PerspectiveProjectionDistortion:
494 case BarrelDistortion:
495 case BarrelInverseDistortion:
499 perror("unknown method given"); /* just fail assertion */
502 /* allocate the array of coefficients needed */
503 coeff = (double *) AcquireQuantumMemory(number_coeff,sizeof(*coeff));
504 if (coeff == (double *) NULL) {
505 (void) ThrowMagickException(exception,GetMagickModule(),
506 ResourceLimitError,"MemoryAllocationFailed",
507 "%s", "GenerateCoefficients");
508 return((double *) NULL);
511 /* zero out coefficients array */
512 for (i=0; i < number_coeff; i++)
517 case AffineDistortion:
521 for each 'value' given
523 Input Arguments are sets of control points...
524 For Distort Images u,v, x,y ...
525 For Sparse Gradients x,y, r,g,b ...
527 if ( number_arguments%cp_size != 0 ||
528 number_arguments < cp_size ) {
529 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
530 "InvalidArgument", "%s : 'require at least %.20g CPs'",
532 coeff=(double *) RelinquishMagickMemory(coeff);
533 return((double *) NULL);
535 /* handle special cases of not enough arguments */
536 if ( number_arguments == cp_size ) {
537 /* Only 1 CP Set Given */
538 if ( cp_values == 0 ) {
539 /* image distortion - translate the image */
541 coeff[2] = arguments[0] - arguments[2];
543 coeff[5] = arguments[1] - arguments[3];
546 /* sparse gradient - use the values directly */
547 for (i=0; i<number_values; i++)
548 coeff[i*3+2] = arguments[cp_values+i];
552 /* 2 or more points (usally 3) given.
553 Solve a least squares simultaneous equation for coefficients.
563 /* create matrix, and a fake vectors matrix */
564 matrix = AcquireMagickMatrix(3UL,3UL);
565 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
566 if (matrix == (double **) NULL || vectors == (double **) NULL)
568 matrix = RelinquishMagickMatrix(matrix, 3UL);
569 vectors = (double **) RelinquishMagickMemory(vectors);
570 coeff = (double *) RelinquishMagickMemory(coeff);
571 (void) ThrowMagickException(exception,GetMagickModule(),
572 ResourceLimitError,"MemoryAllocationFailed",
573 "%s", "DistortCoefficients");
574 return((double *) NULL);
576 /* fake a number_values x3 vectors matrix from coefficients array */
577 for (i=0; i < number_values; i++)
578 vectors[i] = &(coeff[i*3]);
579 /* Add given control point pairs for least squares solving */
580 for (i=0; i < number_arguments; i+=cp_size) {
581 terms[0] = arguments[i+cp_x]; /* x */
582 terms[1] = arguments[i+cp_y]; /* y */
583 terms[2] = 1; /* 1 */
584 LeastSquaresAddTerms(matrix,vectors,terms,
585 &(arguments[i+cp_values]),3UL,number_values);
587 if ( number_arguments == 2*cp_size ) {
588 /* Only two pairs were given, but we need 3 to solve the affine.
589 Fake extra coordinates by rotating p1 around p0 by 90 degrees.
590 x2 = x0 - (y1-y0) y2 = y0 + (x1-x0)
592 terms[0] = arguments[cp_x]
593 - ( arguments[cp_size+cp_y] - arguments[cp_y] ); /* x2 */
594 terms[1] = arguments[cp_y] +
595 + ( arguments[cp_size+cp_x] - arguments[cp_x] ); /* y2 */
596 terms[2] = 1; /* 1 */
597 if ( cp_values == 0 ) {
598 /* Image Distortion - rotate the u,v coordients too */
601 uv2[0] = arguments[0] - arguments[5] + arguments[1]; /* u2 */
602 uv2[1] = arguments[1] + arguments[4] - arguments[0]; /* v2 */
603 LeastSquaresAddTerms(matrix,vectors,terms,uv2,3UL,2UL);
606 /* Sparse Gradient - use values of p0 for linear gradient */
607 LeastSquaresAddTerms(matrix,vectors,terms,
608 &(arguments[cp_values]),3UL,number_values);
611 /* Solve for LeastSquares Coefficients */
612 status=GaussJordanElimination(matrix,vectors,3UL,number_values);
613 matrix = RelinquishMagickMatrix(matrix, 3UL);
614 vectors = (double **) RelinquishMagickMemory(vectors);
615 if ( status == MagickFalse ) {
616 coeff = (double *) RelinquishMagickMemory(coeff);
617 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
618 "InvalidArgument","%s : 'Unsolvable Matrix'",
619 CommandOptionToMnemonic(MagickDistortOptions, *method) );
620 return((double *) NULL);
625 case AffineProjectionDistortion:
628 Arguments: Affine Matrix (forward mapping)
629 Arguments sx, rx, ry, sy, tx, ty
630 Where u = sx*x + ry*y + tx
633 Returns coefficients (in there inverse form) ordered as...
636 AffineProjection Distortion Notes...
637 + Will only work with a 2 number_values for Image Distortion
638 + Can not be used for generating a sparse gradient (interpolation)
641 if (number_arguments != 6) {
642 coeff = (double *) RelinquishMagickMemory(coeff);
643 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
644 "InvalidArgument","%s : 'Needs 6 coeff values'",
645 CommandOptionToMnemonic(MagickDistortOptions, *method) );
646 return((double *) NULL);
648 /* FUTURE: trap test for sx*sy-rx*ry == 0 (determinant = 0, no inverse) */
649 for(i=0; i<6UL; i++ )
650 inverse[i] = arguments[i];
651 AffineArgsToCoefficients(inverse); /* map into coefficents */
652 InvertAffineCoefficients(inverse, coeff); /* invert */
653 *method = AffineDistortion;
657 case ScaleRotateTranslateDistortion:
659 /* Scale, Rotate and Translate Distortion
660 An alternative Affine Distortion
661 Argument options, by number of arguments given:
662 7: x,y, sx,sy, a, nx,ny
669 Where actions are (in order of application)
670 x,y 'center' of transforms (default = image center)
671 sx,sy scale image by this amount (default = 1)
672 a angle of rotation (argument required)
673 nx,ny move 'center' here (default = x,y or no movement)
674 And convert to affine mapping coefficients
676 ScaleRotateTranslate Distortion Notes...
677 + Does not use a set of CPs in any normal way
678 + Will only work with a 2 number_valuesal Image Distortion
679 + Cannot be used for generating a sparse gradient (interpolation)
685 /* set default center, and default scale */
686 x = nx = (double)(image->columns)/2.0 + (double)image->page.x;
687 y = ny = (double)(image->rows)/2.0 + (double)image->page.y;
689 switch ( number_arguments ) {
691 coeff = (double *) RelinquishMagickMemory(coeff);
692 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
693 "InvalidArgument","%s : 'Needs at least 1 argument'",
694 CommandOptionToMnemonic(MagickDistortOptions, *method) );
695 return((double *) NULL);
700 sx = sy = arguments[0];
704 x = nx = arguments[0];
705 y = ny = arguments[1];
706 switch ( number_arguments ) {
711 sx = sy = arguments[2];
720 sx = sy = arguments[2];
733 coeff = (double *) RelinquishMagickMemory(coeff);
734 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
735 "InvalidArgument","%s : 'Too Many Arguments (7 or less)'",
736 CommandOptionToMnemonic(MagickDistortOptions, *method) );
737 return((double *) NULL);
741 /* Trap if sx or sy == 0 -- image is scaled out of existance! */
742 if ( fabs(sx) < MagickEpsilon || fabs(sy) < MagickEpsilon ) {
743 coeff = (double *) RelinquishMagickMemory(coeff);
744 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
745 "InvalidArgument","%s : 'Zero Scale Given'",
746 CommandOptionToMnemonic(MagickDistortOptions, *method) );
747 return((double *) NULL);
749 /* Save the given arguments as an affine distortion */
750 a=DegreesToRadians(a); cosine=cos(a); sine=sin(a);
752 *method = AffineDistortion;
755 coeff[2]=x-nx*coeff[0]-ny*coeff[1];
758 coeff[5]=y-nx*coeff[3]-ny*coeff[4];
761 case PerspectiveDistortion:
763 Perspective Distortion (a ratio of affine distortions)
765 p(x,y) c0*x + c1*y + c2
766 u = ------ = ------------------
767 r(x,y) c6*x + c7*y + 1
769 q(x,y) c3*x + c4*y + c5
770 v = ------ = ------------------
771 r(x,y) c6*x + c7*y + 1
773 c8 = Sign of 'r', or the denominator affine, for the actual image.
774 This determines what part of the distorted image is 'ground'
775 side of the horizon, the other part is 'sky' or invalid.
776 Valid values are +1.0 or -1.0 only.
778 Input Arguments are sets of control points...
779 For Distort Images u,v, x,y ...
780 For Sparse Gradients x,y, r,g,b ...
782 Perspective Distortion Notes...
783 + Can be thought of as ratio of 3 affine transformations
784 + Not separatable: r() or c6 and c7 are used by both equations
785 + All 8 coefficients must be determined simultaniously
786 + Will only work with a 2 number_valuesal Image Distortion
787 + Can not be used for generating a sparse gradient (interpolation)
788 + It is not linear, but is simple to generate an inverse
789 + All lines within an image remain lines.
790 + but distances between points may vary.
804 if ( number_arguments%cp_size != 0 ||
805 number_arguments < cp_size*4 ) {
806 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
807 "InvalidArgument", "%s : 'require at least %.20g CPs'",
808 CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0);
809 coeff=(double *) RelinquishMagickMemory(coeff);
810 return((double *) NULL);
812 /* fake 1x8 vectors matrix directly using the coefficients array */
813 vectors[0] = &(coeff[0]);
814 /* 8x8 least-squares matrix (zeroed) */
815 matrix = AcquireMagickMatrix(8UL,8UL);
816 if (matrix == (double **) NULL) {
817 (void) ThrowMagickException(exception,GetMagickModule(),
818 ResourceLimitError,"MemoryAllocationFailed",
819 "%s", "DistortCoefficients");
820 return((double *) NULL);
822 /* Add control points for least squares solving */
823 for (i=0; i < number_arguments; i+=4) {
824 terms[0]=arguments[i+cp_x]; /* c0*x */
825 terms[1]=arguments[i+cp_y]; /* c1*y */
826 terms[2]=1.0; /* c2*1 */
830 terms[6]=-terms[0]*arguments[i+cp_u]; /* 1/(c6*x) */
831 terms[7]=-terms[1]*arguments[i+cp_u]; /* 1/(c7*y) */
832 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_u]),
838 terms[3]=arguments[i+cp_x]; /* c3*x */
839 terms[4]=arguments[i+cp_y]; /* c4*y */
840 terms[5]=1.0; /* c5*1 */
841 terms[6]=-terms[3]*arguments[i+cp_v]; /* 1/(c6*x) */
842 terms[7]=-terms[4]*arguments[i+cp_v]; /* 1/(c7*y) */
843 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_v]),
846 /* Solve for LeastSquares Coefficients */
847 status=GaussJordanElimination(matrix,vectors,8UL,1UL);
848 matrix = RelinquishMagickMatrix(matrix, 8UL);
849 if ( status == MagickFalse ) {
850 coeff = (double *) RelinquishMagickMemory(coeff);
851 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
852 "InvalidArgument","%s : 'Unsolvable Matrix'",
853 CommandOptionToMnemonic(MagickDistortOptions, *method) );
854 return((double *) NULL);
857 Calculate 9'th coefficient! The ground-sky determination.
858 What is sign of the 'ground' in r() denominator affine function?
859 Just use any valid image coordinate (first control point) in
860 destination for determination of what part of view is 'ground'.
862 coeff[8] = coeff[6]*arguments[cp_x]
863 + coeff[7]*arguments[cp_y] + 1.0;
864 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
868 case PerspectiveProjectionDistortion:
871 Arguments: Perspective Coefficents (forward mapping)
873 if (number_arguments != 8) {
874 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
875 "InvalidArgument", "%s : 'Needs 8 coefficient values'",
876 CommandOptionToMnemonic(MagickDistortOptions, *method));
877 return((double *) NULL);
879 /* FUTURE: trap test c0*c4-c3*c1 == 0 (determinate = 0, no inverse) */
880 InvertPerspectiveCoefficients(arguments, coeff);
882 Calculate 9'th coefficient! The ground-sky determination.
883 What is sign of the 'ground' in r() denominator affine function?
884 Just use any valid image cocodinate in destination for determination.
885 For a forward mapped perspective the images 0,0 coord will map to
886 c2,c5 in the distorted image, so set the sign of denominator of that.
888 coeff[8] = coeff[6]*arguments[2]
889 + coeff[7]*arguments[5] + 1.0;
890 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
891 *method = PerspectiveDistortion;
895 case BilinearForwardDistortion:
896 case BilinearReverseDistortion:
898 /* Bilinear Distortion (Forward mapping)
899 v = c0*x + c1*y + c2*x*y + c3;
900 for each 'value' given
902 This is actually a simple polynomial Distortion! The difference
903 however is when we need to reverse the above equation to generate a
904 BilinearForwardDistortion (see below).
906 Input Arguments are sets of control points...
907 For Distort Images u,v, x,y ...
908 For Sparse Gradients x,y, r,g,b ...
919 /* check the number of arguments */
920 if ( number_arguments%cp_size != 0 ||
921 number_arguments < cp_size*4 ) {
922 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
923 "InvalidArgument", "%s : 'require at least %.20g CPs'",
924 CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0);
925 coeff=(double *) RelinquishMagickMemory(coeff);
926 return((double *) NULL);
928 /* create matrix, and a fake vectors matrix */
929 matrix = AcquireMagickMatrix(4UL,4UL);
930 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
931 if (matrix == (double **) NULL || vectors == (double **) NULL)
933 matrix = RelinquishMagickMatrix(matrix, 4UL);
934 vectors = (double **) RelinquishMagickMemory(vectors);
935 coeff = (double *) RelinquishMagickMemory(coeff);
936 (void) ThrowMagickException(exception,GetMagickModule(),
937 ResourceLimitError,"MemoryAllocationFailed",
938 "%s", "DistortCoefficients");
939 return((double *) NULL);
941 /* fake a number_values x4 vectors matrix from coefficients array */
942 for (i=0; i < number_values; i++)
943 vectors[i] = &(coeff[i*4]);
944 /* Add given control point pairs for least squares solving */
945 for (i=0; i < number_arguments; i+=cp_size) {
946 terms[0] = arguments[i+cp_x]; /* x */
947 terms[1] = arguments[i+cp_y]; /* y */
948 terms[2] = terms[0]*terms[1]; /* x*y */
949 terms[3] = 1; /* 1 */
950 LeastSquaresAddTerms(matrix,vectors,terms,
951 &(arguments[i+cp_values]),4UL,number_values);
953 /* Solve for LeastSquares Coefficients */
954 status=GaussJordanElimination(matrix,vectors,4UL,number_values);
955 matrix = RelinquishMagickMatrix(matrix, 4UL);
956 vectors = (double **) RelinquishMagickMemory(vectors);
957 if ( status == MagickFalse ) {
958 coeff = (double *) RelinquishMagickMemory(coeff);
959 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
960 "InvalidArgument","%s : 'Unsolvable Matrix'",
961 CommandOptionToMnemonic(MagickDistortOptions, *method) );
962 return((double *) NULL);
964 if ( *method == BilinearForwardDistortion ) {
965 /* Bilinear Forward Mapped Distortion
967 The above least-squares solved for coefficents but in the forward
968 direction, due to changes to indexing constants.
970 i = c0*x + c1*y + c2*x*y + c3;
971 j = c4*x + c5*y + c6*x*y + c7;
973 where i,j are in the destination image, NOT the source.
975 Reverse Pixel mapping however needs to use reverse of these
976 functions. It required a full page of algbra to work out the
977 reversed mapping formula, but resolves down to the following...
980 c9 = 2*(c2*c5-c1*c6); // '2*a' in the quadratic formula
982 i = i - c3; j = j - c7;
983 b = c6*i - c2*j + c8; // So that a*y^2 + b*y + c == 0
984 c = c4*i - c0*j; // y = ( -b +- sqrt(bb - 4ac) ) / (2*a)
988 y = ( -b + sqrt(r) ) / c9;
992 x = ( i - c1*y) / ( c1 - c2*y );
994 NB: if 'r' is negative there is no solution!
995 NB: the sign of the sqrt() should be negative if image becomes
996 flipped or flopped, or crosses over itself.
997 NB: techniqually coefficient c5 is not needed, anymore,
998 but kept for completness.
1000 See Anthony Thyssen <A.Thyssen@griffith.edu.au>
1001 or Fred Weinhaus <fmw@alink.net> for more details.
1004 coeff[8] = coeff[0]*coeff[5] - coeff[1]*coeff[4];
1005 coeff[9] = 2*(coeff[2]*coeff[5] - coeff[1]*coeff[6]);
1010 case QuadrilateralDistortion:
1012 /* Map a Quadrilateral to a unit square using BilinearReverse
1013 Then map that unit square back to the final Quadrilateral
1014 using BilinearForward.
1016 Input Arguments are sets of control points...
1017 For Distort Images u,v, x,y ...
1018 For Sparse Gradients x,y, r,g,b ...
1021 /* UNDER CONSTRUCTION */
1026 case PolynomialDistortion:
1028 /* Polynomial Distortion
1030 First two coefficents are used to hole global polynomal information
1031 c0 = Order of the polynimial being created
1032 c1 = number_of_terms in one polynomial equation
1034 Rest of the coefficients map to the equations....
1035 v = c0 + c1*x + c2*y + c3*x*y + c4*x^2 + c5*y^2 + c6*x^3 + ...
1036 for each 'value' (number_values of them) given.
1037 As such total coefficients = 2 + number_terms * number_values
1039 Input Arguments are sets of control points...
1040 For Distort Images order [u,v, x,y] ...
1041 For Sparse Gradients order [x,y, r,g,b] ...
1043 Polynomial Distortion Notes...
1044 + UNDER DEVELOPMENT -- Do not expect this to remain as is.
1045 + Currently polynomial is a reversed mapped distortion.
1046 + Order 1.5 is fudged to map into a bilinear distortion.
1047 though it is not the same order as that distortion.
1055 nterms; /* number of polynomial terms per number_values */
1063 /* first two coefficients hold polynomial order information */
1064 coeff[0] = arguments[0];
1065 coeff[1] = (double) poly_number_terms(arguments[0]);
1066 nterms = (size_t) coeff[1];
1068 /* create matrix, a fake vectors matrix, and least sqs terms */
1069 matrix = AcquireMagickMatrix(nterms,nterms);
1070 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
1071 terms = (double *) AcquireQuantumMemory(nterms, sizeof(*terms));
1072 if (matrix == (double **) NULL ||
1073 vectors == (double **) NULL ||
1074 terms == (double *) NULL )
1076 matrix = RelinquishMagickMatrix(matrix, nterms);
1077 vectors = (double **) RelinquishMagickMemory(vectors);
1078 terms = (double *) RelinquishMagickMemory(terms);
1079 coeff = (double *) RelinquishMagickMemory(coeff);
1080 (void) ThrowMagickException(exception,GetMagickModule(),
1081 ResourceLimitError,"MemoryAllocationFailed",
1082 "%s", "DistortCoefficients");
1083 return((double *) NULL);
1085 /* fake a number_values x3 vectors matrix from coefficients array */
1086 for (i=0; i < number_values; i++)
1087 vectors[i] = &(coeff[2+i*nterms]);
1088 /* Add given control point pairs for least squares solving */
1089 for (i=1; i < number_arguments; i+=cp_size) { /* NB: start = 1 not 0 */
1090 for (j=0; j < (ssize_t) nterms; j++)
1091 terms[j] = poly_basis_fn(j,arguments[i+cp_x],arguments[i+cp_y]);
1092 LeastSquaresAddTerms(matrix,vectors,terms,
1093 &(arguments[i+cp_values]),nterms,number_values);
1095 terms = (double *) RelinquishMagickMemory(terms);
1096 /* Solve for LeastSquares Coefficients */
1097 status=GaussJordanElimination(matrix,vectors,nterms,number_values);
1098 matrix = RelinquishMagickMatrix(matrix, nterms);
1099 vectors = (double **) RelinquishMagickMemory(vectors);
1100 if ( status == MagickFalse ) {
1101 coeff = (double *) RelinquishMagickMemory(coeff);
1102 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1103 "InvalidArgument","%s : 'Unsolvable Matrix'",
1104 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1105 return((double *) NULL);
1112 Args: arc_width rotate top_edge_radius bottom_edge_radius
1113 All but first argument are optional
1114 arc_width The angle over which to arc the image side-to-side
1115 rotate Angle to rotate image from vertical center
1116 top_radius Set top edge of source image at this radius
1117 bottom_radius Set bootom edge to this radius (radial scaling)
1119 By default, if the radii arguments are nor provided the image radius
1120 is calculated so the horizontal center-line is fits the given arc
1123 The output image size is ALWAYS adjusted to contain the whole image,
1124 and an offset is given to position image relative to the 0,0 point of
1125 the origin, allowing users to use relative positioning onto larger
1126 background (via -flatten).
1128 The arguments are converted to these coefficients
1129 c0: angle for center of source image
1130 c1: angle scale for mapping to source image
1131 c2: radius for top of source image
1132 c3: radius scale for mapping source image
1133 c4: centerline of arc within source image
1135 Note the coefficients use a center angle, so asymptotic join is
1136 furthest from both sides of the source image. This also means that
1137 for arc angles greater than 360 the sides of the image will be
1140 Arc Distortion Notes...
1141 + Does not use a set of CPs
1142 + Will only work with Image Distortion
1143 + Can not be used for generating a sparse gradient (interpolation)
1145 if ( number_arguments >= 1 && arguments[0] < MagickEpsilon ) {
1146 coeff = (double *) RelinquishMagickMemory(coeff);
1147 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1148 "InvalidArgument","%s : 'Arc Angle Too Small'",
1149 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1150 return((double *) NULL);
1152 if ( number_arguments >= 3 && arguments[2] < MagickEpsilon ) {
1153 coeff = (double *) RelinquishMagickMemory(coeff);
1154 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1155 "InvalidArgument","%s : 'Outer Radius Too Small'",
1156 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1157 return((double *) NULL);
1159 coeff[0] = -MagickPI2; /* -90, place at top! */
1160 if ( number_arguments >= 1 )
1161 coeff[1] = DegreesToRadians(arguments[0]);
1163 coeff[1] = MagickPI2; /* zero arguments - center is at top */
1164 if ( number_arguments >= 2 )
1165 coeff[0] += DegreesToRadians(arguments[1]);
1166 coeff[0] /= Magick2PI; /* normalize radians */
1167 coeff[0] -= MagickRound(coeff[0]);
1168 coeff[0] *= Magick2PI; /* de-normalize back to radians */
1169 coeff[3] = (double)image->rows-1;
1170 coeff[2] = (double)image->columns/coeff[1] + coeff[3]/2.0;
1171 if ( number_arguments >= 3 ) {
1172 if ( number_arguments >= 4 )
1173 coeff[3] = arguments[2] - arguments[3];
1175 coeff[3] *= arguments[2]/coeff[2];
1176 coeff[2] = arguments[2];
1178 coeff[4] = ((double)image->columns-1.0)/2.0;
1182 case PolarDistortion:
1183 case DePolarDistortion:
1185 /* (De)Polar Distortion (same set of arguments)
1186 Args: Rmax, Rmin, Xcenter,Ycenter, Afrom,Ato
1187 DePolar can also have the extra arguments of Width, Height
1189 Coefficients 0 to 5 is the sanatized version first 6 input args
1190 Coefficient 6 is the angle to coord ratio and visa-versa
1191 Coefficient 7 is the radius to coord ratio and visa-versa
1193 WARNING: It is possible for Radius max<min and/or Angle from>to
1195 if ( number_arguments == 3
1196 || ( number_arguments > 6 && *method == PolarDistortion )
1197 || number_arguments > 8 ) {
1198 (void) ThrowMagickException(exception,GetMagickModule(),
1199 OptionError,"InvalidArgument", "%s : number of arguments",
1200 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1201 coeff=(double *) RelinquishMagickMemory(coeff);
1202 return((double *) NULL);
1204 /* Rmax - if 0 calculate appropriate value */
1205 if ( number_arguments >= 1 )
1206 coeff[0] = arguments[0];
1209 /* Rmin - usally 0 */
1210 coeff[1] = number_arguments >= 2 ? arguments[1] : 0.0;
1212 if ( number_arguments >= 4 ) {
1213 coeff[2] = arguments[2];
1214 coeff[3] = arguments[3];
1216 else { /* center of actual image */
1217 coeff[2] = (double)(image->columns)/2.0+image->page.x;
1218 coeff[3] = (double)(image->rows)/2.0+image->page.y;
1220 /* Angle from,to - about polar center 0 is downward */
1221 coeff[4] = -MagickPI;
1222 if ( number_arguments >= 5 )
1223 coeff[4] = DegreesToRadians(arguments[4]);
1224 coeff[5] = coeff[4];
1225 if ( number_arguments >= 6 )
1226 coeff[5] = DegreesToRadians(arguments[5]);
1227 if ( fabs(coeff[4]-coeff[5]) < MagickEpsilon )
1228 coeff[5] += Magick2PI; /* same angle is a full circle */
1229 /* if radius 0 or negative, its a special value... */
1230 if ( coeff[0] < MagickEpsilon ) {
1231 /* Use closest edge if radius == 0 */
1232 if ( fabs(coeff[0]) < MagickEpsilon ) {
1233 coeff[0]=MagickMin(fabs(coeff[2]-image->page.x),
1234 fabs(coeff[3]-image->page.y));
1235 coeff[0]=MagickMin(coeff[0],
1236 fabs(coeff[2]-image->page.x-image->columns));
1237 coeff[0]=MagickMin(coeff[0],
1238 fabs(coeff[3]-image->page.y-image->rows));
1240 /* furthest diagonal if radius == -1 */
1241 if ( fabs(-1.0-coeff[0]) < MagickEpsilon ) {
1243 rx = coeff[2]-image->page.x;
1244 ry = coeff[3]-image->page.y;
1245 coeff[0] = rx*rx+ry*ry;
1246 ry = coeff[3]-image->page.y-image->rows;
1247 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1248 rx = coeff[2]-image->page.x-image->columns;
1249 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1250 ry = coeff[3]-image->page.y;
1251 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1252 coeff[0] = sqrt(coeff[0]);
1255 /* IF Rmax <= 0 or Rmin < 0 OR Rmax < Rmin, THEN error */
1256 if ( coeff[0] < MagickEpsilon || coeff[1] < -MagickEpsilon
1257 || (coeff[0]-coeff[1]) < MagickEpsilon ) {
1258 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1259 "InvalidArgument", "%s : Invalid Radius",
1260 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1261 coeff=(double *) RelinquishMagickMemory(coeff);
1262 return((double *) NULL);
1264 /* converstion ratios */
1265 if ( *method == PolarDistortion ) {
1266 coeff[6]=(double) image->columns/(coeff[5]-coeff[4]);
1267 coeff[7]=(double) image->rows/(coeff[0]-coeff[1]);
1269 else { /* *method == DePolarDistortion */
1270 coeff[6]=(coeff[5]-coeff[4])/image->columns;
1271 coeff[7]=(coeff[0]-coeff[1])/image->rows;
1275 case Cylinder2PlaneDistortion:
1276 case Plane2CylinderDistortion:
1278 /* 3D Cylinder to/from a Tangential Plane
1280 Projection between a clinder and flat plain from a point on the
1281 center line of the cylinder.
1283 The two surfaces coincide in 3D space at the given centers of
1284 distortion (perpendicular to projection point) on both images.
1287 Coefficents: FOV(radians), Radius, center_x,y, dest_center_x,y
1289 FOV (Field Of View) the angular field of view of the distortion,
1290 across the width of the image, in degrees. The centers are the
1291 points of least distortion in the input and resulting images.
1293 These centers are however determined later.
1295 Coeff 0 is the FOV angle of view of image width in radians
1296 Coeff 1 is calculated radius of cylinder.
1297 Coeff 2,3 center of distortion of input image
1298 Coefficents 4,5 Center of Distortion of dest (determined later)
1300 if ( arguments[0] < MagickEpsilon || arguments[0] > 160.0 ) {
1301 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1302 "InvalidArgument", "%s : Invalid FOV Angle",
1303 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1304 coeff=(double *) RelinquishMagickMemory(coeff);
1305 return((double *) NULL);
1307 coeff[0] = DegreesToRadians(arguments[0]);
1308 if ( *method == Cylinder2PlaneDistortion )
1309 /* image is curved around cylinder, so FOV angle (in radians)
1310 * scales directly to image X coordinate, according to its radius.
1312 coeff[1] = (double) image->columns/coeff[0];
1314 /* radius is distance away from an image with this angular FOV */
1315 coeff[1] = (double) image->columns / ( 2 * tan(coeff[0]/2) );
1317 coeff[2] = (double)(image->columns)/2.0+image->page.x;
1318 coeff[3] = (double)(image->rows)/2.0+image->page.y;
1319 coeff[4] = coeff[2];
1320 coeff[5] = coeff[3]; /* assuming image size is the same */
1323 case BarrelDistortion:
1324 case BarrelInverseDistortion:
1326 /* Barrel Distortion
1327 Rs=(A*Rd^3 + B*Rd^2 + C*Rd + D)*Rd
1328 BarrelInv Distortion
1329 Rs=Rd/(A*Rd^3 + B*Rd^2 + C*Rd + D)
1331 Where Rd is the normalized radius from corner to middle of image
1332 Input Arguments are one of the following forms (number of arguments)...
1337 8: Ax,Bx,Cx,Dx Ay,By,Cy,Dy
1338 10: Ax,Bx,Cx,Dx Ay,By,Cy,Dy X,Y
1340 Returns 10 coefficent values, which are de-normalized (pixel scale)
1341 Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Xc, Yc
1343 /* Radius de-normalization scaling factor */
1345 rscale = 2.0/MagickMin((double) image->columns,(double) image->rows);
1347 /* sanity check number of args must = 3,4,5,6,8,10 or error */
1348 if ( (number_arguments < 3) || (number_arguments == 7) ||
1349 (number_arguments == 9) || (number_arguments > 10) )
1351 coeff=(double *) RelinquishMagickMemory(coeff);
1352 (void) ThrowMagickException(exception,GetMagickModule(),
1353 OptionError,"InvalidArgument", "%s : number of arguments",
1354 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1355 return((double *) NULL);
1357 /* A,B,C,D coefficients */
1358 coeff[0] = arguments[0];
1359 coeff[1] = arguments[1];
1360 coeff[2] = arguments[2];
1361 if ((number_arguments == 3) || (number_arguments == 5) )
1362 coeff[3] = 1.0 - coeff[0] - coeff[1] - coeff[2];
1364 coeff[3] = arguments[3];
1365 /* de-normalize the coefficients */
1366 coeff[0] *= pow(rscale,3.0);
1367 coeff[1] *= rscale*rscale;
1369 /* Y coefficients: as given OR same as X coefficients */
1370 if ( number_arguments >= 8 ) {
1371 coeff[4] = arguments[4] * pow(rscale,3.0);
1372 coeff[5] = arguments[5] * rscale*rscale;
1373 coeff[6] = arguments[6] * rscale;
1374 coeff[7] = arguments[7];
1377 coeff[4] = coeff[0];
1378 coeff[5] = coeff[1];
1379 coeff[6] = coeff[2];
1380 coeff[7] = coeff[3];
1382 /* X,Y Center of Distortion (image coodinates) */
1383 if ( number_arguments == 5 ) {
1384 coeff[8] = arguments[3];
1385 coeff[9] = arguments[4];
1387 else if ( number_arguments == 6 ) {
1388 coeff[8] = arguments[4];
1389 coeff[9] = arguments[5];
1391 else if ( number_arguments == 10 ) {
1392 coeff[8] = arguments[8];
1393 coeff[9] = arguments[9];
1396 /* center of the image provided (image coodinates) */
1397 coeff[8] = (double)image->columns/2.0 + image->page.x;
1398 coeff[9] = (double)image->rows/2.0 + image->page.y;
1402 case ShepardsDistortion:
1404 /* Shepards Distortion input arguments are the coefficents!
1405 Just check the number of arguments is valid!
1406 Args: u1,v1, x1,y1, ...
1407 OR : u1,v1, r1,g1,c1, ...
1409 if ( number_arguments%cp_size != 0 ||
1410 number_arguments < cp_size ) {
1411 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1412 "InvalidArgument", "%s : 'requires CP's (4 numbers each)'",
1413 CommandOptionToMnemonic(MagickDistortOptions, *method));
1414 coeff=(double *) RelinquishMagickMemory(coeff);
1415 return((double *) NULL);
1417 /* User defined weighting power for Shepard's Method */
1418 { const char *artifact=GetImageArtifact(image,"shepards:power");
1419 if ( artifact != (const char *) NULL ) {
1420 coeff[0]=StringToDouble(artifact,(char **) NULL) / 2.0;
1421 if ( coeff[0] < MagickEpsilon ) {
1422 (void) ThrowMagickException(exception,GetMagickModule(),
1423 OptionError,"InvalidArgument","%s", "-define shepards:power" );
1424 coeff=(double *) RelinquishMagickMemory(coeff);
1425 return((double *) NULL);
1429 coeff[0]=1.0; /* Default power of 2 (Inverse Squared) */
1436 /* you should never reach this point */
1437 perror("no method handler"); /* just fail assertion */
1438 return((double *) NULL);
1442 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1446 + D i s t o r t R e s i z e I m a g e %
1450 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1452 % DistortResizeImage() resize image using the equivalent but slower image
1453 % distortion operator. The filter is applied using a EWA cylindrical
1454 % resampling. But like resize the final image size is limited to whole pixels
1455 % with no effects by virtual-pixels on the result.
1457 % Note that images containing a transparency channel will be twice as slow to
1458 % resize as images one without transparency.
1460 % The format of the DistortResizeImage method is:
1462 % Image *AdaptiveResizeImage(const Image *image,const size_t columns,
1463 % const size_t rows,ExceptionInfo *exception)
1465 % A description of each parameter follows:
1467 % o image: the image.
1469 % o columns: the number of columns in the resized image.
1471 % o rows: the number of rows in the resized image.
1473 % o exception: return any errors or warnings in this structure.
1476 MagickExport Image *DistortResizeImage(const Image *image,
1477 const size_t columns,const size_t rows,ExceptionInfo *exception)
1479 #define DistortResizeImageTag "Distort/Image"
1495 Distort resize image.
1497 assert(image != (const Image *) NULL);
1498 assert(image->signature == MagickSignature);
1499 if (image->debug != MagickFalse)
1500 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1501 assert(exception != (ExceptionInfo *) NULL);
1502 assert(exception->signature == MagickSignature);
1503 if ((columns == 0) || (rows == 0))
1504 return((Image *) NULL);
1505 /* Do not short-circuit this resize if final image size is unchanged */
1507 (void) ResetMagickMemory(distort_args,0,12*sizeof(double));
1508 distort_args[4]=(double) image->columns;
1509 distort_args[6]=(double) columns;
1510 distort_args[9]=(double) image->rows;
1511 distort_args[11]=(double) rows;
1513 vp_save=GetImageVirtualPixelMethod(image);
1515 tmp_image=CloneImage(image,0,0,MagickTrue,exception);
1516 if ( tmp_image == (Image *) NULL )
1517 return((Image *) NULL);
1518 (void) SetImageVirtualPixelMethod(tmp_image,TransparentVirtualPixelMethod,
1520 tmp_image->image_info=image->image_info; /* preserve global options */
1522 if (image->alpha_trait != BlendPixelTrait)
1525 Image has not transparency channel, so we free to use it
1527 (void) SetImageAlphaChannel(tmp_image,SetAlphaChannel,exception);
1528 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1529 MagickTrue,exception),
1531 tmp_image=DestroyImage(tmp_image);
1532 if ( resize_image == (Image *) NULL )
1533 return((Image *) NULL);
1535 (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel,exception);
1540 Image has transparency so handle colors and alpha separatly.
1541 Basically we need to separate Virtual-Pixel alpha in the resized
1542 image, so only the actual original images alpha channel is used.
1544 distort alpha channel separately
1549 (void) SetImageAlphaChannel(tmp_image,ExtractAlphaChannel,exception);
1550 (void) SetImageAlphaChannel(tmp_image,OpaqueAlphaChannel,exception);
1551 resize_alpha=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1552 MagickTrue,exception),
1553 tmp_image=DestroyImage(tmp_image);
1554 if (resize_alpha == (Image *) NULL)
1555 return((Image *) NULL);
1557 /* distort the actual image containing alpha + VP alpha */
1558 tmp_image=CloneImage(image,0,0,MagickTrue,exception);
1559 if ( tmp_image == (Image *) NULL )
1560 return((Image *) NULL);
1561 tmp_image->image_info=image->image_info; /* preserve global options */
1562 (void) SetImageVirtualPixelMethod(tmp_image,
1563 TransparentVirtualPixelMethod,exception);
1564 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1565 MagickTrue,exception),
1566 tmp_image=DestroyImage(tmp_image);
1567 if ( resize_image == (Image *) NULL)
1569 resize_alpha=DestroyImage(resize_alpha);
1570 return((Image *) NULL);
1572 /* replace resize images alpha with the separally distorted alpha */
1573 (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel,
1575 (void) SetImageAlphaChannel(resize_alpha,DeactivateAlphaChannel,
1577 (void) CompositeImage(resize_image,resize_alpha,CopyAlphaCompositeOp,
1578 MagickTrue,0,0,exception);
1579 resize_alpha=DestroyImage(resize_alpha);
1581 (void) SetImageVirtualPixelMethod(resize_image,vp_save,exception);
1584 Clean up the results of the Distortion
1586 crop_area.width=columns;
1587 crop_area.height=rows;
1591 tmp_image=resize_image;
1592 resize_image=CropImage(tmp_image,&crop_area,exception);
1593 tmp_image=DestroyImage(tmp_image);
1594 return(resize_image);
1598 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1602 % D i s t o r t I m a g e %
1606 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1608 % DistortImage() distorts an image using various distortion methods, by
1609 % mapping color lookups of the source image to a new destination image
1610 % usally of the same size as the source image, unless 'bestfit' is set to
1613 % If 'bestfit' is enabled, and distortion allows it, the destination image is
1614 % adjusted to ensure the whole source 'image' will just fit within the final
1615 % destination image, which will be sized and offset accordingly. Also in
1616 % many cases the virtual offset of the source image will be taken into
1617 % account in the mapping.
1619 % If the '-verbose' control option has been set print to standard error the
1620 % equicelent '-fx' formula with coefficients for the function, if practical.
1622 % The format of the DistortImage() method is:
1624 % Image *DistortImage(const Image *image,const DistortImageMethod method,
1625 % const size_t number_arguments,const double *arguments,
1626 % MagickBooleanType bestfit, ExceptionInfo *exception)
1628 % A description of each parameter follows:
1630 % o image: the image to be distorted.
1632 % o method: the method of image distortion.
1634 % ArcDistortion always ignores source image offset, and always
1635 % 'bestfit' the destination image with the top left corner offset
1636 % relative to the polar mapping center.
1638 % Affine, Perspective, and Bilinear, do least squares fitting of the
1639 % distrotion when more than the minimum number of control point pairs
1642 % Perspective, and Bilinear, fall back to a Affine distortion when less
1643 % than 4 control point pairs are provided. While Affine distortions
1644 % let you use any number of control point pairs, that is Zero pairs is
1645 % a No-Op (viewport only) distortion, one pair is a translation and
1646 % two pairs of control points do a scale-rotate-translate, without any
1649 % o number_arguments: the number of arguments given.
1651 % o arguments: an array of floating point arguments for this method.
1653 % o bestfit: Attempt to 'bestfit' the size of the resulting image.
1654 % This also forces the resulting image to be a 'layered' virtual
1655 % canvas image. Can be overridden using 'distort:viewport' setting.
1657 % o exception: return any errors or warnings in this structure
1659 % Extra Controls from Image meta-data (artifacts)...
1662 % Output to stderr alternatives, internal coefficents, and FX
1663 % equivalents for the distortion operation (if feasible).
1664 % This forms an extra check of the distortion method, and allows users
1665 % access to the internal constants IM calculates for the distortion.
1667 % o "distort:viewport"
1668 % Directly set the output image canvas area and offest to use for the
1669 % resulting image, rather than use the original images canvas, or a
1670 % calculated 'bestfit' canvas.
1673 % Scale the size of the output canvas by this amount to provide a
1674 % method of Zooming, and for super-sampling the results.
1676 % Other settings that can effect results include
1678 % o 'interpolate' For source image lookups (scale enlargements)
1680 % o 'filter' Set filter to use for area-resampling (scale shrinking).
1681 % Set to 'point' to turn off and use 'interpolate' lookup
1685 MagickExport Image *DistortImage(const Image *image,DistortImageMethod method,
1686 const size_t number_arguments,const double *arguments,
1687 MagickBooleanType bestfit,ExceptionInfo *exception)
1689 #define DistortImageTag "Distort/Image"
1699 geometry; /* geometry of the distorted space viewport */
1704 assert(image != (Image *) NULL);
1705 assert(image->signature == MagickSignature);
1706 if (image->debug != MagickFalse)
1707 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1708 assert(exception != (ExceptionInfo *) NULL);
1709 assert(exception->signature == MagickSignature);
1712 Handle Special Compound Distortions
1714 if ( method == ResizeDistortion )
1716 if ( number_arguments != 2 )
1718 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1719 "InvalidArgument","%s : '%s'","Resize",
1720 "Invalid number of args: 2 only");
1721 return((Image *) NULL);
1723 distort_image=DistortResizeImage(image,(size_t)arguments[0],
1724 (size_t)arguments[1], exception);
1725 return(distort_image);
1729 Convert input arguments (usually as control points for reverse mapping)
1730 into mapping coefficients to apply the distortion.
1732 Note that some distortions are mapped to other distortions,
1733 and as such do not require specific code after this point.
1735 coeff = GenerateCoefficients(image, &method, number_arguments,
1736 arguments, 0, exception);
1737 if ( coeff == (double *) NULL )
1738 return((Image *) NULL);
1741 Determine the size and offset for a 'bestfit' destination.
1742 Usally the four corners of the source image is enough.
1745 /* default output image bounds, when no 'bestfit' is requested */
1746 geometry.width=image->columns;
1747 geometry.height=image->rows;
1751 if ( method == ArcDistortion ) {
1752 bestfit = MagickTrue; /* always calculate a 'best fit' viewport */
1755 /* Work out the 'best fit', (required for ArcDistortion) */
1758 s,d,min,max; /* source, dest coords --mapping--> min, max coords */
1761 fix_bounds = MagickTrue; /* enlarge bounds for VP handling */
1763 s.x=s.y=min.x=max.x=min.y=max.y=0.0; /* keep compiler happy */
1765 /* defines to figure out the bounds of the distorted image */
1766 #define InitalBounds(p) \
1768 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1769 min.x = max.x = p.x; \
1770 min.y = max.y = p.y; \
1772 #define ExpandBounds(p) \
1774 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1775 min.x = MagickMin(min.x,p.x); \
1776 max.x = MagickMax(max.x,p.x); \
1777 min.y = MagickMin(min.y,p.y); \
1778 max.y = MagickMax(max.y,p.y); \
1783 case AffineDistortion:
1784 { double inverse[6];
1785 InvertAffineCoefficients(coeff, inverse);
1786 s.x = (double) image->page.x;
1787 s.y = (double) image->page.y;
1788 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1789 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1791 s.x = (double) image->page.x+image->columns;
1792 s.y = (double) image->page.y;
1793 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1794 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1796 s.x = (double) image->page.x;
1797 s.y = (double) image->page.y+image->rows;
1798 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1799 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1801 s.x = (double) image->page.x+image->columns;
1802 s.y = (double) image->page.y+image->rows;
1803 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1804 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1808 case PerspectiveDistortion:
1809 { double inverse[8], scale;
1810 InvertPerspectiveCoefficients(coeff, inverse);
1811 s.x = (double) image->page.x;
1812 s.y = (double) image->page.y;
1813 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1814 scale=PerceptibleReciprocal(scale);
1815 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1816 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1818 s.x = (double) image->page.x+image->columns;
1819 s.y = (double) image->page.y;
1820 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1821 scale=PerceptibleReciprocal(scale);
1822 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1823 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1825 s.x = (double) image->page.x;
1826 s.y = (double) image->page.y+image->rows;
1827 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1828 scale=PerceptibleReciprocal(scale);
1829 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1830 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1832 s.x = (double) image->page.x+image->columns;
1833 s.y = (double) image->page.y+image->rows;
1834 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1835 scale=PerceptibleReciprocal(scale);
1836 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1837 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1843 /* Forward Map Corners */
1844 a = coeff[0]-coeff[1]/2; ca = cos(a); sa = sin(a);
1848 d.x = (coeff[2]-coeff[3])*ca;
1849 d.y = (coeff[2]-coeff[3])*sa;
1851 a = coeff[0]+coeff[1]/2; ca = cos(a); sa = sin(a);
1855 d.x = (coeff[2]-coeff[3])*ca;
1856 d.y = (coeff[2]-coeff[3])*sa;
1858 /* Orthogonal points along top of arc */
1859 for( a=(double) (ceil((double) ((coeff[0]-coeff[1]/2.0)/MagickPI2))*MagickPI2);
1860 a<(coeff[0]+coeff[1]/2.0); a+=MagickPI2 ) {
1861 ca = cos(a); sa = sin(a);
1867 Convert the angle_to_width and radius_to_height
1868 to appropriate scaling factors, to allow faster processing
1869 in the mapping function.
1871 coeff[1] = (double) (Magick2PI*image->columns/coeff[1]);
1872 coeff[3] = (double)image->rows/coeff[3];
1875 case PolarDistortion:
1877 if (number_arguments < 2)
1878 coeff[2] = coeff[3] = 0.0;
1879 min.x = coeff[2]-coeff[0];
1880 max.x = coeff[2]+coeff[0];
1881 min.y = coeff[3]-coeff[0];
1882 max.y = coeff[3]+coeff[0];
1883 /* should be about 1.0 if Rmin = 0 */
1884 coeff[7]=(double) geometry.height/(coeff[0]-coeff[1]);
1887 case DePolarDistortion:
1889 /* direct calculation as it needs to tile correctly
1890 * for reversibility in a DePolar-Polar cycle */
1891 fix_bounds = MagickFalse;
1892 geometry.x = geometry.y = 0;
1893 geometry.height = (size_t) ceil(coeff[0]-coeff[1]);
1894 geometry.width = (size_t)
1895 ceil((coeff[0]-coeff[1])*(coeff[5]-coeff[4])*0.5);
1896 /* correct scaling factors relative to new size */
1897 coeff[6]=(coeff[5]-coeff[4])/geometry.width; /* changed width */
1898 coeff[7]=(coeff[0]-coeff[1])/geometry.height; /* should be about 1.0 */
1901 case Cylinder2PlaneDistortion:
1903 /* direct calculation so center of distortion is either a pixel
1904 * center, or pixel edge. This allows for reversibility of the
1906 geometry.x = geometry.y = 0;
1907 geometry.width = (size_t) ceil( 2.0*coeff[1]*tan(coeff[0]/2.0) );
1908 geometry.height = (size_t) ceil( 2.0*coeff[3]/cos(coeff[0]/2.0) );
1909 /* correct center of distortion relative to new size */
1910 coeff[4] = (double) geometry.width/2.0;
1911 coeff[5] = (double) geometry.height/2.0;
1912 fix_bounds = MagickFalse;
1915 case Plane2CylinderDistortion:
1917 /* direct calculation center is either pixel center, or pixel edge
1918 * so as to allow reversibility of the image distortion */
1919 geometry.x = geometry.y = 0;
1920 geometry.width = (size_t) ceil(coeff[0]*coeff[1]); /* FOV * radius */
1921 geometry.height = (size_t) (2*coeff[3]); /* input image height */
1922 /* correct center of distortion relative to new size */
1923 coeff[4] = (double) geometry.width/2.0;
1924 coeff[5] = (double) geometry.height/2.0;
1925 fix_bounds = MagickFalse;
1928 case ShepardsDistortion:
1929 case BilinearForwardDistortion:
1930 case BilinearReverseDistortion:
1932 case QuadrilateralDistortion:
1934 case PolynomialDistortion:
1935 case BarrelDistortion:
1936 case BarrelInverseDistortion:
1938 /* no calculated bestfit available for these distortions */
1939 bestfit = MagickFalse;
1940 fix_bounds = MagickFalse;
1944 /* Set the output image geometry to calculated 'bestfit'.
1945 Yes this tends to 'over do' the file image size, ON PURPOSE!
1946 Do not do this for DePolar which needs to be exact for virtual tiling.
1949 geometry.x = (ssize_t) floor(min.x-0.5);
1950 geometry.y = (ssize_t) floor(min.y-0.5);
1951 geometry.width=(size_t) ceil(max.x-geometry.x+0.5);
1952 geometry.height=(size_t) ceil(max.y-geometry.y+0.5);
1955 } /* end bestfit destination image calculations */
1957 /* The user provided a 'viewport' expert option which may
1958 overrides some parts of the current output image geometry.
1959 This also overrides its default 'bestfit' setting.
1961 { const char *artifact=GetImageArtifact(image,"distort:viewport");
1962 viewport_given = MagickFalse;
1963 if ( artifact != (const char *) NULL ) {
1964 MagickStatusType flags=ParseAbsoluteGeometry(artifact,&geometry);
1966 (void) ThrowMagickException(exception,GetMagickModule(),
1967 OptionWarning,"InvalidSetting","'%s' '%s'",
1968 "distort:viewport",artifact);
1970 viewport_given = MagickTrue;
1974 /* Verbose output */
1975 if ( IfStringTrue(GetImageArtifact(image,"verbose")) ) {
1978 char image_gen[MaxTextExtent];
1981 /* Set destination image size and virtual offset */
1982 if ( bestfit || viewport_given ) {
1983 (void) FormatLocaleString(image_gen, MaxTextExtent," -size %.20gx%.20g "
1984 "-page %+.20g%+.20g xc: +insert \\\n",(double) geometry.width,
1985 (double) geometry.height,(double) geometry.x,(double) geometry.y);
1986 lookup="v.p{ xx-v.page.x-.5, yy-v.page.y-.5 }";
1989 image_gen[0] = '\0'; /* no destination to generate */
1990 lookup = "p{ xx-page.x-.5, yy-page.y-.5 }"; /* simplify lookup */
1994 case AffineDistortion:
1998 inverse = (double *) AcquireQuantumMemory(6,sizeof(*inverse));
1999 if (inverse == (double *) NULL) {
2000 coeff = (double *) RelinquishMagickMemory(coeff);
2001 (void) ThrowMagickException(exception,GetMagickModule(),
2002 ResourceLimitError,"MemoryAllocationFailed",
2003 "%s", "DistortImages");
2004 return((Image *) NULL);
2006 InvertAffineCoefficients(coeff, inverse);
2007 CoefficientsToAffineArgs(inverse);
2008 (void) FormatLocaleFile(stderr, "Affine Projection:\n");
2009 (void) FormatLocaleFile(stderr, " -distort AffineProjection \\\n '");
2010 for (i=0; i < 5; i++)
2011 (void) FormatLocaleFile(stderr, "%lf,", inverse[i]);
2012 (void) FormatLocaleFile(stderr, "%lf'\n", inverse[5]);
2013 inverse = (double *) RelinquishMagickMemory(inverse);
2015 (void) FormatLocaleFile(stderr, "Affine Distort, FX Equivelent:\n");
2016 (void) FormatLocaleFile(stderr, "%s", image_gen);
2017 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2018 (void) FormatLocaleFile(stderr, " xx=%+lf*ii %+lf*jj %+lf;\n",
2019 coeff[0], coeff[1], coeff[2]);
2020 (void) FormatLocaleFile(stderr, " yy=%+lf*ii %+lf*jj %+lf;\n",
2021 coeff[3], coeff[4], coeff[5]);
2022 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2027 case PerspectiveDistortion:
2031 inverse = (double *) AcquireQuantumMemory(8,sizeof(*inverse));
2032 if (inverse == (double *) NULL) {
2033 coeff = (double *) RelinquishMagickMemory(coeff);
2034 (void) ThrowMagickException(exception,GetMagickModule(),
2035 ResourceLimitError,"MemoryAllocationFailed",
2036 "%s", "DistortCoefficients");
2037 return((Image *) NULL);
2039 InvertPerspectiveCoefficients(coeff, inverse);
2040 (void) FormatLocaleFile(stderr, "Perspective Projection:\n");
2041 (void) FormatLocaleFile(stderr, " -distort PerspectiveProjection \\\n '");
2043 (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
2044 (void) FormatLocaleFile(stderr, "\n ");
2046 (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
2047 (void) FormatLocaleFile(stderr, "%lf'\n", inverse[7]);
2048 inverse = (double *) RelinquishMagickMemory(inverse);
2050 (void) FormatLocaleFile(stderr, "Perspective Distort, FX Equivelent:\n");
2051 (void) FormatLocaleFile(stderr, "%s", image_gen);
2052 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2053 (void) FormatLocaleFile(stderr, " rr=%+lf*ii %+lf*jj + 1;\n",
2054 coeff[6], coeff[7]);
2055 (void) FormatLocaleFile(stderr, " xx=(%+lf*ii %+lf*jj %+lf)/rr;\n",
2056 coeff[0], coeff[1], coeff[2]);
2057 (void) FormatLocaleFile(stderr, " yy=(%+lf*ii %+lf*jj %+lf)/rr;\n",
2058 coeff[3], coeff[4], coeff[5]);
2059 (void) FormatLocaleFile(stderr, " rr%s0 ? %s : blue' \\\n",
2060 coeff[8] < 0 ? "<" : ">", lookup);
2064 case BilinearForwardDistortion:
2065 (void) FormatLocaleFile(stderr, "BilinearForward Mapping Equations:\n");
2066 (void) FormatLocaleFile(stderr, "%s", image_gen);
2067 (void) FormatLocaleFile(stderr, " i = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
2068 coeff[0], coeff[1], coeff[2], coeff[3]);
2069 (void) FormatLocaleFile(stderr, " j = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
2070 coeff[4], coeff[5], coeff[6], coeff[7]);
2073 (void) FormatLocaleFile(stderr, " c8 = %+lf c9 = 2*a = %+lf;\n",
2074 coeff[8], coeff[9]);
2076 (void) FormatLocaleFile(stderr, "BilinearForward Distort, FX Equivelent:\n");
2077 (void) FormatLocaleFile(stderr, "%s", image_gen);
2078 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
2079 0.5-coeff[3], 0.5-coeff[7]);
2080 (void) FormatLocaleFile(stderr, " bb=%lf*ii %+lf*jj %+lf;\n",
2081 coeff[6], -coeff[2], coeff[8]);
2082 /* Handle Special degenerate (non-quadratic) or trapezoidal case */
2083 if ( coeff[9] != 0 ) {
2084 (void) FormatLocaleFile(stderr, " rt=bb*bb %+lf*(%lf*ii%+lf*jj);\n",
2085 -2*coeff[9], coeff[4], -coeff[0]);
2086 (void) FormatLocaleFile(stderr, " yy=( -bb + sqrt(rt) ) / %lf;\n",
2089 (void) FormatLocaleFile(stderr, " yy=(%lf*ii%+lf*jj)/bb;\n",
2090 -coeff[4], coeff[0]);
2091 (void) FormatLocaleFile(stderr, " xx=(ii %+lf*yy)/(%lf %+lf*yy);\n",
2092 -coeff[1], coeff[0], coeff[2]);
2093 if ( coeff[9] != 0 )
2094 (void) FormatLocaleFile(stderr, " (rt < 0 ) ? red : %s'\n", lookup);
2096 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2099 case BilinearReverseDistortion:
2101 (void) FormatLocaleFile(stderr, "Polynomial Projection Distort:\n");
2102 (void) FormatLocaleFile(stderr, " -distort PolynomialProjection \\\n");
2103 (void) FormatLocaleFile(stderr, " '1.5, %lf, %lf, %lf, %lf,\n",
2104 coeff[3], coeff[0], coeff[1], coeff[2]);
2105 (void) FormatLocaleFile(stderr, " %lf, %lf, %lf, %lf'\n",
2106 coeff[7], coeff[4], coeff[5], coeff[6]);
2108 (void) FormatLocaleFile(stderr, "BilinearReverse Distort, FX Equivelent:\n");
2109 (void) FormatLocaleFile(stderr, "%s", image_gen);
2110 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2111 (void) FormatLocaleFile(stderr, " xx=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
2112 coeff[0], coeff[1], coeff[2], coeff[3]);
2113 (void) FormatLocaleFile(stderr, " yy=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
2114 coeff[4], coeff[5], coeff[6], coeff[7]);
2115 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2118 case PolynomialDistortion:
2120 size_t nterms = (size_t) coeff[1];
2121 (void) FormatLocaleFile(stderr, "Polynomial (order %lg, terms %lu), FX Equivelent\n",
2122 coeff[0],(unsigned long) nterms);
2123 (void) FormatLocaleFile(stderr, "%s", image_gen);
2124 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2125 (void) FormatLocaleFile(stderr, " xx =");
2126 for (i=0; i<(ssize_t) nterms; i++) {
2127 if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n ");
2128 (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i],
2131 (void) FormatLocaleFile(stderr, ";\n yy =");
2132 for (i=0; i<(ssize_t) nterms; i++) {
2133 if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n ");
2134 (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i+nterms],
2137 (void) FormatLocaleFile(stderr, ";\n %s' \\\n", lookup);
2142 (void) FormatLocaleFile(stderr, "Arc Distort, Internal Coefficients:\n");
2143 for ( i=0; i<5; i++ )
2144 (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
2145 (void) FormatLocaleFile(stderr, "Arc Distort, FX Equivelent:\n");
2146 (void) FormatLocaleFile(stderr, "%s", image_gen);
2147 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x; jj=j+page.y;\n");
2148 (void) FormatLocaleFile(stderr, " xx=(atan2(jj,ii)%+lf)/(2*pi);\n",
2150 (void) FormatLocaleFile(stderr, " xx=xx-round(xx);\n");
2151 (void) FormatLocaleFile(stderr, " xx=xx*%lf %+lf;\n",
2152 coeff[1], coeff[4]);
2153 (void) FormatLocaleFile(stderr, " yy=(%lf - hypot(ii,jj)) * %lf;\n",
2154 coeff[2], coeff[3]);
2155 (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n");
2158 case PolarDistortion:
2160 (void) FormatLocaleFile(stderr, "Polar Distort, Internal Coefficents\n");
2161 for ( i=0; i<8; i++ )
2162 (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
2163 (void) FormatLocaleFile(stderr, "Polar Distort, FX Equivelent:\n");
2164 (void) FormatLocaleFile(stderr, "%s", image_gen);
2165 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
2166 -coeff[2], -coeff[3]);
2167 (void) FormatLocaleFile(stderr, " xx=(atan2(ii,jj)%+lf)/(2*pi);\n",
2168 -(coeff[4]+coeff[5])/2 );
2169 (void) FormatLocaleFile(stderr, " xx=xx-round(xx);\n");
2170 (void) FormatLocaleFile(stderr, " xx=xx*2*pi*%lf + v.w/2;\n",
2172 (void) FormatLocaleFile(stderr, " yy=(hypot(ii,jj)%+lf)*%lf;\n",
2173 -coeff[1], coeff[7] );
2174 (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n");
2177 case DePolarDistortion:
2179 (void) FormatLocaleFile(stderr, "DePolar Distort, Internal Coefficents\n");
2180 for ( i=0; i<8; i++ )
2181 (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
2182 (void) FormatLocaleFile(stderr, "DePolar Distort, FX Equivelent:\n");
2183 (void) FormatLocaleFile(stderr, "%s", image_gen);
2184 (void) FormatLocaleFile(stderr, " -fx 'aa=(i+.5)*%lf %+lf;\n", coeff[6], +coeff[4] );
2185 (void) FormatLocaleFile(stderr, " rr=(j+.5)*%lf %+lf;\n", coeff[7], +coeff[1] );
2186 (void) FormatLocaleFile(stderr, " xx=rr*sin(aa) %+lf;\n", coeff[2] );
2187 (void) FormatLocaleFile(stderr, " yy=rr*cos(aa) %+lf;\n", coeff[3] );
2188 (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n");
2191 case Cylinder2PlaneDistortion:
2193 (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, Internal Coefficents\n");
2194 (void) FormatLocaleFile(stderr, " cylinder_radius = %+lf\n", coeff[1]);
2195 (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, FX Equivelent:\n");
2196 (void) FormatLocaleFile(stderr, "%s", image_gen);
2197 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
2198 -coeff[4], -coeff[5]);
2199 (void) FormatLocaleFile(stderr, " aa=atan(ii/%+lf);\n", coeff[1] );
2200 (void) FormatLocaleFile(stderr, " xx=%lf*aa%+lf;\n",
2201 coeff[1], coeff[2] );
2202 (void) FormatLocaleFile(stderr, " yy=jj*cos(aa)%+lf;\n", coeff[3] );
2203 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2206 case Plane2CylinderDistortion:
2208 (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, Internal Coefficents\n");
2209 (void) FormatLocaleFile(stderr, " cylinder_radius = %+lf\n", coeff[1]);
2210 (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, FX Equivelent:\n");
2211 (void) FormatLocaleFile(stderr, "%s", image_gen);
2212 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
2213 -coeff[4], -coeff[5]);
2214 (void) FormatLocaleFile(stderr, " ii=ii/%+lf;\n", coeff[1] );
2215 (void) FormatLocaleFile(stderr, " xx=%lf*tan(ii)%+lf;\n",
2216 coeff[1], coeff[2] );
2217 (void) FormatLocaleFile(stderr, " yy=jj/cos(ii)%+lf;\n",
2219 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2222 case BarrelDistortion:
2223 case BarrelInverseDistortion:
2225 /* NOTE: This does the barrel roll in pixel coords not image coords
2226 ** The internal distortion must do it in image coordinates,
2227 ** so that is what the center coeff (8,9) is given in.
2229 xc = ((double)image->columns-1.0)/2.0 + image->page.x;
2230 yc = ((double)image->rows-1.0)/2.0 + image->page.y;
2231 (void) FormatLocaleFile(stderr, "Barrel%s Distort, FX Equivelent:\n",
2232 method == BarrelDistortion ? "" : "Inv");
2233 (void) FormatLocaleFile(stderr, "%s", image_gen);
2234 if ( fabs(coeff[8]-xc-0.5) < 0.1 && fabs(coeff[9]-yc-0.5) < 0.1 )
2235 (void) FormatLocaleFile(stderr, " -fx 'xc=(w-1)/2; yc=(h-1)/2;\n");
2237 (void) FormatLocaleFile(stderr, " -fx 'xc=%lf; yc=%lf;\n",
2238 coeff[8]-0.5, coeff[9]-0.5);
2239 (void) FormatLocaleFile(stderr,
2240 " ii=i-xc; jj=j-yc; rr=hypot(ii,jj);\n");
2241 (void) FormatLocaleFile(stderr, " ii=ii%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2242 method == BarrelDistortion ? "*" : "/",
2243 coeff[0],coeff[1],coeff[2],coeff[3]);
2244 (void) FormatLocaleFile(stderr, " jj=jj%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2245 method == BarrelDistortion ? "*" : "/",
2246 coeff[4],coeff[5],coeff[6],coeff[7]);
2247 (void) FormatLocaleFile(stderr, " v.p{fx*ii+xc,fy*jj+yc}' \\\n");
2254 /* The user provided a 'scale' expert option will scale the
2255 output image size, by the factor given allowing for super-sampling
2256 of the distorted image space. Any scaling factors must naturally
2257 be halved as a result.
2259 { const char *artifact;
2260 artifact=GetImageArtifact(image,"distort:scale");
2261 output_scaling = 1.0;
2262 if (artifact != (const char *) NULL) {
2263 output_scaling = fabs(StringToDouble(artifact,(char **) NULL));
2264 geometry.width=(size_t) (output_scaling*geometry.width+0.5);
2265 geometry.height=(size_t) (output_scaling*geometry.height+0.5);
2266 geometry.x=(ssize_t) (output_scaling*geometry.x+0.5);
2267 geometry.y=(ssize_t) (output_scaling*geometry.y+0.5);
2268 if ( output_scaling < 0.1 ) {
2269 coeff = (double *) RelinquishMagickMemory(coeff);
2270 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
2271 "InvalidArgument","%s", "-set option:distort:scale" );
2272 return((Image *) NULL);
2274 output_scaling = 1/output_scaling;
2277 #define ScaleFilter(F,A,B,C,D) \
2278 ScaleResampleFilter( (F), \
2279 output_scaling*(A), output_scaling*(B), \
2280 output_scaling*(C), output_scaling*(D) )
2283 Initialize the distort image attributes.
2285 distort_image=CloneImage(image,geometry.width,geometry.height,MagickTrue,
2287 if (distort_image == (Image *) NULL)
2288 return((Image *) NULL);
2289 /* if image is ColorMapped - change it to DirectClass */
2290 if (SetImageStorageClass(distort_image,DirectClass,exception) == MagickFalse)
2292 distort_image=DestroyImage(distort_image);
2293 return((Image *) NULL);
2295 if ((IsPixelInfoGray(&distort_image->background_color) == MagickFalse) &&
2296 (IsGrayColorspace(distort_image->colorspace) != MagickFalse))
2297 (void) SetImageColorspace(distort_image,sRGBColorspace,exception);
2298 if (distort_image->background_color.alpha_trait == BlendPixelTrait)
2299 distort_image->alpha_trait=BlendPixelTrait;
2300 distort_image->page.x=geometry.x;
2301 distort_image->page.y=geometry.y;
2303 { /* ----- MAIN CODE -----
2304 Sample the source image to each pixel in the distort image.
2319 **restrict resample_filter;
2326 GetPixelInfo(distort_image,&zero);
2327 resample_filter=AcquireResampleFilterThreadSet(image,
2328 UndefinedVirtualPixelMethod,MagickFalse,exception);
2329 distort_view=AcquireAuthenticCacheView(distort_image,exception);
2330 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2331 #pragma omp parallel for schedule(static,4) shared(progress,status) \
2332 magick_threads(image,distort_image,distort_image->rows,1)
2334 for (j=0; j < (ssize_t) distort_image->rows; j++)
2337 id = GetOpenMPThreadId();
2340 validity; /* how mathematically valid is this the mapping */
2346 pixel, /* pixel color to assign to distorted image */
2347 invalid; /* the color to assign when distort result is invalid */
2351 s; /* transform destination image x,y to source image x,y */
2359 q=QueueCacheViewAuthenticPixels(distort_view,0,j,distort_image->columns,1,
2361 if (q == (Quantum *) NULL)
2368 /* Define constant scaling vectors for Affine Distortions
2369 Other methods are either variable, or use interpolated lookup
2373 case AffineDistortion:
2374 ScaleFilter( resample_filter[id],
2376 coeff[3], coeff[4] );
2382 /* Initialize default pixel validity
2383 * negative: pixel is invalid output 'matte_color'
2384 * 0.0 to 1.0: antialiased, mix with resample output
2385 * 1.0 or greater: use resampled output.
2389 invalid=distort_image->matte_color;
2390 if (distort_image->colorspace == CMYKColorspace)
2391 ConvertRGBToCMYK(&invalid); /* what about other color spaces? */
2392 for (i=0; i < (ssize_t) distort_image->columns; i++)
2394 /* map pixel coordinate to distortion space coordinate */
2395 d.x = (double) (geometry.x+i+0.5)*output_scaling;
2396 d.y = (double) (geometry.y+j+0.5)*output_scaling;
2397 s = d; /* default is a no-op mapping */
2400 case AffineDistortion:
2402 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2403 s.y=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2404 /* Affine partial derivitives are constant -- set above */
2407 case PerspectiveDistortion:
2410 p,q,r,abs_r,abs_c6,abs_c7,scale;
2411 /* perspective is a ratio of affines */
2412 p=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2413 q=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2414 r=coeff[6]*d.x+coeff[7]*d.y+1.0;
2415 /* Pixel Validity -- is it a 'sky' or 'ground' pixel */
2416 validity = (r*coeff[8] < 0.0) ? 0.0 : 1.0;
2417 /* Determine horizon anti-alias blending */
2419 abs_c6 = fabs(coeff[6]);
2420 abs_c7 = fabs(coeff[7]);
2421 if ( abs_c6 > abs_c7 ) {
2422 if ( abs_r < abs_c6*output_scaling )
2423 validity = 0.5 - coeff[8]*r/(coeff[6]*output_scaling);
2425 else if ( abs_r < abs_c7*output_scaling )
2426 validity = 0.5 - coeff[8]*r/(coeff[7]*output_scaling);
2427 /* Perspective Sampling Point (if valid) */
2428 if ( validity > 0.0 ) {
2429 /* divide by r affine, for perspective scaling */
2433 /* Perspective Partial Derivatives or Scaling Vectors */
2435 ScaleFilter( resample_filter[id],
2436 (r*coeff[0] - p*coeff[6])*scale,
2437 (r*coeff[1] - p*coeff[7])*scale,
2438 (r*coeff[3] - q*coeff[6])*scale,
2439 (r*coeff[4] - q*coeff[7])*scale );
2443 case BilinearReverseDistortion:
2445 /* Reversed Mapped is just a simple polynomial */
2446 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2]*d.x*d.y+coeff[3];
2447 s.y=coeff[4]*d.x+coeff[5]*d.y
2448 +coeff[6]*d.x*d.y+coeff[7];
2449 /* Bilinear partial derivitives of scaling vectors */
2450 ScaleFilter( resample_filter[id],
2451 coeff[0] + coeff[2]*d.y,
2452 coeff[1] + coeff[2]*d.x,
2453 coeff[4] + coeff[6]*d.y,
2454 coeff[5] + coeff[6]*d.x );
2457 case BilinearForwardDistortion:
2459 /* Forward mapped needs reversed polynomial equations
2460 * which unfortunatally requires a square root! */
2462 d.x -= coeff[3]; d.y -= coeff[7];
2463 b = coeff[6]*d.x - coeff[2]*d.y + coeff[8];
2464 c = coeff[4]*d.x - coeff[0]*d.y;
2467 /* Handle Special degenerate (non-quadratic) case
2468 * Currently without horizon anti-alising */
2469 if ( fabs(coeff[9]) < MagickEpsilon )
2472 c = b*b - 2*coeff[9]*c;
2476 s.y = ( -b + sqrt(c) )/coeff[9];
2478 if ( validity > 0.0 )
2479 s.x = ( d.x - coeff[1]*s.y) / ( coeff[0] + coeff[2]*s.y );
2481 /* NOTE: the sign of the square root should be -ve for parts
2482 where the source image becomes 'flipped' or 'mirrored'.
2483 FUTURE: Horizon handling
2484 FUTURE: Scaling factors or Deritives (how?)
2489 case BilinearDistortion:
2490 /* Bilinear mapping of any Quadrilateral to any Quadrilateral */
2491 /* UNDER DEVELOPMENT */
2494 case PolynomialDistortion:
2496 /* multi-ordered polynomial */
2501 nterms=(ssize_t)coeff[1];
2504 du,dv; /* the du,dv vectors from unit dx,dy -- derivatives */
2506 s.x=s.y=du.x=du.y=dv.x=dv.y=0.0;
2507 for(k=0; k < nterms; k++) {
2508 s.x += poly_basis_fn(k,d.x,d.y)*coeff[2+k];
2509 du.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k];
2510 du.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k];
2511 s.y += poly_basis_fn(k,d.x,d.y)*coeff[2+k+nterms];
2512 dv.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k+nterms];
2513 dv.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k+nterms];
2515 ScaleFilter( resample_filter[id], du.x,du.y,dv.x,dv.y );
2520 /* what is the angle and radius in the destination image */
2521 s.x = (double) ((atan2(d.y,d.x) - coeff[0])/Magick2PI);
2522 s.x -= MagickRound(s.x); /* angle */
2523 s.y = hypot(d.x,d.y); /* radius */
2525 /* Arc Distortion Partial Scaling Vectors
2526 Are derived by mapping the perpendicular unit vectors
2527 dR and dA*R*2PI rather than trying to map dx and dy
2528 The results is a very simple orthogonal aligned ellipse.
2530 if ( s.y > MagickEpsilon )
2531 ScaleFilter( resample_filter[id],
2532 (double) (coeff[1]/(Magick2PI*s.y)), 0, 0, coeff[3] );
2534 ScaleFilter( resample_filter[id],
2535 distort_image->columns*2, 0, 0, coeff[3] );
2537 /* now scale the angle and radius for source image lookup point */
2538 s.x = s.x*coeff[1] + coeff[4] + image->page.x +0.5;
2539 s.y = (coeff[2] - s.y) * coeff[3] + image->page.y;
2542 case PolarDistortion:
2543 { /* 2D Cartesain to Polar View */
2546 s.x = atan2(d.x,d.y) - (coeff[4]+coeff[5])/2;
2548 s.x -= MagickRound(s.x);
2549 s.x *= Magick2PI; /* angle - relative to centerline */
2550 s.y = hypot(d.x,d.y); /* radius */
2552 /* Polar Scaling vectors are based on mapping dR and dA vectors
2553 This results in very simple orthogonal scaling vectors
2555 if ( s.y > MagickEpsilon )
2556 ScaleFilter( resample_filter[id],
2557 (double) (coeff[6]/(Magick2PI*s.y)), 0, 0, coeff[7] );
2559 ScaleFilter( resample_filter[id],
2560 distort_image->columns*2, 0, 0, coeff[7] );
2562 /* now finish mapping radius/angle to source x,y coords */
2563 s.x = s.x*coeff[6] + (double)image->columns/2.0 + image->page.x;
2564 s.y = (s.y-coeff[1])*coeff[7] + image->page.y;
2567 case DePolarDistortion:
2568 { /* @D Polar to Carteasain */
2569 /* ignore all destination virtual offsets */
2570 d.x = ((double)i+0.5)*output_scaling*coeff[6]+coeff[4];
2571 d.y = ((double)j+0.5)*output_scaling*coeff[7]+coeff[1];
2572 s.x = d.y*sin(d.x) + coeff[2];
2573 s.y = d.y*cos(d.x) + coeff[3];
2574 /* derivatives are usless - better to use SuperSampling */
2577 case Cylinder2PlaneDistortion:
2578 { /* 3D Cylinder to Tangential Plane */
2580 /* relative to center of distortion */
2581 d.x -= coeff[4]; d.y -= coeff[5];
2582 d.x /= coeff[1]; /* x' = x/r */
2583 ax=atan(d.x); /* aa = atan(x/r) = u/r */
2584 cx=cos(ax); /* cx = cos(atan(x/r)) = 1/sqrt(x^2+u^2) */
2585 s.x = coeff[1]*ax; /* u = r*atan(x/r) */
2586 s.y = d.y*cx; /* v = y*cos(u/r) */
2587 /* derivatives... (see personnal notes) */
2588 ScaleFilter( resample_filter[id],
2589 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2591 if ( i == 0 && j == 0 ) {
2592 fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2593 fprintf(stderr, "phi = %lf\n", (double)(ax * 180.0/MagickPI) );
2594 fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n",
2595 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2598 /* add center of distortion in source */
2599 s.x += coeff[2]; s.y += coeff[3];
2602 case Plane2CylinderDistortion:
2603 { /* 3D Cylinder to Tangential Plane */
2604 /* relative to center of distortion */
2605 d.x -= coeff[4]; d.y -= coeff[5];
2607 /* is pixel valid - horizon of a infinite Virtual-Pixel Plane
2608 * (see Anthony Thyssen's personal note) */
2609 validity = (double) (coeff[1]*MagickPI2 - fabs(d.x))/output_scaling + 0.5;
2611 if ( validity > 0.0 ) {
2613 d.x /= coeff[1]; /* x'= x/r */
2614 cx = 1/cos(d.x); /* cx = 1/cos(x/r) */
2615 tx = tan(d.x); /* tx = tan(x/r) */
2616 s.x = coeff[1]*tx; /* u = r * tan(x/r) */
2617 s.y = d.y*cx; /* v = y / cos(x/r) */
2618 /* derivatives... (see Anthony Thyssen's personal notes) */
2619 ScaleFilter( resample_filter[id],
2620 cx*cx, 0.0, s.y*cx/coeff[1], cx );
2622 /*if ( i == 0 && j == 0 )*/
2623 if ( d.x == 0.5 && d.y == 0.5 ) {
2624 fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2625 fprintf(stderr, "radius = %lf phi = %lf validity = %lf\n",
2626 coeff[1], (double)(d.x * 180.0/MagickPI), validity );
2627 fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n",
2628 cx*cx, 0.0, s.y*cx/coeff[1], cx);
2632 /* add center of distortion in source */
2633 s.x += coeff[2]; s.y += coeff[3];
2636 case BarrelDistortion:
2637 case BarrelInverseDistortion:
2638 { /* Lens Barrel Distionion Correction */
2639 double r,fx,fy,gx,gy;
2640 /* Radial Polynomial Distortion (de-normalized) */
2643 r = sqrt(d.x*d.x+d.y*d.y);
2644 if ( r > MagickEpsilon ) {
2645 fx = ((coeff[0]*r + coeff[1])*r + coeff[2])*r + coeff[3];
2646 fy = ((coeff[4]*r + coeff[5])*r + coeff[6])*r + coeff[7];
2647 gx = ((3*coeff[0]*r + 2*coeff[1])*r + coeff[2])/r;
2648 gy = ((3*coeff[4]*r + 2*coeff[5])*r + coeff[6])/r;
2649 /* adjust functions and scaling for 'inverse' form */
2650 if ( method == BarrelInverseDistortion ) {
2651 fx = 1/fx; fy = 1/fy;
2652 gx *= -fx*fx; gy *= -fy*fy;
2654 /* Set the source pixel to lookup and EWA derivative vectors */
2655 s.x = d.x*fx + coeff[8];
2656 s.y = d.y*fy + coeff[9];
2657 ScaleFilter( resample_filter[id],
2658 gx*d.x*d.x + fx, gx*d.x*d.y,
2659 gy*d.x*d.y, gy*d.y*d.y + fy );
2662 /* Special handling to avoid divide by zero when r==0
2664 ** The source and destination pixels match in this case
2665 ** which was set at the top of the loop using s = d;
2666 ** otherwise... s.x=coeff[8]; s.y=coeff[9];
2668 if ( method == BarrelDistortion )
2669 ScaleFilter( resample_filter[id],
2670 coeff[3], 0, 0, coeff[7] );
2671 else /* method == BarrelInverseDistortion */
2672 /* FUTURE, trap for D==0 causing division by zero */
2673 ScaleFilter( resample_filter[id],
2674 1.0/coeff[3], 0, 0, 1.0/coeff[7] );
2678 case ShepardsDistortion:
2679 { /* Shepards Method, or Inverse Weighted Distance for
2680 displacement around the destination image control points
2681 The input arguments are the coefficents to the function.
2682 This is more of a 'displacement' function rather than an
2683 absolute distortion function.
2685 Note: We can not determine derivatives using shepards method
2686 so only a point sample interpolatation can be used.
2693 denominator = s.x = s.y = 0;
2694 for(i=0; i<number_arguments; i+=4) {
2696 ((double)d.x-arguments[i+2])*((double)d.x-arguments[i+2])
2697 + ((double)d.y-arguments[i+3])*((double)d.y-arguments[i+3]);
2698 weight = pow(weight,coeff[0]); /* shepards power factor */
2699 weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
2701 s.x += (arguments[ i ]-arguments[i+2])*weight;
2702 s.y += (arguments[i+1]-arguments[i+3])*weight;
2703 denominator += weight;
2707 s.x += d.x; /* make it as relative displacement */
2712 break; /* use the default no-op given above */
2714 /* map virtual canvas location back to real image coordinate */
2715 if ( bestfit && method != ArcDistortion ) {
2716 s.x -= image->page.x;
2717 s.y -= image->page.y;
2722 if ( validity <= 0.0 ) {
2723 /* result of distortion is an invalid pixel - don't resample */
2724 SetPixelInfoPixel(distort_image,&invalid,q);
2727 /* resample the source image to find its correct color */
2728 (void) ResamplePixelColor(resample_filter[id],s.x,s.y,&pixel,
2730 /* if validity between 0.0 and 1.0 mix result with invalid pixel */
2731 if ( validity < 1.0 ) {
2732 /* Do a blend of sample color and invalid pixel */
2733 /* should this be a 'Blend', or an 'Over' compose */
2734 CompositePixelInfoBlend(&pixel,validity,&invalid,(1.0-validity),
2737 SetPixelInfoPixel(distort_image,&pixel,q);
2739 q+=GetPixelChannels(distort_image);
2741 sync=SyncCacheViewAuthenticPixels(distort_view,exception);
2742 if (sync == MagickFalse)
2744 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2749 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2750 #pragma omp critical (MagickCore_DistortImage)
2752 proceed=SetImageProgress(image,DistortImageTag,progress++,
2754 if (proceed == MagickFalse)
2758 distort_view=DestroyCacheView(distort_view);
2759 resample_filter=DestroyResampleFilterThreadSet(resample_filter);
2761 if (status == MagickFalse)
2762 distort_image=DestroyImage(distort_image);
2765 /* Arc does not return an offset unless 'bestfit' is in effect
2766 And the user has not provided an overriding 'viewport'.
2768 if ( method == ArcDistortion && !bestfit && !viewport_given ) {
2769 distort_image->page.x = 0;
2770 distort_image->page.y = 0;
2772 coeff = (double *) RelinquishMagickMemory(coeff);
2773 return(distort_image);
2777 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2781 % R o t a t e I m a g e %
2785 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2787 % RotateImage() creates a new image that is a rotated copy of an existing
2788 % one. Positive angles rotate counter-clockwise (right-hand rule), while
2789 % negative angles rotate clockwise. Rotated images are usually larger than
2790 % the originals and have 'empty' triangular corners. X axis. Empty
2791 % triangles left over from shearing the image are filled with the background
2792 % color defined by member 'background_color' of the image. RotateImage
2793 % allocates the memory necessary for the new Image structure and returns a
2794 % pointer to the new image.
2796 % The format of the RotateImage method is:
2798 % Image *RotateImage(const Image *image,const double degrees,
2799 % ExceptionInfo *exception)
2801 % A description of each parameter follows.
2803 % o image: the image.
2805 % o degrees: Specifies the number of degrees to rotate the image.
2807 % o exception: return any errors or warnings in this structure.
2810 MagickExport Image *RotateImage(const Image *image,const double degrees,
2811 ExceptionInfo *exception)
2827 Adjust rotation angle.
2829 assert(image != (Image *) NULL);
2830 assert(image->signature == MagickSignature);
2831 if (image->debug != MagickFalse)
2832 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2833 assert(exception != (ExceptionInfo *) NULL);
2834 assert(exception->signature == MagickSignature);
2836 while (angle < -45.0)
2838 for (rotations=0; angle > 45.0; rotations++)
2841 shear.x=(-tan((double) DegreesToRadians(angle)/2.0));
2842 shear.y=sin((double) DegreesToRadians(angle));
2843 if ((fabs(shear.x) < MagickEpsilon) && (fabs(shear.y) < MagickEpsilon))
2844 return(IntegralRotateImage(image,rotations,exception));
2845 distort_image=CloneImage(image,0,0,MagickTrue,exception);
2846 if (distort_image == (Image *) NULL)
2847 return((Image *) NULL);
2848 (void) SetImageVirtualPixelMethod(distort_image,BackgroundVirtualPixelMethod,
2850 rotate_image=DistortImage(distort_image,ScaleRotateTranslateDistortion,1,
2851 °rees,MagickTrue,exception);
2852 distort_image=DestroyImage(distort_image);
2853 return(rotate_image);
2857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2861 % S p a r s e C o l o r I m a g e %
2865 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2867 % SparseColorImage(), given a set of coordinates, interpolates the colors
2868 % found at those coordinates, across the whole image, using various methods.
2870 % The format of the SparseColorImage() method is:
2872 % Image *SparseColorImage(const Image *image,
2873 % const SparseColorMethod method,const size_t number_arguments,
2874 % const double *arguments,ExceptionInfo *exception)
2876 % A description of each parameter follows:
2878 % o image: the image to be filled in.
2880 % o method: the method to fill in the gradient between the control points.
2882 % The methods used for SparseColor() are often simular to methods
2883 % used for DistortImage(), and even share the same code for determination
2884 % of the function coefficents, though with more dimensions (or resulting
2887 % o number_arguments: the number of arguments given.
2889 % o arguments: array of floating point arguments for this method--
2890 % x,y,color_values-- with color_values given as normalized values.
2892 % o exception: return any errors or warnings in this structure
2895 MagickExport Image *SparseColorImage(const Image *image,
2896 const SparseColorMethod method,const size_t number_arguments,
2897 const double *arguments,ExceptionInfo *exception)
2899 #define SparseColorTag "Distort/SparseColor"
2913 assert(image != (Image *) NULL);
2914 assert(image->signature == MagickSignature);
2915 if (image->debug != MagickFalse)
2916 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2917 assert(exception != (ExceptionInfo *) NULL);
2918 assert(exception->signature == MagickSignature);
2920 /* Determine number of color values needed per control point */
2922 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2924 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2926 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2928 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2929 (image->colorspace == CMYKColorspace))
2931 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2932 (image->alpha_trait == BlendPixelTrait))
2936 Convert input arguments into mapping coefficients, this this case
2937 we are mapping (distorting) colors, rather than coordinates.
2939 { DistortImageMethod
2942 distort_method=(DistortImageMethod) method;
2943 if ( distort_method >= SentinelDistortion )
2944 distort_method = ShepardsDistortion; /* Pretend to be Shepards */
2945 coeff = GenerateCoefficients(image, &distort_method, number_arguments,
2946 arguments, number_colors, exception);
2947 if ( coeff == (double *) NULL )
2948 return((Image *) NULL);
2950 Note some Distort Methods may fall back to other simpler methods,
2951 Currently the only fallback of concern is Bilinear to Affine
2952 (Barycentric), which is alaso sparse_colr method. This also ensures
2953 correct two and one color Barycentric handling.
2955 sparse_method = (SparseColorMethod) distort_method;
2956 if ( distort_method == ShepardsDistortion )
2957 sparse_method = method; /* return non-distort methods to normal */
2958 if ( sparse_method == InverseColorInterpolate )
2959 coeff[0]=0.5; /* sqrt() the squared distance for inverse */
2962 /* Verbose output */
2963 if ( IfStringTrue(GetImageArtifact(image,"verbose")) ) {
2965 switch (sparse_method) {
2966 case BarycentricColorInterpolate:
2968 register ssize_t x=0;
2969 (void) FormatLocaleFile(stderr, "Barycentric Sparse Color:\n");
2970 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2971 (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf' \\\n",
2972 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2973 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2974 (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf' \\\n",
2975 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2976 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2977 (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf' \\\n",
2978 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2979 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2980 (image->colorspace == CMYKColorspace))
2981 (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf' \\\n",
2982 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2983 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2984 (image->alpha_trait == BlendPixelTrait))
2985 (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf' \\\n",
2986 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2989 case BilinearColorInterpolate:
2991 register ssize_t x=0;
2992 (void) FormatLocaleFile(stderr, "Bilinear Sparse Color\n");
2993 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2994 (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2995 coeff[ x ], coeff[x+1],
2996 coeff[x+2], coeff[x+3]),x+=4;
2997 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2998 (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2999 coeff[ x ], coeff[x+1],
3000 coeff[x+2], coeff[x+3]),x+=4;
3001 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3002 (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3003 coeff[ x ], coeff[x+1],
3004 coeff[x+2], coeff[x+3]),x+=4;
3005 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3006 (image->colorspace == CMYKColorspace))
3007 (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3008 coeff[ x ], coeff[x+1],
3009 coeff[x+2], coeff[x+3]),x+=4;
3010 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3011 (image->alpha_trait == BlendPixelTrait))
3012 (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3013 coeff[ x ], coeff[x+1],
3014 coeff[x+2], coeff[x+3]),x+=4;
3018 /* sparse color method is too complex for FX emulation */
3023 /* Generate new image for generated interpolated gradient.
3024 * ASIDE: Actually we could have just replaced the colors of the original
3025 * image, but IM Core policy, is if storage class could change then clone
3029 sparse_image=CloneImage(image,0,0,MagickTrue,exception);
3030 if (sparse_image == (Image *) NULL)
3031 return((Image *) NULL);
3032 if (SetImageStorageClass(sparse_image,DirectClass,exception) == MagickFalse)
3033 { /* if image is ColorMapped - change it to DirectClass */
3034 sparse_image=DestroyImage(sparse_image);
3035 return((Image *) NULL);
3037 { /* ----- MAIN CODE ----- */
3052 sparse_view=AcquireAuthenticCacheView(sparse_image,exception);
3053 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3054 #pragma omp parallel for schedule(static,4) shared(progress,status) \
3055 magick_threads(image,sparse_image,sparse_image->rows,1)
3057 for (j=0; j < (ssize_t) sparse_image->rows; j++)
3063 pixel; /* pixel to assign to distorted image */
3071 q=GetCacheViewAuthenticPixels(sparse_view,0,j,sparse_image->columns,
3073 if (q == (Quantum *) NULL)
3078 GetPixelInfo(sparse_image,&pixel);
3079 for (i=0; i < (ssize_t) image->columns; i++)
3081 GetPixelInfoPixel(image,q,&pixel);
3082 switch (sparse_method)
3084 case BarycentricColorInterpolate:
3086 register ssize_t x=0;
3087 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3088 pixel.red = coeff[x]*i +coeff[x+1]*j
3090 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3091 pixel.green = coeff[x]*i +coeff[x+1]*j
3093 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3094 pixel.blue = coeff[x]*i +coeff[x+1]*j
3096 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3097 (image->colorspace == CMYKColorspace))
3098 pixel.black = coeff[x]*i +coeff[x+1]*j
3100 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3101 (image->alpha_trait == BlendPixelTrait))
3102 pixel.alpha = coeff[x]*i +coeff[x+1]*j
3106 case BilinearColorInterpolate:
3108 register ssize_t x=0;
3109 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3110 pixel.red = coeff[x]*i + coeff[x+1]*j +
3111 coeff[x+2]*i*j + coeff[x+3], x+=4;
3112 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3113 pixel.green = coeff[x]*i + coeff[x+1]*j +
3114 coeff[x+2]*i*j + coeff[x+3], x+=4;
3115 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3116 pixel.blue = coeff[x]*i + coeff[x+1]*j +
3117 coeff[x+2]*i*j + coeff[x+3], x+=4;
3118 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3119 (image->colorspace == CMYKColorspace))
3120 pixel.black = coeff[x]*i + coeff[x+1]*j +
3121 coeff[x+2]*i*j + coeff[x+3], x+=4;
3122 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3123 (image->alpha_trait == BlendPixelTrait))
3124 pixel.alpha = coeff[x]*i + coeff[x+1]*j +
3125 coeff[x+2]*i*j + coeff[x+3], x+=4;
3128 case InverseColorInterpolate:
3129 case ShepardsColorInterpolate:
3130 { /* Inverse (Squared) Distance weights average (IDW) */
3136 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3138 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3140 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3142 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3143 (image->colorspace == CMYKColorspace))
3145 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3146 (image->alpha_trait == BlendPixelTrait))
3149 for(k=0; k<number_arguments; k+=2+number_colors) {
3150 register ssize_t x=(ssize_t) k+2;
3152 ((double)i-arguments[ k ])*((double)i-arguments[ k ])
3153 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
3154 weight = pow(weight,coeff[0]); /* inverse of power factor */
3155 weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
3156 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3157 pixel.red += arguments[x++]*weight;
3158 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3159 pixel.green += arguments[x++]*weight;
3160 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3161 pixel.blue += arguments[x++]*weight;
3162 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3163 (image->colorspace == CMYKColorspace))
3164 pixel.black += arguments[x++]*weight;
3165 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3166 (image->alpha_trait == BlendPixelTrait))
3167 pixel.alpha += arguments[x++]*weight;
3168 denominator += weight;
3170 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3171 pixel.red/=denominator;
3172 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3173 pixel.green/=denominator;
3174 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3175 pixel.blue/=denominator;
3176 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3177 (image->colorspace == CMYKColorspace))
3178 pixel.black/=denominator;
3179 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3180 (image->alpha_trait == BlendPixelTrait))
3181 pixel.alpha/=denominator;
3184 case VoronoiColorInterpolate:
3186 { /* Just use the closest control point you can find! */
3190 minimum = MagickHuge;
3192 for(k=0; k<number_arguments; k+=2+number_colors) {
3194 ((double)i-arguments[ k ])*((double)i-arguments[ k ])
3195 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
3196 if ( distance < minimum ) {
3197 register ssize_t x=(ssize_t) k+2;
3198 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3199 pixel.red=arguments[x++];
3200 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3201 pixel.green=arguments[x++];
3202 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3203 pixel.blue=arguments[x++];
3204 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3205 (image->colorspace == CMYKColorspace))
3206 pixel.black=arguments[x++];
3207 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3208 (image->alpha_trait == BlendPixelTrait))
3209 pixel.alpha=arguments[x++];
3216 /* set the color directly back into the source image */
3217 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3218 pixel.red*=QuantumRange;
3219 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3220 pixel.green*=QuantumRange;
3221 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3222 pixel.blue*=QuantumRange;
3223 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3224 (image->colorspace == CMYKColorspace))
3225 pixel.black*=QuantumRange;
3226 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3227 (image->alpha_trait == BlendPixelTrait))
3228 pixel.alpha*=QuantumRange;
3229 SetPixelInfoPixel(sparse_image,&pixel,q);
3230 q+=GetPixelChannels(sparse_image);
3232 sync=SyncCacheViewAuthenticPixels(sparse_view,exception);
3233 if (sync == MagickFalse)
3235 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3240 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3241 #pragma omp critical (MagickCore_SparseColorImage)
3243 proceed=SetImageProgress(image,SparseColorTag,progress++,image->rows);
3244 if (proceed == MagickFalse)
3248 sparse_view=DestroyCacheView(sparse_view);
3249 if (status == MagickFalse)
3250 sparse_image=DestroyImage(sparse_image);
3252 coeff = (double *) RelinquishMagickMemory(coeff);
3253 return(sparse_image);