root / ImageMagick / trunk / magick / distort.c

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%               DDDD   IIIII  SSSSS  TTTTT   OOO   RRRR   TTTTT               %
%               D   D    I    SS       T    O   O  R   R    T                 %
%               D   D    I     SSS     T    O   O  RRRR     T                 %
%               D   D    I       SS    T    O   O  R R      T                 %
%               DDDD   IIIII  SSSSS    T     OOO   R  R     T                 %
%                                                                             %
%                                                                             %
%                     ImageMagick Image Distortion Methods.                   %
%                                                                             %
%                              Software Design                                %
%                                John Cristy                                  %
%                              Anthony Thyssen                                %
%                                 June 2007                                   %
%                                                                             %
%                                                                             %
%  Copyright 1999-2008 ImageMagick Studio LLC, a non-profit organization      %
%  dedicated to making software imaging solutions freely available.           %
%                                                                             %
%  You may not use this file except in compliance with the License.  You may  %
%  obtain a copy of the License at                                            %
%                                                                             %
%    http://www.imagemagick.org/script/license.php                            %
%                                                                             %
%  Unless required by applicable law or agreed to in writing, software        %
%  distributed under the License is distributed on an "AS IS" BASIS,          %
%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
%  See the License for the specific language governing permissions and        %
%  limitations under the License.                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
*/

/*
  Include declarations.
*/
#include "magick/studio.h"
#include "magick/artifact.h"
#include "magick/cache-view.h"
#include "magick/colorspace-private.h"
#include "magick/composite-private.h"
#include "magick/distort.h"
#include "magick/exception.h"
#include "magick/exception-private.h"
#include "magick/gem.h"
#include "magick/hashmap.h"
#include "magick/image.h"
#include "magick/list.h"
#include "magick/matrix.h"
#include "magick/memory_.h"
#include "magick/monitor-private.h"
#include "magick/pixel.h"
#include "magick/pixel-private.h"
#include "magick/resample.h"
#include "magick/registry.h"
#include "magick/semaphore.h"
#include "magick/string_.h"
#include "magick/token.h"

/*
  Numerous internal routines for image distortions.
*/
static inline double MagickMin(const double x,const double y)
{
  return( x < y ? x : y);
}
static inline double MagickMax(const double x,const double y)
{
  return( x > y ? x : y);
}

static inline void AffineArgsToCoefficients(double *affine)
{
  /* map  external sx,ry,rx,sy,tx,ty  to  internal c0,c2,c4,c1,c3,c5 */
  double tmp[4];  /* note indexes  0 and 5 remain unchanged */
  tmp[0]=affine[1]; tmp[1]=affine[2]; tmp[2]=affine[3]; tmp[3]=affine[4];
  affine[3]=tmp[0]; affine[1]=tmp[1]; affine[4]=tmp[2]; affine[2]=tmp[3];
}
static inline void CoefficientsToAffineArgs(double *coeff)
{
  /* map  internal c0,c1,c2,c3,c4,c5  to  external sx,ry,rx,sy,tx,ty */
  double tmp[4];  /* note indexes 0 and 5 remain unchanged */
  tmp[0]=coeff[3]; tmp[1]=coeff[1]; tmp[2]=coeff[4]; tmp[3]=coeff[2];
  coeff[1]=tmp[0]; coeff[2]=tmp[1]; coeff[3]=tmp[2]; coeff[4]=tmp[3];
}
static void InvertAffineCoefficients(const double *coeff,double *inverse)
{
  /* From "Digital Image Warping" by George Wolberg, page 50 */
  double determinant;

  determinant=1.0/(coeff[0]*coeff[4]-coeff[1]*coeff[3]);
  inverse[0]=determinant*coeff[4];
  inverse[1]=determinant*(-coeff[1]);
  inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[2]*coeff[4]);
  inverse[3]=determinant*(-coeff[3]);
  inverse[4]=determinant*coeff[0];
  inverse[5]=determinant*(coeff[2]*coeff[3]-coeff[0]*coeff[5]);
}

