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