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-2010 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 "magick/studio.h"
44 #include "magick/artifact.h"
45 #include "magick/cache.h"
46 #include "magick/cache-view.h"
47 #include "magick/colorspace-private.h"
48 #include "magick/composite-private.h"
49 #include "magick/distort.h"
50 #include "magick/exception.h"
51 #include "magick/exception-private.h"
52 #include "magick/gem.h"
53 #include "magick/hashmap.h"
54 #include "magick/image.h"
55 #include "magick/list.h"
56 #include "magick/matrix.h"
57 #include "magick/memory_.h"
58 #include "magick/monitor-private.h"
59 #include "magick/option.h"
60 #include "magick/pixel.h"
61 #include "magick/pixel-private.h"
62 #include "magick/resample.h"
63 #include "magick/resample-private.h"
64 #include "magick/registry.h"
65 #include "magick/semaphore.h"
66 #include "magick/string_.h"
67 #include "magick/string-private.h"
68 #include "magick/thread-private.h"
69 #include "magick/token.h"
72 Numerous internal routines for image distortions.
74 static inline double MagickMin(const double x,const double y)
76 return( x < y ? x : y);
78 static inline double MagickMax(const double x,const double y)
80 return( x > y ? x : y);
83 static inline void AffineArgsToCoefficients(double *affine)
85 /* map external sx,ry,rx,sy,tx,ty to internal c0,c2,c4,c1,c3,c5 */
86 double tmp[4]; /* note indexes 0 and 5 remain unchanged */
87 tmp[0]=affine[1]; tmp[1]=affine[2]; tmp[2]=affine[3]; tmp[3]=affine[4];
88 affine[3]=tmp[0]; affine[1]=tmp[1]; affine[4]=tmp[2]; affine[2]=tmp[3];
91 static inline void CoefficientsToAffineArgs(double *coeff)
93 /* map internal c0,c1,c2,c3,c4,c5 to external sx,ry,rx,sy,tx,ty */
94 double tmp[4]; /* note indexes 0 and 5 remain unchanged */
95 tmp[0]=coeff[3]; tmp[1]=coeff[1]; tmp[2]=coeff[4]; tmp[3]=coeff[2];
96 coeff[1]=tmp[0]; coeff[2]=tmp[1]; coeff[3]=tmp[2]; coeff[4]=tmp[3];
98 static void InvertAffineCoefficients(const double *coeff,double *inverse)
100 /* From "Digital Image Warping" by George Wolberg, page 50 */
103 determinant=1.0/(coeff[0]*coeff[4]-coeff[1]*coeff[3]);
104 inverse[0]=determinant*coeff[4];
105 inverse[1]=determinant*(-coeff[1]);
106 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[2]*coeff[4]);
107 inverse[3]=determinant*(-coeff[3]);
108 inverse[4]=determinant*coeff[0];
109 inverse[5]=determinant*(coeff[2]*coeff[3]-coeff[0]*coeff[5]);
112 static void InvertPerspectiveCoefficients(const double *coeff,
115 /* From "Digital Image Warping" by George Wolberg, page 53 */
118 determinant=1.0/(coeff[0]*coeff[4]-coeff[3]*coeff[1]);
119 inverse[0]=determinant*(coeff[4]-coeff[7]*coeff[5]);
120 inverse[1]=determinant*(coeff[7]*coeff[2]-coeff[1]);
121 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[4]*coeff[2]);
122 inverse[3]=determinant*(coeff[6]*coeff[5]-coeff[3]);
123 inverse[4]=determinant*(coeff[0]-coeff[6]*coeff[2]);
124 inverse[5]=determinant*(coeff[3]*coeff[2]-coeff[0]*coeff[5]);
125 inverse[6]=determinant*(coeff[3]*coeff[7]-coeff[6]*coeff[4]);
126 inverse[7]=determinant*(coeff[6]*coeff[1]-coeff[0]*coeff[7]);
129 static inline double MagickRound(double x)
132 Round the fraction to nearest integer.
135 return((double) ((ssize_t) (x+0.5)));
136 return((double) ((ssize_t) (x-0.5)));
140 * Polynomial Term Defining Functions
142 * Order must either be an integer, or 1.5 to produce
143 * the 2 number_valuesal polyminal function...
144 * affine 1 (3) u = c0 + c1*x + c2*y
145 * bilinear 1.5 (4) u = '' + c3*x*y
146 * quadratic 2 (6) u = '' + c4*x*x + c5*y*y
147 * cubic 3 (10) u = '' + c6*x^3 + c7*x*x*y + c8*x*y*y + c9*y^3
148 * quartic 4 (15) u = '' + c10*x^4 + ... + c14*y^4
149 * quintic 5 (21) u = '' + c15*x^5 + ... + c20*y^5
150 * number in parenthesis minimum number of points needed.
151 * Anything beyond quintic, has not been implemented until
152 * a more automated way of determined terms is found.
154 * Note the slight re-ordering of the terms for a quadratic polynomial
155 * which is to allow the use of a bi-linear (order=1.5) polynomial.
156 * All the later polynomials are ordered simply from x^N to y^N
158 static size_t poly_number_terms(double order)
160 /* Return the number of terms for a 2d polynomial */
161 if ( order < 1 || order > 5 ||
162 ( order != floor(order) && (order-1.5) > MagickEpsilon) )
163 return 0; /* invalid polynomial order */
164 return((size_t) floor((order+1)*(order+2)/2));
167 static double poly_basis_fn(ssize_t n, double x, double y)
169 /* Return the result for this polynomial term */
171 case 0: return( 1.0 ); /* constant */
173 case 2: return( y ); /* affine order = 1 terms = 3 */
174 case 3: return( x*y ); /* bilinear order = 1.5 terms = 4 */
175 case 4: return( x*x );
176 case 5: return( y*y ); /* quadratic order = 2 terms = 6 */
177 case 6: return( x*x*x );
178 case 7: return( x*x*y );
179 case 8: return( x*y*y );
180 case 9: return( y*y*y ); /* cubic order = 3 terms = 10 */
181 case 10: return( x*x*x*x );
182 case 11: return( x*x*x*y );
183 case 12: return( x*x*y*y );
184 case 13: return( x*y*y*y );
185 case 14: return( y*y*y*y ); /* quartic order = 4 terms = 15 */
186 case 15: return( x*x*x*x*x );
187 case 16: return( x*x*x*x*y );
188 case 17: return( x*x*x*y*y );
189 case 18: return( x*x*y*y*y );
190 case 19: return( x*y*y*y*y );
191 case 20: return( y*y*y*y*y ); /* quintic order = 5 terms = 21 */
193 return( 0 ); /* should never happen */
195 static const char *poly_basis_str(ssize_t n)
197 /* return the result for this polynomial term */
199 case 0: return(""); /* constant */
200 case 1: return("*ii");
201 case 2: return("*jj"); /* affiine order = 1 terms = 3 */
202 case 3: return("*ii*jj"); /* biiliinear order = 1.5 terms = 4 */
203 case 4: return("*ii*ii");
204 case 5: return("*jj*jj"); /* quadratiic order = 2 terms = 6 */
205 case 6: return("*ii*ii*ii");
206 case 7: return("*ii*ii*jj");
207 case 8: return("*ii*jj*jj");
208 case 9: return("*jj*jj*jj"); /* cubiic order = 3 terms = 10 */
209 case 10: return("*ii*ii*ii*ii");
210 case 11: return("*ii*ii*ii*jj");
211 case 12: return("*ii*ii*jj*jj");
212 case 13: return("*ii*jj*jj*jj");
213 case 14: return("*jj*jj*jj*jj"); /* quartiic order = 4 terms = 15 */
214 case 15: return("*ii*ii*ii*ii*ii");
215 case 16: return("*ii*ii*ii*ii*jj");
216 case 17: return("*ii*ii*ii*jj*jj");
217 case 18: return("*ii*ii*jj*jj*jj");
218 case 19: return("*ii*jj*jj*jj*jj");
219 case 20: return("*jj*jj*jj*jj*jj"); /* quiintiic order = 5 terms = 21 */
221 return( "UNKNOWN" ); /* should never happen */
223 static double poly_basis_dx(ssize_t n, double x, double y)
225 /* polynomial term for x derivative */
227 case 0: return( 0.0 ); /* constant */
228 case 1: return( 1.0 );
229 case 2: return( 0.0 ); /* affine order = 1 terms = 3 */
230 case 3: return( y ); /* bilinear order = 1.5 terms = 4 */
232 case 5: return( 0.0 ); /* quadratic order = 2 terms = 6 */
233 case 6: return( x*x );
234 case 7: return( x*y );
235 case 8: return( y*y );
236 case 9: return( 0.0 ); /* cubic order = 3 terms = 10 */
237 case 10: return( x*x*x );
238 case 11: return( x*x*y );
239 case 12: return( x*y*y );
240 case 13: return( y*y*y );
241 case 14: return( 0.0 ); /* quartic order = 4 terms = 15 */
242 case 15: return( x*x*x*x );
243 case 16: return( x*x*x*y );
244 case 17: return( x*x*y*y );
245 case 18: return( x*y*y*y );
246 case 19: return( y*y*y*y );
247 case 20: return( 0.0 ); /* quintic order = 5 terms = 21 */
249 return( 0.0 ); /* should never happen */
251 static double poly_basis_dy(ssize_t n, double x, double y)
253 /* polynomial term for y derivative */
255 case 0: return( 0.0 ); /* constant */
256 case 1: return( 0.0 );
257 case 2: return( 1.0 ); /* affine order = 1 terms = 3 */
258 case 3: return( x ); /* bilinear order = 1.5 terms = 4 */
259 case 4: return( 0.0 );
260 case 5: return( y ); /* quadratic order = 2 terms = 6 */
261 default: return( poly_basis_dx(n-1,x,y) ); /* weird but true */
263 /* NOTE: the only reason that last is not true for 'quadtratic'
264 is due to the re-arrangement of terms to allow for 'bilinear'
269 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
273 + G e n e r a t e C o e f f i c i e n t s %
277 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
279 % GenerateCoefficients() takes user provided input arguments and generates
280 % the coefficients, needed to apply the specific distortion for either
281 % distorting images (generally using control points) or generating a color
282 % gradient from sparsely separated color points.
284 % The format of the GenerateCoefficients() method is:
286 % Image *GenerateCoefficients(const Image *image,DistortImageMethod method,
287 % const size_t number_arguments,const double *arguments,
288 % size_t number_values, ExceptionInfo *exception)
290 % A description of each parameter follows:
292 % o image: the image to be distorted.
294 % o method: the method of image distortion/ sparse gradient
296 % o number_arguments: the number of arguments given.
298 % o arguments: the arguments for this distortion method.
300 % o number_values: the style and format of given control points, (caller type)
301 % 0: 2 dimensional mapping of control points (Distort)
302 % Format: u,v,x,y where u,v is the 'source' of the
303 % the color to be plotted, for DistortImage()
304 % N: Interpolation of control points with N values (usally r,g,b)
305 % Format: x,y,r,g,b mapping x,y to color values r,g,b
306 % IN future, varible number of values may be given (1 to N)
308 % o exception: return any errors or warnings in this structure
310 % Note that the returned array of double values must be freed by the
311 % calling method using RelinquishMagickMemory(). This however may change in
312 % the future to require a more 'method' specific method.
314 % Because of this this method should not be classed as stable or used
315 % outside other MagickCore library methods.
318 static double *GenerateCoefficients(const Image *image,
319 DistortImageMethod *method,const size_t number_arguments,
320 const double *arguments,size_t number_values,ExceptionInfo *exception)
329 number_coeff, /* number of coefficients to return (array size) */
330 cp_size, /* number floating point numbers per control point */
331 cp_x,cp_y, /* the x,y indexes for control point */
332 cp_values; /* index of values for this control point */
333 /* number_values Number of values given per control point */
335 if ( number_values == 0 ) {
336 /* Image distortion using control points (or other distortion)
337 That is generate a mapping so that x,y->u,v given u,v,x,y
339 number_values = 2; /* special case: two values of u,v */
340 cp_values = 0; /* the values i,j are BEFORE the destination CP x,y */
341 cp_x = 2; /* location of x,y in input control values */
343 /* NOTE: cp_values, also used for later 'reverse map distort' tests */
346 cp_x = 0; /* location of x,y in input control values */
348 cp_values = 2; /* and the other values are after x,y */
349 /* Typically in this case the values are R,G,B color values */
351 cp_size = number_values+2; /* each CP defintion involves this many numbers */
353 /* If not enough control point pairs are found for specific distortions
354 fall back to Affine distortion (allowing 0 to 3 point pairs)
356 if ( number_arguments < 4*cp_size &&
357 ( *method == BilinearForwardDistortion
358 || *method == BilinearReverseDistortion
359 || *method == PerspectiveDistortion
361 *method = AffineDistortion;
365 case AffineDistortion:
366 /* also BarycentricColorInterpolate: */
367 number_coeff=3*number_values;
369 case PolynomialDistortion:
370 /* number of coefficents depend on the given polynomal 'order' */
371 if ( number_arguments <= 1 && (number_arguments-1)%cp_size != 0)
373 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
374 "InvalidArgument","%s : '%s'","Polynomial",
375 "Invalid number of args: order [CPs]...");
376 return((double *) NULL);
378 i = poly_number_terms(arguments[0]);
379 number_coeff = 2 + i*number_values;
381 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
382 "InvalidArgument","%s : '%s'","Polynomial",
383 "Invalid order, should be interger 1 to 5, or 1.5");
384 return((double *) NULL);
386 if ( number_arguments < 1+i*cp_size ) {
387 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
388 "InvalidArgument", "%s : 'require at least %.20g CPs'",
389 "Polynomial", (double) i);
390 return((double *) NULL);
393 case BilinearReverseDistortion:
394 number_coeff=4*number_values;
397 The rest are constants as they are only used for image distorts
399 case BilinearForwardDistortion:
400 number_coeff=10; /* 2*4 coeff plus 2 constants */
401 cp_x = 0; /* Reverse src/dest coords for forward mapping */
406 case QuadraterialDistortion:
407 number_coeff=19; /* BilinearForward + BilinearReverse */
410 case ShepardsDistortion:
411 case VoronoiColorInterpolate:
412 number_coeff=1; /* not used, but provide some type of return */
417 case ScaleRotateTranslateDistortion:
418 case AffineProjectionDistortion:
421 case PolarDistortion:
422 case DePolarDistortion:
425 case PerspectiveDistortion:
426 case PerspectiveProjectionDistortion:
429 case BarrelDistortion:
430 case BarrelInverseDistortion:
433 case UndefinedDistortion:
435 assert(! "Unknown Method Given"); /* just fail assertion */
438 /* allocate the array of coefficients needed */
439 coeff = (double *) AcquireQuantumMemory(number_coeff,sizeof(*coeff));
440 if (coeff == (double *) NULL) {
441 (void) ThrowMagickException(exception,GetMagickModule(),
442 ResourceLimitError,"MemoryAllocationFailed",
443 "%s", "GenerateCoefficients");
444 return((double *) NULL);
447 /* zero out coeffiecents array */
448 for (i=0; i < number_coeff; i++)
453 case AffineDistortion:
457 for each 'value' given
459 Input Arguments are sets of control points...
460 For Distort Images u,v, x,y ...
461 For Sparse Gradients x,y, r,g,b ...
463 if ( number_arguments%cp_size != 0 ||
464 number_arguments < cp_size ) {
465 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
466 "InvalidArgument", "%s : 'require at least %.20g CPs'",
468 coeff=(double *) RelinquishMagickMemory(coeff);
469 return((double *) NULL);
471 /* handle special cases of not enough arguments */
472 if ( number_arguments == cp_size ) {
473 /* Only 1 CP Set Given */
474 if ( cp_values == 0 ) {
475 /* image distortion - translate the image */
477 coeff[2] = arguments[0] - arguments[2];
479 coeff[5] = arguments[1] - arguments[3];
482 /* sparse gradient - use the values directly */
483 for (i=0; i<number_values; i++)
484 coeff[i*3+2] = arguments[cp_values+i];
488 /* 2 or more points (usally 3) given.
489 Solve a least squares simultanious equation for coefficients.
499 /* create matrix, and a fake vectors matrix */
500 matrix = AcquireMagickMatrix(3UL,3UL);
501 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
502 if (matrix == (double **) NULL || vectors == (double **) NULL)
504 matrix = RelinquishMagickMatrix(matrix, 3UL);
505 vectors = (double **) RelinquishMagickMemory(vectors);
506 coeff = (double *) RelinquishMagickMemory(coeff);
507 (void) ThrowMagickException(exception,GetMagickModule(),
508 ResourceLimitError,"MemoryAllocationFailed",
509 "%s", "DistortCoefficients");
510 return((double *) NULL);
512 /* fake a number_values x3 vectors matrix from coefficients array */
513 for (i=0; i < number_values; i++)
514 vectors[i] = &(coeff[i*3]);
515 /* Add given control point pairs for least squares solving */
516 for (i=0; i < number_arguments; i+=cp_size) {
517 terms[0] = arguments[i+cp_x]; /* x */
518 terms[1] = arguments[i+cp_y]; /* y */
519 terms[2] = 1; /* 1 */
520 LeastSquaresAddTerms(matrix,vectors,terms,
521 &(arguments[i+cp_values]),3UL,number_values);
523 if ( number_arguments == 2*cp_size ) {
524 /* Only two pairs were given, but we need 3 to solve the affine.
525 Fake extra coordinates by rotating p1 around p0 by 90 degrees.
526 x2 = x0 - (y1-y0) y2 = y0 + (x1-x0)
528 terms[0] = arguments[cp_x]
529 - ( arguments[cp_size+cp_y] - arguments[cp_y] ); /* x2 */
530 terms[1] = arguments[cp_y] +
531 + ( arguments[cp_size+cp_x] - arguments[cp_x] ); /* y2 */
532 terms[2] = 1; /* 1 */
533 if ( cp_values == 0 ) {
534 /* Image Distortion - rotate the u,v coordients too */
537 uv2[0] = arguments[0] - arguments[5] + arguments[1]; /* u2 */
538 uv2[1] = arguments[1] + arguments[4] - arguments[0]; /* v2 */
539 LeastSquaresAddTerms(matrix,vectors,terms,uv2,3UL,2UL);
542 /* Sparse Gradient - use values of p0 for linear gradient */
543 LeastSquaresAddTerms(matrix,vectors,terms,
544 &(arguments[cp_values]),3UL,number_values);
547 /* Solve for LeastSquares Coefficients */
548 status=GaussJordanElimination(matrix,vectors,3UL,number_values);
549 matrix = RelinquishMagickMatrix(matrix, 3UL);
550 vectors = (double **) RelinquishMagickMemory(vectors);
551 if ( status == MagickFalse ) {
552 coeff = (double *) RelinquishMagickMemory(coeff);
553 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
554 "InvalidArgument","%s : 'Unsolvable Matrix'",
555 MagickOptionToMnemonic(MagickDistortOptions, *method) );
556 return((double *) NULL);
561 case AffineProjectionDistortion:
564 Arguments: Affine Matrix (forward mapping)
565 Arguments sx, rx, ry, sy, tx, ty
566 Where u = sx*x + ry*y + tx
569 Returns coefficients (in there inverse form) ordered as...
572 AffineProjection Distortion Notes...
573 + Will only work with a 2 number_values for Image Distortion
574 + Can not be used for generating a sparse gradient (interpolation)
577 if (number_arguments != 6) {
578 coeff = (double *) RelinquishMagickMemory(coeff);
579 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
580 "InvalidArgument","%s : 'Needs 6 coeff values'",
581 MagickOptionToMnemonic(MagickDistortOptions, *method) );
582 return((double *) NULL);
584 /* FUTURE: trap test for sx*sy-rx*ry == 0 (determinate = 0, no inverse) */
585 for(i=0; i<6UL; i++ )
586 inverse[i] = arguments[i];
587 AffineArgsToCoefficients(inverse); /* map into coefficents */
588 InvertAffineCoefficients(inverse, coeff); /* invert */
589 *method = AffineDistortion;
593 case ScaleRotateTranslateDistortion:
595 /* Scale, Rotate and Translate Distortion
596 An alturnative Affine Distortion
597 Argument options, by number of arguments given:
598 7: x,y, sx,sy, a, nx,ny
605 Where actions are (in order of application)
606 x,y 'center' of transforms (default = image center)
607 sx,sy scale image by this amount (default = 1)
608 a angle of rotation (argument required)
609 nx,ny move 'center' here (default = x,y or no movement)
610 And convert to affine mapping coefficients
612 ScaleRotateTranslate Distortion Notes...
613 + Does not use a set of CPs in any normal way
614 + Will only work with a 2 number_valuesal Image Distortion
615 + Can not be used for generating a sparse gradient (interpolation)
621 /* set default center, and default scale */
622 x = nx = (double)(image->columns)/2.0 + (double)image->page.x;
623 y = ny = (double)(image->rows)/2.0 + (double)image->page.y;
625 switch ( number_arguments ) {
627 coeff = (double *) RelinquishMagickMemory(coeff);
628 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
629 "InvalidArgument","%s : 'Needs at least 1 argument'",
630 MagickOptionToMnemonic(MagickDistortOptions, *method) );
631 return((double *) NULL);
636 sx = sy = arguments[0];
640 x = nx = arguments[0];
641 y = ny = arguments[1];
642 switch ( number_arguments ) {
647 sx = sy = arguments[2];
656 sx = sy = arguments[2];
669 coeff = (double *) RelinquishMagickMemory(coeff);
670 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
671 "InvalidArgument","%s : 'Too Many Arguments (7 or less)'",
672 MagickOptionToMnemonic(MagickDistortOptions, *method) );
673 return((double *) NULL);
677 /* Trap if sx or sy == 0 -- image is scaled out of existance! */
678 if ( fabs(sx) < MagickEpsilon || fabs(sy) < MagickEpsilon ) {
679 coeff = (double *) RelinquishMagickMemory(coeff);
680 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
681 "InvalidArgument","%s : 'Zero Scale Given'",
682 MagickOptionToMnemonic(MagickDistortOptions, *method) );
683 return((double *) NULL);
685 /* Save the given arguments as an affine distortion */
686 a=DegreesToRadians(a); cosine=cos(a); sine=sin(a);
688 *method = AffineDistortion;
691 coeff[2]=x-nx*coeff[0]-ny*coeff[1];
694 coeff[5]=y-nx*coeff[3]-ny*coeff[4];
697 case PerspectiveDistortion:
699 Perspective Distortion (a ratio of affine distortions)
701 p(x,y) c0*x + c1*y + c2
702 u = ------ = ------------------
703 r(x,y) c6*x + c7*y + 1
705 q(x,y) c3*x + c4*y + c5
706 v = ------ = ------------------
707 r(x,y) c6*x + c7*y + 1
709 c8 = Sign of 'r', or the denominator affine, for the actual image.
710 This determines what part of the distorted image is 'ground'
711 side of the horizon, the other part is 'sky' or invalid.
712 Valid values are +1.0 or -1.0 only.
714 Input Arguments are sets of control points...
715 For Distort Images u,v, x,y ...
716 For Sparse Gradients x,y, r,g,b ...
718 Perspective Distortion Notes...
719 + Can be thought of as ratio of 3 affine transformations
720 + Not separatable: r() or c6 and c7 are used by both equations
721 + All 8 coefficients must be determined simultaniously
722 + Will only work with a 2 number_valuesal Image Distortion
723 + Can not be used for generating a sparse gradient (interpolation)
724 + It is not linear, but is simple to generate an inverse
725 + All lines within an image remain lines.
726 + but distances between points may vary.
740 if ( number_arguments%cp_size != 0 ||
741 number_arguments < cp_size*4 ) {
742 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
743 "InvalidArgument", "%s : 'require at least %.20g CPs'",
744 MagickOptionToMnemonic(MagickDistortOptions, *method), 4.0);
745 coeff=(double *) RelinquishMagickMemory(coeff);
746 return((double *) NULL);
748 /* fake 1x8 vectors matrix directly using the coefficients array */
749 vectors[0] = &(coeff[0]);
750 /* 8x8 least-squares matrix (zeroed) */
751 matrix = AcquireMagickMatrix(8UL,8UL);
752 if (matrix == (double **) NULL) {
753 (void) ThrowMagickException(exception,GetMagickModule(),
754 ResourceLimitError,"MemoryAllocationFailed",
755 "%s", "DistortCoefficients");
756 return((double *) NULL);
758 /* Add control points for least squares solving */
759 for (i=0; i < number_arguments; i+=4) {
760 terms[0]=arguments[i+cp_x]; /* c0*x */
761 terms[1]=arguments[i+cp_y]; /* c1*y */
762 terms[2]=1.0; /* c2*1 */
766 terms[6]=-terms[0]*arguments[i+cp_u]; /* 1/(c6*x) */
767 terms[7]=-terms[1]*arguments[i+cp_u]; /* 1/(c7*y) */
768 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_u]),
774 terms[3]=arguments[i+cp_x]; /* c3*x */
775 terms[4]=arguments[i+cp_y]; /* c4*y */
776 terms[5]=1.0; /* c5*1 */
777 terms[6]=-terms[3]*arguments[i+cp_v]; /* 1/(c6*x) */
778 terms[7]=-terms[4]*arguments[i+cp_v]; /* 1/(c7*y) */
779 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_v]),
782 /* Solve for LeastSquares Coefficients */
783 status=GaussJordanElimination(matrix,vectors,8UL,1UL);
784 matrix = RelinquishMagickMatrix(matrix, 8UL);
785 if ( status == MagickFalse ) {
786 coeff = (double *) RelinquishMagickMemory(coeff);
787 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
788 "InvalidArgument","%s : 'Unsolvable Matrix'",
789 MagickOptionToMnemonic(MagickDistortOptions, *method) );
790 return((double *) NULL);
793 Calculate 9'th coefficient! The ground-sky determination.
794 What is sign of the 'ground' in r() denominator affine function?
795 Just use any valid image coordinate (first control point) in
796 destination for determination of what part of view is 'ground'.
798 coeff[8] = coeff[6]*arguments[cp_x]
799 + coeff[7]*arguments[cp_y] + 1.0;
800 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
804 case PerspectiveProjectionDistortion:
807 Arguments: Perspective Coefficents (forward mapping)
809 if (number_arguments != 8) {
810 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
811 "InvalidArgument", "%s : 'Needs 8 coefficient values'",
812 MagickOptionToMnemonic(MagickDistortOptions, *method));
813 return((double *) NULL);
815 /* FUTURE: trap test c0*c4-c3*c1 == 0 (determinate = 0, no inverse) */
816 InvertPerspectiveCoefficients(arguments, coeff);
818 Calculate 9'th coefficient! The ground-sky determination.
819 What is sign of the 'ground' in r() denominator affine function?
820 Just use any valid image cocodinate in destination for determination.
821 For a forward mapped perspective the images 0,0 coord will map to
822 c2,c5 in the distorted image, so set the sign of denominator of that.
824 coeff[8] = coeff[6]*arguments[2]
825 + coeff[7]*arguments[5] + 1.0;
826 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
827 *method = PerspectiveDistortion;
831 case BilinearForwardDistortion:
832 case BilinearReverseDistortion:
834 /* Bilinear Distortion (Forward mapping)
835 v = c0*x + c1*y + c2*x*y + c3;
836 for each 'value' given
838 This is actually a simple polynomial Distortion! The difference
839 however is when we need to reverse the above equation to generate a
840 BilinearForwardDistortion (see below).
842 Input Arguments are sets of control points...
843 For Distort Images u,v, x,y ...
844 For Sparse Gradients x,y, r,g,b ...
855 /* check the number of arguments */
856 if ( number_arguments%cp_size != 0 ||
857 number_arguments < cp_size*4 ) {
858 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
859 "InvalidArgument", "%s : 'require at least %.20g CPs'",
860 MagickOptionToMnemonic(MagickDistortOptions, *method), 4.0);
861 coeff=(double *) RelinquishMagickMemory(coeff);
862 return((double *) NULL);
864 /* create matrix, and a fake vectors matrix */
865 matrix = AcquireMagickMatrix(4UL,4UL);
866 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
867 if (matrix == (double **) NULL || vectors == (double **) NULL)
869 matrix = RelinquishMagickMatrix(matrix, 4UL);
870 vectors = (double **) RelinquishMagickMemory(vectors);
871 coeff = (double *) RelinquishMagickMemory(coeff);
872 (void) ThrowMagickException(exception,GetMagickModule(),
873 ResourceLimitError,"MemoryAllocationFailed",
874 "%s", "DistortCoefficients");
875 return((double *) NULL);
877 /* fake a number_values x4 vectors matrix from coefficients array */
878 for (i=0; i < number_values; i++)
879 vectors[i] = &(coeff[i*4]);
880 /* Add given control point pairs for least squares solving */
881 for (i=0; i < number_arguments; i+=cp_size) {
882 terms[0] = arguments[i+cp_x]; /* x */
883 terms[1] = arguments[i+cp_y]; /* y */
884 terms[2] = terms[0]*terms[1]; /* x*y */
885 terms[3] = 1; /* 1 */
886 LeastSquaresAddTerms(matrix,vectors,terms,
887 &(arguments[i+cp_values]),4UL,number_values);
889 /* Solve for LeastSquares Coefficients */
890 status=GaussJordanElimination(matrix,vectors,4UL,number_values);
891 matrix = RelinquishMagickMatrix(matrix, 4UL);
892 vectors = (double **) RelinquishMagickMemory(vectors);
893 if ( status == MagickFalse ) {
894 coeff = (double *) RelinquishMagickMemory(coeff);
895 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
896 "InvalidArgument","%s : 'Unsolvable Matrix'",
897 MagickOptionToMnemonic(MagickDistortOptions, *method) );
898 return((double *) NULL);
900 if ( *method == BilinearForwardDistortion ) {
901 /* Bilinear Forward Mapped Distortion
903 The above least-squares solved for coefficents but in the forward
904 direction, due to changes to indexing constants.
906 i = c0*x + c1*y + c2*x*y + c3;
907 j = c4*x + c5*y + c6*x*y + c7;
909 where i,j are in the destination image, NOT the source.
911 Reverse Pixel mapping however needs to use reverse of these
912 functions. It required a full page of algbra to work out the
913 reversed mapping formula, but resolves down to the following...
916 c9 = 2*(c2*c5-c1*c6); // '2*a' in the quadratic formula
918 i = i - c3; j = j - c7;
919 b = c6*i - c2*j + c8; // So that a*y^2 + b*y + c == 0
920 c = c4*i - c0*j; // y = ( -b +- sqrt(bb - 4ac) ) / (2*a)
924 y = ( -b + sqrt(r) ) / c9;
928 x = ( i - c1*y) / ( c1 - c2*y );
930 NB: if 'r' is negative there is no solution!
931 NB: the sign of the sqrt() should be negative if image becomes
932 flipped or flopped, or crosses over itself.
933 NB: techniqually coefficient c5 is not needed, anymore,
934 but kept for completness.
936 See Anthony Thyssen <A.Thyssen@griffith.edu.au>
937 or Fred Weinhaus <fmw@alink.net> for more details.
940 coeff[8] = coeff[0]*coeff[5] - coeff[1]*coeff[4];
941 coeff[9] = 2*(coeff[2]*coeff[5] - coeff[1]*coeff[6]);
946 case QuadrilateralDistortion:
948 /* Map a Quadrilateral to a unit square using BilinearReverse
949 Then map that unit square back to the final Quadrilateral
950 using BilinearForward.
952 Input Arguments are sets of control points...
953 For Distort Images u,v, x,y ...
954 For Sparse Gradients x,y, r,g,b ...
957 /* UNDER CONSTRUCTION */
962 case PolynomialDistortion:
964 /* Polynomial Distortion
966 First two coefficents are used to hole global polynomal information
967 c0 = Order of the polynimial being created
968 c1 = number_of_terms in one polynomial equation
970 Rest of the coefficients map to the equations....
971 v = c0 + c1*x + c2*y + c3*x*y + c4*x^2 + c5*y^2 + c6*x^3 + ...
972 for each 'value' (number_values of them) given.
973 As such total coefficients = 2 + number_terms * number_values
975 Input Arguments are sets of control points...
976 For Distort Images order [u,v, x,y] ...
977 For Sparse Gradients order [x,y, r,g,b] ...
979 Polynomial Distortion Notes...
980 + UNDER DEVELOPMENT -- Do not expect this to remain as is.
981 + Currently polynomial is a reversed mapped distortion.
982 + Order 1.5 is fudged to map into a bilinear distortion.
983 though it is not the same order as that distortion.
991 nterms; /* number of polynomial terms per number_values */
999 /* first two coefficients hold polynomial order information */
1000 coeff[0] = arguments[0];
1001 coeff[1] = (double) poly_number_terms(arguments[0]);
1002 nterms = (size_t) coeff[1];
1004 /* create matrix, a fake vectors matrix, and least sqs terms */
1005 matrix = AcquireMagickMatrix(nterms,nterms);
1006 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
1007 terms = (double *) AcquireQuantumMemory(nterms, sizeof(*terms));
1008 if (matrix == (double **) NULL ||
1009 vectors == (double **) NULL ||
1010 terms == (double *) NULL )
1012 matrix = RelinquishMagickMatrix(matrix, nterms);
1013 vectors = (double **) RelinquishMagickMemory(vectors);
1014 terms = (double *) RelinquishMagickMemory(terms);
1015 coeff = (double *) RelinquishMagickMemory(coeff);
1016 (void) ThrowMagickException(exception,GetMagickModule(),
1017 ResourceLimitError,"MemoryAllocationFailed",
1018 "%s", "DistortCoefficients");
1019 return((double *) NULL);
1021 /* fake a number_values x3 vectors matrix from coefficients array */
1022 for (i=0; i < number_values; i++)
1023 vectors[i] = &(coeff[2+i*nterms]);
1024 /* Add given control point pairs for least squares solving */
1025 for (i=1; i < number_arguments; i+=cp_size) { /* NB: start = 1 not 0 */
1026 for (j=0; j < (ssize_t) nterms; j++)
1027 terms[j] = poly_basis_fn(j,arguments[i+cp_x],arguments[i+cp_y]);
1028 LeastSquaresAddTerms(matrix,vectors,terms,
1029 &(arguments[i+cp_values]),nterms,number_values);
1031 terms = (double *) RelinquishMagickMemory(terms);
1032 /* Solve for LeastSquares Coefficients */
1033 status=GaussJordanElimination(matrix,vectors,nterms,number_values);
1034 matrix = RelinquishMagickMatrix(matrix, nterms);
1035 vectors = (double **) RelinquishMagickMemory(vectors);
1036 if ( status == MagickFalse ) {
1037 coeff = (double *) RelinquishMagickMemory(coeff);
1038 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1039 "InvalidArgument","%s : 'Unsolvable Matrix'",
1040 MagickOptionToMnemonic(MagickDistortOptions, *method) );
1041 return((double *) NULL);
1048 Args: arc_width rotate top_edge_radius bottom_edge_radius
1049 All but first argument are optional
1050 arc_width The angle over which to arc the image side-to-side
1051 rotate Angle to rotate image from vertical center
1052 top_radius Set top edge of source image at this radius
1053 bottom_radius Set bootom edge to this radius (radial scaling)
1055 By default, if the radii arguments are nor provided the image radius
1056 is calculated so the horizontal center-line is fits the given arc
1059 The output image size is ALWAYS adjusted to contain the whole image,
1060 and an offset is given to position image relative to the 0,0 point of
1061 the origin, allowing users to use relative positioning onto larger
1062 background (via -flatten).
1064 The arguments are converted to these coefficients
1065 c0: angle for center of source image
1066 c1: angle scale for mapping to source image
1067 c2: radius for top of source image
1068 c3: radius scale for mapping source image
1069 c4: centerline of arc within source image
1071 Note the coefficients use a center angle, so asymptotic join is
1072 furthest from both sides of the source image. This also means that
1073 for arc angles greater than 360 the sides of the image will be
1076 Arc Distortion Notes...
1077 + Does not use a set of CPs
1078 + Will only work with Image Distortion
1079 + Can not be used for generating a sparse gradient (interpolation)
1081 if ( number_arguments >= 1 && arguments[0] < MagickEpsilon ) {
1082 coeff = (double *) RelinquishMagickMemory(coeff);
1083 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1084 "InvalidArgument","%s : 'Arc Angle Too Small'",
1085 MagickOptionToMnemonic(MagickDistortOptions, *method) );
1086 return((double *) NULL);
1088 if ( number_arguments >= 3 && arguments[2] < MagickEpsilon ) {
1089 coeff = (double *) RelinquishMagickMemory(coeff);
1090 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1091 "InvalidArgument","%s : 'Outer Radius Too Small'",
1092 MagickOptionToMnemonic(MagickDistortOptions, *method) );
1093 return((double *) NULL);
1095 coeff[0] = -MagickPI2; /* -90, place at top! */
1096 if ( number_arguments >= 1 )
1097 coeff[1] = DegreesToRadians(arguments[0]);
1099 coeff[1] = MagickPI2; /* zero arguments - center is at top */
1100 if ( number_arguments >= 2 )
1101 coeff[0] += DegreesToRadians(arguments[1]);
1102 coeff[0] /= Magick2PI; /* normalize radians */
1103 coeff[0] -= MagickRound(coeff[0]);
1104 coeff[0] *= Magick2PI; /* de-normalize back to radians */
1105 coeff[3] = (double)image->rows-1;
1106 coeff[2] = (double)image->columns/coeff[1] + coeff[3]/2.0;
1107 if ( number_arguments >= 3 ) {
1108 if ( number_arguments >= 4 )
1109 coeff[3] = arguments[2] - arguments[3];
1111 coeff[3] *= arguments[2]/coeff[2];
1112 coeff[2] = arguments[2];
1114 coeff[4] = ((double)image->columns-1.0)/2.0;
1118 case PolarDistortion:
1119 case DePolarDistortion:
1121 /* (De)Polar Distortion (same set of arguments)
1122 Args: Rmax, Rmin, Xcenter,Ycenter, Afrom,Ato
1123 DePolar can also have the extra arguments of Width, Height
1125 Coefficients 0 to 5 is the sanatized version first 6 input args
1126 Coefficient 6 is the angle to coord ratio and visa-versa
1127 Coefficient 7 is the radius to coord ratio and visa-versa
1129 WARNING: It is posible for Radius max<min and/or Angle from>to
1131 if ( number_arguments == 3
1132 || ( number_arguments > 6 && *method == PolarDistortion )
1133 || number_arguments > 8 ) {
1134 (void) ThrowMagickException(exception,GetMagickModule(),
1135 OptionError,"InvalidArgument", "%s : number of arguments",
1136 MagickOptionToMnemonic(MagickDistortOptions, *method) );
1137 coeff=(double *) RelinquishMagickMemory(coeff);
1138 return((double *) NULL);
1140 /* Rmax - if 0 calculate appropriate value */
1141 if ( number_arguments >= 1 )
1142 coeff[0] = arguments[0];
1145 /* Rmin - usally 0 */
1146 coeff[1] = number_arguments >= 2 ? arguments[1] : 0.0;
1148 if ( number_arguments >= 4 ) {
1149 coeff[2] = arguments[2];
1150 coeff[3] = arguments[3];
1152 else { /* center of actual image */
1153 coeff[2] = (double)(image->columns)/2.0+image->page.x;
1154 coeff[3] = (double)(image->rows)/2.0+image->page.y;
1156 /* Angle from,to - about polar center 0 is downward */
1157 coeff[4] = -MagickPI;
1158 if ( number_arguments >= 5 )
1159 coeff[4] = DegreesToRadians(arguments[4]);
1160 coeff[5] = coeff[4];
1161 if ( number_arguments >= 6 )
1162 coeff[5] = DegreesToRadians(arguments[5]);
1163 if ( fabs(coeff[4]-coeff[5]) < MagickEpsilon )
1164 coeff[5] += Magick2PI; /* same angle is a full circle */
1165 /* if radius 0 or negative, its a special value... */
1166 if ( coeff[0] < MagickEpsilon ) {
1167 /* Use closest edge if radius == 0 */
1168 if ( fabs(coeff[0]) < MagickEpsilon ) {
1169 coeff[0]=MagickMin(fabs(coeff[2]-image->page.x),
1170 fabs(coeff[3]-image->page.y));
1171 coeff[0]=MagickMin(coeff[0],
1172 fabs(coeff[2]-image->page.x-image->columns));
1173 coeff[0]=MagickMin(coeff[0],
1174 fabs(coeff[3]-image->page.y-image->rows));
1176 /* furthest diagonal if radius == -1 */
1177 if ( fabs(-1.0-coeff[0]) < MagickEpsilon ) {
1179 rx = coeff[2]-image->page.x;
1180 ry = coeff[3]-image->page.y;
1181 coeff[0] = rx*rx+ry*ry;
1182 ry = coeff[3]-image->page.y-image->rows;
1183 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1184 rx = coeff[2]-image->page.x-image->columns;
1185 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1186 ry = coeff[3]-image->page.y;
1187 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1188 coeff[0] = sqrt(coeff[0]);
1191 /* IF Rmax <= 0 or Rmin < 0 OR Rmax < Rmin, THEN error */
1192 if ( coeff[0] < MagickEpsilon || coeff[1] < -MagickEpsilon
1193 || (coeff[0]-coeff[1]) < MagickEpsilon ) {
1194 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1195 "InvalidArgument", "%s : Invalid Radius",
1196 MagickOptionToMnemonic(MagickDistortOptions, *method) );
1197 coeff=(double *) RelinquishMagickMemory(coeff);
1198 return((double *) NULL);
1200 /* converstion ratios */
1201 if ( *method == PolarDistortion ) {
1202 coeff[6]=(double) image->columns/(coeff[5]-coeff[4]);
1203 coeff[7]=(double) image->rows/(coeff[0]-coeff[1]);
1205 else { /* *method == DePolarDistortion */
1206 coeff[6]=(coeff[5]-coeff[4])/image->columns;
1207 coeff[7]=(coeff[0]-coeff[1])/image->rows;
1211 case BarrelDistortion:
1212 case BarrelInverseDistortion:
1214 /* Barrel Distortion
1215 Rs=(A*Rd^3 + B*Rd^2 + C*Rd + D)*Rd
1216 BarrelInv Distortion
1217 Rs=Rd/(A*Rd^3 + B*Rd^2 + C*Rd + D)
1219 Where Rd is the normalized radius from corner to middle of image
1220 Input Arguments are one of the following forms (number of arguments)...
1225 8: Ax,Bx,Cx,Dx Ay,By,Cy,Dy
1226 10: Ax,Bx,Cx,Dx Ay,By,Cy,Dy X,Y
1228 Returns 10 coefficent values, which are de-normalized (pixel scale)
1229 Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Xc, Yc
1231 /* Radius de-normalization scaling factor */
1233 rscale = 2.0/MagickMin((double) image->columns,(double) image->rows);
1235 /* sanity check number of args must = 3,4,5,6,8,10 or error */
1236 if ( (number_arguments < 3) || (number_arguments == 7) ||
1237 (number_arguments == 9) || (number_arguments > 10) )
1239 coeff=(double *) RelinquishMagickMemory(coeff);
1240 (void) ThrowMagickException(exception,GetMagickModule(),
1241 OptionError,"InvalidArgument", "%s : number of arguments",
1242 MagickOptionToMnemonic(MagickDistortOptions, *method) );
1243 return((double *) NULL);
1245 /* A,B,C,D coefficients */
1246 coeff[0] = arguments[0];
1247 coeff[1] = arguments[1];
1248 coeff[2] = arguments[2];
1249 if ((number_arguments == 3) || (number_arguments == 5) )
1250 coeff[3] = 1.0 - coeff[0] - coeff[1] - coeff[2];
1252 coeff[3] = arguments[3];
1253 /* de-normalize the coefficients */
1254 coeff[0] *= pow(rscale,3.0);
1255 coeff[1] *= rscale*rscale;
1257 /* Y coefficients: as given OR same as X coefficients */
1258 if ( number_arguments >= 8 ) {
1259 coeff[4] = arguments[4] * pow(rscale,3.0);
1260 coeff[5] = arguments[5] * rscale*rscale;
1261 coeff[6] = arguments[6] * rscale;
1262 coeff[7] = arguments[7];
1265 coeff[4] = coeff[0];
1266 coeff[5] = coeff[1];
1267 coeff[6] = coeff[2];
1268 coeff[7] = coeff[3];
1270 /* X,Y Center of Distortion (image coodinates) */
1271 if ( number_arguments == 5 ) {
1272 coeff[8] = arguments[3];
1273 coeff[9] = arguments[4];
1275 else if ( number_arguments == 6 ) {
1276 coeff[8] = arguments[4];
1277 coeff[9] = arguments[5];
1279 else if ( number_arguments == 10 ) {
1280 coeff[8] = arguments[8];
1281 coeff[9] = arguments[9];
1284 /* center of the image provided (image coodinates) */
1285 coeff[8] = (double)image->columns/2.0 + image->page.x;
1286 coeff[9] = (double)image->rows/2.0 + image->page.y;
1290 case ShepardsDistortion:
1291 case VoronoiColorInterpolate:
1293 /* Shepards Distortion input arguments are the coefficents!
1294 Just check the number of arguments is valid!
1295 Args: u1,v1, x1,y1, ...
1296 OR : u1,v1, r1,g1,c1, ...
1298 if ( number_arguments%cp_size != 0 ||
1299 number_arguments < cp_size ) {
1300 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1301 "InvalidArgument", "%s : 'require at least %.20g CPs'",
1302 MagickOptionToMnemonic(MagickDistortOptions, *method), 1.0);
1303 coeff=(double *) RelinquishMagickMemory(coeff);
1304 return((double *) NULL);
1311 /* you should never reach this point */
1312 assert(! "No Method Handler"); /* just fail assertion */
1313 return((double *) NULL);
1317 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1321 % D i s t o r t I m a g e %
1325 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1327 % DistortImage() distorts an image using various distortion methods, by
1328 % mapping color lookups of the source image to a new destination image
1329 % usally of the same size as the source image, unless 'bestfit' is set to
1332 % If 'bestfit' is enabled, and distortion allows it, the destination image is
1333 % adjusted to ensure the whole source 'image' will just fit within the final
1334 % destination image, which will be sized and offset accordingly. Also in
1335 % many cases the virtual offset of the source image will be taken into
1336 % account in the mapping.
1338 % If the '-verbose' control option has been set print to standard error the
1339 % equicelent '-fx' formula with coefficients for the function, if practical.
1341 % The format of the DistortImage() method is:
1343 % Image *DistortImage(const Image *image,const DistortImageMethod method,
1344 % const size_t number_arguments,const double *arguments,
1345 % MagickBooleanType bestfit, ExceptionInfo *exception)
1347 % A description of each parameter follows:
1349 % o image: the image to be distorted.
1351 % o method: the method of image distortion.
1353 % ArcDistortion always ignores source image offset, and always
1354 % 'bestfit' the destination image with the top left corner offset
1355 % relative to the polar mapping center.
1357 % Affine, Perspective, and Bilinear, do least squares fitting of the
1358 % distrotion when more than the minimum number of control point pairs
1361 % Perspective, and Bilinear, fall back to a Affine distortion when less
1362 % than 4 control point pairs are provided. While Affine distortions
1363 % let you use any number of control point pairs, that is Zero pairs is
1364 % a No-Op (viewport only) distortion, one pair is a translation and
1365 % two pairs of control points do a scale-rotate-translate, without any
1368 % o number_arguments: the number of arguments given.
1370 % o arguments: an array of floating point arguments for this method.
1372 % o bestfit: Attempt to 'bestfit' the size of the resulting image.
1373 % This also forces the resulting image to be a 'layered' virtual
1374 % canvas image. Can be overridden using 'distort:viewport' setting.
1376 % o exception: return any errors or warnings in this structure
1378 % Extra Controls from Image meta-data (artifacts)...
1381 % Output to stderr alternatives, internal coefficents, and FX
1382 % equivelents for the distortion operation (if feasible).
1383 % This forms an extra check of the distortion method, and allows users
1384 % access to the internal constants IM calculates for the distortion.
1386 % o "distort:viewport"
1387 % Directly set the output image canvas area and offest to use for the
1388 % resulting image, rather than use the original images canvas, or a
1389 % calculated 'bestfit' canvas.
1392 % Scale the size of the output canvas by this amount to provide a
1393 % method of Zooming, and for super-sampling the results.
1395 % Other settings that can effect results include
1397 % o 'interpolate' For source image lookups (scale enlargements)
1399 % o 'filter' Set filter to use for area-resampling (scale shrinking).
1400 % Set to 'point' to turn off and use 'interpolate' lookup
1404 MagickExport Image *DistortImage(const Image *image,DistortImageMethod method,
1405 const size_t number_arguments,const double *arguments,
1406 MagickBooleanType bestfit,ExceptionInfo *exception)
1408 #define DistortImageTag "Distort/Image"
1418 geometry; /* geometry of the distorted space viewport */
1423 assert(image != (Image *) NULL);
1424 assert(image->signature == MagickSignature);
1425 if (image->debug != MagickFalse)
1426 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1427 assert(exception != (ExceptionInfo *) NULL);
1428 assert(exception->signature == MagickSignature);
1431 Convert input arguments (usally as control points for reverse mapping)
1432 into mapping coefficients to apply the distortion.
1434 Note that some distortions are mapped to other distortions,
1435 and as such do not require specific code after this point.
1437 coeff = GenerateCoefficients(image, &method, number_arguments,
1438 arguments, 0, exception);
1439 if ( coeff == (double *) NULL )
1440 return((Image *) NULL);
1443 Determine the size and offset for a 'bestfit' destination.
1444 Usally the four corners of the source image is enough.
1447 /* default output image bounds, when no 'bestfit' is requested */
1448 geometry.width=image->columns;
1449 geometry.height=image->rows;
1453 if ( method == ArcDistortion ) {
1454 /* always use the 'best fit' viewport */
1455 bestfit = MagickTrue;
1458 /* Work out the 'best fit', (required for ArcDistortion) */
1463 s.x=s.y=min.x=max.x=min.y=max.y=0.0; /* keep compiler happy */
1465 /* defines to figure out the bounds of the distorted image */
1466 #define InitalBounds(p) \
1468 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1469 min.x = max.x = p.x; \
1470 min.y = max.y = p.y; \
1472 #define ExpandBounds(p) \
1474 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1475 min.x = MagickMin(min.x,p.x); \
1476 max.x = MagickMax(max.x,p.x); \
1477 min.y = MagickMin(min.y,p.y); \
1478 max.y = MagickMax(max.y,p.y); \
1483 case AffineDistortion:
1484 { double inverse[6];
1485 InvertAffineCoefficients(coeff, inverse);
1486 s.x = (double) image->page.x;
1487 s.y = (double) image->page.y;
1488 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1489 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1491 s.x = (double) image->page.x+image->columns;
1492 s.y = (double) image->page.y;
1493 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1494 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1496 s.x = (double) image->page.x;
1497 s.y = (double) image->page.y+image->rows;
1498 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1499 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1501 s.x = (double) image->page.x+image->columns;
1502 s.y = (double) image->page.y+image->rows;
1503 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1504 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1508 case PerspectiveDistortion:
1509 { double inverse[8], scale;
1510 InvertPerspectiveCoefficients(coeff, inverse);
1511 s.x = (double) image->page.x;
1512 s.y = (double) image->page.y;
1513 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1514 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1515 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1516 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1518 s.x = (double) image->page.x+image->columns;
1519 s.y = (double) image->page.y;
1520 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1521 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1522 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1523 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1525 s.x = (double) image->page.x;
1526 s.y = (double) image->page.y+image->rows;
1527 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1528 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1529 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1530 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1532 s.x = (double) image->page.x+image->columns;
1533 s.y = (double) image->page.y+image->rows;
1534 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1535 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1536 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1537 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1543 /* Forward Map Corners */
1544 a = coeff[0]-coeff[1]/2; ca = cos(a); sa = sin(a);
1548 d.x = (coeff[2]-coeff[3])*ca;
1549 d.y = (coeff[2]-coeff[3])*sa;
1551 a = coeff[0]+coeff[1]/2; ca = cos(a); sa = sin(a);
1555 d.x = (coeff[2]-coeff[3])*ca;
1556 d.y = (coeff[2]-coeff[3])*sa;
1558 /* Orthogonal points along top of arc */
1559 for( a=(double) (ceil((double) ((coeff[0]-coeff[1]/2.0)/MagickPI2))*MagickPI2);
1560 a<(coeff[0]+coeff[1]/2.0); a+=MagickPI2 ) {
1561 ca = cos(a); sa = sin(a);
1567 Convert the angle_to_width and radius_to_height
1568 to appropriate scaling factors, to allow faster processing
1569 in the mapping function.
1571 coeff[1] = (double) (Magick2PI*image->columns/coeff[1]);
1572 coeff[3] = (double)image->rows/coeff[3];
1575 case PolarDistortion:
1577 if (number_arguments < 2)
1578 coeff[2] = coeff[3] = 0.0;
1579 min.x = coeff[2]-coeff[0];
1580 max.x = coeff[2]+coeff[0];
1581 min.y = coeff[3]-coeff[0];
1582 max.y = coeff[3]+coeff[0];
1583 /* should be about 1.0 if Rmin = 0 */
1584 coeff[7]=(double) geometry.height/(coeff[0]-coeff[1]);
1587 case DePolarDistortion:
1589 /* direct calculation as it needs to tile correctly
1590 * for reversibility in a DePolar-Polar cycle */
1591 geometry.x = geometry.y = 0;
1592 geometry.height = (size_t) ceil(coeff[0]-coeff[1]);
1593 geometry.width = (size_t)
1594 ceil((coeff[0]-coeff[1])*(coeff[5]-coeff[4])*0.5);
1597 case ShepardsDistortion:
1598 case BilinearForwardDistortion:
1599 case BilinearReverseDistortion:
1601 case QuadrilateralDistortion:
1603 case PolynomialDistortion:
1604 case BarrelDistortion:
1605 case BarrelInverseDistortion:
1607 /* no bestfit available for this distortion */
1608 bestfit = MagickFalse;
1612 /* Set the output image geometry to calculated 'bestfit'.
1613 Yes this tends to 'over do' the file image size, ON PURPOSE!
1614 Do not do this for DePolar which needs to be exact for virtual tiling.
1616 if ( bestfit && method != DePolarDistortion ) {
1617 geometry.x = (ssize_t) floor(min.x-0.5);
1618 geometry.y = (ssize_t) floor(min.y-0.5);
1619 geometry.width=(size_t) ceil(max.x-geometry.x+0.5);
1620 geometry.height=(size_t) ceil(max.y-geometry.y+0.5);
1623 /* Now that we have a new size lets some distortions to it exactly
1624 This is for correct handling of Depolar and its virtual tile handling
1626 if ( method == DePolarDistortion ) {
1627 coeff[6]=(coeff[5]-coeff[4])/geometry.width; /* changed width */
1628 coeff[7]=(coeff[0]-coeff[1])/geometry.height; /* should be about 1.0 */
1632 /* The user provided a 'viewport' expert option which may
1633 overrides some parts of the current output image geometry.
1634 For ArcDistortion, this also overrides its default 'bestfit' setting.
1636 { const char *artifact=GetImageArtifact(image,"distort:viewport");
1637 viewport_given = MagickFalse;
1638 if ( artifact != (const char *) NULL ) {
1639 (void) ParseAbsoluteGeometry(artifact,&geometry);
1640 viewport_given = MagickTrue;
1644 /* Verbose output */
1645 if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
1648 char image_gen[MaxTextExtent];
1651 /* Set destination image size and virtual offset */
1652 if ( bestfit || viewport_given ) {
1653 (void) FormatMagickString(image_gen, MaxTextExtent," -size %.20gx%.20g "
1654 "-page %+.20g%+.20g xc: +insert \\\n",(double) geometry.width,
1655 (double) geometry.height,(double) geometry.x,(double) geometry.y);
1656 lookup="v.p{ xx-v.page.x-.5, yy-v.page.x-.5 }";
1659 image_gen[0] = '\0'; /* no destination to generate */
1660 lookup = "p{ xx-page.x-.5, yy-page.x-.5 }"; /* simplify lookup */
1664 case AffineDistortion:
1668 inverse = (double *) AcquireQuantumMemory(6,sizeof(*inverse));
1669 if (inverse == (double *) NULL) {
1670 coeff = (double *) RelinquishMagickMemory(coeff);
1671 (void) ThrowMagickException(exception,GetMagickModule(),
1672 ResourceLimitError,"MemoryAllocationFailed",
1673 "%s", "DistortImages");
1674 return((Image *) NULL);
1676 InvertAffineCoefficients(coeff, inverse);
1677 CoefficientsToAffineArgs(inverse);
1678 fprintf(stderr, "Affine Projection:\n");
1679 fprintf(stderr, " -distort AffineProjection \\\n '");
1681 fprintf(stderr, "%lf,", inverse[i]);
1682 fprintf(stderr, "%lf'\n", inverse[5]);
1683 inverse = (double *) RelinquishMagickMemory(inverse);
1685 fprintf(stderr, "Affine Distort, FX Equivelent:\n");
1686 fprintf(stderr, "%s", image_gen);
1687 fprintf(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1688 fprintf(stderr, " xx=%+lf*ii %+lf*jj %+lf;\n",
1689 coeff[0], coeff[1], coeff[2]);
1690 fprintf(stderr, " yy=%+lf*ii %+lf*jj %+lf;\n",
1691 coeff[3], coeff[4], coeff[5]);
1692 fprintf(stderr, " %s'\n", lookup);
1697 case PerspectiveDistortion:
1701 inverse = (double *) AcquireQuantumMemory(8,sizeof(*inverse));
1702 if (inverse == (double *) NULL) {
1703 coeff = (double *) RelinquishMagickMemory(coeff);
1704 (void) ThrowMagickException(exception,GetMagickModule(),
1705 ResourceLimitError,"MemoryAllocationFailed",
1706 "%s", "DistortCoefficients");
1707 return((Image *) NULL);
1709 InvertPerspectiveCoefficients(coeff, inverse);
1710 fprintf(stderr, "Perspective Projection:\n");
1711 fprintf(stderr, " -distort PerspectiveProjection \\\n '");
1713 fprintf(stderr, "%lf, ", inverse[i]);
1714 fprintf(stderr, "\n ");
1716 fprintf(stderr, "%lf, ", inverse[i]);
1717 fprintf(stderr, "%lf'\n", inverse[7]);
1718 inverse = (double *) RelinquishMagickMemory(inverse);
1720 fprintf(stderr, "Perspective Distort, FX Equivelent:\n");
1721 fprintf(stderr, "%s", image_gen);
1722 fprintf(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1723 fprintf(stderr, " rr=%+lf*ii %+lf*jj + 1;\n",
1724 coeff[6], coeff[7]);
1725 fprintf(stderr, " xx=(%+lf*ii %+lf*jj %+lf)/rr;\n",
1726 coeff[0], coeff[1], coeff[2]);
1727 fprintf(stderr, " yy=(%+lf*ii %+lf*jj %+lf)/rr;\n",
1728 coeff[3], coeff[4], coeff[5]);
1729 fprintf(stderr, " rr%s0 ? %s : blue'\n",
1730 coeff[8] < 0 ? "<" : ">", lookup);
1734 case BilinearForwardDistortion:
1735 fprintf(stderr, "BilinearForward Mapping Equations:\n");
1736 fprintf(stderr, "%s", image_gen);
1737 fprintf(stderr, " i = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
1738 coeff[0], coeff[1], coeff[2], coeff[3]);
1739 fprintf(stderr, " j = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
1740 coeff[4], coeff[5], coeff[6], coeff[7]);
1743 fprintf(stderr, " c8 = %+lf c9 = 2*a = %+lf;\n",
1744 coeff[8], coeff[9]);
1746 fprintf(stderr, "BilinearForward Distort, FX Equivelent:\n");
1747 fprintf(stderr, "%s", image_gen);
1748 fprintf(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
1749 0.5-coeff[3], 0.5-coeff[7]);
1750 fprintf(stderr, " bb=%lf*ii %+lf*jj %+lf;\n",
1751 coeff[6], -coeff[2], coeff[8]);
1752 /* Handle Special degenerate (non-quadratic) or trapezoidal case */
1753 if ( coeff[9] != 0 ) {
1754 fprintf(stderr, " rt=bb*bb %+lf*(%lf*ii%+lf*jj);\n",
1755 -2*coeff[9], coeff[4], -coeff[0]);
1756 fprintf(stderr, " yy=( -bb + sqrt(rt) ) / %lf;\n",
1759 fprintf(stderr, " yy=(%lf*ii%+lf*jj)/bb;\n",
1760 -coeff[4], coeff[0]);
1761 fprintf(stderr, " xx=(ii %+lf*yy)/(%lf %+lf*yy);\n",
1762 -coeff[1], coeff[0], coeff[2]);
1763 if ( coeff[9] != 0 )
1764 fprintf(stderr, " (rt < 0 ) ? red : %s'\n", lookup);
1766 fprintf(stderr, " %s'\n", lookup);
1769 case BilinearReverseDistortion:
1771 fprintf(stderr, "Polynomial Projection Distort:\n");
1772 fprintf(stderr, " -distort PolynomialProjection \\\n");
1773 fprintf(stderr, " '1.5, %lf, %lf, %lf, %lf,\n",
1774 coeff[3], coeff[0], coeff[1], coeff[2]);
1775 fprintf(stderr, " %lf, %lf, %lf, %lf'\n",
1776 coeff[7], coeff[4], coeff[5], coeff[6]);
1778 fprintf(stderr, "BilinearReverse Distort, FX Equivelent:\n");
1779 fprintf(stderr, "%s", image_gen);
1780 fprintf(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1781 fprintf(stderr, " xx=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
1782 coeff[0], coeff[1], coeff[2], coeff[3]);
1783 fprintf(stderr, " yy=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
1784 coeff[4], coeff[5], coeff[6], coeff[7]);
1785 fprintf(stderr, " %s'\n", lookup);
1788 case PolynomialDistortion:
1790 size_t nterms = (size_t) coeff[1];
1791 fprintf(stderr, "Polynomial (order %lg, terms %lu), FX Equivelent\n",
1792 coeff[0],(unsigned long) nterms);
1793 fprintf(stderr, "%s", image_gen);
1794 fprintf(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1795 fprintf(stderr, " xx =");
1796 for (i=0; i<(ssize_t) nterms; i++) {
1797 if ( i != 0 && i%4 == 0 ) fprintf(stderr, "\n ");
1798 fprintf(stderr, " %+lf%s", coeff[2+i],
1801 fprintf(stderr, ";\n yy =");
1802 for (i=0; i<(ssize_t) nterms; i++) {
1803 if ( i != 0 && i%4 == 0 ) fprintf(stderr, "\n ");
1804 fprintf(stderr, " %+lf%s", coeff[2+i+nterms],
1807 fprintf(stderr, ";\n %s'\n", lookup);
1812 fprintf(stderr, "Arc Distort, Internal Coefficients:\n");
1813 for ( i=0; i<5; i++ )
1814 fprintf(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
1815 fprintf(stderr, "Arc Distort, FX Equivelent:\n");
1816 fprintf(stderr, "%s", image_gen);
1817 fprintf(stderr, " -fx 'ii=i+page.x; jj=j+page.y;\n");
1818 fprintf(stderr, " xx=(atan2(jj,ii)%+lf)/(2*pi);\n",
1820 fprintf(stderr, " xx=xx-round(xx);\n");
1821 fprintf(stderr, " xx=xx*%lf %+lf;\n",
1822 coeff[1], coeff[4]);
1823 fprintf(stderr, " yy=(%lf - hypot(ii,jj)) * %lf;\n",
1824 coeff[2], coeff[3]);
1825 fprintf(stderr, " v.p{xx-.5,yy-.5}'\n");
1828 case PolarDistortion:
1830 fprintf(stderr, "Polar Distort, Internal Coefficents\n");
1831 for ( i=0; i<8; i++ )
1832 fprintf(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
1833 fprintf(stderr, "Polar Distort, FX Equivelent:\n");
1834 fprintf(stderr, "%s", image_gen);
1835 fprintf(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
1836 -coeff[2], -coeff[3]);
1837 fprintf(stderr, " xx=(atan2(ii,jj)%+lf)/(2*pi);\n",
1838 -(coeff[4]+coeff[5])/2 );
1839 fprintf(stderr, " xx=xx-round(xx);\n");
1840 fprintf(stderr, " xx=xx*2*pi*%lf + v.w/2;\n",
1842 fprintf(stderr, " yy=(hypot(ii,jj)%+lf)*%lf;\n",
1843 -coeff[1], coeff[7] );
1844 fprintf(stderr, " v.p{xx-.5,yy-.5}'\n");
1847 case DePolarDistortion:
1849 fprintf(stderr, "DePolar Distort, Internal Coefficents\n");
1850 for ( i=0; i<8; i++ )
1851 fprintf(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
1852 fprintf(stderr, "DePolar Distort, FX Equivelent:\n");
1853 fprintf(stderr, "%s", image_gen);
1854 fprintf(stderr, " -fx 'aa=(i+.5)*%lf %+lf;\n", coeff[6], -coeff[4] );
1855 fprintf(stderr, " rr=(j+.5)*%lf %+lf;\n", coeff[7], +coeff[1] );
1856 fprintf(stderr, " xx=rr*sin(aa) %+lf;\n", coeff[2] );
1857 fprintf(stderr, " yy=rr*cos(aa) %+lf;\n", coeff[3] );
1858 fprintf(stderr, " v.p{xx-.5,yy-.5}'\n");
1861 case BarrelDistortion:
1862 case BarrelInverseDistortion:
1864 /* NOTE: This does the barrel roll in pixel coords not image coords
1865 ** The internal distortion must do it in image coordinates,
1866 ** so that is what the center coeff (8,9) is given in.
1868 xc = ((double)image->columns-1.0)/2.0 + image->page.x;
1869 yc = ((double)image->rows-1.0)/2.0 + image->page.y;
1870 fprintf(stderr, "Barrel%s Distort, FX Equivelent:\n",
1871 method == BarrelDistortion ? "" : "Inv");
1872 fprintf(stderr, "%s", image_gen);
1873 if ( fabs(coeff[8]-xc-0.5) < 0.1 && fabs(coeff[9]-yc-0.5) < 0.1 )
1874 fprintf(stderr, " -fx 'xc=(w-1)/2; yc=(h-1)/2;\n");
1876 fprintf(stderr, " -fx 'xc=%lf; yc=%lf;\n",
1877 coeff[8]-0.5, coeff[9]-0.5);
1879 " ii=i-xc; jj=j-yc; rr=hypot(ii,jj);\n");
1880 fprintf(stderr, " ii=ii%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
1881 method == BarrelDistortion ? "*" : "/",
1882 coeff[0],coeff[1],coeff[2],coeff[3]);
1883 fprintf(stderr, " jj=jj%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
1884 method == BarrelDistortion ? "*" : "/",
1885 coeff[4],coeff[5],coeff[6],coeff[7]);
1886 fprintf(stderr, " v.p{fx*ii+xc,fy*jj+yc}'\n");
1893 /* The user provided a 'scale' expert option will scale the
1894 output image size, by the factor given allowing for super-sampling
1895 of the distorted image space. Any scaling factors must naturally
1896 be halved as a result.
1898 { const char *artifact;
1899 artifact=GetImageArtifact(image,"distort:scale");
1900 output_scaling = 1.0;
1901 if (artifact != (const char *) NULL) {
1902 output_scaling = fabs(StringToDouble(artifact));
1903 geometry.width *= output_scaling;
1904 geometry.height *= output_scaling;
1905 geometry.x *= output_scaling;
1906 geometry.y *= output_scaling;
1907 if ( output_scaling < 0.1 ) {
1908 coeff = (double *) RelinquishMagickMemory(coeff);
1909 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1910 "InvalidArgument","%s", "-set option:distort:scale" );
1911 return((Image *) NULL);
1913 output_scaling = 1/output_scaling;
1916 #define ScaleFilter(F,A,B,C,D) \
1917 ScaleResampleFilter( (F), \
1918 output_scaling*(A), output_scaling*(B), \
1919 output_scaling*(C), output_scaling*(D) )
1922 Initialize the distort image attributes.
1924 distort_image=CloneImage(image,geometry.width,geometry.height,MagickTrue,
1926 if (distort_image == (Image *) NULL)
1927 return((Image *) NULL);
1928 /* if image is ColorMapped - change it to DirectClass */
1929 if (SetImageStorageClass(distort_image,DirectClass) == MagickFalse)
1931 InheritException(exception,&distort_image->exception);
1932 distort_image=DestroyImage(distort_image);
1933 return((Image *) NULL);
1935 distort_image->page.x=geometry.x;
1936 distort_image->page.y=geometry.y;
1937 if (distort_image->background_color.opacity != OpaqueOpacity)
1938 distort_image->matte=MagickTrue;
1940 { /* ----- MAIN CODE -----
1941 Sample the source image to each pixel in the distort image.
1956 **restrict resample_filter;
1963 GetMagickPixelPacket(distort_image,&zero);
1964 resample_filter=AcquireResampleFilterThreadSet(image,
1965 UndefinedVirtualPixelMethod,MagickFalse,exception);
1966 distort_view=AcquireCacheView(distort_image);
1967 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1968 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
1970 for (j=0; j < (ssize_t) distort_image->rows; j++)
1973 id = GetOpenMPThreadId();
1976 validity; /* how mathematically valid is this the mapping */
1982 pixel, /* pixel color to assign to distorted image */
1983 invalid; /* the color to assign when distort result is invalid */
1987 s; /* transform destination image x,y to source image x,y */
1989 register IndexPacket
1995 register PixelPacket
1998 q=QueueCacheViewAuthenticPixels(distort_view,0,j,distort_image->columns,1,
2000 if (q == (PixelPacket *) NULL)
2005 indexes=GetCacheViewAuthenticIndexQueue(distort_view);
2008 /* Define constant scaling vectors for Affine Distortions
2009 Other methods are either variable, or use interpolated lookup
2013 case AffineDistortion:
2014 ScaleFilter( resample_filter[id],
2016 coeff[3], coeff[4] );
2022 /* Initialize default pixel validity
2023 * negative: pixel is invalid output 'matte_color'
2024 * 0.0 to 1.0: antialiased, mix with resample output
2025 * 1.0 or greater: use resampled output.
2029 GetMagickPixelPacket(distort_image,&invalid);
2030 SetMagickPixelPacket(distort_image,&distort_image->matte_color,
2031 (IndexPacket *) NULL, &invalid);
2032 if (distort_image->colorspace == CMYKColorspace)
2033 ConvertRGBToCMYK(&invalid); /* what about other color spaces? */
2035 for (i=0; i < (ssize_t) distort_image->columns; i++)
2037 /* map pixel coordinate to distortion space coordinate */
2038 d.x = (double) (geometry.x+i+0.5)*output_scaling;
2039 d.y = (double) (geometry.y+j+0.5)*output_scaling;
2040 s = d; /* default is a no-op mapping */
2043 case AffineDistortion:
2045 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2046 s.y=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2047 /* Affine partial derivitives are constant -- set above */
2050 case PerspectiveDistortion:
2053 p,q,r,abs_r,abs_c6,abs_c7,scale;
2054 /* perspective is a ratio of affines */
2055 p=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2056 q=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2057 r=coeff[6]*d.x+coeff[7]*d.y+1.0;
2058 /* Pixel Validity -- is it a 'sky' or 'ground' pixel */
2059 validity = (r*coeff[8] < 0.0) ? 0.0 : 1.0;
2060 /* Determine horizon anti-alias blending */
2062 abs_c6 = fabs(coeff[6]);
2063 abs_c7 = fabs(coeff[7]);
2064 if ( abs_c6 > abs_c7 ) {
2065 if ( abs_r < abs_c6*output_scaling )
2066 validity = 0.5 - coeff[8]*r/(coeff[6]*output_scaling);
2068 else if ( abs_r < abs_c7*output_scaling )
2069 validity = 0.5 - coeff[8]*r/(coeff[7]*output_scaling);
2070 /* Perspective Sampling Point (if valid) */
2071 if ( validity > 0.0 ) {
2072 /* divide by r affine, for perspective scaling */
2076 /* Perspective Partial Derivatives or Scaling Vectors */
2078 ScaleFilter( resample_filter[id],
2079 (r*coeff[0] - p*coeff[6])*scale,
2080 (r*coeff[1] - p*coeff[7])*scale,
2081 (r*coeff[3] - q*coeff[6])*scale,
2082 (r*coeff[4] - q*coeff[7])*scale );
2086 case BilinearReverseDistortion:
2088 /* Reversed Mapped is just a simple polynomial */
2089 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2]*d.x*d.y+coeff[3];
2090 s.y=coeff[4]*d.x+coeff[5]*d.y
2091 +coeff[6]*d.x*d.y+coeff[7];
2092 /* Bilinear partial derivitives of scaling vectors */
2093 ScaleFilter( resample_filter[id],
2094 coeff[0] + coeff[2]*d.y,
2095 coeff[1] + coeff[2]*d.x,
2096 coeff[4] + coeff[6]*d.y,
2097 coeff[5] + coeff[6]*d.x );
2100 case BilinearForwardDistortion:
2102 /* Forward mapped needs reversed polynomial equations
2103 * which unfortunatally requires a square root! */
2105 d.x -= coeff[3]; d.y -= coeff[7];
2106 b = coeff[6]*d.x - coeff[2]*d.y + coeff[8];
2107 c = coeff[4]*d.x - coeff[0]*d.y;
2110 /* Handle Special degenerate (non-quadratic) case */
2111 if ( fabs(coeff[9]) < MagickEpsilon )
2114 c = b*b - 2*coeff[9]*c;
2118 s.y = ( -b + sqrt(c) )/coeff[9];
2120 if ( validity > 0.0 )
2121 s.x = ( d.x - coeff[1]*s.y) / ( coeff[0] + coeff[2]*s.y );
2123 /* NOTE: the sign of the square root should be -ve for parts
2124 where the source image becomes 'flipped' or 'mirrored'.
2125 FUTURE: Horizon handling
2126 FUTURE: Scaling factors or Deritives (how?)
2131 case BilinearDistortion:
2132 /* Bilinear mapping of any Quadrilateral to any Quadrilateral */
2133 /* UNDER DEVELOPMENT */
2136 case PolynomialDistortion:
2138 /* multi-ordered polynomial */
2142 nterms=(ssize_t)coeff[1];
2145 du,dv; /* the du,dv vectors from unit dx,dy -- derivatives */
2147 s.x=s.y=du.x=du.y=dv.x=dv.y=0.0;
2148 for(k=0; k < nterms; k++) {
2149 s.x += poly_basis_fn(k,d.x,d.y)*coeff[2+k];
2150 du.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k];
2151 du.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k];
2152 s.y += poly_basis_fn(k,d.x,d.y)*coeff[2+k+nterms];
2153 dv.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k+nterms];
2154 dv.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k+nterms];
2156 ScaleFilter( resample_filter[id], du.x,du.y,dv.x,dv.y );
2161 /* what is the angle and radius in the destination image */
2162 s.x = (double) ((atan2(d.y,d.x) - coeff[0])/Magick2PI);
2163 s.x -= MagickRound(s.x); /* angle */
2164 s.y = hypot(d.x,d.y); /* radius */
2166 /* Arc Distortion Partial Scaling Vectors
2167 Are derived by mapping the perpendicular unit vectors
2168 dR and dA*R*2PI rather than trying to map dx and dy
2169 The results is a very simple orthogonal aligned ellipse.
2171 if ( s.y > MagickEpsilon )
2172 ScaleFilter( resample_filter[id],
2173 (double) (coeff[1]/(Magick2PI*s.y)), 0, 0, coeff[3] );
2175 ScaleFilter( resample_filter[id],
2176 distort_image->columns*2, 0, 0, coeff[3] );
2178 /* now scale the angle and radius for source image lookup point */
2179 s.x = s.x*coeff[1] + coeff[4] + image->page.x +0.5;
2180 s.y = (coeff[2] - s.y) * coeff[3] + image->page.y;
2183 case PolarDistortion:
2184 { /* Rect/Cartesain/Cylinder to Polar View */
2187 s.x = atan2(d.x,d.y) - (coeff[4]+coeff[5])/2;
2189 s.x -= MagickRound(s.x);
2190 s.x *= Magick2PI; /* angle - relative to centerline */
2191 s.y = hypot(d.x,d.y); /* radius */
2193 /* Polar Scaling vectors are based on mapping dR and dA vectors
2194 This results in very simple orthogonal scaling vectors
2196 if ( s.y > MagickEpsilon )
2197 ScaleFilter( resample_filter[id],
2198 (double) (coeff[6]/(Magick2PI*s.y)), 0, 0, coeff[7] );
2200 ScaleFilter( resample_filter[id],
2201 distort_image->columns*2, 0, 0, coeff[7] );
2203 /* now finish mapping radius/angle to source x,y coords */
2204 s.x = s.x*coeff[6] + (double)image->columns/2.0 + image->page.x;
2205 s.y = (s.y-coeff[1])*coeff[7] + image->page.y;
2208 case DePolarDistortion:
2209 { /* Polar to Cylindical */
2210 /* ignore all destination virtual offsets */
2211 d.x = ((double)i+0.5)*output_scaling*coeff[6]-coeff[4];
2212 d.y = ((double)j+0.5)*output_scaling*coeff[7]+coeff[1];
2213 s.x = d.y*sin(d.x) + coeff[2];
2214 s.y = d.y*cos(d.x) + coeff[3];
2215 /* derivatives are usless - better to use SuperSampling */
2218 case BarrelDistortion:
2219 case BarrelInverseDistortion:
2221 double r,fx,fy,gx,gy;
2222 /* Radial Polynomial Distortion (de-normalized) */
2225 r = sqrt(d.x*d.x+d.y*d.y);
2226 if ( r > MagickEpsilon ) {
2227 fx = ((coeff[0]*r + coeff[1])*r + coeff[2])*r + coeff[3];
2228 fy = ((coeff[4]*r + coeff[5])*r + coeff[6])*r + coeff[7];
2229 gx = ((3*coeff[0]*r + 2*coeff[1])*r + coeff[2])/r;
2230 gy = ((3*coeff[4]*r + 2*coeff[5])*r + coeff[6])/r;
2231 /* adjust functions and scaling for 'inverse' form */
2232 if ( method == BarrelInverseDistortion ) {
2233 fx = 1/fx; fy = 1/fy;
2234 gx *= -fx*fx; gy *= -fy*fy;
2236 /* Set the source pixel to lookup and EWA derivative vectors */
2237 s.x = d.x*fx + coeff[8];
2238 s.y = d.y*fy + coeff[9];
2239 ScaleFilter( resample_filter[id],
2240 gx*d.x*d.x + fx, gx*d.x*d.y,
2241 gy*d.x*d.y, gy*d.y*d.y + fy );
2244 /* Special handling to avoid divide by zero when r==0
2246 ** The source and destination pixels match in this case
2247 ** which was set at the top of the loop using s = d;
2248 ** otherwise... s.x=coeff[8]; s.y=coeff[9];
2250 if ( method == BarrelDistortion )
2251 ScaleFilter( resample_filter[id],
2252 coeff[3], 0, 0, coeff[7] );
2253 else /* method == BarrelInverseDistortion */
2254 /* FUTURE, trap for D==0 causing division by zero */
2255 ScaleFilter( resample_filter[id],
2256 1.0/coeff[3], 0, 0, 1.0/coeff[7] );
2260 case ShepardsDistortion:
2261 { /* Shepards Method, or Inverse Weighted Distance for
2262 displacement around the destination image control points
2263 The input arguments are the coefficents to the function.
2264 This is more of a 'displacement' function rather than an
2265 absolute distortion function.
2272 denominator = s.x = s.y = 0;
2273 for(i=0; i<number_arguments; i+=4) {
2275 ((double)d.x-arguments[i+2])*((double)d.x-arguments[i+2])
2276 + ((double)d.y-arguments[i+3])*((double)d.y-arguments[i+3]);
2282 s.x += (arguments[ i ]-arguments[i+2])*weight;
2283 s.y += (arguments[i+1]-arguments[i+3])*weight;
2284 denominator += weight;
2291 /* We can not determine derivatives using shepards method
2292 only color interpolatation, not area-resampling */
2296 break; /* use the default no-op given above */
2298 /* map virtual canvas location back to real image coordinate */
2299 if ( bestfit && method != ArcDistortion ) {
2300 s.x -= image->page.x;
2301 s.y -= image->page.y;
2306 if ( validity <= 0.0 ) {
2307 /* result of distortion is an invalid pixel - don't resample */
2308 SetPixelPacket(distort_image,&invalid,q,indexes);
2311 /* resample the source image to find its correct color */
2312 (void) ResamplePixelColor(resample_filter[id],s.x,s.y,&pixel);
2313 /* if validity between 0.0 and 1.0 mix result with invalid pixel */
2314 if ( validity < 1.0 ) {
2315 /* Do a blend of sample color and invalid pixel */
2316 /* should this be a 'Blend', or an 'Over' compose */
2317 MagickPixelCompositeBlend(&pixel,validity,&invalid,(1.0-validity),
2320 SetPixelPacket(distort_image,&pixel,q,indexes);
2325 sync=SyncCacheViewAuthenticPixels(distort_view,exception);
2326 if (sync == MagickFalse)
2328 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2333 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2334 #pragma omp critical (MagickCore_DistortImage)
2336 proceed=SetImageProgress(image,DistortImageTag,progress++,
2338 if (proceed == MagickFalse)
2342 fprintf(stderr, "\n");
2345 distort_view=DestroyCacheView(distort_view);
2346 resample_filter=DestroyResampleFilterThreadSet(resample_filter);
2348 if (status == MagickFalse)
2349 distort_image=DestroyImage(distort_image);
2352 /* Arc does not return an offset unless 'bestfit' is in effect
2353 And the user has not provided an overriding 'viewport'.
2355 if ( method == ArcDistortion && !bestfit && !viewport_given ) {
2356 distort_image->page.x = 0;
2357 distort_image->page.y = 0;
2359 coeff = (double *) RelinquishMagickMemory(coeff);
2360 return(distort_image);
2364 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2368 % S p a r s e C o l o r I m a g e %
2372 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2374 % SparseColorImage(), given a set of coordinates, interpolates the colors
2375 % found at those coordinates, across the whole image, using various methods.
2377 % The format of the SparseColorImage() method is:
2379 % Image *SparseColorImage(const Image *image,const ChannelType channel,
2380 % const SparseColorMethod method,const size_t number_arguments,
2381 % const double *arguments,ExceptionInfo *exception)
2383 % A description of each parameter follows:
2385 % o image: the image to be filled in.
2387 % o channel: Specify which color values (in RGBKA sequence) are being set.
2388 % This also determines the number of color_values in above.
2390 % o method: the method to fill in the gradient between the control points.
2392 % The methods used for SparseColor() are often simular to methods
2393 % used for DistortImage(), and even share the same code for determination
2394 % of the function coefficents, though with more dimensions (or resulting
2397 % o number_arguments: the number of arguments given.
2399 % o arguments: array of floating point arguments for this method--
2400 % x,y,color_values-- with color_values given as normalized values.
2402 % o exception: return any errors or warnings in this structure
2405 MagickExport Image *SparseColorImage(const Image *image,
2406 const ChannelType channel,const SparseColorMethod method,
2407 const size_t number_arguments,const double *arguments,
2408 ExceptionInfo *exception)
2410 #define SparseColorTag "Distort/SparseColor"
2427 assert(image != (Image *) NULL);
2428 assert(image->signature == MagickSignature);
2429 if (image->debug != MagickFalse)
2430 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2431 assert(exception != (ExceptionInfo *) NULL);
2432 assert(exception->signature == MagickSignature);
2434 /* Determine number of color values needed per control point */
2436 if ( channel & RedChannel ) number_colors++;
2437 if ( channel & GreenChannel ) number_colors++;
2438 if ( channel & BlueChannel ) number_colors++;
2439 if ( channel & IndexChannel ) number_colors++;
2440 if ( channel & OpacityChannel ) number_colors++;
2443 Convert input arguments into mapping coefficients to apply the distortion.
2444 Note some Methods may fall back to other simpler methods.
2446 distort_method=(DistortImageMethod) method;
2447 coeff = GenerateCoefficients(image, &distort_method, number_arguments,
2448 arguments, number_colors, exception);
2449 if ( coeff == (double *) NULL )
2450 return((Image *) NULL);
2452 /* Verbose output */
2453 if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
2456 case BarycentricColorInterpolate:
2458 register ssize_t x=0;
2459 fprintf(stderr, "Barycentric Sparse Color:\n");
2460 if ( channel & RedChannel )
2461 fprintf(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf' \\\n",
2462 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2463 if ( channel & GreenChannel )
2464 fprintf(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf' \\\n",
2465 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2466 if ( channel & BlueChannel )
2467 fprintf(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf' \\\n",
2468 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2469 if ( channel & IndexChannel )
2470 fprintf(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf' \\\n",
2471 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2472 if ( channel & OpacityChannel )
2473 fprintf(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf' \\\n",
2474 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2477 case BilinearColorInterpolate:
2479 register ssize_t x=0;
2480 fprintf(stderr, "Bilinear Sparse Color\n");
2481 if ( channel & RedChannel )
2482 fprintf(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2483 coeff[ x ], coeff[x+1],
2484 coeff[x+2], coeff[x+3]),x+=4;
2485 if ( channel & GreenChannel )
2486 fprintf(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2487 coeff[ x ], coeff[x+1],
2488 coeff[x+2], coeff[x+3]),x+=4;
2489 if ( channel & BlueChannel )
2490 fprintf(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2491 coeff[ x ], coeff[x+1],
2492 coeff[x+2], coeff[x+3]),x+=4;
2493 if ( channel & IndexChannel )
2494 fprintf(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2495 coeff[ x ], coeff[x+1],
2496 coeff[x+2], coeff[x+3]),x+=4;
2497 if ( channel & OpacityChannel )
2498 fprintf(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2499 coeff[ x ], coeff[x+1],
2500 coeff[x+2], coeff[x+3]),x+=4;
2504 /* unknown, or which are too complex for FX alturnatives */
2509 /* Generate new image for generated interpolated gradient.
2510 * ASIDE: Actually we could have just replaced the colors of the original
2511 * image, but IM core policy, is if storage class could change then clone
2515 sparse_image=CloneImage(image,image->columns,image->rows,MagickTrue,
2517 if (sparse_image == (Image *) NULL)
2518 return((Image *) NULL);
2519 if (SetImageStorageClass(sparse_image,DirectClass) == MagickFalse)
2520 { /* if image is ColorMapped - change it to DirectClass */
2521 InheritException(exception,&image->exception);
2522 sparse_image=DestroyImage(sparse_image);
2523 return((Image *) NULL);
2525 { /* ----- MAIN CODE ----- */
2540 GetMagickPixelPacket(sparse_image,&zero);
2541 sparse_view=AcquireCacheView(sparse_image);
2542 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2543 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2545 for (j=0; j < (ssize_t) sparse_image->rows; j++)
2551 pixel; /* pixel to assign to distorted image */
2553 register IndexPacket
2559 register PixelPacket
2562 q=QueueCacheViewAuthenticPixels(sparse_view,0,j,sparse_image->columns,
2564 if (q == (PixelPacket *) NULL)
2569 /* FUTURE: get pixel from source image - so channel can replace parts */
2570 indexes=GetCacheViewAuthenticIndexQueue(sparse_view);
2572 for (i=0; i < (ssize_t) sparse_image->columns; i++)
2576 case BarycentricColorInterpolate:
2578 register ssize_t x=0;
2579 if ( channel & RedChannel )
2580 pixel.red = coeff[x]*i +coeff[x+1]*j
2582 if ( channel & GreenChannel )
2583 pixel.green = coeff[x]*i +coeff[x+1]*j
2585 if ( channel & BlueChannel )
2586 pixel.blue = coeff[x]*i +coeff[x+1]*j
2588 if ( channel & IndexChannel )
2589 pixel.index = coeff[x]*i +coeff[x+1]*j
2591 if ( channel & OpacityChannel )
2592 pixel.opacity = coeff[x]*i +coeff[x+1]*j
2596 case BilinearColorInterpolate:
2598 register ssize_t x=0;
2599 if ( channel & RedChannel )
2600 pixel.red = coeff[x]*i + coeff[x+1]*j +
2601 coeff[x+2]*i*j + coeff[x+3], x+=4;
2602 if ( channel & GreenChannel )
2603 pixel.green = coeff[x]*i + coeff[x+1]*j +
2604 coeff[x+2]*i*j + coeff[x+3], x+=4;
2605 if ( channel & BlueChannel )
2606 pixel.blue = coeff[x]*i + coeff[x+1]*j +
2607 coeff[x+2]*i*j + coeff[x+3], x+=4;
2608 if ( channel & IndexChannel )
2609 pixel.index = coeff[x]*i + coeff[x+1]*j +
2610 coeff[x+2]*i*j + coeff[x+3], x+=4;
2611 if ( channel & OpacityChannel )
2612 pixel.opacity = coeff[x]*i + coeff[x+1]*j +
2613 coeff[x+2]*i*j + coeff[x+3], x+=4;
2616 case ShepardsColorInterpolate:
2617 { /* Shepards Method, uses its own input arguments as coefficients.
2624 if ( channel & RedChannel ) pixel.red = 0.0;
2625 if ( channel & GreenChannel ) pixel.green = 0.0;
2626 if ( channel & BlueChannel ) pixel.blue = 0.0;
2627 if ( channel & IndexChannel ) pixel.index = 0.0;
2628 if ( channel & OpacityChannel ) pixel.opacity = 0.0;
2630 for(k=0; k<number_arguments; k+=2+number_colors) {
2631 register ssize_t x=(ssize_t) k+2;
2633 ((double)i-arguments[ k ])*((double)i-arguments[ k ])
2634 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
2639 if ( channel & RedChannel )
2640 pixel.red += arguments[x++]*weight;
2641 if ( channel & GreenChannel )
2642 pixel.green += arguments[x++]*weight;
2643 if ( channel & BlueChannel )
2644 pixel.blue += arguments[x++]*weight;
2645 if ( channel & IndexChannel )
2646 pixel.index += arguments[x++]*weight;
2647 if ( channel & OpacityChannel )
2648 pixel.opacity += arguments[x++]*weight;
2649 denominator += weight;
2651 if ( channel & RedChannel ) pixel.red /= denominator;
2652 if ( channel & GreenChannel ) pixel.green /= denominator;
2653 if ( channel & BlueChannel ) pixel.blue /= denominator;
2654 if ( channel & IndexChannel ) pixel.index /= denominator;
2655 if ( channel & OpacityChannel ) pixel.opacity /= denominator;
2658 case VoronoiColorInterpolate:
2660 { /* Just use the closest control point you can find! */
2664 minimum = MagickHuge;
2666 for(k=0; k<number_arguments; k+=2+number_colors) {
2668 ((double)i-arguments[ k ])*((double)i-arguments[ k ])
2669 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
2670 if ( distance < minimum ) {
2671 register ssize_t x=(ssize_t) k+2;
2672 if ( channel & RedChannel ) pixel.red = arguments[x++];
2673 if ( channel & GreenChannel ) pixel.green = arguments[x++];
2674 if ( channel & BlueChannel ) pixel.blue = arguments[x++];
2675 if ( channel & IndexChannel ) pixel.index = arguments[x++];
2676 if ( channel & OpacityChannel ) pixel.opacity = arguments[x++];
2683 /* set the color directly back into the source image */
2684 if ( channel & RedChannel ) pixel.red *= QuantumRange;
2685 if ( channel & GreenChannel ) pixel.green *= QuantumRange;
2686 if ( channel & BlueChannel ) pixel.blue *= QuantumRange;
2687 if ( channel & IndexChannel ) pixel.index *= QuantumRange;
2688 if ( channel & OpacityChannel ) pixel.opacity *= QuantumRange;
2689 SetPixelPacket(sparse_image,&pixel,q,indexes);
2693 sync=SyncCacheViewAuthenticPixels(sparse_view,exception);
2694 if (sync == MagickFalse)
2696 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2701 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2702 #pragma omp critical (MagickCore_SparseColorImage)
2704 proceed=SetImageProgress(image,SparseColorTag,progress++,image->rows);
2705 if (proceed == MagickFalse)
2709 sparse_view=DestroyCacheView(sparse_view);
2710 if (status == MagickFalse)
2711 sparse_image=DestroyImage(sparse_image);
2713 coeff = (double *) RelinquishMagickMemory(coeff);
2714 return(sparse_image);