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