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