static void InvertPerspectiveCoefficients(const double *coeff,
  double *inverse)
{
  /* From "Digital Image Warping" by George Wolberg, page 53 */
  double determinant;

  determinant=1.0/(coeff[0]*coeff[4]-coeff[3]*coeff[1]);
  inverse[0]=determinant*(coeff[4]-coeff[7]*coeff[5]);
  inverse[1]=determinant*(coeff[7]*coeff[2]-coeff[1]);
  inverse[2]=determinant*(coeff[1]*coeff[5]-coeff[4]*coeff[2]);
  inverse[3]=determinant*(coeff[6]*coeff[5]-coeff[3]);
  inverse[4]=determinant*(coeff[0]-coeff[6]*coeff[2]);
  inverse[5]=determinant*(coeff[3]*coeff[2]-coeff[0]*coeff[5]);
  inverse[6]=determinant*(coeff[3]*coeff[7]-coeff[6]*coeff[4]);
  inverse[7]=determinant*(coeff[6]*coeff[1]-coeff[0]*coeff[7]);
}

static inline double MagickRound(double x)
{
  /* round the fraction to nearest integer */
  if (x >= 0.0)
    return((double) ((long) (x+0.5)));
  return((double) ((long) (x-0.5)));
}

static unsigned long poly_number_terms(double order)
{
  /* Return the number of terms for a 2d polynomial
     Order must either be an integer, or 1.5 to produce
     the 2 number_valuesal polyminal function...
        affine     1 (3)      u = c0 + c1*x + c2*y
        bilinear  1.5 (4)     u = '' + c3*x*y
        quadratic  2 (6)      u = '' + c4*x*x + c5*y*y
        cubic      3 (10)     u = '' + c6*x^3 + c7*x*x*y + c8*x*y*y + c9*y^3
        quartic    4 (15)     u = '' + c10*x^4 + ... + c14*y^4
        quintic    5 (21)     u = '' + c15*x^5 + ... + c20*y^5
     number in parenthesis minimum number of points needed.
     Anything beyond quintic, has not been implemented until
     a more automated way of determined terms is found.
   */
  if ( order < 1 || order > 5 ||
       ( order != floor(order) && (order-1.5) > MagickEpsilon) )
    return 0; /* invalid polynomial order */
  return ((unsigned long) floor((order+1)*(order+2)/2));
}

