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

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