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

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