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