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