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