MagickCore  6.7.5
distort.c
Go to the documentation of this file.
00001 /*
00002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00003 %                                                                             %
00004 %                                                                             %
00005 %                                                                             %
00006 %               DDDD   IIIII  SSSSS  TTTTT   OOO   RRRR   TTTTT               %
00007 %               D   D    I    SS       T    O   O  R   R    T                 %
00008 %               D   D    I     SSS     T    O   O  RRRR     T                 %
00009 %               D   D    I       SS    T    O   O  R R      T                 %
00010 %               DDDD   IIIII  SSSSS    T     OOO   R  R     T                 %
00011 %                                                                             %
00012 %                                                                             %
00013 %                     MagickCore Image Distortion Methods                     %
00014 %                                                                             %
00015 %                              Software Design                                %
00016 %                                John Cristy                                  %
00017 %                              Anthony Thyssen                                %
00018 %                                 June 2007                                   %
00019 %                                                                             %
00020 %                                                                             %
00021 %  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
00022 %  dedicated to making software imaging solutions freely available.           %
00023 %                                                                             %
00024 %  You may not use this file except in compliance with the License.  You may  %
00025 %  obtain a copy of the License at                                            %
00026 %                                                                             %
00027 %    http://www.imagemagick.org/script/license.php                            %
00028 %                                                                             %
00029 %  Unless required by applicable law or agreed to in writing, software        %
00030 %  distributed under the License is distributed on an "AS IS" BASIS,          %
00031 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
00032 %  See the License for the specific language governing permissions and        %
00033 %  limitations under the License.                                             %
00034 %                                                                             %
00035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00036 %
00037 %
00038 */
00039 
00040 /*
00041   Include declarations.
00042 */
00043 #include "MagickCore/studio.h"
00044 #include "MagickCore/artifact.h"
00045 #include "MagickCore/cache.h"
00046 #include "MagickCore/cache-view.h"
00047 #include "MagickCore/colorspace-private.h"
00048 #include "MagickCore/composite-private.h"
00049 #include "MagickCore/distort.h"
00050 #include "MagickCore/exception.h"
00051 #include "MagickCore/exception-private.h"
00052 #include "MagickCore/gem.h"
00053 #include "MagickCore/hashmap.h"
00054 #include "MagickCore/image.h"
00055 #include "MagickCore/list.h"
00056 #include "MagickCore/matrix.h"
00057 #include "MagickCore/matrix-private.h"
00058 #include "MagickCore/memory_.h"
00059 #include "MagickCore/monitor-private.h"
00060 #include "MagickCore/option.h"
00061 #include "MagickCore/pixel.h"
00062 #include "MagickCore/pixel-accessor.h"
00063 #include "MagickCore/resample.h"
00064 #include "MagickCore/resample-private.h"
00065 #include "MagickCore/registry.h"
00066 #include "MagickCore/semaphore.h"
00067 #include "MagickCore/shear.h"
00068 #include "MagickCore/string_.h"
00069 #include "MagickCore/string-private.h"
00070 #include "MagickCore/thread-private.h"
00071 #include "MagickCore/token.h"
00072 #include "MagickCore/transform.h"
00073 
00074 /*
00075   Numerous internal routines for image distortions.
00076 */
00077 static inline double MagickMin(const double x,const double y)
00078 {
00079   return( x < y ? x : y);
00080 }
00081 static inline double MagickMax(const double x,const double y)
00082 {
00083   return( x > y ? x : y);
00084 }
00085 
00086 static inline void AffineArgsToCoefficients(double *affine)
00087 {
00088   /* map  external sx,ry,rx,sy,tx,ty  to  internal c0,c2,c4,c1,c3,c5 */
00089   double tmp[4];  /* note indexes  0 and 5 remain unchanged */
00090   tmp[0]=affine[1]; tmp[1]=affine[2]; tmp[2]=affine[3]; tmp[3]=affine[4];
00091   affine[3]=tmp[0]; affine[1]=tmp[1]; affine[4]=tmp[2]; affine[2]=tmp[3];
00092 }
00093 
00094 static inline void CoefficientsToAffineArgs(double *coeff)
00095 {
00096   /* map  internal c0,c1,c2,c3,c4,c5  to  external sx,ry,rx,sy,tx,ty */
00097   double tmp[4];  /* note indexes 0 and 5 remain unchanged */
00098   tmp[0]=coeff[3]; tmp[1]=coeff[1]; tmp[2]=coeff[4]; tmp[3]=coeff[2];
00099   coeff[1]=tmp[0]; coeff[2]=tmp[1]; coeff[3]=tmp[2]; coeff[4]=tmp[3];
00100 }
00101 static void InvertAffineCoefficients(const double *coeff,double *inverse)
00102 {
00103   /* From "Digital Image Warping" by George Wolberg, page 50 */
00104   double determinant;
00105 
00106   determinant=1.0/(coeff[0]*coeff[4]-coeff[1]*coeff[3]);
00107   inverse[0]=determinant*coeff[4];
00108   inverse[1]=determinant*(-coeff[1]);
00109   inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[2]*coeff[4]);
00110   inverse[3]=determinant*(-coeff[3]);
00111   inverse[4]=determinant*coeff[0];
00112   inverse[5]=determinant*(coeff[2]*coeff[3]-coeff[0]*coeff[5]);
00113 }
00114 
00115 static void InvertPerspectiveCoefficients(const double *coeff,
00116   double *inverse)
00117 {
00118   /* From "Digital Image Warping" by George Wolberg, page 53 */
00119   double determinant;
00120 
00121   determinant=1.0/(coeff[0]*coeff[4]-coeff[3]*coeff[1]);
00122   inverse[0]=determinant*(coeff[4]-coeff[7]*coeff[5]);
00123   inverse[1]=determinant*(coeff[7]*coeff[2]-coeff[1]);
00124   inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[4]*coeff[2]);
00125   inverse[3]=determinant*(coeff[6]*coeff[5]-coeff[3]);
00126   inverse[4]=determinant*(coeff[0]-coeff[6]*coeff[2]);
00127   inverse[5]=determinant*(coeff[3]*coeff[2]-coeff[0]*coeff[5]);
00128   inverse[6]=determinant*(coeff[3]*coeff[7]-coeff[6]*coeff[4]);
00129   inverse[7]=determinant*(coeff[6]*coeff[1]-coeff[0]*coeff[7]);
00130 }
00131 
00132 static inline double MagickRound(double x)
00133 {
00134   /*
00135     Round the fraction to nearest integer.
00136   */
00137   if (x >= 0.0)
00138     return((double) ((ssize_t) (x+0.5)));
00139   return((double) ((ssize_t) (x-0.5)));
00140 }
00141 
00142 /*
00143  * Polynomial Term Defining Functions
00144  *
00145  * Order must either be an integer, or 1.5 to produce
00146  * the 2 number_valuesal polynomial function...
00147  *    affine     1   (3)      u = c0 + c1*x + c2*y
00148  *    bilinear   1.5 (4)      u = '' + c3*x*y
00149  *    quadratic  2   (6)      u = '' + c4*x*x + c5*y*y
00150  *    cubic      3   (10)     u = '' + c6*x^3 + c7*x*x*y + c8*x*y*y + c9*y^3
00151  *    quartic    4   (15)     u = '' + c10*x^4 + ... + c14*y^4
00152  *    quintic    5   (21)     u = '' + c15*x^5 + ... + c20*y^5
00153  * number in parenthesis minimum number of points needed.
00154  * Anything beyond quintic, has not been implemented until
00155  * a more automated way of determining terms is found.
00156 
00157  * Note the slight re-ordering of the terms for a quadratic polynomial
00158  * which is to allow the use of a bi-linear (order=1.5) polynomial.
00159  * All the later polynomials are ordered simply from x^N to y^N
00160  */
00161 static size_t poly_number_terms(double order)
00162 {
00163  /* Return the number of terms for a 2d polynomial */
00164   if ( order < 1 || order > 5 ||
00165        ( order != floor(order) && (order-1.5) > MagickEpsilon) )
00166     return 0; /* invalid polynomial order */
00167   return((size_t) floor((order+1)*(order+2)/2));
00168 }
00169 
00170 static double poly_basis_fn(ssize_t n, double x, double y)
00171 {
00172   /* Return the result for this polynomial term */
00173   switch(n) {
00174     case  0:  return( 1.0 ); /* constant */
00175     case  1:  return(  x  );
00176     case  2:  return(  y  ); /* affine          order = 1   terms = 3 */
00177     case  3:  return( x*y ); /* bilinear        order = 1.5 terms = 4 */
00178     case  4:  return( x*x );
00179     case  5:  return( y*y ); /* quadratic       order = 2   terms = 6 */
00180     case  6:  return( x*x*x );
00181     case  7:  return( x*x*y );
00182     case  8:  return( x*y*y );
00183     case  9:  return( y*y*y ); /* cubic         order = 3   terms = 10 */
00184     case 10:  return( x*x*x*x );
00185     case 11:  return( x*x*x*y );
00186     case 12:  return( x*x*y*y );
00187     case 13:  return( x*y*y*y );
00188     case 14:  return( y*y*y*y ); /* quartic     order = 4   terms = 15 */
00189     case 15:  return( x*x*x*x*x );
00190     case 16:  return( x*x*x*x*y );
00191     case 17:  return( x*x*x*y*y );
00192     case 18:  return( x*x*y*y*y );
00193     case 19:  return( x*y*y*y*y );
00194     case 20:  return( y*y*y*y*y ); /* quintic   order = 5   terms = 21 */
00195   }
00196   return( 0 ); /* should never happen */
00197 }
00198 static const char *poly_basis_str(ssize_t n)
00199 {
00200   /* return the result for this polynomial term */
00201   switch(n) {
00202     case  0:  return(""); /* constant */
00203     case  1:  return("*ii");
00204     case  2:  return("*jj"); /* affine                order = 1   terms = 3 */
00205     case  3:  return("*ii*jj"); /* bilinear           order = 1.5 terms = 4 */
00206     case  4:  return("*ii*ii");
00207     case  5:  return("*jj*jj"); /* quadratic          order = 2   terms = 6 */
00208     case  6:  return("*ii*ii*ii");
00209     case  7:  return("*ii*ii*jj");
00210     case  8:  return("*ii*jj*jj");
00211     case  9:  return("*jj*jj*jj"); /* cubic           order = 3   terms = 10 */
00212     case 10:  return("*ii*ii*ii*ii");
00213     case 11:  return("*ii*ii*ii*jj");
00214     case 12:  return("*ii*ii*jj*jj");
00215     case 13:  return("*ii*jj*jj*jj");
00216     case 14:  return("*jj*jj*jj*jj"); /* quartic      order = 4   terms = 15 */
00217     case 15:  return("*ii*ii*ii*ii*ii");
00218     case 16:  return("*ii*ii*ii*ii*jj");
00219     case 17:  return("*ii*ii*ii*jj*jj");
00220     case 18:  return("*ii*ii*jj*jj*jj");
00221     case 19:  return("*ii*jj*jj*jj*jj");
00222     case 20:  return("*jj*jj*jj*jj*jj"); /* quintic   order = 5   terms = 21 */
00223   }
00224   return( "UNKNOWN" ); /* should never happen */
00225 }
00226 static double poly_basis_dx(ssize_t n, double x, double y)
00227 {
00228   /* polynomial term for x derivative */
00229   switch(n) {
00230     case  0:  return( 0.0 ); /* constant */
00231     case  1:  return( 1.0 );
00232     case  2:  return( 0.0 ); /* affine      order = 1   terms = 3 */
00233     case  3:  return(  y  ); /* bilinear    order = 1.5 terms = 4 */
00234     case  4:  return(  x  );
00235     case  5:  return( 0.0 ); /* quadratic   order = 2   terms = 6 */
00236     case  6:  return( x*x );
00237     case  7:  return( x*y );
00238     case  8:  return( y*y );
00239     case  9:  return( 0.0 ); /* cubic       order = 3   terms = 10 */
00240     case 10:  return( x*x*x );
00241     case 11:  return( x*x*y );
00242     case 12:  return( x*y*y );
00243     case 13:  return( y*y*y );
00244     case 14:  return( 0.0 ); /* quartic     order = 4   terms = 15 */
00245     case 15:  return( x*x*x*x );
00246     case 16:  return( x*x*x*y );
00247     case 17:  return( x*x*y*y );
00248     case 18:  return( x*y*y*y );
00249     case 19:  return( y*y*y*y );
00250     case 20:  return( 0.0 ); /* quintic     order = 5   terms = 21 */
00251   }
00252   return( 0.0 ); /* should never happen */
00253 }
00254 static double poly_basis_dy(ssize_t n, double x, double y)
00255 {
00256   /* polynomial term for y derivative */
00257   switch(n) {
00258     case  0:  return( 0.0 ); /* constant */
00259     case  1:  return( 0.0 );
00260     case  2:  return( 1.0 ); /* affine      order = 1   terms = 3 */
00261     case  3:  return(  x  ); /* bilinear    order = 1.5 terms = 4 */
00262     case  4:  return( 0.0 );
00263     case  5:  return(  y  ); /* quadratic   order = 2   terms = 6 */
00264     default:  return( poly_basis_dx(n-1,x,y) ); /* weird but true */
00265   }
00266   /* NOTE: the only reason that last is not true for 'quadratic'
00267      is due to the re-arrangement of terms to allow for 'bilinear'
00268   */
00269 }
00270 
00271 /*
00272 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00273 %                                                                             %
00274 %                                                                             %
00275 %                                                                             %
00276 %     A f f i n e T r a n s f o r m I m a g e                                 %
00277 %                                                                             %
00278 %                                                                             %
00279 %                                                                             %
00280 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00281 %
00282 %  AffineTransformImage() transforms an image as dictated by the affine matrix.
00283 %  It allocates the memory necessary for the new Image structure and returns
00284 %  a pointer to the new image.
00285 %
00286 %  The format of the AffineTransformImage method is:
00287 %
00288 %      Image *AffineTransformImage(const Image *image,
00289 %        AffineMatrix *affine_matrix,ExceptionInfo *exception)
00290 %
00291 %  A description of each parameter follows:
00292 %
00293 %    o image: the image.
00294 %
00295 %    o affine_matrix: the affine matrix.
00296 %
00297 %    o exception: return any errors or warnings in this structure.
00298 %
00299 */
00300 MagickExport Image *AffineTransformImage(const Image *image,
00301   const AffineMatrix *affine_matrix,ExceptionInfo *exception)
00302 {
00303   double
00304     distort[6];
00305 
00306   Image
00307     *deskew_image;
00308 
00309   /*
00310     Affine transform image.
00311   */
00312   assert(image->signature == MagickSignature);
00313   if (image->debug != MagickFalse)
00314     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00315   assert(affine_matrix != (AffineMatrix *) NULL);
00316   assert(exception != (ExceptionInfo *) NULL);
00317   assert(exception->signature == MagickSignature);
00318   distort[0]=affine_matrix->sx;
00319   distort[1]=affine_matrix->rx;
00320   distort[2]=affine_matrix->ry;
00321   distort[3]=affine_matrix->sy;
00322   distort[4]=affine_matrix->tx;
00323   distort[5]=affine_matrix->ty;
00324   deskew_image=DistortImage(image,AffineProjectionDistortion,6,distort,
00325     MagickTrue,exception);
00326   return(deskew_image);
00327 }
00328 
00329 /*
00330 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00331 %                                                                             %
00332 %                                                                             %
00333 %                                                                             %
00334 +   G e n e r a t e C o e f f i c i e n t s                                   %
00335 %                                                                             %
00336 %                                                                             %
00337 %                                                                             %
00338 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00339 %
00340 %  GenerateCoefficients() takes user provided input arguments and generates
00341 %  the coefficients, needed to apply the specific distortion for either
00342 %  distorting images (generally using control points) or generating a color
00343 %  gradient from sparsely separated color points.
00344 %
00345 %  The format of the GenerateCoefficients() method is:
00346 %
00347 %    Image *GenerateCoefficients(const Image *image,DistortImageMethod method,
00348 %        const size_t number_arguments,const double *arguments,
00349 %        size_t number_values, ExceptionInfo *exception)
00350 %
00351 %  A description of each parameter follows:
00352 %
00353 %    o image: the image to be distorted.
00354 %
00355 %    o method: the method of image distortion/ sparse gradient
00356 %
00357 %    o number_arguments: the number of arguments given.
00358 %
00359 %    o arguments: the arguments for this distortion method.
00360 %
00361 %    o number_values: the style and format of given control points, (caller type)
00362 %         0: 2 dimensional mapping of control points (Distort)
00363 %            Format:  u,v,x,y  where u,v is the 'source' of the
00364 %            the color to be plotted, for DistortImage()
00365 %         N: Interpolation of control points with N values (usally r,g,b)
00366 %            Format: x,y,r,g,b    mapping x,y to color values r,g,b
00367 %            IN future, variable number of values may be given (1 to N)
00368 %
00369 %    o exception: return any errors or warnings in this structure
00370 %
00371 %  Note that the returned array of double values must be freed by the
00372 %  calling method using RelinquishMagickMemory().  This however may change in
00373 %  the future to require a more 'method' specific method.
00374 %
00375 %  Because of this this method should not be classed as stable or used
00376 %  outside other MagickCore library methods.
00377 */
00378 
00379 static double *GenerateCoefficients(const Image *image,
00380   DistortImageMethod *method,const size_t number_arguments,
00381   const double *arguments,size_t number_values,ExceptionInfo *exception)
00382 {
00383   double
00384     *coeff;
00385 
00386   register size_t
00387     i;
00388 
00389   size_t
00390     number_coeff, /* number of coefficients to return (array size) */
00391     cp_size,      /* number floating point numbers per control point */
00392     cp_x,cp_y,    /* the x,y indexes for control point */
00393     cp_values;    /* index of values for this control point */
00394     /* number_values   Number of values given per control point */
00395 
00396   if ( number_values == 0 ) {
00397     /* Image distortion using control points (or other distortion)
00398        That is generate a mapping so that   x,y->u,v   given  u,v,x,y
00399     */
00400     number_values = 2;   /* special case: two values of u,v */
00401     cp_values = 0;       /* the values i,j are BEFORE the destination CP x,y */
00402     cp_x = 2;            /* location of x,y in input control values */
00403     cp_y = 3;
00404     /* NOTE: cp_values, also used for later 'reverse map distort' tests */
00405   }
00406   else {
00407     cp_x = 0;            /* location of x,y in input control values */
00408     cp_y = 1;
00409     cp_values = 2;       /* and the other values are after x,y */
00410     /* Typically in this case the values are R,G,B color values */
00411   }
00412   cp_size = number_values+2; /* each CP defintion involves this many numbers */
00413 
00414   /* If not enough control point pairs are found for specific distortions
00415      fall back to Affine distortion (allowing 0 to 3 point pairs)
00416   */
00417   if ( number_arguments < 4*cp_size &&
00418        (  *method == BilinearForwardDistortion
00419        || *method == BilinearReverseDistortion
00420        || *method == PerspectiveDistortion
00421        ) )
00422     *method = AffineDistortion;
00423 
00424   number_coeff=0;
00425   switch (*method) {
00426     case AffineDistortion:
00427     /* also BarycentricColorInterpolate: */
00428       number_coeff=3*number_values;
00429       break;
00430     case PolynomialDistortion:
00431       /* number of coefficents depend on the given polynomal 'order' */
00432       if ( number_arguments <= 1 && (number_arguments-1)%cp_size != 0)
00433       {
00434         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00435                    "InvalidArgument","%s : '%s'","Polynomial",
00436                    "Invalid number of args: order [CPs]...");
00437         return((double *) NULL);
00438       }
00439       i = poly_number_terms(arguments[0]);
00440       number_coeff = 2 + i*number_values;
00441       if ( i == 0 ) {
00442         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00443                    "InvalidArgument","%s : '%s'","Polynomial",
00444                    "Invalid order, should be interger 1 to 5, or 1.5");
00445         return((double *) NULL);
00446       }
00447       if ( number_arguments < 1+i*cp_size ) {
00448         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00449                "InvalidArgument", "%s : 'require at least %.20g CPs'",
00450                "Polynomial", (double) i);
00451         return((double *) NULL);
00452       }
00453       break;
00454     case BilinearReverseDistortion:
00455       number_coeff=4*number_values;
00456       break;
00457     /*
00458       The rest are constants as they are only used for image distorts
00459     */
00460     case BilinearForwardDistortion:
00461       number_coeff=10; /* 2*4 coeff plus 2 constants */
00462       cp_x = 0;        /* Reverse src/dest coords for forward mapping */
00463       cp_y = 1;
00464       cp_values = 2;
00465       break;
00466 #if 0
00467     case QuadraterialDistortion:
00468       number_coeff=19; /* BilinearForward + BilinearReverse */
00469 #endif
00470       break;
00471     case ShepardsDistortion:
00472       number_coeff=1;  /* not used, but provide some type of return */
00473       break;
00474     case ArcDistortion:
00475       number_coeff=5;
00476       break;
00477     case ScaleRotateTranslateDistortion:
00478     case AffineProjectionDistortion:
00479     case Plane2CylinderDistortion:
00480     case Cylinder2PlaneDistortion:
00481       number_coeff=6;
00482       break;
00483     case PolarDistortion:
00484     case DePolarDistortion:
00485       number_coeff=8;
00486       break;
00487     case PerspectiveDistortion:
00488     case PerspectiveProjectionDistortion:
00489       number_coeff=9;
00490       break;
00491     case BarrelDistortion:
00492     case BarrelInverseDistortion:
00493       number_coeff=10;
00494       break;
00495     default:
00496       assert(! "Unknown Method Given"); /* just fail assertion */
00497   }
00498 
00499   /* allocate the array of coefficients needed */
00500   coeff = (double *) AcquireQuantumMemory(number_coeff,sizeof(*coeff));
00501   if (coeff == (double *) NULL) {
00502     (void) ThrowMagickException(exception,GetMagickModule(),
00503                   ResourceLimitError,"MemoryAllocationFailed",
00504                   "%s", "GenerateCoefficients");
00505     return((double *) NULL);
00506   }
00507 
00508   /* zero out coefficients array */
00509   for (i=0; i < number_coeff; i++)
00510     coeff[i] = 0.0;
00511 
00512   switch (*method)
00513   {
00514     case AffineDistortion:
00515     {
00516       /* Affine Distortion
00517            v =  c0*x + c1*y + c2
00518          for each 'value' given
00519 
00520          Input Arguments are sets of control points...
00521          For Distort Images    u,v, x,y  ...
00522          For Sparse Gradients  x,y, r,g,b  ...
00523       */
00524       if ( number_arguments%cp_size != 0 ||
00525            number_arguments < cp_size ) {
00526         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00527                "InvalidArgument", "%s : 'require at least %.20g CPs'",
00528                "Affine", 1.0);
00529         coeff=(double *) RelinquishMagickMemory(coeff);
00530         return((double *) NULL);
00531       }
00532       /* handle special cases of not enough arguments */
00533       if ( number_arguments == cp_size ) {
00534         /* Only 1 CP Set Given */
00535         if ( cp_values == 0 ) {
00536           /* image distortion - translate the image */
00537           coeff[0] = 1.0;
00538           coeff[2] = arguments[0] - arguments[2];
00539           coeff[4] = 1.0;
00540           coeff[5] = arguments[1] - arguments[3];
00541         }
00542         else {
00543           /* sparse gradient - use the values directly */
00544           for (i=0; i<number_values; i++)
00545             coeff[i*3+2] = arguments[cp_values+i];
00546         }
00547       }
00548       else {
00549         /* 2 or more points (usally 3) given.
00550            Solve a least squares simultaneous equation for coefficients.
00551         */
00552         double
00553           **matrix,
00554           **vectors,
00555           terms[3];
00556 
00557         MagickBooleanType
00558           status;
00559 
00560         /* create matrix, and a fake vectors matrix */
00561         matrix = AcquireMagickMatrix(3UL,3UL);
00562         vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
00563         if (matrix == (double **) NULL || vectors == (double **) NULL)
00564         {
00565           matrix  = RelinquishMagickMatrix(matrix, 3UL);
00566           vectors = (double **) RelinquishMagickMemory(vectors);
00567           coeff   = (double *) RelinquishMagickMemory(coeff);
00568           (void) ThrowMagickException(exception,GetMagickModule(),
00569                   ResourceLimitError,"MemoryAllocationFailed",
00570                   "%s", "DistortCoefficients");
00571           return((double *) NULL);
00572         }
00573         /* fake a number_values x3 vectors matrix from coefficients array */
00574         for (i=0; i < number_values; i++)
00575           vectors[i] = &(coeff[i*3]);
00576         /* Add given control point pairs for least squares solving */
00577         for (i=0; i < number_arguments; i+=cp_size) {
00578           terms[0] = arguments[i+cp_x];  /* x */
00579           terms[1] = arguments[i+cp_y];  /* y */
00580           terms[2] = 1;                  /* 1 */
00581           LeastSquaresAddTerms(matrix,vectors,terms,
00582                    &(arguments[i+cp_values]),3UL,number_values);
00583         }
00584         if ( number_arguments == 2*cp_size ) {
00585           /* Only two pairs were given, but we need 3 to solve the affine.
00586              Fake extra coordinates by rotating p1 around p0 by 90 degrees.
00587                x2 = x0 - (y1-y0)   y2 = y0 + (x1-x0)
00588            */
00589           terms[0] = arguments[cp_x]
00590                    - ( arguments[cp_size+cp_y] - arguments[cp_y] ); /* x2 */
00591           terms[1] = arguments[cp_y] +
00592                    + ( arguments[cp_size+cp_x] - arguments[cp_x] ); /* y2 */
00593           terms[2] = 1;                                             /* 1 */
00594           if ( cp_values == 0 ) {
00595             /* Image Distortion - rotate the u,v coordients too */
00596             double
00597               uv2[2];
00598             uv2[0] = arguments[0] - arguments[5] + arguments[1];   /* u2 */
00599             uv2[1] = arguments[1] + arguments[4] - arguments[0];   /* v2 */
00600             LeastSquaresAddTerms(matrix,vectors,terms,uv2,3UL,2UL);
00601           }
00602           else {
00603             /* Sparse Gradient - use values of p0 for linear gradient */
00604             LeastSquaresAddTerms(matrix,vectors,terms,
00605                   &(arguments[cp_values]),3UL,number_values);
00606           }
00607         }
00608         /* Solve for LeastSquares Coefficients */
00609         status=GaussJordanElimination(matrix,vectors,3UL,number_values);
00610         matrix = RelinquishMagickMatrix(matrix, 3UL);
00611         vectors = (double **) RelinquishMagickMemory(vectors);
00612         if ( status == MagickFalse ) {
00613           coeff = (double *) RelinquishMagickMemory(coeff);
00614           (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00615               "InvalidArgument","%s : 'Unsolvable Matrix'",
00616               CommandOptionToMnemonic(MagickDistortOptions, *method) );
00617           return((double *) NULL);
00618         }
00619       }
00620       return(coeff);
00621     }
00622     case AffineProjectionDistortion:
00623     {
00624       /*
00625         Arguments: Affine Matrix (forward mapping)
00626         Arguments  sx, rx, ry, sy, tx, ty
00627         Where      u = sx*x + ry*y + tx
00628                    v = rx*x + sy*y + ty
00629 
00630         Returns coefficients (in there inverse form) ordered as...
00631              sx ry tx  rx sy ty
00632 
00633         AffineProjection Distortion Notes...
00634            + Will only work with a 2 number_values for Image Distortion
00635            + Can not be used for generating a sparse gradient (interpolation)
00636       */
00637       double inverse[8];
00638       if (number_arguments != 6) {
00639         coeff = (double *) RelinquishMagickMemory(coeff);
00640         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00641               "InvalidArgument","%s : 'Needs 6 coeff values'",
00642               CommandOptionToMnemonic(MagickDistortOptions, *method) );
00643         return((double *) NULL);
00644       }
00645       /* FUTURE: trap test for sx*sy-rx*ry == 0 (determinant = 0, no inverse) */
00646       for(i=0; i<6UL; i++ )
00647         inverse[i] = arguments[i];
00648       AffineArgsToCoefficients(inverse); /* map into coefficents */
00649       InvertAffineCoefficients(inverse, coeff); /* invert */
00650       *method = AffineDistortion;
00651 
00652       return(coeff);
00653     }
00654     case ScaleRotateTranslateDistortion:
00655     {
00656       /* Scale, Rotate and Translate Distortion
00657          An alternative Affine Distortion
00658          Argument options, by number of arguments given:
00659            7: x,y, sx,sy, a, nx,ny
00660            6: x,y,   s,   a, nx,ny
00661            5: x,y, sx,sy, a
00662            4: x,y,   s,   a
00663            3: x,y,        a
00664            2:        s,   a
00665            1:             a
00666          Where actions are (in order of application)
00667             x,y     'center' of transforms     (default = image center)
00668             sx,sy   scale image by this amount (default = 1)
00669             a       angle of rotation          (argument required)
00670             nx,ny   move 'center' here         (default = x,y or no movement)
00671          And convert to affine mapping coefficients
00672 
00673          ScaleRotateTranslate Distortion Notes...
00674            + Does not use a set of CPs in any normal way
00675            + Will only work with a 2 number_valuesal Image Distortion
00676            + Cannot be used for generating a sparse gradient (interpolation)
00677       */
00678       double
00679         cosine, sine,
00680         x,y,sx,sy,a,nx,ny;
00681 
00682       /* set default center, and default scale */
00683       x = nx = (double)(image->columns)/2.0 + (double)image->page.x;
00684       y = ny = (double)(image->rows)/2.0    + (double)image->page.y;
00685       sx = sy = 1.0;
00686       switch ( number_arguments ) {
00687       case 0:
00688         coeff = (double *) RelinquishMagickMemory(coeff);
00689         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00690               "InvalidArgument","%s : 'Needs at least 1 argument'",
00691               CommandOptionToMnemonic(MagickDistortOptions, *method) );
00692         return((double *) NULL);
00693       case 1:
00694         a = arguments[0];
00695         break;
00696       case 2:
00697         sx = sy = arguments[0];
00698         a = arguments[1];
00699         break;
00700       default:
00701         x = nx = arguments[0];
00702         y = ny = arguments[1];
00703         switch ( number_arguments ) {
00704         case 3:
00705           a = arguments[2];
00706           break;
00707         case 4:
00708           sx = sy = arguments[2];
00709           a = arguments[3];
00710           break;
00711         case 5:
00712           sx = arguments[2];
00713           sy = arguments[3];
00714           a = arguments[4];
00715           break;
00716         case 6:
00717           sx = sy = arguments[2];
00718           a = arguments[3];
00719           nx = arguments[4];
00720           ny = arguments[5];
00721           break;
00722         case 7:
00723           sx = arguments[2];
00724           sy = arguments[3];
00725           a = arguments[4];
00726           nx = arguments[5];
00727           ny = arguments[6];
00728           break;
00729         default:
00730           coeff = (double *) RelinquishMagickMemory(coeff);
00731           (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00732               "InvalidArgument","%s : 'Too Many Arguments (7 or less)'",
00733               CommandOptionToMnemonic(MagickDistortOptions, *method) );
00734           return((double *) NULL);
00735         }
00736         break;
00737       }
00738       /* Trap if sx or sy == 0 -- image is scaled out of existance! */
00739       if ( fabs(sx) < MagickEpsilon || fabs(sy) < MagickEpsilon ) {
00740         coeff = (double *) RelinquishMagickMemory(coeff);
00741         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00742               "InvalidArgument","%s : 'Zero Scale Given'",
00743               CommandOptionToMnemonic(MagickDistortOptions, *method) );
00744         return((double *) NULL);
00745       }
00746       /* Save the given arguments as an affine distortion */
00747       a=DegreesToRadians(a); cosine=cos(a); sine=sin(a);
00748 
00749       *method = AffineDistortion;
00750       coeff[0]=cosine/sx;
00751       coeff[1]=sine/sx;
00752       coeff[2]=x-nx*coeff[0]-ny*coeff[1];
00753       coeff[3]=(-sine)/sy;
00754       coeff[4]=cosine/sy;
00755       coeff[5]=y-nx*coeff[3]-ny*coeff[4];
00756       return(coeff);
00757     }
00758     case PerspectiveDistortion:
00759     { /*
00760          Perspective Distortion (a ratio of affine distortions)
00761 
00762                 p(x,y)    c0*x + c1*y + c2
00763             u = ------ = ------------------
00764                 r(x,y)    c6*x + c7*y + 1
00765 
00766                 q(x,y)    c3*x + c4*y + c5
00767             v = ------ = ------------------
00768                 r(x,y)    c6*x + c7*y + 1
00769 
00770            c8 = Sign of 'r', or the denominator affine, for the actual image.
00771                 This determines what part of the distorted image is 'ground'
00772                 side of the horizon, the other part is 'sky' or invalid.
00773                 Valid values are  +1.0  or  -1.0  only.
00774 
00775          Input Arguments are sets of control points...
00776          For Distort Images    u,v, x,y  ...
00777          For Sparse Gradients  x,y, r,g,b  ...
00778 
00779          Perspective Distortion Notes...
00780            + Can be thought of as ratio of  3 affine transformations
00781            + Not separatable: r() or c6 and c7 are used by both equations
00782            + All 8 coefficients must be determined simultaniously
00783            + Will only work with a 2 number_valuesal Image Distortion
00784            + Can not be used for generating a sparse gradient (interpolation)
00785            + It is not linear, but is simple to generate an inverse
00786            + All lines within an image remain lines.
00787            + but distances between points may vary.
00788       */
00789       double
00790         **matrix,
00791         *vectors[1],
00792         terms[8];
00793 
00794       size_t
00795         cp_u = cp_values,
00796         cp_v = cp_values+1;
00797 
00798       MagickBooleanType
00799         status;
00800 
00801       if ( number_arguments%cp_size != 0 ||
00802            number_arguments < cp_size*4 ) {
00803         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00804               "InvalidArgument", "%s : 'require at least %.20g CPs'",
00805               CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0);
00806         coeff=(double *) RelinquishMagickMemory(coeff);
00807         return((double *) NULL);
00808       }
00809       /* fake 1x8 vectors matrix directly using the coefficients array */
00810       vectors[0] = &(coeff[0]);
00811       /* 8x8 least-squares matrix (zeroed) */
00812       matrix = AcquireMagickMatrix(8UL,8UL);
00813       if (matrix == (double **) NULL) {
00814         (void) ThrowMagickException(exception,GetMagickModule(),
00815                   ResourceLimitError,"MemoryAllocationFailed",
00816                   "%s", "DistortCoefficients");
00817         return((double *) NULL);
00818       }
00819       /* Add control points for least squares solving */
00820       for (i=0; i < number_arguments; i+=4) {
00821         terms[0]=arguments[i+cp_x];            /*   c0*x   */
00822         terms[1]=arguments[i+cp_y];            /*   c1*y   */
00823         terms[2]=1.0;                          /*   c2*1   */
00824         terms[3]=0.0;
00825         terms[4]=0.0;
00826         terms[5]=0.0;
00827         terms[6]=-terms[0]*arguments[i+cp_u];  /* 1/(c6*x) */
00828         terms[7]=-terms[1]*arguments[i+cp_u];  /* 1/(c7*y) */
00829         LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_u]),
00830             8UL,1UL);
00831 
00832         terms[0]=0.0;
00833         terms[1]=0.0;
00834         terms[2]=0.0;
00835         terms[3]=arguments[i+cp_x];           /*   c3*x   */
00836         terms[4]=arguments[i+cp_y];           /*   c4*y   */
00837         terms[5]=1.0;                         /*   c5*1   */
00838         terms[6]=-terms[3]*arguments[i+cp_v]; /* 1/(c6*x) */
00839         terms[7]=-terms[4]*arguments[i+cp_v]; /* 1/(c7*y) */
00840         LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_v]),
00841             8UL,1UL);
00842       }
00843       /* Solve for LeastSquares Coefficients */
00844       status=GaussJordanElimination(matrix,vectors,8UL,1UL);
00845       matrix = RelinquishMagickMatrix(matrix, 8UL);
00846       if ( status == MagickFalse ) {
00847         coeff = (double *) RelinquishMagickMemory(coeff);
00848         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00849             "InvalidArgument","%s : 'Unsolvable Matrix'",
00850             CommandOptionToMnemonic(MagickDistortOptions, *method) );
00851         return((double *) NULL);
00852       }
00853       /*
00854         Calculate 9'th coefficient! The ground-sky determination.
00855         What is sign of the 'ground' in r() denominator affine function?
00856         Just use any valid image coordinate (first control point) in
00857         destination for determination of what part of view is 'ground'.
00858       */
00859       coeff[8] = coeff[6]*arguments[cp_x]
00860                       + coeff[7]*arguments[cp_y] + 1.0;
00861       coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
00862 
00863       return(coeff);
00864     }
00865     case PerspectiveProjectionDistortion:
00866     {
00867       /*
00868         Arguments: Perspective Coefficents (forward mapping)
00869       */
00870       if (number_arguments != 8) {
00871         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00872               "InvalidArgument", "%s : 'Needs 8 coefficient values'",
00873               CommandOptionToMnemonic(MagickDistortOptions, *method));
00874         return((double *) NULL);
00875       }
00876       /* FUTURE: trap test  c0*c4-c3*c1 == 0  (determinate = 0, no inverse) */
00877       InvertPerspectiveCoefficients(arguments, coeff);
00878       /*
00879         Calculate 9'th coefficient! The ground-sky determination.
00880         What is sign of the 'ground' in r() denominator affine function?
00881         Just use any valid image cocodinate in destination for determination.
00882         For a forward mapped perspective the images 0,0 coord will map to
00883         c2,c5 in the distorted image, so set the sign of denominator of that.
00884       */
00885       coeff[8] = coeff[6]*arguments[2]
00886                            + coeff[7]*arguments[5] + 1.0;
00887       coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
00888       *method = PerspectiveDistortion;
00889 
00890       return(coeff);
00891     }
00892     case BilinearForwardDistortion:
00893     case BilinearReverseDistortion:
00894     {
00895       /* Bilinear Distortion (Forward mapping)
00896             v = c0*x + c1*y + c2*x*y + c3;
00897          for each 'value' given
00898 
00899          This is actually a simple polynomial Distortion!  The difference
00900          however is when we need to reverse the above equation to generate a
00901          BilinearForwardDistortion (see below).
00902 
00903          Input Arguments are sets of control points...
00904          For Distort Images    u,v, x,y  ...
00905          For Sparse Gradients  x,y, r,g,b  ...
00906 
00907       */
00908       double
00909         **matrix,
00910         **vectors,
00911         terms[4];
00912 
00913       MagickBooleanType
00914         status;
00915 
00916       /* check the number of arguments */
00917       if ( number_arguments%cp_size != 0 ||
00918            number_arguments < cp_size*4 ) {
00919         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00920               "InvalidArgument", "%s : 'require at least %.20g CPs'",
00921               CommandOptionToMnemonic(MagickDistortOptions, *method), 4.0);
00922         coeff=(double *) RelinquishMagickMemory(coeff);
00923         return((double *) NULL);
00924       }
00925       /* create matrix, and a fake vectors matrix */
00926       matrix = AcquireMagickMatrix(4UL,4UL);
00927       vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
00928       if (matrix == (double **) NULL || vectors == (double **) NULL)
00929       {
00930         matrix  = RelinquishMagickMatrix(matrix, 4UL);
00931         vectors = (double **) RelinquishMagickMemory(vectors);
00932         coeff   = (double *) RelinquishMagickMemory(coeff);
00933         (void) ThrowMagickException(exception,GetMagickModule(),
00934                 ResourceLimitError,"MemoryAllocationFailed",
00935                 "%s", "DistortCoefficients");
00936         return((double *) NULL);
00937       }
00938       /* fake a number_values x4 vectors matrix from coefficients array */
00939       for (i=0; i < number_values; i++)
00940         vectors[i] = &(coeff[i*4]);
00941       /* Add given control point pairs for least squares solving */
00942       for (i=0; i < number_arguments; i+=cp_size) {
00943         terms[0] = arguments[i+cp_x];   /*  x  */
00944         terms[1] = arguments[i+cp_y];   /*  y  */
00945         terms[2] = terms[0]*terms[1];   /* x*y */
00946         terms[3] = 1;                   /*  1  */
00947         LeastSquaresAddTerms(matrix,vectors,terms,
00948              &(arguments[i+cp_values]),4UL,number_values);
00949       }
00950       /* Solve for LeastSquares Coefficients */
00951       status=GaussJordanElimination(matrix,vectors,4UL,number_values);
00952       matrix  = RelinquishMagickMatrix(matrix, 4UL);
00953       vectors = (double **) RelinquishMagickMemory(vectors);
00954       if ( status == MagickFalse ) {
00955         coeff = (double *) RelinquishMagickMemory(coeff);
00956         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
00957             "InvalidArgument","%s : 'Unsolvable Matrix'",
00958             CommandOptionToMnemonic(MagickDistortOptions, *method) );
00959         return((double *) NULL);
00960       }
00961       if ( *method == BilinearForwardDistortion ) {
00962          /* Bilinear Forward Mapped Distortion
00963 
00964          The above least-squares solved for coefficents but in the forward
00965          direction, due to changes to indexing constants.
00966 
00967             i = c0*x + c1*y + c2*x*y + c3;
00968             j = c4*x + c5*y + c6*x*y + c7;
00969 
00970          where i,j are in the destination image, NOT the source.
00971 
00972          Reverse Pixel mapping however needs to use reverse of these
00973          functions.  It required a full page of algbra to work out the
00974          reversed mapping formula, but resolves down to the following...
00975 
00976             c8 = c0*c5-c1*c4;
00977             c9 = 2*(c2*c5-c1*c6);   // '2*a' in the quadratic formula
00978 
00979             i = i - c3;   j = j - c7;
00980             b = c6*i - c2*j + c8;   // So that   a*y^2 + b*y + c == 0
00981             c = c4*i -  c0*j;       // y = ( -b +- sqrt(bb - 4ac) ) / (2*a)
00982 
00983             r = b*b - c9*(c+c);
00984             if ( c9 != 0 )
00985               y = ( -b + sqrt(r) ) / c9;
00986             else
00987               y = -c/b;
00988 
00989             x = ( i - c1*y) / ( c1 - c2*y );
00990 
00991          NB: if 'r' is negative there is no solution!
00992          NB: the sign of the sqrt() should be negative if image becomes
00993              flipped or flopped, or crosses over itself.
00994          NB: techniqually coefficient c5 is not needed, anymore,
00995              but kept for completness.
00996 
00997          See Anthony Thyssen <A.Thyssen@griffith.edu.au>
00998          or  Fred Weinhaus <fmw@alink.net>  for more details.
00999 
01000          */
01001          coeff[8] = coeff[0]*coeff[5] - coeff[1]*coeff[4];
01002          coeff[9] = 2*(coeff[2]*coeff[5] - coeff[1]*coeff[6]);
01003       }
01004       return(coeff);
01005     }
01006 #if 0
01007     case QuadrilateralDistortion:
01008     {
01009       /* Map a Quadrilateral to a unit square using BilinearReverse
01010          Then map that unit square back to the final Quadrilateral
01011          using BilinearForward.
01012 
01013          Input Arguments are sets of control points...
01014          For Distort Images    u,v, x,y  ...
01015          For Sparse Gradients  x,y, r,g,b  ...
01016 
01017       */
01018       /* UNDER CONSTRUCTION */
01019       return(coeff);
01020     }
01021 #endif
01022 
01023     case PolynomialDistortion:
01024     {
01025       /* Polynomial Distortion
01026 
01027          First two coefficents are used to hole global polynomal information
01028            c0 = Order of the polynimial being created
01029            c1 = number_of_terms in one polynomial equation
01030 
01031          Rest of the coefficients map to the equations....
01032             v = c0 + c1*x + c2*y + c3*x*y + c4*x^2 + c5*y^2 + c6*x^3 + ...
01033          for each 'value' (number_values of them) given.
01034          As such total coefficients =  2 + number_terms * number_values
01035 
01036          Input Arguments are sets of control points...
01037          For Distort Images    order  [u,v, x,y] ...
01038          For Sparse Gradients  order  [x,y, r,g,b] ...
01039 
01040          Polynomial Distortion Notes...
01041            + UNDER DEVELOPMENT -- Do not expect this to remain as is.
01042            + Currently polynomial is a reversed mapped distortion.
01043            + Order 1.5 is fudged to map into a bilinear distortion.
01044              though it is not the same order as that distortion.
01045       */
01046       double
01047         **matrix,
01048         **vectors,
01049         *terms;
01050 
01051       size_t
01052         nterms;   /* number of polynomial terms per number_values */
01053 
01054       register ssize_t
01055         j;
01056 
01057       MagickBooleanType
01058         status;
01059 
01060       /* first two coefficients hold polynomial order information */
01061       coeff[0] = arguments[0];
01062       coeff[1] = (double) poly_number_terms(arguments[0]);
01063       nterms = (size_t) coeff[1];
01064 
01065       /* create matrix, a fake vectors matrix, and least sqs terms */
01066       matrix = AcquireMagickMatrix(nterms,nterms);
01067       vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
01068       terms = (double *) AcquireQuantumMemory(nterms, sizeof(*terms));
01069       if (matrix  == (double **) NULL ||
01070           vectors == (double **) NULL ||
01071           terms   == (double *) NULL )
01072       {
01073         matrix  = RelinquishMagickMatrix(matrix, nterms);
01074         vectors = (double **) RelinquishMagickMemory(vectors);
01075         terms   = (double *) RelinquishMagickMemory(terms);
01076         coeff   = (double *) RelinquishMagickMemory(coeff);
01077         (void) ThrowMagickException(exception,GetMagickModule(),
01078                 ResourceLimitError,"MemoryAllocationFailed",
01079                 "%s", "DistortCoefficients");
01080         return((double *) NULL);
01081       }
01082       /* fake a number_values x3 vectors matrix from coefficients array */
01083       for (i=0; i < number_values; i++)
01084         vectors[i] = &(coeff[2+i*nterms]);
01085       /* Add given control point pairs for least squares solving */
01086       for (i=1; i < number_arguments; i+=cp_size) { /* NB: start = 1 not 0 */
01087         for (j=0; j < (ssize_t) nterms; j++)
01088           terms[j] = poly_basis_fn(j,arguments[i+cp_x],arguments[i+cp_y]);
01089         LeastSquaresAddTerms(matrix,vectors,terms,
01090              &(arguments[i+cp_values]),nterms,number_values);
01091       }
01092       terms = (double *) RelinquishMagickMemory(terms);
01093       /* Solve for LeastSquares Coefficients */
01094       status=GaussJordanElimination(matrix,vectors,nterms,number_values);
01095       matrix  = RelinquishMagickMatrix(matrix, nterms);
01096       vectors = (double **) RelinquishMagickMemory(vectors);
01097       if ( status == MagickFalse ) {
01098         coeff = (double *) RelinquishMagickMemory(coeff);
01099         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
01100             "InvalidArgument","%s : 'Unsolvable Matrix'",
01101             CommandOptionToMnemonic(MagickDistortOptions, *method) );
01102         return((double *) NULL);
01103       }
01104       return(coeff);
01105     }
01106     case ArcDistortion:
01107     {
01108       /* Arc Distortion
01109          Args: arc_width  rotate  top_edge_radius  bottom_edge_radius
01110          All but first argument are optional
01111             arc_width      The angle over which to arc the image side-to-side
01112             rotate         Angle to rotate image from vertical center
01113             top_radius     Set top edge of source image at this radius
01114             bottom_radius  Set bootom edge to this radius (radial scaling)
01115 
01116          By default, if the radii arguments are nor provided the image radius
01117          is calculated so the horizontal center-line is fits the given arc
01118          without scaling.
01119 
01120          The output image size is ALWAYS adjusted to contain the whole image,
01121          and an offset is given to position image relative to the 0,0 point of
01122          the origin, allowing users to use relative positioning onto larger
01123          background (via -flatten).
01124 
01125          The arguments are converted to these coefficients
01126             c0: angle for center of source image
01127             c1: angle scale for mapping to source image
01128             c2: radius for top of source image
01129             c3: radius scale for mapping source image
01130             c4: centerline of arc within source image
01131 
01132          Note the coefficients use a center angle, so asymptotic join is
01133          furthest from both sides of the source image. This also means that
01134          for arc angles greater than 360 the sides of the image will be
01135          trimmed equally.
01136 
01137          Arc Distortion Notes...
01138            + Does not use a set of CPs
01139            + Will only work with Image Distortion
01140            + Can not be used for generating a sparse gradient (interpolation)
01141       */
01142       if ( number_arguments >= 1 && arguments[0] < MagickEpsilon ) {
01143         coeff = (double *) RelinquishMagickMemory(coeff);
01144         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
01145             "InvalidArgument","%s : 'Arc Angle Too Small'",
01146             CommandOptionToMnemonic(MagickDistortOptions, *method) );
01147         return((double *) NULL);
01148       }
01149       if ( number_arguments >= 3 && arguments[2] < MagickEpsilon ) {
01150         coeff = (double *) RelinquishMagickMemory(coeff);
01151         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
01152             "InvalidArgument","%s : 'Outer Radius Too Small'",
01153             CommandOptionToMnemonic(MagickDistortOptions, *method) );
01154         return((double *) NULL);
01155       }
01156       coeff[0] = -MagickPI2;   /* -90, place at top! */
01157       if ( number_arguments >= 1 )
01158         coeff[1] = DegreesToRadians(arguments[0]);
01159       else
01160         coeff[1] = MagickPI2;   /* zero arguments - center is at top */
01161       if ( number_arguments >= 2 )
01162         coeff[0] += DegreesToRadians(arguments[1]);
01163       coeff[0] /= Magick2PI;  /* normalize radians */
01164       coeff[0] -= MagickRound(coeff[0]);
01165       coeff[0] *= Magick2PI;  /* de-normalize back to radians */
01166       coeff[3] = (double)image->rows-1;
01167       coeff[2] = (double)image->columns/coeff[1] + coeff[3]/2.0;
01168       if ( number_arguments >= 3 ) {
01169         if ( number_arguments >= 4 )
01170           coeff[3] = arguments[2] - arguments[3];
01171         else
01172           coeff[3] *= arguments[2]/coeff[2];
01173         coeff[2] = arguments[2];
01174       }
01175       coeff[4] = ((double)image->columns-1.0)/2.0;
01176 
01177       return(coeff);
01178     }
01179     case PolarDistortion:
01180     case DePolarDistortion:
01181     {
01182       /* (De)Polar Distortion   (same set of arguments)
01183          Args:  Rmax, Rmin,  Xcenter,Ycenter,  Afrom,Ato
01184          DePolar can also have the extra arguments of Width, Height
01185 
01186          Coefficients 0 to 5 is the sanatized version first 6 input args
01187          Coefficient 6  is the angle to coord ratio  and visa-versa
01188          Coefficient 7  is the radius to coord ratio and visa-versa
01189 
01190          WARNING: It is possible for  Radius max<min  and/or  Angle from>to
01191       */
01192       if ( number_arguments == 3
01193           || ( number_arguments > 6 && *method == PolarDistortion )
01194           || number_arguments > 8 ) {
01195           (void) ThrowMagickException(exception,GetMagickModule(),
01196             OptionError,"InvalidArgument", "%s : number of arguments",
01197             CommandOptionToMnemonic(MagickDistortOptions, *method) );
01198         coeff=(double *) RelinquishMagickMemory(coeff);
01199         return((double *) NULL);
01200       }
01201       /* Rmax -  if 0 calculate appropriate value */
01202       if ( number_arguments >= 1 )
01203         coeff[0] = arguments[0];
01204       else
01205         coeff[0] = 0.0;
01206       /* Rmin  - usally 0 */
01207       coeff[1] = number_arguments >= 2 ? arguments[1] : 0.0;
01208       /* Center X,Y */
01209       if ( number_arguments >= 4 ) {
01210         coeff[2] = arguments[2];
01211         coeff[3] = arguments[3];
01212       }
01213       else { /* center of actual image */
01214         coeff[2] = (double)(image->columns)/2.0+image->page.x;
01215         coeff[3] = (double)(image->rows)/2.0+image->page.y;
01216       }
01217       /* Angle from,to - about polar center 0 is downward */
01218       coeff[4] = -MagickPI;
01219       if ( number_arguments >= 5 )
01220         coeff[4] = DegreesToRadians(arguments[4]);
01221       coeff[5] = coeff[4];
01222       if ( number_arguments >= 6 )
01223         coeff[5] = DegreesToRadians(arguments[5]);
01224       if ( fabs(coeff[4]-coeff[5]) < MagickEpsilon )
01225         coeff[5] += Magick2PI; /* same angle is a full circle */
01226       /* if radius 0 or negative,  its a special value... */
01227       if ( coeff[0] < MagickEpsilon ) {
01228         /* Use closest edge  if radius == 0 */
01229         if ( fabs(coeff[0]) < MagickEpsilon ) {
01230           coeff[0]=MagickMin(fabs(coeff[2]-image->page.x),
01231                              fabs(coeff[3]-image->page.y));
01232           coeff[0]=MagickMin(coeff[0],
01233                        fabs(coeff[2]-image->page.x-image->columns));
01234           coeff[0]=MagickMin(coeff[0],
01235                        fabs(coeff[3]-image->page.y-image->rows));
01236         }
01237         /* furthest diagonal if radius == -1 */
01238         if ( fabs(-1.0-coeff[0]) < MagickEpsilon ) {
01239           double rx,ry;
01240           rx = coeff[2]-image->page.x;
01241           ry = coeff[3]-image->page.y;
01242           coeff[0] = rx*rx+ry*ry;
01243           ry = coeff[3]-image->page.y-image->rows;
01244           coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
01245           rx = coeff[2]-image->page.x-image->columns;
01246           coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
01247           ry = coeff[3]-image->page.y;
01248           coeff[0] = MagickMax(coeff[0],rx*rx+ry*ry);
01249           coeff[0] = sqrt(coeff[0]);
01250         }
01251       }
01252       /* IF Rmax <= 0 or Rmin < 0 OR Rmax < Rmin, THEN error */
01253       if ( coeff[0] < MagickEpsilon || coeff[1] < -MagickEpsilon
01254            || (coeff[0]-coeff[1]) < MagickEpsilon ) {
01255         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
01256             "InvalidArgument", "%s : Invalid Radius",
01257             CommandOptionToMnemonic(MagickDistortOptions, *method) );
01258         coeff=(double *) RelinquishMagickMemory(coeff);
01259         return((double *) NULL);
01260       }
01261       /* converstion ratios */
01262       if ( *method == PolarDistortion ) {
01263         coeff[6]=(double) image->columns/(coeff[5]-coeff[4]);
01264         coeff[7]=(double) image->rows/(coeff[0]-coeff[1]);
01265       }
01266       else { /* *method == DePolarDistortion */
01267         coeff[6]=(coeff[5]-coeff[4])/image->columns;
01268         coeff[7]=(coeff[0]-coeff[1])/image->rows;
01269       }
01270       return(coeff);
01271     }
01272     case Cylinder2PlaneDistortion:
01273     case Plane2CylinderDistortion:
01274     {
01275       /* 3D Cylinder to/from a Tangential Plane
01276 
01277          Projection between a clinder and flat plain from a point on the
01278          center line of the cylinder.
01279 
01280          The two surfaces coincide in 3D space at the given centers of
01281          distortion (perpendicular to projection point) on both images.
01282 
01283          Args:  FOV_arc_width
01284          Coefficents: FOV(radians), Radius, center_x,y, dest_center_x,y
01285 
01286          FOV (Field Of View) the angular field of view of the distortion,
01287          across the width of the image, in degrees.  The centers are the
01288          points of least distortion in the input and resulting images.
01289 
01290          These centers are however determined later.
01291 
01292          Coeff 0 is the FOV angle of view of image width in radians
01293          Coeff 1 is calculated radius of cylinder.
01294          Coeff 2,3  center of distortion of input image
01295          Coefficents 4,5 Center of Distortion of dest (determined later)
01296       */
01297       if ( arguments[0] < MagickEpsilon || arguments[0] > 160.0 ) {
01298         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
01299             "InvalidArgument", "%s : Invalid FOV Angle",
01300             CommandOptionToMnemonic(MagickDistortOptions, *method) );
01301         coeff=(double *) RelinquishMagickMemory(coeff);
01302         return((double *) NULL);
01303       }
01304       coeff[0] = DegreesToRadians(arguments[0]);
01305       if ( *method == Cylinder2PlaneDistortion )
01306         /* image is curved around cylinder, so FOV angle (in radians)
01307          * scales directly to image X coordinate, according to its radius.
01308          */
01309         coeff[1] = (double) image->columns/coeff[0];
01310       else
01311         /* radius is distance away from an image with this angular FOV */
01312         coeff[1] = (double) image->columns / ( 2 * tan(coeff[0]/2) );
01313 
01314       coeff[2] = (double)(image->columns)/2.0+image->page.x;
01315       coeff[3] = (double)(image->rows)/2.0+image->page.y;
01316       coeff[4] = coeff[2];
01317       coeff[5] = coeff[3]; /* assuming image size is the same */
01318       return(coeff);
01319     }
01320     case BarrelDistortion:
01321     case BarrelInverseDistortion:
01322     {
01323       /* Barrel Distortion
01324            Rs=(A*Rd^3 + B*Rd^2 + C*Rd + D)*Rd
01325          BarrelInv Distortion
01326            Rs=Rd/(A*Rd^3 + B*Rd^2 + C*Rd + D)
01327 
01328         Where Rd is the normalized radius from corner to middle of image
01329         Input Arguments are one of the following forms (number of arguments)...
01330             3:  A,B,C
01331             4:  A,B,C,D
01332             5:  A,B,C    X,Y
01333             6:  A,B,C,D  X,Y
01334             8:  Ax,Bx,Cx,Dx  Ay,By,Cy,Dy
01335            10:  Ax,Bx,Cx,Dx  Ay,By,Cy,Dy   X,Y
01336 
01337         Returns 10 coefficent values, which are de-normalized (pixel scale)
01338           Ax, Bx, Cx, Dx,   Ay, By, Cy, Dy,    Xc, Yc
01339       */
01340       /* Radius de-normalization scaling factor */
01341       double
01342         rscale = 2.0/MagickMin((double) image->columns,(double) image->rows);
01343 
01344       /* sanity check  number of args must = 3,4,5,6,8,10 or error */
01345       if ( (number_arguments  < 3) || (number_arguments == 7) ||
01346            (number_arguments == 9) || (number_arguments > 10) )
01347         {
01348           coeff=(double *) RelinquishMagickMemory(coeff);
01349           (void) ThrowMagickException(exception,GetMagickModule(),
01350             OptionError,"InvalidArgument", "%s : number of arguments",
01351             CommandOptionToMnemonic(MagickDistortOptions, *method) );
01352           return((double *) NULL);
01353         }
01354       /* A,B,C,D coefficients */
01355       coeff[0] = arguments[0];
01356       coeff[1] = arguments[1];
01357       coeff[2] = arguments[2];
01358       if ((number_arguments == 3) || (number_arguments == 5) )
01359         coeff[3] = 1.0 - coeff[0] - coeff[1] - coeff[2];
01360       else
01361         coeff[3] = arguments[3];
01362       /* de-normalize the coefficients */
01363       coeff[0] *= pow(rscale,3.0);
01364       coeff[1] *= rscale*rscale;
01365       coeff[2] *= rscale;
01366       /* Y coefficients: as given OR same as X coefficients */
01367       if ( number_arguments >= 8 ) {
01368         coeff[4] = arguments[4] * pow(rscale,3.0);
01369         coeff[5] = arguments[5] * rscale*rscale;
01370         coeff[6] = arguments[6] * rscale;
01371         coeff[7] = arguments[7];
01372       }
01373       else {
01374         coeff[4] = coeff[0];
01375         coeff[5] = coeff[1];
01376         coeff[6] = coeff[2];
01377         coeff[7] = coeff[3];
01378       }
01379       /* X,Y Center of Distortion (image coodinates) */
01380       if ( number_arguments == 5 )  {
01381         coeff[8] = arguments[3];
01382         coeff[9] = arguments[4];
01383       }
01384       else if ( number_arguments == 6 ) {
01385         coeff[8] = arguments[4];
01386         coeff[9] = arguments[5];
01387       }
01388       else if ( number_arguments == 10 ) {
01389         coeff[8] = arguments[8];
01390         coeff[9] = arguments[9];
01391       }
01392       else {
01393         /* center of the image provided (image coodinates) */
01394         coeff[8] = (double)image->columns/2.0 + image->page.x;
01395         coeff[9] = (double)image->rows/2.0    + image->page.y;
01396       }
01397       return(coeff);
01398     }
01399     case ShepardsDistortion:
01400     {
01401       /* Shepards Distortion  input arguments are the coefficents!
01402          Just check the number of arguments is valid!
01403          Args:  u1,v1, x1,y1, ...
01404           OR :  u1,v1, r1,g1,c1, ...
01405       */
01406       if ( number_arguments%cp_size != 0 ||
01407            number_arguments < cp_size ) {
01408         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
01409               "InvalidArgument", "%s : 'require at least %.20g CPs'",
01410               CommandOptionToMnemonic(MagickDistortOptions, *method), 1.0);
01411         coeff=(double *) RelinquishMagickMemory(coeff);
01412         return((double *) NULL);
01413       }
01414       return(coeff);
01415     }
01416     default:
01417       break;
01418   }
01419   /* you should never reach this point */
01420   assert(! "No Method Handler"); /* just fail assertion */
01421   return((double *) NULL);
01422 }
01423 
01424 /*
01425 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01426 %                                                                             %
01427 %                                                                             %
01428 %                                                                             %
01429 +   D i s t o r t R e s i z e I m a g e                                       %
01430 %                                                                             %
01431 %                                                                             %
01432 %                                                                             %
01433 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01434 %
01435 %  DistortResizeImage() resize image using the equivalent but slower image
01436 %  distortion operator.  The filter is applied using a EWA cylindrical
01437 %  resampling. But like resize the final image size is limited to whole pixels
01438 %  with no effects by virtual-pixels on the result.
01439 %
01440 %  Note that images containing a transparency channel will be twice as slow to
01441 %  resize as images one without transparency.
01442 %
01443 %  The format of the DistortResizeImage method is:
01444 %
01445 %      Image *AdaptiveResizeImage(const Image *image,const size_t columns,
01446 %        const size_t rows,ExceptionInfo *exception)
01447 %
01448 %  A description of each parameter follows:
01449 %
01450 %    o image: the image.
01451 %
01452 %    o columns: the number of columns in the resized image.
01453 %
01454 %    o rows: the number of rows in the resized image.
01455 %
01456 %    o exception: return any errors or warnings in this structure.
01457 %
01458 */
01459 MagickExport Image *DistortResizeImage(const Image *image,
01460   const size_t columns,const size_t rows,ExceptionInfo *exception)
01461 {
01462 #define DistortResizeImageTag  "Distort/Image"
01463 
01464   Image
01465     *resize_image,
01466     *tmp_image;
01467 
01468   RectangleInfo
01469     crop_area;
01470 
01471   double
01472     distort_args[12];
01473 
01474   VirtualPixelMethod
01475     vp_save;
01476 
01477   /*
01478     Distort resize image.
01479   */
01480   assert(image != (const Image *) NULL);
01481   assert(image->signature == MagickSignature);
01482   if (image->debug != MagickFalse)
01483     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01484   assert(exception != (ExceptionInfo *) NULL);
01485   assert(exception->signature == MagickSignature);
01486   if ((columns == 0) || (rows == 0))
01487     return((Image *) NULL);
01488   /* Do not short-circuit this resize if final image size is unchanged */
01489 
01490 
01491   (void) ResetMagickMemory(distort_args,0,12*sizeof(double));
01492   distort_args[4]=(double) image->columns;
01493   distort_args[6]=(double) columns;
01494   distort_args[9]=(double) image->rows;
01495   distort_args[11]=(double) rows;
01496 
01497   vp_save=GetImageVirtualPixelMethod(image);
01498 
01499   tmp_image=CloneImage(image,0,0,MagickTrue,exception);
01500   if ( tmp_image == (Image *) NULL )
01501     return((Image *) NULL);
01502   (void) SetImageVirtualPixelMethod(tmp_image,TransparentVirtualPixelMethod,
01503     exception);
01504 
01505   if (image->matte == MagickFalse)
01506     {
01507       /*
01508         Image has not transparency channel, so we free to use it
01509       */
01510       (void) SetImageAlphaChannel(tmp_image,SetAlphaChannel,exception);
01511       resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
01512         MagickTrue,exception),
01513 
01514       tmp_image=DestroyImage(tmp_image);
01515       if ( resize_image == (Image *) NULL )
01516         return((Image *) NULL);
01517 
01518       (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel,exception);
01519     }
01520   else
01521     {
01522       Image
01523         *resize_alpha;
01524 
01525       /*
01526         Image has transparency so handle colors and alpha separatly.
01527         Basically we need to separate Virtual-Pixel alpha in the resized
01528         image, so only the actual original images alpha channel is used.
01529 
01530         distort alpha channel separately
01531       */
01532       (void) SetImageAlphaChannel(tmp_image,ExtractAlphaChannel,exception);
01533       (void) SetImageAlphaChannel(tmp_image,OpaqueAlphaChannel,exception);
01534       resize_alpha=DistortImage(tmp_image,AffineDistortion,12,distort_args,
01535         MagickTrue,exception),
01536       tmp_image=DestroyImage(tmp_image);
01537       if (resize_alpha == (Image *) NULL)
01538         return((Image *) NULL);
01539 
01540       /* distort the actual image containing alpha + VP alpha */
01541       tmp_image=CloneImage(image,0,0,MagickTrue,exception);
01542       if ( tmp_image == (Image *) NULL )
01543         return((Image *) NULL);
01544       (void) SetImageVirtualPixelMethod(tmp_image,
01545         TransparentVirtualPixelMethod,exception);
01546       resize_image=DistortImage(tmp_image,AffineDistortion,12,distort_args,
01547         MagickTrue,exception),
01548       tmp_image=DestroyImage(tmp_image);
01549       if ( resize_image == (Image *) NULL)
01550         {
01551           resize_alpha=DestroyImage(resize_alpha);
01552           return((Image *) NULL);
01553         }
01554       /* replace resize images alpha with the separally distorted alpha */
01555       (void) SetImageAlphaChannel(resize_image,DeactivateAlphaChannel,
01556         exception);
01557       (void) SetImageAlphaChannel(resize_alpha,DeactivateAlphaChannel,
01558         exception);
01559       (void) CompositeImage(resize_image,CopyAlphaCompositeOp,resize_alpha,
01560         0,0,exception);
01561       resize_alpha=DestroyImage(resize_alpha);
01562     }
01563   (void) SetImageVirtualPixelMethod(resize_image,vp_save,exception);
01564 
01565   /*
01566     Clean up the results of the Distortion
01567   */
01568   crop_area.width=columns;
01569   crop_area.height=rows;
01570   crop_area.x=0;
01571   crop_area.y=0;
01572 
01573   tmp_image=resize_image;
01574   resize_image=CropImage(tmp_image,&crop_area,exception);
01575   tmp_image=DestroyImage(tmp_image);
01576 
01577   if ( resize_image == (Image *) NULL )
01578     return((Image *) NULL);
01579 
01580   return(resize_image);
01581 }
01582 
01583 /*
01584 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01585 %                                                                             %
01586 %                                                                             %
01587 %                                                                             %
01588 %   D i s t o r t I m a g e                                                   %
01589 %                                                                             %
01590 %                                                                             %
01591 %                                                                             %
01592 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01593 %
01594 %  DistortImage() distorts an image using various distortion methods, by
01595 %  mapping color lookups of the source image to a new destination image
01596 %  usally of the same size as the source image, unless 'bestfit' is set to
01597 %  true.
01598 %
01599 %  If 'bestfit' is enabled, and distortion allows it, the destination image is
01600 %  adjusted to ensure the whole source 'image' will just fit within the final
01601 %  destination image, which will be sized and offset accordingly.  Also in
01602 %  many cases the virtual offset of the source image will be taken into
01603 %  account in the mapping.
01604 %
01605 %  If the '-verbose' control option has been set print to standard error the
01606 %  equicelent '-fx' formula with coefficients for the function, if practical.
01607 %
01608 %  The format of the DistortImage() method is:
01609 %
01610 %      Image *DistortImage(const Image *image,const DistortImageMethod method,
01611 %        const size_t number_arguments,const double *arguments,
01612 %        MagickBooleanType bestfit, ExceptionInfo *exception)
01613 %
01614 %  A description of each parameter follows:
01615 %
01616 %    o image: the image to be distorted.
01617 %
01618 %    o method: the method of image distortion.
01619 %
01620 %        ArcDistortion always ignores source image offset, and always
01621 %        'bestfit' the destination image with the top left corner offset
01622 %        relative to the polar mapping center.
01623 %
01624 %        Affine, Perspective, and Bilinear, do least squares fitting of the
01625 %        distrotion when more than the minimum number of control point pairs
01626 %        are provided.
01627 %
01628 %        Perspective, and Bilinear, fall back to a Affine distortion when less
01629 %        than 4 control point pairs are provided.  While Affine distortions
01630 %        let you use any number of control point pairs, that is Zero pairs is
01631 %        a No-Op (viewport only) distortion, one pair is a translation and
01632 %        two pairs of control points do a scale-rotate-translate, without any
01633 %        shearing.
01634 %
01635 %    o number_arguments: the number of arguments given.
01636 %
01637 %    o arguments: an array of floating point arguments for this method.
01638 %
01639 %    o bestfit: Attempt to 'bestfit' the size of the resulting image.
01640 %        This also forces the resulting image to be a 'layered' virtual
01641 %        canvas image.  Can be overridden using 'distort:viewport' setting.
01642 %
01643 %    o exception: return any errors or warnings in this structure
01644 %
01645 %  Extra Controls from Image meta-data (artifacts)...
01646 %
01647 %    o "verbose"
01648 %        Output to stderr alternatives, internal coefficents, and FX
01649 %        equivalents for the distortion operation (if feasible).
01650 %        This forms an extra check of the distortion method, and allows users
01651 %        access to the internal constants IM calculates for the distortion.
01652 %
01653 %    o "distort:viewport"
01654 %        Directly set the output image canvas area and offest to use for the
01655 %        resulting image, rather than use the original images canvas, or a
01656 %        calculated 'bestfit' canvas.
01657 %
01658 %    o "distort:scale"
01659 %        Scale the size of the output canvas by this amount to provide a
01660 %        method of Zooming, and for super-sampling the results.
01661 %
01662 %  Other settings that can effect results include
01663 %
01664 %    o 'interpolate' For source image lookups (scale enlargements)
01665 %
01666 %    o 'filter'      Set filter to use for area-resampling (scale shrinking).
01667 %                    Set to 'point' to turn off and use 'interpolate' lookup
01668 %                    instead
01669 %
01670 */
01671 MagickExport Image *DistortImage(const Image *image,DistortImageMethod method,
01672   const size_t number_arguments,const double *arguments,
01673   MagickBooleanType bestfit,ExceptionInfo *exception)
01674 {
01675 #define DistortImageTag  "Distort/Image"
01676 
01677   double
01678     *coeff,
01679     output_scaling;
01680 
01681   Image
01682     *distort_image;
01683 
01684   RectangleInfo
01685     geometry;  /* geometry of the distorted space viewport */
01686 
01687   MagickBooleanType
01688     viewport_given;
01689 
01690   assert(image != (Image *) NULL);
01691   assert(image->signature == MagickSignature);
01692   if (image->debug != MagickFalse)
01693     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
01694   assert(exception != (ExceptionInfo *) NULL);
01695   assert(exception->signature == MagickSignature);
01696 
01697 
01698   /*
01699     Handle Special Compound Distortions
01700   */
01701   if ( method == ResizeDistortion )
01702     {
01703       if ( number_arguments != 2 )
01704         {
01705           (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
01706                     "InvalidArgument","%s : '%s'","Resize",
01707                     "Invalid number of args: 2 only");
01708           return((Image *) NULL);
01709         }
01710       distort_image=DistortResizeImage(image,(size_t)arguments[0],
01711          (size_t)arguments[1], exception);
01712       return(distort_image);
01713     }
01714 
01715   /*
01716     Convert input arguments (usually as control points for reverse mapping)
01717     into mapping coefficients to apply the distortion.
01718 
01719     Note that some distortions are mapped to other distortions,
01720     and as such do not require specific code after this point.
01721   */
01722   coeff = GenerateCoefficients(image, &method, number_arguments,
01723       arguments, 0, exception);
01724   if ( coeff == (double *) NULL )
01725     return((Image *) NULL);
01726 
01727   /*
01728     Determine the size and offset for a 'bestfit' destination.
01729     Usally the four corners of the source image is enough.
01730   */
01731 
01732   /* default output image bounds, when no 'bestfit' is requested */
01733   geometry.width=image->columns;
01734   geometry.height=image->rows;
01735   geometry.x=0;
01736   geometry.y=0;
01737 
01738   if ( method == ArcDistortion ) {
01739     bestfit = MagickTrue;  /* always calculate a 'best fit' viewport */
01740   }
01741 
01742   /* Work out the 'best fit', (required for ArcDistortion) */
01743   if ( bestfit ) {
01744     PointInfo
01745       s,d,min,max;  /* source, dest coords --mapping--> min, max coords */
01746 
01747     MagickBooleanType
01748       fix_bounds = MagickTrue;   /* enlarge bounds for VP handling */
01749 
01750     s.x=s.y=min.x=max.x=min.y=max.y=0.0;   /* keep compiler happy */
01751 
01752 /* defines to figure out the bounds of the distorted image */
01753 #define InitalBounds(p) \
01754 { \
01755   /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
01756   min.x = max.x = p.x; \
01757   min.y = max.y = p.y; \
01758 }
01759 #define ExpandBounds(p) \
01760 { \
01761   /* printf("%lg,%lg -> %lg,%lg\n", s.x,s.y, d.x,d.y); */ \
01762   min.x = MagickMin(min.x,p.x); \
01763   max.x = MagickMax(max.x,p.x); \
01764   min.y = MagickMin(min.y,p.y); \
01765   max.y = MagickMax(max.y,p.y); \
01766 }
01767 
01768     switch (method)
01769     {
01770       case AffineDistortion:
01771       { double inverse[6];
01772         InvertAffineCoefficients(coeff, inverse);
01773         s.x = (double) image->page.x;
01774         s.y = (double) image->page.y;
01775         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
01776         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
01777         InitalBounds(d);
01778         s.x = (double) image->page.x+image->columns;
01779         s.y = (double) image->page.y;
01780         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
01781         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
01782         ExpandBounds(d);
01783         s.x = (double) image->page.x;
01784         s.y = (double) image->page.y+image->rows;
01785         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
01786         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
01787         ExpandBounds(d);
01788         s.x = (double) image->page.x+image->columns;
01789         s.y = (double) image->page.y+image->rows;
01790         d.x = inverse[0]*s.x+inverse[1]*s.y+inverse[2];
01791         d.y = inverse[3]*s.x+inverse[4]*s.y+inverse[5];
01792         ExpandBounds(d);
01793         break;
01794       }
01795       case PerspectiveDistortion:
01796       { double inverse[8], scale;
01797         InvertPerspectiveCoefficients(coeff, inverse);
01798         s.x = (double) image->page.x;
01799         s.y = (double) image->page.y;
01800         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
01801         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
01802         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
01803         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
01804         InitalBounds(d);
01805         s.x = (double) image->page.x+image->columns;
01806         s.y = (double) image->page.y;
01807         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
01808         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
01809         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
01810         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
01811         ExpandBounds(d);
01812         s.x = (double) image->page.x;
01813         s.y = (double) image->page.y+image->rows;
01814         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
01815         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
01816         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
01817         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
01818         ExpandBounds(d);
01819         s.x = (double) image->page.x+image->columns;
01820         s.y = (double) image->page.y+image->rows;
01821         scale=inverse[6]*s.x+inverse[7]*s.y+1.0;
01822         scale=1.0/(  (fabs(scale) <= MagickEpsilon) ? 1.0 : scale );
01823         d.x = scale*(inverse[0]*s.x+inverse[1]*s.y+inverse[2]);
01824         d.y = scale*(inverse[3]*s.x+inverse[4]*s.y+inverse[5]);
01825         ExpandBounds(d);
01826         break;
01827       }
01828       case ArcDistortion:
01829       { double a, ca, sa;
01830         /* Forward Map Corners */
01831         a = coeff[0]-coeff[1]/2; ca = cos(a); sa = sin(a);
01832         d.x = coeff[2]*ca;
01833         d.y = coeff[2]*sa;
01834         InitalBounds(d);
01835         d.x = (coeff[2]-coeff[3])*ca;
01836         d.y = (coeff[2]-coeff[3])*sa;
01837         ExpandBounds(d);
01838         a = coeff[0]+coeff[1]/2; ca = cos(a); sa = sin(a);
01839         d.x = coeff[2]*ca;
01840         d.y = coeff[2]*sa;
01841         ExpandBounds(d);
01842         d.x = (coeff[2]-coeff[3])*ca;
01843         d.y = (coeff[2]-coeff[3])*sa;
01844         ExpandBounds(d);
01845         /* Orthogonal points along top of arc */
01846         for( a=(double) (ceil((double) ((coeff[0]-coeff[1]/2.0)/MagickPI2))*MagickPI2);
01847                a<(coeff[0]+coeff[1]/2.0); a+=MagickPI2 ) {
01848           ca = cos(a); sa = sin(a);
01849           d.x = coeff[2]*ca;
01850           d.y = coeff[2]*sa;
01851           ExpandBounds(d);
01852         }
01853         /*
01854           Convert the angle_to_width and radius_to_height
01855           to appropriate scaling factors, to allow faster processing
01856           in the mapping function.
01857         */
01858         coeff[1] = (double) (Magick2PI*image->columns/coeff[1]);
01859         coeff[3] = (double)image->rows/coeff[3];
01860         break;
01861       }
01862       case PolarDistortion:
01863       {
01864         if (number_arguments < 2)
01865           coeff[2] = coeff[3] = 0.0;
01866         min.x = coeff[2]-coeff[0];
01867         max.x = coeff[2]+coeff[0];
01868         min.y = coeff[3]-coeff[0];
01869         max.y = coeff[3]+coeff[0];
01870         /* should be about 1.0 if Rmin = 0 */
01871         coeff[7]=(double) geometry.height/(coeff[0]-coeff[1]);
01872         break;
01873       }
01874       case DePolarDistortion:
01875       {
01876         /* direct calculation as it needs to tile correctly
01877          * for reversibility in a DePolar-Polar cycle */
01878         fix_bounds = MagickFalse;
01879         geometry.x = geometry.y = 0;
01880         geometry.height = (size_t) ceil(coeff[0]-coeff[1]);
01881         geometry.width = (size_t)
01882                   ceil((coeff[0]-coeff[1])*(coeff[5]-coeff[4])*0.5);
01883         /* correct scaling factors relative to new size */
01884         coeff[6]=(coeff[5]-coeff[4])/geometry.width; /* changed width */
01885         coeff[7]=(coeff[0]-coeff[1])/geometry.height; /* should be about 1.0 */
01886         break;
01887       }
01888       case Cylinder2PlaneDistortion:
01889       {
01890         /* direct calculation so center of distortion is either a pixel
01891          * center, or pixel edge. This allows for reversibility of the
01892          * distortion */
01893         geometry.x = geometry.y = 0;
01894         geometry.width = (size_t) ceil( 2.0*coeff[1]*tan(coeff[0]/2.0) );
01895         geometry.height = (size_t) ceil( 2.0*coeff[3]/cos(coeff[0]/2.0) );
01896         /* correct center of distortion relative to new size */
01897         coeff[4] = (double) geometry.width/2.0;
01898         coeff[5] = (double) geometry.height/2.0;
01899         fix_bounds = MagickFalse;
01900         break;
01901       }
01902       case Plane2CylinderDistortion:
01903       {
01904         /* direct calculation center is either pixel center, or pixel edge
01905          * so as to allow reversibility of the image distortion */
01906         geometry.x = geometry.y = 0;
01907         geometry.width = (size_t) ceil(coeff[0]*coeff[1]);  /* FOV * radius */
01908         geometry.height = (size_t) (2*coeff[3]);              /* input image height */
01909         /* correct center of distortion relative to new size */
01910         coeff[4] = (double) geometry.width/2.0;
01911         coeff[5] = (double) geometry.height/2.0;
01912         fix_bounds = MagickFalse;
01913         break;
01914       }
01915       case ShepardsDistortion:
01916       case BilinearForwardDistortion:
01917       case BilinearReverseDistortion:
01918 #if 0
01919       case QuadrilateralDistortion:
01920 #endif
01921       case PolynomialDistortion:
01922       case BarrelDistortion:
01923       case BarrelInverseDistortion:
01924       default:
01925         /* no calculated bestfit available for these distortions */
01926         bestfit = MagickFalse;
01927         fix_bounds = MagickFalse;
01928         break;
01929     }
01930 
01931     /* Set the output image geometry to calculated 'bestfit'.
01932        Yes this tends to 'over do' the file image size, ON PURPOSE!
01933        Do not do this for DePolar which needs to be exact for virtual tiling.
01934     */
01935     if ( fix_bounds ) {
01936       geometry.x = (ssize_t) floor(min.x-0.5);
01937       geometry.y = (ssize_t) floor(min.y-0.5);
01938       geometry.width=(size_t) ceil(max.x-geometry.x+0.5);
01939       geometry.height=(size_t) ceil(max.y-geometry.y+0.5);
01940     }
01941 
01942   }  /* end bestfit destination image calculations */
01943 
01944   /* The user provided a 'viewport' expert option which may
01945      overrides some parts of the current output image geometry.
01946      This also overrides its default 'bestfit' setting.
01947   */
01948   { const char *artifact=GetImageArtifact(image,"distort:viewport");
01949     viewport_given = MagickFalse;
01950     if ( artifact != (const char *) NULL ) {
01951       (void) ParseAbsoluteGeometry(artifact,&geometry);
01952       viewport_given = MagickTrue;
01953     }
01954   }
01955 
01956   /* Verbose output */
01957   if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
01958     register ssize_t
01959        i;
01960     char image_gen[MaxTextExtent];
01961     const char *lookup;
01962 
01963     /* Set destination image size and virtual offset */
01964     if ( bestfit || viewport_given ) {
01965       (void) FormatLocaleString(image_gen, MaxTextExtent,"  -size %.20gx%.20g "
01966         "-page %+.20g%+.20g xc: +insert \\\n",(double) geometry.width,
01967         (double) geometry.height,(double) geometry.x,(double) geometry.y);
01968       lookup="v.p{ xx-v.page.x-.5, yy-v.page.y-.5 }";
01969     }
01970     else {
01971       image_gen[0] = '\0';             /* no destination to generate */
01972       lookup = "p{ xx-page.x-.5, yy-page.y-.5 }"; /* simplify lookup */
01973     }
01974 
01975     switch (method) {
01976       case AffineDistortion:
01977       {
01978         double *inverse;
01979 
01980         inverse = (double *) AcquireQuantumMemory(6,sizeof(*inverse));
01981         if (inverse == (double *) NULL) {
01982           coeff = (double *) RelinquishMagickMemory(coeff);
01983           (void) ThrowMagickException(exception,GetMagickModule(),
01984                   ResourceLimitError,"MemoryAllocationFailed",
01985                   "%s", "DistortImages");
01986           return((Image *) NULL);
01987         }
01988         InvertAffineCoefficients(coeff, inverse);
01989         CoefficientsToAffineArgs(inverse);
01990         (void) FormatLocaleFile(stderr, "Affine Projection:\n");
01991         (void) FormatLocaleFile(stderr, "  -distort AffineProjection \\\n      '");
01992         for (i=0; i < 5; i++)
01993           (void) FormatLocaleFile(stderr, "%lf,", inverse[i]);
01994         (void) FormatLocaleFile(stderr, "%lf'\n", inverse[5]);
01995         inverse = (double *) RelinquishMagickMemory(inverse);
01996 
01997         (void) FormatLocaleFile(stderr, "Affine Distort, FX Equivelent:\n");
01998         (void) FormatLocaleFile(stderr, "%s", image_gen);
01999         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
02000         (void) FormatLocaleFile(stderr, "       xx=%+lf*ii %+lf*jj %+lf;\n",
02001             coeff[0], coeff[1], coeff[2]);
02002         (void) FormatLocaleFile(stderr, "       yy=%+lf*ii %+lf*jj %+lf;\n",
02003             coeff[3], coeff[4], coeff[5]);
02004         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
02005 
02006         break;
02007       }
02008 
02009       case PerspectiveDistortion:
02010       {
02011         double *inverse;
02012 
02013         inverse = (double *) AcquireQuantumMemory(8,sizeof(*inverse));
02014         if (inverse == (double *) NULL) {
02015           coeff = (double *) RelinquishMagickMemory(coeff);
02016           (void) ThrowMagickException(exception,GetMagickModule(),
02017                   ResourceLimitError,"MemoryAllocationFailed",
02018                   "%s", "DistortCoefficients");
02019           return((Image *) NULL);
02020         }
02021         InvertPerspectiveCoefficients(coeff, inverse);
02022         (void) FormatLocaleFile(stderr, "Perspective Projection:\n");
02023         (void) FormatLocaleFile(stderr, "  -distort PerspectiveProjection \\\n      '");
02024         for (i=0; i<4; i++)
02025           (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
02026         (void) FormatLocaleFile(stderr, "\n       ");
02027         for (; i<7; i++)
02028           (void) FormatLocaleFile(stderr, "%lf, ", inverse[i]);
02029         (void) FormatLocaleFile(stderr, "%lf'\n", inverse[7]);
02030         inverse = (double *) RelinquishMagickMemory(inverse);
02031 
02032         (void) FormatLocaleFile(stderr, "Perspective Distort, FX Equivelent:\n");
02033         (void) FormatLocaleFile(stderr, "%s", image_gen);
02034         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
02035         (void) FormatLocaleFile(stderr, "       rr=%+lf*ii %+lf*jj + 1;\n",
02036             coeff[6], coeff[7]);
02037         (void) FormatLocaleFile(stderr, "       xx=(%+lf*ii %+lf*jj %+lf)/rr;\n",
02038             coeff[0], coeff[1], coeff[2]);
02039         (void) FormatLocaleFile(stderr, "       yy=(%+lf*ii %+lf*jj %+lf)/rr;\n",
02040             coeff[3], coeff[4], coeff[5]);
02041         (void) FormatLocaleFile(stderr, "       rr%s0 ? %s : blue' \\\n",
02042             coeff[8] < 0 ? "<" : ">", lookup);
02043         break;
02044       }
02045 
02046       case BilinearForwardDistortion:
02047         (void) FormatLocaleFile(stderr, "BilinearForward Mapping Equations:\n");
02048         (void) FormatLocaleFile(stderr, "%s", image_gen);
02049         (void) FormatLocaleFile(stderr, "    i = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
02050             coeff[0], coeff[1], coeff[2], coeff[3]);
02051         (void) FormatLocaleFile(stderr, "    j = %+lf*x %+lf*y %+lf*x*y %+lf;\n",
02052             coeff[4], coeff[5], coeff[6], coeff[7]);
02053 #if 0
02054         /* for debugging */
02055         (void) FormatLocaleFile(stderr, "   c8 = %+lf  c9 = 2*a = %+lf;\n",
02056             coeff[8], coeff[9]);
02057 #endif
02058         (void) FormatLocaleFile(stderr, "BilinearForward Distort, FX Equivelent:\n");
02059         (void) FormatLocaleFile(stderr, "%s", image_gen);
02060         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
02061             0.5-coeff[3], 0.5-coeff[7]);
02062         (void) FormatLocaleFile(stderr, "       bb=%lf*ii %+lf*jj %+lf;\n",
02063             coeff[6], -coeff[2], coeff[8]);
02064         /* Handle Special degenerate (non-quadratic) or trapezoidal case */
02065         if ( coeff[9] != 0 ) {
02066           (void) FormatLocaleFile(stderr, "       rt=bb*bb %+lf*(%lf*ii%+lf*jj);\n",
02067               -2*coeff[9],  coeff[4], -coeff[0]);
02068           (void) FormatLocaleFile(stderr, "       yy=( -bb + sqrt(rt) ) / %lf;\n",
02069                coeff[9]);
02070         } else
02071           (void) FormatLocaleFile(stderr, "       yy=(%lf*ii%+lf*jj)/bb;\n",
02072                 -coeff[4], coeff[0]);
02073         (void) FormatLocaleFile(stderr, "       xx=(ii %+lf*yy)/(%lf %+lf*yy);\n",
02074              -coeff[1], coeff[0], coeff[2]);
02075         if ( coeff[9] != 0 )
02076           (void) FormatLocaleFile(stderr, "       (rt < 0 ) ? red : %s'\n", lookup);
02077         else
02078           (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
02079         break;
02080 
02081       case BilinearReverseDistortion:
02082 #if 0
02083         (void) FormatLocaleFile(stderr, "Polynomial Projection Distort:\n");
02084         (void) FormatLocaleFile(stderr, "  -distort PolynomialProjection \\\n");
02085         (void) FormatLocaleFile(stderr, "      '1.5, %lf, %lf, %lf, %lf,\n",
02086             coeff[3], coeff[0], coeff[1], coeff[2]);
02087         (void) FormatLocaleFile(stderr, "            %lf, %lf, %lf, %lf'\n",
02088             coeff[7], coeff[4], coeff[5], coeff[6]);
02089 #endif
02090         (void) FormatLocaleFile(stderr, "BilinearReverse Distort, FX Equivelent:\n");
02091         (void) FormatLocaleFile(stderr, "%s", image_gen);
02092         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
02093         (void) FormatLocaleFile(stderr, "       xx=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
02094             coeff[0], coeff[1], coeff[2], coeff[3]);
02095         (void) FormatLocaleFile(stderr, "       yy=%+lf*ii %+lf*jj %+lf*ii*jj %+lf;\n",
02096             coeff[4], coeff[5], coeff[6], coeff[7]);
02097         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
02098         break;
02099 
02100       case PolynomialDistortion:
02101       {
02102         size_t nterms = (size_t) coeff[1];
02103         (void) FormatLocaleFile(stderr, "Polynomial (order %lg, terms %lu), FX Equivelent\n",
02104           coeff[0],(unsigned long) nterms);
02105         (void) FormatLocaleFile(stderr, "%s", image_gen);
02106         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x+0.5; jj=j+page.y+0.5;\n");
02107         (void) FormatLocaleFile(stderr, "       xx =");
02108         for (i=0; i<(ssize_t) nterms; i++) {
02109           if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n         ");
02110           (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i],
02111                poly_basis_str(i));
02112         }
02113         (void) FormatLocaleFile(stderr, ";\n       yy =");
02114         for (i=0; i<(ssize_t) nterms; i++) {
02115           if ( i != 0 && i%4 == 0 ) (void) FormatLocaleFile(stderr, "\n         ");
02116           (void) FormatLocaleFile(stderr, " %+lf%s", coeff[2+i+nterms],
02117                poly_basis_str(i));
02118         }
02119         (void) FormatLocaleFile(stderr, ";\n       %s' \\\n", lookup);
02120         break;
02121       }
02122       case ArcDistortion:
02123       {
02124         (void) FormatLocaleFile(stderr, "Arc Distort, Internal Coefficients:\n");
02125         for ( i=0; i<5; i++ )
02126           (void) FormatLocaleFile(stderr, "  c%.20g = %+lf\n", (double) i, coeff[i]);
02127         (void) FormatLocaleFile(stderr, "Arc Distort, FX Equivelent:\n");
02128         (void) FormatLocaleFile(stderr, "%s", image_gen);
02129         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x; jj=j+page.y;\n");
02130         (void) FormatLocaleFile(stderr, "       xx=(atan2(jj,ii)%+lf)/(2*pi);\n",
02131                                   -coeff[0]);
02132         (void) FormatLocaleFile(stderr, "       xx=xx-round(xx);\n");
02133         (void) FormatLocaleFile(stderr, "       xx=xx*%lf %+lf;\n",
02134                             coeff[1], coeff[4]);
02135         (void) FormatLocaleFile(stderr, "       yy=(%lf - hypot(ii,jj)) * %lf;\n",
02136                             coeff[2], coeff[3]);
02137         (void) FormatLocaleFile(stderr, "       v.p{xx-.5,yy-.5}' \\\n");
02138         break;
02139       }
02140       case PolarDistortion:
02141       {
02142         (void) FormatLocaleFile(stderr, "Polar Distort, Internal Coefficents\n");
02143         for ( i=0; i<8; i++ )
02144           (void) FormatLocaleFile(stderr, "  c%.20g = %+lf\n", (double) i, coeff[i]);
02145         (void) FormatLocaleFile(stderr, "Polar Distort, FX Equivelent:\n");
02146         (void) FormatLocaleFile(stderr, "%s", image_gen);
02147         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf; jj=j+page.y%+lf;\n",
02148                          -coeff[2], -coeff[3]);
02149         (void) FormatLocaleFile(stderr, "       xx=(atan2(ii,jj)%+lf)/(2*pi);\n",
02150                          -(coeff[4]+coeff[5])/2 );
02151         (void) FormatLocaleFile(stderr, "       xx=xx-round(xx);\n");
02152         (void) FormatLocaleFile(stderr, "       xx=xx*2*pi*%lf + v.w/2;\n",
02153                          coeff[6] );
02154         (void) FormatLocaleFile(stderr, "       yy=(hypot(ii,jj)%+lf)*%lf;\n",
02155                          -coeff[1], coeff[7] );
02156         (void) FormatLocaleFile(stderr, "       v.p{xx-.5,yy-.5}' \\\n");
02157         break;
02158       }
02159       case DePolarDistortion:
02160       {
02161         (void) FormatLocaleFile(stderr, "DePolar Distort, Internal Coefficents\n");
02162         for ( i=0; i<8; i++ )
02163           (void) FormatLocaleFile(stderr, "  c%.20g = %+lf\n", (double) i, coeff[i]);
02164         (void) FormatLocaleFile(stderr, "DePolar Distort, FX Equivelent:\n");
02165         (void) FormatLocaleFile(stderr, "%s", image_gen);
02166         (void) FormatLocaleFile(stderr, "  -fx 'aa=(i+.5)*%lf %+lf;\n", coeff[6], -coeff[4] );
02167         (void) FormatLocaleFile(stderr, "       rr=(j+.5)*%lf %+lf;\n", coeff[7], +coeff[1] );
02168         (void) FormatLocaleFile(stderr, "       xx=rr*sin(aa) %+lf;\n", coeff[2] );
02169         (void) FormatLocaleFile(stderr, "       yy=rr*cos(aa) %+lf;\n", coeff[3] );
02170         (void) FormatLocaleFile(stderr, "       v.p{xx-.5,yy-.5}' \\\n");
02171         break;
02172       }
02173       case Cylinder2PlaneDistortion:
02174       {
02175         (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, Internal Coefficents\n");
02176         (void) FormatLocaleFile(stderr, "  cylinder_radius = %+lf\n", coeff[1]);
02177         (void) FormatLocaleFile(stderr, "Cylinder to Plane Distort, FX Equivelent:\n");
02178         (void) FormatLocaleFile(stderr, "%s", image_gen);
02179         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
02180                          -coeff[4], -coeff[5]);
02181         (void) FormatLocaleFile(stderr, "       aa=atan(ii/%+lf);\n", coeff[1] );
02182         (void) FormatLocaleFile(stderr, "       xx=%lf*aa%+lf;\n",
02183                          coeff[1], coeff[2] );
02184         (void) FormatLocaleFile(stderr, "       yy=jj*cos(aa)%+lf;\n", coeff[3] );
02185         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
02186         break;
02187       }
02188       case Plane2CylinderDistortion:
02189       {
02190         (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, Internal Coefficents\n");
02191         (void) FormatLocaleFile(stderr, "  cylinder_radius = %+lf\n", coeff[1]);
02192         (void) FormatLocaleFile(stderr, "Plane to Cylinder Distort, FX Equivelent:\n");
02193         (void) FormatLocaleFile(stderr, "%s", image_gen);
02194         (void) FormatLocaleFile(stderr, "  -fx 'ii=i+page.x%+lf+0.5; jj=j+page.y%+lf+0.5;\n",
02195                          -coeff[4], -coeff[5]);
02196         (void) FormatLocaleFile(stderr, "       ii=ii/%+lf;\n", coeff[1] );
02197         (void) FormatLocaleFile(stderr, "       xx=%lf*tan(ii)%+lf;\n",
02198                          coeff[1], coeff[2] );
02199         (void) FormatLocaleFile(stderr, "       yy=jj/cos(ii)%+lf;\n",
02200                          coeff[3] );
02201         (void) FormatLocaleFile(stderr, "       %s' \\\n", lookup);
02202         break;
02203       }
02204       case BarrelDistortion:
02205       case BarrelInverseDistortion:
02206       { double xc,yc;
02207         /* NOTE: This does the barrel roll in pixel coords not image coords
02208         ** The internal distortion must do it in image coordinates,
02209         ** so that is what the center coeff (8,9) is given in.
02210         */
02211         xc = ((double)image->columns-1.0)/2.0 + image->page.x;
02212         yc = ((double)image->rows-1.0)/2.0    + image->page.y;
02213         (void) FormatLocaleFile(stderr, "Barrel%s Distort, FX Equivelent:\n",
02214              method == BarrelDistortion ? "" : "Inv");
02215         (void) FormatLocaleFile(stderr, "%s", image_gen);
02216         if ( fabs(coeff[8]-xc-0.5) < 0.1 && fabs(coeff[9]-yc-0.5) < 0.1 )
02217           (void) FormatLocaleFile(stderr, "  -fx 'xc=(w-1)/2;  yc=(h-1)/2;\n");
02218         else
02219           (void) FormatLocaleFile(stderr, "  -fx 'xc=%lf;  yc=%lf;\n",
02220                coeff[8]-0.5, coeff[9]-0.5);
02221         (void) FormatLocaleFile(stderr,
02222              "       ii=i-xc;  jj=j-yc;  rr=hypot(ii,jj);\n");
02223         (void) FormatLocaleFile(stderr, "       ii=ii%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
02224              method == BarrelDistortion ? "*" : "/",
02225              coeff[0],coeff[1],coeff[2],coeff[3]);
02226         (void) FormatLocaleFile(stderr, "       jj=jj%s(%lf*rr*rr*rr %+lf*rr*rr %+lf*rr %+lf);\n",
02227              method == BarrelDistortion ? "*" : "/",
02228              coeff[4],coeff[5],coeff[6],coeff[7]);
02229         (void) FormatLocaleFile(stderr, "       v.p{fx*ii+xc,fy*jj+yc}' \\\n");
02230       }
02231       default:
02232         break;
02233     }
02234   }
02235 
02236   /* The user provided a 'scale' expert option will scale the
02237      output image size, by the factor given allowing for super-sampling
02238      of the distorted image space.  Any scaling factors must naturally
02239      be halved as a result.
02240   */
02241   { const char *artifact;
02242     artifact=GetImageArtifact(image,"distort:scale");
02243     output_scaling = 1.0;
02244     if (artifact != (const char *) NULL) {
02245       output_scaling = fabs(StringToDouble(artifact,(char **) NULL));
02246       geometry.width  *= (size_t) output_scaling;
02247       geometry.height *= (size_t) output_scaling;
02248       geometry.x      *= (ssize_t) output_scaling;
02249       geometry.y      *= (ssize_t) output_scaling;
02250       if ( output_scaling < 0.1 ) {
02251         coeff = (double *) RelinquishMagickMemory(coeff);
02252         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
02253                 "InvalidArgument","%s", "-set option:distort:scale" );
02254         return((Image *) NULL);
02255       }
02256       output_scaling = 1/output_scaling;
02257     }
02258   }
02259 #define ScaleFilter(F,A,B,C,D) \
02260     ScaleResampleFilter( (F), \
02261       output_scaling*(A), output_scaling*(B), \
02262       output_scaling*(C), output_scaling*(D) )
02263 
02264   /*
02265     Initialize the distort image attributes.
02266   */
02267   distort_image=CloneImage(image,geometry.width,geometry.height,MagickTrue,
02268     exception);
02269   if (distort_image == (Image *) NULL)
02270     return((Image *) NULL);
02271   /* if image is ColorMapped - change it to DirectClass */
02272   if (SetImageStorageClass(distort_image,DirectClass,exception) == MagickFalse)
02273     {
02274       distort_image=DestroyImage(distort_image);
02275       return((Image *) NULL);
02276     }
02277   distort_image->page.x=geometry.x;
02278   distort_image->page.y=geometry.y;
02279   if (distort_image->background_color.matte != MagickFalse)
02280     distort_image->matte=MagickTrue;
02281 
02282   { /* ----- MAIN CODE -----
02283        Sample the source image to each pixel in the distort image.
02284      */
02285     CacheView
02286       *distort_view;
02287 
02288     MagickBooleanType
02289       status;
02290 
02291     MagickOffsetType
02292       progress;
02293 
02294     PixelInfo
02295       zero;
02296 
02297     ResampleFilter
02298       **restrict resample_filter;
02299 
02300     ssize_t
02301       j;
02302 
02303     status=MagickTrue;
02304     progress=0;
02305     GetPixelInfo(distort_image,&zero);
02306     resample_filter=AcquireResampleFilterThreadSet(image,
02307       UndefinedVirtualPixelMethod,MagickFalse,exception);
02308     distort_view=AcquireCacheView(distort_image);
02309 #if defined(MAGICKCORE_OPENMP_SUPPORT)
02310   #pragma omp parallel for schedule(static,4) shared(progress,status)
02311 #endif
02312     for (j=0; j < (ssize_t) distort_image->rows; j++)
02313     {
02314       const int
02315         id = GetOpenMPThreadId();
02316 
02317       double
02318         validity;  /* how mathematically valid is this the mapping */
02319 
02320       MagickBooleanType
02321         sync;
02322 
02323       PixelInfo
02324         pixel,    /* pixel color to assign to distorted image */
02325         invalid;  /* the color to assign when distort result is invalid */
02326 
02327       PointInfo
02328         d,
02329         s;  /* transform destination image x,y  to source image x,y */
02330 
02331       register ssize_t
02332         i;
02333 
02334       register Quantum
02335         *restrict q;
02336 
02337       q=QueueCacheViewAuthenticPixels(distort_view,0,j,distort_image->columns,1,
02338         exception);
02339       if (q == (Quantum *) NULL)
02340         {
02341           status=MagickFalse;
02342           continue;
02343         }
02344       pixel=zero;
02345 
02346       /* Define constant scaling vectors for Affine Distortions
02347         Other methods are either variable, or use interpolated lookup
02348       */
02349       switch (method)
02350       {
02351         case AffineDistortion:
02352           ScaleFilter( resample_filter[id],
02353             coeff[0], coeff[1],
02354             coeff[3], coeff[4] );
02355           break;
02356         default:
02357           break;
02358       }
02359 
02360       /* Initialize default pixel validity
02361       *    negative:         pixel is invalid  output 'matte_color'
02362       *    0.0 to 1.0:       antialiased, mix with resample output
02363       *    1.0 or greater:   use resampled output.
02364       */
02365       validity = 1.0;
02366 
02367       invalid=distort_image->matte_color;
02368       if (distort_image->colorspace == CMYKColorspace)
02369         ConvertRGBToCMYK(&invalid);   /* what about other color spaces? */
02370       for (i=0; i < (ssize_t) distort_image->columns; i++)
02371       {
02372         /* map pixel coordinate to distortion space coordinate */
02373         d.x = (double) (geometry.x+i+0.5)*output_scaling;
02374         d.y = (double) (geometry.y+j+0.5)*output_scaling;
02375         s = d;  /* default is a no-op mapping */
02376         switch (method)
02377         {
02378           case AffineDistortion:
02379           {
02380             s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
02381             s.y=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
02382             /* Affine partial derivitives are constant -- set above */
02383             break;
02384           }
02385           case PerspectiveDistortion:
02386           {
02387             double
02388               p,q,r,abs_r,abs_c6,abs_c7,scale;
02389             /* perspective is a ratio of affines */
02390             p=coeff[0]*d.x+coeff[1]*d.y+coeff[2];
02391             q=coeff[3]*d.x+coeff[4]*d.y+coeff[5];
02392             r=coeff[6]*d.x+coeff[7]*d.y+1.0;
02393             /* Pixel Validity -- is it a 'sky' or 'ground' pixel */
02394             validity = (r*coeff[8] < 0.0) ? 0.0 : 1.0;
02395             /* Determine horizon anti-alias blending */
02396             abs_r = fabs(r)*2;
02397             abs_c6 = fabs(coeff[6]);
02398             abs_c7 = fabs(coeff[7]);
02399             if ( abs_c6 > abs_c7 ) {
02400               if ( abs_r < abs_c6*output_scaling )
02401                 validity = 0.5 - coeff[8]*r/(coeff[6]*output_scaling);
02402             }
02403             else if ( abs_r < abs_c7*output_scaling )
02404               validity = 0.5 - coeff[8]*r/(coeff[7]*output_scaling);
02405             /* Perspective Sampling Point (if valid) */
02406             if ( validity > 0.0 ) {
02407               /* divide by r affine, for perspective scaling */
02408               scale = 1.0/r;
02409               s.x = p*scale;
02410               s.y = q*scale;
02411               /* Perspective Partial Derivatives or Scaling Vectors */
02412               scale *= scale;
02413               ScaleFilter( resample_filter[id],
02414                 (r*coeff[0] - p*coeff[6])*scale,
02415                 (r*coeff[1] - p*coeff[7])*scale,
02416                 (r*coeff[3] - q*coeff[6])*scale,
02417                 (r*coeff[4] - q*coeff[7])*scale );
02418             }
02419             break;
02420           }
02421           case BilinearReverseDistortion:
02422           {
02423             /* Reversed Mapped is just a simple polynomial */
02424             s.x=coeff[0]*d.x+coeff[1]*d.y+coeff[2]*d.x*d.y+coeff[3];
02425             s.y=coeff[4]*d.x+coeff[5]*d.y
02426                     +coeff[6]*d.x*d.y+coeff[7];
02427             /* Bilinear partial derivitives of scaling vectors */
02428             ScaleFilter( resample_filter[id],
02429                 coeff[0] + coeff[2]*d.y,
02430                 coeff[1] + coeff[2]*d.x,
02431                 coeff[4] + coeff[6]*d.y,
02432                 coeff[5] + coeff[6]*d.x );
02433             break;
02434           }
02435           case BilinearForwardDistortion:
02436           {
02437             /* Forward mapped needs reversed polynomial equations
02438              * which unfortunatally requires a square root!  */
02439             double b,c;
02440             d.x -= coeff[3];  d.y -= coeff[7];
02441             b = coeff[6]*d.x - coeff[2]*d.y + coeff[8];
02442             c = coeff[4]*d.x - coeff[0]*d.y;
02443 
02444             validity = 1.0;
02445             /* Handle Special degenerate (non-quadratic) case
02446              * Currently without horizon anti-alising */
02447             if ( fabs(coeff[9]) < MagickEpsilon )
02448               s.y =  -c/b;
02449             else {
02450               c = b*b - 2*coeff[9]*c;
02451               if ( c < 0.0 )
02452                 validity = 0.0;
02453               else
02454                 s.y = ( -b + sqrt(c) )/coeff[9];
02455             }
02456             if ( validity > 0.0 )
02457               s.x = ( d.x - coeff[1]*s.y) / ( coeff[0] + coeff[2]*s.y );
02458 
02459             /* NOTE: the sign of the square root should be -ve for parts
02460                      where the source image becomes 'flipped' or 'mirrored'.
02461                FUTURE: Horizon handling
02462                FUTURE: Scaling factors or Deritives (how?)
02463             */
02464             break;
02465           }
02466 #if 0
02467           case BilinearDistortion:
02468             /* Bilinear mapping of any Quadrilateral to any Quadrilateral */
02469             /* UNDER DEVELOPMENT */
02470             break;
02471 #endif
02472           case PolynomialDistortion:
02473           {
02474             /* multi-ordered polynomial */
02475             register ssize_t
02476               k;
02477 
02478             ssize_t
02479               nterms=(ssize_t)coeff[1];
02480 
02481             PointInfo
02482               du,dv; /* the du,dv vectors from unit dx,dy -- derivatives */
02483 
02484             s.x=s.y=du.x=du.y=dv.x=dv.y=0.0;
02485             for(k=0; k < nterms; k++) {
02486               s.x  += poly_basis_fn(k,d.x,d.y)*coeff[2+k];
02487               du.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k];
02488               du.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k];
02489               s.y  += poly_basis_fn(k,d.x,d.y)*coeff[2+k+nterms];
02490               dv.x += poly_basis_dx(k,d.x,d.y)*coeff[2+k+nterms];
02491               dv.y += poly_basis_dy(k,d.x,d.y)*coeff[2+k+nterms];
02492             }
02493             ScaleFilter( resample_filter[id], du.x,du.y,dv.x,dv.y );
02494             break;
02495           }
02496           case ArcDistortion:
02497           {
02498             /* what is the angle and radius in the destination image */
02499             s.x  = (double) ((atan2(d.y,d.x) - coeff[0])/Magick2PI);
02500             s.x -= MagickRound(s.x);     /* angle */
02501             s.y  = hypot(d.x,d.y);       /* radius */
02502 
02503             /* Arc Distortion Partial Scaling Vectors
02504               Are derived by mapping the perpendicular unit vectors
02505               dR  and  dA*R*2PI  rather than trying to map dx and dy
02506               The results is a very simple orthogonal aligned ellipse.
02507             */
02508             if ( s.y > MagickEpsilon )
02509               ScaleFilter( resample_filter[id],
02510                   (double) (coeff[1]/(Magick2PI*s.y)), 0, 0, coeff[3] );
02511             else
02512               ScaleFilter( resample_filter[id],
02513                   distort_image->columns*2, 0, 0, coeff[3] );
02514 
02515             /* now scale the angle and radius for source image lookup point */
02516             s.x = s.x*coeff[1] + coeff[4] + image->page.x +0.5;
02517             s.y = (coeff[2] - s.y) * coeff[3] + image->page.y;
02518             break;
02519           }
02520           case PolarDistortion:
02521           { /* 2D Cartesain to Polar View */
02522             d.x -= coeff[2];
02523             d.y -= coeff[3];
02524             s.x  = atan2(d.x,d.y) - (coeff[4]+coeff[5])/2;
02525             s.x /= Magick2PI;
02526             s.x -= MagickRound(s.x);
02527             s.x *= Magick2PI;       /* angle - relative to centerline */
02528             s.y  = hypot(d.x,d.y);  /* radius */
02529 
02530             /* Polar Scaling vectors are based on mapping dR and dA vectors
02531                This results in very simple orthogonal scaling vectors
02532             */
02533             if ( s.y > MagickEpsilon )
02534               ScaleFilter( resample_filter[id],
02535                 (double) (coeff[6]/(Magick2PI*s.y)), 0, 0, coeff[7] );
02536             else
02537               ScaleFilter( resample_filter[id],
02538                   distort_image->columns*2, 0, 0, coeff[7] );
02539 
02540             /* now finish mapping radius/angle to source x,y coords */
02541             s.x = s.x*coeff[6] + (double)image->columns/2.0 + image->page.x;
02542             s.y = (s.y-coeff[1])*coeff[7] + image->page.y;
02543             break;
02544           }
02545           case DePolarDistortion:
02546           { /* @D Polar to Carteasain  */
02547             /* ignore all destination virtual offsets */
02548             d.x = ((double)i+0.5)*output_scaling*coeff[6]-coeff[4];
02549             d.y = ((double)j+0.5)*output_scaling*coeff[7]+coeff[1];
02550             s.x = d.y*sin(d.x) + coeff[2];
02551             s.y = d.y*cos(d.x) + coeff[3];
02552             /* derivatives are usless - better to use SuperSampling */
02553             break;
02554           }
02555           case Cylinder2PlaneDistortion:
02556           { /* 3D Cylinder to Tangential Plane */
02557             double ax, cx;
02558             /* relative to center of distortion */
02559             d.x -= coeff[4]; d.y -= coeff[5];
02560             d.x /= coeff[1];        /* x' = x/r */
02561             ax=atan(d.x);           /* aa = atan(x/r) = u/r  */
02562             cx=cos(ax);             /* cx = cos(atan(x/r)) = 1/sqrt(x^2+u^2) */
02563             s.x = coeff[1]*ax;      /* u  = r*atan(x/r) */
02564             s.y = d.y*cx;           /* v  = y*cos(u/r) */
02565             /* derivatives... (see personnal notes) */
02566             ScaleFilter( resample_filter[id],
02567                   1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
02568 #if 0
02569 if ( i == 0 && j == 0 ) {
02570   fprintf(stderr, "x=%lf  y=%lf  u=%lf  v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
02571   fprintf(stderr, "phi = %lf\n", (double)(ax * 180.0/MagickPI) );
02572   fprintf(stderr, "du/dx=%lf  du/dx=%lf  dv/dx=%lf  dv/dy=%lf\n",
02573                 1.0/(1.0+d.x*d.x), 0.0, -d.x*s.y*cx*cx/coeff[1], s.y/d.y );
02574   fflush(stderr); }
02575 #endif
02576             /* add center of distortion in source */
02577             s.x += coeff[2]; s.y += coeff[3];
02578             break;
02579           }
02580           case Plane2CylinderDistortion:
02581           { /* 3D Cylinder to Tangential Plane */
02582             /* relative to center of distortion */
02583             d.x -= coeff[4]; d.y -= coeff[5];
02584 
02585             /* is pixel valid - horizon of a infinite Virtual-Pixel Plane
02586              * (see Anthony Thyssen's personal note) */
02587             validity = (double) (coeff[1]*MagickPI2 - fabs(d.x))/output_scaling + 0.5;
02588 
02589             if ( validity > 0.0 ) {
02590               double cx,tx;
02591               d.x /= coeff[1];           /* x'= x/r */
02592               cx = 1/cos(d.x);           /* cx = 1/cos(x/r) */
02593               tx = tan(d.x);             /* tx = tan(x/r) */
02594               s.x = coeff[1]*tx;         /* u = r * tan(x/r) */
02595               s.y = d.y*cx;              /* v = y / cos(x/r) */
02596               /* derivatives...  (see Anthony Thyssen's personal notes) */
02597               ScaleFilter( resample_filter[id],
02598                     cx*cx, 0.0, s.y*cx/coeff[1], cx );
02599 #if 1
02600 /*if ( i == 0 && j == 0 )*/
02601 if ( d.x == 0.5 && d.y == 0.5 ) {
02602   fprintf(stderr, "x=%lf  y=%lf  u=%lf  v=%lf\n", d.x*coeff[1], d.y, s.x, s.y);
02603   fprintf(stderr, "radius = %lf  phi = %lf  validity = %lf\n",
02604       coeff[1],  (double)(d.x * 180.0/MagickPI), validity );
02605   fprintf(stderr, "du/dx=%lf  du/dx=%lf  dv/dx=%lf  dv/dy=%lf\n",
02606       cx*cx, 0.0, s.y*cx/coeff[1], cx);
02607   fflush(stderr); }
02608 #endif
02609             }
02610             /* add center of distortion in source */
02611             s.x += coeff[2]; s.y += coeff[3];
02612             break;
02613           }
02614           case BarrelDistortion:
02615           case BarrelInverseDistortion:
02616           { /* Lens Barrel Distionion Correction */
02617             double r,fx,fy,gx,gy;
02618             /* Radial Polynomial Distortion (de-normalized) */
02619             d.x -= coeff[8];
02620             d.y -= coeff[9];
02621             r = sqrt(d.x*d.x+d.y*d.y);
02622             if ( r > MagickEpsilon ) {
02623               fx = ((coeff[0]*r + coeff[1])*r + coeff[2])*r + coeff[3];
02624               fy = ((coeff[4]*r + coeff[5])*r + coeff[6])*r + coeff[7];
02625               gx = ((3*coeff[0]*r + 2*coeff[1])*r + coeff[2])/r;
02626               gy = ((3*coeff[4]*r + 2*coeff[5])*r + coeff[6])/r;
02627               /* adjust functions and scaling for 'inverse' form */
02628               if ( method == BarrelInverseDistortion ) {
02629                 fx = 1/fx;  fy = 1/fy;
02630                 gx *= -fx*fx;  gy *= -fy*fy;
02631               }
02632               /* Set the source pixel to lookup and EWA derivative vectors */
02633               s.x = d.x*fx + coeff[8];
02634               s.y = d.y*fy + coeff[9];
02635               ScaleFilter( resample_filter[id],
02636                   gx*d.x*d.x + fx, gx*d.x*d.y,
02637                   gy*d.x*d.y,      gy*d.y*d.y + fy );
02638             }
02639             else {
02640               /* Special handling to avoid divide by zero when r==0
02641               **
02642               ** The source and destination pixels match in this case
02643               ** which was set at the top of the loop using  s = d;
02644               ** otherwise...   s.x=coeff[8]; s.y=coeff[9];
02645               */
02646               if ( method == BarrelDistortion )
02647                 ScaleFilter( resample_filter[id],
02648                      coeff[3], 0, 0, coeff[7] );
02649               else /* method == BarrelInverseDistortion */
02650                 /* FUTURE, trap for D==0 causing division by zero */
02651                 ScaleFilter( resample_filter[id],
02652                      1.0/coeff[3], 0, 0, 1.0/coeff[7] );
02653             }
02654             break;
02655           }
02656           case ShepardsDistortion:
02657           { /* Shepards Method, or Inverse Weighted Distance for
02658               displacement around the destination image control points
02659               The input arguments are the coefficents to the function.
02660               This is more of a 'displacement' function rather than an
02661               absolute distortion function.
02662             */
02663             size_t
02664               i;
02665             double
02666               denominator;
02667 
02668             denominator = s.x = s.y = 0;
02669             for(i=0; i<number_arguments; i+=4) {
02670               double weight =
02671                   ((double)d.x-arguments[i+2])*((double)d.x-arguments[i+2])
02672                 + ((double)d.y-arguments[i+3])*((double)d.y-arguments[i+3]);
02673               if ( weight != 0 )
02674                 weight = 1/weight;
02675               else
02676                 weight = 1;
02677 
02678               s.x += (arguments[ i ]-arguments[i+2])*weight;
02679               s.y += (arguments[i+1]-arguments[i+3])*weight;
02680               denominator += weight;
02681             }
02682             s.x /= denominator;
02683             s.y /= denominator;
02684             s.x += d.x;
02685             s.y += d.y;
02686 
02687             /* We can not determine derivatives using shepards method
02688                only color interpolatation, not area-resampling */
02689             break;
02690           }
02691           default:
02692             break; /* use the default no-op given above */
02693         }
02694         /* map virtual canvas location back to real image coordinate */
02695         if ( bestfit && method != ArcDistortion ) {
02696           s.x -= image->page.x;
02697           s.y -= image->page.y;
02698         }
02699         s.x -= 0.5;
02700         s.y -= 0.5;
02701 
02702         if ( validity <= 0.0 ) {
02703           /* result of distortion is an invalid pixel - don't resample */
02704           SetPixelInfoPixel(distort_image,&invalid,q);
02705         }
02706         else {
02707           /* resample the source image to find its correct color */
02708           (void) ResamplePixelColor(resample_filter[id],s.x,s.y,&pixel);
02709           /* if validity between 0.0 and 1.0 mix result with invalid pixel */
02710           if ( validity < 1.0 ) {
02711             /* Do a blend of sample color and invalid pixel */
02712             /* should this be a 'Blend', or an 'Over' compose */
02713             CompositePixelInfoBlend(&pixel,validity,&invalid,(1.0-validity),
02714               &pixel);
02715           }
02716           SetPixelInfoPixel(distort_image,&pixel,q);
02717         }
02718         q+=GetPixelChannels(distort_image);
02719       }
02720       sync=SyncCacheViewAuthenticPixels(distort_view,exception);
02721       if (sync == MagickFalse)
02722         status=MagickFalse;
02723       if (image->progress_monitor != (MagickProgressMonitor) NULL)
02724         {
02725           MagickBooleanType
02726             proceed;
02727 
02728 #if defined(MAGICKCORE_OPENMP_SUPPORT)
02729   #pragma omp critical (MagickCore_DistortImage)
02730 #endif
02731           proceed=SetImageProgress(image,DistortImageTag,progress++,
02732             image->rows);
02733           if (proceed == MagickFalse)
02734             status=MagickFalse;
02735         }
02736     }
02737     distort_view=DestroyCacheView(distort_view);
02738     resample_filter=DestroyResampleFilterThreadSet(resample_filter);
02739 
02740     if (status == MagickFalse)
02741       distort_image=DestroyImage(distort_image);
02742   }
02743 
02744   /* Arc does not return an offset unless 'bestfit' is in effect
02745      And the user has not provided an overriding 'viewport'.
02746    */
02747   if ( method == ArcDistortion && !bestfit && !viewport_given ) {
02748     distort_image->page.x = 0;
02749     distort_image->page.y = 0;
02750   }
02751   coeff = (double *) RelinquishMagickMemory(coeff);
02752   return(distort_image);
02753 }
02754 
02755 /*
02756 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02757 %                                                                             %
02758 %                                                                             %
02759 %                                                                             %
02760 %   R o t a t e I m a g e                                                     %
02761 %                                                                             %
02762 %                                                                             %
02763 %                                                                             %
02764 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02765 %
02766 %  RotateImage() creates a new image that is a rotated copy of an existing
02767 %  one.  Positive angles rotate counter-clockwise (right-hand rule), while
02768 %  negative angles rotate clockwise.  Rotated images are usually larger than
02769 %  the originals and have 'empty' triangular corners.  X axis.  Empty
02770 %  triangles left over from shearing the image are filled with the background
02771 %  color defined by member 'background_color' of the image.  RotateImage
02772 %  allocates the memory necessary for the new Image structure and returns a
02773 %  pointer to the new image.
02774 %
02775 %  The format of the RotateImage method is:
02776 %
02777 %      Image *RotateImage(const Image *image,const double degrees,
02778 %        ExceptionInfo *exception)
02779 %
02780 %  A description of each parameter follows.
02781 %
02782 %    o image: the image.
02783 %
02784 %    o degrees: Specifies the number of degrees to rotate the image.
02785 %
02786 %    o exception: return any errors or warnings in this structure.
02787 %
02788 */
02789 MagickExport Image *RotateImage(const Image *image,const double degrees,
02790   ExceptionInfo *exception)
02791 {
02792   Image
02793     *distort_image,
02794     *rotate_image;
02795 
02796   MagickRealType
02797     angle;
02798 
02799   PointInfo
02800     shear;
02801 
02802   size_t
02803     rotations;
02804 
02805   /*
02806     Adjust rotation angle.
02807   */
02808   assert(image != (Image *) NULL);
02809   assert(image->signature == MagickSignature);
02810   if (image->debug != MagickFalse)
02811     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
02812   assert(exception != (ExceptionInfo *) NULL);
02813   assert(exception->signature == MagickSignature);
02814   angle=degrees;
02815   while (angle < -45.0)
02816     angle+=360.0;
02817   for (rotations=0; angle > 45.0; rotations++)
02818     angle-=90.0;
02819   rotations%=4;
02820   shear.x=(-tan((double) DegreesToRadians(angle)/2.0));
02821   shear.y=sin((double) DegreesToRadians(angle));
02822   if ((fabs(shear.x) < MagickEpsilon) && (fabs(shear.y) < MagickEpsilon))
02823     return(IntegralRotateImage(image,rotations,exception));
02824   distort_image=CloneImage(image,0,0,MagickTrue,exception);
02825   if (distort_image == (Image *) NULL)
02826     return((Image *) NULL);
02827   (void) SetImageVirtualPixelMethod(distort_image,BackgroundVirtualPixelMethod,
02828     exception);
02829   rotate_image=DistortImage(distort_image,ScaleRotateTranslateDistortion,1,
02830     &degrees,MagickTrue,exception);
02831   distort_image=DestroyImage(distort_image);
02832   return(rotate_image);
02833 }
02834 
02835 /*
02836 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02837 %                                                                             %
02838 %                                                                             %
02839 %                                                                             %
02840 %   S p a r s e C o l o r I m a g e                                           %
02841 %                                                                             %
02842 %                                                                             %
02843 %                                                                             %
02844 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
02845 %
02846 %  SparseColorImage(), given a set of coordinates, interpolates the colors
02847 %  found at those coordinates, across the whole image, using various methods.
02848 %
02849 %  The format of the SparseColorImage() method is:
02850 %
02851 %      Image *SparseColorImage(const Image *image,
02852 %        const SparseColorMethod method,const size_t number_arguments,
02853 %        const double *arguments,ExceptionInfo *exception)
02854 %
02855 %  A description of each parameter follows:
02856 %
02857 %    o image: the image to be filled in.
02858 %
02859 %    o method: the method to fill in the gradient between the control points.
02860 %
02861 %        The methods used for SparseColor() are often simular to methods
02862 %        used for DistortImage(), and even share the same code for determination
02863 %        of the function coefficents, though with more dimensions (or resulting
02864 %        values).
02865 %
02866 %    o number_arguments: the number of arguments given.
02867 %
02868 %    o arguments: array of floating point arguments for this method--
02869 %        x,y,color_values-- with color_values given as normalized values.
02870 %
02871 %    o exception: return any errors or warnings in this structure
02872 %
02873 */
02874 MagickExport Image *SparseColorImage(const Image *image,
02875   const SparseColorMethod method,const size_t number_arguments,
02876   const double *arguments,ExceptionInfo *exception)
02877 {
02878 #define SparseColorTag  "Distort/SparseColor"
02879 
02880   SparseColorMethod
02881     sparse_method;
02882 
02883   double
02884     *coeff;
02885 
02886   Image
02887     *sparse_image;
02888 
02889   size_t
02890     number_colors;
02891 
02892   assert(image != (Image *) NULL);
02893   assert(image->signature == MagickSignature);
02894   if (image->debug != MagickFalse)
02895     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
02896   assert(exception != (ExceptionInfo *) NULL);
02897   assert(exception->signature == MagickSignature);
02898 
02899   /* Determine number of color values needed per control point */
02900   number_colors=0;
02901   if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
02902     number_colors++;
02903   if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
02904     number_colors++;
02905   if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
02906     number_colors++;
02907   if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
02908       (image->colorspace == CMYKColorspace))
02909     number_colors++;
02910   if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
02911       (image->matte != MagickFalse))
02912     number_colors++;
02913 
02914   /*
02915     Convert input arguments into mapping coefficients, this this case
02916     we are mapping (distorting) colors, rather than coordinates.
02917   */
02918   { DistortImageMethod
02919       distort_method;
02920 
02921     distort_method=(DistortImageMethod) method;
02922     if ( distort_method >= SentinelDistortion )
02923       distort_method = ShepardsDistortion; /* Pretend to be Shepards */
02924     coeff = GenerateCoefficients(image, &distort_method, number_arguments,
02925                 arguments, number_colors, exception);
02926     if ( coeff == (double *) NULL )
02927       return((Image *) NULL);
02928     /*
02929       Note some Distort Methods may fall back to other simpler methods,
02930       Currently the only fallback of concern is Bilinear to Affine
02931       (Barycentric), which is alaso sparse_colr method.  This also ensures
02932       correct two and one color Barycentric handling.
02933     */
02934     sparse_method = (SparseColorMethod) distort_method;
02935     if ( distort_method == ShepardsDistortion )
02936       sparse_method = method;   /* return non-distiort methods to normal */
02937   }
02938 
02939   /* Verbose output */
02940   if ( GetImageArtifact(image,"verbose") != (const char *) NULL ) {
02941 
02942     switch (sparse_method) {
02943       case BarycentricColorInterpolate:
02944       {
02945         register ssize_t x=0;
02946         (void) FormatLocaleFile(stderr, "Barycentric Sparse Color:\n");
02947         if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
02948           (void) FormatLocaleFile(stderr, "  -channel R -fx '%+lf*i %+lf*j %+lf' \\\n",
02949               coeff[x], coeff[x+1], coeff[x+2]),x+=3;
02950         if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
02951           (void) FormatLocaleFile(stderr, "  -channel G -fx '%+lf*i %+lf*j %+lf' \\\n",
02952               coeff[x], coeff[x+1], coeff[x+2]),x+=3;
02953         if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
02954           (void) FormatLocaleFile(stderr, "  -channel B -fx '%+lf*i %+lf*j %+lf' \\\n",
02955               coeff[x], coeff[x+1], coeff[x+2]),x+=3;
02956         if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
02957             (image->colorspace == CMYKColorspace))
02958           (void) FormatLocaleFile(stderr, "  -channel K -fx '%+lf*i %+lf*j %+lf' \\\n",
02959               coeff[x], coeff[x+1], coeff[x+2]),x+=3;
02960         if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
02961             (image->matte != MagickFalse))
02962           (void) FormatLocaleFile(stderr, "  -channel A -fx '%+lf*i %+lf*j %+lf' \\\n",
02963               coeff[x], coeff[x+1], coeff[x+2]),x+=3;
02964         break;
02965       }
02966       case BilinearColorInterpolate:
02967       {
02968         register ssize_t x=0;
02969         (void) FormatLocaleFile(stderr, "Bilinear Sparse Color\n");
02970         if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
02971           (void) FormatLocaleFile(stderr, "   -channel R -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
02972               coeff[ x ], coeff[x+1],
02973               coeff[x+2], coeff[x+3]),x+=4;
02974         if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
02975           (void) FormatLocaleFile(stderr, "   -channel G -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
02976               coeff[ x ], coeff[x+1],
02977               coeff[x+2], coeff[x+3]),x+=4;
02978         if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
02979           (void) FormatLocaleFile(stderr, "   -channel B -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
02980               coeff[ x ], coeff[x+1],
02981               coeff[x+2], coeff[x+3]),x+=4;
02982         if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
02983             (image->colorspace == CMYKColorspace))
02984           (void) FormatLocaleFile(stderr, "   -channel K -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
02985               coeff[ x ], coeff[x+1],
02986               coeff[x+2], coeff[x+3]),x+=4;
02987         if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
02988             (image->matte != MagickFalse))
02989           (void) FormatLocaleFile(stderr, "   -channel A -fx '%+lf*i %+lf*j %+lf*i*j %+lf;\n",
02990               coeff[ x ], coeff[x+1],
02991               coeff[x+2], coeff[x+3]),x+=4;
02992         break;
02993       }
02994       default:
02995         /* sparse color method is too complex for FX emulation */
02996         break;
02997     }
02998   }
02999 
03000   /* Generate new image for generated interpolated gradient.
03001    * ASIDE: Actually we could have just replaced the colors of the original
03002    * image, but IM Core policy, is if storage class could change then clone
03003    * the image.
03004    */
03005 
03006   sparse_image=CloneImage(image,0,0,MagickTrue,exception);
03007   if (sparse_image == (Image *) NULL)
03008     return((Image *) NULL);
03009   if (SetImageStorageClass(sparse_image,DirectClass,exception) == MagickFalse)
03010     { /* if image is ColorMapped - change it to DirectClass */
03011       sparse_image=DestroyImage(sparse_image);
03012       return((Image *) NULL);
03013     }
03014   { /* ----- MAIN CODE ----- */
03015     CacheView
03016       *sparse_view;
03017 
03018     MagickBooleanType
03019       status;
03020 
03021     MagickOffsetType
03022       progress;
03023 
03024     ssize_t
03025       j;
03026 
03027     status=MagickTrue;
03028     progress=0;
03029     sparse_view=AcquireCacheView(sparse_image);
03030 #if defined(MAGICKCORE_OPENMP_SUPPORT)
03031   #pragma omp parallel for schedule(static,4) shared(progress,status)
03032 #endif
03033     for (j=0; j < (ssize_t) sparse_image->rows; j++)
03034     {
03035       MagickBooleanType
03036         sync;
03037 
03038       PixelInfo
03039         pixel;    /* pixel to assign to distorted image */
03040 
03041       register ssize_t
03042         i;
03043 
03044       register Quantum
03045         *restrict q;
03046 
03047       q=GetCacheViewAuthenticPixels(sparse_view,0,j,sparse_image->columns,
03048         1,exception);
03049       if (q == (Quantum *) NULL)
03050         {
03051           status=MagickFalse;
03052           continue;
03053         }
03054       GetPixelInfo(sparse_image,&pixel);
03055       for (i=0; i < (ssize_t) image->columns; i++)
03056       {
03057         GetPixelInfoPixel(image,q,&pixel);
03058         switch (sparse_method)
03059         {
03060           case BarycentricColorInterpolate:
03061           {
03062             register ssize_t x=0;
03063             if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
03064               pixel.red     = coeff[x]*i +coeff[x+1]*j
03065                               +coeff[x+2], x+=3;
03066             if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
03067               pixel.green   = coeff[x]*i +coeff[x+1]*j
03068                               +coeff[x+2], x+=3;
03069             if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
03070               pixel.blue    = coeff[x]*i +coeff[x+1]*j
03071                               +coeff[x+2], x+=3;
03072             if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
03073                 (image->colorspace == CMYKColorspace))
03074               pixel.black   = coeff[x]*i +coeff[x+1]*j
03075                               +coeff[x+2], x+=3;
03076             if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
03077                 (image->matte != MagickFalse))
03078               pixel.alpha = coeff[x]*i +coeff[x+1]*j
03079                               +coeff[x+2], x+=3;
03080             break;
03081           }
03082           case BilinearColorInterpolate:
03083           {
03084             register ssize_t x=0;
03085             if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
03086               pixel.red     = coeff[x]*i     + coeff[x+1]*j +
03087                               coeff[x+2]*i*j + coeff[x+3], x+=4;
03088             if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
03089               pixel.green   = coeff[x]*i     + coeff[x+1]*j +
03090                               coeff[x+2]*i*j + coeff[x+3], x+=4;
03091             if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
03092               pixel.blue    = coeff[x]*i     + coeff[x+1]*j +
03093                               coeff[x+2]*i*j + coeff[x+3], x+=4;
03094             if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
03095                 (image->colorspace == CMYKColorspace))
03096               pixel.black   = coeff[x]*i     + coeff[x+1]*j +
03097                               coeff[x+2]*i*j + coeff[x+3], x+=4;
03098             if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
03099                 (image->matte != MagickFalse))
03100               pixel.alpha = coeff[x]*i     + coeff[x+1]*j +
03101                               coeff[x+2]*i*j + coeff[x+3], x+=4;
03102             break;
03103           }
03104           case InverseColorInterpolate:
03105           case ShepardsColorInterpolate:
03106           { /* Inverse (Squared) Distance weights average (IDW) */
03107             size_t
03108               k;
03109             double
03110               denominator;
03111 
03112             if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
03113               pixel.red=0.0;
03114             if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
03115               pixel.green=0.0;
03116             if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
03117               pixel.blue=0.0;
03118             if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
03119                 (image->colorspace == CMYKColorspace))
03120               pixel.black=0.0;
03121             if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
03122                 (image->matte != MagickFalse))
03123               pixel.alpha=0.0;
03124             denominator = 0.0;
03125             for(k=0; k<number_arguments; k+=2+number_colors) {
03126               register ssize_t x=(ssize_t) k+2;
03127               double weight =
03128                   ((double)i-arguments[ k ])*((double)i-arguments[ k ])
03129                 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
03130               if ( method == InverseColorInterpolate )
03131                 weight = sqrt(weight);  /* inverse, not inverse squared */
03132               weight = ( weight < 1.0 ) ? 1.0 : 1.0/weight;
03133               if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
03134                 pixel.red     += arguments[x++]*weight;
03135               if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
03136                 pixel.green   += arguments[x++]*weight;
03137               if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
03138                 pixel.blue    += arguments[x++]*weight;
03139               if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
03140                   (image->colorspace == CMYKColorspace))
03141                 pixel.black   += arguments[x++]*weight;
03142               if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
03143                   (image->matte != MagickFalse))
03144                 pixel.alpha += arguments[x++]*weight;
03145               denominator += weight;
03146             }
03147             if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
03148               pixel.red/=denominator;
03149             if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
03150               pixel.green/=denominator;
03151             if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
03152               pixel.blue/=denominator;
03153             if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
03154                 (image->colorspace == CMYKColorspace))
03155               pixel.black/=denominator;
03156             if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
03157                 (image->matte != MagickFalse))
03158               pixel.alpha/=denominator;
03159             break;
03160           }
03161           case VoronoiColorInterpolate:
03162           default:
03163           { /* Just use the closest control point you can find! */
03164             size_t
03165               k;
03166             double
03167               minimum = MagickHuge;
03168 
03169             for(k=0; k<number_arguments; k+=2+number_colors) {
03170               double distance =
03171                   ((double)i-arguments[ k ])*((double)i-arguments[ k ])
03172                 + ((double)j-arguments[k+1])*((double)j-arguments[k+1]);
03173               if ( distance < minimum ) {
03174                 register ssize_t x=(ssize_t) k+2;
03175                 if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
03176                   pixel.red=arguments[x++];
03177                 if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
03178                   pixel.green=arguments[x++];
03179                 if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
03180                   pixel.blue=arguments[x++];
03181                 if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
03182                     (image->colorspace == CMYKColorspace))
03183                   pixel.black=arguments[x++];
03184                 if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
03185                     (image->matte != MagickFalse))
03186                   pixel.alpha=arguments[x++];
03187                 minimum = distance;
03188               }
03189             }
03190             break;
03191           }
03192         }
03193         /* set the color directly back into the source image */
03194         if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
03195           pixel.red*=QuantumRange;
03196         if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
03197           pixel.green*=QuantumRange;
03198         if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
03199           pixel.blue*=QuantumRange;
03200         if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
03201             (image->colorspace == CMYKColorspace))
03202           pixel.black*=QuantumRange;
03203         if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
03204             (image->matte != MagickFalse))
03205           pixel.alpha*=QuantumRange;
03206         SetPixelInfoPixel(sparse_image,&pixel,q);
03207         q+=GetPixelChannels(sparse_image);
03208       }
03209       sync=SyncCacheViewAuthenticPixels(sparse_view,exception);
03210       if (sync == MagickFalse)
03211         status=MagickFalse;
03212       if (image->progress_monitor != (MagickProgressMonitor) NULL)
03213         {
03214           MagickBooleanType
03215             proceed;
03216 
03217 #if defined(MAGICKCORE_OPENMP_SUPPORT)
03218   #pragma omp critical (MagickCore_SparseColorImage)
03219 #endif
03220           proceed=SetImageProgress(image,SparseColorTag,progress++,image->rows);
03221           if (proceed == MagickFalse)
03222             status=MagickFalse;
03223         }
03224     }
03225     sparse_view=DestroyCacheView(sparse_view);
03226     if (status == MagickFalse)
03227       sparse_image=DestroyImage(sparse_image);
03228   }
03229   coeff = (double *) RelinquishMagickMemory(coeff);
03230   return(sparse_image);
03231 }