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