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