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