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

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