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