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-2011 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/colorspace-private.h"
48 #include "MagickCore/composite-private.h"
49 #include "MagickCore/distort.h"
50 #include "MagickCore/exception.h"
51 #include "MagickCore/exception-private.h"
52 #include "MagickCore/gem.h"
53 #include "MagickCore/hashmap.h"
54 #include "MagickCore/image.h"
55 #include "MagickCore/list.h"
56 #include "MagickCore/matrix.h"
57 #include "MagickCore/matrix-private.h"
58 #include "MagickCore/memory_.h"
59 #include "MagickCore/monitor-private.h"
60 #include "MagickCore/option.h"
61 #include "MagickCore/pixel.h"
62 #include "MagickCore/pixel-accessor.h"
63 #include "MagickCore/resample.h"
64 #include "MagickCore/resample-private.h"
65 #include "MagickCore/registry.h"
66 #include "MagickCore/semaphore.h"
67 #include "MagickCore/string_.h"
68 #include "MagickCore/string-private.h"
69 #include "MagickCore/thread-private.h"
70 #include "MagickCore/token.h"
71 #include "MagickCore/transform.h"
74 Numerous internal routines for image distortions.
76 static inline double MagickMin(const double x,const double y)
78 return( x < y ? x : y);
80 static inline double MagickMax(const double x,const double y)
82 return( x > y ? x : y);
85 static inline void AffineArgsToCoefficients(double *affine)
87 /* map external sx,ry,rx,sy,tx,ty to internal c0,c2,c4,c1,c3,c5 */
88 double tmp[4]; /* note indexes 0 and 5 remain unchanged */
89 tmp[0]=affine[1]; tmp[1]=affine[2]; tmp[2]=affine[3]; tmp[3]=affine[4];
90 affine[3]=tmp[0]; affine[1]=tmp[1]; affine[4]=tmp[2]; affine[2]=tmp[3];
93 static inline void CoefficientsToAffineArgs(double *coeff)
95 /* map internal c0,c1,c2,c3,c4,c5 to external sx,ry,rx,sy,tx,ty */
96 double tmp[4]; /* note indexes 0 and 5 remain unchanged */
97 tmp[0]=coeff[3]; tmp[1]=coeff[1]; tmp[2]=coeff[4]; tmp[3]=coeff[2];
98 coeff[1]=tmp[0]; coeff[2]=tmp[1]; coeff[3]=tmp[2]; coeff[4]=tmp[3];
100 static void InvertAffineCoefficients(const double *coeff,double *inverse)
102 /* From "Digital Image Warping" by George Wolberg, page 50 */
105 determinant=1.0/(coeff[0]*coeff[4]-coeff[1]*coeff[3]);
106 inverse[0]=determinant*coeff[4];
107 inverse[1]=determinant*(-coeff[1]);
108 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[2]*coeff[4]);
109 inverse[3]=determinant*(-coeff[3]);
110 inverse[4]=determinant*coeff[0];
111 inverse[5]=determinant*(coeff[2]*coeff[3]-coeff[0]*coeff[5]);
114 static void InvertPerspectiveCoefficients(const double *coeff,
117 /* From "Digital Image Warping" by George Wolberg, page 53 */
120 determinant=1.0/(coeff[0]*coeff[4]-coeff[3]*coeff[1]);
121 inverse[0]=determinant*(coeff[4]-coeff[7]*coeff[5]);
122 inverse[1]=determinant*(coeff[7]*coeff[2]-coeff[1]);
123 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[4]*coeff[2]);
124 inverse[3]=determinant*(coeff[6]*coeff[5]-coeff[3]);
125 inverse[4]=determinant*(coeff[0]-coeff[6]*coeff[2]);
126 inverse[5]=determinant*(coeff[3]*coeff[2]-coeff[0]*coeff[5]);
127 inverse[6]=determinant*(coeff[3]*coeff[7]-coeff[6]*coeff[4]);
128 inverse[7]=determinant*(coeff[6]*coeff[1]-coeff[0]*coeff[7]);
131 static inline double MagickRound(double x)
134 Round the fraction to nearest integer.
137 return((double) ((ssize_t) (x+0.5)));
138 return((double) ((ssize_t) (x-0.5)));
142 * Polynomial Term Defining Functions
144 * Order must either be an integer, or 1.5 to produce
145 * the 2 number_valuesal polynomial function...
146 * affine 1 (3) u = c0 + c1*x + c2*y
147 * bilinear 1.5 (4) u = '' + c3*x*y
148 * quadratic 2 (6) u = '' + c4*x*x + c5*y*y
149 * cubic 3 (10) u = '' + c6*x^3 + c7*x*x*y + c8*x*y*y + c9*y^3
150 * quartic 4 (15) u = '' + c10*x^4 + ... + c14*y^4
151 * quintic 5 (21) u = '' + c15*x^5 + ... + c20*y^5
152 * number in parenthesis minimum number of points needed.
153 * Anything beyond quintic, has not been implemented until
154 * a more automated way of determining terms is found.
156 * Note the slight re-ordering of the terms for a quadratic polynomial
157 * which is to allow the use of a bi-linear (order=1.5) polynomial.
158 * All the later polynomials are ordered simply from x^N to y^N
160 static size_t poly_number_terms(double order)
162 /* Return the number of terms for a 2d polynomial */
163 if ( order < 1 || order > 5 ||
164 ( order != floor(order) && (order-1.5) > MagickEpsilon) )
165 return 0; /* invalid polynomial order */
166 return((size_t) floor((order+1)*(order+2)/2));
169 static double poly_basis_fn(ssize_t n, double x, double y)
171 /* Return the result for this polynomial term */
173 case 0: return( 1.0 ); /* constant */
175 case 2: return( y ); /* affine order = 1 terms = 3 */
176 case 3: return( x*y ); /* bilinear order = 1.5 terms = 4 */
177 case 4: return( x*x );
178 case 5: return( y*y ); /* quadratic order = 2 terms = 6 */
179 case 6: return( x*x*x );
180 case 7: return( x*x*y );
181 case 8: return( x*y*y );
182 case 9: return( y*y*y ); /* cubic order = 3 terms = 10 */
183 case 10: return( x*x*x*x );
184 case 11: return( x*x*x*y );
185 case 12: return( x*x*y*y );
186 case 13: return( x*y*y*y );
187 case 14: return( y*y*y*y ); /* quartic order = 4 terms = 15 */
188 case 15: return( x*x*x*x*x );
189 case 16: return( x*x*x*x*y );
190 case 17: return( x*x*x*y*y );
191 case 18: return( x*x*y*y*y );
192 case 19: return( x*y*y*y*y );
193 case 20: return( y*y*y*y*y ); /* quintic order = 5 terms = 21 */
195 return( 0 ); /* should never happen */
197 static const char *poly_basis_str(ssize_t n)
199 /* return the result for this polynomial term */
201 case 0: return(""); /* constant */
202 case 1: return("*ii");
203 case 2: return("*jj"); /* affine order = 1 terms = 3 */
204 case 3: return("*ii*jj"); /* bilinear order = 1.5 terms = 4 */
205 case 4: return("*ii*ii");
206 case 5: return("*jj*jj"); /* quadratic order = 2 terms = 6 */
207 case 6: return("*ii*ii*ii");
208 case 7: return("*ii*ii*jj");
209 case 8: return("*ii*jj*jj");
210 case 9: return("*jj*jj*jj"); /* cubic order = 3 terms = 10 */
211 case 10: return("*ii*ii*ii*ii");
212 case 11: return("*ii*ii*ii*jj");
213 case 12: return("*ii*ii*jj*jj");
214 case 13: return("*ii*jj*jj*jj");
215 case 14: return("*jj*jj*jj*jj"); /* quartic order = 4 terms = 15 */
216 case 15: return("*ii*ii*ii*ii*ii");
217 case 16: return("*ii*ii*ii*ii*jj");
218 case 17: return("*ii*ii*ii*jj*jj");
219 case 18: return("*ii*ii*jj*jj*jj");
220 case 19: return("*ii*jj*jj*jj*jj");
221 case 20: return("*jj*jj*jj*jj*jj"); /* quintic order = 5 terms = 21 */
223 return( "UNKNOWN" ); /* should never happen */
225 static double poly_basis_dx(ssize_t n, double x, double y)
227 /* polynomial term for x derivative */
229 case 0: return( 0.0 ); /* constant */
230 case 1: return( 1.0 );
231 case 2: return( 0.0 ); /* affine order = 1 terms = 3 */
232 case 3: return( y ); /* bilinear order = 1.5 terms = 4 */
234 case 5: return( 0.0 ); /* quadratic order = 2 terms = 6 */
235 case 6: return( x*x );
236 case 7: return( x*y );
237 case 8: return( y*y );
238 case 9: return( 0.0 ); /* cubic order = 3 terms = 10 */
239 case 10: return( x*x*x );
240 case 11: return( x*x*y );
241 case 12: return( x*y*y );
242 case 13: return( y*y*y );
243 case 14: return( 0.0 ); /* quartic order = 4 terms = 15 */
244 case 15: return( x*x*x*x );
245 case 16: return( x*x*x*y );
246 case 17: return( x*x*y*y );
247 case 18: return( x*y*y*y );
248 case 19: return( y*y*y*y );
249 case 20: return( 0.0 ); /* quintic order = 5 terms = 21 */
251 return( 0.0 ); /* should never happen */
253 static double poly_basis_dy(ssize_t n, double x, double y)
255 /* polynomial term for y derivative */
257 case 0: return( 0.0 ); /* constant */
258 case 1: return( 0.0 );
259 case 2: return( 1.0 ); /* affine order = 1 terms = 3 */
260 case 3: return( x ); /* bilinear order = 1.5 terms = 4 */
261 case 4: return( 0.0 );
262 case 5: return( y ); /* quadratic order = 2 terms = 6 */
263 default: return( poly_basis_dx(n-1,x,y) ); /* weird but true */
265 /* NOTE: the only reason that last is not true for 'quadratic'
266 is due to the re-arrangement of terms to allow for 'bilinear'
271 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
275 + G e n e r a t e C o e f f i c i e n t s %
279 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
281 % GenerateCoefficients() takes user provided input arguments and generates
282 % the coefficients, needed to apply the specific distortion for either
283 % distorting images (generally using control points) or generating a color
284 % gradient from sparsely separated color points.
286 % The format of the GenerateCoefficients() method is:
288 % Image *GenerateCoefficients(const Image *image,DistortImageMethod method,
289 % const size_t number_arguments,const double *arguments,
290 % size_t number_values, ExceptionInfo *exception)
292 % A description of each parameter follows:
294 % o image: the image to be distorted.
296 % o method: the method of image distortion/ sparse gradient
298 % o number_arguments: the number of arguments given.
300 % o arguments: the arguments for this distortion method.
302 % o number_values: the style and format of given control points, (caller type)
303 % 0: 2 dimensional mapping of control points (Distort)
304 % Format: u,v,x,y where u,v is the 'source' of the
305 % the color to be plotted, for DistortImage()
306 % N: Interpolation of control points with N values (usally r,g,b)
307 % Format: x,y,r,g,b mapping x,y to color values r,g,b
308 % IN future, variable number of values may be given (1 to N)
310 % o exception: return any errors or warnings in this structure
312 % Note that the returned array of double values must be freed by the
313 % calling method using RelinquishMagickMemory(). This however may change in
314 % the future to require a more 'method' specific method.
316 % Because of this this method should not be classed as stable or used
317 % outside other MagickCore library methods.
320 static double *GenerateCoefficients(const Image *image,
321 DistortImageMethod *method,const size_t number_arguments,
322 const double *arguments,size_t number_values,ExceptionInfo *exception)
331 number_coeff, /* number of coefficients to return (array size) */
332 cp_size, /* number floating point numbers per control point */
333 cp_x,cp_y, /* the x,y indexes for control point */
334 cp_values; /* index of values for this control point */
335 /* number_values Number of values given per control point */
337 if ( number_values == 0 ) {
338 /* Image distortion using control points (or other distortion)
339 That is generate a mapping so that x,y->u,v given u,v,x,y
341 number_values = 2; /* special case: two values of u,v */
342 cp_values = 0; /* the values i,j are BEFORE the destination CP x,y */
343 cp_x = 2; /* location of x,y in input control values */
345 /* NOTE: cp_values, also used for later 'reverse map distort' tests */
348 cp_x = 0; /* location of x,y in input control values */
350 cp_values = 2; /* and the other values are after x,y */
351 /* Typically in this case the values are R,G,B color values */
353 cp_size = number_values+2; /* each CP defintion involves this many numbers */
355 /* If not enough control point pairs are found for specific distortions
356 fall back to Affine distortion (allowing 0 to 3 point pairs)
358 if ( number_arguments < 4*cp_size &&
359 ( *method == BilinearForwardDistortion
360 || *method == BilinearReverseDistortion
361 || *method == PerspectiveDistortion
363 *method = AffineDistortion;
367 case AffineDistortion:
368 /* also BarycentricColorInterpolate: */
369 number_coeff=3*number_values;
371 case PolynomialDistortion:
372 /* number of coefficents depend on the given polynomal 'order' */
373 if ( number_arguments <= 1 && (number_arguments-1)%cp_size != 0)
375 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
376 "InvalidArgument","%s : '%s'","Polynomial",
377 "Invalid number of args: order [CPs]...");
378 return((double *) NULL);
380 i = poly_number_terms(arguments[0]);
381 number_coeff = 2 + i*number_values;
383 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
384 "InvalidArgument","%s : '%s'","Polynomial",
385 "Invalid order, should be interger 1 to 5, or 1.5");
386 return((double *) NULL);
388 if ( number_arguments < 1+i*cp_size ) {
389 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
390 "InvalidArgument", "%s : 'require at least %.20g CPs'",
391 "Polynomial", (double) i);
392 return((double *) NULL);
395 case BilinearReverseDistortion:
396 number_coeff=4*number_values;
399 The rest are constants as they are only used for image distorts
401 case BilinearForwardDistortion:
402 number_coeff=10; /* 2*4 coeff plus 2 constants */
403 cp_x = 0; /* Reverse src/dest coords for forward mapping */
408 case QuadraterialDistortion:
409 number_coeff=19; /* BilinearForward + BilinearReverse */
412 case ShepardsDistortion:
413 number_coeff=1; /* not used, but provide some type of return */
418 case ScaleRotateTranslateDistortion:
419 case AffineProjectionDistortion:
420 case Plane2CylinderDistortion:
421 case Cylinder2PlaneDistortion:
424 case PolarDistortion:
425 case DePolarDistortion:
428 case PerspectiveDistortion:
429 case PerspectiveProjectionDistortion:
432 case BarrelDistortion:
433 case BarrelInverseDistortion:
437 assert(! "Unknown Method Given"); /* just fail assertion */
440 /* allocate the array of coefficients needed */
441 coeff = (double *) AcquireQuantumMemory(number_coeff,sizeof(*coeff));
442 if (coeff == (double *) NULL) {
443 (void) ThrowMagickException(exception,GetMagickModule(),
444 ResourceLimitError,"MemoryAllocationFailed",
445 "%s", "GenerateCoefficients");
446 return((double *) NULL);
449 /* zero out coefficients array */
450 for (i=0; i < number_coeff; i++)
455 case AffineDistortion:
459 for each 'value' given
461 Input Arguments are sets of control points...
462 For Distort Images u,v, x,y ...
463 For Sparse Gradients x,y, r,g,b ...
465 if ( number_arguments%cp_size != 0 ||
466 number_arguments < cp_size ) {
467 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
468 "InvalidArgument", "%s : 'require at least %.20g CPs'",
470 coeff=(double *) RelinquishMagickMemory(coeff);
471 return((double *) NULL);
473 /* handle special cases of not enough arguments */
474 if ( number_arguments == cp_size ) {
475 /* Only 1 CP Set Given */
476 if ( cp_values == 0 ) {
477 /* image distortion - translate the image */
479 coeff[2] = arguments[0] - arguments[2];
481 coeff[5] = arguments[1] - arguments[3];
484 /* sparse gradient - use the values directly */
485 for (i=0; i<number_values; i++)
486 coeff[i*3+2] = arguments[cp_values+i];
490 /* 2 or more points (usally 3) given.
491 Solve a least squares simultaneous equation for coefficients.
501 /* create matrix, and a fake vectors matrix */
502 matrix = AcquireMagickMatrix(3UL,3UL);
503 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
504 if (matrix == (double **) NULL || vectors == (double **) NULL)
506 matrix = RelinquishMagickMatrix(matrix, 3UL);
507 vectors = (double **) RelinquishMagickMemory(vectors);
508 coeff = (double *) RelinquishMagickMemory(coeff);
509 (void) ThrowMagickException(exception,GetMagickModule(),
510 ResourceLimitError,"MemoryAllocationFailed",
511 "%s", "DistortCoefficients");
512 return((double *) NULL);
514 /* fake a number_values x3 vectors matrix from coefficients array */
515 for (i=0; i < number_values; i++)
516 vectors[i] = &(coeff[i*3]);
517 /* Add given control point pairs for least squares solving */
518 for (i=0; i < number_arguments; i+=cp_size) {
519 terms[0] = arguments[i+cp_x]; /* x */
520 terms[1] = arguments[i+cp_y]; /* y */
521 terms[2] = 1; /* 1 */
522 LeastSquaresAddTerms(matrix,vectors,terms,
523 &(arguments[i+cp_values]),3UL,number_values);
525 if ( number_arguments == 2*cp_size ) {
526 /* Only two pairs were given, but we need 3 to solve the affine.
527 Fake extra coordinates by rotating p1 around p0 by 90 degrees.
528 x2 = x0 - (y1-y0) y2 = y0 + (x1-x0)
530 terms[0] = arguments[cp_x]
531 - ( arguments[cp_size+cp_y] - arguments[cp_y] ); /* x2 */
532 terms[1] = arguments[cp_y] +
533 + ( arguments[cp_size+cp_x] - arguments[cp_x] ); /* y2 */
534 terms[2] = 1; /* 1 */
535 if ( cp_values == 0 ) {
536 /* Image Distortion - rotate the u,v coordients too */
539 uv2[0] = arguments[0] - arguments[5] + arguments[1]; /* u2 */
540 uv2[1] = arguments[1] + arguments[4] - arguments[0]; /* v2 */
541 LeastSquaresAddTerms(matrix,vectors,terms,uv2,3UL,2UL);
544 /* Sparse Gradient - use values of p0 for linear gradient */
545 LeastSquaresAddTerms(matrix,vectors,terms,
546 &(arguments[cp_values]),3UL,number_values);
549 /* Solve for LeastSquares Coefficients */
550 status=GaussJordanElimination(matrix,vectors,3UL,number_values);
551 matrix = RelinquishMagickMatrix(matrix, 3UL);
552 vectors = (double **) RelinquishMagickMemory(vectors);
553 if ( status == MagickFalse ) {
554 coeff = (double *) RelinquishMagickMemory(coeff);
555 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
556 "InvalidArgument","%s : 'Unsolvable Matrix'",
557 CommandOptionToMnemonic(MagickDistortOptions, *method) );
558 return((double *) NULL);
563 case AffineProjectionDistortion:
566 Arguments: Affine Matrix (forward mapping)
567 Arguments sx, rx, ry, sy, tx, ty
568 Where u = sx*x + ry*y + tx
571 Returns coefficients (in there inverse form) ordered as...
574 AffineProjection Distortion Notes...
575 + Will only work with a 2 number_values for Image Distortion
576 + Can not be used for generating a sparse gradient (interpolation)
579 if (number_arguments != 6) {
580 coeff = (double *) RelinquishMagickMemory(coeff);
581 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
582 "InvalidArgument","%s : 'Needs 6 coeff values'",
583 CommandOptionToMnemonic(MagickDistortOptions, *method) );
584 return((double *) NULL);
586 /* FUTURE: trap test for sx*sy-rx*ry == 0 (determinant = 0, no inverse) */
587 for(i=0; i<6UL; i++ )
588 inverse[i] = arguments[i];
589 AffineArgsToCoefficients(inverse); /* map into coefficents */
590 InvertAffineCoefficients(inverse, coeff); /* invert */
591 *method = AffineDistortion;
595 case ScaleRotateTranslateDistortion:
597 /* Scale, Rotate and Translate Distortion
598 An alternative Affine Distortion
599 Argument options, by number of arguments given:
600 7: x,y, sx,sy, a, nx,ny
607 Where actions are (in order of application)
608 x,y 'center' of transforms (default = image center)
609 sx,sy scale image by this amount (default = 1)
610 a angle of rotation (argument required)
611 nx,ny move 'center' here (default = x,y or no movement)
612 And convert to affine mapping coefficients
614 ScaleRotateTranslate Distortion Notes...
615 + Does not use a set of CPs in any normal way
616 + Will only work with a 2 number_valuesal Image Distortion
617 + Cannot be used for generating a sparse gradient (interpolation)
623 /* set default center, and default scale */
624 x = nx = (double)(image->columns)/2.0 + (double)image->page.x;
625 y = ny = (double)(image->rows)/2.0 + (double)image->page.y;
627 switch ( number_arguments ) {
629 coeff = (double *) RelinquishMagickMemory(coeff);
630 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
631 "InvalidArgument","%s : 'Needs at least 1 argument'",
632 CommandOptionToMnemonic(MagickDistortOptions, *method) );
633 return((double *) NULL);
638 sx = sy = arguments[0];
642 x = nx = arguments[0];
643 y = ny = arguments[1];
644 switch ( number_arguments ) {
649 sx = sy = arguments[2];
658 sx = sy = arguments[2];
671 coeff = (double *) RelinquishMagickMemory(coeff);
672 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
673 "InvalidArgument","%s : 'Too Many Arguments (7 or less)'",
674 CommandOptionToMnemonic(MagickDistortOptions, *method) );
675 return((double *) NULL);
679 /* Trap if sx or sy == 0 -- image is scaled out of existance! */
680 if ( fabs(sx) < MagickEpsilon || fabs(sy) < MagickEpsilon ) {
681 coeff = (double *) RelinquishMagickMemory(coeff);
682 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
683 "InvalidArgument","%s : 'Zero Scale Given'",
684 CommandOptionToMnemonic(MagickDistortOptions, *method) );
685 return((double *) NULL);
687 /* Save the given arguments as an affine distortion */
688 a=DegreesToRadians(a); cosine=cos(a); sine=sin(a);
690 *method = AffineDistortion;
693 coeff[2]=x-nx*coeff[0]-ny*coeff[1];
696 coeff[5]=y-nx*coeff[3]-ny*coeff[4];
699 case PerspectiveDistortion:
701 Perspective Distortion (a ratio of affine distortions)
703 p(x,y) c0*x + c1*y + c2
704 u = ------ = ------------------
705 r(x,y) c6*x + c7*y + 1
707 q(x,y) c3*x + c4*y + c5
708 v = ------ = ------------------
709 r(x,y) c6*x + c7*y + 1
711 c8 = Sign of 'r', or the denominator affine, for the actual image.
712 This determines what part of the distorted image is 'ground'
713 side of the horizon, the other part is 'sky' or invalid.
714 Valid values are +1.0 or -1.0 only.
716 Input Arguments are sets of control points...
717 For Distort Images u,v, x,y ...
718 For Sparse Gradients x,y, r,g,b ...
720 Perspective Distortion Notes...
721 + Can be thought of as ratio of 3 affine transformations
722 + Not separatable: r() or c6 and c7 are used by both equations
723 + All 8 coefficients must be determined simultaniously
724 + Will only work with a 2 number_valuesal Image Distortion
725 + Can not be used for generating a sparse gradient (interpolation)
726 + It is not linear, but is simple to generate an inverse
727 + All lines within an image remain lines.
728 + but distances between points may vary.
742 if ( number_arguments%cp_size != 0 ||
743 number_arguments < cp_size*4 ) {
744 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
745 "InvalidArgument", "%s : 'require at least %.20g CPs'",
746 CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0);
747 coeff=(double *) RelinquishMagickMemory(coeff);
748 return((double *) NULL);
750 /* fake 1x8 vectors matrix directly using the coefficients array */
751 vectors[0] = &(coeff[0]);
752 /* 8x8 least-squares matrix (zeroed) */
753 matrix = AcquireMagickMatrix(8UL,8UL);
754 if (matrix == (double **) NULL) {
755 (void) ThrowMagickException(exception,GetMagickModule(),
756 ResourceLimitError,"MemoryAllocationFailed",
757 "%s", "DistortCoefficients");
758 return((double *) NULL);
760 /* Add control points for least squares solving */
761 for (i=0; i < number_arguments; i+=4) {
762 terms[0]=arguments[i+cp_x]; /* c0*x */
763 terms[1]=arguments[i+cp_y]; /* c1*y */
764 terms[2]=1.0; /* c2*1 */
768 terms[6]=-terms[0]*arguments[i+cp_u]; /* 1/(c6*x) */
769 terms[7]=-terms[1]*arguments[i+cp_u]; /* 1/(c7*y) */
770 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_u]),
776 terms[3]=arguments[i+cp_x]; /* c3*x */
777 terms[4]=arguments[i+cp_y]; /* c4*y */
778 terms[5]=1.0; /* c5*1 */
779 terms[6]=-terms[3]*arguments[i+cp_v]; /* 1/(c6*x) */
780 terms[7]=-terms[4]*arguments[i+cp_v]; /* 1/(c7*y) */
781 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_v]),
784 /* Solve for LeastSquares Coefficients */
785 status=GaussJordanElimination(matrix,vectors,8UL,1UL);
786 matrix = RelinquishMagickMatrix(matrix, 8UL);
787 if ( status == MagickFalse ) {
788 coeff = (double *) RelinquishMagickMemory(coeff);
789 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
790 "InvalidArgument","%s : 'Unsolvable Matrix'",
791 CommandOptionToMnemonic(MagickDistortOptions, *method) );
792 return((double *) NULL);
795 Calculate 9'th coefficient! The ground-sky determination.
796 What is sign of the 'ground' in r() denominator affine function?
797 Just use any valid image coordinate (first control point) in
798 destination for determination of what part of view is 'ground'.
800 coeff[8] = coeff[6]*arguments[cp_x]
801 + coeff[7]*arguments[cp_y] + 1.0;
802 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
806 case PerspectiveProjectionDistortion:
809 Arguments: Perspective Coefficents (forward mapping)
811 if (number_arguments != 8) {
812 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
813 "InvalidArgument", "%s : 'Needs 8 coefficient values'",
814 CommandOptionToMnemonic(MagickDistortOptions, *method));
815 return((double *) NULL);
817 /* FUTURE: trap test c0*c4-c3*c1 == 0 (determinate = 0, no inverse) */
818 InvertPerspectiveCoefficients(arguments, coeff);
820 Calculate 9'th coefficient! The ground-sky determination.
821 What is sign of the 'ground' in r() denominator affine function?
822 Just use any valid image cocodinate in destination for determination.
823 For a forward mapped perspective the images 0,0 coord will map to
824 c2,c5 in the distorted image, so set the sign of denominator of that.
826 coeff[8] = coeff[6]*arguments[2]
827 + coeff[7]*arguments[5] + 1.0;
828 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
829 *method = PerspectiveDistortion;
833 case BilinearForwardDistortion:
834 case BilinearReverseDistortion:
836 /* Bilinear Distortion (Forward mapping)
837 v = c0*x + c1*y + c2*x*y + c3;
838 for each 'value' given
840 This is actually a simple polynomial Distortion! The difference
841 however is when we need to reverse the above equation to generate a
842 BilinearForwardDistortion (see below).
844 Input Arguments are sets of control points...
845 For Distort Images u,v, x,y ...
846 For Sparse Gradients x,y, r,g,b ...
857 /* check the number of arguments */
858 if ( number_arguments%cp_size != 0 ||
859 number_arguments < cp_size*4 ) {
860 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
861 "InvalidArgument", "%s : 'require at least %.20g CPs'",
862 CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0);
863 coeff=(double *) RelinquishMagickMemory(coeff);
864 return((double *) NULL);
866 /* create matrix, and a fake vectors matrix */
867 matrix = AcquireMagickMatrix(4UL,4UL);
868 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
869 if (matrix == (double **) NULL || vectors == (double **) NULL)
871 matrix = RelinquishMagickMatrix(matrix, 4UL);
872 vectors = (double **) RelinquishMagickMemory(vectors);
873 coeff = (double *) RelinquishMagickMemory(coeff);
874 (void) ThrowMagickException(exception,GetMagickModule(),
875 ResourceLimitError,"MemoryAllocationFailed",
876 "%s", "DistortCoefficients");
877 return((double *) NULL);
879 /* fake a number_values x4 vectors matrix from coefficients array */
880 for (i=0; i < number_values; i++)
881 vectors[i] = &(coeff[i*4]);
882 /* Add given control point pairs for least squares solving */
883 for (i=0; i < number_arguments; i+=cp_size) {
884 terms[0] = arguments[i+cp_x]; /* x */
885 terms[1] = arguments[i+cp_y]; /* y */
886 terms[2] = terms[0]*terms[1]; /* x*y */
887 terms[3] = 1; /* 1 */
888 LeastSquaresAddTerms(matrix,vectors,terms,
889 &(arguments[i+cp_values]),4UL,number_values);
891 /* Solve for LeastSquares Coefficients */
892 status=GaussJordanElimination(matrix,vectors,4UL,number_values);
893 matrix = RelinquishMagickMatrix(matrix, 4UL);
894 vectors = (double **) RelinquishMagickMemory(vectors);
895 if ( status == MagickFalse ) {
896 coeff = (double *) RelinquishMagickMemory(coeff);
897 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
898 "InvalidArgument","%s : 'Unsolvable Matrix'",
899 CommandOptionToMnemonic(MagickDistortOptions, *method) );
900 return((double *) NULL);
902 if ( *method == BilinearForwardDistortion ) {
903 /* Bilinear Forward Mapped Distortion
905 The above least-squares solved for coefficents but in the forward
906 direction, due to changes to indexing constants.
908 i = c0*x + c1*y + c2*x*y + c3;
909 j = c4*x + c5*y + c6*x*y + c7;
911 where i,j are in the destination image, NOT the source.
913 Reverse Pixel mapping however needs to use reverse of these
914 functions. It required a full page of algbra to work out the
915 reversed mapping formula, but resolves down to the following...
918 c9 = 2*(c2*c5-c1*c6); // '2*a' in the quadratic formula
920 i = i - c3; j = j - c7;
921 b = c6*i - c2*j + c8; // So that a*y^2 + b*y + c == 0
922 c = c4*i - c0*j; // y = ( -b +- sqrt(bb - 4ac) ) / (2*a)
926 y = ( -b + sqrt(r) ) / c9;
930 x = ( i - c1*y) / ( c1 - c2*y );
932 NB: if 'r' is negative there is no solution!
933 NB: the sign of the sqrt() should be negative if image becomes
934 flipped or flopped, or crosses over itself.
935 NB: techniqually coefficient c5 is not needed, anymore,
936 but kept for completness.
938 See Anthony Thyssen <A.Thyssen@griffith.edu.au>
939 or Fred Weinhaus <fmw@alink.net> for more details.
942 coeff[8] = coeff[0]*coeff[5] - coeff[1]*coeff[4];
943 coeff[9] = 2*(coeff[2]*coeff[5] - coeff[1]*coeff[6]);
948 case QuadrilateralDistortion:
950 /* Map a Quadrilateral to a unit square using BilinearReverse
951 Then map that unit square back to the final Quadrilateral
952 using BilinearForward.
954 Input Arguments are sets of control points...
955 For Distort Images u,v, x,y ...
956 For Sparse Gradients x,y, r,g,b ...
959 /* UNDER CONSTRUCTION */
964 case PolynomialDistortion:
966 /* Polynomial Distortion
968 First two coefficents are used to hole global polynomal information
969 c0 = Order of the polynimial being created
970 c1 = number_of_terms in one polynomial equation
972 Rest of the coefficients map to the equations....
973 v = c0 + c1*x + c2*y + c3*x*y + c4*x^2 + c5*y^2 + c6*x^3 + ...
974 for each 'value' (number_values of them) given.
975 As such total coefficients = 2 + number_terms * number_values
977 Input Arguments are sets of control points...
978 For Distort Images order [u,v, x,y] ...
979 For Sparse Gradients order [x,y, r,g,b] ...
981 Polynomial Distortion Notes...
982 + UNDER DEVELOPMENT -- Do not expect this to remain as is.
983 + Currently polynomial is a reversed mapped distortion.
984 + Order 1.5 is fudged to map into a bilinear distortion.
985 though it is not the same order as that distortion.
993 nterms; /* number of polynomial terms per number_values */
1001 /* first two coefficients hold polynomial order information */
1002 coeff[0] = arguments[0];
1003 coeff[1] = (double) poly_number_terms(arguments[0]);
1004 nterms = (size_t) coeff[1];
1006 /* create matrix, a fake vectors matrix, and least sqs terms */
1007 matrix = AcquireMagickMatrix(nterms,nterms);
1008 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
1009 terms = (double *) AcquireQuantumMemory(nterms, sizeof(*terms));
1010 if (matrix == (double **) NULL ||
1011 vectors == (double **) NULL ||
1012 terms == (double *) NULL )
1014 matrix = RelinquishMagickMatrix(matrix, nterms);
1015 vectors = (double **) RelinquishMagickMemory(vectors);
1016 terms = (double *) RelinquishMagickMemory(terms);
1017 coeff = (double *) RelinquishMagickMemory(coeff);
1018 (void) ThrowMagickException(exception,GetMagickModule(),
1019 ResourceLimitError,"MemoryAllocationFailed",
1020 "%s", "DistortCoefficients");
1021 return((double *) NULL);
1023 /* fake a number_values x3 vectors matrix from coefficients array */
1024 for (i=0; i < number_values; i++)
1025 vectors[i] = &(coeff[2+i*nterms]);
1026 /* Add given control point pairs for least squares solving */
1027 for (i=1; i < number_arguments; i+=cp_size) { /* NB: start = 1 not 0 */
1028 for (j=0; j < (ssize_t) nterms; j++)
1029 terms[j] = poly_basis_fn(j,arguments[i+cp_x],arguments[i+cp_y]);
1030 LeastSquaresAddTerms(matrix,vectors,terms,
1031 &(arguments[i+cp_values]),nterms,number_values);
1033 terms = (double *) RelinquishMagickMemory(terms);
1034 /* Solve for LeastSquares Coefficients */
1035 status=GaussJordanElimination(matrix,vectors,nterms,number_values);
1036 matrix = RelinquishMagickMatrix(matrix, nterms);
1037 vectors = (double **) RelinquishMagickMemory(vectors);
1038 if ( status == MagickFalse ) {
1039 coeff = (double *) RelinquishMagickMemory(coeff);
1040 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1041 "InvalidArgument","%s : 'Unsolvable Matrix'",
1042 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1043 return((double *) NULL);
1050 Args: arc_width rotate top_edge_radius bottom_edge_radius
1051 All but first argument are optional
1052 arc_width The angle over which to arc the image side-to-side
1053 rotate Angle to rotate image from vertical center
1054 top_radius Set top edge of source image at this radius
1055 bottom_radius Set bootom edge to this radius (radial scaling)
1057 By default, if the radii arguments are nor provided the image radius
1058 is calculated so the horizontal center-line is fits the given arc
1061 The output image size is ALWAYS adjusted to contain the whole image,
1062 and an offset is given to position image relative to the 0,0 point of
1063 the origin, allowing users to use relative positioning onto larger
1064 background (via -flatten).
1066 The arguments are converted to these coefficients
1067 c0: angle for center of source image
1068 c1: angle scale for mapping to source image
1069 c2: radius for top of source image
1070 c3: radius scale for mapping source image
1071 c4: centerline of arc within source image
1073 Note the coefficients use a center angle, so asymptotic join is
1074 furthest from both sides of the source image. This also means that
1075 for arc angles greater than 360 the sides of the image will be
1078 Arc Distortion Notes...
1079 + Does not use a set of CPs
1080 + Will only work with Image Distortion
1081 + Can not be used for generating a sparse gradient (interpolation)
1083 if ( number_arguments >= 1 && arguments[0] < MagickEpsilon ) {
1084 coeff = (double *) RelinquishMagickMemory(coeff);
1085 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1086 "InvalidArgument","%s : 'Arc Angle Too Small'",
1087 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1088 return((double *) NULL);
1090 if ( number_arguments >= 3 && arguments[2] < MagickEpsilon ) {
1091 coeff = (double *) RelinquishMagickMemory(coeff);
1092 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1093 "InvalidArgument","%s : 'Outer Radius Too Small'",
1094 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1095 return((double *) NULL);
1097 coeff[0] = -MagickPI2; /* -90, place at top! */
1098 if ( number_arguments >= 1 )
1099 coeff[1] = DegreesToRadians(arguments[0]);
1101 coeff[1] = MagickPI2; /* zero arguments - center is at top */
1102 if ( number_arguments >= 2 )
1103 coeff[0] += DegreesToRadians(arguments[1]);
1104 coeff[0] /= Magick2PI; /* normalize radians */
1105 coeff[0] -= MagickRound(coeff[0]);
1106 coeff[0] *= Magick2PI; /* de-normalize back to radians */
1107 coeff[3] = (double)image->rows-1;
1108 coeff[2] = (double)image->columns/coeff[1] + coeff[3]/2.0;
1109 if ( number_arguments >= 3 ) {
1110 if ( number_arguments >= 4 )
1111 coeff[3] = arguments[2] - arguments[3];
1113 coeff[3] *= arguments[2]/coeff[2];
1114 coeff[2] = arguments[2];
1116 coeff[4] = ((double)image->columns-1.0)/2.0;
1120 case PolarDistortion:
1121 case DePolarDistortion:
1123 /* (De)Polar Distortion (same set of arguments)
1124 Args: Rmax, Rmin, Xcenter,Ycenter, Afrom,Ato
1125 DePolar can also have the extra arguments of Width, Height
1127 Coefficients 0 to 5 is the sanatized version first 6 input args
1128 Coefficient 6 is the angle to coord ratio and visa-versa
1129 Coefficient 7 is the radius to coord ratio and visa-versa
1131 WARNING: It is possible for Radius max<min and/or Angle from>to
1133 if ( number_arguments == 3
1134 || ( number_arguments > 6 && *method == PolarDistortion )
1135 || number_arguments > 8 ) {
1136 (void) ThrowMagickException(exception,GetMagickModule(),
1137 OptionError,"InvalidArgument", "%s : number of arguments",
1138 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1139 coeff=(double *) RelinquishMagickMemory(coeff);
1140 return((double *) NULL);
1142 /* Rmax - if 0 calculate appropriate value */
1143 if ( number_arguments >= 1 )
1144 coeff[0] = arguments[0];
1147 /* Rmin - usally 0 */
1148 coeff[1] = number_arguments >= 2 ? arguments[1] : 0.0;
1150 if ( number_arguments >= 4 ) {
1151 coeff[2] = arguments[2];
1152 coeff[3] = arguments[3];
1154 else { /* center of actual image */
1155 coeff[2] = (double)(image->columns)/2.0+image->page.x;
1156 coeff[3] = (double)(image->rows)/2.0+image->page.y;
1158 /* Angle from,to - about polar center 0 is downward */
1159 coeff[4] = -MagickPI;
1160 if ( number_arguments >= 5 )
1161 coeff[4] = DegreesToRadians(arguments[4]);
1162 coeff[5] = coeff[4];
1163 if ( number_arguments >= 6 )
1164 coeff[5] = DegreesToRadians(arguments[5]);
1165 if ( fabs(coeff[4]-coeff[5]) < MagickEpsilon )
1166 coeff[5] += Magick2PI; /* same angle is a full circle */
1167 /* if radius 0 or negative, its a special value... */
1168 if ( coeff[0] < MagickEpsilon ) {
1169 /* Use closest edge if radius == 0 */
1170 if ( fabs(coeff[0]) < MagickEpsilon ) {
1171 coeff[0]=MagickMin(fabs(coeff[2]-image->page.x),
1172 fabs(coeff[3]-image->page.y));
1173 coeff[0]=MagickMin(coeff[0],
1174 fabs(coeff[2]-image->page.x-image->columns));
1175 coeff[0]=MagickMin(coeff[0],
1176 fabs(coeff[3]-image->page.y-image->rows));
1178 /* furthest diagonal if radius == -1 */
1179 if ( fabs(-1.0-coeff[0]) < MagickEpsilon ) {
1181 rx = coeff[2]-image->page.x;
1182 ry = coeff[3]-image->page.y;
1183 coeff[0] = rx*rx+ry*ry;
1184 ry = coeff[3]-image->page.y-image->rows;
1185 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1186 rx = coeff[2]-image->page.x-image->columns;
1187 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1188 ry = coeff[3]-image->page.y;
1189 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1190 coeff[0] = sqrt(coeff[0]);
1193 /* IF Rmax <= 0 or Rmin < 0 OR Rmax < Rmin, THEN error */
1194 if ( coeff[0] < MagickEpsilon || coeff[1] < -MagickEpsilon
1195 || (coeff[0]-coeff[1]) < MagickEpsilon ) {
1196 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1197 "InvalidArgument", "%s : Invalid Radius",
1198 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1199 coeff=(double *) RelinquishMagickMemory(coeff);
1200 return((double *) NULL);
1202 /* converstion ratios */
1203 if ( *method == PolarDistortion ) {
1204 coeff[6]=(double) image->columns/(coeff[5]-coeff[4]);
1205 coeff[7]=(double) image->rows/(coeff[0]-coeff[1]);
1207 else { /* *method == DePolarDistortion */
1208 coeff[6]=(coeff[5]-coeff[4])/image->columns;
1209 coeff[7]=(coeff[0]-coeff[1])/image->rows;
1213 case Cylinder2PlaneDistortion:
1214 case Plane2CylinderDistortion:
1216 /* 3D Cylinder to/from a Tangential Plane
1218 Projection between a clinder and flat plain from a point on the
1219 center line of the cylinder.
1221 The two surfaces coincide in 3D space at the given centers of
1222 distortion (perpendicular to projection point) on both images.
1225 Coefficents: FOV(radians), Radius, center_x,y, dest_center_x,y
1227 FOV (Field Of View) the angular field of view of the distortion,
1228 across the width of the image, in degrees. The centers are the
1229 points of least distortion in the input and resulting images.
1231 These centers are however determined later.
1233 Coeff 0 is the FOV angle of view of image width in radians
1234 Coeff 1 is calculated radius of cylinder.
1235 Coeff 2,3 center of distortion of input image
1236 Coefficents 4,5 Center of Distortion of dest (determined later)
1238 if ( arguments[0] < MagickEpsilon || arguments[0] > 160.0 ) {
1239 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1240 "InvalidArgument", "%s : Invalid FOV Angle",
1241 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1242 coeff=(double *) RelinquishMagickMemory(coeff);
1243 return((double *) NULL);
1245 coeff[0] = DegreesToRadians(arguments[0]);
1246 if ( *method == Cylinder2PlaneDistortion )
1247 /* image is curved around cylinder, so FOV angle (in radians)
1248 * scales directly to image X coordinate, according to its radius.
1250 coeff[1] = (double) image->columns/coeff[0];
1252 /* radius is distance away from an image with this angular FOV */
1253 coeff[1] = (double) image->columns / ( 2 * tan(coeff[0]/2) );
1255 coeff[2] = (double)(image->columns)/2.0+image->page.x;
1256 coeff[3] = (double)(image->rows)/2.0+image->page.y;
1257 coeff[4] = coeff[2];
1258 coeff[5] = coeff[3]; /* assuming image size is the same */
1261 case BarrelDistortion:
1262 case BarrelInverseDistortion:
1264 /* Barrel Distortion
1265 Rs=(A*Rd^3 + B*Rd^2 + C*Rd + D)*Rd
1266 BarrelInv Distortion
1267 Rs=Rd/(A*Rd^3 + B*Rd^2 + C*Rd + D)
1269 Where Rd is the normalized radius from corner to middle of image
1270 Input Arguments are one of the following forms (number of arguments)...
1275 8: Ax,Bx,Cx,Dx Ay,By,Cy,Dy
1276 10: Ax,Bx,Cx,Dx Ay,By,Cy,Dy X,Y
1278 Returns 10 coefficent values, which are de-normalized (pixel scale)
1279 Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Xc, Yc
1281 /* Radius de-normalization scaling factor */
1283 rscale = 2.0/MagickMin((double) image->columns,(double) image->rows);
1285 /* sanity check number of args must = 3,4,5,6,8,10 or error */
1286 if ( (number_arguments < 3) || (number_arguments == 7) ||
1287 (number_arguments == 9) || (number_arguments > 10) )
1289 coeff=(double *) RelinquishMagickMemory(coeff);
1290 (void) ThrowMagickException(exception,GetMagickModule(),
1291 OptionError,"InvalidArgument", "%s : number of arguments",
1292 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1293 return((double *) NULL);
1295 /* A,B,C,D coefficients */
1296 coeff[0] = arguments[0];
1297 coeff[1] = arguments[1];
1298 coeff[2] = arguments[2];
1299 if ((number_arguments == 3) || (number_arguments == 5) )
1300 coeff[3] = 1.0 - coeff[0] - coeff[1] - coeff[2];
1302 coeff[3] = arguments[3];
1303 /* de-normalize the coefficients */
1304 coeff[0] *= pow(rscale,3.0);
1305 coeff[1] *= rscale*rscale;
1307 /* Y coefficients: as given OR same as X coefficients */
1308 if ( number_arguments >= 8 ) {
1309 coeff[4] = arguments[4] * pow(rscale,3.0);
1310 coeff[5] = arguments[5] * rscale*rscale;
1311 coeff[6] = arguments[6] * rscale;
1312 coeff[7] = arguments[7];
1315 coeff[4] = coeff[0];
1316 coeff[5] = coeff[1];
1317 coeff[6] = coeff[2];
1318 coeff[7] = coeff[3];
1320 /* X,Y Center of Distortion (image coodinates) */
1321 if ( number_arguments == 5 ) {
1322 coeff[8] = arguments[3];
1323 coeff[9] = arguments[4];
1325 else if ( number_arguments == 6 ) {
1326 coeff[8] = arguments[4];
1327 coeff[9] = arguments[5];
1329 else if ( number_arguments == 10 ) {
1330 coeff[8] = arguments[8];
1331 coeff[9] = arguments[9];
1334 /* center of the image provided (image coodinates) */
1335 coeff[8] = (double)image->columns/2.0 + image->page.x;
1336 coeff[9] = (double)image->rows/2.0 + image->page.y;
1340 case ShepardsDistortion:
1342 /* Shepards Distortion input arguments are the coefficents!
1343 Just check the number of arguments is valid!
1344 Args: u1,v1, x1,y1, ...
1345 OR : u1,v1, r1,g1,c1, ...
1347 if ( number_arguments%cp_size != 0 ||
1348 number_arguments < cp_size ) {
1349 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1350 "InvalidArgument", "%s : 'require at least %.20g CPs'",
1351 CommandOptionToMnemonic(MagickDistortOptions, *method), 1.0);
1352 coeff=(double *) RelinquishMagickMemory(coeff);
1353 return((double *) NULL);
1360 /* you should never reach this point */
1361 assert(! "No Method Handler"); /* just fail assertion */
1362 return((double *) NULL);
1366 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1370 + D i s t o r t R e s i z e I m a g e %
1374 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1376 % DistortResizeImage() resize image using the equivalent but slower image
1377 % distortion operator. The filter is applied using a EWA cylindrical
1378 % resampling. But like resize the final image size is limited to whole pixels
1379 % with no effects by virtual-pixels on the result.
1381 % Note that images containing a transparency channel will be twice as slow to
1382 % resize as images one without transparency.
1384 % The format of the DistortResizeImage method is:
1386 % Image *AdaptiveResizeImage(const Image *image,const size_t columns,
1387 % const size_t rows,ExceptionInfo *exception)
1389 % A description of each parameter follows:
1391 % o image: the image.
1393 % o columns: the number of columns in the resized image.
1395 % o rows: the number of rows in the resized image.
1397 % o exception: return any errors or warnings in this structure.
1400 MagickExport Image *DistortResizeImage(const Image *image,
1401 const size_t columns,const size_t rows,ExceptionInfo *exception)
1403 #define DistortResizeImageTag "Distort/Image"
1419 Distort resize image.
1421 assert(image != (const Image *) NULL);
1422 assert(image->signature == MagickSignature);
1423 if (image->debug != MagickFalse)
1424 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1425 assert(exception != (ExceptionInfo *) NULL);
1426 assert(exception->signature == MagickSignature);
1427 if ((columns == 0) || (rows == 0))
1428 return((Image *) NULL);
1429 /* Do not short-circuit this resize if final image size is unchanged */
1431 (void) SetImageVirtualPixelMethod(image,TransparentVirtualPixelMethod);
1433 (void) ResetMagickMemory(distort_args,0,12*sizeof(double));
1434 distort_args[4]=(double) image->columns;
1435 distort_args[6]=(double) columns;
1436 distort_args[9]=(double) image->rows;
1437 distort_args[11]=(double) rows;
1439 vp_save=GetImageVirtualPixelMethod(image);
1441 tmp_image=CloneImage(image,0,0,MagickTrue,exception);
1442 if ( tmp_image == (Image *) NULL )
1443 return((Image *) NULL);
1444 (void) SetImageVirtualPixelMethod(tmp_image,TransparentVirtualPixelMethod);
1446 if (image->matte == MagickFalse)
1449 Image has not transparency channel, so we free to use it
1451 (void) SetImageAlphaChannel(tmp_image,SetAlphaChannel,exception);
1452 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1453 MagickTrue,exception),
1455 tmp_image=DestroyImage(tmp_image);
1456 if ( resize_image == (Image *) NULL )
1457 return((Image *) NULL);
1459 (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel,exception);
1470 Image has transparency so handle colors and alpha separatly.
1471 Basically we need to separate Virtual-Pixel alpha in the resized
1472 image, so only the actual original images alpha channel is used.
1474 distort alpha channel separately
1476 channel_mask=SetPixelChannelMask(tmp_image,AlphaChannel);
1477 (void) SeparateImage(tmp_image);
1478 SetPixelChannelMap(tmp_image,channel_mask);
1479 (void) SetImageAlphaChannel(tmp_image,OpaqueAlphaChannel,exception);
1480 resize_alpha=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1481 MagickTrue,exception),
1482 tmp_image=DestroyImage(tmp_image);
1483 if (resize_alpha == (Image *) NULL)
1484 return((Image *) NULL);
1486 /* distort the actual image containing alpha + VP alpha */
1487 tmp_image=CloneImage(image,0,0,MagickTrue,exception);
1488 if ( tmp_image == (Image *) NULL )
1489 return((Image *) NULL);
1490 (void) SetImageVirtualPixelMethod(tmp_image,
1491 TransparentVirtualPixelMethod);
1492 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1493 MagickTrue,exception),
1494 tmp_image=DestroyImage(tmp_image);
1495 if ( resize_image == (Image *) NULL)
1497 resize_alpha=DestroyImage(resize_alpha);
1498 return((Image *) NULL);
1500 /* replace resize images alpha with the separally distorted alpha */
1501 (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel,exception);
1502 (void) SetImageAlphaChannel(resize_alpha,DeactivateAlphaChannel,exception);
1503 (void) CompositeImage(resize_image,CopyOpacityCompositeOp,resize_alpha,
1505 resize_alpha=DestroyImage(resize_alpha);
1507 (void) SetImageVirtualPixelMethod(resize_image,vp_save);
1510 Clean up the results of the Distortion
1512 crop_area.width=columns;
1513 crop_area.height=rows;
1517 tmp_image=resize_image;
1518 resize_image=CropImage(tmp_image,&crop_area,exception);
1519 tmp_image=DestroyImage(tmp_image);
1521 if ( resize_image == (Image *) NULL )
1522 return((Image *) NULL);
1524 return(resize_image);
1528 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1532 % D i s t o r t I m a g e %
1536 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1538 % DistortImage() distorts an image using various distortion methods, by
1539 % mapping color lookups of the source image to a new destination image
1540 % usally of the same size as the source image, unless 'bestfit' is set to
1543 % If 'bestfit' is enabled, and distortion allows it, the destination image is
1544 % adjusted to ensure the whole source 'image' will just fit within the final
1545 % destination image, which will be sized and offset accordingly. Also in
1546 % many cases the virtual offset of the source image will be taken into
1547 % account in the mapping.
1549 % If the '-verbose' control option has been set print to standard error the
1550 % equicelent '-fx' formula with coefficients for the function, if practical.
1552 % The format of the DistortImage() method is:
1554 % Image *DistortImage(const Image *image,const DistortImageMethod method,
1555 % const size_t number_arguments,const double *arguments,
1556 % MagickBooleanType bestfit, ExceptionInfo *exception)
1558 % A description of each parameter follows:
1560 % o image: the image to be distorted.
1562 % o method: the method of image distortion.
1564 % ArcDistortion always ignores source image offset, and always
1565 % 'bestfit' the destination image with the top left corner offset
1566 % relative to the polar mapping center.
1568 % Affine, Perspective, and Bilinear, do least squares fitting of the
1569 % distrotion when more than the minimum number of control point pairs
1572 % Perspective, and Bilinear, fall back to a Affine distortion when less
1573 % than 4 control point pairs are provided. While Affine distortions
1574 % let you use any number of control point pairs, that is Zero pairs is
1575 % a No-Op (viewport only) distortion, one pair is a translation and
1576 % two pairs of control points do a scale-rotate-translate, without any
1579 % o number_arguments: the number of arguments given.
1581 % o arguments: an array of floating point arguments for this method.
1583 % o bestfit: Attempt to 'bestfit' the size of the resulting image.
1584 % This also forces the resulting image to be a 'layered' virtual
1585 % canvas image. Can be overridden using 'distort:viewport' setting.
1587 % o exception: return any errors or warnings in this structure
1589 % Extra Controls from Image meta-data (artifacts)...
1592 % Output to stderr alternatives, internal coefficents, and FX
1593 % equivalents for the distortion operation (if feasible).
1594 % This forms an extra check of the distortion method, and allows users
1595 % access to the internal constants IM calculates for the distortion.
1597 % o "distort:viewport"
1598 % Directly set the output image canvas area and offest to use for the
1599 % resulting image, rather than use the original images canvas, or a
1600 % calculated 'bestfit' canvas.
1603 % Scale the size of the output canvas by this amount to provide a
1604 % method of Zooming, and for super-sampling the results.
1606 % Other settings that can effect results include
1608 % o 'interpolate' For source image lookups (scale enlargements)
1610 % o 'filter' Set filter to use for area-resampling (scale shrinking).
1611 % Set to 'point' to turn off and use 'interpolate' lookup
1616 MagickExport Image *DistortImage(const Image *image,DistortImageMethod method,
1617 const size_t number_arguments,const double *arguments,
1618 MagickBooleanType bestfit,ExceptionInfo *exception)
1620 #define DistortImageTag "Distort/Image"
1630 geometry; /* geometry of the distorted space viewport */
1635 assert(image != (Image *) NULL);
1636 assert(image->signature == MagickSignature);
1637 if (image->debug != MagickFalse)
1638 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1639 assert(exception != (ExceptionInfo *) NULL);
1640 assert(exception->signature == MagickSignature);
1644 Handle Special Compound Distortions
1646 if ( method == ResizeDistortion )
1648 if ( number_arguments != 2 )
1650 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1651 "InvalidArgument","%s : '%s'","Resize",
1652 "Invalid number of args: 2 only");
1653 return((Image *) NULL);
1655 distort_image=DistortResizeImage(image,(size_t)arguments[0],
1656 (size_t)arguments[1], exception);
1657 return(distort_image);
1661 Convert input arguments (usually as control points for reverse mapping)
1662 into mapping coefficients to apply the distortion.
1664 Note that some distortions are mapped to other distortions,
1665 and as such do not require specific code after this point.
1667 coeff = GenerateCoefficients(image, &method, number_arguments,
1668 arguments, 0, exception);
1669 if ( coeff == (double *) NULL )
1670 return((Image *) NULL);
1673 Determine the size and offset for a 'bestfit' destination.
1674 Usally the four corners of the source image is enough.
1677 /* default output image bounds, when no 'bestfit' is requested */
1678 geometry.width=image->columns;
1679 geometry.height=image->rows;
1683 if ( method == ArcDistortion ) {
1684 bestfit = MagickTrue; /* always calculate a 'best fit' viewport */
1687 /* Work out the 'best fit', (required for ArcDistortion) */
1690 s,d,min,max; /* source, dest coords --mapping--> min, max coords */
1693 fix_bounds = MagickTrue; /* enlarge bounds for VP handling */
1695 s.x=s.y=min.x=max.x=min.y=max.y=0.0; /* keep compiler happy */
1697 /* defines to figure out the bounds of the distorted image */
1698 #define InitalBounds(p) \
1700 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1701 min.x = max.x = p.x; \
1702 min.y = max.y = p.y; \
1704 #define ExpandBounds(p) \
1706 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1707 min.x = MagickMin(min.x,p.x); \
1708 max.x = MagickMax(max.x,p.x); \
1709 min.y = MagickMin(min.y,p.y); \
1710 max.y = MagickMax(max.y,p.y); \
1715 case AffineDistortion:
1716 { double inverse[6];
1717 InvertAffineCoefficients(coeff, inverse);
1718 s.x = (double) image->page.x;
1719 s.y = (double) image->page.y;
1720 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1721 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1723 s.x = (double) image->page.x+image->columns;
1724 s.y = (double) image->page.y;
1725 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1726 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1728 s.x = (double) image->page.x;
1729 s.y = (double) image->page.y+image->rows;
1730 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1731 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1733 s.x = (double) image->page.x+image->columns;
1734 s.y = (double) image->page.y+image->rows;
1735 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1736 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1740 case PerspectiveDistortion:
1741 { double inverse[8], scale;
1742 InvertPerspectiveCoefficients(coeff, inverse);
1743 s.x = (double) image->page.x;
1744 s.y = (double) image->page.y;
1745 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1746 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1747 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1748 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1750 s.x = (double) image->page.x+image->columns;
1751 s.y = (double) image->page.y;
1752 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1753 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1754 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1755 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1757 s.x = (double) image->page.x;
1758 s.y = (double) image->page.y+image->rows;
1759 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1760 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1761 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1762 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1764 s.x = (double) image->page.x+image->columns;
1765 s.y = (double) image->page.y+image->rows;
1766 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1767 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1768 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1769 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1775 /* Forward Map Corners */
1776 a = coeff[0]-coeff[1]/2; ca = cos(a); sa = sin(a);
1780 d.x = (coeff[2]-coeff[3])*ca;
1781 d.y = (coeff[2]-coeff[3])*sa;
1783 a = coeff[0]+coeff[1]/2; ca = cos(a); sa = sin(a);
1787 d.x = (coeff[2]-coeff[3])*ca;
1788 d.y = (coeff[2]-coeff[3])*sa;
1790 /* Orthogonal points along top of arc */
1791 for( a=(double) (ceil((double) ((coeff[0]-coeff[1]/2.0)/MagickPI2))*MagickPI2);
1792 a<(coeff[0]+coeff[1]/2.0); a+=MagickPI2 ) {
1793 ca = cos(a); sa = sin(a);
1799 Convert the angle_to_width and radius_to_height
1800 to appropriate scaling factors, to allow faster processing
1801 in the mapping function.
1803 coeff[1] = (double) (Magick2PI*image->columns/coeff[1]);
1804 coeff[3] = (double)image->rows/coeff[3];
1807 case PolarDistortion:
1809 if (number_arguments < 2)
1810 coeff[2] = coeff[3] = 0.0;
1811 min.x = coeff[2]-coeff[0];
1812 max.x = coeff[2]+coeff[0];
1813 min.y = coeff[3]-coeff[0];
1814 max.y = coeff[3]+coeff[0];
1815 /* should be about 1.0 if Rmin = 0 */
1816 coeff[7]=(double) geometry.height/(coeff[0]-coeff[1]);
1819 case DePolarDistortion:
1821 /* direct calculation as it needs to tile correctly
1822 * for reversibility in a DePolar-Polar cycle */
1823 fix_bounds = MagickFalse;
1824 geometry.x = geometry.y = 0;
1825 geometry.height = (size_t) ceil(coeff[0]-coeff[1]);
1826 geometry.width = (size_t)
1827 ceil((coeff[0]-coeff[1])*(coeff[5]-coeff[4])*0.5);
1828 /* correct scaling factors relative to new size */
1829 coeff[6]=(coeff[5]-coeff[4])/geometry.width; /* changed width */
1830 coeff[7]=(coeff[0]-coeff[1])/geometry.height; /* should be about 1.0 */
1833 case Cylinder2PlaneDistortion:
1835 /* direct calculation so center of distortion is either a pixel
1836 * center, or pixel edge. This allows for reversibility of the
1838 geometry.x = geometry.y = 0;
1839 geometry.width = (size_t) ceil( 2.0*coeff[1]*tan(coeff[0]/2.0) );
1840 geometry.height = (size_t) ceil( 2.0*coeff[3]/cos(coeff[0]/2.0) );
1841 /* correct center of distortion relative to new size */
1842 coeff[4] = (double) geometry.width/2.0;
1843 coeff[5] = (double) geometry.height/2.0;
1844 fix_bounds = MagickFalse;
1847 case Plane2CylinderDistortion:
1849 /* direct calculation center is either pixel center, or pixel edge
1850 * so as to allow reversibility of the image distortion */
1851 geometry.x = geometry.y = 0;
1852 geometry.width = (size_t) ceil(coeff[0]*coeff[1]); /* FOV * radius */
1853 geometry.height = (size_t) (2*coeff[3]); /* input image height */
1854 /* correct center of distortion relative to new size */
1855 coeff[4] = (double) geometry.width/2.0;
1856 coeff[5] = (double) geometry.height/2.0;
1857 fix_bounds = MagickFalse;
1860 case ShepardsDistortion:
1861 case BilinearForwardDistortion:
1862 case BilinearReverseDistortion:
1864 case QuadrilateralDistortion:
1866 case PolynomialDistortion:
1867 case BarrelDistortion:
1868 case BarrelInverseDistortion:
1870 /* no calculated bestfit available for these distortions */
1871 bestfit = MagickFalse;
1872 fix_bounds = MagickFalse;
1876 /* Set the output image geometry to calculated 'bestfit'.
1877 Yes this tends to 'over do' the file image size, ON PURPOSE!
1878 Do not do this for DePolar which needs to be exact for virtual tiling.
1881 geometry.x = (ssize_t) floor(min.x-0.5);
1882 geometry.y = (ssize_t) floor(min.y-0.5);
1883 geometry.width=(size_t) ceil(max.x-geometry.x+0.5);
1884 geometry.height=(size_t) ceil(max.y-geometry.y+0.5);
1887 } /* end bestfit destination image calculations */
1889 /* The user provided a 'viewport' expert option which may
1890 overrides some parts of the current output image geometry.
1891 This also overrides its default 'bestfit' setting.
1893 { const char *artifact=GetImageArtifact(image,"distort:viewport");
1894 viewport_given = MagickFalse;
1895 if ( artifact != (const char *) NULL ) {
1896 (void) ParseAbsoluteGeometry(artifact,&geometry);
1897 viewport_given = MagickTrue;
1901 /* Verbose output */
1902 if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
1905 char image_gen[MaxTextExtent];
1908 /* Set destination image size and virtual offset */
1909 if ( bestfit || viewport_given ) {
1910 (void) FormatLocaleString(image_gen, MaxTextExtent," -size %.20gx%.20g "
1911 "-page %+.20g%+.20g xc: +insert \\\n",(double) geometry.width,
1912 (double) geometry.height,(double) geometry.x,(double) geometry.y);
1913 lookup="v.p{ xx-v.page.x-.5, yy-v.page.y-.5 }";
1916 image_gen[0] = '\0'; /* no destination to generate */
1917 lookup = "p{ xx-page.x-.5, yy-page.y-.5 }"; /* simplify lookup */
1921 case AffineDistortion:
1925 inverse = (double *) AcquireQuantumMemory(6,sizeof(*inverse));
1926 if (inverse == (double *) NULL) {
1927 coeff = (double *) RelinquishMagickMemory(coeff);
1928 (void) ThrowMagickException(exception,GetMagickModule(),
1929 ResourceLimitError,"MemoryAllocationFailed",
1930 "%s", "DistortImages");
1931 return((Image *) NULL);
1933 InvertAffineCoefficients(coeff, inverse);
1934 CoefficientsToAffineArgs(inverse);
1935 (void) FormatLocaleFile(stderr, "Affine Projection:\n");
1936 (void) FormatLocaleFile(stderr, " -distort AffineProjection \\\n '");
1937 for (i=0; i < 5; i++)
1938 (void) FormatLocaleFile(stderr, "%lf,", inverse[i]);
1939 (void) FormatLocaleFile(stderr, "%lf'\n", inverse[5]);
1940 inverse = (double *) RelinquishMagickMemory(inverse);
1942 (void) FormatLocaleFile(stderr, "Affine Distort, FX Equivelent:\n");
1943 (void) FormatLocaleFile(stderr, "%s", image_gen);
1944 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1945 (void) FormatLocaleFile(stderr, " xx=%+lf*ii %+lf*jj %+lf;\n",
1946 coeff[0], coeff[1], coeff[2]);
1947 (void) FormatLocaleFile(stderr, " yy=%+lf*ii %+lf*jj %+lf;\n",
1948 coeff[3], coeff[4], coeff[5]);
1949 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
1954 case PerspectiveDistortion:
1958 inverse = (double *) AcquireQuantumMemory(8,sizeof(*inverse));
1959 if (inverse == (double *) NULL) {
1960 coeff = (double *) RelinquishMagickMemory(coeff);
1961 (void) ThrowMagickException(exception,GetMagickModule(),
1962 ResourceLimitError,"MemoryAllocationFailed",
1963 "%s", "DistortCoefficients");
1964 return((Image *) NULL);
1966 InvertPerspectiveCoefficients(coeff, inverse);
1967 (void) FormatLocaleFile(stderr, "Perspective Projection:\n");
1968 (void) FormatLocaleFile(stderr, " -distort PerspectiveProjection \\\n '");
1970 (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
1971 (void) FormatLocaleFile(stderr, "\n ");
1973 (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
1974 (void) FormatLocaleFile(stderr, "%lf'\n", inverse[7]);
1975 inverse = (double *) RelinquishMagickMemory(inverse);
1977 (void) FormatLocaleFile(stderr, "Perspective Distort, FX Equivelent:\n");
1978 (void) FormatLocaleFile(stderr, "%s", image_gen);
1979 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1980 (void) FormatLocaleFile(stderr, " rr=%+lf*ii %+lf*jj + 1;\n",
1981 coeff[6], coeff[7]);
1982 (void) FormatLocaleFile(stderr, " xx=(%+lf*ii %+lf*jj %+lf)/rr;\n",
1983 coeff[0], coeff[1], coeff[2]);
1984 (void) FormatLocaleFile(stderr, " yy=(%+lf*ii %+lf*jj %+lf)/rr;\n",
1985 coeff[3], coeff[4], coeff[5]);
1986 (void) FormatLocaleFile(stderr, " rr%s0 ? %s : blue' \\\n",
1987 coeff[8] < 0 ? "<" : ">", lookup);
1991 case BilinearForwardDistortion:
1992 (void) FormatLocaleFile(stderr, "BilinearForward Mapping Equations:\n");
1993 (void) FormatLocaleFile(stderr, "%s", image_gen);
1994 (void) FormatLocaleFile(stderr, " i = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
1995 coeff[0], coeff[1], coeff[2], coeff[3]);
1996 (void) FormatLocaleFile(stderr, " j = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
1997 coeff[4], coeff[5], coeff[6], coeff[7]);
2000 (void) FormatLocaleFile(stderr, " c8 = %+lf c9 = 2*a = %+lf;\n",
2001 coeff[8], coeff[9]);
2003 (void) FormatLocaleFile(stderr, "BilinearForward Distort, FX Equivelent:\n");
2004 (void) FormatLocaleFile(stderr, "%s", image_gen);
2005 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
2006 0.5-coeff[3], 0.5-coeff[7]);
2007 (void) FormatLocaleFile(stderr, " bb=%lf*ii %+lf*jj %+lf;\n",
2008 coeff[6], -coeff[2], coeff[8]);
2009 /* Handle Special degenerate (non-quadratic) or trapezoidal case */
2010 if ( coeff[9] != 0 ) {
2011 (void) FormatLocaleFile(stderr, " rt=bb*bb %+lf*(%lf*ii%+lf*jj);\n",
2012 -2*coeff[9], coeff[4], -coeff[0]);
2013 (void) FormatLocaleFile(stderr, " yy=( -bb + sqrt(rt) ) / %lf;\n",
2016 (void) FormatLocaleFile(stderr, " yy=(%lf*ii%+lf*jj)/bb;\n",
2017 -coeff[4], coeff[0]);
2018 (void) FormatLocaleFile(stderr, " xx=(ii %+lf*yy)/(%lf %+lf*yy);\n",
2019 -coeff[1], coeff[0], coeff[2]);
2020 if ( coeff[9] != 0 )
2021 (void) FormatLocaleFile(stderr, " (rt < 0 ) ? red : %s'\n", lookup);
2023 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2026 case BilinearReverseDistortion:
2028 (void) FormatLocaleFile(stderr, "Polynomial Projection Distort:\n");
2029 (void) FormatLocaleFile(stderr, " -distort PolynomialProjection \\\n");
2030 (void) FormatLocaleFile(stderr, " '1.5, %lf, %lf, %lf, %lf,\n",
2031 coeff[3], coeff[0], coeff[1], coeff[2]);
2032 (void) FormatLocaleFile(stderr, " %lf, %lf, %lf, %lf'\n",
2033 coeff[7], coeff[4], coeff[5], coeff[6]);
2035 (void) FormatLocaleFile(stderr, "BilinearReverse Distort, FX Equivelent:\n");
2036 (void) FormatLocaleFile(stderr, "%s", image_gen);
2037 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2038 (void) FormatLocaleFile(stderr, " xx=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
2039 coeff[0], coeff[1], coeff[2], coeff[3]);
2040 (void) FormatLocaleFile(stderr, " yy=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
2041 coeff[4], coeff[5], coeff[6], coeff[7]);
2042 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2045 case PolynomialDistortion:
2047 size_t nterms = (size_t) coeff[1];
2048 (void) FormatLocaleFile(stderr, "Polynomial (order %lg, terms %lu), FX Equivelent\n",
2049 coeff[0],(unsigned long) nterms);
2050 (void) FormatLocaleFile(stderr, "%s", image_gen);
2051 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2052 (void) FormatLocaleFile(stderr, " xx =");
2053 for (i=0; i<(ssize_t) nterms; i++) {
2054 if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n ");
2055 (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i],
2058 (void) FormatLocaleFile(stderr, ";\n yy =");
2059 for (i=0; i<(ssize_t) nterms; i++) {
2060 if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n ");
2061 (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i+nterms],
2064 (void) FormatLocaleFile(stderr, ";\n %s' \\\n", lookup);
2069 (void) FormatLocaleFile(stderr, "Arc Distort, Internal Coefficients:\n");
2070 for ( i=0; i<5; i++ )
2071 (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
2072 (void) FormatLocaleFile(stderr, "Arc Distort, FX Equivelent:\n");
2073 (void) FormatLocaleFile(stderr, "%s", image_gen);
2074 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x; jj=j+page.y;\n");
2075 (void) FormatLocaleFile(stderr, " xx=(atan2(jj,ii)%+lf)/(2*pi);\n",
2077 (void) FormatLocaleFile(stderr, " xx=xx-round(xx);\n");
2078 (void) FormatLocaleFile(stderr, " xx=xx*%lf %+lf;\n",
2079 coeff[1], coeff[4]);
2080 (void) FormatLocaleFile(stderr, " yy=(%lf - hypot(ii,jj)) * %lf;\n",
2081 coeff[2], coeff[3]);
2082 (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n");
2085 case PolarDistortion:
2087 (void) FormatLocaleFile(stderr, "Polar Distort, Internal Coefficents\n");
2088 for ( i=0; i<8; i++ )
2089 (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
2090 (void) FormatLocaleFile(stderr, "Polar Distort, FX Equivelent:\n");
2091 (void) FormatLocaleFile(stderr, "%s", image_gen);
2092 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
2093 -coeff[2], -coeff[3]);
2094 (void) FormatLocaleFile(stderr, " xx=(atan2(ii,jj)%+lf)/(2*pi);\n",
2095 -(coeff[4]+coeff[5])/2 );
2096 (void) FormatLocaleFile(stderr, " xx=xx-round(xx);\n");
2097 (void) FormatLocaleFile(stderr, " xx=xx*2*pi*%lf + v.w/2;\n",
2099 (void) FormatLocaleFile(stderr, " yy=(hypot(ii,jj)%+lf)*%lf;\n",
2100 -coeff[1], coeff[7] );
2101 (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n");
2104 case DePolarDistortion:
2106 (void) FormatLocaleFile(stderr, "DePolar Distort, Internal Coefficents\n");
2107 for ( i=0; i<8; i++ )
2108 (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
2109 (void) FormatLocaleFile(stderr, "DePolar Distort, FX Equivelent:\n");
2110 (void) FormatLocaleFile(stderr, "%s", image_gen);
2111 (void) FormatLocaleFile(stderr, " -fx 'aa=(i+.5)*%lf %+lf;\n", coeff[6], -coeff[4] );
2112 (void) FormatLocaleFile(stderr, " rr=(j+.5)*%lf %+lf;\n", coeff[7], +coeff[1] );
2113 (void) FormatLocaleFile(stderr, " xx=rr*sin(aa) %+lf;\n", coeff[2] );
2114 (void) FormatLocaleFile(stderr, " yy=rr*cos(aa) %+lf;\n", coeff[3] );
2115 (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n");
2118 case Cylinder2PlaneDistortion:
2120 (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, Internal Coefficents\n");
2121 (void) FormatLocaleFile(stderr, " cylinder_radius = %+lf\n", coeff[1]);
2122 (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, FX Equivelent:\n");
2123 (void) FormatLocaleFile(stderr, "%s", image_gen);
2124 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
2125 -coeff[4], -coeff[5]);
2126 (void) FormatLocaleFile(stderr, " aa=atan(ii/%+lf);\n", coeff[1] );
2127 (void) FormatLocaleFile(stderr, " xx=%lf*aa%+lf;\n",
2128 coeff[1], coeff[2] );
2129 (void) FormatLocaleFile(stderr, " yy=jj*cos(aa)%+lf;\n", coeff[3] );
2130 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2133 case Plane2CylinderDistortion:
2135 (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, Internal Coefficents\n");
2136 (void) FormatLocaleFile(stderr, " cylinder_radius = %+lf\n", coeff[1]);
2137 (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, FX Equivelent:\n");
2138 (void) FormatLocaleFile(stderr, "%s", image_gen);
2139 (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
2140 -coeff[4], -coeff[5]);
2141 (void) FormatLocaleFile(stderr, " ii=ii/%+lf;\n", coeff[1] );
2142 (void) FormatLocaleFile(stderr, " xx=%lf*tan(ii)%+lf;\n",
2143 coeff[1], coeff[2] );
2144 (void) FormatLocaleFile(stderr, " yy=jj/cos(ii)%+lf;\n",
2146 (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2149 case BarrelDistortion:
2150 case BarrelInverseDistortion:
2152 /* NOTE: This does the barrel roll in pixel coords not image coords
2153 ** The internal distortion must do it in image coordinates,
2154 ** so that is what the center coeff (8,9) is given in.
2156 xc = ((double)image->columns-1.0)/2.0 + image->page.x;
2157 yc = ((double)image->rows-1.0)/2.0 + image->page.y;
2158 (void) FormatLocaleFile(stderr, "Barrel%s Distort, FX Equivelent:\n",
2159 method == BarrelDistortion ? "" : "Inv");
2160 (void) FormatLocaleFile(stderr, "%s", image_gen);
2161 if ( fabs(coeff[8]-xc-0.5) < 0.1 && fabs(coeff[9]-yc-0.5) < 0.1 )
2162 (void) FormatLocaleFile(stderr, " -fx 'xc=(w-1)/2; yc=(h-1)/2;\n");
2164 (void) FormatLocaleFile(stderr, " -fx 'xc=%lf; yc=%lf;\n",
2165 coeff[8]-0.5, coeff[9]-0.5);
2166 (void) FormatLocaleFile(stderr,
2167 " ii=i-xc; jj=j-yc; rr=hypot(ii,jj);\n");
2168 (void) FormatLocaleFile(stderr, " ii=ii%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2169 method == BarrelDistortion ? "*" : "/",
2170 coeff[0],coeff[1],coeff[2],coeff[3]);
2171 (void) FormatLocaleFile(stderr, " jj=jj%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2172 method == BarrelDistortion ? "*" : "/",
2173 coeff[4],coeff[5],coeff[6],coeff[7]);
2174 (void) FormatLocaleFile(stderr, " v.p{fx*ii+xc,fy*jj+yc}' \\\n");
2181 /* The user provided a 'scale' expert option will scale the
2182 output image size, by the factor given allowing for super-sampling
2183 of the distorted image space. Any scaling factors must naturally
2184 be halved as a result.
2186 { const char *artifact;
2187 artifact=GetImageArtifact(image,"distort:scale");
2188 output_scaling = 1.0;
2189 if (artifact != (const char *) NULL) {
2190 output_scaling = fabs(InterpretLocaleValue(artifact,(char **) NULL));
2191 geometry.width *= (size_t) output_scaling;
2192 geometry.height *= (size_t) output_scaling;
2193 geometry.x *= (ssize_t) output_scaling;
2194 geometry.y *= (ssize_t) output_scaling;
2195 if ( output_scaling < 0.1 ) {
2196 coeff = (double *) RelinquishMagickMemory(coeff);
2197 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
2198 "InvalidArgument","%s", "-set option:distort:scale" );
2199 return((Image *) NULL);
2201 output_scaling = 1/output_scaling;
2204 #define ScaleFilter(F,A,B,C,D) \
2205 ScaleResampleFilter( (F), \
2206 output_scaling*(A), output_scaling*(B), \
2207 output_scaling*(C), output_scaling*(D) )
2210 Initialize the distort image attributes.
2212 distort_image=CloneImage(image,geometry.width,geometry.height,MagickTrue,
2214 if (distort_image == (Image *) NULL)
2215 return((Image *) NULL);
2216 /* if image is ColorMapped - change it to DirectClass */
2217 if (SetImageStorageClass(distort_image,DirectClass,exception) == MagickFalse)
2219 distort_image=DestroyImage(distort_image);
2220 return((Image *) NULL);
2222 distort_image->page.x=geometry.x;
2223 distort_image->page.y=geometry.y;
2224 if (distort_image->background_color.alpha != OpaqueAlpha)
2225 distort_image->matte=MagickTrue;
2227 { /* ----- MAIN CODE -----
2228 Sample the source image to each pixel in the distort image.
2243 **restrict resample_filter;
2250 GetPixelInfo(distort_image,&zero);
2251 resample_filter=AcquireResampleFilterThreadSet(image,
2252 UndefinedVirtualPixelMethod,MagickFalse,exception);
2253 distort_view=AcquireCacheView(distort_image);
2254 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2255 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2257 for (j=0; j < (ssize_t) distort_image->rows; j++)
2260 id = GetOpenMPThreadId();
2263 validity; /* how mathematically valid is this the mapping */
2269 pixel, /* pixel color to assign to distorted image */
2270 invalid; /* the color to assign when distort result is invalid */
2274 s; /* transform destination image x,y to source image x,y */
2282 q=QueueCacheViewAuthenticPixels(distort_view,0,j,distort_image->columns,1,
2284 if (q == (Quantum *) NULL)
2291 /* Define constant scaling vectors for Affine Distortions
2292 Other methods are either variable, or use interpolated lookup
2296 case AffineDistortion:
2297 ScaleFilter( resample_filter[id],
2299 coeff[3], coeff[4] );
2305 /* Initialize default pixel validity
2306 * negative: pixel is invalid output 'matte_color'
2307 * 0.0 to 1.0: antialiased, mix with resample output
2308 * 1.0 or greater: use resampled output.
2312 GetPixelInfo(distort_image,&invalid);
2313 SetPixelInfoPacket(distort_image,&distort_image->matte_color,&invalid);
2314 if (distort_image->colorspace == CMYKColorspace)
2315 ConvertRGBToCMYK(&invalid); /* what about other color spaces? */
2317 for (i=0; i < (ssize_t) distort_image->columns; i++)
2319 /* map pixel coordinate to distortion space coordinate */
2320 d.x = (double) (geometry.x+i+0.5)*output_scaling;
2321 d.y = (double) (geometry.y+j+0.5)*output_scaling;
2322 s = d; /* default is a no-op mapping */
2325 case AffineDistortion:
2327 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2328 s.y=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2329 /* Affine partial derivitives are constant -- set above */
2332 case PerspectiveDistortion:
2335 p,q,r,abs_r,abs_c6,abs_c7,scale;
2336 /* perspective is a ratio of affines */
2337 p=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2338 q=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2339 r=coeff[6]*d.x+coeff[7]*d.y+1.0;
2340 /* Pixel Validity -- is it a 'sky' or 'ground' pixel */
2341 validity = (r*coeff[8] < 0.0) ? 0.0 : 1.0;
2342 /* Determine horizon anti-alias blending */
2344 abs_c6 = fabs(coeff[6]);
2345 abs_c7 = fabs(coeff[7]);
2346 if ( abs_c6 > abs_c7 ) {
2347 if ( abs_r < abs_c6*output_scaling )
2348 validity = 0.5 - coeff[8]*r/(coeff[6]*output_scaling);
2350 else if ( abs_r < abs_c7*output_scaling )
2351 validity = 0.5 - coeff[8]*r/(coeff[7]*output_scaling);
2352 /* Perspective Sampling Point (if valid) */
2353 if ( validity > 0.0 ) {
2354 /* divide by r affine, for perspective scaling */
2358 /* Perspective Partial Derivatives or Scaling Vectors */
2360 ScaleFilter( resample_filter[id],
2361 (r*coeff[0] - p*coeff[6])*scale,
2362 (r*coeff[1] - p*coeff[7])*scale,
2363 (r*coeff[3] - q*coeff[6])*scale,
2364 (r*coeff[4] - q*coeff[7])*scale );
2368 case BilinearReverseDistortion:
2370 /* Reversed Mapped is just a simple polynomial */
2371 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2]*d.x*d.y+coeff[3];
2372 s.y=coeff[4]*d.x+coeff[5]*d.y
2373 +coeff[6]*d.x*d.y+coeff[7];
2374 /* Bilinear partial derivitives of scaling vectors */
2375 ScaleFilter( resample_filter[id],
2376 coeff[0] + coeff[2]*d.y,
2377 coeff[1] + coeff[2]*d.x,
2378 coeff[4] + coeff[6]*d.y,
2379 coeff[5] + coeff[6]*d.x );
2382 case BilinearForwardDistortion:
2384 /* Forward mapped needs reversed polynomial equations
2385 * which unfortunatally requires a square root! */
2387 d.x -= coeff[3]; d.y -= coeff[7];
2388 b = coeff[6]*d.x - coeff[2]*d.y + coeff[8];
2389 c = coeff[4]*d.x - coeff[0]*d.y;
2392 /* Handle Special degenerate (non-quadratic) case
2393 * Currently without horizon anti-alising */
2394 if ( fabs(coeff[9]) < MagickEpsilon )
2397 c = b*b - 2*coeff[9]*c;
2401 s.y = ( -b + sqrt(c) )/coeff[9];
2403 if ( validity > 0.0 )
2404 s.x = ( d.x - coeff[1]*s.y) / ( coeff[0] + coeff[2]*s.y );
2406 /* NOTE: the sign of the square root should be -ve for parts
2407 where the source image becomes 'flipped' or 'mirrored'.
2408 FUTURE: Horizon handling
2409 FUTURE: Scaling factors or Deritives (how?)
2414 case BilinearDistortion:
2415 /* Bilinear mapping of any Quadrilateral to any Quadrilateral */
2416 /* UNDER DEVELOPMENT */
2419 case PolynomialDistortion:
2421 /* multi-ordered polynomial */
2426 nterms=(ssize_t)coeff[1];
2429 du,dv; /* the du,dv vectors from unit dx,dy -- derivatives */
2431 s.x=s.y=du.x=du.y=dv.x=dv.y=0.0;
2432 for(k=0; k < nterms; k++) {
2433 s.x += poly_basis_fn(k,d.x,d.y)*coeff[2+k];
2434 du.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k];
2435 du.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k];
2436 s.y += poly_basis_fn(k,d.x,d.y)*coeff[2+k+nterms];
2437 dv.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k+nterms];
2438 dv.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k+nterms];
2440 ScaleFilter( resample_filter[id], du.x,du.y,dv.x,dv.y );
2445 /* what is the angle and radius in the destination image */
2446 s.x = (double) ((atan2(d.y,d.x) - coeff[0])/Magick2PI);
2447 s.x -= MagickRound(s.x); /* angle */
2448 s.y = hypot(d.x,d.y); /* radius */
2450 /* Arc Distortion Partial Scaling Vectors
2451 Are derived by mapping the perpendicular unit vectors
2452 dR and dA*R*2PI rather than trying to map dx and dy
2453 The results is a very simple orthogonal aligned ellipse.
2455 if ( s.y > MagickEpsilon )
2456 ScaleFilter( resample_filter[id],
2457 (double) (coeff[1]/(Magick2PI*s.y)), 0, 0, coeff[3] );
2459 ScaleFilter( resample_filter[id],
2460 distort_image->columns*2, 0, 0, coeff[3] );
2462 /* now scale the angle and radius for source image lookup point */
2463 s.x = s.x*coeff[1] + coeff[4] + image->page.x +0.5;
2464 s.y = (coeff[2] - s.y) * coeff[3] + image->page.y;
2467 case PolarDistortion:
2468 { /* 2D Cartesain to Polar View */
2471 s.x = atan2(d.x,d.y) - (coeff[4]+coeff[5])/2;
2473 s.x -= MagickRound(s.x);
2474 s.x *= Magick2PI; /* angle - relative to centerline */
2475 s.y = hypot(d.x,d.y); /* radius */
2477 /* Polar Scaling vectors are based on mapping dR and dA vectors
2478 This results in very simple orthogonal scaling vectors
2480 if ( s.y > MagickEpsilon )
2481 ScaleFilter( resample_filter[id],
2482 (double) (coeff[6]/(Magick2PI*s.y)), 0, 0, coeff[7] );
2484 ScaleFilter( resample_filter[id],
2485 distort_image->columns*2, 0, 0, coeff[7] );
2487 /* now finish mapping radius/angle to source x,y coords */
2488 s.x = s.x*coeff[6] + (double)image->columns/2.0 + image->page.x;
2489 s.y = (s.y-coeff[1])*coeff[7] + image->page.y;
2492 case DePolarDistortion:
2493 { /* @D Polar to Carteasain */
2494 /* ignore all destination virtual offsets */
2495 d.x = ((double)i+0.5)*output_scaling*coeff[6]-coeff[4];
2496 d.y = ((double)j+0.5)*output_scaling*coeff[7]+coeff[1];
2497 s.x = d.y*sin(d.x) + coeff[2];
2498 s.y = d.y*cos(d.x) + coeff[3];
2499 /* derivatives are usless - better to use SuperSampling */
2502 case Cylinder2PlaneDistortion:
2503 { /* 3D Cylinder to Tangential Plane */
2505 /* relative to center of distortion */
2506 d.x -= coeff[4]; d.y -= coeff[5];
2507 d.x /= coeff[1]; /* x' = x/r */
2508 ax=atan(d.x); /* aa = atan(x/r) = u/r */
2509 cx=cos(ax); /* cx = cos(atan(x/r)) = 1/sqrt(x^2+u^2) */
2510 s.x = coeff[1]*ax; /* u = r*atan(x/r) */
2511 s.y = d.y*cx; /* v = y*cos(u/r) */
2512 /* derivatives... (see personnal notes) */
2513 ScaleFilter( resample_filter[id],
2514 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2516 if ( i == 0 && j == 0 ) {
2517 fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2518 fprintf(stderr, "phi = %lf\n", (double)(ax * 180.0/MagickPI) );
2519 fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n",
2520 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2523 /* add center of distortion in source */
2524 s.x += coeff[2]; s.y += coeff[3];
2527 case Plane2CylinderDistortion:
2528 { /* 3D Cylinder to Tangential Plane */
2529 /* relative to center of distortion */
2530 d.x -= coeff[4]; d.y -= coeff[5];
2532 /* is pixel valid - horizon of a infinite Virtual-Pixel Plane
2533 * (see Anthony Thyssen's personal note) */
2534 validity = (double) (coeff[1]*MagickPI2 - fabs(d.x))/output_scaling + 0.5;
2536 if ( validity > 0.0 ) {
2538 d.x /= coeff[1]; /* x'= x/r */
2539 cx = 1/cos(d.x); /* cx = 1/cos(x/r) */
2540 tx = tan(d.x); /* tx = tan(x/r) */
2541 s.x = coeff[1]*tx; /* u = r * tan(x/r) */
2542 s.y = d.y*cx; /* v = y / cos(x/r) */
2543 /* derivatives... (see Anthony Thyssen's personal notes) */
2544 ScaleFilter( resample_filter[id],
2545 cx*cx, 0.0, s.y*cx/coeff[1], cx );
2547 /*if ( i == 0 && j == 0 )*/
2548 if ( d.x == 0.5 && d.y == 0.5 ) {
2549 fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2550 fprintf(stderr, "radius = %lf phi = %lf validity = %lf\n",
2551 coeff[1], (double)(d.x * 180.0/MagickPI), validity );
2552 fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n",
2553 cx*cx, 0.0, s.y*cx/coeff[1], cx);
2557 /* add center of distortion in source */
2558 s.x += coeff[2]; s.y += coeff[3];
2561 case BarrelDistortion:
2562 case BarrelInverseDistortion:
2563 { /* Lens Barrel Distionion Correction */
2564 double r,fx,fy,gx,gy;
2565 /* Radial Polynomial Distortion (de-normalized) */
2568 r = sqrt(d.x*d.x+d.y*d.y);
2569 if ( r > MagickEpsilon ) {
2570 fx = ((coeff[0]*r + coeff[1])*r + coeff[2])*r + coeff[3];
2571 fy = ((coeff[4]*r + coeff[5])*r + coeff[6])*r + coeff[7];
2572 gx = ((3*coeff[0]*r + 2*coeff[1])*r + coeff[2])/r;
2573 gy = ((3*coeff[4]*r + 2*coeff[5])*r + coeff[6])/r;
2574 /* adjust functions and scaling for 'inverse' form */
2575 if ( method == BarrelInverseDistortion ) {
2576 fx = 1/fx; fy = 1/fy;
2577 gx *= -fx*fx; gy *= -fy*fy;
2579 /* Set the source pixel to lookup and EWA derivative vectors */
2580 s.x = d.x*fx + coeff[8];
2581 s.y = d.y*fy + coeff[9];
2582 ScaleFilter( resample_filter[id],
2583 gx*d.x*d.x + fx, gx*d.x*d.y,
2584 gy*d.x*d.y, gy*d.y*d.y + fy );
2587 /* Special handling to avoid divide by zero when r==0
2589 ** The source and destination pixels match in this case
2590 ** which was set at the top of the loop using s = d;
2591 ** otherwise... s.x=coeff[8]; s.y=coeff[9];
2593 if ( method == BarrelDistortion )
2594 ScaleFilter( resample_filter[id],
2595 coeff[3], 0, 0, coeff[7] );
2596 else /* method == BarrelInverseDistortion */
2597 /* FUTURE, trap for D==0 causing division by zero */
2598 ScaleFilter( resample_filter[id],
2599 1.0/coeff[3], 0, 0, 1.0/coeff[7] );
2603 case ShepardsDistortion:
2604 { /* Shepards Method, or Inverse Weighted Distance for
2605 displacement around the destination image control points
2606 The input arguments are the coefficents to the function.
2607 This is more of a 'displacement' function rather than an
2608 absolute distortion function.
2615 denominator = s.x = s.y = 0;
2616 for(i=0; i<number_arguments; i+=4) {
2618 ((double)d.x-arguments[i+2])*((double)d.x-arguments[i+2])
2619 + ((double)d.y-arguments[i+3])*((double)d.y-arguments[i+3]);
2625 s.x += (arguments[ i ]-arguments[i+2])*weight;
2626 s.y += (arguments[i+1]-arguments[i+3])*weight;
2627 denominator += weight;
2634 /* We can not determine derivatives using shepards method
2635 only color interpolatation, not area-resampling */
2639 break; /* use the default no-op given above */
2641 /* map virtual canvas location back to real image coordinate */
2642 if ( bestfit && method != ArcDistortion ) {
2643 s.x -= image->page.x;
2644 s.y -= image->page.y;
2649 if ( validity <= 0.0 ) {
2650 /* result of distortion is an invalid pixel - don't resample */
2651 SetPixelPixelInfo(distort_image,&invalid,q);
2654 /* resample the source image to find its correct color */
2655 (void) ResamplePixelColor(resample_filter[id],s.x,s.y,&pixel);
2656 /* if validity between 0.0 and 1.0 mix result with invalid pixel */
2657 if ( validity < 1.0 ) {
2658 /* Do a blend of sample color and invalid pixel */
2659 /* should this be a 'Blend', or an 'Over' compose */
2660 CompositePixelInfoBlend(&pixel,validity,&invalid,(1.0-validity),
2663 SetPixelPixelInfo(distort_image,&pixel,q);
2665 q+=GetPixelChannels(distort_image);
2667 sync=SyncCacheViewAuthenticPixels(distort_view,exception);
2668 if (sync == MagickFalse)
2670 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2675 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2676 #pragma omp critical (MagickCore_DistortImage)
2678 proceed=SetImageProgress(image,DistortImageTag,progress++,
2680 if (proceed == MagickFalse)
2684 distort_view=DestroyCacheView(distort_view);
2685 resample_filter=DestroyResampleFilterThreadSet(resample_filter);
2687 if (status == MagickFalse)
2688 distort_image=DestroyImage(distort_image);
2691 /* Arc does not return an offset unless 'bestfit' is in effect
2692 And the user has not provided an overriding 'viewport'.
2694 if ( method == ArcDistortion && !bestfit && !viewport_given ) {
2695 distort_image->page.x = 0;
2696 distort_image->page.y = 0;
2698 coeff = (double *) RelinquishMagickMemory(coeff);
2699 return(distort_image);
2703 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2707 % S p a r s e C o l o r I m a g e %
2711 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2713 % SparseColorImage(), given a set of coordinates, interpolates the colors
2714 % found at those coordinates, across the whole image, using various methods.
2716 % The format of the SparseColorImage() method is:
2718 % Image *SparseColorImage(const Image *image,
2719 % const SparseColorMethod method,const size_t number_arguments,
2720 % const double *arguments,ExceptionInfo *exception)
2722 % A description of each parameter follows:
2724 % o image: the image to be filled in.
2726 % o method: the method to fill in the gradient between the control points.
2728 % The methods used for SparseColor() are often simular to methods
2729 % used for DistortImage(), and even share the same code for determination
2730 % of the function coefficents, though with more dimensions (or resulting
2733 % o number_arguments: the number of arguments given.
2735 % o arguments: array of floating point arguments for this method--
2736 % x,y,color_values-- with color_values given as normalized values.
2738 % o exception: return any errors or warnings in this structure
2741 MagickExport Image *SparseColorImage(const Image *image,
2742 const SparseColorMethod method,const size_t number_arguments,
2743 const double *arguments,ExceptionInfo *exception)
2745 #define SparseColorTag "Distort/SparseColor"
2759 assert(image != (Image *) NULL);
2760 assert(image->signature == MagickSignature);
2761 if (image->debug != MagickFalse)
2762 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2763 assert(exception != (ExceptionInfo *) NULL);
2764 assert(exception->signature == MagickSignature);
2766 /* Determine number of color values needed per control point */
2768 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2770 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2772 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2774 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2775 (image->colorspace == CMYKColorspace))
2777 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2778 (image->matte != MagickFalse))
2782 Convert input arguments into mapping coefficients, this this case
2783 we are mapping (distorting) colors, rather than coordinates.
2785 { DistortImageMethod
2788 distort_method=(DistortImageMethod) method;
2789 if ( distort_method >= SentinelDistortion )
2790 distort_method = ShepardsDistortion; /* Pretend to be Shepards */
2791 coeff = GenerateCoefficients(image, &distort_method, number_arguments,
2792 arguments, number_colors, exception);
2793 if ( coeff == (double *) NULL )
2794 return((Image *) NULL);
2796 Note some Distort Methods may fall back to other simpler methods,
2797 Currently the only fallback of concern is Bilinear to Affine
2798 (Barycentric), which is alaso sparse_colr method. This also ensures
2799 correct two and one color Barycentric handling.
2801 sparse_method = (SparseColorMethod) distort_method;
2802 if ( distort_method == ShepardsDistortion )
2803 sparse_method = method; /* return non-distiort methods to normal */
2806 /* Verbose output */
2807 if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
2809 switch (sparse_method) {
2810 case BarycentricColorInterpolate:
2812 register ssize_t x=0;
2813 (void) FormatLocaleFile(stderr, "Barycentric Sparse Color:\n");
2814 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2815 (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf' \\\n",
2816 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2817 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2818 (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf' \\\n",
2819 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2820 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2821 (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf' \\\n",
2822 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2823 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2824 (image->colorspace == CMYKColorspace))
2825 (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf' \\\n",
2826 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2827 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2828 (image->matte != MagickFalse))
2829 (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf' \\\n",
2830 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2833 case BilinearColorInterpolate:
2835 register ssize_t x=0;
2836 (void) FormatLocaleFile(stderr, "Bilinear Sparse Color\n");
2837 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2838 (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2839 coeff[ x ], coeff[x+1],
2840 coeff[x+2], coeff[x+3]),x+=4;
2841 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2842 (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2843 coeff[ x ], coeff[x+1],
2844 coeff[x+2], coeff[x+3]),x+=4;
2845 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2846 (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2847 coeff[ x ], coeff[x+1],
2848 coeff[x+2], coeff[x+3]),x+=4;
2849 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2850 (image->colorspace == CMYKColorspace))
2851 (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2852 coeff[ x ], coeff[x+1],
2853 coeff[x+2], coeff[x+3]),x+=4;
2854 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2855 (image->matte != MagickFalse))
2856 (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2857 coeff[ x ], coeff[x+1],
2858 coeff[x+2], coeff[x+3]),x+=4;
2862 /* sparse color method is too complex for FX emulation */
2867 /* Generate new image for generated interpolated gradient.
2868 * ASIDE: Actually we could have just replaced the colors of the original
2869 * image, but IM Core policy, is if storage class could change then clone
2873 sparse_image=CloneImage(image,0,0,MagickTrue,exception);
2874 if (sparse_image == (Image *) NULL)
2875 return((Image *) NULL);
2876 if (SetImageStorageClass(sparse_image,DirectClass,exception) == MagickFalse)
2877 { /* if image is ColorMapped - change it to DirectClass */
2878 sparse_image=DestroyImage(sparse_image);
2879 return((Image *) NULL);
2881 { /* ----- MAIN CODE ----- */
2896 sparse_view=AcquireCacheView(sparse_image);
2897 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2898 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2900 for (j=0; j < (ssize_t) sparse_image->rows; j++)
2906 pixel; /* pixel to assign to distorted image */
2914 q=GetCacheViewAuthenticPixels(sparse_view,0,j,sparse_image->columns,
2916 if (q == (Quantum *) NULL)
2921 GetPixelInfo(sparse_image,&pixel);
2922 for (i=0; i < (ssize_t) image->columns; i++)
2924 SetPixelInfo(image,q,&pixel);
2925 switch (sparse_method)
2927 case BarycentricColorInterpolate:
2929 register ssize_t x=0;
2930 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2931 pixel.red = coeff[x]*i +coeff[x+1]*j
2933 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2934 pixel.green = coeff[x]*i +coeff[x+1]*j
2936 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2937 pixel.blue = coeff[x]*i +coeff[x+1]*j
2939 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2940 (image->colorspace == CMYKColorspace))
2941 pixel.black = coeff[x]*i +coeff[x+1]*j
2943 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2944 (image->matte != MagickFalse))
2945 pixel.alpha = coeff[x]*i +coeff[x+1]*j
2949 case BilinearColorInterpolate:
2951 register ssize_t x=0;
2952 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2953 pixel.red = coeff[x]*i + coeff[x+1]*j +
2954 coeff[x+2]*i*j + coeff[x+3], x+=4;
2955 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2956 pixel.green = coeff[x]*i + coeff[x+1]*j +
2957 coeff[x+2]*i*j + coeff[x+3], x+=4;
2958 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2959 pixel.blue = coeff[x]*i + coeff[x+1]*j +
2960 coeff[x+2]*i*j + coeff[x+3], x+=4;
2961 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2962 (image->colorspace == CMYKColorspace))
2963 pixel.black = coeff[x]*i + coeff[x+1]*j +
2964 coeff[x+2]*i*j + coeff[x+3], x+=4;
2965 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2966 (image->matte != MagickFalse))
2967 pixel.alpha = coeff[x]*i + coeff[x+1]*j +
2968 coeff[x+2]*i*j + coeff[x+3], x+=4;
2971 case InverseColorInterpolate:
2972 case ShepardsColorInterpolate:
2973 { /* Inverse (Squared) Distance weights average (IDW) */
2979 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2981 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2983 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2985 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2986 (image->colorspace == CMYKColorspace))
2988 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2989 (image->matte != MagickFalse))
2992 for(k=0; k<number_arguments; k+=2+number_colors) {
2993 register ssize_t x=(ssize_t) k+2;
2995 ((double)i-arguments[ k ])*((double)i-arguments[ k ])
2996 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
2997 if ( method == InverseColorInterpolate )
2998 weight = sqrt(weight); /* inverse, not inverse squared */
2999 weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
3000 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3001 pixel.red += arguments[x++]*weight;
3002 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3003 pixel.green += arguments[x++]*weight;
3004 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3005 pixel.blue += arguments[x++]*weight;
3006 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3007 (image->colorspace == CMYKColorspace))
3008 pixel.black += arguments[x++]*weight;
3009 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3010 (image->matte != MagickFalse))
3011 pixel.alpha += arguments[x++]*weight;
3012 denominator += weight;
3014 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3015 pixel.red/=denominator;
3016 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3017 pixel.green/=denominator;
3018 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3019 pixel.blue/=denominator;
3020 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3021 (image->colorspace == CMYKColorspace))
3022 pixel.black/=denominator;
3023 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3024 (image->matte != MagickFalse))
3025 pixel.alpha/=denominator;
3028 case VoronoiColorInterpolate:
3030 { /* Just use the closest control point you can find! */
3034 minimum = MagickHuge;
3036 for(k=0; k<number_arguments; k+=2+number_colors) {
3038 ((double)i-arguments[ k ])*((double)i-arguments[ k ])
3039 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
3040 if ( distance < minimum ) {
3041 register ssize_t x=(ssize_t) k+2;
3042 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3043 pixel.red=arguments[x++];
3044 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3045 pixel.green=arguments[x++];
3046 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3047 pixel.blue=arguments[x++];
3048 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3049 (image->colorspace == CMYKColorspace))
3050 pixel.black=arguments[x++];
3051 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3052 (image->matte != MagickFalse))
3053 pixel.alpha=arguments[x++];
3060 /* set the color directly back into the source image */
3061 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3062 pixel.red*=QuantumRange;
3063 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3064 pixel.green*=QuantumRange;
3065 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3066 pixel.blue*=QuantumRange;
3067 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3068 (image->colorspace == CMYKColorspace))
3069 pixel.black*=QuantumRange;
3070 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3071 (image->matte != MagickFalse))
3072 pixel.alpha*=QuantumRange;
3073 SetPixelPixelInfo(sparse_image,&pixel,q);
3074 q+=GetPixelChannels(sparse_image);
3076 sync=SyncCacheViewAuthenticPixels(sparse_view,exception);
3077 if (sync == MagickFalse)
3079 if (image->progress_monitor != (MagickProgressMonitor) NULL)
3084 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3085 #pragma omp critical (MagickCore_SparseColorImage)
3087 proceed=SetImageProgress(image,SparseColorTag,progress++,image->rows);
3088 if (proceed == MagickFalse)
3092 sparse_view=DestroyCacheView(sparse_view);
3093 if (status == MagickFalse)
3094 sparse_image=DestroyImage(sparse_image);
3096 coeff = (double *) RelinquishMagickMemory(coeff);
3097 return(sparse_image);