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);
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,exception);
1502       (void) SetImageAlphaChannel(resize_alpha,DeactivateAlphaChannel,exception);
1503       (void) CompositeImage(resize_image,CopyOpacityCompositeOp,resize_alpha,
1504         0,0);
1505       resize_alpha=DestroyImage(resize_alpha);
1506     }
1507   (void) SetImageVirtualPixelMethod(resize_image,vp_save);
1508
1509   /*
1510     Clean up the results of the Distortion
1511   */
1512   crop_area.width=columns;
1513   crop_area.height=rows;
1514   crop_area.x=0;
1515   crop_area.y=0;
1516
1517   tmp_image=resize_image;
1518   resize_image=CropImage(tmp_image,&crop_area,exception);
1519   tmp_image=DestroyImage(tmp_image);
1520
1521   if ( resize_image == (Image *) NULL )
1522     return((Image *) NULL);
1523
1524   return(resize_image);
1525 }
1526 \f
1527 /*
1528 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1529 %                                                                             %
1530 %                                                                             %
1531 %                                                                             %
1532 %   D i s t o r t I m a g e                                                   %
1533 %                                                                             %
1534 %                                                                             %
1535 %                                                                             %
1536 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1537 %
1538 %  DistortImage() distorts an image using various distortion methods, by
1539 %  mapping color lookups of the source image to a new destination image
1540 %  usally of the same size as the source image, unless 'bestfit' is set to
1541 %  true.
1542 %
1543 %  If 'bestfit' is enabled, and distortion allows it, the destination image is
1544 %  adjusted to ensure the whole source 'image' will just fit within the final
1545 %  destination image, which will be sized and offset accordingly.  Also in
1546 %  many cases the virtual offset of the source image will be taken into
1547 %  account in the mapping.
1548 %
1549 %  If the '-verbose' control option has been set print to standard error the
1550 %  equicelent '-fx' formula with coefficients for the function, if practical.
1551 %
1552 %  The format of the DistortImage() method is:
1553 %
1554 %      Image *DistortImage(const Image *image,const DistortImageMethod method,
1555 %        const size_t number_arguments,const double *arguments,
1556 %        MagickBooleanType bestfit, ExceptionInfo *exception)
1557 %
1558 %  A description of each parameter follows:
1559 %
1560 %    o image: the image to be distorted.
1561 %
1562 %    o method: the method of image distortion.
1563 %
1564 %        ArcDistortion always ignores source image offset, and always
1565 %        'bestfit' the destination image with the top left corner offset
1566 %        relative to the polar mapping center.
1567 %
1568 %        Affine, Perspective, and Bilinear, do least squares fitting of the
1569 %        distrotion when more than the minimum number of control point pairs
1570 %        are provided.
1571 %
1572 %        Perspective, and Bilinear, fall back to a Affine distortion when less
1573 %        than 4 control point pairs are provided.  While Affine distortions
1574 %        let you use any number of control point pairs, that is Zero pairs is
1575 %        a No-Op (viewport only) distortion, one pair is a translation and
1576 %        two pairs of control points do a scale-rotate-translate, without any
1577 %        shearing.
1578 %
1579 %    o number_arguments: the number of arguments given.
1580 %
1581 %    o arguments: an array of floating point arguments for this method.
1582 %
1583 %    o bestfit: Attempt to 'bestfit' the size of the resulting image.
1584 %        This also forces the resulting image to be a 'layered' virtual
1585 %        canvas image.  Can be overridden using 'distort:viewport' setting.
1586 %
1587 %    o exception: return any errors or warnings in this structure
1588 %
1589 %  Extra Controls from Image meta-data (artifacts)...
1590 %
1591 %    o "verbose"
1592 %        Output to stderr alternatives, internal coefficents, and FX
1593 %        equivalents for the distortion operation (if feasible).
1594 %        This forms an extra check of the distortion method, and allows users
1595 %        access to the internal constants IM calculates for the distortion.
1596 %
1597 %    o "distort:viewport"
1598 %        Directly set the output image canvas area and offest to use for the
1599 %        resulting image, rather than use the original images canvas, or a
1600 %        calculated 'bestfit' canvas.
1601 %
1602 %    o "distort:scale"
1603 %        Scale the size of the output canvas by this amount to provide a
1604 %        method of Zooming, and for super-sampling the results.
1605 %
1606 %  Other settings that can effect results include
1607 %
1608 %    o 'interpolate' For source image lookups (scale enlargements)
1609 %
1610 %    o 'filter'      Set filter to use for area-resampling (scale shrinking).
1611 %                    Set to 'point' to turn off and use 'interpolate' lookup
1612 %                    instead
1613 %
1614 */
1615
1616 MagickExport Image *DistortImage(const Image *image,DistortImageMethod method,
1617   const size_t number_arguments,const double *arguments,
1618   MagickBooleanType bestfit,ExceptionInfo *exception)
1619 {
1620 #define DistortImageTag  "Distort/Image"
1621
1622   double
1623     *coeff,
1624     output_scaling;
1625
1626   Image
1627     *distort_image;
1628
1629   RectangleInfo
1630     geometry;  /* geometry of the distorted space viewport */
1631
1632   MagickBooleanType
1633     viewport_given;
1634
1635   assert(image != (Image *) NULL);
1636   assert(image->signature == MagickSignature);
1637   if (image->debug != MagickFalse)
1638     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1639   assert(exception != (ExceptionInfo *) NULL);
1640   assert(exception->signature == MagickSignature);
1641
1642
1643   /*
1644     Handle Special Compound Distortions
1645   */
1646   if ( method == ResizeDistortion )
1647     {
1648       if ( number_arguments != 2 )
1649         {
1650           (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1651                     "InvalidArgument","%s : '%s'","Resize",
1652                     "Invalid number of args: 2 only");
1653           return((Image *) NULL);
1654         }
1655       distort_image=DistortResizeImage(image,(size_t)arguments[0],
1656          (size_t)arguments[1], exception);
1657       return(distort_image);
1658     }
1659
1660   /*
1661     Convert input arguments (usually as control points for reverse mapping)
1662     into mapping coefficients to apply the distortion.
1663
1664     Note that some distortions are mapped to other distortions,
1665     and as such do not require specific code after this point.
1666   */
1667   coeff = GenerateCoefficients(image, &method, number_arguments,
1668       arguments, 0, exception);
1669   if ( coeff == (double *) NULL )
1670     return((Image *) NULL);
1671
1672   /*
1673     Determine the size and offset for a 'bestfit' destination.
1674     Usally the four corners of the source image is enough.
1675   */
1676
1677   /* default output image bounds, when no 'bestfit' is requested */
1678   geometry.width=image->columns;
1679   geometry.height=image->rows;
1680   geometry.x=0;
1681   geometry.y=0;
1682
1683   if ( method == ArcDistortion ) {
1684     bestfit = MagickTrue;  /* always calculate a 'best fit' viewport */
1685   }
1686
1687   /* Work out the 'best fit', (required for ArcDistortion) */
1688   if ( bestfit ) {
1689     PointInfo
1690       s,d,min,max;  /* source, dest coords --mapping--> min, max coords */
1691
1692     MagickBooleanType
1693       fix_bounds = MagickTrue;   /* enlarge bounds for VP handling */
1694
1695     s.x=s.y=min.x=max.x=min.y=max.y=0.0;   /* keep compiler happy */
1696
1697 /* defines to figure out the bounds of the distorted image */
1698 #define InitalBounds(p) \
1699 { \
1700   /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1701   min.x = max.x = p.x; \
1702   min.y = max.y = p.y; \
1703 }
1704 #define ExpandBounds(p) \
1705 { \
1706   /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1707   min.x = MagickMin(min.x,p.x); \
1708   max.x = MagickMax(max.x,p.x); \
1709   min.y = MagickMin(min.y,p.y); \
1710   max.y = MagickMax(max.y,p.y); \
1711 }
1712
1713     switch (method)
1714     {
1715       case AffineDistortion:
1716       { double inverse[6];
1717         InvertAffineCoefficients(coeff, inverse);
1718         s.x = (double) image->page.x;
1719         s.y = (double) image->page.y;
1720         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1721         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1722         InitalBounds(d);
1723         s.x = (double) image->page.x+image->columns;
1724         s.y = (double) image->page.y;
1725         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1726         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1727         ExpandBounds(d);
1728         s.x = (double) image->page.x;
1729         s.y = (double) image->page.y+image->rows;
1730         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1731         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1732         ExpandBounds(d);
1733         s.x = (double) image->page.x+image->columns;
1734         s.y = (double) image->page.y+image->rows;
1735         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1736         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1737         ExpandBounds(d);
1738         break;
1739       }
1740       case PerspectiveDistortion:
1741       { double inverse[8], scale;
1742         InvertPerspectiveCoefficients(coeff, inverse);
1743         s.x = (double) image->page.x;
1744         s.y = (double) image->page.y;
1745         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1746         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1747         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1748         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1749         InitalBounds(d);
1750         s.x = (double) image->page.x+image->columns;
1751         s.y = (double) image->page.y;
1752         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1753         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1754         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1755         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1756         ExpandBounds(d);
1757         s.x = (double) image->page.x;
1758         s.y = (double) image->page.y+image->rows;
1759         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1760         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1761         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1762         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1763         ExpandBounds(d);
1764         s.x = (double) image->page.x+image->columns;
1765         s.y = (double) image->page.y+image->rows;
1766         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1767         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1768         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1769         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1770         ExpandBounds(d);
1771         break;
1772       }
1773       case ArcDistortion:
1774       { double a, ca, sa;
1775         /* Forward Map Corners */
1776         a = coeff[0]-coeff[1]/2; ca = cos(a); sa = sin(a);
1777         d.x = coeff[2]*ca;
1778         d.y = coeff[2]*sa;
1779         InitalBounds(d);
1780         d.x = (coeff[2]-coeff[3])*ca;
1781         d.y = (coeff[2]-coeff[3])*sa;
1782         ExpandBounds(d);
1783         a = coeff[0]+coeff[1]/2; ca = cos(a); sa = sin(a);
1784         d.x = coeff[2]*ca;
1785         d.y = coeff[2]*sa;
1786         ExpandBounds(d);
1787         d.x = (coeff[2]-coeff[3])*ca;
1788         d.y = (coeff[2]-coeff[3])*sa;
1789         ExpandBounds(d);
1790         /* Orthogonal points along top of arc */
1791         for( a=(double) (ceil((double) ((coeff[0]-coeff[1]/2.0)/MagickPI2))*MagickPI2);
1792                a<(coeff[0]+coeff[1]/2.0); a+=MagickPI2 ) {
1793           ca = cos(a); sa = sin(a);
1794           d.x = coeff[2]*ca;
1795           d.y = coeff[2]*sa;
1796           ExpandBounds(d);
1797         }
1798         /*
1799           Convert the angle_to_width and radius_to_height
1800           to appropriate scaling factors, to allow faster processing
1801           in the mapping function.
1802         */
1803         coeff[1] = (double) (Magick2PI*image->columns/coeff[1]);
1804         coeff[3] = (double)image->rows/coeff[3];
1805         break;
1806       }
1807       case PolarDistortion:
1808       {
1809         if (number_arguments < 2)
1810           coeff[2] = coeff[3] = 0.0;
1811         min.x = coeff[2]-coeff[0];
1812         max.x = coeff[2]+coeff[0];
1813         min.y = coeff[3]-coeff[0];
1814         max.y = coeff[3]+coeff[0];
1815         /* should be about 1.0 if Rmin = 0 */
1816         coeff[7]=(double) geometry.height/(coeff[0]-coeff[1]);
1817         break;
1818       }
1819       case DePolarDistortion:
1820       {
1821         /* direct calculation as it needs to tile correctly
1822          * for reversibility in a DePolar-Polar cycle */
1823         fix_bounds = MagickFalse;
1824         geometry.x = geometry.y = 0;
1825         geometry.height = (size_t) ceil(coeff[0]-coeff[1]);
1826         geometry.width = (size_t)
1827                   ceil((coeff[0]-coeff[1])*(coeff[5]-coeff[4])*0.5);
1828         /* correct scaling factors relative to new size */
1829         coeff[6]=(coeff[5]-coeff[4])/geometry.width; /* changed width */
1830         coeff[7]=(coeff[0]-coeff[1])/geometry.height; /* should be about 1.0 */
1831         break;
1832       }
1833       case Cylinder2PlaneDistortion:
1834       {
1835         /* direct calculation so center of distortion is either a pixel
1836          * center, or pixel edge. This allows for reversibility of the
1837          * distortion */
1838         geometry.x = geometry.y = 0;
1839         geometry.width = (size_t) ceil( 2.0*coeff[1]*tan(coeff[0]/2.0) );
1840         geometry.height = (size_t) ceil( 2.0*coeff[3]/cos(coeff[0]/2.0) );
1841         /* correct center of distortion relative to new size */
1842         coeff[4] = (double) geometry.width/2.0;
1843         coeff[5] = (double) geometry.height/2.0;
1844         fix_bounds = MagickFalse;
1845         break;
1846       }
1847       case Plane2CylinderDistortion:
1848       {
1849         /* direct calculation center is either pixel center, or pixel edge
1850          * so as to allow reversibility of the image distortion */
1851         geometry.x = geometry.y = 0;
1852         geometry.width = (size_t) ceil(coeff[0]*coeff[1]);  /* FOV * radius */
1853         geometry.height = (size_t) (2*coeff[3]);              /* input image height */
1854         /* correct center of distortion relative to new size */
1855         coeff[4] = (double) geometry.width/2.0;
1856         coeff[5] = (double) geometry.height/2.0;
1857         fix_bounds = MagickFalse;
1858         break;
1859       }
1860       case ShepardsDistortion:
1861       case BilinearForwardDistortion:
1862       case BilinearReverseDistortion:
1863 #if 0
1864       case QuadrilateralDistortion:
1865 #endif
1866       case PolynomialDistortion:
1867       case BarrelDistortion:
1868       case BarrelInverseDistortion:
1869       default:
1870         /* no calculated bestfit available for these distortions */
1871         bestfit = MagickFalse;
1872         fix_bounds = MagickFalse;
1873         break;
1874     }
1875
1876     /* Set the output image geometry to calculated 'bestfit'.
1877        Yes this tends to 'over do' the file image size, ON PURPOSE!
1878        Do not do this for DePolar which needs to be exact for virtual tiling.
1879     */
1880     if ( fix_bounds ) {
1881       geometry.x = (ssize_t) floor(min.x-0.5);
1882       geometry.y = (ssize_t) floor(min.y-0.5);
1883       geometry.width=(size_t) ceil(max.x-geometry.x+0.5);
1884       geometry.height=(size_t) ceil(max.y-geometry.y+0.5);
1885     }
1886
1887   }  /* end bestfit destination image calculations */
1888
1889   /* The user provided a 'viewport' expert option which may
1890      overrides some parts of the current output image geometry.
1891      This also overrides its default 'bestfit' setting.
1892   */
1893   { const char *artifact=GetImageArtifact(image,"distort:viewport");
1894     viewport_given = MagickFalse;
1895     if ( artifact != (const char *) NULL ) {
1896       (void) ParseAbsoluteGeometry(artifact,&geometry);
1897       viewport_given = MagickTrue;
1898     }
1899   }
1900
1901   /* Verbose output */
1902   if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
1903     register ssize_t
1904        i;
1905     char image_gen[MaxTextExtent];
1906     const char *lookup;
1907
1908     /* Set destination image size and virtual offset */
1909     if ( bestfit || viewport_given ) {
1910       (void) FormatLocaleString(image_gen, MaxTextExtent,"  -size %.20gx%.20g "
1911         "-page %+.20g%+.20g xc: +insert \\\n",(double) geometry.width,
1912         (double) geometry.height,(double) geometry.x,(double) geometry.y);
1913       lookup="v.p{ xx-v.page.x-.5, yy-v.page.y-.5 }";
1914     }
1915     else {
1916       image_gen[0] = '\0';             /* no destination to generate */
1917       lookup = "p{ xx-page.x-.5, yy-page.y-.5 }"; /* simplify lookup */
1918     }
1919
1920     switch (method) {
1921       case AffineDistortion:
1922       {
1923         double *inverse;
1924
1925         inverse = (double *) AcquireQuantumMemory(6,sizeof(*inverse));
1926         if (inverse == (double *) NULL) {
1927           coeff = (double *) RelinquishMagickMemory(coeff);
1928           (void) ThrowMagickException(exception,GetMagickModule(),
1929                   ResourceLimitError,"MemoryAllocationFailed",
1930                   "%s", "DistortImages");
1931           return((Image *) NULL);
1932         }
1933         InvertAffineCoefficients(coeff, inverse);
1934         CoefficientsToAffineArgs(inverse);
1935         (void) FormatLocaleFile(stderr, "Affine Projection:\n");
1936         (void) FormatLocaleFile(stderr, "  -distort AffineProjection \\\n      '");
1937         for (i=0; i < 5; i++)
1938           (void) FormatLocaleFile(stderr, "%lf,", inverse[i]);
1939         (void) FormatLocaleFile(stderr, "%lf'\n", inverse[5]);
1940         inverse = (double *) RelinquishMagickMemory(inverse);
1941
1942         (void) FormatLocaleFile(stderr, "Affine Distort, FX Equivelent:\n");
1943         (void) FormatLocaleFile(stderr, "%s", image_gen);
1944         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1945         (void) FormatLocaleFile(stderr, "       xx=%+lf*ii %+lf*jj %+lf;\n",
1946             coeff[0], coeff[1], coeff[2]);
1947         (void) FormatLocaleFile(stderr, "       yy=%+lf*ii %+lf*jj %+lf;\n",
1948             coeff[3], coeff[4], coeff[5]);
1949         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
1950
1951         break;
1952       }
1953
1954       case PerspectiveDistortion:
1955       {
1956         double *inverse;
1957
1958         inverse = (double *) AcquireQuantumMemory(8,sizeof(*inverse));
1959         if (inverse == (double *) NULL) {
1960           coeff = (double *) RelinquishMagickMemory(coeff);
1961           (void) ThrowMagickException(exception,GetMagickModule(),
1962                   ResourceLimitError,"MemoryAllocationFailed",
1963                   "%s", "DistortCoefficients");
1964           return((Image *) NULL);
1965         }
1966         InvertPerspectiveCoefficients(coeff, inverse);
1967         (void) FormatLocaleFile(stderr, "Perspective Projection:\n");
1968         (void) FormatLocaleFile(stderr, "  -distort PerspectiveProjection \\\n      '");
1969         for (i=0; i<4; i++)
1970           (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
1971         (void) FormatLocaleFile(stderr, "\n       ");
1972         for (; i<7; i++)
1973           (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
1974         (void) FormatLocaleFile(stderr, "%lf'\n", inverse[7]);
1975         inverse = (double *) RelinquishMagickMemory(inverse);
1976
1977         (void) FormatLocaleFile(stderr, "Perspective Distort, FX Equivelent:\n");
1978         (void) FormatLocaleFile(stderr, "%s", image_gen);
1979         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1980         (void) FormatLocaleFile(stderr, "       rr=%+lf*ii %+lf*jj + 1;\n",
1981             coeff[6], coeff[7]);
1982         (void) FormatLocaleFile(stderr, "       xx=(%+lf*ii %+lf*jj %+lf)/rr;\n",
1983             coeff[0], coeff[1], coeff[2]);
1984         (void) FormatLocaleFile(stderr, "       yy=(%+lf*ii %+lf*jj %+lf)/rr;\n",
1985             coeff[3], coeff[4], coeff[5]);
1986         (void) FormatLocaleFile(stderr, "       rr%s0 ? %s : blue' \\\n",
1987             coeff[8] < 0 ? "<" : ">", lookup);
1988         break;
1989       }
1990
1991       case BilinearForwardDistortion:
1992         (void) FormatLocaleFile(stderr, "BilinearForward Mapping Equations:\n");
1993         (void) FormatLocaleFile(stderr, "%s", image_gen);
1994         (void) FormatLocaleFile(stderr, "    i = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
1995             coeff[0], coeff[1], coeff[2], coeff[3]);
1996         (void) FormatLocaleFile(stderr, "    j = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
1997             coeff[4], coeff[5], coeff[6], coeff[7]);
1998 #if 0
1999         /* for debugging */
2000         (void) FormatLocaleFile(stderr, "   c8 = %+lf  c9 = 2*a = %+lf;\n",
2001             coeff[8], coeff[9]);
2002 #endif
2003         (void) FormatLocaleFile(stderr, "BilinearForward Distort, FX Equivelent:\n");
2004         (void) FormatLocaleFile(stderr, "%s", image_gen);
2005         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
2006             0.5-coeff[3], 0.5-coeff[7]);
2007         (void) FormatLocaleFile(stderr, "       bb=%lf*ii %+lf*jj %+lf;\n",
2008             coeff[6], -coeff[2], coeff[8]);
2009         /* Handle Special degenerate (non-quadratic) or trapezoidal case */
2010         if ( coeff[9] != 0 ) {
2011           (void) FormatLocaleFile(stderr, "       rt=bb*bb %+lf*(%lf*ii%+lf*jj);\n",
2012               -2*coeff[9],  coeff[4], -coeff[0]);
2013           (void) FormatLocaleFile(stderr, "       yy=( -bb + sqrt(rt) ) / %lf;\n",
2014                coeff[9]);
2015         } else
2016           (void) FormatLocaleFile(stderr, "       yy=(%lf*ii%+lf*jj)/bb;\n",
2017                 -coeff[4], coeff[0]);
2018         (void) FormatLocaleFile(stderr, "       xx=(ii %+lf*yy)/(%lf %+lf*yy);\n",
2019              -coeff[1], coeff[0], coeff[2]);
2020         if ( coeff[9] != 0 )
2021           (void) FormatLocaleFile(stderr, "       (rt < 0 ) ? red : %s'\n", lookup);
2022         else
2023           (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
2024         break;
2025
2026       case BilinearReverseDistortion:
2027 #if 0
2028         (void) FormatLocaleFile(stderr, "Polynomial Projection Distort:\n");
2029         (void) FormatLocaleFile(stderr, "  -distort PolynomialProjection \\\n");
2030         (void) FormatLocaleFile(stderr, "      '1.5, %lf, %lf, %lf, %lf,\n",
2031             coeff[3], coeff[0], coeff[1], coeff[2]);
2032         (void) FormatLocaleFile(stderr, "            %lf, %lf, %lf, %lf'\n",
2033             coeff[7], coeff[4], coeff[5], coeff[6]);
2034 #endif
2035         (void) FormatLocaleFile(stderr, "BilinearReverse Distort, FX Equivelent:\n");
2036         (void) FormatLocaleFile(stderr, "%s", image_gen);
2037         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2038         (void) FormatLocaleFile(stderr, "       xx=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
2039             coeff[0], coeff[1], coeff[2], coeff[3]);
2040         (void) FormatLocaleFile(stderr, "       yy=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
2041             coeff[4], coeff[5], coeff[6], coeff[7]);
2042         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
2043         break;
2044
2045       case PolynomialDistortion:
2046       {
2047         size_t nterms = (size_t) coeff[1];
2048         (void) FormatLocaleFile(stderr, "Polynomial (order %lg, terms %lu), FX Equivelent\n",
2049           coeff[0],(unsigned long) nterms);
2050         (void) FormatLocaleFile(stderr, "%s", image_gen);
2051         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2052         (void) FormatLocaleFile(stderr, "       xx =");
2053         for (i=0; i<(ssize_t) nterms; i++) {
2054           if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n         ");
2055           (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i],
2056                poly_basis_str(i));
2057         }
2058         (void) FormatLocaleFile(stderr, ";\n       yy =");
2059         for (i=0; i<(ssize_t) nterms; i++) {
2060           if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n         ");
2061           (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i+nterms],
2062                poly_basis_str(i));
2063         }
2064         (void) FormatLocaleFile(stderr, ";\n       %s' \\\n", lookup);
2065         break;
2066       }
2067       case ArcDistortion:
2068       {
2069         (void) FormatLocaleFile(stderr, "Arc Distort, Internal Coefficients:\n");
2070         for ( i=0; i<5; i++ )
2071           (void) FormatLocaleFile(stderr, "  c%.20g = %+lf\n", (double) i, coeff[i]);
2072         (void) FormatLocaleFile(stderr, "Arc Distort, FX Equivelent:\n");
2073         (void) FormatLocaleFile(stderr, "%s", image_gen);
2074         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x; jj=j+page.y;\n");
2075         (void) FormatLocaleFile(stderr, "       xx=(atan2(jj,ii)%+lf)/(2*pi);\n",
2076                                   -coeff[0]);
2077         (void) FormatLocaleFile(stderr, "       xx=xx-round(xx);\n");
2078         (void) FormatLocaleFile(stderr, "       xx=xx*%lf %+lf;\n",
2079                             coeff[1], coeff[4]);
2080         (void) FormatLocaleFile(stderr, "       yy=(%lf - hypot(ii,jj)) * %lf;\n",
2081                             coeff[2], coeff[3]);
2082         (void) FormatLocaleFile(stderr, "       v.p{xx-.5,yy-.5}' \\\n");
2083         break;
2084       }
2085       case PolarDistortion:
2086       {
2087         (void) FormatLocaleFile(stderr, "Polar Distort, Internal Coefficents\n");
2088         for ( i=0; i<8; i++ )
2089           (void) FormatLocaleFile(stderr, "  c%.20g = %+lf\n", (double) i, coeff[i]);
2090         (void) FormatLocaleFile(stderr, "Polar Distort, FX Equivelent:\n");
2091         (void) FormatLocaleFile(stderr, "%s", image_gen);
2092         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
2093                          -coeff[2], -coeff[3]);
2094         (void) FormatLocaleFile(stderr, "       xx=(atan2(ii,jj)%+lf)/(2*pi);\n",
2095                          -(coeff[4]+coeff[5])/2 );
2096         (void) FormatLocaleFile(stderr, "       xx=xx-round(xx);\n");
2097         (void) FormatLocaleFile(stderr, "       xx=xx*2*pi*%lf + v.w/2;\n",
2098                          coeff[6] );
2099         (void) FormatLocaleFile(stderr, "       yy=(hypot(ii,jj)%+lf)*%lf;\n",
2100                          -coeff[1], coeff[7] );
2101         (void) FormatLocaleFile(stderr, "       v.p{xx-.5,yy-.5}' \\\n");
2102         break;
2103       }
2104       case DePolarDistortion:
2105       {
2106         (void) FormatLocaleFile(stderr, "DePolar Distort, Internal Coefficents\n");
2107         for ( i=0; i<8; i++ )
2108           (void) FormatLocaleFile(stderr, "  c%.20g = %+lf\n", (double) i, coeff[i]);
2109         (void) FormatLocaleFile(stderr, "DePolar Distort, FX Equivelent:\n");
2110         (void) FormatLocaleFile(stderr, "%s", image_gen);
2111         (void) FormatLocaleFile(stderr, "  -fx 'aa=(i+.5)*%lf %+lf;\n", coeff[6], -coeff[4] );
2112         (void) FormatLocaleFile(stderr, "       rr=(j+.5)*%lf %+lf;\n", coeff[7], +coeff[1] );
2113         (void) FormatLocaleFile(stderr, "       xx=rr*sin(aa) %+lf;\n", coeff[2] );
2114         (void) FormatLocaleFile(stderr, "       yy=rr*cos(aa) %+lf;\n", coeff[3] );
2115         (void) FormatLocaleFile(stderr, "       v.p{xx-.5,yy-.5}' \\\n");
2116         break;
2117       }
2118       case Cylinder2PlaneDistortion:
2119       {
2120         (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, Internal Coefficents\n");
2121         (void) FormatLocaleFile(stderr, "  cylinder_radius = %+lf\n", coeff[1]);
2122         (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, FX Equivelent:\n");
2123         (void) FormatLocaleFile(stderr, "%s", image_gen);
2124         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
2125                          -coeff[4], -coeff[5]);
2126         (void) FormatLocaleFile(stderr, "       aa=atan(ii/%+lf);\n", coeff[1] );
2127         (void) FormatLocaleFile(stderr, "       xx=%lf*aa%+lf;\n",
2128                          coeff[1], coeff[2] );
2129         (void) FormatLocaleFile(stderr, "       yy=jj*cos(aa)%+lf;\n", coeff[3] );
2130         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
2131         break;
2132       }
2133       case Plane2CylinderDistortion:
2134       {
2135         (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, Internal Coefficents\n");
2136         (void) FormatLocaleFile(stderr, "  cylinder_radius = %+lf\n", coeff[1]);
2137         (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, FX Equivelent:\n");
2138         (void) FormatLocaleFile(stderr, "%s", image_gen);
2139         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
2140                          -coeff[4], -coeff[5]);
2141         (void) FormatLocaleFile(stderr, "       ii=ii/%+lf;\n", coeff[1] );
2142         (void) FormatLocaleFile(stderr, "       xx=%lf*tan(ii)%+lf;\n",
2143                          coeff[1], coeff[2] );
2144         (void) FormatLocaleFile(stderr, "       yy=jj/cos(ii)%+lf;\n",
2145                          coeff[3] );
2146         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
2147         break;
2148       }
2149       case BarrelDistortion:
2150       case BarrelInverseDistortion:
2151       { double xc,yc;
2152         /* NOTE: This does the barrel roll in pixel coords not image coords
2153         ** The internal distortion must do it in image coordinates,
2154         ** so that is what the center coeff (8,9) is given in.
2155         */
2156         xc = ((double)image->columns-1.0)/2.0 + image->page.x;
2157         yc = ((double)image->rows-1.0)/2.0    + image->page.y;
2158         (void) FormatLocaleFile(stderr, "Barrel%s Distort, FX Equivelent:\n",
2159              method == BarrelDistortion ? "" : "Inv");
2160         (void) FormatLocaleFile(stderr, "%s", image_gen);
2161         if ( fabs(coeff[8]-xc-0.5) < 0.1 && fabs(coeff[9]-yc-0.5) < 0.1 )
2162           (void) FormatLocaleFile(stderr, "  -fx 'xc=(w-1)/2;  yc=(h-1)/2;\n");
2163         else
2164           (void) FormatLocaleFile(stderr, "  -fx 'xc=%lf;  yc=%lf;\n",
2165                coeff[8]-0.5, coeff[9]-0.5);
2166         (void) FormatLocaleFile(stderr,
2167              "       ii=i-xc;  jj=j-yc;  rr=hypot(ii,jj);\n");
2168         (void) FormatLocaleFile(stderr, "       ii=ii%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2169              method == BarrelDistortion ? "*" : "/",
2170              coeff[0],coeff[1],coeff[2],coeff[3]);
2171         (void) FormatLocaleFile(stderr, "       jj=jj%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2172              method == BarrelDistortion ? "*" : "/",
2173              coeff[4],coeff[5],coeff[6],coeff[7]);
2174         (void) FormatLocaleFile(stderr, "       v.p{fx*ii+xc,fy*jj+yc}' \\\n");
2175       }
2176       default:
2177         break;
2178     }
2179   }
2180
2181   /* The user provided a 'scale' expert option will scale the
2182      output image size, by the factor given allowing for super-sampling
2183      of the distorted image space.  Any scaling factors must naturally
2184      be halved as a result.
2185   */
2186   { const char *artifact;
2187     artifact=GetImageArtifact(image,"distort:scale");
2188     output_scaling = 1.0;
2189     if (artifact != (const char *) NULL) {
2190       output_scaling = fabs(InterpretLocaleValue(artifact,(char **) NULL));
2191       geometry.width  *= output_scaling;
2192       geometry.height *= output_scaling;
2193       geometry.x      *= output_scaling;
2194       geometry.y      *= output_scaling;
2195       if ( output_scaling < 0.1 ) {
2196         coeff = (double *) RelinquishMagickMemory(coeff);
2197         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
2198                 "InvalidArgument","%s", "-set option:distort:scale" );
2199         return((Image *) NULL);
2200       }
2201       output_scaling = 1/output_scaling;
2202     }
2203   }
2204 #define ScaleFilter(F,A,B,C,D) \
2205     ScaleResampleFilter( (F), \
2206       output_scaling*(A), output_scaling*(B), \
2207       output_scaling*(C), output_scaling*(D) )
2208
2209   /*
2210     Initialize the distort image attributes.
2211   */
2212   distort_image=CloneImage(image,geometry.width,geometry.height,MagickTrue,
2213     exception);
2214   if (distort_image == (Image *) NULL)
2215     return((Image *) NULL);
2216   /* if image is ColorMapped - change it to DirectClass */
2217   if (SetImageStorageClass(distort_image,DirectClass,exception) == MagickFalse)
2218     {
2219       distort_image=DestroyImage(distort_image);
2220       return((Image *) NULL);
2221     }
2222   distort_image->page.x=geometry.x;
2223   distort_image->page.y=geometry.y;
2224   if (distort_image->background_color.alpha != OpaqueAlpha)
2225     distort_image->matte=MagickTrue;
2226
2227   { /* ----- MAIN CODE -----
2228        Sample the source image to each pixel in the distort image.
2229      */
2230     CacheView
2231       *distort_view;
2232
2233     MagickBooleanType
2234       status;
2235
2236     MagickOffsetType
2237       progress;
2238
2239     PixelInfo
2240       zero;
2241
2242     ResampleFilter
2243       **restrict resample_filter;
2244
2245     ssize_t
2246       j;
2247
2248     status=MagickTrue;
2249     progress=0;
2250     GetPixelInfo(distort_image,&zero);
2251     resample_filter=AcquireResampleFilterThreadSet(image,
2252       UndefinedVirtualPixelMethod,MagickFalse,exception);
2253     distort_view=AcquireCacheView(distort_image);
2254 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2255   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2256 #endif
2257     for (j=0; j < (ssize_t) distort_image->rows; j++)
2258     {
2259       const int
2260         id = GetOpenMPThreadId();
2261
2262       double
2263         validity;  /* how mathematically valid is this the mapping */
2264
2265       MagickBooleanType
2266         sync;
2267
2268       PixelInfo
2269         pixel,    /* pixel color to assign to distorted image */
2270         invalid;  /* the color to assign when distort result is invalid */
2271
2272       PointInfo
2273         d,
2274         s;  /* transform destination image x,y  to source image x,y */
2275
2276       register ssize_t
2277         i;
2278
2279       register Quantum
2280         *restrict q;
2281
2282       q=QueueCacheViewAuthenticPixels(distort_view,0,j,distort_image->columns,1,
2283         exception);
2284       if (q == (Quantum *) NULL)
2285         {
2286           status=MagickFalse;
2287           continue;
2288         }
2289       pixel=zero;
2290
2291       /* Define constant scaling vectors for Affine Distortions
2292         Other methods are either variable, or use interpolated lookup
2293       */
2294       switch (method)
2295       {
2296         case AffineDistortion:
2297           ScaleFilter( resample_filter[id],
2298             coeff[0], coeff[1],
2299             coeff[3], coeff[4] );
2300           break;
2301         default:
2302           break;
2303       }
2304
2305       /* Initialize default pixel validity
2306       *    negative:         pixel is invalid  output 'matte_color'
2307       *    0.0 to 1.0:       antialiased, mix with resample output
2308       *    1.0 or greater:   use resampled output.
2309       */
2310       validity = 1.0;
2311
2312       GetPixelInfo(distort_image,&invalid);
2313       SetPixelInfoPacket(distort_image,&distort_image->matte_color,&invalid);
2314       if (distort_image->colorspace == CMYKColorspace)
2315         ConvertRGBToCMYK(&invalid);   /* what about other color spaces? */
2316
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 }