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