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 "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"
70 #include "magick/transform.h"
73 Numerous internal routines for image distortions.
75 static inline double MagickMin(const double x,const double y)
77 return( x < y ? x : y);
79 static inline double MagickMax(const double x,const double y)
81 return( x > y ? x : y);
84 static inline void AffineArgsToCoefficients(double *affine)
86 /* map external sx,ry,rx,sy,tx,ty to internal c0,c2,c4,c1,c3,c5 */
87 double tmp[4]; /* note indexes 0 and 5 remain unchanged */
88 tmp[0]=affine[1]; tmp[1]=affine[2]; tmp[2]=affine[3]; tmp[3]=affine[4];
89 affine[3]=tmp[0]; affine[1]=tmp[1]; affine[4]=tmp[2]; affine[2]=tmp[3];
92 static inline void CoefficientsToAffineArgs(double *coeff)
94 /* map internal c0,c1,c2,c3,c4,c5 to external sx,ry,rx,sy,tx,ty */
95 double tmp[4]; /* note indexes 0 and 5 remain unchanged */
96 tmp[0]=coeff[3]; tmp[1]=coeff[1]; tmp[2]=coeff[4]; tmp[3]=coeff[2];
97 coeff[1]=tmp[0]; coeff[2]=tmp[1]; coeff[3]=tmp[2]; coeff[4]=tmp[3];
99 static void InvertAffineCoefficients(const double *coeff,double *inverse)
101 /* From "Digital Image Warping" by George Wolberg, page 50 */
104 determinant=1.0/(coeff[0]*coeff[4]-coeff[1]*coeff[3]);
105 inverse[0]=determinant*coeff[4];
106 inverse[1]=determinant*(-coeff[1]);
107 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[2]*coeff[4]);
108 inverse[3]=determinant*(-coeff[3]);
109 inverse[4]=determinant*coeff[0];
110 inverse[5]=determinant*(coeff[2]*coeff[3]-coeff[0]*coeff[5]);
113 static void InvertPerspectiveCoefficients(const double *coeff,
116 /* From "Digital Image Warping" by George Wolberg, page 53 */
119 determinant=1.0/(coeff[0]*coeff[4]-coeff[3]*coeff[1]);
120 inverse[0]=determinant*(coeff[4]-coeff[7]*coeff[5]);
121 inverse[1]=determinant*(coeff[7]*coeff[2]-coeff[1]);
122 inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[4]*coeff[2]);
123 inverse[3]=determinant*(coeff[6]*coeff[5]-coeff[3]);
124 inverse[4]=determinant*(coeff[0]-coeff[6]*coeff[2]);
125 inverse[5]=determinant*(coeff[3]*coeff[2]-coeff[0]*coeff[5]);
126 inverse[6]=determinant*(coeff[3]*coeff[7]-coeff[6]*coeff[4]);
127 inverse[7]=determinant*(coeff[6]*coeff[1]-coeff[0]*coeff[7]);
130 static inline double MagickRound(double x)
133 Round the fraction to nearest integer.
136 return((double) ((ssize_t) (x+0.5)));
137 return((double) ((ssize_t) (x-0.5)));
141 * Polynomial Term Defining Functions
143 * Order must either be an integer, or 1.5 to produce
144 * the 2 number_valuesal polynomial function...
145 * affine 1 (3) u = c0 + c1*x + c2*y
146 * bilinear 1.5 (4) u = '' + c3*x*y
147 * quadratic 2 (6) u = '' + c4*x*x + c5*y*y
148 * cubic 3 (10) u = '' + c6*x^3 + c7*x*x*y + c8*x*y*y + c9*y^3
149 * quartic 4 (15) u = '' + c10*x^4 + ... + c14*y^4
150 * quintic 5 (21) u = '' + c15*x^5 + ... + c20*y^5
151 * number in parenthesis minimum number of points needed.
152 * Anything beyond quintic, has not been implemented until
153 * a more automated way of determining terms is found.
155 * Note the slight re-ordering of the terms for a quadratic polynomial
156 * which is to allow the use of a bi-linear (order=1.5) polynomial.
157 * All the later polynomials are ordered simply from x^N to y^N
159 static size_t poly_number_terms(double order)
161 /* Return the number of terms for a 2d polynomial */
162 if ( order < 1 || order > 5 ||
163 ( order != floor(order) && (order-1.5) > MagickEpsilon) )
164 return 0; /* invalid polynomial order */
165 return((size_t) floor((order+1)*(order+2)/2));
168 static double poly_basis_fn(ssize_t n, double x, double y)
170 /* Return the result for this polynomial term */
172 case 0: return( 1.0 ); /* constant */
174 case 2: return( y ); /* affine order = 1 terms = 3 */
175 case 3: return( x*y ); /* bilinear order = 1.5 terms = 4 */
176 case 4: return( x*x );
177 case 5: return( y*y ); /* quadratic order = 2 terms = 6 */
178 case 6: return( x*x*x );
179 case 7: return( x*x*y );
180 case 8: return( x*y*y );
181 case 9: return( y*y*y ); /* cubic order = 3 terms = 10 */
182 case 10: return( x*x*x*x );
183 case 11: return( x*x*x*y );
184 case 12: return( x*x*y*y );
185 case 13: return( x*y*y*y );
186 case 14: return( y*y*y*y ); /* quartic order = 4 terms = 15 */
187 case 15: return( x*x*x*x*x );
188 case 16: return( x*x*x*x*y );
189 case 17: return( x*x*x*y*y );
190 case 18: return( x*x*y*y*y );
191 case 19: return( x*y*y*y*y );
192 case 20: return( y*y*y*y*y ); /* quintic order = 5 terms = 21 */
194 return( 0 ); /* should never happen */
196 static const char *poly_basis_str(ssize_t n)
198 /* return the result for this polynomial term */
200 case 0: return(""); /* constant */
201 case 1: return("*ii");
202 case 2: return("*jj"); /* affine order = 1 terms = 3 */
203 case 3: return("*ii*jj"); /* bilinear order = 1.5 terms = 4 */
204 case 4: return("*ii*ii");
205 case 5: return("*jj*jj"); /* quadratic order = 2 terms = 6 */
206 case 6: return("*ii*ii*ii");
207 case 7: return("*ii*ii*jj");
208 case 8: return("*ii*jj*jj");
209 case 9: return("*jj*jj*jj"); /* cubic order = 3 terms = 10 */
210 case 10: return("*ii*ii*ii*ii");
211 case 11: return("*ii*ii*ii*jj");
212 case 12: return("*ii*ii*jj*jj");
213 case 13: return("*ii*jj*jj*jj");
214 case 14: return("*jj*jj*jj*jj"); /* quartic order = 4 terms = 15 */
215 case 15: return("*ii*ii*ii*ii*ii");
216 case 16: return("*ii*ii*ii*ii*jj");
217 case 17: return("*ii*ii*ii*jj*jj");
218 case 18: return("*ii*ii*jj*jj*jj");
219 case 19: return("*ii*jj*jj*jj*jj");
220 case 20: return("*jj*jj*jj*jj*jj"); /* quintic order = 5 terms = 21 */
222 return( "UNKNOWN" ); /* should never happen */
224 static double poly_basis_dx(ssize_t n, double x, double y)
226 /* polynomial term for x derivative */
228 case 0: return( 0.0 ); /* constant */
229 case 1: return( 1.0 );
230 case 2: return( 0.0 ); /* affine order = 1 terms = 3 */
231 case 3: return( y ); /* bilinear order = 1.5 terms = 4 */
233 case 5: return( 0.0 ); /* quadratic order = 2 terms = 6 */
234 case 6: return( x*x );
235 case 7: return( x*y );
236 case 8: return( y*y );
237 case 9: return( 0.0 ); /* cubic order = 3 terms = 10 */
238 case 10: return( x*x*x );
239 case 11: return( x*x*y );
240 case 12: return( x*y*y );
241 case 13: return( y*y*y );
242 case 14: return( 0.0 ); /* quartic order = 4 terms = 15 */
243 case 15: return( x*x*x*x );
244 case 16: return( x*x*x*y );
245 case 17: return( x*x*y*y );
246 case 18: return( x*y*y*y );
247 case 19: return( y*y*y*y );
248 case 20: return( 0.0 ); /* quintic order = 5 terms = 21 */
250 return( 0.0 ); /* should never happen */
252 static double poly_basis_dy(ssize_t n, double x, double y)
254 /* polynomial term for y derivative */
256 case 0: return( 0.0 ); /* constant */
257 case 1: return( 0.0 );
258 case 2: return( 1.0 ); /* affine order = 1 terms = 3 */
259 case 3: return( x ); /* bilinear order = 1.5 terms = 4 */
260 case 4: return( 0.0 );
261 case 5: return( y ); /* quadratic order = 2 terms = 6 */
262 default: return( poly_basis_dx(n-1,x,y) ); /* weird but true */
264 /* NOTE: the only reason that last is not true for 'quadratic'
265 is due to the re-arrangement of terms to allow for 'bilinear'
270 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
274 + G e n e r a t e C o e f f i c i e n t s %
278 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
280 % GenerateCoefficients() takes user provided input arguments and generates
281 % the coefficients, needed to apply the specific distortion for either
282 % distorting images (generally using control points) or generating a color
283 % gradient from sparsely separated color points.
285 % The format of the GenerateCoefficients() method is:
287 % Image *GenerateCoefficients(const Image *image,DistortImageMethod method,
288 % const size_t number_arguments,const double *arguments,
289 % size_t number_values, ExceptionInfo *exception)
291 % A description of each parameter follows:
293 % o image: the image to be distorted.
295 % o method: the method of image distortion/ sparse gradient
297 % o number_arguments: the number of arguments given.
299 % o arguments: the arguments for this distortion method.
301 % o number_values: the style and format of given control points, (caller type)
302 % 0: 2 dimensional mapping of control points (Distort)
303 % Format: u,v,x,y where u,v is the 'source' of the
304 % the color to be plotted, for DistortImage()
305 % N: Interpolation of control points with N values (usally r,g,b)
306 % Format: x,y,r,g,b mapping x,y to color values r,g,b
307 % IN future, variable number of values may be given (1 to N)
309 % o exception: return any errors or warnings in this structure
311 % Note that the returned array of double values must be freed by the
312 % calling method using RelinquishMagickMemory(). This however may change in
313 % the future to require a more 'method' specific method.
315 % Because of this this method should not be classed as stable or used
316 % outside other MagickCore library methods.
319 static double *GenerateCoefficients(const Image *image,
320 DistortImageMethod *method,const size_t number_arguments,
321 const double *arguments,size_t number_values,ExceptionInfo *exception)
330 number_coeff, /* number of coefficients to return (array size) */
331 cp_size, /* number floating point numbers per control point */
332 cp_x,cp_y, /* the x,y indexes for control point */
333 cp_values; /* index of values for this control point */
334 /* number_values Number of values given per control point */
336 if ( number_values == 0 ) {
337 /* Image distortion using control points (or other distortion)
338 That is generate a mapping so that x,y->u,v given u,v,x,y
340 number_values = 2; /* special case: two values of u,v */
341 cp_values = 0; /* the values i,j are BEFORE the destination CP x,y */
342 cp_x = 2; /* location of x,y in input control values */
344 /* NOTE: cp_values, also used for later 'reverse map distort' tests */
347 cp_x = 0; /* location of x,y in input control values */
349 cp_values = 2; /* and the other values are after x,y */
350 /* Typically in this case the values are R,G,B color values */
352 cp_size = number_values+2; /* each CP defintion involves this many numbers */
354 /* If not enough control point pairs are found for specific distortions
355 fall back to Affine distortion (allowing 0 to 3 point pairs)
357 if ( number_arguments < 4*cp_size &&
358 ( *method == BilinearForwardDistortion
359 || *method == BilinearReverseDistortion
360 || *method == PerspectiveDistortion
362 *method = AffineDistortion;
366 case AffineDistortion:
367 /* also BarycentricColorInterpolate: */
368 number_coeff=3*number_values;
370 case PolynomialDistortion:
371 /* number of coefficents depend on the given polynomal 'order' */
372 if ( number_arguments <= 1 && (number_arguments-1)%cp_size != 0)
374 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
375 "InvalidArgument","%s : '%s'","Polynomial",
376 "Invalid number of args: order [CPs]...");
377 return((double *) NULL);
379 i = poly_number_terms(arguments[0]);
380 number_coeff = 2 + i*number_values;
382 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
383 "InvalidArgument","%s : '%s'","Polynomial",
384 "Invalid order, should be interger 1 to 5, or 1.5");
385 return((double *) NULL);
387 if ( number_arguments < 1+i*cp_size ) {
388 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
389 "InvalidArgument", "%s : 'require at least %.20g CPs'",
390 "Polynomial", (double) i);
391 return((double *) NULL);
394 case BilinearReverseDistortion:
395 number_coeff=4*number_values;
398 The rest are constants as they are only used for image distorts
400 case BilinearForwardDistortion:
401 number_coeff=10; /* 2*4 coeff plus 2 constants */
402 cp_x = 0; /* Reverse src/dest coords for forward mapping */
407 case QuadraterialDistortion:
408 number_coeff=19; /* BilinearForward + BilinearReverse */
411 case ShepardsDistortion:
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:
434 assert(! "Unknown Method Given"); /* just fail assertion */
437 /* allocate the array of coefficients needed */
438 coeff = (double *) AcquireQuantumMemory(number_coeff,sizeof(*coeff));
439 if (coeff == (double *) NULL) {
440 (void) ThrowMagickException(exception,GetMagickModule(),
441 ResourceLimitError,"MemoryAllocationFailed",
442 "%s", "GenerateCoefficients");
443 return((double *) NULL);
446 /* zero out coefficients array */
447 for (i=0; i < number_coeff; i++)
452 case AffineDistortion:
456 for each 'value' given
458 Input Arguments are sets of control points...
459 For Distort Images u,v, x,y ...
460 For Sparse Gradients x,y, r,g,b ...
462 if ( number_arguments%cp_size != 0 ||
463 number_arguments < cp_size ) {
464 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
465 "InvalidArgument", "%s : 'require at least %.20g CPs'",
467 coeff=(double *) RelinquishMagickMemory(coeff);
468 return((double *) NULL);
470 /* handle special cases of not enough arguments */
471 if ( number_arguments == cp_size ) {
472 /* Only 1 CP Set Given */
473 if ( cp_values == 0 ) {
474 /* image distortion - translate the image */
476 coeff[2] = arguments[0] - arguments[2];
478 coeff[5] = arguments[1] - arguments[3];
481 /* sparse gradient - use the values directly */
482 for (i=0; i<number_values; i++)
483 coeff[i*3+2] = arguments[cp_values+i];
487 /* 2 or more points (usally 3) given.
488 Solve a least squares simultaneous equation for coefficients.
498 /* create matrix, and a fake vectors matrix */
499 matrix = AcquireMagickMatrix(3UL,3UL);
500 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
501 if (matrix == (double **) NULL || vectors == (double **) NULL)
503 matrix = RelinquishMagickMatrix(matrix, 3UL);
504 vectors = (double **) RelinquishMagickMemory(vectors);
505 coeff = (double *) RelinquishMagickMemory(coeff);
506 (void) ThrowMagickException(exception,GetMagickModule(),
507 ResourceLimitError,"MemoryAllocationFailed",
508 "%s", "DistortCoefficients");
509 return((double *) NULL);
511 /* fake a number_values x3 vectors matrix from coefficients array */
512 for (i=0; i < number_values; i++)
513 vectors[i] = &(coeff[i*3]);
514 /* Add given control point pairs for least squares solving */
515 for (i=0; i < number_arguments; i+=cp_size) {
516 terms[0] = arguments[i+cp_x]; /* x */
517 terms[1] = arguments[i+cp_y]; /* y */
518 terms[2] = 1; /* 1 */
519 LeastSquaresAddTerms(matrix,vectors,terms,
520 &(arguments[i+cp_values]),3UL,number_values);
522 if ( number_arguments == 2*cp_size ) {
523 /* Only two pairs were given, but we need 3 to solve the affine.
524 Fake extra coordinates by rotating p1 around p0 by 90 degrees.
525 x2 = x0 - (y1-y0) y2 = y0 + (x1-x0)
527 terms[0] = arguments[cp_x]
528 - ( arguments[cp_size+cp_y] - arguments[cp_y] ); /* x2 */
529 terms[1] = arguments[cp_y] +
530 + ( arguments[cp_size+cp_x] - arguments[cp_x] ); /* y2 */
531 terms[2] = 1; /* 1 */
532 if ( cp_values == 0 ) {
533 /* Image Distortion - rotate the u,v coordients too */
536 uv2[0] = arguments[0] - arguments[5] + arguments[1]; /* u2 */
537 uv2[1] = arguments[1] + arguments[4] - arguments[0]; /* v2 */
538 LeastSquaresAddTerms(matrix,vectors,terms,uv2,3UL,2UL);
541 /* Sparse Gradient - use values of p0 for linear gradient */
542 LeastSquaresAddTerms(matrix,vectors,terms,
543 &(arguments[cp_values]),3UL,number_values);
546 /* Solve for LeastSquares Coefficients */
547 status=GaussJordanElimination(matrix,vectors,3UL,number_values);
548 matrix = RelinquishMagickMatrix(matrix, 3UL);
549 vectors = (double **) RelinquishMagickMemory(vectors);
550 if ( status == MagickFalse ) {
551 coeff = (double *) RelinquishMagickMemory(coeff);
552 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
553 "InvalidArgument","%s : 'Unsolvable Matrix'",
554 CommandOptionToMnemonic(MagickDistortOptions, *method) );
555 return((double *) NULL);
560 case AffineProjectionDistortion:
563 Arguments: Affine Matrix (forward mapping)
564 Arguments sx, rx, ry, sy, tx, ty
565 Where u = sx*x + ry*y + tx
568 Returns coefficients (in there inverse form) ordered as...
571 AffineProjection Distortion Notes...
572 + Will only work with a 2 number_values for Image Distortion
573 + Can not be used for generating a sparse gradient (interpolation)
576 if (number_arguments != 6) {
577 coeff = (double *) RelinquishMagickMemory(coeff);
578 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
579 "InvalidArgument","%s : 'Needs 6 coeff values'",
580 CommandOptionToMnemonic(MagickDistortOptions, *method) );
581 return((double *) NULL);
583 /* FUTURE: trap test for sx*sy-rx*ry == 0 (determinant = 0, no inverse) */
584 for(i=0; i<6UL; i++ )
585 inverse[i] = arguments[i];
586 AffineArgsToCoefficients(inverse); /* map into coefficents */
587 InvertAffineCoefficients(inverse, coeff); /* invert */
588 *method = AffineDistortion;
592 case ScaleRotateTranslateDistortion:
594 /* Scale, Rotate and Translate Distortion
595 An alternative Affine Distortion
596 Argument options, by number of arguments given:
597 7: x,y, sx,sy, a, nx,ny
604 Where actions are (in order of application)
605 x,y 'center' of transforms (default = image center)
606 sx,sy scale image by this amount (default = 1)
607 a angle of rotation (argument required)
608 nx,ny move 'center' here (default = x,y or no movement)
609 And convert to affine mapping coefficients
611 ScaleRotateTranslate Distortion Notes...
612 + Does not use a set of CPs in any normal way
613 + Will only work with a 2 number_valuesal Image Distortion
614 + Cannot be used for generating a sparse gradient (interpolation)
620 /* set default center, and default scale */
621 x = nx = (double)(image->columns)/2.0 + (double)image->page.x;
622 y = ny = (double)(image->rows)/2.0 + (double)image->page.y;
624 switch ( number_arguments ) {
626 coeff = (double *) RelinquishMagickMemory(coeff);
627 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
628 "InvalidArgument","%s : 'Needs at least 1 argument'",
629 CommandOptionToMnemonic(MagickDistortOptions, *method) );
630 return((double *) NULL);
635 sx = sy = arguments[0];
639 x = nx = arguments[0];
640 y = ny = arguments[1];
641 switch ( number_arguments ) {
646 sx = sy = arguments[2];
655 sx = sy = arguments[2];
668 coeff = (double *) RelinquishMagickMemory(coeff);
669 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
670 "InvalidArgument","%s : 'Too Many Arguments (7 or less)'",
671 CommandOptionToMnemonic(MagickDistortOptions, *method) );
672 return((double *) NULL);
676 /* Trap if sx or sy == 0 -- image is scaled out of existance! */
677 if ( fabs(sx) < MagickEpsilon || fabs(sy) < MagickEpsilon ) {
678 coeff = (double *) RelinquishMagickMemory(coeff);
679 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
680 "InvalidArgument","%s : 'Zero Scale Given'",
681 CommandOptionToMnemonic(MagickDistortOptions, *method) );
682 return((double *) NULL);
684 /* Save the given arguments as an affine distortion */
685 a=DegreesToRadians(a); cosine=cos(a); sine=sin(a);
687 *method = AffineDistortion;
690 coeff[2]=x-nx*coeff[0]-ny*coeff[1];
693 coeff[5]=y-nx*coeff[3]-ny*coeff[4];
696 case PerspectiveDistortion:
698 Perspective Distortion (a ratio of affine distortions)
700 p(x,y) c0*x + c1*y + c2
701 u = ------ = ------------------
702 r(x,y) c6*x + c7*y + 1
704 q(x,y) c3*x + c4*y + c5
705 v = ------ = ------------------
706 r(x,y) c6*x + c7*y + 1
708 c8 = Sign of 'r', or the denominator affine, for the actual image.
709 This determines what part of the distorted image is 'ground'
710 side of the horizon, the other part is 'sky' or invalid.
711 Valid values are +1.0 or -1.0 only.
713 Input Arguments are sets of control points...
714 For Distort Images u,v, x,y ...
715 For Sparse Gradients x,y, r,g,b ...
717 Perspective Distortion Notes...
718 + Can be thought of as ratio of 3 affine transformations
719 + Not separatable: r() or c6 and c7 are used by both equations
720 + All 8 coefficients must be determined simultaniously
721 + Will only work with a 2 number_valuesal Image Distortion
722 + Can not be used for generating a sparse gradient (interpolation)
723 + It is not linear, but is simple to generate an inverse
724 + All lines within an image remain lines.
725 + but distances between points may vary.
739 if ( number_arguments%cp_size != 0 ||
740 number_arguments < cp_size*4 ) {
741 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
742 "InvalidArgument", "%s : 'require at least %.20g CPs'",
743 CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0);
744 coeff=(double *) RelinquishMagickMemory(coeff);
745 return((double *) NULL);
747 /* fake 1x8 vectors matrix directly using the coefficients array */
748 vectors[0] = &(coeff[0]);
749 /* 8x8 least-squares matrix (zeroed) */
750 matrix = AcquireMagickMatrix(8UL,8UL);
751 if (matrix == (double **) NULL) {
752 (void) ThrowMagickException(exception,GetMagickModule(),
753 ResourceLimitError,"MemoryAllocationFailed",
754 "%s", "DistortCoefficients");
755 return((double *) NULL);
757 /* Add control points for least squares solving */
758 for (i=0; i < number_arguments; i+=4) {
759 terms[0]=arguments[i+cp_x]; /* c0*x */
760 terms[1]=arguments[i+cp_y]; /* c1*y */
761 terms[2]=1.0; /* c2*1 */
765 terms[6]=-terms[0]*arguments[i+cp_u]; /* 1/(c6*x) */
766 terms[7]=-terms[1]*arguments[i+cp_u]; /* 1/(c7*y) */
767 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_u]),
773 terms[3]=arguments[i+cp_x]; /* c3*x */
774 terms[4]=arguments[i+cp_y]; /* c4*y */
775 terms[5]=1.0; /* c5*1 */
776 terms[6]=-terms[3]*arguments[i+cp_v]; /* 1/(c6*x) */
777 terms[7]=-terms[4]*arguments[i+cp_v]; /* 1/(c7*y) */
778 LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_v]),
781 /* Solve for LeastSquares Coefficients */
782 status=GaussJordanElimination(matrix,vectors,8UL,1UL);
783 matrix = RelinquishMagickMatrix(matrix, 8UL);
784 if ( status == MagickFalse ) {
785 coeff = (double *) RelinquishMagickMemory(coeff);
786 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
787 "InvalidArgument","%s : 'Unsolvable Matrix'",
788 CommandOptionToMnemonic(MagickDistortOptions, *method) );
789 return((double *) NULL);
792 Calculate 9'th coefficient! The ground-sky determination.
793 What is sign of the 'ground' in r() denominator affine function?
794 Just use any valid image coordinate (first control point) in
795 destination for determination of what part of view is 'ground'.
797 coeff[8] = coeff[6]*arguments[cp_x]
798 + coeff[7]*arguments[cp_y] + 1.0;
799 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
803 case PerspectiveProjectionDistortion:
806 Arguments: Perspective Coefficents (forward mapping)
808 if (number_arguments != 8) {
809 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
810 "InvalidArgument", "%s : 'Needs 8 coefficient values'",
811 CommandOptionToMnemonic(MagickDistortOptions, *method));
812 return((double *) NULL);
814 /* FUTURE: trap test c0*c4-c3*c1 == 0 (determinate = 0, no inverse) */
815 InvertPerspectiveCoefficients(arguments, coeff);
817 Calculate 9'th coefficient! The ground-sky determination.
818 What is sign of the 'ground' in r() denominator affine function?
819 Just use any valid image cocodinate in destination for determination.
820 For a forward mapped perspective the images 0,0 coord will map to
821 c2,c5 in the distorted image, so set the sign of denominator of that.
823 coeff[8] = coeff[6]*arguments[2]
824 + coeff[7]*arguments[5] + 1.0;
825 coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
826 *method = PerspectiveDistortion;
830 case BilinearForwardDistortion:
831 case BilinearReverseDistortion:
833 /* Bilinear Distortion (Forward mapping)
834 v = c0*x + c1*y + c2*x*y + c3;
835 for each 'value' given
837 This is actually a simple polynomial Distortion! The difference
838 however is when we need to reverse the above equation to generate a
839 BilinearForwardDistortion (see below).
841 Input Arguments are sets of control points...
842 For Distort Images u,v, x,y ...
843 For Sparse Gradients x,y, r,g,b ...
854 /* check the number of arguments */
855 if ( number_arguments%cp_size != 0 ||
856 number_arguments < cp_size*4 ) {
857 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
858 "InvalidArgument", "%s : 'require at least %.20g CPs'",
859 CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0);
860 coeff=(double *) RelinquishMagickMemory(coeff);
861 return((double *) NULL);
863 /* create matrix, and a fake vectors matrix */
864 matrix = AcquireMagickMatrix(4UL,4UL);
865 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
866 if (matrix == (double **) NULL || vectors == (double **) NULL)
868 matrix = RelinquishMagickMatrix(matrix, 4UL);
869 vectors = (double **) RelinquishMagickMemory(vectors);
870 coeff = (double *) RelinquishMagickMemory(coeff);
871 (void) ThrowMagickException(exception,GetMagickModule(),
872 ResourceLimitError,"MemoryAllocationFailed",
873 "%s", "DistortCoefficients");
874 return((double *) NULL);
876 /* fake a number_values x4 vectors matrix from coefficients array */
877 for (i=0; i < number_values; i++)
878 vectors[i] = &(coeff[i*4]);
879 /* Add given control point pairs for least squares solving */
880 for (i=0; i < number_arguments; i+=cp_size) {
881 terms[0] = arguments[i+cp_x]; /* x */
882 terms[1] = arguments[i+cp_y]; /* y */
883 terms[2] = terms[0]*terms[1]; /* x*y */
884 terms[3] = 1; /* 1 */
885 LeastSquaresAddTerms(matrix,vectors,terms,
886 &(arguments[i+cp_values]),4UL,number_values);
888 /* Solve for LeastSquares Coefficients */
889 status=GaussJordanElimination(matrix,vectors,4UL,number_values);
890 matrix = RelinquishMagickMatrix(matrix, 4UL);
891 vectors = (double **) RelinquishMagickMemory(vectors);
892 if ( status == MagickFalse ) {
893 coeff = (double *) RelinquishMagickMemory(coeff);
894 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
895 "InvalidArgument","%s : 'Unsolvable Matrix'",
896 CommandOptionToMnemonic(MagickDistortOptions, *method) );
897 return((double *) NULL);
899 if ( *method == BilinearForwardDistortion ) {
900 /* Bilinear Forward Mapped Distortion
902 The above least-squares solved for coefficents but in the forward
903 direction, due to changes to indexing constants.
905 i = c0*x + c1*y + c2*x*y + c3;
906 j = c4*x + c5*y + c6*x*y + c7;
908 where i,j are in the destination image, NOT the source.
910 Reverse Pixel mapping however needs to use reverse of these
911 functions. It required a full page of algbra to work out the
912 reversed mapping formula, but resolves down to the following...
915 c9 = 2*(c2*c5-c1*c6); // '2*a' in the quadratic formula
917 i = i - c3; j = j - c7;
918 b = c6*i - c2*j + c8; // So that a*y^2 + b*y + c == 0
919 c = c4*i - c0*j; // y = ( -b +- sqrt(bb - 4ac) ) / (2*a)
923 y = ( -b + sqrt(r) ) / c9;
927 x = ( i - c1*y) / ( c1 - c2*y );
929 NB: if 'r' is negative there is no solution!
930 NB: the sign of the sqrt() should be negative if image becomes
931 flipped or flopped, or crosses over itself.
932 NB: techniqually coefficient c5 is not needed, anymore,
933 but kept for completness.
935 See Anthony Thyssen <A.Thyssen@griffith.edu.au>
936 or Fred Weinhaus <fmw@alink.net> for more details.
939 coeff[8] = coeff[0]*coeff[5] - coeff[1]*coeff[4];
940 coeff[9] = 2*(coeff[2]*coeff[5] - coeff[1]*coeff[6]);
945 case QuadrilateralDistortion:
947 /* Map a Quadrilateral to a unit square using BilinearReverse
948 Then map that unit square back to the final Quadrilateral
949 using BilinearForward.
951 Input Arguments are sets of control points...
952 For Distort Images u,v, x,y ...
953 For Sparse Gradients x,y, r,g,b ...
956 /* UNDER CONSTRUCTION */
961 case PolynomialDistortion:
963 /* Polynomial Distortion
965 First two coefficents are used to hole global polynomal information
966 c0 = Order of the polynimial being created
967 c1 = number_of_terms in one polynomial equation
969 Rest of the coefficients map to the equations....
970 v = c0 + c1*x + c2*y + c3*x*y + c4*x^2 + c5*y^2 + c6*x^3 + ...
971 for each 'value' (number_values of them) given.
972 As such total coefficients = 2 + number_terms * number_values
974 Input Arguments are sets of control points...
975 For Distort Images order [u,v, x,y] ...
976 For Sparse Gradients order [x,y, r,g,b] ...
978 Polynomial Distortion Notes...
979 + UNDER DEVELOPMENT -- Do not expect this to remain as is.
980 + Currently polynomial is a reversed mapped distortion.
981 + Order 1.5 is fudged to map into a bilinear distortion.
982 though it is not the same order as that distortion.
990 nterms; /* number of polynomial terms per number_values */
998 /* first two coefficients hold polynomial order information */
999 coeff[0] = arguments[0];
1000 coeff[1] = (double) poly_number_terms(arguments[0]);
1001 nterms = (size_t) coeff[1];
1003 /* create matrix, a fake vectors matrix, and least sqs terms */
1004 matrix = AcquireMagickMatrix(nterms,nterms);
1005 vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
1006 terms = (double *) AcquireQuantumMemory(nterms, sizeof(*terms));
1007 if (matrix == (double **) NULL ||
1008 vectors == (double **) NULL ||
1009 terms == (double *) NULL )
1011 matrix = RelinquishMagickMatrix(matrix, nterms);
1012 vectors = (double **) RelinquishMagickMemory(vectors);
1013 terms = (double *) RelinquishMagickMemory(terms);
1014 coeff = (double *) RelinquishMagickMemory(coeff);
1015 (void) ThrowMagickException(exception,GetMagickModule(),
1016 ResourceLimitError,"MemoryAllocationFailed",
1017 "%s", "DistortCoefficients");
1018 return((double *) NULL);
1020 /* fake a number_values x3 vectors matrix from coefficients array */
1021 for (i=0; i < number_values; i++)
1022 vectors[i] = &(coeff[2+i*nterms]);
1023 /* Add given control point pairs for least squares solving */
1024 for (i=1; i < number_arguments; i+=cp_size) { /* NB: start = 1 not 0 */
1025 for (j=0; j < (ssize_t) nterms; j++)
1026 terms[j] = poly_basis_fn(j,arguments[i+cp_x],arguments[i+cp_y]);
1027 LeastSquaresAddTerms(matrix,vectors,terms,
1028 &(arguments[i+cp_values]),nterms,number_values);
1030 terms = (double *) RelinquishMagickMemory(terms);
1031 /* Solve for LeastSquares Coefficients */
1032 status=GaussJordanElimination(matrix,vectors,nterms,number_values);
1033 matrix = RelinquishMagickMatrix(matrix, nterms);
1034 vectors = (double **) RelinquishMagickMemory(vectors);
1035 if ( status == MagickFalse ) {
1036 coeff = (double *) RelinquishMagickMemory(coeff);
1037 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1038 "InvalidArgument","%s : 'Unsolvable Matrix'",
1039 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1040 return((double *) NULL);
1047 Args: arc_width rotate top_edge_radius bottom_edge_radius
1048 All but first argument are optional
1049 arc_width The angle over which to arc the image side-to-side
1050 rotate Angle to rotate image from vertical center
1051 top_radius Set top edge of source image at this radius
1052 bottom_radius Set bootom edge to this radius (radial scaling)
1054 By default, if the radii arguments are nor provided the image radius
1055 is calculated so the horizontal center-line is fits the given arc
1058 The output image size is ALWAYS adjusted to contain the whole image,
1059 and an offset is given to position image relative to the 0,0 point of
1060 the origin, allowing users to use relative positioning onto larger
1061 background (via -flatten).
1063 The arguments are converted to these coefficients
1064 c0: angle for center of source image
1065 c1: angle scale for mapping to source image
1066 c2: radius for top of source image
1067 c3: radius scale for mapping source image
1068 c4: centerline of arc within source image
1070 Note the coefficients use a center angle, so asymptotic join is
1071 furthest from both sides of the source image. This also means that
1072 for arc angles greater than 360 the sides of the image will be
1075 Arc Distortion Notes...
1076 + Does not use a set of CPs
1077 + Will only work with Image Distortion
1078 + Can not be used for generating a sparse gradient (interpolation)
1080 if ( number_arguments >= 1 && arguments[0] < MagickEpsilon ) {
1081 coeff = (double *) RelinquishMagickMemory(coeff);
1082 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1083 "InvalidArgument","%s : 'Arc Angle Too Small'",
1084 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1085 return((double *) NULL);
1087 if ( number_arguments >= 3 && arguments[2] < MagickEpsilon ) {
1088 coeff = (double *) RelinquishMagickMemory(coeff);
1089 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1090 "InvalidArgument","%s : 'Outer Radius Too Small'",
1091 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1092 return((double *) NULL);
1094 coeff[0] = -MagickPI2; /* -90, place at top! */
1095 if ( number_arguments >= 1 )
1096 coeff[1] = DegreesToRadians(arguments[0]);
1098 coeff[1] = MagickPI2; /* zero arguments - center is at top */
1099 if ( number_arguments >= 2 )
1100 coeff[0] += DegreesToRadians(arguments[1]);
1101 coeff[0] /= Magick2PI; /* normalize radians */
1102 coeff[0] -= MagickRound(coeff[0]);
1103 coeff[0] *= Magick2PI; /* de-normalize back to radians */
1104 coeff[3] = (double)image->rows-1;
1105 coeff[2] = (double)image->columns/coeff[1] + coeff[3]/2.0;
1106 if ( number_arguments >= 3 ) {
1107 if ( number_arguments >= 4 )
1108 coeff[3] = arguments[2] - arguments[3];
1110 coeff[3] *= arguments[2]/coeff[2];
1111 coeff[2] = arguments[2];
1113 coeff[4] = ((double)image->columns-1.0)/2.0;
1117 case PolarDistortion:
1118 case DePolarDistortion:
1120 /* (De)Polar Distortion (same set of arguments)
1121 Args: Rmax, Rmin, Xcenter,Ycenter, Afrom,Ato
1122 DePolar can also have the extra arguments of Width, Height
1124 Coefficients 0 to 5 is the sanatized version first 6 input args
1125 Coefficient 6 is the angle to coord ratio and visa-versa
1126 Coefficient 7 is the radius to coord ratio and visa-versa
1128 WARNING: It is posible for Radius max<min and/or Angle from>to
1130 if ( number_arguments == 3
1131 || ( number_arguments > 6 && *method == PolarDistortion )
1132 || number_arguments > 8 ) {
1133 (void) ThrowMagickException(exception,GetMagickModule(),
1134 OptionError,"InvalidArgument", "%s : number of arguments",
1135 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1136 coeff=(double *) RelinquishMagickMemory(coeff);
1137 return((double *) NULL);
1139 /* Rmax - if 0 calculate appropriate value */
1140 if ( number_arguments >= 1 )
1141 coeff[0] = arguments[0];
1144 /* Rmin - usally 0 */
1145 coeff[1] = number_arguments >= 2 ? arguments[1] : 0.0;
1147 if ( number_arguments >= 4 ) {
1148 coeff[2] = arguments[2];
1149 coeff[3] = arguments[3];
1151 else { /* center of actual image */
1152 coeff[2] = (double)(image->columns)/2.0+image->page.x;
1153 coeff[3] = (double)(image->rows)/2.0+image->page.y;
1155 /* Angle from,to - about polar center 0 is downward */
1156 coeff[4] = -MagickPI;
1157 if ( number_arguments >= 5 )
1158 coeff[4] = DegreesToRadians(arguments[4]);
1159 coeff[5] = coeff[4];
1160 if ( number_arguments >= 6 )
1161 coeff[5] = DegreesToRadians(arguments[5]);
1162 if ( fabs(coeff[4]-coeff[5]) < MagickEpsilon )
1163 coeff[5] += Magick2PI; /* same angle is a full circle */
1164 /* if radius 0 or negative, its a special value... */
1165 if ( coeff[0] < MagickEpsilon ) {
1166 /* Use closest edge if radius == 0 */
1167 if ( fabs(coeff[0]) < MagickEpsilon ) {
1168 coeff[0]=MagickMin(fabs(coeff[2]-image->page.x),
1169 fabs(coeff[3]-image->page.y));
1170 coeff[0]=MagickMin(coeff[0],
1171 fabs(coeff[2]-image->page.x-image->columns));
1172 coeff[0]=MagickMin(coeff[0],
1173 fabs(coeff[3]-image->page.y-image->rows));
1175 /* furthest diagonal if radius == -1 */
1176 if ( fabs(-1.0-coeff[0]) < MagickEpsilon ) {
1178 rx = coeff[2]-image->page.x;
1179 ry = coeff[3]-image->page.y;
1180 coeff[0] = rx*rx+ry*ry;
1181 ry = coeff[3]-image->page.y-image->rows;
1182 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1183 rx = coeff[2]-image->page.x-image->columns;
1184 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1185 ry = coeff[3]-image->page.y;
1186 coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
1187 coeff[0] = sqrt(coeff[0]);
1190 /* IF Rmax <= 0 or Rmin < 0 OR Rmax < Rmin, THEN error */
1191 if ( coeff[0] < MagickEpsilon || coeff[1] < -MagickEpsilon
1192 || (coeff[0]-coeff[1]) < MagickEpsilon ) {
1193 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1194 "InvalidArgument", "%s : Invalid Radius",
1195 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1196 coeff=(double *) RelinquishMagickMemory(coeff);
1197 return((double *) NULL);
1199 /* converstion ratios */
1200 if ( *method == PolarDistortion ) {
1201 coeff[6]=(double) image->columns/(coeff[5]-coeff[4]);
1202 coeff[7]=(double) image->rows/(coeff[0]-coeff[1]);
1204 else { /* *method == DePolarDistortion */
1205 coeff[6]=(coeff[5]-coeff[4])/image->columns;
1206 coeff[7]=(coeff[0]-coeff[1])/image->rows;
1210 case BarrelDistortion:
1211 case BarrelInverseDistortion:
1213 /* Barrel Distortion
1214 Rs=(A*Rd^3 + B*Rd^2 + C*Rd + D)*Rd
1215 BarrelInv Distortion
1216 Rs=Rd/(A*Rd^3 + B*Rd^2 + C*Rd + D)
1218 Where Rd is the normalized radius from corner to middle of image
1219 Input Arguments are one of the following forms (number of arguments)...
1224 8: Ax,Bx,Cx,Dx Ay,By,Cy,Dy
1225 10: Ax,Bx,Cx,Dx Ay,By,Cy,Dy X,Y
1227 Returns 10 coefficent values, which are de-normalized (pixel scale)
1228 Ax, Bx, Cx, Dx, Ay, By, Cy, Dy, Xc, Yc
1230 /* Radius de-normalization scaling factor */
1232 rscale = 2.0/MagickMin((double) image->columns,(double) image->rows);
1234 /* sanity check number of args must = 3,4,5,6,8,10 or error */
1235 if ( (number_arguments < 3) || (number_arguments == 7) ||
1236 (number_arguments == 9) || (number_arguments > 10) )
1238 coeff=(double *) RelinquishMagickMemory(coeff);
1239 (void) ThrowMagickException(exception,GetMagickModule(),
1240 OptionError,"InvalidArgument", "%s : number of arguments",
1241 CommandOptionToMnemonic(MagickDistortOptions, *method) );
1242 return((double *) NULL);
1244 /* A,B,C,D coefficients */
1245 coeff[0] = arguments[0];
1246 coeff[1] = arguments[1];
1247 coeff[2] = arguments[2];
1248 if ((number_arguments == 3) || (number_arguments == 5) )
1249 coeff[3] = 1.0 - coeff[0] - coeff[1] - coeff[2];
1251 coeff[3] = arguments[3];
1252 /* de-normalize the coefficients */
1253 coeff[0] *= pow(rscale,3.0);
1254 coeff[1] *= rscale*rscale;
1256 /* Y coefficients: as given OR same as X coefficients */
1257 if ( number_arguments >= 8 ) {
1258 coeff[4] = arguments[4] * pow(rscale,3.0);
1259 coeff[5] = arguments[5] * rscale*rscale;
1260 coeff[6] = arguments[6] * rscale;
1261 coeff[7] = arguments[7];
1264 coeff[4] = coeff[0];
1265 coeff[5] = coeff[1];
1266 coeff[6] = coeff[2];
1267 coeff[7] = coeff[3];
1269 /* X,Y Center of Distortion (image coodinates) */
1270 if ( number_arguments == 5 ) {
1271 coeff[8] = arguments[3];
1272 coeff[9] = arguments[4];
1274 else if ( number_arguments == 6 ) {
1275 coeff[8] = arguments[4];
1276 coeff[9] = arguments[5];
1278 else if ( number_arguments == 10 ) {
1279 coeff[8] = arguments[8];
1280 coeff[9] = arguments[9];
1283 /* center of the image provided (image coodinates) */
1284 coeff[8] = (double)image->columns/2.0 + image->page.x;
1285 coeff[9] = (double)image->rows/2.0 + image->page.y;
1289 case ShepardsDistortion:
1291 /* Shepards Distortion input arguments are the coefficents!
1292 Just check the number of arguments is valid!
1293 Args: u1,v1, x1,y1, ...
1294 OR : u1,v1, r1,g1,c1, ...
1296 if ( number_arguments%cp_size != 0 ||
1297 number_arguments < cp_size ) {
1298 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1299 "InvalidArgument", "%s : 'require at least %.20g CPs'",
1300 CommandOptionToMnemonic(MagickDistortOptions, *method), 1.0);
1301 coeff=(double *) RelinquishMagickMemory(coeff);
1302 return((double *) NULL);
1309 /* you should never reach this point */
1310 assert(! "No Method Handler"); /* just fail assertion */
1311 return((double *) NULL);
1315 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1319 + D i s t o r t R e s i z e I m a g e %
1323 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1325 % DistortResizeImage() resize image using the equivalent but slower image
1326 % distortion operator. The filter is applied using a EWA cylindrical
1327 % resampling. But like resize the final image size is limited to whole pixels
1328 % with no effects by virtual-pixels on the result.
1330 % Note that images containing a transparency channel will be twice as slow to
1331 % resize as images one without transparency.
1333 % The format of the DistortResizeImage method is:
1335 % Image *AdaptiveResizeImage(const Image *image,const size_t columns,
1336 % const size_t rows,ExceptionInfo *exception)
1338 % A description of each parameter follows:
1340 % o image: the image.
1342 % o columns: the number of columns in the resized image.
1344 % o rows: the number of rows in the resized image.
1346 % o exception: return any errors or warnings in this structure.
1349 MagickExport Image *DistortResizeImage(const Image *image,
1350 const size_t columns,const size_t rows,ExceptionInfo *exception)
1352 #define DistortResizeImageTag "Distort/Image"
1368 Distort resize image.
1370 assert(image != (const Image *) NULL);
1371 assert(image->signature == MagickSignature);
1372 if (image->debug != MagickFalse)
1373 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1374 assert(exception != (ExceptionInfo *) NULL);
1375 assert(exception->signature == MagickSignature);
1376 if ((columns == 0) || (rows == 0))
1377 return((Image *) NULL);
1378 /* Do not short-circuit this resize if final image size is unchanged */
1380 (void) SetImageVirtualPixelMethod(image,TransparentVirtualPixelMethod);
1382 (void) ResetMagickMemory(distort_args,0,12*sizeof(double));
1383 distort_args[4]=(double) image->columns;
1384 distort_args[6]=(double) columns;
1385 distort_args[9]=(double) image->rows;
1386 distort_args[11]=(double) rows;
1388 vp_save=GetImageVirtualPixelMethod(image);
1390 tmp_image=CloneImage(image,0,0,MagickTrue,exception);
1391 if ( tmp_image == (Image *) NULL )
1392 return((Image *) NULL);
1393 (void) SetImageVirtualPixelMethod(tmp_image,TransparentVirtualPixelMethod);
1395 if (image->matte == MagickFalse)
1398 Image has not transparency channel, so we free to use it
1400 (void) SetImageAlphaChannel(tmp_image,SetAlphaChannel);
1401 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1402 MagickTrue,exception),
1404 tmp_image=DestroyImage(tmp_image);
1405 if ( resize_image == (Image *) NULL )
1406 return((Image *) NULL);
1408 (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel);
1409 InheritException(exception,&image->exception);
1414 Image has transparency so handle colors and alpha separatly.
1415 Basically we need to separate Virtual-Pixel alpha in the resized
1416 image, so only the actual original images alpha channel is used.
1421 /* distort alpha channel separatally */
1422 (void) SeparateImageChannel(tmp_image,TrueAlphaChannel);
1423 (void) SetImageAlphaChannel(tmp_image,OpaqueAlphaChannel);
1424 resize_alpha=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1425 MagickTrue,exception),
1426 tmp_image=DestroyImage(tmp_image);
1427 if ( resize_alpha == (Image *) NULL )
1428 return((Image *) NULL);
1430 /* distort the actual image containing alpha + VP alpha */
1431 tmp_image=CloneImage(image,0,0,MagickTrue,exception);
1432 if ( tmp_image == (Image *) NULL )
1433 return((Image *) NULL);
1434 (void) SetImageVirtualPixelMethod(tmp_image,
1435 TransparentVirtualPixelMethod);
1436 resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1437 MagickTrue,exception),
1438 tmp_image=DestroyImage(tmp_image);
1439 if ( resize_image == (Image *) NULL)
1441 resize_alpha=DestroyImage(resize_alpha);
1442 return((Image *) NULL);
1445 /* replace resize images alpha with the separally distorted alpha */
1446 (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel);
1447 (void) SetImageAlphaChannel(resize_alpha,DeactivateAlphaChannel);
1448 (void) CompositeImage(resize_image,CopyOpacityCompositeOp,resize_alpha,
1450 InheritException(exception,&resize_image->exception);
1451 resize_alpha=DestroyImage(resize_alpha);
1453 (void) SetImageVirtualPixelMethod(resize_image,vp_save);
1456 Clean up the results of the Distortion
1458 crop_area.width=columns;
1459 crop_area.height=rows;
1463 tmp_image=resize_image;
1464 resize_image=CropImage(tmp_image,&crop_area,exception);
1465 tmp_image=DestroyImage(tmp_image);
1467 if ( resize_image == (Image *) NULL )
1468 return((Image *) NULL);
1470 return(resize_image);
1474 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1478 % D i s t o r t I m a g e %
1482 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1484 % DistortImage() distorts an image using various distortion methods, by
1485 % mapping color lookups of the source image to a new destination image
1486 % usally of the same size as the source image, unless 'bestfit' is set to
1489 % If 'bestfit' is enabled, and distortion allows it, the destination image is
1490 % adjusted to ensure the whole source 'image' will just fit within the final
1491 % destination image, which will be sized and offset accordingly. Also in
1492 % many cases the virtual offset of the source image will be taken into
1493 % account in the mapping.
1495 % If the '-verbose' control option has been set print to standard error the
1496 % equicelent '-fx' formula with coefficients for the function, if practical.
1498 % The format of the DistortImage() method is:
1500 % Image *DistortImage(const Image *image,const DistortImageMethod method,
1501 % const size_t number_arguments,const double *arguments,
1502 % MagickBooleanType bestfit, ExceptionInfo *exception)
1504 % A description of each parameter follows:
1506 % o image: the image to be distorted.
1508 % o method: the method of image distortion.
1510 % ArcDistortion always ignores source image offset, and always
1511 % 'bestfit' the destination image with the top left corner offset
1512 % relative to the polar mapping center.
1514 % Affine, Perspective, and Bilinear, do least squares fitting of the
1515 % distrotion when more than the minimum number of control point pairs
1518 % Perspective, and Bilinear, fall back to a Affine distortion when less
1519 % than 4 control point pairs are provided. While Affine distortions
1520 % let you use any number of control point pairs, that is Zero pairs is
1521 % a No-Op (viewport only) distortion, one pair is a translation and
1522 % two pairs of control points do a scale-rotate-translate, without any
1525 % o number_arguments: the number of arguments given.
1527 % o arguments: an array of floating point arguments for this method.
1529 % o bestfit: Attempt to 'bestfit' the size of the resulting image.
1530 % This also forces the resulting image to be a 'layered' virtual
1531 % canvas image. Can be overridden using 'distort:viewport' setting.
1533 % o exception: return any errors or warnings in this structure
1535 % Extra Controls from Image meta-data (artifacts)...
1538 % Output to stderr alternatives, internal coefficents, and FX
1539 % equivalents for the distortion operation (if feasible).
1540 % This forms an extra check of the distortion method, and allows users
1541 % access to the internal constants IM calculates for the distortion.
1543 % o "distort:viewport"
1544 % Directly set the output image canvas area and offest to use for the
1545 % resulting image, rather than use the original images canvas, or a
1546 % calculated 'bestfit' canvas.
1549 % Scale the size of the output canvas by this amount to provide a
1550 % method of Zooming, and for super-sampling the results.
1552 % Other settings that can effect results include
1554 % o 'interpolate' For source image lookups (scale enlargements)
1556 % o 'filter' Set filter to use for area-resampling (scale shrinking).
1557 % Set to 'point' to turn off and use 'interpolate' lookup
1562 MagickExport Image *DistortImage(const Image *image,DistortImageMethod method,
1563 const size_t number_arguments,const double *arguments,
1564 MagickBooleanType bestfit,ExceptionInfo *exception)
1566 #define DistortImageTag "Distort/Image"
1576 geometry; /* geometry of the distorted space viewport */
1581 assert(image != (Image *) NULL);
1582 assert(image->signature == MagickSignature);
1583 if (image->debug != MagickFalse)
1584 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1585 assert(exception != (ExceptionInfo *) NULL);
1586 assert(exception->signature == MagickSignature);
1590 Handle Special Compound Distortions (in-direct distortions)
1592 if ( method == ResizeDistortion )
1594 if ( number_arguments != 2 )
1596 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1597 "InvalidArgument","%s : '%s'","Resize",
1598 "Invalid number of args: 2 only");
1599 return((Image *) NULL);
1601 distort_image=DistortResizeImage(image,(size_t)arguments[0],
1602 (size_t)arguments[1], exception);
1603 return(distort_image);
1607 Convert input arguments (usally as control points for reverse mapping)
1608 into mapping coefficients to apply the distortion.
1610 Note that some distortions are mapped to other distortions,
1611 and as such do not require specific code after this point.
1613 coeff = GenerateCoefficients(image, &method, number_arguments,
1614 arguments, 0, exception);
1615 if ( coeff == (double *) NULL )
1616 return((Image *) NULL);
1619 Determine the size and offset for a 'bestfit' destination.
1620 Usally the four corners of the source image is enough.
1623 /* default output image bounds, when no 'bestfit' is requested */
1624 geometry.width=image->columns;
1625 geometry.height=image->rows;
1629 if ( method == ArcDistortion ) {
1630 /* always use the 'best fit' viewport */
1631 bestfit = MagickTrue;
1634 /* Work out the 'best fit', (required for ArcDistortion) */
1639 s.x=s.y=min.x=max.x=min.y=max.y=0.0; /* keep compiler happy */
1641 /* defines to figure out the bounds of the distorted image */
1642 #define InitalBounds(p) \
1644 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1645 min.x = max.x = p.x; \
1646 min.y = max.y = p.y; \
1648 #define ExpandBounds(p) \
1650 /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1651 min.x = MagickMin(min.x,p.x); \
1652 max.x = MagickMax(max.x,p.x); \
1653 min.y = MagickMin(min.y,p.y); \
1654 max.y = MagickMax(max.y,p.y); \
1659 case AffineDistortion:
1660 { double inverse[6];
1661 InvertAffineCoefficients(coeff, inverse);
1662 s.x = (double) image->page.x;
1663 s.y = (double) image->page.y;
1664 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1665 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1667 s.x = (double) image->page.x+image->columns;
1668 s.y = (double) image->page.y;
1669 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1670 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1672 s.x = (double) image->page.x;
1673 s.y = (double) image->page.y+image->rows;
1674 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1675 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1677 s.x = (double) image->page.x+image->columns;
1678 s.y = (double) image->page.y+image->rows;
1679 d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1680 d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1684 case PerspectiveDistortion:
1685 { double inverse[8], scale;
1686 InvertPerspectiveCoefficients(coeff, inverse);
1687 s.x = (double) image->page.x;
1688 s.y = (double) image->page.y;
1689 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1690 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1691 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1692 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1694 s.x = (double) image->page.x+image->columns;
1695 s.y = (double) image->page.y;
1696 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1697 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1698 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1699 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1701 s.x = (double) image->page.x;
1702 s.y = (double) image->page.y+image->rows;
1703 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1704 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1705 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1706 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1708 s.x = (double) image->page.x+image->columns;
1709 s.y = (double) image->page.y+image->rows;
1710 scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1711 scale=1.0/( (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1712 d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1713 d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1719 /* Forward Map Corners */
1720 a = coeff[0]-coeff[1]/2; ca = cos(a); sa = sin(a);
1724 d.x = (coeff[2]-coeff[3])*ca;
1725 d.y = (coeff[2]-coeff[3])*sa;
1727 a = coeff[0]+coeff[1]/2; ca = cos(a); sa = sin(a);
1731 d.x = (coeff[2]-coeff[3])*ca;
1732 d.y = (coeff[2]-coeff[3])*sa;
1734 /* Orthogonal points along top of arc */
1735 for( a=(double) (ceil((double) ((coeff[0]-coeff[1]/2.0)/MagickPI2))*MagickPI2);
1736 a<(coeff[0]+coeff[1]/2.0); a+=MagickPI2 ) {
1737 ca = cos(a); sa = sin(a);
1743 Convert the angle_to_width and radius_to_height
1744 to appropriate scaling factors, to allow faster processing
1745 in the mapping function.
1747 coeff[1] = (double) (Magick2PI*image->columns/coeff[1]);
1748 coeff[3] = (double)image->rows/coeff[3];
1751 case PolarDistortion:
1753 if (number_arguments < 2)
1754 coeff[2] = coeff[3] = 0.0;
1755 min.x = coeff[2]-coeff[0];
1756 max.x = coeff[2]+coeff[0];
1757 min.y = coeff[3]-coeff[0];
1758 max.y = coeff[3]+coeff[0];
1759 /* should be about 1.0 if Rmin = 0 */
1760 coeff[7]=(double) geometry.height/(coeff[0]-coeff[1]);
1763 case DePolarDistortion:
1765 /* direct calculation as it needs to tile correctly
1766 * for reversibility in a DePolar-Polar cycle */
1767 geometry.x = geometry.y = 0;
1768 geometry.height = (size_t) ceil(coeff[0]-coeff[1]);
1769 geometry.width = (size_t)
1770 ceil((coeff[0]-coeff[1])*(coeff[5]-coeff[4])*0.5);
1773 case ShepardsDistortion:
1774 case BilinearForwardDistortion:
1775 case BilinearReverseDistortion:
1777 case QuadrilateralDistortion:
1779 case PolynomialDistortion:
1780 case BarrelDistortion:
1781 case BarrelInverseDistortion:
1783 /* no bestfit available for this distortion */
1784 bestfit = MagickFalse;
1788 /* Set the output image geometry to calculated 'bestfit'.
1789 Yes this tends to 'over do' the file image size, ON PURPOSE!
1790 Do not do this for DePolar which needs to be exact for virtual tiling.
1792 if ( bestfit && method != DePolarDistortion ) {
1793 geometry.x = (ssize_t) floor(min.x-0.5);
1794 geometry.y = (ssize_t) floor(min.y-0.5);
1795 geometry.width=(size_t) ceil(max.x-geometry.x+0.5);
1796 geometry.height=(size_t) ceil(max.y-geometry.y+0.5);
1799 /* Now that we have a new size lets some distortions to it exactly
1800 This is for correct handling of Depolar and its virtual tile handling
1802 if ( method == DePolarDistortion ) {
1803 coeff[6]=(coeff[5]-coeff[4])/geometry.width; /* changed width */
1804 coeff[7]=(coeff[0]-coeff[1])/geometry.height; /* should be about 1.0 */
1808 /* The user provided a 'viewport' expert option which may
1809 overrides some parts of the current output image geometry.
1810 For ArcDistortion, this also overrides its default 'bestfit' setting.
1812 { const char *artifact=GetImageArtifact(image,"distort:viewport");
1813 viewport_given = MagickFalse;
1814 if ( artifact != (const char *) NULL ) {
1815 (void) ParseAbsoluteGeometry(artifact,&geometry);
1816 viewport_given = MagickTrue;
1820 /* Verbose output */
1821 if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
1824 char image_gen[MaxTextExtent];
1827 /* Set destination image size and virtual offset */
1828 if ( bestfit || viewport_given ) {
1829 (void) FormatLocaleString(image_gen, MaxTextExtent," -size %.20gx%.20g "
1830 "-page %+.20g%+.20g xc: +insert \\\n",(double) geometry.width,
1831 (double) geometry.height,(double) geometry.x,(double) geometry.y);
1832 lookup="v.p{ xx-v.page.x-.5, yy-v.page.x-.5 }";
1835 image_gen[0] = '\0'; /* no destination to generate */
1836 lookup = "p{ xx-page.x-.5, yy-page.x-.5 }"; /* simplify lookup */
1840 case AffineDistortion:
1844 inverse = (double *) AcquireQuantumMemory(6,sizeof(*inverse));
1845 if (inverse == (double *) NULL) {
1846 coeff = (double *) RelinquishMagickMemory(coeff);
1847 (void) ThrowMagickException(exception,GetMagickModule(),
1848 ResourceLimitError,"MemoryAllocationFailed",
1849 "%s", "DistortImages");
1850 return((Image *) NULL);
1852 InvertAffineCoefficients(coeff, inverse);
1853 CoefficientsToAffineArgs(inverse);
1854 FormatLocaleFile(stderr, "Affine Projection:\n");
1855 FormatLocaleFile(stderr, " -distort AffineProjection \\\n '");
1856 for (i=0; i < 5; i++)
1857 FormatLocaleFile(stderr, "%lf,", inverse[i]);
1858 FormatLocaleFile(stderr, "%lf'\n", inverse[5]);
1859 inverse = (double *) RelinquishMagickMemory(inverse);
1861 FormatLocaleFile(stderr, "Affine Distort, FX Equivelent:\n");
1862 FormatLocaleFile(stderr, "%s", image_gen);
1863 FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1864 FormatLocaleFile(stderr, " xx=%+lf*ii %+lf*jj %+lf;\n",
1865 coeff[0], coeff[1], coeff[2]);
1866 FormatLocaleFile(stderr, " yy=%+lf*ii %+lf*jj %+lf;\n",
1867 coeff[3], coeff[4], coeff[5]);
1868 FormatLocaleFile(stderr, " %s'\n", lookup);
1873 case PerspectiveDistortion:
1877 inverse = (double *) AcquireQuantumMemory(8,sizeof(*inverse));
1878 if (inverse == (double *) NULL) {
1879 coeff = (double *) RelinquishMagickMemory(coeff);
1880 (void) ThrowMagickException(exception,GetMagickModule(),
1881 ResourceLimitError,"MemoryAllocationFailed",
1882 "%s", "DistortCoefficients");
1883 return((Image *) NULL);
1885 InvertPerspectiveCoefficients(coeff, inverse);
1886 FormatLocaleFile(stderr, "Perspective Projection:\n");
1887 FormatLocaleFile(stderr, " -distort PerspectiveProjection \\\n '");
1889 FormatLocaleFile(stderr, "%lf, ", inverse[i]);
1890 FormatLocaleFile(stderr, "\n ");
1892 FormatLocaleFile(stderr, "%lf, ", inverse[i]);
1893 FormatLocaleFile(stderr, "%lf'\n", inverse[7]);
1894 inverse = (double *) RelinquishMagickMemory(inverse);
1896 FormatLocaleFile(stderr, "Perspective Distort, FX Equivelent:\n");
1897 FormatLocaleFile(stderr, "%s", image_gen);
1898 FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1899 FormatLocaleFile(stderr, " rr=%+lf*ii %+lf*jj + 1;\n",
1900 coeff[6], coeff[7]);
1901 FormatLocaleFile(stderr, " xx=(%+lf*ii %+lf*jj %+lf)/rr;\n",
1902 coeff[0], coeff[1], coeff[2]);
1903 FormatLocaleFile(stderr, " yy=(%+lf*ii %+lf*jj %+lf)/rr;\n",
1904 coeff[3], coeff[4], coeff[5]);
1905 FormatLocaleFile(stderr, " rr%s0 ? %s : blue'\n",
1906 coeff[8] < 0 ? "<" : ">", lookup);
1910 case BilinearForwardDistortion:
1911 FormatLocaleFile(stderr, "BilinearForward Mapping Equations:\n");
1912 FormatLocaleFile(stderr, "%s", image_gen);
1913 FormatLocaleFile(stderr, " i = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
1914 coeff[0], coeff[1], coeff[2], coeff[3]);
1915 FormatLocaleFile(stderr, " j = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
1916 coeff[4], coeff[5], coeff[6], coeff[7]);
1919 FormatLocaleFile(stderr, " c8 = %+lf c9 = 2*a = %+lf;\n",
1920 coeff[8], coeff[9]);
1922 FormatLocaleFile(stderr, "BilinearForward Distort, FX Equivelent:\n");
1923 FormatLocaleFile(stderr, "%s", image_gen);
1924 FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
1925 0.5-coeff[3], 0.5-coeff[7]);
1926 FormatLocaleFile(stderr, " bb=%lf*ii %+lf*jj %+lf;\n",
1927 coeff[6], -coeff[2], coeff[8]);
1928 /* Handle Special degenerate (non-quadratic) or trapezoidal case */
1929 if ( coeff[9] != 0 ) {
1930 FormatLocaleFile(stderr, " rt=bb*bb %+lf*(%lf*ii%+lf*jj);\n",
1931 -2*coeff[9], coeff[4], -coeff[0]);
1932 FormatLocaleFile(stderr, " yy=( -bb + sqrt(rt) ) / %lf;\n",
1935 FormatLocaleFile(stderr, " yy=(%lf*ii%+lf*jj)/bb;\n",
1936 -coeff[4], coeff[0]);
1937 FormatLocaleFile(stderr, " xx=(ii %+lf*yy)/(%lf %+lf*yy);\n",
1938 -coeff[1], coeff[0], coeff[2]);
1939 if ( coeff[9] != 0 )
1940 FormatLocaleFile(stderr, " (rt < 0 ) ? red : %s'\n", lookup);
1942 FormatLocaleFile(stderr, " %s'\n", lookup);
1945 case BilinearReverseDistortion:
1947 FormatLocaleFile(stderr, "Polynomial Projection Distort:\n");
1948 FormatLocaleFile(stderr, " -distort PolynomialProjection \\\n");
1949 FormatLocaleFile(stderr, " '1.5, %lf, %lf, %lf, %lf,\n",
1950 coeff[3], coeff[0], coeff[1], coeff[2]);
1951 FormatLocaleFile(stderr, " %lf, %lf, %lf, %lf'\n",
1952 coeff[7], coeff[4], coeff[5], coeff[6]);
1954 FormatLocaleFile(stderr, "BilinearReverse Distort, FX Equivelent:\n");
1955 FormatLocaleFile(stderr, "%s", image_gen);
1956 FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1957 FormatLocaleFile(stderr, " xx=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
1958 coeff[0], coeff[1], coeff[2], coeff[3]);
1959 FormatLocaleFile(stderr, " yy=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
1960 coeff[4], coeff[5], coeff[6], coeff[7]);
1961 FormatLocaleFile(stderr, " %s'\n", lookup);
1964 case PolynomialDistortion:
1966 size_t nterms = (size_t) coeff[1];
1967 FormatLocaleFile(stderr, "Polynomial (order %lg, terms %lu), FX Equivelent\n",
1968 coeff[0],(unsigned long) nterms);
1969 FormatLocaleFile(stderr, "%s", image_gen);
1970 FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1971 FormatLocaleFile(stderr, " xx =");
1972 for (i=0; i<(ssize_t) nterms; i++) {
1973 if ( i != 0 && i%4 == 0 ) FormatLocaleFile(stderr, "\n ");
1974 FormatLocaleFile(stderr, " %+lf%s", coeff[2+i],
1977 FormatLocaleFile(stderr, ";\n yy =");
1978 for (i=0; i<(ssize_t) nterms; i++) {
1979 if ( i != 0 && i%4 == 0 ) FormatLocaleFile(stderr, "\n ");
1980 FormatLocaleFile(stderr, " %+lf%s", coeff[2+i+nterms],
1983 FormatLocaleFile(stderr, ";\n %s'\n", lookup);
1988 FormatLocaleFile(stderr, "Arc Distort, Internal Coefficients:\n");
1989 for ( i=0; i<5; i++ )
1990 FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
1991 FormatLocaleFile(stderr, "Arc Distort, FX Equivelent:\n");
1992 FormatLocaleFile(stderr, "%s", image_gen);
1993 FormatLocaleFile(stderr, " -fx 'ii=i+page.x; jj=j+page.y;\n");
1994 FormatLocaleFile(stderr, " xx=(atan2(jj,ii)%+lf)/(2*pi);\n",
1996 FormatLocaleFile(stderr, " xx=xx-round(xx);\n");
1997 FormatLocaleFile(stderr, " xx=xx*%lf %+lf;\n",
1998 coeff[1], coeff[4]);
1999 FormatLocaleFile(stderr, " yy=(%lf - hypot(ii,jj)) * %lf;\n",
2000 coeff[2], coeff[3]);
2001 FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}'\n");
2004 case PolarDistortion:
2006 FormatLocaleFile(stderr, "Polar Distort, Internal Coefficents\n");
2007 for ( i=0; i<8; i++ )
2008 FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
2009 FormatLocaleFile(stderr, "Polar Distort, FX Equivelent:\n");
2010 FormatLocaleFile(stderr, "%s", image_gen);
2011 FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
2012 -coeff[2], -coeff[3]);
2013 FormatLocaleFile(stderr, " xx=(atan2(ii,jj)%+lf)/(2*pi);\n",
2014 -(coeff[4]+coeff[5])/2 );
2015 FormatLocaleFile(stderr, " xx=xx-round(xx);\n");
2016 FormatLocaleFile(stderr, " xx=xx*2*pi*%lf + v.w/2;\n",
2018 FormatLocaleFile(stderr, " yy=(hypot(ii,jj)%+lf)*%lf;\n",
2019 -coeff[1], coeff[7] );
2020 FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}'\n");
2023 case DePolarDistortion:
2025 FormatLocaleFile(stderr, "DePolar Distort, Internal Coefficents\n");
2026 for ( i=0; i<8; i++ )
2027 FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
2028 FormatLocaleFile(stderr, "DePolar Distort, FX Equivelent:\n");
2029 FormatLocaleFile(stderr, "%s", image_gen);
2030 FormatLocaleFile(stderr, " -fx 'aa=(i+.5)*%lf %+lf;\n", coeff[6], -coeff[4] );
2031 FormatLocaleFile(stderr, " rr=(j+.5)*%lf %+lf;\n", coeff[7], +coeff[1] );
2032 FormatLocaleFile(stderr, " xx=rr*sin(aa) %+lf;\n", coeff[2] );
2033 FormatLocaleFile(stderr, " yy=rr*cos(aa) %+lf;\n", coeff[3] );
2034 FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}'\n");
2037 case BarrelDistortion:
2038 case BarrelInverseDistortion:
2040 /* NOTE: This does the barrel roll in pixel coords not image coords
2041 ** The internal distortion must do it in image coordinates,
2042 ** so that is what the center coeff (8,9) is given in.
2044 xc = ((double)image->columns-1.0)/2.0 + image->page.x;
2045 yc = ((double)image->rows-1.0)/2.0 + image->page.y;
2046 FormatLocaleFile(stderr, "Barrel%s Distort, FX Equivelent:\n",
2047 method == BarrelDistortion ? "" : "Inv");
2048 FormatLocaleFile(stderr, "%s", image_gen);
2049 if ( fabs(coeff[8]-xc-0.5) < 0.1 && fabs(coeff[9]-yc-0.5) < 0.1 )
2050 FormatLocaleFile(stderr, " -fx 'xc=(w-1)/2; yc=(h-1)/2;\n");
2052 FormatLocaleFile(stderr, " -fx 'xc=%lf; yc=%lf;\n",
2053 coeff[8]-0.5, coeff[9]-0.5);
2054 FormatLocaleFile(stderr,
2055 " ii=i-xc; jj=j-yc; rr=hypot(ii,jj);\n");
2056 FormatLocaleFile(stderr, " ii=ii%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2057 method == BarrelDistortion ? "*" : "/",
2058 coeff[0],coeff[1],coeff[2],coeff[3]);
2059 FormatLocaleFile(stderr, " jj=jj%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2060 method == BarrelDistortion ? "*" : "/",
2061 coeff[4],coeff[5],coeff[6],coeff[7]);
2062 FormatLocaleFile(stderr, " v.p{fx*ii+xc,fy*jj+yc}'\n");
2069 /* The user provided a 'scale' expert option will scale the
2070 output image size, by the factor given allowing for super-sampling
2071 of the distorted image space. Any scaling factors must naturally
2072 be halved as a result.
2074 { const char *artifact;
2075 artifact=GetImageArtifact(image,"distort:scale");
2076 output_scaling = 1.0;
2077 if (artifact != (const char *) NULL) {
2078 output_scaling = fabs(InterpretLocaleValue(artifact,(char **) NULL));
2079 geometry.width *= output_scaling;
2080 geometry.height *= output_scaling;
2081 geometry.x *= output_scaling;
2082 geometry.y *= output_scaling;
2083 if ( output_scaling < 0.1 ) {
2084 coeff = (double *) RelinquishMagickMemory(coeff);
2085 (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
2086 "InvalidArgument","%s", "-set option:distort:scale" );
2087 return((Image *) NULL);
2089 output_scaling = 1/output_scaling;
2092 #define ScaleFilter(F,A,B,C,D) \
2093 ScaleResampleFilter( (F), \
2094 output_scaling*(A), output_scaling*(B), \
2095 output_scaling*(C), output_scaling*(D) )
2098 Initialize the distort image attributes.
2100 distort_image=CloneImage(image,geometry.width,geometry.height,MagickTrue,
2102 if (distort_image == (Image *) NULL)
2103 return((Image *) NULL);
2104 /* if image is ColorMapped - change it to DirectClass */
2105 if (SetImageStorageClass(distort_image,DirectClass) == MagickFalse)
2107 InheritException(exception,&distort_image->exception);
2108 distort_image=DestroyImage(distort_image);
2109 return((Image *) NULL);
2111 distort_image->page.x=geometry.x;
2112 distort_image->page.y=geometry.y;
2113 if (distort_image->background_color.opacity != OpaqueOpacity)
2114 distort_image->matte=MagickTrue;
2116 { /* ----- MAIN CODE -----
2117 Sample the source image to each pixel in the distort image.
2132 **restrict resample_filter;
2139 GetMagickPixelPacket(distort_image,&zero);
2140 resample_filter=AcquireResampleFilterThreadSet(image,
2141 UndefinedVirtualPixelMethod,MagickFalse,exception);
2142 distort_view=AcquireCacheView(distort_image);
2143 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2144 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2146 for (j=0; j < (ssize_t) distort_image->rows; j++)
2149 id = GetOpenMPThreadId();
2152 validity; /* how mathematically valid is this the mapping */
2158 pixel, /* pixel color to assign to distorted image */
2159 invalid; /* the color to assign when distort result is invalid */
2163 s; /* transform destination image x,y to source image x,y */
2165 register IndexPacket
2171 register PixelPacket
2174 q=QueueCacheViewAuthenticPixels(distort_view,0,j,distort_image->columns,1,
2176 if (q == (PixelPacket *) NULL)
2181 indexes=GetCacheViewAuthenticIndexQueue(distort_view);
2184 /* Define constant scaling vectors for Affine Distortions
2185 Other methods are either variable, or use interpolated lookup
2189 case AffineDistortion:
2190 ScaleFilter( resample_filter[id],
2192 coeff[3], coeff[4] );
2198 /* Initialize default pixel validity
2199 * negative: pixel is invalid output 'matte_color'
2200 * 0.0 to 1.0: antialiased, mix with resample output
2201 * 1.0 or greater: use resampled output.
2205 GetMagickPixelPacket(distort_image,&invalid);
2206 SetMagickPixelPacket(distort_image,&distort_image->matte_color,
2207 (IndexPacket *) NULL, &invalid);
2208 if (distort_image->colorspace == CMYKColorspace)
2209 ConvertRGBToCMYK(&invalid); /* what about other color spaces? */
2211 for (i=0; i < (ssize_t) distort_image->columns; i++)
2213 /* map pixel coordinate to distortion space coordinate */
2214 d.x = (double) (geometry.x+i+0.5)*output_scaling;
2215 d.y = (double) (geometry.y+j+0.5)*output_scaling;
2216 s = d; /* default is a no-op mapping */
2219 case AffineDistortion:
2221 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2222 s.y=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2223 /* Affine partial derivitives are constant -- set above */
2226 case PerspectiveDistortion:
2229 p,q,r,abs_r,abs_c6,abs_c7,scale;
2230 /* perspective is a ratio of affines */
2231 p=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2232 q=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2233 r=coeff[6]*d.x+coeff[7]*d.y+1.0;
2234 /* Pixel Validity -- is it a 'sky' or 'ground' pixel */
2235 validity = (r*coeff[8] < 0.0) ? 0.0 : 1.0;
2236 /* Determine horizon anti-alias blending */
2238 abs_c6 = fabs(coeff[6]);
2239 abs_c7 = fabs(coeff[7]);
2240 if ( abs_c6 > abs_c7 ) {
2241 if ( abs_r < abs_c6*output_scaling )
2242 validity = 0.5 - coeff[8]*r/(coeff[6]*output_scaling);
2244 else if ( abs_r < abs_c7*output_scaling )
2245 validity = 0.5 - coeff[8]*r/(coeff[7]*output_scaling);
2246 /* Perspective Sampling Point (if valid) */
2247 if ( validity > 0.0 ) {
2248 /* divide by r affine, for perspective scaling */
2252 /* Perspective Partial Derivatives or Scaling Vectors */
2254 ScaleFilter( resample_filter[id],
2255 (r*coeff[0] - p*coeff[6])*scale,
2256 (r*coeff[1] - p*coeff[7])*scale,
2257 (r*coeff[3] - q*coeff[6])*scale,
2258 (r*coeff[4] - q*coeff[7])*scale );
2262 case BilinearReverseDistortion:
2264 /* Reversed Mapped is just a simple polynomial */
2265 s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2]*d.x*d.y+coeff[3];
2266 s.y=coeff[4]*d.x+coeff[5]*d.y
2267 +coeff[6]*d.x*d.y+coeff[7];
2268 /* Bilinear partial derivitives of scaling vectors */
2269 ScaleFilter( resample_filter[id],
2270 coeff[0] + coeff[2]*d.y,
2271 coeff[1] + coeff[2]*d.x,
2272 coeff[4] + coeff[6]*d.y,
2273 coeff[5] + coeff[6]*d.x );
2276 case BilinearForwardDistortion:
2278 /* Forward mapped needs reversed polynomial equations
2279 * which unfortunatally requires a square root! */
2281 d.x -= coeff[3]; d.y -= coeff[7];
2282 b = coeff[6]*d.x - coeff[2]*d.y + coeff[8];
2283 c = coeff[4]*d.x - coeff[0]*d.y;
2286 /* Handle Special degenerate (non-quadratic) case */
2287 if ( fabs(coeff[9]) < MagickEpsilon )
2290 c = b*b - 2*coeff[9]*c;
2294 s.y = ( -b + sqrt(c) )/coeff[9];
2296 if ( validity > 0.0 )
2297 s.x = ( d.x - coeff[1]*s.y) / ( coeff[0] + coeff[2]*s.y );
2299 /* NOTE: the sign of the square root should be -ve for parts
2300 where the source image becomes 'flipped' or 'mirrored'.
2301 FUTURE: Horizon handling
2302 FUTURE: Scaling factors or Deritives (how?)
2307 case BilinearDistortion:
2308 /* Bilinear mapping of any Quadrilateral to any Quadrilateral */
2309 /* UNDER DEVELOPMENT */
2312 case PolynomialDistortion:
2314 /* multi-ordered polynomial */
2319 nterms=(ssize_t)coeff[1];
2322 du,dv; /* the du,dv vectors from unit dx,dy -- derivatives */
2324 s.x=s.y=du.x=du.y=dv.x=dv.y=0.0;
2325 for(k=0; k < nterms; k++) {
2326 s.x += poly_basis_fn(k,d.x,d.y)*coeff[2+k];
2327 du.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k];
2328 du.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k];
2329 s.y += poly_basis_fn(k,d.x,d.y)*coeff[2+k+nterms];
2330 dv.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k+nterms];
2331 dv.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k+nterms];
2333 ScaleFilter( resample_filter[id], du.x,du.y,dv.x,dv.y );
2338 /* what is the angle and radius in the destination image */
2339 s.x = (double) ((atan2(d.y,d.x) - coeff[0])/Magick2PI);
2340 s.x -= MagickRound(s.x); /* angle */
2341 s.y = hypot(d.x,d.y); /* radius */
2343 /* Arc Distortion Partial Scaling Vectors
2344 Are derived by mapping the perpendicular unit vectors
2345 dR and dA*R*2PI rather than trying to map dx and dy
2346 The results is a very simple orthogonal aligned ellipse.
2348 if ( s.y > MagickEpsilon )
2349 ScaleFilter( resample_filter[id],
2350 (double) (coeff[1]/(Magick2PI*s.y)), 0, 0, coeff[3] );
2352 ScaleFilter( resample_filter[id],
2353 distort_image->columns*2, 0, 0, coeff[3] );
2355 /* now scale the angle and radius for source image lookup point */
2356 s.x = s.x*coeff[1] + coeff[4] + image->page.x +0.5;
2357 s.y = (coeff[2] - s.y) * coeff[3] + image->page.y;
2360 case PolarDistortion:
2361 { /* Rect/Cartesain/Cylinder to Polar View */
2364 s.x = atan2(d.x,d.y) - (coeff[4]+coeff[5])/2;
2366 s.x -= MagickRound(s.x);
2367 s.x *= Magick2PI; /* angle - relative to centerline */
2368 s.y = hypot(d.x,d.y); /* radius */
2370 /* Polar Scaling vectors are based on mapping dR and dA vectors
2371 This results in very simple orthogonal scaling vectors
2373 if ( s.y > MagickEpsilon )
2374 ScaleFilter( resample_filter[id],
2375 (double) (coeff[6]/(Magick2PI*s.y)), 0, 0, coeff[7] );
2377 ScaleFilter( resample_filter[id],
2378 distort_image->columns*2, 0, 0, coeff[7] );
2380 /* now finish mapping radius/angle to source x,y coords */
2381 s.x = s.x*coeff[6] + (double)image->columns/2.0 + image->page.x;
2382 s.y = (s.y-coeff[1])*coeff[7] + image->page.y;
2385 case DePolarDistortion:
2386 { /* Polar to Cylindical */
2387 /* ignore all destination virtual offsets */
2388 d.x = ((double)i+0.5)*output_scaling*coeff[6]-coeff[4];
2389 d.y = ((double)j+0.5)*output_scaling*coeff[7]+coeff[1];
2390 s.x = d.y*sin(d.x) + coeff[2];
2391 s.y = d.y*cos(d.x) + coeff[3];
2392 /* derivatives are usless - better to use SuperSampling */
2395 case BarrelDistortion:
2396 case BarrelInverseDistortion:
2398 double r,fx,fy,gx,gy;
2399 /* Radial Polynomial Distortion (de-normalized) */
2402 r = sqrt(d.x*d.x+d.y*d.y);
2403 if ( r > MagickEpsilon ) {
2404 fx = ((coeff[0]*r + coeff[1])*r + coeff[2])*r + coeff[3];
2405 fy = ((coeff[4]*r + coeff[5])*r + coeff[6])*r + coeff[7];
2406 gx = ((3*coeff[0]*r + 2*coeff[1])*r + coeff[2])/r;
2407 gy = ((3*coeff[4]*r + 2*coeff[5])*r + coeff[6])/r;
2408 /* adjust functions and scaling for 'inverse' form */
2409 if ( method == BarrelInverseDistortion ) {
2410 fx = 1/fx; fy = 1/fy;
2411 gx *= -fx*fx; gy *= -fy*fy;
2413 /* Set the source pixel to lookup and EWA derivative vectors */
2414 s.x = d.x*fx + coeff[8];
2415 s.y = d.y*fy + coeff[9];
2416 ScaleFilter( resample_filter[id],
2417 gx*d.x*d.x + fx, gx*d.x*d.y,
2418 gy*d.x*d.y, gy*d.y*d.y + fy );
2421 /* Special handling to avoid divide by zero when r==0
2423 ** The source and destination pixels match in this case
2424 ** which was set at the top of the loop using s = d;
2425 ** otherwise... s.x=coeff[8]; s.y=coeff[9];
2427 if ( method == BarrelDistortion )
2428 ScaleFilter( resample_filter[id],
2429 coeff[3], 0, 0, coeff[7] );
2430 else /* method == BarrelInverseDistortion */
2431 /* FUTURE, trap for D==0 causing division by zero */
2432 ScaleFilter( resample_filter[id],
2433 1.0/coeff[3], 0, 0, 1.0/coeff[7] );
2437 case ShepardsDistortion:
2438 { /* Shepards Method, or Inverse Weighted Distance for
2439 displacement around the destination image control points
2440 The input arguments are the coefficents to the function.
2441 This is more of a 'displacement' function rather than an
2442 absolute distortion function.
2449 denominator = s.x = s.y = 0;
2450 for(i=0; i<number_arguments; i+=4) {
2452 ((double)d.x-arguments[i+2])*((double)d.x-arguments[i+2])
2453 + ((double)d.y-arguments[i+3])*((double)d.y-arguments[i+3]);
2459 s.x += (arguments[ i ]-arguments[i+2])*weight;
2460 s.y += (arguments[i+1]-arguments[i+3])*weight;
2461 denominator += weight;
2468 /* We can not determine derivatives using shepards method
2469 only color interpolatation, not area-resampling */
2473 break; /* use the default no-op given above */
2475 /* map virtual canvas location back to real image coordinate */
2476 if ( bestfit && method != ArcDistortion ) {
2477 s.x -= image->page.x;
2478 s.y -= image->page.y;
2483 if ( validity <= 0.0 ) {
2484 /* result of distortion is an invalid pixel - don't resample */
2485 SetPixelPacket(distort_image,&invalid,q,indexes);
2488 /* resample the source image to find its correct color */
2489 (void) ResamplePixelColor(resample_filter[id],s.x,s.y,&pixel);
2490 /* if validity between 0.0 and 1.0 mix result with invalid pixel */
2491 if ( validity < 1.0 ) {
2492 /* Do a blend of sample color and invalid pixel */
2493 /* should this be a 'Blend', or an 'Over' compose */
2494 MagickPixelCompositeBlend(&pixel,validity,&invalid,(1.0-validity),
2497 SetPixelPacket(distort_image,&pixel,q,indexes);
2502 sync=SyncCacheViewAuthenticPixels(distort_view,exception);
2503 if (sync == MagickFalse)
2505 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2510 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2511 #pragma omp critical (MagickCore_DistortImage)
2513 proceed=SetImageProgress(image,DistortImageTag,progress++,
2515 if (proceed == MagickFalse)
2519 distort_view=DestroyCacheView(distort_view);
2520 resample_filter=DestroyResampleFilterThreadSet(resample_filter);
2522 if (status == MagickFalse)
2523 distort_image=DestroyImage(distort_image);
2526 /* Arc does not return an offset unless 'bestfit' is in effect
2527 And the user has not provided an overriding 'viewport'.
2529 if ( method == ArcDistortion && !bestfit && !viewport_given ) {
2530 distort_image->page.x = 0;
2531 distort_image->page.y = 0;
2533 coeff = (double *) RelinquishMagickMemory(coeff);
2534 return(distort_image);
2538 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2542 % S p a r s e C o l o r I m a g e %
2546 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2548 % SparseColorImage(), given a set of coordinates, interpolates the colors
2549 % found at those coordinates, across the whole image, using various methods.
2551 % The format of the SparseColorImage() method is:
2553 % Image *SparseColorImage(const Image *image,const ChannelType channel,
2554 % const SparseColorMethod method,const size_t number_arguments,
2555 % const double *arguments,ExceptionInfo *exception)
2557 % A description of each parameter follows:
2559 % o image: the image to be filled in.
2561 % o channel: Specify which color values (in RGBKA sequence) are being set.
2562 % This also determines the number of color_values in above.
2564 % o method: the method to fill in the gradient between the control points.
2566 % The methods used for SparseColor() are often simular to methods
2567 % used for DistortImage(), and even share the same code for determination
2568 % of the function coefficents, though with more dimensions (or resulting
2571 % o number_arguments: the number of arguments given.
2573 % o arguments: array of floating point arguments for this method--
2574 % x,y,color_values-- with color_values given as normalized values.
2576 % o exception: return any errors or warnings in this structure
2579 MagickExport Image *SparseColorImage(const Image *image,
2580 const ChannelType channel,const SparseColorMethod method,
2581 const size_t number_arguments,const double *arguments,
2582 ExceptionInfo *exception)
2584 #define SparseColorTag "Distort/SparseColor"
2598 assert(image != (Image *) NULL);
2599 assert(image->signature == MagickSignature);
2600 if (image->debug != MagickFalse)
2601 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2602 assert(exception != (ExceptionInfo *) NULL);
2603 assert(exception->signature == MagickSignature);
2605 /* Determine number of color values needed per control point */
2607 if ( channel & RedChannel ) number_colors++;
2608 if ( channel & GreenChannel ) number_colors++;
2609 if ( channel & BlueChannel ) number_colors++;
2610 if ( channel & IndexChannel ) number_colors++;
2611 if ( channel & OpacityChannel ) number_colors++;
2614 Convert input arguments into mapping coefficients, this this case
2615 we are mapping (distorting) colors, rather than coordinates.
2617 { DistortImageMethod
2620 distort_method=(DistortImageMethod) method;
2621 if ( distort_method >= SentinelDistortion )
2622 distort_method = ShepardsDistortion; /* Pretend to be Shepards */
2623 coeff = GenerateCoefficients(image, &distort_method, number_arguments,
2624 arguments, number_colors, exception);
2625 if ( coeff == (double *) NULL )
2626 return((Image *) NULL);
2628 Note some Distort Methods may fall back to other simpler methods,
2629 Currently the only fallback of concern is Bilinear to Affine
2630 (Barycentric), which is alaso sparse_colr method. This also ensures
2631 correct two and one color Barycentric handling.
2633 sparse_method = (SparseColorMethod) distort_method;
2634 if ( distort_method == ShepardsDistortion )
2635 sparse_method = method; /* return non-distiort methods to normal */
2638 /* Verbose output */
2639 if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
2641 switch (sparse_method) {
2642 case BarycentricColorInterpolate:
2644 register ssize_t x=0;
2645 FormatLocaleFile(stderr, "Barycentric Sparse Color:\n");
2646 if ( channel & RedChannel )
2647 FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf' \\\n",
2648 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2649 if ( channel & GreenChannel )
2650 FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf' \\\n",
2651 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2652 if ( channel & BlueChannel )
2653 FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf' \\\n",
2654 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2655 if ( channel & IndexChannel )
2656 FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf' \\\n",
2657 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2658 if ( channel & OpacityChannel )
2659 FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf' \\\n",
2660 coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2663 case BilinearColorInterpolate:
2665 register ssize_t x=0;
2666 FormatLocaleFile(stderr, "Bilinear Sparse Color\n");
2667 if ( channel & RedChannel )
2668 FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2669 coeff[ x ], coeff[x+1],
2670 coeff[x+2], coeff[x+3]),x+=4;
2671 if ( channel & GreenChannel )
2672 FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2673 coeff[ x ], coeff[x+1],
2674 coeff[x+2], coeff[x+3]),x+=4;
2675 if ( channel & BlueChannel )
2676 FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2677 coeff[ x ], coeff[x+1],
2678 coeff[x+2], coeff[x+3]),x+=4;
2679 if ( channel & IndexChannel )
2680 FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2681 coeff[ x ], coeff[x+1],
2682 coeff[x+2], coeff[x+3]),x+=4;
2683 if ( channel & OpacityChannel )
2684 FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2685 coeff[ x ], coeff[x+1],
2686 coeff[x+2], coeff[x+3]),x+=4;
2690 /* sparse color method is too complex for FX emulation */
2695 /* Generate new image for generated interpolated gradient.
2696 * ASIDE: Actually we could have just replaced the colors of the original
2697 * image, but IM Core policy, is if storage class could change then clone
2701 sparse_image=CloneImage(image,0,0,MagickTrue,exception);
2702 if (sparse_image == (Image *) NULL)
2703 return((Image *) NULL);
2704 if (SetImageStorageClass(sparse_image,DirectClass) == MagickFalse)
2705 { /* if image is ColorMapped - change it to DirectClass */
2706 InheritException(exception,&image->exception);
2707 sparse_image=DestroyImage(sparse_image);
2708 return((Image *) NULL);
2710 { /* ----- MAIN CODE ----- */
2725 sparse_view=AcquireCacheView(sparse_image);
2726 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2727 #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2729 for (j=0; j < (ssize_t) sparse_image->rows; j++)
2735 pixel; /* pixel to assign to distorted image */
2737 register IndexPacket
2743 register PixelPacket
2746 q=GetCacheViewAuthenticPixels(sparse_view,0,j,sparse_image->columns,
2748 if (q == (PixelPacket *) NULL)
2753 indexes=GetCacheViewAuthenticIndexQueue(sparse_view);
2754 GetMagickPixelPacket(sparse_image,&pixel);
2755 for (i=0; i < (ssize_t) image->columns; i++)
2757 SetMagickPixelPacket(image,q,indexes,&pixel);
2758 switch (sparse_method)
2760 case BarycentricColorInterpolate:
2762 register ssize_t x=0;
2763 if ( channel & RedChannel )
2764 pixel.red = coeff[x]*i +coeff[x+1]*j
2766 if ( channel & GreenChannel )
2767 pixel.green = coeff[x]*i +coeff[x+1]*j
2769 if ( channel & BlueChannel )
2770 pixel.blue = coeff[x]*i +coeff[x+1]*j
2772 if ( channel & IndexChannel )
2773 pixel.index = coeff[x]*i +coeff[x+1]*j
2775 if ( channel & OpacityChannel )
2776 pixel.opacity = coeff[x]*i +coeff[x+1]*j
2780 case BilinearColorInterpolate:
2782 register ssize_t x=0;
2783 if ( channel & RedChannel )
2784 pixel.red = coeff[x]*i + coeff[x+1]*j +
2785 coeff[x+2]*i*j + coeff[x+3], x+=4;
2786 if ( channel & GreenChannel )
2787 pixel.green = coeff[x]*i + coeff[x+1]*j +
2788 coeff[x+2]*i*j + coeff[x+3], x+=4;
2789 if ( channel & BlueChannel )
2790 pixel.blue = coeff[x]*i + coeff[x+1]*j +
2791 coeff[x+2]*i*j + coeff[x+3], x+=4;
2792 if ( channel & IndexChannel )
2793 pixel.index = coeff[x]*i + coeff[x+1]*j +
2794 coeff[x+2]*i*j + coeff[x+3], x+=4;
2795 if ( channel & OpacityChannel )
2796 pixel.opacity = coeff[x]*i + coeff[x+1]*j +
2797 coeff[x+2]*i*j + coeff[x+3], x+=4;
2800 case InverseColorInterpolate:
2801 case ShepardsColorInterpolate:
2802 { /* Inverse (Squared) Distance weights average (IDW) */
2808 if ( channel & RedChannel ) pixel.red = 0.0;
2809 if ( channel & GreenChannel ) pixel.green = 0.0;
2810 if ( channel & BlueChannel ) pixel.blue = 0.0;
2811 if ( channel & IndexChannel ) pixel.index = 0.0;
2812 if ( channel & OpacityChannel ) pixel.opacity = 0.0;
2814 for(k=0; k<number_arguments; k+=2+number_colors) {
2815 register ssize_t x=(ssize_t) k+2;
2817 ((double)i-arguments[ k ])*((double)i-arguments[ k ])
2818 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
2819 if ( method == InverseColorInterpolate )
2820 weight = sqrt(weight); /* inverse, not inverse squared */
2821 weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
2822 if ( channel & RedChannel )
2823 pixel.red += arguments[x++]*weight;
2824 if ( channel & GreenChannel )
2825 pixel.green += arguments[x++]*weight;
2826 if ( channel & BlueChannel )
2827 pixel.blue += arguments[x++]*weight;
2828 if ( channel & IndexChannel )
2829 pixel.index += arguments[x++]*weight;
2830 if ( channel & OpacityChannel )
2831 pixel.opacity += arguments[x++]*weight;
2832 denominator += weight;
2834 if ( channel & RedChannel ) pixel.red /= denominator;
2835 if ( channel & GreenChannel ) pixel.green /= denominator;
2836 if ( channel & BlueChannel ) pixel.blue /= denominator;
2837 if ( channel & IndexChannel ) pixel.index /= denominator;
2838 if ( channel & OpacityChannel ) pixel.opacity /= denominator;
2841 case VoronoiColorInterpolate:
2843 { /* Just use the closest control point you can find! */
2847 minimum = MagickHuge;
2849 for(k=0; k<number_arguments; k+=2+number_colors) {
2851 ((double)i-arguments[ k ])*((double)i-arguments[ k ])
2852 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
2853 if ( distance < minimum ) {
2854 register ssize_t x=(ssize_t) k+2;
2855 if ( channel & RedChannel ) pixel.red = arguments[x++];
2856 if ( channel & GreenChannel ) pixel.green = arguments[x++];
2857 if ( channel & BlueChannel ) pixel.blue = arguments[x++];
2858 if ( channel & IndexChannel ) pixel.index = arguments[x++];
2859 if ( channel & OpacityChannel ) pixel.opacity = arguments[x++];
2866 /* set the color directly back into the source image */
2867 if ( channel & RedChannel ) pixel.red *= QuantumRange;
2868 if ( channel & GreenChannel ) pixel.green *= QuantumRange;
2869 if ( channel & BlueChannel ) pixel.blue *= QuantumRange;
2870 if ( channel & IndexChannel ) pixel.index *= QuantumRange;
2871 if ( channel & OpacityChannel ) pixel.opacity *= QuantumRange;
2872 SetPixelPacket(sparse_image,&pixel,q,indexes);
2876 sync=SyncCacheViewAuthenticPixels(sparse_view,exception);
2877 if (sync == MagickFalse)
2879 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2884 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2885 #pragma omp critical (MagickCore_SparseColorImage)
2887 proceed=SetImageProgress(image,SparseColorTag,progress++,image->rows);
2888 if (proceed == MagickFalse)
2892 sparse_view=DestroyCacheView(sparse_view);
2893 if (status == MagickFalse)
2894 sparse_image=DestroyImage(sparse_image);
2896 coeff = (double *) RelinquishMagickMemory(coeff);
2897 return(sparse_image);