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

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