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