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