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