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,exception);
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,exception);
1459     }
1460   else
1461     {
1462       ChannelType
1463         channel_mask;
1464
1465       Image
1466         *resize_alpha;
1467
1468       /*
1469         Image has transparency so handle colors and alpha separatly.
1470         Basically we need to separate Virtual-Pixel alpha in the resized
1471         image, so only the actual original images alpha channel is used.
1472
1473         distort alpha channel separately
1474       */
1475       channel_mask=SetPixelChannelMask(tmp_image,AlphaChannel);
1476       (void) SeparateImage(tmp_image);
1477       SetPixelChannelMap(tmp_image,channel_mask);
1478       (void) SetImageAlphaChannel(tmp_image,OpaqueAlphaChannel,exception);
1479       resize_alpha=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1480         MagickTrue,exception),
1481       tmp_image=DestroyImage(tmp_image);
1482       if (resize_alpha == (Image *) NULL)
1483         return((Image *) NULL);
1484
1485       /* distort the actual image containing alpha + VP alpha */
1486       tmp_image=CloneImage(image,0,0,MagickTrue,exception);
1487       if ( tmp_image == (Image *) NULL )
1488         return((Image *) NULL);
1489       (void) SetImageVirtualPixelMethod(tmp_image,
1490         TransparentVirtualPixelMethod);
1491       resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
1492         MagickTrue,exception),
1493       tmp_image=DestroyImage(tmp_image);
1494       if ( resize_image == (Image *) NULL)
1495         {
1496           resize_alpha=DestroyImage(resize_alpha);
1497           return((Image *) NULL);
1498         }
1499       /* replace resize images alpha with the separally distorted alpha */
1500       (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel,exception);
1501       (void) SetImageAlphaChannel(resize_alpha,DeactivateAlphaChannel,exception);
1502       (void) CompositeImage(resize_image,CopyOpacityCompositeOp,resize_alpha,
1503         0,0);
1504       resize_alpha=DestroyImage(resize_alpha);
1505     }
1506   (void) SetImageVirtualPixelMethod(resize_image,vp_save);
1507
1508   /*
1509     Clean up the results of the Distortion
1510   */
1511   crop_area.width=columns;
1512   crop_area.height=rows;
1513   crop_area.x=0;
1514   crop_area.y=0;
1515
1516   tmp_image=resize_image;
1517   resize_image=CropImage(tmp_image,&crop_area,exception);
1518   tmp_image=DestroyImage(tmp_image);
1519
1520   if ( resize_image == (Image *) NULL )
1521     return((Image *) NULL);
1522
1523   return(resize_image);
1524 }
1525 \f
1526 /*
1527 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1528 %                                                                             %
1529 %                                                                             %
1530 %                                                                             %
1531 %   D i s t o r t I m a g e                                                   %
1532 %                                                                             %
1533 %                                                                             %
1534 %                                                                             %
1535 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1536 %
1537 %  DistortImage() distorts an image using various distortion methods, by
1538 %  mapping color lookups of the source image to a new destination image
1539 %  usally of the same size as the source image, unless 'bestfit' is set to
1540 %  true.
1541 %
1542 %  If 'bestfit' is enabled, and distortion allows it, the destination image is
1543 %  adjusted to ensure the whole source 'image' will just fit within the final
1544 %  destination image, which will be sized and offset accordingly.  Also in
1545 %  many cases the virtual offset of the source image will be taken into
1546 %  account in the mapping.
1547 %
1548 %  If the '-verbose' control option has been set print to standard error the
1549 %  equicelent '-fx' formula with coefficients for the function, if practical.
1550 %
1551 %  The format of the DistortImage() method is:
1552 %
1553 %      Image *DistortImage(const Image *image,const DistortImageMethod method,
1554 %        const size_t number_arguments,const double *arguments,
1555 %        MagickBooleanType bestfit, ExceptionInfo *exception)
1556 %
1557 %  A description of each parameter follows:
1558 %
1559 %    o image: the image to be distorted.
1560 %
1561 %    o method: the method of image distortion.
1562 %
1563 %        ArcDistortion always ignores source image offset, and always
1564 %        'bestfit' the destination image with the top left corner offset
1565 %        relative to the polar mapping center.
1566 %
1567 %        Affine, Perspective, and Bilinear, do least squares fitting of the
1568 %        distrotion when more than the minimum number of control point pairs
1569 %        are provided.
1570 %
1571 %        Perspective, and Bilinear, fall back to a Affine distortion when less
1572 %        than 4 control point pairs are provided.  While Affine distortions
1573 %        let you use any number of control point pairs, that is Zero pairs is
1574 %        a No-Op (viewport only) distortion, one pair is a translation and
1575 %        two pairs of control points do a scale-rotate-translate, without any
1576 %        shearing.
1577 %
1578 %    o number_arguments: the number of arguments given.
1579 %
1580 %    o arguments: an array of floating point arguments for this method.
1581 %
1582 %    o bestfit: Attempt to 'bestfit' the size of the resulting image.
1583 %        This also forces the resulting image to be a 'layered' virtual
1584 %        canvas image.  Can be overridden using 'distort:viewport' setting.
1585 %
1586 %    o exception: return any errors or warnings in this structure
1587 %
1588 %  Extra Controls from Image meta-data (artifacts)...
1589 %
1590 %    o "verbose"
1591 %        Output to stderr alternatives, internal coefficents, and FX
1592 %        equivalents for the distortion operation (if feasible).
1593 %        This forms an extra check of the distortion method, and allows users
1594 %        access to the internal constants IM calculates for the distortion.
1595 %
1596 %    o "distort:viewport"
1597 %        Directly set the output image canvas area and offest to use for the
1598 %        resulting image, rather than use the original images canvas, or a
1599 %        calculated 'bestfit' canvas.
1600 %
1601 %    o "distort:scale"
1602 %        Scale the size of the output canvas by this amount to provide a
1603 %        method of Zooming, and for super-sampling the results.
1604 %
1605 %  Other settings that can effect results include
1606 %
1607 %    o 'interpolate' For source image lookups (scale enlargements)
1608 %
1609 %    o 'filter'      Set filter to use for area-resampling (scale shrinking).
1610 %                    Set to 'point' to turn off and use 'interpolate' lookup
1611 %                    instead
1612 %
1613 */
1614
1615 MagickExport Image *DistortImage(const Image *image,DistortImageMethod method,
1616   const size_t number_arguments,const double *arguments,
1617   MagickBooleanType bestfit,ExceptionInfo *exception)
1618 {
1619 #define DistortImageTag  "Distort/Image"
1620
1621   double
1622     *coeff,
1623     output_scaling;
1624
1625   Image
1626     *distort_image;
1627
1628   RectangleInfo
1629     geometry;  /* geometry of the distorted space viewport */
1630
1631   MagickBooleanType
1632     viewport_given;
1633
1634   assert(image != (Image *) NULL);
1635   assert(image->signature == MagickSignature);
1636   if (image->debug != MagickFalse)
1637     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1638   assert(exception != (ExceptionInfo *) NULL);
1639   assert(exception->signature == MagickSignature);
1640
1641
1642   /*
1643     Handle Special Compound Distortions
1644   */
1645   if ( method == ResizeDistortion )
1646     {
1647       if ( number_arguments != 2 )
1648         {
1649           (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
1650                     "InvalidArgument","%s : '%s'","Resize",
1651                     "Invalid number of args: 2 only");
1652           return((Image *) NULL);
1653         }
1654       distort_image=DistortResizeImage(image,(size_t)arguments[0],
1655          (size_t)arguments[1], exception);
1656       return(distort_image);
1657     }
1658
1659   /*
1660     Convert input arguments (usually as control points for reverse mapping)
1661     into mapping coefficients to apply the distortion.
1662
1663     Note that some distortions are mapped to other distortions,
1664     and as such do not require specific code after this point.
1665   */
1666   coeff = GenerateCoefficients(image, &method, number_arguments,
1667       arguments, 0, exception);
1668   if ( coeff == (double *) NULL )
1669     return((Image *) NULL);
1670
1671   /*
1672     Determine the size and offset for a 'bestfit' destination.
1673     Usally the four corners of the source image is enough.
1674   */
1675
1676   /* default output image bounds, when no 'bestfit' is requested */
1677   geometry.width=image->columns;
1678   geometry.height=image->rows;
1679   geometry.x=0;
1680   geometry.y=0;
1681
1682   if ( method == ArcDistortion ) {
1683     bestfit = MagickTrue;  /* always calculate a 'best fit' viewport */
1684   }
1685
1686   /* Work out the 'best fit', (required for ArcDistortion) */
1687   if ( bestfit ) {
1688     PointInfo
1689       s,d,min,max;  /* source, dest coords --mapping--> min, max coords */
1690
1691     MagickBooleanType
1692       fix_bounds = MagickTrue;   /* enlarge bounds for VP handling */
1693
1694     s.x=s.y=min.x=max.x=min.y=max.y=0.0;   /* keep compiler happy */
1695
1696 /* defines to figure out the bounds of the distorted image */
1697 #define InitalBounds(p) \
1698 { \
1699   /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1700   min.x = max.x = p.x; \
1701   min.y = max.y = p.y; \
1702 }
1703 #define ExpandBounds(p) \
1704 { \
1705   /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
1706   min.x = MagickMin(min.x,p.x); \
1707   max.x = MagickMax(max.x,p.x); \
1708   min.y = MagickMin(min.y,p.y); \
1709   max.y = MagickMax(max.y,p.y); \
1710 }
1711
1712     switch (method)
1713     {
1714       case AffineDistortion:
1715       { double inverse[6];
1716         InvertAffineCoefficients(coeff, inverse);
1717         s.x = (double) image->page.x;
1718         s.y = (double) image->page.y;
1719         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1720         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1721         InitalBounds(d);
1722         s.x = (double) image->page.x+image->columns;
1723         s.y = (double) image->page.y;
1724         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1725         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1726         ExpandBounds(d);
1727         s.x = (double) image->page.x;
1728         s.y = (double) image->page.y+image->rows;
1729         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1730         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1731         ExpandBounds(d);
1732         s.x = (double) image->page.x+image->columns;
1733         s.y = (double) image->page.y+image->rows;
1734         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
1735         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
1736         ExpandBounds(d);
1737         break;
1738       }
1739       case PerspectiveDistortion:
1740       { double inverse[8], scale;
1741         InvertPerspectiveCoefficients(coeff, inverse);
1742         s.x = (double) image->page.x;
1743         s.y = (double) image->page.y;
1744         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1745         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1746         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1747         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1748         InitalBounds(d);
1749         s.x = (double) image->page.x+image->columns;
1750         s.y = (double) image->page.y;
1751         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1752         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1753         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1754         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1755         ExpandBounds(d);
1756         s.x = (double) image->page.x;
1757         s.y = (double) image->page.y+image->rows;
1758         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1759         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1760         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1761         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1762         ExpandBounds(d);
1763         s.x = (double) image->page.x+image->columns;
1764         s.y = (double) image->page.y+image->rows;
1765         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
1766         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
1767         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
1768         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
1769         ExpandBounds(d);
1770         break;
1771       }
1772       case ArcDistortion:
1773       { double a, ca, sa;
1774         /* Forward Map Corners */
1775         a = coeff[0]-coeff[1]/2; ca = cos(a); sa = sin(a);
1776         d.x = coeff[2]*ca;
1777         d.y = coeff[2]*sa;
1778         InitalBounds(d);
1779         d.x = (coeff[2]-coeff[3])*ca;
1780         d.y = (coeff[2]-coeff[3])*sa;
1781         ExpandBounds(d);
1782         a = coeff[0]+coeff[1]/2; ca = cos(a); sa = sin(a);
1783         d.x = coeff[2]*ca;
1784         d.y = coeff[2]*sa;
1785         ExpandBounds(d);
1786         d.x = (coeff[2]-coeff[3])*ca;
1787         d.y = (coeff[2]-coeff[3])*sa;
1788         ExpandBounds(d);
1789         /* Orthogonal points along top of arc */
1790         for( a=(double) (ceil((double) ((coeff[0]-coeff[1]/2.0)/MagickPI2))*MagickPI2);
1791                a<(coeff[0]+coeff[1]/2.0); a+=MagickPI2 ) {
1792           ca = cos(a); sa = sin(a);
1793           d.x = coeff[2]*ca;
1794           d.y = coeff[2]*sa;
1795           ExpandBounds(d);
1796         }
1797         /*
1798           Convert the angle_to_width and radius_to_height
1799           to appropriate scaling factors, to allow faster processing
1800           in the mapping function.
1801         */
1802         coeff[1] = (double) (Magick2PI*image->columns/coeff[1]);
1803         coeff[3] = (double)image->rows/coeff[3];
1804         break;
1805       }
1806       case PolarDistortion:
1807       {
1808         if (number_arguments < 2)
1809           coeff[2] = coeff[3] = 0.0;
1810         min.x = coeff[2]-coeff[0];
1811         max.x = coeff[2]+coeff[0];
1812         min.y = coeff[3]-coeff[0];
1813         max.y = coeff[3]+coeff[0];
1814         /* should be about 1.0 if Rmin = 0 */
1815         coeff[7]=(double) geometry.height/(coeff[0]-coeff[1]);
1816         break;
1817       }
1818       case DePolarDistortion:
1819       {
1820         /* direct calculation as it needs to tile correctly
1821          * for reversibility in a DePolar-Polar cycle */
1822         fix_bounds = MagickFalse;
1823         geometry.x = geometry.y = 0;
1824         geometry.height = (size_t) ceil(coeff[0]-coeff[1]);
1825         geometry.width = (size_t)
1826                   ceil((coeff[0]-coeff[1])*(coeff[5]-coeff[4])*0.5);
1827         /* correct scaling factors relative to new size */
1828         coeff[6]=(coeff[5]-coeff[4])/geometry.width; /* changed width */
1829         coeff[7]=(coeff[0]-coeff[1])/geometry.height; /* should be about 1.0 */
1830         break;
1831       }
1832       case Cylinder2PlaneDistortion:
1833       {
1834         /* direct calculation so center of distortion is either a pixel
1835          * center, or pixel edge. This allows for reversibility of the
1836          * distortion */
1837         geometry.x = geometry.y = 0;
1838         geometry.width = (size_t) ceil( 2.0*coeff[1]*tan(coeff[0]/2.0) );
1839         geometry.height = (size_t) ceil( 2.0*coeff[3]/cos(coeff[0]/2.0) );
1840         /* correct center of distortion relative to new size */
1841         coeff[4] = (double) geometry.width/2.0;
1842         coeff[5] = (double) geometry.height/2.0;
1843         fix_bounds = MagickFalse;
1844         break;
1845       }
1846       case Plane2CylinderDistortion:
1847       {
1848         /* direct calculation center is either pixel center, or pixel edge
1849          * so as to allow reversibility of the image distortion */
1850         geometry.x = geometry.y = 0;
1851         geometry.width = (size_t) ceil(coeff[0]*coeff[1]);  /* FOV * radius */
1852         geometry.height = (size_t) (2*coeff[3]);              /* input image height */
1853         /* correct center of distortion relative to new size */
1854         coeff[4] = (double) geometry.width/2.0;
1855         coeff[5] = (double) geometry.height/2.0;
1856         fix_bounds = MagickFalse;
1857         break;
1858       }
1859       case ShepardsDistortion:
1860       case BilinearForwardDistortion:
1861       case BilinearReverseDistortion:
1862 #if 0
1863       case QuadrilateralDistortion:
1864 #endif
1865       case PolynomialDistortion:
1866       case BarrelDistortion:
1867       case BarrelInverseDistortion:
1868       default:
1869         /* no calculated bestfit available for these distortions */
1870         bestfit = MagickFalse;
1871         fix_bounds = MagickFalse;
1872         break;
1873     }
1874
1875     /* Set the output image geometry to calculated 'bestfit'.
1876        Yes this tends to 'over do' the file image size, ON PURPOSE!
1877        Do not do this for DePolar which needs to be exact for virtual tiling.
1878     */
1879     if ( fix_bounds ) {
1880       geometry.x = (ssize_t) floor(min.x-0.5);
1881       geometry.y = (ssize_t) floor(min.y-0.5);
1882       geometry.width=(size_t) ceil(max.x-geometry.x+0.5);
1883       geometry.height=(size_t) ceil(max.y-geometry.y+0.5);
1884     }
1885
1886   }  /* end bestfit destination image calculations */
1887
1888   /* The user provided a 'viewport' expert option which may
1889      overrides some parts of the current output image geometry.
1890      This also overrides its default 'bestfit' setting.
1891   */
1892   { const char *artifact=GetImageArtifact(image,"distort:viewport");
1893     viewport_given = MagickFalse;
1894     if ( artifact != (const char *) NULL ) {
1895       (void) ParseAbsoluteGeometry(artifact,&geometry);
1896       viewport_given = MagickTrue;
1897     }
1898   }
1899
1900   /* Verbose output */
1901   if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
1902     register ssize_t
1903        i;
1904     char image_gen[MaxTextExtent];
1905     const char *lookup;
1906
1907     /* Set destination image size and virtual offset */
1908     if ( bestfit || viewport_given ) {
1909       (void) FormatLocaleString(image_gen, MaxTextExtent,"  -size %.20gx%.20g "
1910         "-page %+.20g%+.20g xc: +insert \\\n",(double) geometry.width,
1911         (double) geometry.height,(double) geometry.x,(double) geometry.y);
1912       lookup="v.p{ xx-v.page.x-.5, yy-v.page.y-.5 }";
1913     }
1914     else {
1915       image_gen[0] = '\0';             /* no destination to generate */
1916       lookup = "p{ xx-page.x-.5, yy-page.y-.5 }"; /* simplify lookup */
1917     }
1918
1919     switch (method) {
1920       case AffineDistortion:
1921       {
1922         double *inverse;
1923
1924         inverse = (double *) AcquireQuantumMemory(6,sizeof(*inverse));
1925         if (inverse == (double *) NULL) {
1926           coeff = (double *) RelinquishMagickMemory(coeff);
1927           (void) ThrowMagickException(exception,GetMagickModule(),
1928                   ResourceLimitError,"MemoryAllocationFailed",
1929                   "%s", "DistortImages");
1930           return((Image *) NULL);
1931         }
1932         InvertAffineCoefficients(coeff, inverse);
1933         CoefficientsToAffineArgs(inverse);
1934         (void) FormatLocaleFile(stderr, "Affine Projection:\n");
1935         (void) FormatLocaleFile(stderr, "  -distort AffineProjection \\\n      '");
1936         for (i=0; i < 5; i++)
1937           (void) FormatLocaleFile(stderr, "%lf,", inverse[i]);
1938         (void) FormatLocaleFile(stderr, "%lf'\n", inverse[5]);
1939         inverse = (double *) RelinquishMagickMemory(inverse);
1940
1941         (void) FormatLocaleFile(stderr, "Affine Distort, FX Equivelent:\n");
1942         (void) FormatLocaleFile(stderr, "%s", image_gen);
1943         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1944         (void) FormatLocaleFile(stderr, "       xx=%+lf*ii %+lf*jj %+lf;\n",
1945             coeff[0], coeff[1], coeff[2]);
1946         (void) FormatLocaleFile(stderr, "       yy=%+lf*ii %+lf*jj %+lf;\n",
1947             coeff[3], coeff[4], coeff[5]);
1948         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
1949
1950         break;
1951       }
1952
1953       case PerspectiveDistortion:
1954       {
1955         double *inverse;
1956
1957         inverse = (double *) AcquireQuantumMemory(8,sizeof(*inverse));
1958         if (inverse == (double *) NULL) {
1959           coeff = (double *) RelinquishMagickMemory(coeff);
1960           (void) ThrowMagickException(exception,GetMagickModule(),
1961                   ResourceLimitError,"MemoryAllocationFailed",
1962                   "%s", "DistortCoefficients");
1963           return((Image *) NULL);
1964         }
1965         InvertPerspectiveCoefficients(coeff, inverse);
1966         (void) FormatLocaleFile(stderr, "Perspective Projection:\n");
1967         (void) FormatLocaleFile(stderr, "  -distort PerspectiveProjection \\\n      '");
1968         for (i=0; i<4; i++)
1969           (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
1970         (void) FormatLocaleFile(stderr, "\n       ");
1971         for (; i<7; i++)
1972           (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
1973         (void) FormatLocaleFile(stderr, "%lf'\n", inverse[7]);
1974         inverse = (double *) RelinquishMagickMemory(inverse);
1975
1976         (void) FormatLocaleFile(stderr, "Perspective Distort, FX Equivelent:\n");
1977         (void) FormatLocaleFile(stderr, "%s", image_gen);
1978         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
1979         (void) FormatLocaleFile(stderr, "       rr=%+lf*ii %+lf*jj + 1;\n",
1980             coeff[6], coeff[7]);
1981         (void) FormatLocaleFile(stderr, "       xx=(%+lf*ii %+lf*jj %+lf)/rr;\n",
1982             coeff[0], coeff[1], coeff[2]);
1983         (void) FormatLocaleFile(stderr, "       yy=(%+lf*ii %+lf*jj %+lf)/rr;\n",
1984             coeff[3], coeff[4], coeff[5]);
1985         (void) FormatLocaleFile(stderr, "       rr%s0 ? %s : blue' \\\n",
1986             coeff[8] < 0 ? "<" : ">", lookup);
1987         break;
1988       }
1989
1990       case BilinearForwardDistortion:
1991         (void) FormatLocaleFile(stderr, "BilinearForward Mapping Equations:\n");
1992         (void) FormatLocaleFile(stderr, "%s", image_gen);
1993         (void) FormatLocaleFile(stderr, "    i = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
1994             coeff[0], coeff[1], coeff[2], coeff[3]);
1995         (void) FormatLocaleFile(stderr, "    j = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
1996             coeff[4], coeff[5], coeff[6], coeff[7]);
1997 #if 0
1998         /* for debugging */
1999         (void) FormatLocaleFile(stderr, "   c8 = %+lf  c9 = 2*a = %+lf;\n",
2000             coeff[8], coeff[9]);
2001 #endif
2002         (void) FormatLocaleFile(stderr, "BilinearForward Distort, FX Equivelent:\n");
2003         (void) FormatLocaleFile(stderr, "%s", image_gen);
2004         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
2005             0.5-coeff[3], 0.5-coeff[7]);
2006         (void) FormatLocaleFile(stderr, "       bb=%lf*ii %+lf*jj %+lf;\n",
2007             coeff[6], -coeff[2], coeff[8]);
2008         /* Handle Special degenerate (non-quadratic) or trapezoidal case */
2009         if ( coeff[9] != 0 ) {
2010           (void) FormatLocaleFile(stderr, "       rt=bb*bb %+lf*(%lf*ii%+lf*jj);\n",
2011               -2*coeff[9],  coeff[4], -coeff[0]);
2012           (void) FormatLocaleFile(stderr, "       yy=( -bb + sqrt(rt) ) / %lf;\n",
2013                coeff[9]);
2014         } else
2015           (void) FormatLocaleFile(stderr, "       yy=(%lf*ii%+lf*jj)/bb;\n",
2016                 -coeff[4], coeff[0]);
2017         (void) FormatLocaleFile(stderr, "       xx=(ii %+lf*yy)/(%lf %+lf*yy);\n",
2018              -coeff[1], coeff[0], coeff[2]);
2019         if ( coeff[9] != 0 )
2020           (void) FormatLocaleFile(stderr, "       (rt < 0 ) ? red : %s'\n", lookup);
2021         else
2022           (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
2023         break;
2024
2025       case BilinearReverseDistortion:
2026 #if 0
2027         (void) FormatLocaleFile(stderr, "Polynomial Projection Distort:\n");
2028         (void) FormatLocaleFile(stderr, "  -distort PolynomialProjection \\\n");
2029         (void) FormatLocaleFile(stderr, "      '1.5, %lf, %lf, %lf, %lf,\n",
2030             coeff[3], coeff[0], coeff[1], coeff[2]);
2031         (void) FormatLocaleFile(stderr, "            %lf, %lf, %lf, %lf'\n",
2032             coeff[7], coeff[4], coeff[5], coeff[6]);
2033 #endif
2034         (void) FormatLocaleFile(stderr, "BilinearReverse Distort, FX Equivelent:\n");
2035         (void) FormatLocaleFile(stderr, "%s", image_gen);
2036         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2037         (void) FormatLocaleFile(stderr, "       xx=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
2038             coeff[0], coeff[1], coeff[2], coeff[3]);
2039         (void) FormatLocaleFile(stderr, "       yy=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
2040             coeff[4], coeff[5], coeff[6], coeff[7]);
2041         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
2042         break;
2043
2044       case PolynomialDistortion:
2045       {
2046         size_t nterms = (size_t) coeff[1];
2047         (void) FormatLocaleFile(stderr, "Polynomial (order %lg, terms %lu), FX Equivelent\n",
2048           coeff[0],(unsigned long) nterms);
2049         (void) FormatLocaleFile(stderr, "%s", image_gen);
2050         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
2051         (void) FormatLocaleFile(stderr, "       xx =");
2052         for (i=0; i<(ssize_t) nterms; i++) {
2053           if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n         ");
2054           (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i],
2055                poly_basis_str(i));
2056         }
2057         (void) FormatLocaleFile(stderr, ";\n       yy =");
2058         for (i=0; i<(ssize_t) nterms; i++) {
2059           if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n         ");
2060           (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i+nterms],
2061                poly_basis_str(i));
2062         }
2063         (void) FormatLocaleFile(stderr, ";\n       %s' \\\n", lookup);
2064         break;
2065       }
2066       case ArcDistortion:
2067       {
2068         (void) FormatLocaleFile(stderr, "Arc Distort, Internal Coefficients:\n");
2069         for ( i=0; i<5; i++ )
2070           (void) FormatLocaleFile(stderr, "  c%.20g = %+lf\n", (double) i, coeff[i]);
2071         (void) FormatLocaleFile(stderr, "Arc Distort, FX Equivelent:\n");
2072         (void) FormatLocaleFile(stderr, "%s", image_gen);
2073         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x; jj=j+page.y;\n");
2074         (void) FormatLocaleFile(stderr, "       xx=(atan2(jj,ii)%+lf)/(2*pi);\n",
2075                                   -coeff[0]);
2076         (void) FormatLocaleFile(stderr, "       xx=xx-round(xx);\n");
2077         (void) FormatLocaleFile(stderr, "       xx=xx*%lf %+lf;\n",
2078                             coeff[1], coeff[4]);
2079         (void) FormatLocaleFile(stderr, "       yy=(%lf - hypot(ii,jj)) * %lf;\n",
2080                             coeff[2], coeff[3]);
2081         (void) FormatLocaleFile(stderr, "       v.p{xx-.5,yy-.5}' \\\n");
2082         break;
2083       }
2084       case PolarDistortion:
2085       {
2086         (void) FormatLocaleFile(stderr, "Polar Distort, Internal Coefficents\n");
2087         for ( i=0; i<8; i++ )
2088           (void) FormatLocaleFile(stderr, "  c%.20g = %+lf\n", (double) i, coeff[i]);
2089         (void) FormatLocaleFile(stderr, "Polar Distort, FX Equivelent:\n");
2090         (void) FormatLocaleFile(stderr, "%s", image_gen);
2091         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
2092                          -coeff[2], -coeff[3]);
2093         (void) FormatLocaleFile(stderr, "       xx=(atan2(ii,jj)%+lf)/(2*pi);\n",
2094                          -(coeff[4]+coeff[5])/2 );
2095         (void) FormatLocaleFile(stderr, "       xx=xx-round(xx);\n");
2096         (void) FormatLocaleFile(stderr, "       xx=xx*2*pi*%lf + v.w/2;\n",
2097                          coeff[6] );
2098         (void) FormatLocaleFile(stderr, "       yy=(hypot(ii,jj)%+lf)*%lf;\n",
2099                          -coeff[1], coeff[7] );
2100         (void) FormatLocaleFile(stderr, "       v.p{xx-.5,yy-.5}' \\\n");
2101         break;
2102       }
2103       case DePolarDistortion:
2104       {
2105         (void) FormatLocaleFile(stderr, "DePolar Distort, Internal Coefficents\n");
2106         for ( i=0; i<8; i++ )
2107           (void) FormatLocaleFile(stderr, "  c%.20g = %+lf\n", (double) i, coeff[i]);
2108         (void) FormatLocaleFile(stderr, "DePolar Distort, FX Equivelent:\n");
2109         (void) FormatLocaleFile(stderr, "%s", image_gen);
2110         (void) FormatLocaleFile(stderr, "  -fx 'aa=(i+.5)*%lf %+lf;\n", coeff[6], -coeff[4] );
2111         (void) FormatLocaleFile(stderr, "       rr=(j+.5)*%lf %+lf;\n", coeff[7], +coeff[1] );
2112         (void) FormatLocaleFile(stderr, "       xx=rr*sin(aa) %+lf;\n", coeff[2] );
2113         (void) FormatLocaleFile(stderr, "       yy=rr*cos(aa) %+lf;\n", coeff[3] );
2114         (void) FormatLocaleFile(stderr, "       v.p{xx-.5,yy-.5}' \\\n");
2115         break;
2116       }
2117       case Cylinder2PlaneDistortion:
2118       {
2119         (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, Internal Coefficents\n");
2120         (void) FormatLocaleFile(stderr, "  cylinder_radius = %+lf\n", coeff[1]);
2121         (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, FX Equivelent:\n");
2122         (void) FormatLocaleFile(stderr, "%s", image_gen);
2123         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
2124                          -coeff[4], -coeff[5]);
2125         (void) FormatLocaleFile(stderr, "       aa=atan(ii/%+lf);\n", coeff[1] );
2126         (void) FormatLocaleFile(stderr, "       xx=%lf*aa%+lf;\n",
2127                          coeff[1], coeff[2] );
2128         (void) FormatLocaleFile(stderr, "       yy=jj*cos(aa)%+lf;\n", coeff[3] );
2129         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
2130         break;
2131       }
2132       case Plane2CylinderDistortion:
2133       {
2134         (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, Internal Coefficents\n");
2135         (void) FormatLocaleFile(stderr, "  cylinder_radius = %+lf\n", coeff[1]);
2136         (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, FX Equivelent:\n");
2137         (void) FormatLocaleFile(stderr, "%s", image_gen);
2138         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
2139                          -coeff[4], -coeff[5]);
2140         (void) FormatLocaleFile(stderr, "       ii=ii/%+lf;\n", coeff[1] );
2141         (void) FormatLocaleFile(stderr, "       xx=%lf*tan(ii)%+lf;\n",
2142                          coeff[1], coeff[2] );
2143         (void) FormatLocaleFile(stderr, "       yy=jj/cos(ii)%+lf;\n",
2144                          coeff[3] );
2145         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
2146         break;
2147       }
2148       case BarrelDistortion:
2149       case BarrelInverseDistortion:
2150       { double xc,yc;
2151         /* NOTE: This does the barrel roll in pixel coords not image coords
2152         ** The internal distortion must do it in image coordinates,
2153         ** so that is what the center coeff (8,9) is given in.
2154         */
2155         xc = ((double)image->columns-1.0)/2.0 + image->page.x;
2156         yc = ((double)image->rows-1.0)/2.0    + image->page.y;
2157         (void) FormatLocaleFile(stderr, "Barrel%s Distort, FX Equivelent:\n",
2158              method == BarrelDistortion ? "" : "Inv");
2159         (void) FormatLocaleFile(stderr, "%s", image_gen);
2160         if ( fabs(coeff[8]-xc-0.5) < 0.1 && fabs(coeff[9]-yc-0.5) < 0.1 )
2161           (void) FormatLocaleFile(stderr, "  -fx 'xc=(w-1)/2;  yc=(h-1)/2;\n");
2162         else
2163           (void) FormatLocaleFile(stderr, "  -fx 'xc=%lf;  yc=%lf;\n",
2164                coeff[8]-0.5, coeff[9]-0.5);
2165         (void) FormatLocaleFile(stderr,
2166              "       ii=i-xc;  jj=j-yc;  rr=hypot(ii,jj);\n");
2167         (void) FormatLocaleFile(stderr, "       ii=ii%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2168              method == BarrelDistortion ? "*" : "/",
2169              coeff[0],coeff[1],coeff[2],coeff[3]);
2170         (void) FormatLocaleFile(stderr, "       jj=jj%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
2171              method == BarrelDistortion ? "*" : "/",
2172              coeff[4],coeff[5],coeff[6],coeff[7]);
2173         (void) FormatLocaleFile(stderr, "       v.p{fx*ii+xc,fy*jj+yc}' \\\n");
2174       }
2175       default:
2176         break;
2177     }
2178   }
2179
2180   /* The user provided a 'scale' expert option will scale the
2181      output image size, by the factor given allowing for super-sampling
2182      of the distorted image space.  Any scaling factors must naturally
2183      be halved as a result.
2184   */
2185   { const char *artifact;
2186     artifact=GetImageArtifact(image,"distort:scale");
2187     output_scaling = 1.0;
2188     if (artifact != (const char *) NULL) {
2189       output_scaling = fabs(InterpretLocaleValue(artifact,(char **) NULL));
2190       geometry.width  *= output_scaling;
2191       geometry.height *= output_scaling;
2192       geometry.x      *= output_scaling;
2193       geometry.y      *= output_scaling;
2194       if ( output_scaling < 0.1 ) {
2195         coeff = (double *) RelinquishMagickMemory(coeff);
2196         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
2197                 "InvalidArgument","%s", "-set option:distort:scale" );
2198         return((Image *) NULL);
2199       }
2200       output_scaling = 1/output_scaling;
2201     }
2202   }
2203 #define ScaleFilter(F,A,B,C,D) \
2204     ScaleResampleFilter( (F), \
2205       output_scaling*(A), output_scaling*(B), \
2206       output_scaling*(C), output_scaling*(D) )
2207
2208   /*
2209     Initialize the distort image attributes.
2210   */
2211   distort_image=CloneImage(image,geometry.width,geometry.height,MagickTrue,
2212     exception);
2213   if (distort_image == (Image *) NULL)
2214     return((Image *) NULL);
2215   /* if image is ColorMapped - change it to DirectClass */
2216   if (SetImageStorageClass(distort_image,DirectClass,exception) == MagickFalse)
2217     {
2218       distort_image=DestroyImage(distort_image);
2219       return((Image *) NULL);
2220     }
2221   distort_image->page.x=geometry.x;
2222   distort_image->page.y=geometry.y;
2223   if (distort_image->background_color.alpha != OpaqueAlpha)
2224     distort_image->matte=MagickTrue;
2225
2226   { /* ----- MAIN CODE -----
2227        Sample the source image to each pixel in the distort image.
2228      */
2229     CacheView
2230       *distort_view;
2231
2232     MagickBooleanType
2233       status;
2234
2235     MagickOffsetType
2236       progress;
2237
2238     PixelInfo
2239       zero;
2240
2241     ResampleFilter
2242       **restrict resample_filter;
2243
2244     ssize_t
2245       j;
2246
2247     status=MagickTrue;
2248     progress=0;
2249     GetPixelInfo(distort_image,&zero);
2250     resample_filter=AcquireResampleFilterThreadSet(image,
2251       UndefinedVirtualPixelMethod,MagickFalse,exception);
2252     distort_view=AcquireCacheView(distort_image);
2253 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2254   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2255 #endif
2256     for (j=0; j < (ssize_t) distort_image->rows; j++)
2257     {
2258       const int
2259         id = GetOpenMPThreadId();
2260
2261       double
2262         validity;  /* how mathematically valid is this the mapping */
2263
2264       MagickBooleanType
2265         sync;
2266
2267       PixelInfo
2268         pixel,    /* pixel color to assign to distorted image */
2269         invalid;  /* the color to assign when distort result is invalid */
2270
2271       PointInfo
2272         d,
2273         s;  /* transform destination image x,y  to source image x,y */
2274
2275       register ssize_t
2276         i;
2277
2278       register Quantum
2279         *restrict q;
2280
2281       q=QueueCacheViewAuthenticPixels(distort_view,0,j,distort_image->columns,1,
2282         exception);
2283       if (q == (const Quantum *) NULL)
2284         {
2285           status=MagickFalse;
2286           continue;
2287         }
2288       pixel=zero;
2289
2290       /* Define constant scaling vectors for Affine Distortions
2291         Other methods are either variable, or use interpolated lookup
2292       */
2293       switch (method)
2294       {
2295         case AffineDistortion:
2296           ScaleFilter( resample_filter[id],
2297             coeff[0], coeff[1],
2298             coeff[3], coeff[4] );
2299           break;
2300         default:
2301           break;
2302       }
2303
2304       /* Initialize default pixel validity
2305       *    negative:         pixel is invalid  output 'matte_color'
2306       *    0.0 to 1.0:       antialiased, mix with resample output
2307       *    1.0 or greater:   use resampled output.
2308       */
2309       validity = 1.0;
2310
2311       GetPixelInfo(distort_image,&invalid);
2312       SetPixelInfoPacket(distort_image,&distort_image->matte_color,&invalid);
2313       if (distort_image->colorspace == CMYKColorspace)
2314         ConvertRGBToCMYK(&invalid);   /* what about other color spaces? */
2315
2316       for (i=0; i < (ssize_t) distort_image->columns; i++)
2317       {
2318         /* map pixel coordinate to distortion space coordinate */
2319         d.x = (double) (geometry.x+i+0.5)*output_scaling;
2320         d.y = (double) (geometry.y+j+0.5)*output_scaling;
2321         s = d;  /* default is a no-op mapping */
2322         switch (method)
2323         {
2324           case AffineDistortion:
2325           {
2326             s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2327             s.y=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2328             /* Affine partial derivitives are constant -- set above */
2329             break;
2330           }
2331           case PerspectiveDistortion:
2332           {
2333             double
2334               p,q,r,abs_r,abs_c6,abs_c7,scale;
2335             /* perspective is a ratio of affines */
2336             p=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
2337             q=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
2338             r=coeff[6]*d.x+coeff[7]*d.y+1.0;
2339             /* Pixel Validity -- is it a 'sky' or 'ground' pixel */
2340             validity = (r*coeff[8] < 0.0) ? 0.0 : 1.0;
2341             /* Determine horizon anti-alias blending */
2342             abs_r = fabs(r)*2;
2343             abs_c6 = fabs(coeff[6]);
2344             abs_c7 = fabs(coeff[7]);
2345             if ( abs_c6 > abs_c7 ) {
2346               if ( abs_r < abs_c6*output_scaling )
2347                 validity = 0.5 - coeff[8]*r/(coeff[6]*output_scaling);
2348             }
2349             else if ( abs_r < abs_c7*output_scaling )
2350               validity = 0.5 - coeff[8]*r/(coeff[7]*output_scaling);
2351             /* Perspective Sampling Point (if valid) */
2352             if ( validity > 0.0 ) {
2353               /* divide by r affine, for perspective scaling */
2354               scale = 1.0/r;
2355               s.x = p*scale;
2356               s.y = q*scale;
2357               /* Perspective Partial Derivatives or Scaling Vectors */
2358               scale *= scale;
2359               ScaleFilter( resample_filter[id],
2360                 (r*coeff[0] - p*coeff[6])*scale,
2361                 (r*coeff[1] - p*coeff[7])*scale,
2362                 (r*coeff[3] - q*coeff[6])*scale,
2363                 (r*coeff[4] - q*coeff[7])*scale );
2364             }
2365             break;
2366           }
2367           case BilinearReverseDistortion:
2368           {
2369             /* Reversed Mapped is just a simple polynomial */
2370             s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2]*d.x*d.y+coeff[3];
2371             s.y=coeff[4]*d.x+coeff[5]*d.y
2372                     +coeff[6]*d.x*d.y+coeff[7];
2373             /* Bilinear partial derivitives of scaling vectors */
2374             ScaleFilter( resample_filter[id],
2375                 coeff[0] + coeff[2]*d.y,
2376                 coeff[1] + coeff[2]*d.x,
2377                 coeff[4] + coeff[6]*d.y,
2378                 coeff[5] + coeff[6]*d.x );
2379             break;
2380           }
2381           case BilinearForwardDistortion:
2382           {
2383             /* Forward mapped needs reversed polynomial equations
2384              * which unfortunatally requires a square root!  */
2385             double b,c;
2386             d.x -= coeff[3];  d.y -= coeff[7];
2387             b = coeff[6]*d.x - coeff[2]*d.y + coeff[8];
2388             c = coeff[4]*d.x - coeff[0]*d.y;
2389
2390             validity = 1.0;
2391             /* Handle Special degenerate (non-quadratic) case
2392              * Currently without horizon anti-alising */
2393             if ( fabs(coeff[9]) < MagickEpsilon )
2394               s.y =  -c/b;
2395             else {
2396               c = b*b - 2*coeff[9]*c;
2397               if ( c < 0.0 )
2398                 validity = 0.0;
2399               else
2400                 s.y = ( -b + sqrt(c) )/coeff[9];
2401             }
2402             if ( validity > 0.0 )
2403               s.x = ( d.x - coeff[1]*s.y) / ( coeff[0] + coeff[2]*s.y );
2404
2405             /* NOTE: the sign of the square root should be -ve for parts
2406                      where the source image becomes 'flipped' or 'mirrored'.
2407                FUTURE: Horizon handling
2408                FUTURE: Scaling factors or Deritives (how?)
2409             */
2410             break;
2411           }
2412 #if 0
2413           case BilinearDistortion:
2414             /* Bilinear mapping of any Quadrilateral to any Quadrilateral */
2415             /* UNDER DEVELOPMENT */
2416             break;
2417 #endif
2418           case PolynomialDistortion:
2419           {
2420             /* multi-ordered polynomial */
2421             register ssize_t
2422               k;
2423
2424             ssize_t
2425               nterms=(ssize_t)coeff[1];
2426
2427             PointInfo
2428               du,dv; /* the du,dv vectors from unit dx,dy -- derivatives */
2429
2430             s.x=s.y=du.x=du.y=dv.x=dv.y=0.0;
2431             for(k=0; k < nterms; k++) {
2432               s.x  += poly_basis_fn(k,d.x,d.y)*coeff[2+k];
2433               du.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k];
2434               du.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k];
2435               s.y  += poly_basis_fn(k,d.x,d.y)*coeff[2+k+nterms];
2436               dv.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k+nterms];
2437               dv.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k+nterms];
2438             }
2439             ScaleFilter( resample_filter[id], du.x,du.y,dv.x,dv.y );
2440             break;
2441           }
2442           case ArcDistortion:
2443           {
2444             /* what is the angle and radius in the destination image */
2445             s.x  = (double) ((atan2(d.y,d.x) - coeff[0])/Magick2PI);
2446             s.x -= MagickRound(s.x);     /* angle */
2447             s.y  = hypot(d.x,d.y);       /* radius */
2448
2449             /* Arc Distortion Partial Scaling Vectors
2450               Are derived by mapping the perpendicular unit vectors
2451               dR  and  dA*R*2PI  rather than trying to map dx and dy
2452               The results is a very simple orthogonal aligned ellipse.
2453             */
2454             if ( s.y > MagickEpsilon )
2455               ScaleFilter( resample_filter[id],
2456                   (double) (coeff[1]/(Magick2PI*s.y)), 0, 0, coeff[3] );
2457             else
2458               ScaleFilter( resample_filter[id],
2459                   distort_image->columns*2, 0, 0, coeff[3] );
2460
2461             /* now scale the angle and radius for source image lookup point */
2462             s.x = s.x*coeff[1] + coeff[4] + image->page.x +0.5;
2463             s.y = (coeff[2] - s.y) * coeff[3] + image->page.y;
2464             break;
2465           }
2466           case PolarDistortion:
2467           { /* 2D Cartesain to Polar View */
2468             d.x -= coeff[2];
2469             d.y -= coeff[3];
2470             s.x  = atan2(d.x,d.y) - (coeff[4]+coeff[5])/2;
2471             s.x /= Magick2PI;
2472             s.x -= MagickRound(s.x);
2473             s.x *= Magick2PI;       /* angle - relative to centerline */
2474             s.y  = hypot(d.x,d.y);  /* radius */
2475
2476             /* Polar Scaling vectors are based on mapping dR and dA vectors
2477                This results in very simple orthogonal scaling vectors
2478             */
2479             if ( s.y > MagickEpsilon )
2480               ScaleFilter( resample_filter[id],
2481                 (double) (coeff[6]/(Magick2PI*s.y)), 0, 0, coeff[7] );
2482             else
2483               ScaleFilter( resample_filter[id],
2484                   distort_image->columns*2, 0, 0, coeff[7] );
2485
2486             /* now finish mapping radius/angle to source x,y coords */
2487             s.x = s.x*coeff[6] + (double)image->columns/2.0 + image->page.x;
2488             s.y = (s.y-coeff[1])*coeff[7] + image->page.y;
2489             break;
2490           }
2491           case DePolarDistortion:
2492           { /* @D Polar to Carteasain  */
2493             /* ignore all destination virtual offsets */
2494             d.x = ((double)i+0.5)*output_scaling*coeff[6]-coeff[4];
2495             d.y = ((double)j+0.5)*output_scaling*coeff[7]+coeff[1];
2496             s.x = d.y*sin(d.x) + coeff[2];
2497             s.y = d.y*cos(d.x) + coeff[3];
2498             /* derivatives are usless - better to use SuperSampling */
2499             break;
2500           }
2501           case Cylinder2PlaneDistortion:
2502           { /* 3D Cylinder to Tangential Plane */
2503             double ax, cx;
2504             /* relative to center of distortion */
2505             d.x -= coeff[4]; d.y -= coeff[5];
2506             d.x /= coeff[1];        /* x' = x/r */
2507             ax=atan(d.x);           /* aa = atan(x/r) = u/r  */
2508             cx=cos(ax);             /* cx = cos(atan(x/r)) = 1/sqrt(x^2+u^2) */
2509             s.x = coeff[1]*ax;      /* u  = r*atan(x/r) */
2510             s.y = d.y*cx;           /* v  = y*cos(u/r) */
2511             /* derivatives... (see personnal notes) */
2512             ScaleFilter( resample_filter[id],
2513                   1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2514 #if 0
2515 if ( i == 0 && j == 0 ) {
2516   fprintf(stderr, "x=%lf  y=%lf  u=%lf  v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2517   fprintf(stderr, "phi = %lf\n", (double)(ax * 180.0/MagickPI) );
2518   fprintf(stderr, "du/dx=%lf  du/dx=%lf  dv/dx=%lf  dv/dy=%lf\n",
2519                 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
2520   fflush(stderr); }
2521 #endif
2522             /* add center of distortion in source */
2523             s.x += coeff[2]; s.y += coeff[3];
2524             break;
2525           }
2526           case Plane2CylinderDistortion:
2527           { /* 3D Cylinder to Tangential Plane */
2528             /* relative to center of distortion */
2529             d.x -= coeff[4]; d.y -= coeff[5];
2530
2531             /* is pixel valid - horizon of a infinite Virtual-Pixel Plane
2532              * (see Anthony Thyssen's personal note) */
2533             validity = (double) (coeff[1]*MagickPI2 - fabs(d.x))/output_scaling + 0.5;
2534
2535             if ( validity > 0.0 ) {
2536               double cx,tx;
2537               d.x /= coeff[1];           /* x'= x/r */
2538               cx = 1/cos(d.x);           /* cx = 1/cos(x/r) */
2539               tx = tan(d.x);             /* tx = tan(x/r) */
2540               s.x = coeff[1]*tx;         /* u = r * tan(x/r) */
2541               s.y = d.y*cx;              /* v = y / cos(x/r) */
2542               /* derivatives...  (see Anthony Thyssen's personal notes) */
2543               ScaleFilter( resample_filter[id],
2544                     cx*cx, 0.0, s.y*cx/coeff[1], cx );
2545 #if 1
2546 /*if ( i == 0 && j == 0 )*/
2547 if ( d.x == 0.5 && d.y == 0.5 ) {
2548   fprintf(stderr, "x=%lf  y=%lf  u=%lf  v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
2549   fprintf(stderr, "radius = %lf  phi = %lf  validity = %lf\n",
2550       coeff[1],  (double)(d.x * 180.0/MagickPI), validity );
2551   fprintf(stderr, "du/dx=%lf  du/dx=%lf  dv/dx=%lf  dv/dy=%lf\n",
2552       cx*cx, 0.0, s.y*cx/coeff[1], cx);
2553   fflush(stderr); }
2554 #endif
2555             }
2556             /* add center of distortion in source */
2557             s.x += coeff[2]; s.y += coeff[3];
2558             break;
2559           }
2560           case BarrelDistortion:
2561           case BarrelInverseDistortion:
2562           { /* Lens Barrel Distionion Correction */
2563             double r,fx,fy,gx,gy;
2564             /* Radial Polynomial Distortion (de-normalized) */
2565             d.x -= coeff[8];
2566             d.y -= coeff[9];
2567             r = sqrt(d.x*d.x+d.y*d.y);
2568             if ( r > MagickEpsilon ) {
2569               fx = ((coeff[0]*r + coeff[1])*r + coeff[2])*r + coeff[3];
2570               fy = ((coeff[4]*r + coeff[5])*r + coeff[6])*r + coeff[7];
2571               gx = ((3*coeff[0]*r + 2*coeff[1])*r + coeff[2])/r;
2572               gy = ((3*coeff[4]*r + 2*coeff[5])*r + coeff[6])/r;
2573               /* adjust functions and scaling for 'inverse' form */
2574               if ( method == BarrelInverseDistortion ) {
2575                 fx = 1/fx;  fy = 1/fy;
2576                 gx *= -fx*fx;  gy *= -fy*fy;
2577               }
2578               /* Set the source pixel to lookup and EWA derivative vectors */
2579               s.x = d.x*fx + coeff[8];
2580               s.y = d.y*fy + coeff[9];
2581               ScaleFilter( resample_filter[id],
2582                   gx*d.x*d.x + fx, gx*d.x*d.y,
2583                   gy*d.x*d.y,      gy*d.y*d.y + fy );
2584             }
2585             else {
2586               /* Special handling to avoid divide by zero when r==0
2587               **
2588               ** The source and destination pixels match in this case
2589               ** which was set at the top of the loop using  s = d;
2590               ** otherwise...   s.x=coeff[8]; s.y=coeff[9];
2591               */
2592               if ( method == BarrelDistortion )
2593                 ScaleFilter( resample_filter[id],
2594                      coeff[3], 0, 0, coeff[7] );
2595               else /* method == BarrelInverseDistortion */
2596                 /* FUTURE, trap for D==0 causing division by zero */
2597                 ScaleFilter( resample_filter[id],
2598                      1.0/coeff[3], 0, 0, 1.0/coeff[7] );
2599             }
2600             break;
2601           }
2602           case ShepardsDistortion:
2603           { /* Shepards Method, or Inverse Weighted Distance for
2604               displacement around the destination image control points
2605               The input arguments are the coefficents to the function.
2606               This is more of a 'displacement' function rather than an
2607               absolute distortion function.
2608             */
2609             size_t
2610               i;
2611             double
2612               denominator;
2613
2614             denominator = s.x = s.y = 0;
2615             for(i=0; i<number_arguments; i+=4) {
2616               double weight =
2617                   ((double)d.x-arguments[i+2])*((double)d.x-arguments[i+2])
2618                 + ((double)d.y-arguments[i+3])*((double)d.y-arguments[i+3]);
2619               if ( weight != 0 )
2620                 weight = 1/weight;
2621               else
2622                 weight = 1;
2623
2624               s.x += (arguments[ i ]-arguments[i+2])*weight;
2625               s.y += (arguments[i+1]-arguments[i+3])*weight;
2626               denominator += weight;
2627             }
2628             s.x /= denominator;
2629             s.y /= denominator;
2630             s.x += d.x;
2631             s.y += d.y;
2632
2633             /* We can not determine derivatives using shepards method
2634                only color interpolatation, not area-resampling */
2635             break;
2636           }
2637           default:
2638             break; /* use the default no-op given above */
2639         }
2640         /* map virtual canvas location back to real image coordinate */
2641         if ( bestfit && method != ArcDistortion ) {
2642           s.x -= image->page.x;
2643           s.y -= image->page.y;
2644         }
2645         s.x -= 0.5;
2646         s.y -= 0.5;
2647
2648         if ( validity <= 0.0 ) {
2649           /* result of distortion is an invalid pixel - don't resample */
2650           SetPixelPixelInfo(distort_image,&invalid,q);
2651         }
2652         else {
2653           /* resample the source image to find its correct color */
2654           (void) ResamplePixelColor(resample_filter[id],s.x,s.y,&pixel);
2655           /* if validity between 0.0 and 1.0 mix result with invalid pixel */
2656           if ( validity < 1.0 ) {
2657             /* Do a blend of sample color and invalid pixel */
2658             /* should this be a 'Blend', or an 'Over' compose */
2659             CompositePixelInfoBlend(&pixel,validity,&invalid,(1.0-validity),
2660               &pixel);
2661           }
2662           SetPixelPixelInfo(distort_image,&pixel,q);
2663         }
2664         q+=GetPixelChannels(distort_image);
2665       }
2666       sync=SyncCacheViewAuthenticPixels(distort_view,exception);
2667       if (sync == MagickFalse)
2668         status=MagickFalse;
2669       if (image->progress_monitor != (MagickProgressMonitor) NULL)
2670         {
2671           MagickBooleanType
2672             proceed;
2673
2674 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2675   #pragma omp critical (MagickCore_DistortImage)
2676 #endif
2677           proceed=SetImageProgress(image,DistortImageTag,progress++,
2678             image->rows);
2679           if (proceed == MagickFalse)
2680             status=MagickFalse;
2681         }
2682     }
2683     distort_view=DestroyCacheView(distort_view);
2684     resample_filter=DestroyResampleFilterThreadSet(resample_filter);
2685
2686     if (status == MagickFalse)
2687       distort_image=DestroyImage(distort_image);
2688   }
2689
2690   /* Arc does not return an offset unless 'bestfit' is in effect
2691      And the user has not provided an overriding 'viewport'.
2692    */
2693   if ( method == ArcDistortion && !bestfit && !viewport_given ) {
2694     distort_image->page.x = 0;
2695     distort_image->page.y = 0;
2696   }
2697   coeff = (double *) RelinquishMagickMemory(coeff);
2698   return(distort_image);
2699 }
2700 \f
2701 /*
2702 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2703 %                                                                             %
2704 %                                                                             %
2705 %                                                                             %
2706 %   S p a r s e C o l o r I m a g e                                           %
2707 %                                                                             %
2708 %                                                                             %
2709 %                                                                             %
2710 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2711 %
2712 %  SparseColorImage(), given a set of coordinates, interpolates the colors
2713 %  found at those coordinates, across the whole image, using various methods.
2714 %
2715 %  The format of the SparseColorImage() method is:
2716 %
2717 %      Image *SparseColorImage(const Image *image,
2718 %        const SparseColorMethod method,const size_t number_arguments,
2719 %        const double *arguments,ExceptionInfo *exception)
2720 %
2721 %  A description of each parameter follows:
2722 %
2723 %    o image: the image to be filled in.
2724 %
2725 %    o method: the method to fill in the gradient between the control points.
2726 %
2727 %        The methods used for SparseColor() are often simular to methods
2728 %        used for DistortImage(), and even share the same code for determination
2729 %        of the function coefficents, though with more dimensions (or resulting
2730 %        values).
2731 %
2732 %    o number_arguments: the number of arguments given.
2733 %
2734 %    o arguments: array of floating point arguments for this method--
2735 %        x,y,color_values-- with color_values given as normalized values.
2736 %
2737 %    o exception: return any errors or warnings in this structure
2738 %
2739 */
2740 MagickExport Image *SparseColorImage(const Image *image,
2741   const SparseColorMethod method,const size_t number_arguments,
2742   const double *arguments,ExceptionInfo *exception)
2743 {
2744 #define SparseColorTag  "Distort/SparseColor"
2745
2746   SparseColorMethod
2747     sparse_method;
2748
2749   double
2750     *coeff;
2751
2752   Image
2753     *sparse_image;
2754
2755   size_t
2756     number_colors;
2757
2758   assert(image != (Image *) NULL);
2759   assert(image->signature == MagickSignature);
2760   if (image->debug != MagickFalse)
2761     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2762   assert(exception != (ExceptionInfo *) NULL);
2763   assert(exception->signature == MagickSignature);
2764
2765   /* Determine number of color values needed per control point */
2766   number_colors=0;
2767   if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2768     number_colors++;
2769   if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2770     number_colors++;
2771   if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2772     number_colors++;
2773   if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2774       (image->colorspace == CMYKColorspace))
2775     number_colors++;
2776   if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2777       (image->matte != MagickFalse))
2778     number_colors++;
2779
2780   /*
2781     Convert input arguments into mapping coefficients, this this case
2782     we are mapping (distorting) colors, rather than coordinates.
2783   */
2784   { DistortImageMethod
2785       distort_method;
2786
2787     distort_method=(DistortImageMethod) method;
2788     if ( distort_method >= SentinelDistortion )
2789       distort_method = ShepardsDistortion; /* Pretend to be Shepards */
2790     coeff = GenerateCoefficients(image, &distort_method, number_arguments,
2791                 arguments, number_colors, exception);
2792     if ( coeff == (double *) NULL )
2793       return((Image *) NULL);
2794     /*
2795       Note some Distort Methods may fall back to other simpler methods,
2796       Currently the only fallback of concern is Bilinear to Affine
2797       (Barycentric), which is alaso sparse_colr method.  This also ensures
2798       correct two and one color Barycentric handling.
2799     */
2800     sparse_method = (SparseColorMethod) distort_method;
2801     if ( distort_method == ShepardsDistortion )
2802       sparse_method = method;   /* return non-distiort methods to normal */
2803   }
2804
2805   /* Verbose output */
2806   if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
2807
2808     switch (sparse_method) {
2809       case BarycentricColorInterpolate:
2810       {
2811         register ssize_t x=0;
2812         (void) FormatLocaleFile(stderr, "Barycentric Sparse Color:\n");
2813         if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2814           (void) FormatLocaleFile(stderr, "  -channel R -fx '%+lf*i %+lf*j %+lf' \\\n",
2815               coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2816         if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2817           (void) FormatLocaleFile(stderr, "  -channel G -fx '%+lf*i %+lf*j %+lf' \\\n",
2818               coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2819         if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2820           (void) FormatLocaleFile(stderr, "  -channel B -fx '%+lf*i %+lf*j %+lf' \\\n",
2821               coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2822         if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2823             (image->colorspace == CMYKColorspace))
2824           (void) FormatLocaleFile(stderr, "  -channel K -fx '%+lf*i %+lf*j %+lf' \\\n",
2825               coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2826         if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2827             (image->matte != MagickFalse))
2828           (void) FormatLocaleFile(stderr, "  -channel A -fx '%+lf*i %+lf*j %+lf' \\\n",
2829               coeff[x], coeff[x+1], coeff[x+2]),x+=3;
2830         break;
2831       }
2832       case BilinearColorInterpolate:
2833       {
2834         register ssize_t x=0;
2835         (void) FormatLocaleFile(stderr, "Bilinear Sparse Color\n");
2836         if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2837           (void) FormatLocaleFile(stderr, "   -channel R -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2838               coeff[ x ], coeff[x+1],
2839               coeff[x+2], coeff[x+3]),x+=4;
2840         if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2841           (void) FormatLocaleFile(stderr, "   -channel G -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2842               coeff[ x ], coeff[x+1],
2843               coeff[x+2], coeff[x+3]),x+=4;
2844         if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2845           (void) FormatLocaleFile(stderr, "   -channel B -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2846               coeff[ x ], coeff[x+1],
2847               coeff[x+2], coeff[x+3]),x+=4;
2848         if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2849             (image->colorspace == CMYKColorspace))
2850           (void) FormatLocaleFile(stderr, "   -channel K -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2851               coeff[ x ], coeff[x+1],
2852               coeff[x+2], coeff[x+3]),x+=4;
2853         if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2854             (image->matte != MagickFalse))
2855           (void) FormatLocaleFile(stderr, "   -channel A -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
2856               coeff[ x ], coeff[x+1],
2857               coeff[x+2], coeff[x+3]),x+=4;
2858         break;
2859       }
2860       default:
2861         /* sparse color method is too complex for FX emulation */
2862         break;
2863     }
2864   }
2865
2866   /* Generate new image for generated interpolated gradient.
2867    * ASIDE: Actually we could have just replaced the colors of the original
2868    * image, but IM Core policy, is if storage class could change then clone
2869    * the image.
2870    */
2871
2872   sparse_image=CloneImage(image,0,0,MagickTrue,exception);
2873   if (sparse_image == (Image *) NULL)
2874     return((Image *) NULL);
2875   if (SetImageStorageClass(sparse_image,DirectClass,exception) == MagickFalse)
2876     { /* if image is ColorMapped - change it to DirectClass */
2877       sparse_image=DestroyImage(sparse_image);
2878       return((Image *) NULL);
2879     }
2880   { /* ----- MAIN CODE ----- */
2881     CacheView
2882       *sparse_view;
2883
2884     MagickBooleanType
2885       status;
2886
2887     MagickOffsetType
2888       progress;
2889
2890     ssize_t
2891       j;
2892
2893     status=MagickTrue;
2894     progress=0;
2895     sparse_view=AcquireCacheView(sparse_image);
2896 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2897   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
2898 #endif
2899     for (j=0; j < (ssize_t) sparse_image->rows; j++)
2900     {
2901       MagickBooleanType
2902         sync;
2903
2904       PixelInfo
2905         pixel;    /* pixel to assign to distorted image */
2906
2907       register ssize_t
2908         i;
2909
2910       register Quantum
2911         *restrict q;
2912
2913       q=GetCacheViewAuthenticPixels(sparse_view,0,j,sparse_image->columns,
2914         1,exception);
2915       if (q == (const Quantum *) NULL)
2916         {
2917           status=MagickFalse;
2918           continue;
2919         }
2920       GetPixelInfo(sparse_image,&pixel);
2921       for (i=0; i < (ssize_t) image->columns; i++)
2922       {
2923         SetPixelInfo(image,q,&pixel);
2924         switch (sparse_method)
2925         {
2926           case BarycentricColorInterpolate:
2927           {
2928             register ssize_t x=0;
2929             if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2930               pixel.red     = coeff[x]*i +coeff[x+1]*j
2931                               +coeff[x+2], x+=3;
2932             if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2933               pixel.green   = coeff[x]*i +coeff[x+1]*j
2934                               +coeff[x+2], x+=3;
2935             if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2936               pixel.blue    = coeff[x]*i +coeff[x+1]*j
2937                               +coeff[x+2], x+=3;
2938             if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2939                 (image->colorspace == CMYKColorspace))
2940               pixel.black   = coeff[x]*i +coeff[x+1]*j
2941                               +coeff[x+2], x+=3;
2942             if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2943                 (image->matte != MagickFalse))
2944               pixel.alpha = coeff[x]*i +coeff[x+1]*j
2945                               +coeff[x+2], x+=3;
2946             break;
2947           }
2948           case BilinearColorInterpolate:
2949           {
2950             register ssize_t x=0;
2951             if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2952               pixel.red     = coeff[x]*i     + coeff[x+1]*j +
2953                               coeff[x+2]*i*j + coeff[x+3], x+=4;
2954             if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2955               pixel.green   = coeff[x]*i     + coeff[x+1]*j +
2956                               coeff[x+2]*i*j + coeff[x+3], x+=4;
2957             if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2958               pixel.blue    = coeff[x]*i     + coeff[x+1]*j +
2959                               coeff[x+2]*i*j + coeff[x+3], x+=4;
2960             if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2961                 (image->colorspace == CMYKColorspace))
2962               pixel.black   = coeff[x]*i     + coeff[x+1]*j +
2963                               coeff[x+2]*i*j + coeff[x+3], x+=4;
2964             if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2965                 (image->matte != MagickFalse))
2966               pixel.alpha = coeff[x]*i     + coeff[x+1]*j +
2967                               coeff[x+2]*i*j + coeff[x+3], x+=4;
2968             break;
2969           }
2970           case InverseColorInterpolate:
2971           case ShepardsColorInterpolate:
2972           { /* Inverse (Squared) Distance weights average (IDW) */
2973             size_t
2974               k;
2975             double
2976               denominator;
2977
2978             if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
2979               pixel.red=0.0;
2980             if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
2981               pixel.green=0.0;
2982             if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
2983               pixel.blue=0.0;
2984             if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
2985                 (image->colorspace == CMYKColorspace))
2986               pixel.black=0.0;
2987             if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
2988                 (image->matte != MagickFalse))
2989               pixel.alpha=0.0;
2990             denominator = 0.0;
2991             for(k=0; k<number_arguments; k+=2+number_colors) {
2992               register ssize_t x=(ssize_t) k+2;
2993               double weight =
2994                   ((double)i-arguments[ k ])*((double)i-arguments[ k ])
2995                 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
2996               if ( method == InverseColorInterpolate )
2997                 weight = sqrt(weight);  /* inverse, not inverse squared */
2998               weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
2999               if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3000                 pixel.red     += arguments[x++]*weight;
3001               if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3002                 pixel.green   += arguments[x++]*weight;
3003               if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3004                 pixel.blue    += arguments[x++]*weight;
3005               if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3006                   (image->colorspace == CMYKColorspace))
3007                 pixel.black   += arguments[x++]*weight;
3008               if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3009                   (image->matte != MagickFalse))
3010                 pixel.alpha += arguments[x++]*weight;
3011               denominator += weight;
3012             }
3013             if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3014               pixel.red/=denominator;
3015             if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3016               pixel.green/=denominator;
3017             if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3018               pixel.blue/=denominator;
3019             if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3020                 (image->colorspace == CMYKColorspace))
3021               pixel.black/=denominator;
3022             if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3023                 (image->matte != MagickFalse))
3024               pixel.alpha/=denominator;
3025             break;
3026           }
3027           case VoronoiColorInterpolate:
3028           default:
3029           { /* Just use the closest control point you can find! */
3030             size_t
3031               k;
3032             double
3033               minimum = MagickHuge;
3034
3035             for(k=0; k<number_arguments; k+=2+number_colors) {
3036               double distance =
3037                   ((double)i-arguments[ k ])*((double)i-arguments[ k ])
3038                 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
3039               if ( distance < minimum ) {
3040                 register ssize_t x=(ssize_t) k+2;
3041                 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3042                   pixel.red=arguments[x++];
3043                 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3044                   pixel.green=arguments[x++];
3045                 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3046                   pixel.blue=arguments[x++];
3047                 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3048                     (image->colorspace == CMYKColorspace))
3049                   pixel.black=arguments[x++];
3050                 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3051                     (image->matte != MagickFalse))
3052                   pixel.alpha=arguments[x++];
3053                 minimum = distance;
3054               }
3055             }
3056             break;
3057           }
3058         }
3059         /* set the color directly back into the source image */
3060         if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
3061           pixel.red*=QuantumRange;
3062         if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
3063           pixel.green*=QuantumRange;
3064         if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
3065           pixel.blue*=QuantumRange;
3066         if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
3067             (image->colorspace == CMYKColorspace))
3068           pixel.black*=QuantumRange;
3069         if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
3070             (image->matte != MagickFalse))
3071           pixel.alpha*=QuantumRange;
3072         SetPixelPixelInfo(sparse_image,&pixel,q);
3073         q+=GetPixelChannels(sparse_image);
3074       }
3075       sync=SyncCacheViewAuthenticPixels(sparse_view,exception);
3076       if (sync == MagickFalse)
3077         status=MagickFalse;
3078       if (image->progress_monitor != (MagickProgressMonitor) NULL)
3079         {
3080           MagickBooleanType
3081             proceed;
3082
3083 #if defined(MAGICKCORE_OPENMP_SUPPORT)
3084   #pragma omp critical (MagickCore_SparseColorImage)
3085 #endif
3086           proceed=SetImageProgress(image,SparseColorTag,progress++,image->rows);
3087           if (proceed == MagickFalse)
3088             status=MagickFalse;
3089         }
3090     }
3091     sparse_view=DestroyCacheView(sparse_view);
3092     if (status == MagickFalse)
3093       sparse_image=DestroyImage(sparse_image);
3094   }
3095   coeff = (double *) RelinquishMagickMemory(coeff);
3096   return(sparse_image);
3097 }