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

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