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