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

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