static double poly_basis_fn(long n, double x, double y)
{
  /* return the result for this polynomial term */
  switch(n) {
    case  0:  return( 1.0 ); /* constant */
    case  1:  return(  x  );
    case  2:  return(  y  ); /* affine      order = 1   terms = 3 */
    case  3:  return( x*y ); /* bilinear    order = 1.5 terms = 4 */
    case  4:  return( x*x );
    case  5:  return( y*y ); /* quadratic   order = 2   terms = 6 */
    case  6:  return( x*x*x );
    case  7:  return( x*x*y );
    case  8:  return( x*y*y );
    case  9:  return( y*y*y ); /* cubic       order = 3   terms = 10 */
    case 10:  return( x*x*x*x );
    case 11:  return( x*x*x*y );
    case 12:  return( x*x*y*y );
    case 13:  return( x*y*y*y );
    case 14:  return( y*y*y*y ); /* quartic     order = 4   terms = 15 */
    case 15:  return( x*x*x*x*x );
    case 16:  return( x*x*x*x*y );
    case 17:  return( x*x*x*y*y );
    case 18:  return( x*x*y*y*y );
    case 19:  return( x*y*y*y*y );
    case 20:  return( y*y*y*y*y ); /* quintic     order = 5   terms = 21 */
  }
  return( 0 ); /* should never happen */
}
static const char *poly_basis_str(long n)
{
  /* return the result for this polynomial term */
  switch(n) {
    case  0:  return(""); /* constant */
    case  1:  return("*ii");
    case  2:  return("*jj"); /* affiine      order = 1   terms = 3 */
    case  3:  return("*ii*jj"); /* biiliinear    order = 1.5 terms = 4 */
    case  4:  return("*ii*ii");
    case  5:  return("*jj*jj"); /* quadratiic   order = 2   terms = 6 */
    case  6:  return("*ii*ii*ii");
    case  7:  return("*ii*ii*jj");
    case  8:  return("*ii*jj*jj");
    case  9:  return("*jj*jj*jj"); /* cubiic       order = 3   terms = 10 */
    case 10:  return("*ii*ii*ii*ii");
    case 11:  return("*ii*ii*ii*jj");
    case 12:  return("*ii*ii*jj*jj");
    case 13:  return("*ii*jj*jj*jj");
    case 14:  return("*jj*jj*jj*jj"); /* quartiic     order = 4   terms = 15 */
    case 15:  return("*ii*ii*ii*ii*ii");
    case 16:  return("*ii*ii*ii*ii*jj");
    case 17:  return("*ii*ii*ii*jj*jj");
    case 18:  return("*ii*ii*jj*jj*jj");
    case 19:  return("*ii*jj*jj*jj*jj");
    case 20:  return("*jj*jj*jj*jj*jj"); /* quiintiic     order = 5   terms = 21 */
  }
  return( "UNKNOWN" ); /* should never happen */
}
static double poly_basis_dx(long n, double x, double y)
{
  /* polynomial term for x derivative */
  switch(n) {
    case  0:  return( 0.0 ); /* constant */
    case  1:  return( 1.0 );
    case  2:  return( 0.0 ); /* affine      order = 1   terms = 3 */
    case  3:  return(  y  ); /* bilinear    order = 1.5 terms = 4 */
    case  4:  return(  x  );
    case  5:  return( 0.0 ); /* quadratic   order = 2   terms = 6 */
    case  6:  return( x*x );
    case  7:  return( x*y );
    case  8:  return( y*y );
    case  9:  return( 0.0 ); /* cubic       order = 3   terms = 10 */
    case 10:  return( x*x*x );
    case 11:  return( x*x*y );
    case 12:  return( x*y*y );
    case 13:  return( y*y*y );
    case 14:  return( 0.0 ); /* quartic     order = 4   terms = 15 */
    case 15:  return( x*x*x*x );
    case 16:  return( x*x*x*y );
    case 17:  return( x*x*y*y );
    case 18:  return( x*y*y*y );
    case 19:  return( y*y*y*y );
    case 20:  return( 0.0 ); /* quintic     order = 5   terms = 21 */
  }
  return( 0.0 ); /* should never happen */
}
static double poly_basis_dy(long n, double x, double y)
{
  /* polynomial term for y derivative */
  switch(n) {
    case  0:  return( 0.0 ); /* constant */
    case  1:  return( 0.0 );
    case  2:  return( 1.0 ); /* affine      order = 1   terms = 3 */
    case  3:  return(  x  ); /* bilinear    order = 1.5 terms = 4 */
    case  4:  return( 0.0 );
    case  5:  return(  y  ); /* quadratic   order = 2   terms = 6 */
    default:  return( poly_basis_dx(n-1,x,y) ); /* weird but true */
  }
  /* NOTE: the only reason that last is not true for 'quadtratic'
     is due to the re-arrangement of terms to allow for 'bilinear'
  */
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   G e n e r a t e C o e f f i c i e n t s                                   %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GenerateCoefficients() takes user provided input arguments and generates
%  the coefficients, needed to apply the specific distortion for either
%  distorting images (generally using control points) or generating a color
%  gradient from sparsely separated color points.
%
%  The format of the GenerateCoefficients() method is:
%
%    Image *GenerateCoefficients(const Image *image, DistortImageMethod method,
%        const unsigned long number_arguments,const double *arguments,
%        unsigned long number_values, ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the image to be distorted.
%
%    o method: the method of image distortion/ sparse gradient
%
%    o number_arguments: the number of arguments given.
%
%    o arguments: the arguments for this distortion method.
%
%    o number_values: the style and format of given control points, (caller type)
%         0: 2 dimensional mapping of control points (Distort)
%            Format:  u,v,x,y  where u,v is the 'source' of the
%            the color to be plotted, for DistortImage()
%         N: Interpolation of control points with N values (usally r,g,b)
%            Format: x,y,r,g,b    mapping x,y to color values r,g,b
%            IN future, varible number of values may be given (1 to N)
%
%    o exception: Return any errors or warnings in this structure
%
%  Note that the returned array of double values must be freed by the
%  calling method using RelinquishMagickMemory().  This however may change in
%  the future to require a more 'method' specific method.
%
%  Because of this this method should not be classed as stable or used
%  outside other MagickCore library methods.
*/

static double *GenerateCoefficients(const Image *image,
  DistortImageMethod *method,const unsigned long number_arguments,
  const double *arguments,unsigned long number_values,ExceptionInfo *exception)
{
  double
    *coeff;

  register unsigned long
    i;

  unsigned long
    number_coeff, /* number of coefficients to return (array size) */
    cp_size,      /* number floating point numbers per control point */
    cp_x,cp_y,    /* the x,y indexes for control point */
    cp_values;    /* index of values for this control point */
    /* number_values   Number of values given per control point */

  if ( number_values == 0 ) {
    /* Image distortion using control points (or other distortion)
       That is generate a mapping so that   x,y->u,v   given  u,v,x,y
    */
    number_values = 2;   /* special case: two values of u,v */
    cp_values = 0;       /* with the values are BEFORE the CP x,y */
    cp_x = 2;            /* location of x,y in input control values */
    cp_y = 3;
    /* NOTE: cp_values, also used for later 'reverse map distort' tests */
  }
  else {
    cp_x = 0;            /* location of x,y in input control values */
    cp_y = 1;
    cp_values = 2;       /* and the other values are after x,y */
    /* Typically in this case the values are R,G,B color values */
  }
  cp_size = number_values+2; /* each CP defintion involves this many numbers */


  /* If not enough control point pairs are found for specific distortions
    fall back to Affine distortion (allowing 0 to 3 point pairs)
  */
  if ( number_arguments < 4*cp_size &&
       (  *method == BilinearDistortion
       || *method == PerspectiveDistortion
       ) )
    *method = AffineDistortion;

  switch (*method) {
    case AffineDistortion:
    /* also BarycentricColorInterpolate: */
      number_coeff=3*number_values;
      break;
    case BilinearDistortion:
      number_coeff=4*number_values;
      break;
    case PolynomialDistortion:
      /* number of coefficents depend on the given polynomal 'order' */
      if ( number_arguments <= 1 && (number_arguments-1)%cp_size != 0)
      {
        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
                   "InvalidArgument","%s : '%s'","Polynomial",
                   "Invalid number of args: order [CPs]...");
        return((double *) NULL);
      }
      i = poly_number_terms(arguments[0]);
      number_coeff = 2 + i*number_values;
      if ( i == 0 ) {
        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
                   "InvalidArgument","%s : '%s'","Polynomial",
                   "Invalid order, should be 1 to 5, or 1.5");
        return((double *) NULL);
      }
      if ( number_arguments < 1+i*cp_size ) {
        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
               "InvalidArgument", "%s : 'require at least %ld CPs'",
               "Polynomial", i);
        return((double *) NULL);
      }
      break;
    /* The rest are constants as they are only used for image distorts */
    case ShepardsDistortion:
    case VoronoiColorInterpolate:
      number_coeff=1;  /* may not be used, but provide some type of return */
      break;
    case ArcDistortion:
      number_coeff=5;
      break;
    case ScaleRotateTranslateDistortion:
    case AffineProjectionDistortion:
      number_coeff=6;
      break;
    case PolarDistortion:
    case DePolarDistortion:
      number_coeff=8;
      number_coeff=8;
      break;
    case PerspectiveDistortion:
    case PerspectiveProjectionDistortion:
      number_coeff=9;
      break;
    case BarrelDistortion:
    case BarrelInverseDistortion:
      number_coeff=10;
      break;
    case UndefinedDistortion:
    default:
      assert(! "Unknown Method Given"); /* just fail assertion */
  }

  /* allocate the array of coefficients needed */
  coeff = (double *) AcquireQuantumMemory(number_coeff,sizeof(*coeff));
  if (coeff == (double *) NULL) {
    (void) ThrowMagickException(exception,GetMagickModule(),
                  ResourceLimitError,"MemoryAllocationFailed",
                  "%s", "GenerateCoefficients");
    return((double *) NULL);
  }

  /* zero out coeffiecents array */
  for (i=0; i < number_coeff; i++)
    coeff[i] = 0.0;

  switch (*method)
  {
    case AffineDistortion:
    {
      /* Affine Distortion
           v =  c0*x + c1*y + c2
         for each 'value' given

         Input Arguments are sets of control points...
         For Distort Images    u,v, x,y  ...
         For Sparse Gradients  x,y, r,g,b  ...
      */
      if ( number_arguments%cp_size != 0 ||
           number_arguments < cp_size ) {
        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
               "InvalidArgument", "%s : 'require at least %ld CPs'",
               "Affine", 1L);
        coeff=(double *) RelinquishMagickMemory(coeff);
        return((double *) NULL);
      }
      /* handle special cases of not enough arguments */
      if ( number_arguments == cp_size ) {
        /* Only 1 CP Set Given */
        if ( cp_values == 0 ) {
          /* image distortion - translate the image */
          coeff[0] = 1.0;
          coeff[2] = arguments[0] - arguments[2];
          coeff[4] = 1.0;
          coeff[5] = arguments[1] - arguments[3];
        }
        else {
          /* sparse gradient - use the values directly */
          for (i=0; i<number_values; i++)
            coeff[i*3+2] = arguments[cp_values+i];
        }
      }
      else {
        /* 2 or more points (usally 3) given.
           Solve a least squares simultanious equation for coefficients.
        */
        double
          **matrix,
          **vectors,
          terms[3];

        MagickBooleanType
          status;

        /* create matrix, and a fake vectors matrix */
        matrix = AcquireMagickMatrix(3UL,3UL);
        vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
        if (matrix == (double **) NULL || vectors == (double **) NULL)
        {
          matrix  = RelinquishMagickMatrix(matrix, 3UL);
          vectors = (double **) RelinquishMagickMemory(vectors);
          coeff   = (double *) RelinquishMagickMemory(coeff);
          (void) ThrowMagickException(exception,GetMagickModule(),
                  ResourceLimitError,"MemoryAllocationFailed",
                  "%s", "DistortCoefficients");
          return((double *) NULL);
        }
        /* fake a number_values x3 vectors matrix from coefficients array */
        for (i=0; i < number_values; i++)
          vectors[i] = &(coeff[i*3]);
        /* Add given control point pairs for least squares solving */
        for (i=0; i < number_arguments; i+=cp_size) {
          terms[0] = arguments[i+cp_x];  /* x */
          terms[1] = arguments[i+cp_y];  /* y */
          terms[2] = 1;                  /* 1 */
          LeastSquaresAddTerms(matrix,vectors,terms,
                   &(arguments[i+cp_values]),3UL,number_values);
        }
        if ( number_arguments == 2*cp_size ) {
          /* Only two pairs were given, but we need 3 to solve the affine.
             Fake extra coordinates by rotating p1 around p0 by 90 degrees.
               x2 = x0 - (y1-y0)   y2 = y0 + (x1-x0)
           */
          terms[0] = arguments[cp_x]
                   - ( arguments[cp_size+cp_y] - arguments[cp_y] ); /* x2 */
          terms[1] = arguments[cp_y] +
                   + ( arguments[cp_size+cp_x] - arguments[cp_x] ); /* y2 */
          terms[2] = 1;                                             /* 1 */
          if ( cp_values == 0 ) {
            /* Image Distortion - rotate the u,v coordients too */
            double
              uv2[2];
            uv2[0] = arguments[0] - arguments[5] + arguments[1];   /* u2 */
            uv2[1] = arguments[1] + arguments[4] - arguments[0];   /* v2 */
            LeastSquaresAddTerms(matrix,vectors,terms,uv2,3UL,2UL);
          }
          else {
            /* Sparse Gradient - use values of p0 for linear gradient */
            LeastSquaresAddTerms(matrix,vectors,terms,
                  &(arguments[cp_values]),3UL,number_values);
          }
        }
        /* Solve for LeastSquares Coefficients */
        status=GaussJordanElimination(matrix,vectors,3UL,number_values);
        matrix = RelinquishMagickMatrix(matrix, 3UL);
        vectors = (double **) RelinquishMagickMemory(vectors);
        if ( status == MagickFalse ) {
          coeff = (double *) RelinquishMagickMemory(coeff);
          (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
                      "InvalidArgument","%s : '%s'","Affine",
                      "Unsolvable Matrix");
          return((double *) NULL);
        }
      }
      return(coeff);
    }
    case AffineProjectionDistortion:
    {
      /*
        Arguments: Affine Matrix (forward mapping)
        Arguments  sx, rx, ry, sy, tx, ty
        Where      u = sx*x + ry*y + tx
                   v = rx*x + sy*y + ty

        Returns coefficients (in there inverse form) ordered as...
             sx ry tx  rx sy ty

        AffineProjection Distortion Notes...
           + Will only work with a 2 number_values for Image Distortion
           + Can not be used for generating a sparse gradient (interpolation)
      */
      double inverse[8];
      if (number_arguments != 6) {
        coeff = (double *) RelinquishMagickMemory(coeff);
        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
                     "InvalidArgument","%s : '%s'","AffineProjection",
                     "Needs 6 coeff values");
        return((double *) NULL);
      }
      /* FUTURE: trap test for sx*sy-rx*ry == 0 (determinate = 0, no inverse) */
      for(i=0; i<6UL; i++ )
        inverse[i] = arguments[i];
      AffineArgsToCoefficients(inverse); /* map into coefficents */
      InvertAffineCoefficients(inverse, coeff); /* invert */
      *method = AffineDistortion;

      return(coeff);
    }
    case ScaleRotateTranslateDistortion:
    {
      /* Scale, Rotate and Translate Distortion
         An alturnative Affine Distortion
         Argument options, by number of arguments given:
           7: x,y, sx,sy, a, nx,ny
           6: x,y,   s,   a, nx,ny
           5: x,y, sx,sy, a
           4: x,y,   s,   a
           3: x,y,        a
           2:        s,   a
           1:             a
         Where actions are (in order of application)
            x,y     'center' of transforms     (default = image center)
            sx,sy   scale image by this amount (default = 1)
            a       angle of rotation          (argument required)
            nx,ny   move 'center' here         (default = no movement)
         And convert to affine mapping coefficients

         ScaleRotateTranslate Distortion Notes...
           + Does not use a set of CPs in any normal way
           + Will only work with a 2 number_valuesal Image Distortion
           + Can not be used for generating a sparse gradient (interpolation)
      */
      double
        cosine, sine,
        x,y,sx,sy,a,nx,ny;

      /* set default center, and default scale */
      x = nx = (double)(image->columns-1.0)/2.0 + image->page.x;
      y = ny = (double)(image->rows-1.0)/2.0    + image->page.y;
      sx = sy = 1.0;
      switch ( number_arguments ) {
      case 0:
        coeff = (double *) RelinquishMagickMemory(coeff);
        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
                     "InvalidArgument","%s : '%s'", "ScaleTranslateRotate",
                     "Needs at least 1 argument");
        return((double *) NULL);
      case 1:
        a = arguments[0];
        break;
      case 2:
        sx = sy = arguments[0];
        a = arguments[1];
        break;
      default:
        x = nx = arguments[0];
        y = ny = arguments[1];
        switch ( number_arguments ) {
        case 3:
          a = arguments[2];
          break;
        case 4:
          sx = sy = arguments[2];
          a = arguments[3];
          break;
        case 5:
          sx = arguments[2];
          sy = arguments[3];
          a = arguments[4];
          break;
        case 6:
          sx = sy = arguments[2];
          a = arguments[3];
          nx = arguments[4];
          ny = arguments[5];
          break;
        case 7:
          sx = arguments[2];
          sy = arguments[3];
          a = arguments[4];
          nx = arguments[5];
          ny = arguments[6];
          break;
        default:
          coeff = (double *) RelinquishMagickMemory(coeff);
          (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
                     "InvalidArgument","%s : '%s'", "ScaleTranslateRotate",
                     "Too Many Arguments (7 or less)");
          return((double *) NULL);
        }
        break;
      }
      /* Trap if sx or sy == 0 -- image is scaled out of existance! */
      if ( fabs(sx) < MagickEpsilon || fabs(sy) < MagickEpsilon ) {
        coeff = (double *) RelinquishMagickMemory(coeff);
        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
                     "InvalidArgument","%s : '%s'", "ScaleTranslateRotate",
                     "Zero Scale Given");
        return((double *) NULL);
      }
      /* Save the given arguments as an affine distortion */
      a=DegreesToRadians(a); cosine=cos(a); sine=sin(a);

      *method = AffineDistortion;
      coeff[0]=cosine/sx;
      coeff[1]=sine/sx;
      coeff[2]=x-nx*coeff[0]-ny*coeff[1];
      coeff[3]=(-sine)/sy;
      coeff[4]=cosine/sy;
      coeff[5]=y-nx*coeff[3]-ny*coeff[4];
      return(coeff);
    }
    case PerspectiveDistortion:
    { /*
         Perspective Distortion (a ratio of affine distortions)

                p(x,y)    c0*x + c1*y + c2
            u = ------ = ------------------
                r(x,y)    c6*x + c7*y + 1

                q(x,y)    c3*x + c4*y + c5
            v = ------ = ------------------
                r(x,y)    c6*x + c7*y + 1

           c8 = Sign of 'r', or the denominator affine, for the actual image.
                This determines what part of the distorted image is 'ground'
                side of the horizon, the other part is 'sky' or invalid.
                Valid values are  +1.0  or  -1.0  only.

         Input Arguments are sets of control points...
         For Distort Images    u,v, x,y  ...
         For Sparse Gradients  x,y, r,g,b  ...

         Perspective Distortion Notes...
           + Can be thought of as ratio of  3 affine transformations
           + Not separatable: r() or c6 and c7 are used by both equations
           + All 8 coefficients must be determined simultaniously
           + Will only work with a 2 number_valuesal Image Distortion
           + Can not be used for generating a sparse gradient (interpolation)
           + It is not linear, but is simple to generate an inverse
           + All lines within an image remain lines.
           + but distances between points may vary.
      */
      double
        **matrix,
        *vectors[1],
        terms[8];

      unsigned long
        cp_u = cp_values,
        cp_v = cp_values+1;

      MagickBooleanType
        status;

      /* fake 1x8 vectors matrix directly using the coefficients array */
      vectors[0] = &(coeff[0]);
      /* 8x8 least-squares matrix (zeroed) */
      matrix = AcquireMagickMatrix(8UL,8UL);
      if (matrix == (double **) NULL) {
        (void) ThrowMagickException(exception,GetMagickModule(),
                  ResourceLimitError,"MemoryAllocationFailed",
                  "%s", "DistortCoefficients");
        return((double *) NULL);
      }
      /* Add control points for least squares solving */
      for (i=0; i < number_arguments; i+=4) {
        terms[0]=arguments[i+cp_x];            /*   c0*x   */
        terms[1]=arguments[i+cp_y];            /*   c1*y   */
        terms[2]=1.0;                          /*   c2*1   */
        terms[3]=0.0;
        terms[4]=0.0;
        terms[5]=0.0;
        terms[6]=-terms[0]*arguments[i+cp_u];  /* 1/(c6*x) */
        terms[7]=-terms[1]*arguments[i+cp_u];  /* 1/(c7*y) */
        LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_u]),
            8UL,1UL);

        terms[0]=0.0;
        terms[1]=0.0;
        terms[2]=0.0;
        terms[3]=arguments[i+cp_x];           /*   c3*x   */
        terms[4]=arguments[i+cp_y];           /*   c4*y   */
        terms[5]=1.0;                         /*   c5*1   */
        terms[6]=-terms[3]*arguments[i+cp_v]; /* 1/(c6*x) */
        terms[7]=-terms[4]*arguments[i+cp_v]; /* 1/(c7*y) */
        LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+cp_v]),
            8UL,1UL);
      }
      /* Solve for LeastSquares Coefficients */
      status=GaussJordanElimination(matrix,vectors,8UL,1UL);
      matrix = RelinquishMagickMatrix(matrix, 8UL);
      if ( status == MagickFalse ) {
        coeff = (double *) RelinquishMagickMemory(coeff);
        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
                    "InvalidArgument","%s : '%s'","Perspective",
                    "Unsolvable Matrix");
        return((double *) NULL);
      }
      /*
        Calculate 9'th coefficient! The ground-sky determination.
        What is sign of the 'ground' in r() denominator affine function?
        Just use any valid image coordinate (first control point) in
        destination for determination of what part of view is 'ground'.
      */
      coeff[8] = coeff[6]*arguments[cp_x]
                      + coeff[7]*arguments[cp_y] + 1.0;
      coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;

      return(coeff);
    }
    case PerspectiveProjectionDistortion:
    {
      /*
        Arguments: Perspective Coefficents (forward mapping)
      */
      if (number_arguments != 8) {
        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
                     "InvalidArgument","%s : '%s'","PerspectiveProjection",
                     "Needs 8 coefficient values");
        return((double *) NULL);
      }
      /* FUTURE: trap test  c0*c4-c3*c1 == 0  (determinate = 0, no inverse) */
      InvertPerspectiveCoefficients(arguments, coeff);
      /*
        Calculate 9'th coefficient! The ground-sky determination.
        What is sign of the 'ground' in r() denominator affine function?
        Just use any valid image cocodinate in destination for determination.
        For a forward mapped perspective the images 0,0 coord will map to
        c2,c5 in the distorted image, so set the sign of denominator of that.
      */
      coeff[8] = coeff[6]*arguments[2]
                           + coeff[7]*arguments[5] + 1.0;
      coeff[8] = (coeff[8] < 0.0) ? -1.0 : +1.0;
      *method = PerspectiveDistortion;

      return(coeff);
    }
    case BilinearDistortion:
    {
      /* Bilinear Distortion
            v = c0*x + c1*y + c2*x*y + c3;
         for each 'value' given

         Input Arguments are sets of control points...
         For Distort Images    u,v, x,y  ...
         For Sparse Gradients  x,y, r,g,b  ...

         Bilinear Distortion Notes...
           + This a reversed mapped distortion
           + Can map source orthogonal rectangles to non-planar quadraterials
           + may be used for generating forward mapped 'grid' interpolation
      */
      double
        **matrix,
        **vectors,
        terms[4];

      MagickBooleanType
        status;

      /* create matrix, and a fake vectors matrix */
      matrix = AcquireMagickMatrix(4UL,4UL);
      vectors = (double **) AcquireQuantumMemory(number_values,sizeof(*vectors));
      if (matrix == (double **) NULL || vectors == (double **) NULL)
      {
        matrix  = RelinquishMagickMatrix(matrix, 4UL);
        vectors = (double **) RelinquishMagickMemory(vectors);
        coeff   = (double *) RelinquishMagickMemory(coeff);
        (void) ThrowMagickException(exception,GetMagickModule(),
                ResourceLimitError,"MemoryAllocationFailed",
                "%s", "DistortCoefficients");
        return((double *) NULL);
      }
      /* fake a number_values x3 vectors matrix from coefficients array */
      for (i=0; i < number_values; i++)
        vectors[i] = &(coeff[i*4]);
      /* Add given control point pairs for least squares solving */
      for (i=0; i < number_arguments; i+=cp_size) {
        terms[0] = arguments[i+cp_x];   /*  x  */
        terms[1] = arguments[i+cp_y];   /*  y  */
        terms[2] = terms[0]*terms[1];   /* x*y */
        terms[3] = 1;                   /*