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