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