MagickCore  7.0.7
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,12*sizeof(double));
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,TransparentVirtualPixelMethod, exception);
1548  resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1549  MagickTrue,exception),
1550  tmp_image=DestroyImage(tmp_image);
1551  if ( resize_image == (Image *) NULL)
1552  {
1553  resize_alpha=DestroyImage(resize_alpha);
1554  return((Image *) NULL);
1555  }
1556  /* replace resize images alpha with the separally distorted alpha */
1557  (void) SetImageAlphaChannel(resize_image,OffAlphaChannel,exception);
1558  (void) SetImageAlphaChannel(resize_alpha,OffAlphaChannel,exception);
1559  (void) CompositeImage(resize_image,resize_alpha,CopyAlphaCompositeOp,
1560  MagickTrue,0,0,exception);
1561  resize_alpha=DestroyImage(resize_alpha);
1562  }
1563  (void) SetImageVirtualPixelMethod(resize_image,vp_save,exception);
1564 
1565  /*
1566  Clean up the results of the Distortion
1567  */
1568  crop_area.width=columns;
1569  crop_area.height=rows;
1570  crop_area.x=0;
1571  crop_area.y=0;
1572 
1573  tmp_image=resize_image;
1574  resize_image=CropImage(tmp_image,&crop_area,exception);
1575  tmp_image=DestroyImage(tmp_image);
1576  if (resize_image != (Image *) NULL)
1577  {
1578  resize_image->alpha_trait=image->alpha_trait;
1579  resize_image->compose=image->compose;
1580  resize_image->page.width=0;
1581  resize_image->page.height=0;
1582  }
1583  return(resize_image);
1584 }
1585 
1586 /*
1587 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1588 % %
1589 % %
1590 % %
1591 % D i s t o r t I m a g e %
1592 % %
1593 % %
1594 % %
1595 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1596 %
1597 % DistortImage() distorts an image using various distortion methods, by
1598 % mapping color lookups of the source image to a new destination image
1599 % usally of the same size as the source image, unless 'bestfit' is set to
1600 % true.
1601 %
1602 % If 'bestfit' is enabled, and distortion allows it, the destination image is
1603 % adjusted to ensure the whole source 'image' will just fit within the final
1604 % destination image, which will be sized and offset accordingly. Also in
1605 % many cases the virtual offset of the source image will be taken into
1606 % account in the mapping.
1607 %
1608 % If the '-verbose' control option has been set print to standard error the
1609 % equicelent '-fx' formula with coefficients for the function, if practical.
1610 %
1611 % The format of the DistortImage() method is:
1612 %
1613 % Image *DistortImage(const Image *image,const DistortMethod method,
1614 % const size_t number_arguments,const double *arguments,
1615 % MagickBooleanType bestfit, ExceptionInfo *exception)
1616 %
1617 % A description of each parameter follows:
1618 %
1619 % o image: the image to be distorted.
1620 %
1621 % o method: the method of image distortion.
1622 %
1623 % ArcDistortion always ignores source image offset, and always
1624 % 'bestfit' the destination image with the top left corner offset
1625 % relative to the polar mapping center.
1626 %
1627 % Affine, Perspective, and Bilinear, do least squares fitting of the
1628 % distrotion when more than the minimum number of control point pairs
1629 % are provided.
1630 %
1631 % Perspective, and Bilinear, fall back to a Affine distortion when less
1632 % than 4 control point pairs are provided. While Affine distortions
1633 % let you use any number of control point pairs, that is Zero pairs is
1634 % a No-Op (viewport only) distortion, one pair is a translation and
1635 % two pairs of control points do a scale-rotate-translate, without any
1636 % shearing.
1637 %
1638 % o number_arguments: the number of arguments given.
1639 %
1640 % o arguments: an array of floating point arguments for this method.
1641 %
1642 % o bestfit: Attempt to 'bestfit' the size of the resulting image.
1643 % This also forces the resulting image to be a 'layered' virtual
1644 % canvas image. Can be overridden using 'distort:viewport' setting.
1645 %
1646 % o exception: return any errors or warnings in this structure
1647 %
1648 % Extra Controls from Image meta-data (artifacts)...
1649 %
1650 % o "verbose"
1651 % Output to stderr alternatives, internal coefficents, and FX
1652 % equivalents for the distortion operation (if feasible).
1653 % This forms an extra check of the distortion method, and allows users
1654 % access to the internal constants IM calculates for the distortion.
1655 %
1656 % o "distort:viewport"
1657 % Directly set the output image canvas area and offest to use for the
1658 % resulting image, rather than use the original images canvas, or a
1659 % calculated 'bestfit' canvas.
1660 %
1661 % o "distort:scale"
1662 % Scale the size of the output canvas by this amount to provide a
1663 % method of Zooming, and for super-sampling the results.
1664 %
1665 % Other settings that can effect results include
1666 %
1667 % o 'interpolate' For source image lookups (scale enlargements)
1668 %
1669 % o 'filter' Set filter to use for area-resampling (scale shrinking).
1670 % Set to 'point' to turn off and use 'interpolate' lookup
1671 % instead
1672 %
1673 */
1675  const size_t number_arguments,const double *arguments,
1676  MagickBooleanType bestfit,ExceptionInfo *exception)
1677 {
1678 #define DistortImageTag "Distort/Image"
1679 
1680  double
1681  *coeff,
1682  output_scaling;
1683 
1684  Image
1685  *distort_image;
1686 
1688  geometry; /* geometry of the distorted space viewport */
1689 
1691  viewport_given;
1692 
1693  assert(image != (Image *) NULL);
1694  assert(image->signature == MagickCoreSignature);
1695  if (image->debug != MagickFalse)
1696  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1697  assert(exception != (ExceptionInfo *) NULL);
1698  assert(exception->signature == MagickCoreSignature);
1699 
1700  /*
1701  Handle Special Compound Distortions
1702  */
1703  if ( method == ResizeDistortion )
1704  {
1705  if ( number_arguments != 2 )
1706  {
1708  "InvalidArgument","%s : '%s'","Resize",
1709  "Invalid number of args: 2 only");
1710  return((Image *) NULL);
1711  }
1712  distort_image=DistortResizeImage(image,(size_t)arguments[0],
1713  (size_t)arguments[1], exception);
1714  return(distort_image);
1715  }
1716 
1717  /*
1718  Convert input arguments (usually as control points for reverse mapping)
1719  into mapping coefficients to apply the distortion.
1720 
1721  Note that some distortions are mapped to other distortions,
1722  and as such do not require specific code after this point.
1723  */
1724  coeff = GenerateCoefficients(image, &method, number_arguments,
1725  arguments, 0, exception);
1726  if ( coeff == (double *) NULL )
1727  return((Image *) NULL);
1728 
1729  /*
1730  Determine the size and offset for a 'bestfit' destination.
1731  Usally the four corners of the source image is enough.
1732  */
1733 
1734  /* default output image bounds, when no 'bestfit' is requested */
1735  geometry.width=image->columns;
1736  geometry.height=image->rows;
1737  geometry.x=0;
1738  geometry.y=0;
1739 
1740  if ( method == ArcDistortion ) {
1741  bestfit = MagickTrue; /* always calculate a 'best fit' viewport */
1742  }
1743 
1744  /* Work out the 'best fit', (required for ArcDistortion) */
1745  if ( bestfit ) {
1746  PointInfo
1747  s,d,min,max; /* source, dest coords --mapping--> min, max coords */
1748 
1750  fix_bounds = MagickTrue; /* enlarge bounds for VP handling */
1751 
1752  s.x=s.y=min.x=max.x=min.y=max.y=0.0; /* keep compiler happy */
1753 
1754 /* defines to figure out the bounds of the distorted image */
1755 #define InitalBounds(p) \
1756 { \
1757  /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1758  min.x = max.x = p.x; \
1759  min.y = max.y = p.y; \
1760 }
1761 #define ExpandBounds(p) \
1762 { \
1763  /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1764  min.x = MagickMin(min.x,p.x); \
1765  max.x = MagickMax(max.x,p.x); \
1766  min.y = MagickMin(min.y,p.y); \
1767  max.y = MagickMax(max.y,p.y); \
1768 }
1769 
1770  switch (method)
1771  {
1772  case AffineDistortion:
1773  { double inverse[6];
1774  InvertAffineCoefficients(coeff, inverse);
1775  s.x = (double) image->page.x;
1776  s.y = (double) image->page.y;
1777  d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1778  d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1779  InitalBounds(d);
1780  s.x = (double) image->page.x+image->columns;
1781  s.y = (double) image->page.y;
1782  d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1783  d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1784  ExpandBounds(d);
1785  s.x = (double) image->page.x;
1786  s.y = (double) image->page.y+image->rows;
1787  d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1788  d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1789  ExpandBounds(d);
1790  s.x = (double) image->page.x+image->columns;
1791  s.y = (double) image->page.y+image->rows;
1792  d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1793  d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1794  ExpandBounds(d);
1795  break;
1796  }
1797  case PerspectiveDistortion:
1798  { double inverse[8], scale;
1799  InvertPerspectiveCoefficients(coeff, inverse);
1800  s.x = (double) image->page.x;
1801  s.y = (double) image->page.y;
1802  scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1803  scale=PerceptibleReciprocal(scale);
1804  d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1805  d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1806  InitalBounds(d);
1807  s.x = (double) image->page.x+image->columns;
1808  s.y = (double) image->page.y;
1809  scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1810  scale=PerceptibleReciprocal(scale);
1811  d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1812  d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1813  ExpandBounds(d);
1814  s.x = (double) image->page.x;
1815  s.y = (double) image->page.y+image->rows;
1816  scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1817  scale=PerceptibleReciprocal(scale);
1818  d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1819  d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1820  ExpandBounds(d);
1821  s.x = (double) image->page.x+image->columns;
1822  s.y = (double) image->page.y+image->rows;
1823  scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1824  scale=PerceptibleReciprocal(scale);
1825  d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1826  d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1827  ExpandBounds(d);
1828  break;
1829  }
1830  case ArcDistortion:
1831  { double a, ca, sa;
1832  /* Forward Map Corners */
1833  a = coeff[0]-coeff[1]/2; ca = cos(a); sa = sin(a);
1834  d.x = coeff[2]*ca;
1835  d.y = coeff[2]*sa;
1836  InitalBounds(d);
1837  d.x = (coeff[2]-coeff[3])*ca;
1838  d.y = (coeff[2]-coeff[3])*sa;
1839  ExpandBounds(d);
1840  a = coeff[0]+coeff[1]/2; ca = cos(a); sa = sin(a);
1841  d.x = coeff[2]*ca;
1842  d.y = coeff[2]*sa;
1843  ExpandBounds(d);
1844  d.x = (coeff[2]-coeff[3])*ca;
1845  d.y = (coeff[2]-coeff[3])*sa;
1846  ExpandBounds(d);
1847  /* Orthogonal points along top of arc */
1848  for( a=(double) (ceil((double) ((coeff[0]-coeff[1]/2.0)/MagickPI2))*MagickPI2);
1849  a<(coeff[0]+coeff[1]/2.0); a+=MagickPI2 ) {
1850  ca = cos(a); sa = sin(a);
1851  d.x = coeff[2]*ca;
1852  d.y = coeff[2]*sa;
1853  ExpandBounds(d);
1854  }
1855  /*
1856  Convert the angle_to_width and radius_to_height
1857  to appropriate scaling factors, to allow faster processing
1858  in the mapping function.
1859  */
1860  coeff[1] = (double) (Magick2PI*image->columns/coeff[1]);
1861  coeff[3] = (double)image->rows/coeff[3];
1862  break;
1863  }
1864  case PolarDistortion:
1865  {
1866  if (number_arguments < 2)
1867  coeff[2] = coeff[3] = 0.0;
1868  min.x = coeff[2]-coeff[0];
1869  max.x = coeff[2]+coeff[0];
1870  min.y = coeff[3]-coeff[0];
1871  max.y = coeff[3]+coeff[0];
1872  /* should be about 1.0 if Rmin = 0 */
1873  coeff[7]=(double) geometry.height/(coeff[0]-coeff[1]);
1874  break;
1875  }
1876  case DePolarDistortion:
1877  {
1878  /* direct calculation as it needs to tile correctly
1879  * for reversibility in a DePolar-Polar cycle */
1880  fix_bounds = MagickFalse;
1881  geometry.x = geometry.y = 0;
1882  geometry.height = (size_t) ceil(coeff[0]-coeff[1]);
1883  geometry.width = (size_t)
1884  ceil((coeff[0]-coeff[1])*(coeff[5]-coeff[4])*0.5);
1885  /* correct scaling factors relative to new size */
1886  coeff[6]=(coeff[5]-coeff[4])/geometry.width; /* changed width */
1887  coeff[7]=(coeff[0]-coeff[1])/geometry.height; /* should be about 1.0 */
1888  break;
1889  }
1891  {
1892  /* direct calculation so center of distortion is either a pixel
1893  * center, or pixel edge. This allows for reversibility of the
1894  * distortion */
1895  geometry.x = geometry.y = 0;
1896  geometry.width = (size_t) ceil( 2.0*coeff[1]*tan(coeff[0]/2.0) );
1897  geometry.height = (size_t) ceil( 2.0*coeff[3]/cos(coeff[0]/2.0) );
1898  /* correct center of distortion relative to new size */
1899  coeff[4] = (double) geometry.width/2.0;
1900  coeff[5] = (double) geometry.height/2.0;
1901  fix_bounds = MagickFalse;
1902  break;
1903  }
1905  {
1906  /* direct calculation center is either pixel center, or pixel edge
1907  * so as to allow reversibility of the image distortion */
1908  geometry.x = geometry.y = 0;
1909  geometry.width = (size_t) ceil(coeff[0]*coeff[1]); /* FOV * radius */
1910  geometry.height = (size_t) (2*coeff[3]); /* input image height */
1911  /* correct center of distortion relative to new size */
1912  coeff[4] = (double) geometry.width/2.0;
1913  coeff[5] = (double) geometry.height/2.0;
1914  fix_bounds = MagickFalse;
1915  break;
1916  }
1917  case ShepardsDistortion:
1920 #if 0
1921  case QuadrilateralDistortion:
1922 #endif
1923  case PolynomialDistortion:
1924  case BarrelDistortion:
1926  default:
1927  /* no calculated bestfit available for these distortions */
1928  bestfit = MagickFalse;
1929  fix_bounds = MagickFalse;
1930  break;
1931  }
1932 
1933  /* Set the output image geometry to calculated 'bestfit'.
1934  Yes this tends to 'over do' the file image size, ON PURPOSE!
1935  Do not do this for DePolar which needs to be exact for virtual tiling.
1936  */
1937  if ( fix_bounds ) {
1938  geometry.x = (ssize_t) floor(min.x-0.5);
1939  geometry.y = (ssize_t) floor(min.y-0.5);
1940  geometry.width=(size_t) ceil(max.x-geometry.x+0.5);
1941  geometry.height=(size_t) ceil(max.y-geometry.y+0.5);
1942  }
1943 
1944  } /* end bestfit destination image calculations */
1945 
1946  /* The user provided a 'viewport' expert option which may
1947  overrides some parts of the current output image geometry.
1948  This also overrides its default 'bestfit' setting.
1949  */
1950  { const char *artifact=GetImageArtifact(image,"distort:viewport");
1951  viewport_given = MagickFalse;
1952  if ( artifact != (const char *) NULL ) {
1953  MagickStatusType flags=ParseAbsoluteGeometry(artifact,&geometry);
1954  if (flags==NoValue)
1955  (void) ThrowMagickException(exception,GetMagickModule(),
1956  OptionWarning,"InvalidSetting","'%s' '%s'",
1957  "distort:viewport",artifact);
1958  else
1959  viewport_given = MagickTrue;
1960  }
1961  }
1962 
1963  /* Verbose output */
1964  if (IsStringTrue(GetImageArtifact(image,"verbose")) != MagickFalse) {
1965  register ssize_t
1966  i;
1967  char image_gen[MagickPathExtent];
1968  const char *lookup;
1969 
1970  /* Set destination image size and virtual offset */
1971  if ( bestfit || viewport_given ) {
1972  (void) FormatLocaleString(image_gen, MagickPathExtent," -size %.20gx%.20g "
1973  "-page %+.20g%+.20g xc: +insert \\\n",(double) geometry.width,
1974  (double) geometry.height,(double) geometry.x,(double) geometry.y);
1975  lookup="v.p{ xx-v.page.x-.5, yy-v.page.y-.5 }";
1976  }
1977  else {
1978  image_gen[0] = '\0'; /* no destination to generate */
1979  lookup = "p{ xx-page.x-.5, yy-page.y-.5 }"; /* simplify lookup */
1980  }
1981 
1982  switch (method) {
1983  case AffineDistortion:
1984  {
1985  double *inverse;
1986 
1987  inverse = (double *) AcquireQuantumMemory(6,sizeof(*inverse));
1988  if (inverse == (double *) NULL) {
1989  coeff = (double *) RelinquishMagickMemory(coeff);
1990  (void) ThrowMagickException(exception,GetMagickModule(),
1991  ResourceLimitError,"MemoryAllocationFailed",
1992  "%s", "DistortImages");
1993  return((Image *) NULL);
1994  }
1995  InvertAffineCoefficients(coeff, inverse);
1996  CoefficientsToAffineArgs(inverse);
1997  (void) FormatLocaleFile(stderr, "Affine Projection:\n");
1998  (void) FormatLocaleFile(stderr, " -distort AffineProjection \\\n '");
1999  for (i=0; i < 5; i++)
2000  (void) FormatLocaleFile(stderr, "%lf,", inverse[i]);
2001  (void) FormatLocaleFile(stderr, "%lf'\n", inverse[5]);
2002  inverse = (double *) RelinquishMagickMemory(inverse);
2003 
2004  (void) FormatLocaleFile(stderr, "Affine Distort, FX Equivelent:\n");
2005  (void) FormatLocaleFile(stderr, "%s", image_gen);
2006  (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2007  (void) FormatLocaleFile(stderr, " xx=%+lf*ii %+lf*jj %+lf;\n",
2008  coeff[0], coeff[1], coeff[2]);
2009  (void) FormatLocaleFile(stderr, " yy=%+lf*ii %+lf*jj %+lf;\n",
2010  coeff[3], coeff[4], coeff[5]);
2011  (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2012 
2013  break;
2014  }
2015 
2016  case PerspectiveDistortion:
2017  {
2018  double *inverse;
2019 
2020  inverse = (double *) AcquireQuantumMemory(8,sizeof(*inverse));
2021  if (inverse == (double *) NULL) {
2022  coeff = (double *) RelinquishMagickMemory(coeff);
2023  (void) ThrowMagickException(exception,GetMagickModule(),
2024  ResourceLimitError,"MemoryAllocationFailed",
2025  "%s", "DistortCoefficients");
2026  return((Image *) NULL);
2027  }
2028  InvertPerspectiveCoefficients(coeff, inverse);
2029  (void) FormatLocaleFile(stderr, "Perspective Projection:\n");
2030  (void) FormatLocaleFile(stderr, " -distort PerspectiveProjection \\\n '");
2031  for (i=0; i<4; i++)
2032  (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
2033  (void) FormatLocaleFile(stderr, "\n ");
2034  for (; i<7; i++)
2035  (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
2036  (void) FormatLocaleFile(stderr, "%lf'\n", inverse[7]);
2037  inverse = (double *) RelinquishMagickMemory(inverse);
2038 
2039  (void) FormatLocaleFile(stderr, "Perspective Distort, FX Equivelent:\n");
2040  (void) FormatLocaleFile(stderr, "%s", image_gen);
2041  (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2042  (void) FormatLocaleFile(stderr, " rr=%+lf*ii %+lf*jj + 1;\n",
2043  coeff[6], coeff[7]);
2044  (void) FormatLocaleFile(stderr, " xx=(%+lf*ii %+lf*jj %+lf)/rr;\n",
2045  coeff[0], coeff[1], coeff[2]);
2046  (void) FormatLocaleFile(stderr, " yy=(%+lf*ii %+lf*jj %+lf)/rr;\n",
2047  coeff[3], coeff[4], coeff[5]);
2048  (void) FormatLocaleFile(stderr, " rr%s0 ? %s : blue' \\\n",
2049  coeff[8] < 0 ? "<" : ">", lookup);
2050  break;
2051  }
2052 
2054  (void) FormatLocaleFile(stderr, "BilinearForward Mapping Equations:\n");
2055  (void) FormatLocaleFile(stderr, "%s", image_gen);
2056  (void) FormatLocaleFile(stderr, " i = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
2057  coeff[0], coeff[1], coeff[2], coeff[3]);
2058  (void) FormatLocaleFile(stderr, " j = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
2059  coeff[4], coeff[5], coeff[6], coeff[7]);
2060 #if 0
2061  /* for debugging */
2062  (void) FormatLocaleFile(stderr, " c8 = %+lf c9 = 2*a = %+lf;\n",
2063  coeff[8], coeff[9]);
2064 #endif
2065  (void) FormatLocaleFile(stderr, "BilinearForward Distort, FX Equivelent:\n");
2066  (void) FormatLocaleFile(stderr, "%s", image_gen);
2067  (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
2068  0.5-coeff[3], 0.5-coeff[7]);
2069  (void) FormatLocaleFile(stderr, " bb=%lf*ii %+lf*jj %+lf;\n",
2070  coeff[6], -coeff[2], coeff[8]);
2071  /* Handle Special degenerate (non-quadratic) or trapezoidal case */
2072  if ( coeff[9] != 0 ) {
2073  (void) FormatLocaleFile(stderr, " rt=bb*bb %+lf*(%lf*ii%+lf*jj);\n",
2074  -2*coeff[9], coeff[4], -coeff[0]);
2075  (void) FormatLocaleFile(stderr, " yy=( -bb + sqrt(rt) ) / %lf;\n",
2076  coeff[9]);
2077  } else
2078  (void) FormatLocaleFile(stderr, " yy=(%lf*ii%+lf*jj)/bb;\n",
2079  -coeff[4], coeff[0]);
2080  (void) FormatLocaleFile(stderr, " xx=(ii %+lf*yy)/(%lf %+lf*yy);\n",
2081  -coeff[1], coeff[0], coeff[2]);
2082  if ( coeff[9] != 0 )
2083  (void) FormatLocaleFile(stderr, " (rt < 0 ) ? red : %s'\n", lookup);
2084  else
2085  (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2086  break;
2087 
2089 #if 0
2090  (void) FormatLocaleFile(stderr, "Polynomial Projection Distort:\n");
2091  (void) FormatLocaleFile(stderr, " -distort PolynomialProjection \\\n");
2092  (void) FormatLocaleFile(stderr, " '1.5, %lf, %lf, %lf, %lf,\n",
2093  coeff[3], coeff[0], coeff[1], coeff[2]);
2094  (void) FormatLocaleFile(stderr, " %lf, %lf, %lf, %lf'\n",
2095  coeff[7], coeff[4], coeff[5], coeff[6]);
2096 #endif
2097  (void) FormatLocaleFile(stderr, "BilinearReverse Distort, FX Equivelent:\n");
2098  (void) FormatLocaleFile(stderr, "%s", image_gen);
2099  (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2100  (void) FormatLocaleFile(stderr, " xx=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
2101  coeff[0], coeff[1], coeff[2], coeff[3]);
2102  (void) FormatLocaleFile(stderr, " yy=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
2103  coeff[4], coeff[5], coeff[6], coeff[7]);
2104  (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2105  break;
2106 
2107  case PolynomialDistortion:
2108  {
2109  size_t nterms = (size_t) coeff[1];
2110  (void) FormatLocaleFile(stderr, "Polynomial (order %lg, terms %lu), FX Equivelent\n",
2111  coeff[0],(unsigned long) nterms);
2112  (void) FormatLocaleFile(stderr, "%s", image_gen);
2113  (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2114  (void) FormatLocaleFile(stderr, " xx =");
2115  for (i=0; i<(ssize_t) nterms; i++) {
2116  if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n ");
2117  (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i],
2118  poly_basis_str(i));
2119  }
2120  (void) FormatLocaleFile(stderr, ";\n yy =");
2121  for (i=0; i<(ssize_t) nterms; i++) {
2122  if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n ");
2123  (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i+nterms],
2124  poly_basis_str(i));
2125  }
2126  (void) FormatLocaleFile(stderr, ";\n %s' \\\n", lookup);
2127  break;
2128  }
2129  case ArcDistortion:
2130  {
2131  (void) FormatLocaleFile(stderr, "Arc Distort, Internal Coefficients:\n");
2132  for ( i=0; i<5; i++ )
2133  (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
2134  (void) FormatLocaleFile(stderr, "Arc Distort, FX Equivelent:\n");
2135  (void) FormatLocaleFile(stderr, "%s", image_gen);
2136  (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x; jj=j+page.y;\n");
2137  (void) FormatLocaleFile(stderr, " xx=(atan2(jj,ii)%+lf)/(2*pi);\n",
2138  -coeff[0]);
2139  (void) FormatLocaleFile(stderr, " xx=xx-round(xx);\n");
2140  (void) FormatLocaleFile(stderr, " xx=xx*%lf %+lf;\n",
2141  coeff[1], coeff[4]);
2142  (void) FormatLocaleFile(stderr, " yy=(%lf - hypot(ii,jj)) * %lf;\n",
2143  coeff[2], coeff[3]);
2144  (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n");
2145  break;
2146  }
2147  case PolarDistortion:
2148  {
2149  (void) FormatLocaleFile(stderr, "Polar Distort, Internal Coefficents\n");
2150  for ( i=0; i<8; i++ )
2151  (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
2152  (void) FormatLocaleFile(stderr, "Polar Distort, FX Equivelent:\n");
2153  (void) FormatLocaleFile(stderr, "%s", image_gen);
2154  (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
2155  -coeff[2], -coeff[3]);
2156  (void) FormatLocaleFile(stderr, " xx=(atan2(ii,jj)%+lf)/(2*pi);\n",
2157  -(coeff[4]+coeff[5])/2 );
2158  (void) FormatLocaleFile(stderr, " xx=xx-round(xx);\n");
2159  (void) FormatLocaleFile(stderr, " xx=xx*2*pi*%lf + v.w/2;\n",
2160  coeff[6] );
2161  (void) FormatLocaleFile(stderr, " yy=(hypot(ii,jj)%+lf)*%lf;\n",
2162  -coeff[1], coeff[7] );
2163  (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n");
2164  break;
2165  }
2166  case DePolarDistortion:
2167  {
2168  (void) FormatLocaleFile(stderr, "DePolar Distort, Internal Coefficents\n");
2169  for ( i=0; i<8; i++ )
2170  (void) FormatLocaleFile(stderr, " c%.20g = %+lf\n", (double) i, coeff[i]);
2171  (void) FormatLocaleFile(stderr, "DePolar Distort, FX Equivelent:\n");
2172  (void) FormatLocaleFile(stderr, "%s", image_gen);
2173  (void) FormatLocaleFile(stderr, " -fx 'aa=(i+.5)*%lf %+lf;\n", coeff[6], +coeff[4] );
2174  (void) FormatLocaleFile(stderr, " rr=(j+.5)*%lf %+lf;\n", coeff[7], +coeff[1] );
2175  (void) FormatLocaleFile(stderr, " xx=rr*sin(aa) %+lf;\n", coeff[2] );
2176  (void) FormatLocaleFile(stderr, " yy=rr*cos(aa) %+lf;\n", coeff[3] );
2177  (void) FormatLocaleFile(stderr, " v.p{xx-.5,yy-.5}' \\\n");
2178  break;
2179  }
2181  {
2182  (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, Internal Coefficents\n");
2183  (void) FormatLocaleFile(stderr, " cylinder_radius = %+lf\n", coeff[1]);
2184  (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, FX Equivelent:\n");
2185  (void) FormatLocaleFile(stderr, "%s", image_gen);
2186  (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
2187  -coeff[4], -coeff[5]);
2188  (void) FormatLocaleFile(stderr, " aa=atan(ii/%+lf);\n", coeff[1] );
2189  (void) FormatLocaleFile(stderr, " xx=%lf*aa%+lf;\n",
2190  coeff[1], coeff[2] );
2191  (void) FormatLocaleFile(stderr, " yy=jj*cos(aa)%+lf;\n", coeff[3] );
2192  (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2193  break;
2194  }
2196  {
2197  (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, Internal Coefficents\n");
2198  (void) FormatLocaleFile(stderr, " cylinder_radius = %+lf\n", coeff[1]);
2199  (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, FX Equivelent:\n");
2200  (void) FormatLocaleFile(stderr, "%s", image_gen);
2201  (void) FormatLocaleFile(stderr, " -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
2202  -coeff[4], -coeff[5]);
2203  (void) FormatLocaleFile(stderr, " ii=ii/%+lf;\n", coeff[1] );
2204  (void) FormatLocaleFile(stderr, " xx=%lf*tan(ii)%+lf;\n",
2205  coeff[1], coeff[2] );
2206  (void) FormatLocaleFile(stderr, " yy=jj/cos(ii)%+lf;\n",
2207  coeff[3] );
2208  (void) FormatLocaleFile(stderr, " %s' \\\n", lookup);
2209  break;
2210  }
2211  case BarrelDistortion:
2213  { double xc,yc;
2214  /* NOTE: This does the barrel roll in pixel coords not image coords
2215  ** The internal distortion must do it in image coordinates,
2216  ** so that is what the center coeff (8,9) is given in.
2217  */
2218  xc = ((double)image->columns-1.0)/2.0 + image->page.x;
2219  yc = ((double)image->rows-1.0)/2.0 + image->page.y;
2220  (void) FormatLocaleFile(stderr, "Barrel%s Distort, FX Equivelent:\n",
2221  method == BarrelDistortion ? "" : "Inv");
2222  (void) FormatLocaleFile(stderr, "%s", image_gen);
2223  if ( fabs(coeff[8]-xc-0.5) < 0.1 && fabs(coeff[9]-yc-0.5) < 0.1 )
2224  (void) FormatLocaleFile(stderr, " -fx 'xc=(w-1)/2; yc=(h-1)/2;\n");
2225  else
2226  (void) FormatLocaleFile(stderr, " -fx 'xc=%lf; yc=%lf;\n",
2227  coeff[8]-0.5, coeff[9]-0.5);
2228  (void) FormatLocaleFile(stderr,
2229  " ii=i-xc; jj=j-yc; rr=hypot(ii,jj);\n");
2230  (void) FormatLocaleFile(stderr, " ii=ii%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2231  method == BarrelDistortion ? "*" : "/",
2232  coeff[0],coeff[1],coeff[2],coeff[3]);
2233  (void) FormatLocaleFile(stderr, " jj=jj%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2234  method == BarrelDistortion ? "*" : "/",
2235  coeff[4],coeff[5],coeff[6],coeff[7]);
2236  (void) FormatLocaleFile(stderr, " v.p{fx*ii+xc,fy*jj+yc}' \\\n");
2237  }
2238  default:
2239  break;
2240  }
2241  }
2242 
2243  /* The user provided a 'scale' expert option will scale the
2244  output image size, by the factor given allowing for super-sampling
2245  of the distorted image space. Any scaling factors must naturally
2246  be halved as a result.
2247  */
2248  { const char *artifact;
2249  artifact=GetImageArtifact(image,"distort:scale");
2250  output_scaling = 1.0;
2251  if (artifact != (const char *) NULL) {
2252  output_scaling = fabs(StringToDouble(artifact,(char **) NULL));
2253  geometry.width=(size_t) (output_scaling*geometry.width+0.5);
2254  geometry.height=(size_t) (output_scaling*geometry.height+0.5);
2255  geometry.x=(ssize_t) (output_scaling*geometry.x+0.5);
2256  geometry.y=(ssize_t) (output_scaling*geometry.y+0.5);
2257  if ( output_scaling < 0.1 ) {
2258  coeff = (double *) RelinquishMagickMemory(coeff);
2260  "InvalidArgument","%s", "-set option:distort:scale" );
2261  return((Image *) NULL);
2262  }
2263  output_scaling = 1/output_scaling;
2264  }
2265  }
2266 #define ScaleFilter(F,A,B,C,D) \
2267  ScaleResampleFilter( (F), \
2268  output_scaling*(A), output_scaling*(B), \
2269  output_scaling*(C), output_scaling*(D) )
2270 
2271  /*
2272  Initialize the distort image attributes.
2273  */
2274  distort_image=CloneImage(image,geometry.width,geometry.height,MagickTrue,
2275  exception);
2276  if (distort_image == (Image *) NULL)
2277  {
2278  coeff=(double *) RelinquishMagickMemory(coeff);
2279  return((Image *) NULL);
2280  }
2281  /* if image is ColorMapped - change it to DirectClass */
2282  if (SetImageStorageClass(distort_image,DirectClass,exception) == MagickFalse)
2283  {
2284  coeff=(double *) RelinquishMagickMemory(coeff);
2285  distort_image=DestroyImage(distort_image);
2286  return((Image *) NULL);
2287  }
2288  if ((IsPixelInfoGray(&distort_image->background_color) == MagickFalse) &&
2289  (IsGrayColorspace(distort_image->colorspace) != MagickFalse))
2290  (void) SetImageColorspace(distort_image,sRGBColorspace,exception);
2291  if (distort_image->background_color.alpha_trait != UndefinedPixelTrait)
2292  distort_image->alpha_trait=BlendPixelTrait;
2293  distort_image->page.x=geometry.x;
2294  distort_image->page.y=geometry.y;
2295 
2296  { /* ----- MAIN CODE -----
2297  Sample the source image to each pixel in the distort image.
2298  */
2299  CacheView
2300  *distort_view;
2301 
2303  status;
2304 
2306  progress;
2307 
2308  PixelInfo
2309  zero;
2310 
2312  **magick_restrict resample_filter;
2313 
2314  ssize_t
2315  j;
2316 
2317  status=MagickTrue;
2318  progress=0;
2319  GetPixelInfo(distort_image,&zero);
2320  resample_filter=AcquireResampleFilterThreadSet(image,
2322  distort_view=AcquireAuthenticCacheView(distort_image,exception);
2323 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2324  #pragma omp parallel for schedule(static) shared(progress,status) \
2325  magick_number_threads(image,distort_image,distort_image->rows,1)
2326 #endif
2327  for (j=0; j < (ssize_t) distort_image->rows; j++)
2328  {
2329  const int
2330  id = GetOpenMPThreadId();
2331 
2332  double
2333  validity; /* how mathematically valid is this the mapping */
2334 
2336  sync;
2337 
2338  PixelInfo
2339  pixel, /* pixel color to assign to distorted image */
2340  invalid; /* the color to assign when distort result is invalid */
2341 
2342  PointInfo
2343  d,
2344  s; /* transform destination image x,y to source image x,y */
2345 
2346  register ssize_t
2347  i;
2348 
2349  register Quantum
2350  *magick_restrict q;
2351 
2352  q=QueueCacheViewAuthenticPixels(distort_view,0,j,distort_image->columns,1,
2353  exception);
2354  if (q == (Quantum *) NULL)
2355  {
2356  status=MagickFalse;
2357  continue;
2358  }
2359  pixel=zero;
2360 
2361  /* Define constant scaling vectors for Affine Distortions
2362  Other methods are either variable, or use interpolated lookup
2363  */
2364  switch (method)
2365  {
2366  case AffineDistortion:
2367  ScaleFilter( resample_filter[id],
2368  coeff[0], coeff[1],
2369  coeff[3], coeff[4] );
2370  break;
2371  default:
2372  break;
2373  }
2374 
2375  /* Initialize default pixel validity
2376  * negative: pixel is invalid output 'matte_color'
2377  * 0.0 to 1.0: antialiased, mix with resample output
2378  * 1.0 or greater: use resampled output.
2379  */
2380  validity = 1.0;
2381 
2382  ConformPixelInfo(distort_image,&distort_image->matte_color,&invalid,
2383  exception);
2384  for (i=0; i < (ssize_t) distort_image->columns; i++)
2385  {
2386  /* map pixel coordinate to distortion space coordinate */
2387  d.x = (double) (geometry.x+i+0.5)*output_scaling;
2388  d.y = (double) (geometry.y+j+0.5)*output_scaling;
2389  s = d; /* default is a no-op mapping */
2390  switch (method)
2391  {
2392  case AffineDistortion:
2393  {
2394  s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2395  s.y=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2396  /* Affine partial derivitives are constant -- set above */
2397  break;
2398  }
2399  case PerspectiveDistortion:
2400  {
2401  double
2402  p,q,r,abs_r,abs_c6,abs_c7,scale;
2403  /* perspective is a ratio of affines */
2404  p=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2405  q=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2406  r=coeff[6]*d.x+coeff[7]*d.y+1.0;
2407  /* Pixel Validity -- is it a 'sky' or 'ground' pixel */
2408  validity = (r*coeff[8] < 0.0) ? 0.0 : 1.0;
2409  /* Determine horizon anti-alias blending */
2410  abs_r = fabs(r)*2;
2411  abs_c6 = fabs(coeff[6]);
2412  abs_c7 = fabs(coeff[7]);
2413  if ( abs_c6 > abs_c7 ) {
2414  if ( abs_r < abs_c6*output_scaling )
2415  validity = 0.5 - coeff[8]*r/(coeff[6]*output_scaling);
2416  }
2417  else if ( abs_r < abs_c7*output_scaling )
2418  validity = 0.5 - coeff[8]*r/(coeff[7]*output_scaling);
2419  /* Perspective Sampling Point (if valid) */
2420  if ( validity > 0.0 ) {
2421  /* divide by r affine, for perspective scaling */
2422  scale = 1.0/r;
2423  s.x = p*scale;
2424  s.y = q*scale;
2425  /* Perspective Partial Derivatives or Scaling Vectors */
2426  scale *= scale;
2427  ScaleFilter( resample_filter[id],
2428  (r*coeff[0] - p*coeff[6])*scale,
2429  (r*coeff[1] - p*coeff[7])*scale,
2430  (r*coeff[3] - q*coeff[6])*scale,
2431  (r*coeff[4] - q*coeff[7])*scale );
2432  }
2433  break;
2434  }
2436  {
2437  /* Reversed Mapped is just a simple polynomial */
2438  s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2]*d.x*d.y+coeff[3];
2439  s.y=coeff[4]*d.x+coeff[5]*d.y
2440  +coeff[6]*d.x*d.y+coeff[7];
2441  /* Bilinear partial derivitives of scaling vectors */
2442  ScaleFilter( resample_filter[id],
2443  coeff[0] + coeff[2]*d.y,
2444  coeff[1] + coeff[2]*d.x,
2445  coeff[4] + coeff[6]*d.y,
2446  coeff[5] + coeff[6]*d.x );
2447  break;
2448  }
2450  {
2451  /* Forward mapped needs reversed polynomial equations
2452  * which unfortunatally requires a square root! */
2453  double b,c;
2454  d.x -= coeff[3]; d.y -= coeff[7];
2455  b = coeff[6]*d.x - coeff[2]*d.y + coeff[8];
2456  c = coeff[4]*d.x - coeff[0]*d.y;
2457 
2458  validity = 1.0;
2459  /* Handle Special degenerate (non-quadratic) case
2460  * Currently without horizon anti-alising */
2461  if ( fabs(coeff[9]) < MagickEpsilon )
2462  s.y = -c/b;
2463  else {
2464  c = b*b - 2*coeff[9]*c;
2465  if ( c < 0.0 )
2466  validity = 0.0;
2467  else
2468  s.y = ( -b + sqrt(c) )/coeff[9];
2469  }
2470  if ( validity > 0.0 )
2471  s.x = ( d.x - coeff[1]*s.y) / ( coeff[0] + coeff[2]*s.y );
2472 
2473  /* NOTE: the sign of the square root should be -ve for parts
2474  where the source image becomes 'flipped' or 'mirrored'.
2475  FUTURE: Horizon handling
2476  FUTURE: Scaling factors or Deritives (how?)
2477  */
2478  break;
2479  }
2480 #if 0
2481  case BilinearDistortion:
2482  /* Bilinear mapping of any Quadrilateral to any Quadrilateral */
2483  /* UNDER DEVELOPMENT */
2484  break;
2485 #endif
2486  case PolynomialDistortion:
2487  {
2488  /* multi-ordered polynomial */
2489  register ssize_t
2490  k;
2491 
2492  ssize_t
2493  nterms=(ssize_t)coeff[1];
2494 
2495  PointInfo
2496  du,dv; /* the du,dv vectors from unit dx,dy -- derivatives */
2497 
2498  s.x=s.y=du.x=du.y=dv.x=dv.y=0.0;
2499  for(k=0; k < nterms; k++) {
2500  s.x += poly_basis_fn(k,d.x,d.y)*coeff[2+k];
2501  du.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k];
2502  du.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k];
2503  s.y += poly_basis_fn(k,d.x,d.y)*coeff[2+k+nterms];
2504  dv.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k+nterms];
2505  dv.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k+nterms];
2506  }
2507  ScaleFilter( resample_filter[id], du.x,du.y,dv.x,dv.y );
2508  break;
2509  }
2510  case ArcDistortion:
2511  {
2512  /* what is the angle and radius in the destination image */
2513  s.x = (double) ((atan2(d.y,d.x) - coeff[0])/Magick2PI);
2514  s.x -= MagickRound(s.x); /* angle */
2515  s.y = hypot(d.x,d.y); /* radius */
2516 
2517  /* Arc Distortion Partial Scaling Vectors
2518  Are derived by mapping the perpendicular unit vectors
2519  dR and dA*R*2PI rather than trying to map dx and dy
2520  The results is a very simple orthogonal aligned ellipse.
2521  */
2522  if ( s.y > MagickEpsilon )
2523  ScaleFilter( resample_filter[id],
2524  (double) (coeff[1]/(Magick2PI*s.y)), 0, 0, coeff[3] );
2525  else
2526  ScaleFilter( resample_filter[id],
2527  distort_image->columns*2, 0, 0, coeff[3] );
2528 
2529  /* now scale the angle and radius for source image lookup point */
2530  s.x = s.x*coeff[1] + coeff[4] + image->page.x +0.5;
2531  s.y = (coeff[2] - s.y) * coeff[3] + image->page.y;
2532  break;
2533  }
2534  case PolarDistortion:
2535  { /* 2D Cartesain to Polar View */
2536  d.x -= coeff[2];
2537  d.y -= coeff[3];
2538  s.x = atan2(d.x,d.y) - (coeff[4]+coeff[5])/2;
2539  s.x /= Magick2PI;
2540  s.x -= MagickRound(s.x);
2541  s.x *= Magick2PI; /* angle - relative to centerline */
2542  s.y = hypot(d.x,d.y); /* radius */
2543 
2544  /* Polar Scaling vectors are based on mapping dR and dA vectors
2545  This results in very simple orthogonal scaling vectors
2546  */
2547  if ( s.y > MagickEpsilon )
2548  ScaleFilter( resample_filter[id],
2549  (double) (coeff[6]/(Magick2PI*s.y)), 0, 0, coeff[7] );
2550  else
2551  ScaleFilter( resample_filter[id],
2552  distort_image->columns*2, 0, 0, coeff[7] );
2553 
2554  /* now finish mapping radius/angle to source x,y coords */
2555  s.x = s.x*coeff[6] + (double)image->columns/2.0 + image->page.x;
2556  s.y = (s.y-coeff[1])*coeff[7] + image->page.y;
2557  break;
2558  }
2559  case DePolarDistortion:
2560  { /* @D Polar to Carteasain */
2561  /* ignore all destination virtual offsets */
2562  d.x = ((double)i+0.5)*output_scaling*coeff[6]+coeff[4];
2563  d.y = ((double)j+0.5)*output_scaling*coeff[7]+coeff[1];
2564  s.x = d.y*sin(d.x) + coeff[2];
2565  s.y = d.y*cos(d.x) + coeff[3];
2566  /* derivatives are usless - better to use SuperSampling */
2567  break;
2568  }
2570  { /* 3D Cylinder to Tangential Plane */
2571  double ax, cx;
2572  /* relative to center of distortion */
2573  d.x -= coeff[4]; d.y -= coeff[5];
2574  d.x /= coeff[1]; /* x' = x/r */
2575  ax=atan(d.x); /* aa = atan(x/r) = u/r */
2576  cx=cos(ax); /* cx = cos(atan(x/r)) = 1/sqrt(x^2+u^2) */
2577  s.x = coeff[1]*ax; /* u = r*atan(x/r) */
2578  s.y = d.y*cx; /* v = y*cos(u/r) */
2579  /* derivatives... (see personnal notes) */
2580  ScaleFilter( resample_filter[id],
2581  1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2582 #if 0
2583 if ( i == 0 && j == 0 ) {
2584  fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2585  fprintf(stderr, "phi = %lf\n", (double)(ax * 180.0/MagickPI) );
2586  fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n",
2587  1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2588  fflush(stderr); }
2589 #endif
2590  /* add center of distortion in source */
2591  s.x += coeff[2]; s.y += coeff[3];
2592  break;
2593  }
2595  { /* 3D Cylinder to Tangential Plane */
2596  /* relative to center of distortion */
2597  d.x -= coeff[4]; d.y -= coeff[5];
2598 
2599  /* is pixel valid - horizon of a infinite Virtual-Pixel Plane
2600  * (see Anthony Thyssen's personal note) */
2601  validity = (double) (coeff[1]*MagickPI2 - fabs(d.x))/output_scaling + 0.5;
2602 
2603  if ( validity > 0.0 ) {
2604  double cx,tx;
2605  d.x /= coeff[1]; /* x'= x/r */
2606  cx = 1/cos(d.x); /* cx = 1/cos(x/r) */
2607  tx = tan(d.x); /* tx = tan(x/r) */
2608  s.x = coeff[1]*tx; /* u = r * tan(x/r) */
2609  s.y = d.y*cx; /* v = y / cos(x/r) */
2610  /* derivatives... (see Anthony Thyssen's personal notes) */
2611  ScaleFilter( resample_filter[id],
2612  cx*cx, 0.0, s.y*cx/coeff[1], cx );
2613 #if 0
2614 /*if ( i == 0 && j == 0 )*/
2615 if ( d.x == 0.5 && d.y == 0.5 ) {
2616  fprintf(stderr, "x=%lf y=%lf u=%lf v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2617  fprintf(stderr, "radius = %lf phi = %lf validity = %lf\n",
2618  coeff[1], (double)(d.x * 180.0/MagickPI), validity );
2619  fprintf(stderr, "du/dx=%lf du/dx=%lf dv/dx=%lf dv/dy=%lf\n",
2620  cx*cx, 0.0, s.y*cx/coeff[1], cx);
2621  fflush(stderr); }
2622 #endif
2623  }
2624  /* add center of distortion in source */
2625  s.x += coeff[2]; s.y += coeff[3];
2626  break;
2627  }
2628  case BarrelDistortion:
2630  { /* Lens Barrel Distionion Correction */
2631  double r,fx,fy,gx,gy;
2632  /* Radial Polynomial Distortion (de-normalized) */
2633  d.x -= coeff[8];
2634  d.y -= coeff[9];
2635  r = sqrt(d.x*d.x+d.y*d.y);
2636  if ( r > MagickEpsilon ) {
2637  fx = ((coeff[0]*r + coeff[1])*r + coeff[2])*r + coeff[3];
2638  fy = ((coeff[4]*r + coeff[5])*r + coeff[6])*r + coeff[7];
2639  gx = ((3*coeff[0]*r + 2*coeff[1])*r + coeff[2])/r;
2640  gy = ((3*coeff[4]*r + 2*coeff[5])*r + coeff[6])/r;
2641  /* adjust functions and scaling for 'inverse' form */
2642  if ( method == BarrelInverseDistortion ) {
2643  fx = 1/fx; fy = 1/fy;
2644  gx *= -fx*fx; gy *= -fy*fy;
2645  }
2646  /* Set the source pixel to lookup and EWA derivative vectors */
2647  s.x = d.x*fx + coeff[8];
2648  s.y = d.y*fy + coeff[9];
2649  ScaleFilter( resample_filter[id],
2650  gx*d.x*d.x + fx, gx*d.x*d.y,
2651  gy*d.x*d.y, gy*d.y*d.y + fy );
2652  }
2653  else {
2654  /* Special handling to avoid divide by zero when r==0
2655  **
2656  ** The source and destination pixels match in this case
2657  ** which was set at the top of the loop using s = d;
2658  ** otherwise... s.x=coeff[8]; s.y=coeff[9];
2659  */
2660  if ( method == BarrelDistortion )
2661  ScaleFilter( resample_filter[id],
2662  coeff[3], 0, 0, coeff[7] );
2663  else /* method == BarrelInverseDistortion */
2664  /* FUTURE, trap for D==0 causing division by zero */
2665  ScaleFilter( resample_filter[id],
2666  1.0/coeff[3], 0, 0, 1.0/coeff[7] );
2667  }
2668  break;
2669  }
2670  case ShepardsDistortion:
2671  { /* Shepards Method, or Inverse Weighted Distance for
2672  displacement around the destination image control points
2673  The input arguments are the coefficents to the function.
2674  This is more of a 'displacement' function rather than an
2675  absolute distortion function.
2676 
2677  Note: We can not determine derivatives using shepards method
2678  so only a point sample interpolatation can be used.
2679  */
2680  size_t
2681  i;
2682  double
2683  denominator;
2684 
2685  denominator = s.x = s.y = 0;
2686  for(i=0; i<number_arguments; i+=4) {
2687  double weight =
2688  ((double)d.x-arguments[i+2])*((double)d.x-arguments[i+2])
2689  + ((double)d.y-arguments[i+3])*((double)d.y-arguments[i+3]);
2690  weight = pow(weight,coeff[0]); /* shepards power factor */
2691  weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
2692 
2693  s.x += (arguments[ i ]-arguments[i+2])*weight;
2694  s.y += (arguments[i+1]-arguments[i+3])*weight;
2695  denominator += weight;
2696  }
2697  s.x /= denominator;
2698  s.y /= denominator;
2699  s.x += d.x; /* make it as relative displacement */
2700  s.y += d.y;
2701  break;
2702  }
2703  default:
2704  break; /* use the default no-op given above */
2705  }
2706  /* map virtual canvas location back to real image coordinate */
2707  if ( bestfit && method != ArcDistortion ) {
2708  s.x -= image->page.x;
2709  s.y -= image->page.y;
2710  }
2711  s.x -= 0.5;
2712  s.y -= 0.5;
2713 
2714  if ( validity <= 0.0 ) {
2715  /* result of distortion is an invalid pixel - don't resample */
2716  SetPixelViaPixelInfo(distort_image,&invalid,q);
2717  }
2718  else {
2719  /* resample the source image to find its correct color */
2720  (void) ResamplePixelColor(resample_filter[id],s.x,s.y,&pixel,
2721  exception);
2722  /* if validity between 0.0 and 1.0 mix result with invalid pixel */
2723  if ( validity < 1.0 ) {
2724  /* Do a blend of sample color and invalid pixel */
2725  /* should this be a 'Blend', or an 'Over' compose */
2726  CompositePixelInfoBlend(&pixel,validity,&invalid,(1.0-validity),
2727  &pixel);
2728  }
2729  SetPixelViaPixelInfo(distort_image,&pixel,q);
2730  }
2731  q+=GetPixelChannels(distort_image);
2732  }
2733  sync=SyncCacheViewAuthenticPixels(distort_view,exception);
2734  if (sync == MagickFalse)
2735  status=MagickFalse;
2736  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2737  {
2739  proceed;
2740 
2741 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2742  #pragma omp critical (MagickCore_DistortImage)
2743 #endif
2744  proceed=SetImageProgress(image,DistortImageTag,progress++,
2745  image->rows);
2746  if (proceed == MagickFalse)
2747  status=MagickFalse;
2748  }
2749  }
2750  distort_view=DestroyCacheView(distort_view);
2751  resample_filter=DestroyResampleFilterThreadSet(resample_filter);
2752 
2753  if (status == MagickFalse)
2754  distort_image=DestroyImage(distort_image);
2755  }
2756 
2757  /* Arc does not return an offset unless 'bestfit' is in effect
2758  And the user has not provided an overriding 'viewport'.
2759  */
2760  if ( method == ArcDistortion && !bestfit && !viewport_given ) {
2761  distort_image->page.x = 0;
2762  distort_image->page.y = 0;
2763  }
2764  coeff=(double *) RelinquishMagickMemory(coeff);
2765  return(distort_image);
2766 }
2767 
2768 /*
2769 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2770 % %
2771 % %
2772 % %
2773 % R o t a t e I m a g e %
2774 % %
2775 % %
2776 % %
2777 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2778 %
2779 % RotateImage() creates a new image that is a rotated copy of an existing
2780 % one. Positive angles rotate counter-clockwise (right-hand rule), while
2781 % negative angles rotate clockwise. Rotated images are usually larger than
2782 % the originals and have 'empty' triangular corners. X axis. Empty
2783 % triangles left over from shearing the image are filled with the background
2784 % color defined by member 'background_color' of the image. RotateImage
2785 % allocates the memory necessary for the new Image structure and returns a
2786 % pointer to the new image.
2787 %
2788 % The format of the RotateImage method is:
2789 %
2790 % Image *RotateImage(const Image *image,const double degrees,
2791 % ExceptionInfo *exception)
2792 %
2793 % A description of each parameter follows.
2794 %
2795 % o image: the image.
2796 %
2797 % o degrees: Specifies the number of degrees to rotate the image.
2798 %
2799 % o exception: return any errors or warnings in this structure.
2800 %
2801 */
2802 MagickExport Image *RotateImage(const Image *image,const double degrees,
2803  ExceptionInfo *exception)
2804 {
2805  Image
2806  *distort_image,
2807  *rotate_image;
2808 
2809  double
2810  angle;
2811 
2812  PointInfo
2813  shear;
2814 
2815  size_t
2816  rotations;
2817 
2818  /*
2819  Adjust rotation angle.
2820  */
2821  assert(image != (Image *) NULL);
2822  assert(image->signature == MagickCoreSignature);
2823  if (image->debug != MagickFalse)
2824  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2825  assert(exception != (ExceptionInfo *) NULL);
2826  assert(exception->signature == MagickCoreSignature);
2827  angle=fmod(degrees,360.0);
2828  while (angle < -45.0)
2829  angle+=360.0;
2830  for (rotations=0; angle > 45.0; rotations++)
2831  angle-=90.0;
2832  rotations%=4;
2833  shear.x=(-tan((double) DegreesToRadians(angle)/2.0));
2834  shear.y=sin((double) DegreesToRadians(angle));
2835  if ((fabs(shear.x) < MagickEpsilon) && (fabs(shear.y) < MagickEpsilon))
2836  return(IntegralRotateImage(image,rotations,exception));
2837  distort_image=CloneImage(image,0,0,MagickTrue,exception);
2838  if (distort_image == (Image *) NULL)
2839  return((Image *) NULL);
2841  exception);
2842  rotate_image=DistortImage(distort_image,ScaleRotateTranslateDistortion,1,
2843  &degrees,MagickTrue,exception);
2844  distort_image=DestroyImage(distort_image);
2845  return(rotate_image);
2846 }
2847 
2848 /*
2849 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2850 % %
2851 % %
2852 % %
2853 % S p a r s e C o l o r I m a g e %
2854 % %
2855 % %
2856 % %
2857 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2858 %
2859 % SparseColorImage(), given a set of coordinates, interpolates the colors
2860 % found at those coordinates, across the whole image, using various methods.
2861 %
2862 % The format of the SparseColorImage() method is:
2863 %
2864 % Image *SparseColorImage(const Image *image,
2865 % const SparseColorMethod method,const size_t number_arguments,
2866 % const double *arguments,ExceptionInfo *exception)
2867 %
2868 % A description of each parameter follows:
2869 %
2870 % o image: the image to be filled in.
2871 %
2872 % o method: the method to fill in the gradient between the control points.
2873 %
2874 % The methods used for SparseColor() are often simular to methods
2875 % used for DistortImage(), and even share the same code for determination
2876 % of the function coefficents, though with more dimensions (or resulting
2877 % values).
2878 %
2879 % o number_arguments: the number of arguments given.
2880 %
2881 % o arguments: array of floating point arguments for this method--
2882 % x,y,color_values-- with color_values given as normalized values.
2883 %
2884 % o exception: return any errors or warnings in this structure
2885 %
2886 */
2888  const SparseColorMethod method,const size_t number_arguments,
2889  const double *arguments,ExceptionInfo *exception)
2890 {
2891 #define SparseColorTag "Distort/SparseColor"
2892 
2894  sparse_method;
2895 
2896  double
2897  *coeff;
2898 
2899  Image
2900  *sparse_image;
2901 
2902  size_t
2903  number_colors;
2904 
2905  assert(image != (Image *) NULL);
2906  assert(image->signature == MagickCoreSignature);
2907  if (image->debug != MagickFalse)
2908  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2909  assert(exception != (ExceptionInfo *) NULL);
2910  assert(exception->signature == MagickCoreSignature);
2911 
2912  /* Determine number of color values needed per control point */
2913  number_colors=0;
2914  if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2915  number_colors++;
2916  if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2917  number_colors++;
2918  if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2919  number_colors++;
2920  if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2921  (image->colorspace == CMYKColorspace))
2922  number_colors++;
2923  if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2924  (image->alpha_trait != UndefinedPixelTrait))
2925  number_colors++;
2926 
2927  /*
2928  Convert input arguments into mapping coefficients, this this case
2929  we are mapping (distorting) colors, rather than coordinates.
2930  */
2931  { DistortMethod
2932  distort_method;
2933 
2934  distort_method=(DistortMethod) method;
2935  if ( distort_method >= SentinelDistortion )
2936  distort_method = ShepardsDistortion; /* Pretend to be Shepards */
2937  coeff = GenerateCoefficients(image, &distort_method, number_arguments,
2938  arguments, number_colors, exception);
2939  if ( coeff == (double *) NULL )
2940  return((Image *) NULL);
2941  /*
2942  Note some Distort Methods may fall back to other simpler methods,
2943  Currently the only fallback of concern is Bilinear to Affine
2944  (Barycentric), which is alaso sparse_colr method. This also ensures
2945  correct two and one color Barycentric handling.
2946  */
2947  sparse_method = (SparseColorMethod) distort_method;
2948  if ( distort_method == ShepardsDistortion )
2949  sparse_method = method; /* return non-distort methods to normal */
2950  if ( sparse_method == InverseColorInterpolate )
2951  coeff[0]=0.5; /* sqrt() the squared distance for inverse */
2952  }
2953 
2954  /* Verbose output */
2955  if (IsStringTrue(GetImageArtifact(image,"verbose")) != MagickFalse) {
2956 
2957  switch (sparse_method) {
2959  {
2960  register ssize_t x=0;
2961  (void) FormatLocaleFile(stderr, "Barycentric Sparse Color:\n");
2962  if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2963  (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf' \\\n",
2964  coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2965  if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2966  (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf' \\\n",
2967  coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2968  if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2969  (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf' \\\n",
2970  coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2971  if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2972  (image->colorspace == CMYKColorspace))
2973  (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf' \\\n",
2974  coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2975  if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2976  (image->alpha_trait != UndefinedPixelTrait))
2977  (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf' \\\n",
2978  coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2979  break;
2980  }
2982  {
2983  register ssize_t x=0;
2984  (void) FormatLocaleFile(stderr, "Bilinear Sparse Color\n");
2985  if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2986  (void) FormatLocaleFile(stderr, " -channel R -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2987  coeff[ x ], coeff[x+1],
2988  coeff[x+2], coeff[x+3]),x+=4;
2989  if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2990  (void) FormatLocaleFile(stderr, " -channel G -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2991  coeff[ x ], coeff[x+1],
2992  coeff[x+2], coeff[x+3]),x+=4;
2993  if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2994  (void) FormatLocaleFile(stderr, " -channel B -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2995  coeff[ x ], coeff[x+1],
2996  coeff[x+2], coeff[x+3]),x+=4;
2997  if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2998  (image->colorspace == CMYKColorspace))
2999  (void) FormatLocaleFile(stderr, " -channel K -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3000  coeff[ x ], coeff[x+1],
3001  coeff[x+2], coeff[x+3]),x+=4;
3002  if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3003  (image->alpha_trait != UndefinedPixelTrait))
3004  (void) FormatLocaleFile(stderr, " -channel A -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
3005  coeff[ x ], coeff[x+1],
3006  coeff[x+2], coeff[x+3]),x+=4;
3007  break;
3008  }
3009  default:
3010  /* sparse color method is too complex for FX emulation */
3011  break;
3012  }
3013  }
3014 
3015  /* Generate new image for generated interpolated gradient.
3016  * ASIDE: Actually we could have just replaced the colors of the original
3017  * image, but IM Core policy, is if storage class could change then clone
3018  * the image.
3019  */
3020 
3021  sparse_image=CloneImage(image,0,0,MagickTrue,exception);
3022  if (sparse_image == (Image *) NULL)
3023  return((Image *) NULL);
3024  if (SetImageStorageClass(sparse_image,DirectClass,exception) == MagickFalse)
3025  { /* if image is ColorMapped - change it to DirectClass */
3026  sparse_image=DestroyImage(sparse_image);
3027  return((Image *) NULL);
3028  }
3029  { /* ----- MAIN CODE ----- */
3030  CacheView
3031  *sparse_view;
3032 
3034  status;
3035 
3037  progress;
3038 
3039  ssize_t
3040  j;
3041 
3042  status=MagickTrue;
3043  progress=0;
3044  sparse_view=AcquireAuthenticCacheView(sparse_image,exception);
3045 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3046  #pragma omp parallel for schedule(static) shared(progress,status) \
3047  magick_number_threads(image,sparse_image,sparse_image->rows,1)
3048 #endif
3049  for (j=0; j < (ssize_t) sparse_image->rows; j++)
3050  {
3052  sync;
3053 
3054  PixelInfo
3055  pixel; /* pixel to assign to distorted image */
3056 
3057  register ssize_t
3058  i;
3059 
3060  register Quantum
3061  *magick_restrict q;
3062 
3063  q=GetCacheViewAuthenticPixels(sparse_view,0,j,sparse_image->columns,
3064  1,exception);
3065  if (q == (Quantum *) NULL)
3066  {
3067  status=MagickFalse;
3068  continue;
3069  }
3070  GetPixelInfo(sparse_image,&pixel);
3071  for (i=0; i < (ssize_t) image->columns; i++)
3072  {
3073  GetPixelInfoPixel(image,q,&pixel);
3074  switch (sparse_method)
3075  {
3077  {
3078  register ssize_t x=0;
3079  if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3080  pixel.red = coeff[x]*i +coeff[x+1]*j
3081  +coeff[x+2], x+=3;
3082  if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3083  pixel.green = coeff[x]*i +coeff[x+1]*j
3084  +coeff[x+2], x+=3;
3085  if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3086  pixel.blue = coeff[x]*i +coeff[x+1]*j
3087  +coeff[x+2], x+=3;
3088  if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3089  (image->colorspace == CMYKColorspace))
3090  pixel.black = coeff[x]*i +coeff[x+1]*j
3091  +coeff[x+2], x+=3;
3092  if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3093  (image->alpha_trait != UndefinedPixelTrait))
3094  pixel.alpha = coeff[x]*i +coeff[x+1]*j
3095  +coeff[x+2], x+=3;
3096  break;
3097  }
3099  {
3100  register ssize_t x=0;
3101  if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3102  pixel.red = coeff[x]*i + coeff[x+1]*j +
3103  coeff[x+2]*i*j + coeff[x+3], x+=4;
3104  if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3105  pixel.green = coeff[x]*i + coeff[x+1]*j +
3106  coeff[x+2]*i*j + coeff[x+3], x+=4;
3107  if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3108  pixel.blue = coeff[x]*i + coeff[x+1]*j +
3109  coeff[x+2]*i*j + coeff[x+3], x+=4;
3110  if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3111  (image->colorspace == CMYKColorspace))
3112  pixel.black = coeff[x]*i + coeff[x+1]*j +
3113  coeff[x+2]*i*j + coeff[x+3], x+=4;
3114  if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3115  (image->alpha_trait != UndefinedPixelTrait))
3116  pixel.alpha = coeff[x]*i + coeff[x+1]*j +
3117  coeff[x+2]*i*j + coeff[x+3], x+=4;
3118  break;
3119  }
3122  { /* Inverse (Squared) Distance weights average (IDW) */
3123  size_t
3124  k;
3125  double
3126  denominator;
3127 
3128  if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3129  pixel.red=0.0;
3130  if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3131  pixel.green=0.0;
3132  if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3133  pixel.blue=0.0;
3134  if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3135  (image->colorspace == CMYKColorspace))
3136  pixel.black=0.0;
3137  if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3138  (image->alpha_trait != UndefinedPixelTrait))
3139  pixel.alpha=0.0;
3140  denominator = 0.0;
3141  for(k=0; k<number_arguments; k+=2+number_colors) {
3142  register ssize_t x=(ssize_t) k+2;
3143  double weight =
3144  ((double)i-arguments[ k ])*((double)i-arguments[ k ])
3145  + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
3146  weight = pow(weight,coeff[0]); /* inverse of power factor */
3147  weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
3148  if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3149  pixel.red += arguments[x++]*weight;
3150  if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3151  pixel.green += arguments[x++]*weight;
3152  if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3153  pixel.blue += arguments[x++]*weight;
3154  if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3155  (image->colorspace == CMYKColorspace))
3156  pixel.black += arguments[x++]*weight;
3157  if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3158  (image->alpha_trait != UndefinedPixelTrait))
3159  pixel.alpha += arguments[x++]*weight;
3160  denominator += weight;
3161  }
3162  if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3163  pixel.red/=denominator;
3164  if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3165  pixel.green/=denominator;
3166  if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3167  pixel.blue/=denominator;
3168  if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3169  (image->colorspace == CMYKColorspace))
3170  pixel.black/=denominator;
3171  if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3172  (image->alpha_trait != UndefinedPixelTrait))
3173  pixel.alpha/=denominator;
3174  break;
3175  }
3177  {
3178  size_t
3179  k;
3180 
3181  double
3182  minimum = MagickMaximumValue;
3183 
3184  /*
3185  Just use the closest control point you can find!
3186  */
3187  for(k=0; k<number_arguments; k+=2+number_colors) {
3188  double distance =
3189  fabs((double)i-arguments[ k ])
3190  + fabs((double)j-arguments[k+1]);
3191  if ( distance < minimum ) {
3192  register ssize_t x=(ssize_t) k+2;
3193  if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3194  pixel.red=arguments[x++];
3195  if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3196  pixel.green=arguments[x++];
3197  if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3198  pixel.blue=arguments[x++];
3199  if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3200  (image->colorspace == CMYKColorspace))
3201  pixel.black=arguments[x++];
3202  if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3203  (image->alpha_trait != UndefinedPixelTrait))
3204  pixel.alpha=arguments[x++];
3205  minimum = distance;
3206  }
3207  }
3208  break;
3209  }
3211  default:
3212  {
3213  size_t
3214  k;
3215 
3216  double
3217  minimum = MagickMaximumValue;
3218 
3219  /*
3220  Just use the closest control point you can find!
3221  */
3222  for (k=0; k<number_arguments; k+=2+number_colors) {
3223  double distance =
3224  ((double)i-arguments[ k ])*((double)i-arguments[ k ])
3225  + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
3226  if ( distance < minimum ) {
3227  register ssize_t x=(ssize_t) k+2;
3228  if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3229  pixel.red=arguments[x++];
3230  if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3231  pixel.green=arguments[x++];
3232  if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3233  pixel.blue=arguments[x++];
3234  if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3235  (image->colorspace == CMYKColorspace))
3236  pixel.black=arguments[x++];
3237  if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3238  (image->alpha_trait != UndefinedPixelTrait))
3239  pixel.alpha=arguments[x++];
3240  minimum = distance;
3241  }
3242  }
3243  break;
3244  }
3245  }
3246  /* set the color directly back into the source image */
3247  if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3248  pixel.red=(MagickRealType) ClampPixel(QuantumRange*pixel.red);
3249  if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3251  if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3253  if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3254  (image->colorspace == CMYKColorspace))
3256  if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3257  (image->alpha_trait != UndefinedPixelTrait))
3259  SetPixelViaPixelInfo(sparse_image,&pixel,q);
3260  q+=GetPixelChannels(sparse_image);
3261  }
3262  sync=SyncCacheViewAuthenticPixels(sparse_view,exception);
3263  if (sync == MagickFalse)
3264  status=MagickFalse;
3265  if (image->progress_monitor != (MagickProgressMonitor) NULL)
3266  {
3268  proceed;
3269 
3270 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3271  #pragma omp critical (MagickCore_SparseColorImage)
3272 #endif
3273  proceed=SetImageProgress(image,SparseColorTag,progress++,image->rows);
3274  if (proceed == MagickFalse)
3275  status=MagickFalse;
3276  }
3277  }
3278  sparse_view=DestroyCacheView(sparse_view);
3279  if (status == MagickFalse)
3280  sparse_image=DestroyImage(sparse_image);
3281  }
3282  coeff = (double *) RelinquishMagickMemory(coeff);
3283  return(sparse_image);
3284 }
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:176
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:2887
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:188
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:472
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:188
#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:539
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:1581
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:377
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:2669
#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:1498
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:188
#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:3409
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:1058
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:2802
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:2602
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:188
#define MagickMin(x, y)
Definition: image-private.h:27
MagickExport void * RelinquishMagickMemory(void *memory)
Definition: memory.c:1027
MagickRealType green
Definition: pixel.h:188
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:1674
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:1183
MagickExport Image * CloneImage(const Image *image, const size_t columns, const size_t rows, const MagickBooleanType detach, ExceptionInfo *exception)
Definition: image.c:800
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)