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