granicus.if.org Git - imagemagick/blob - magick/distort.c

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