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