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