Changeset 11612 for ImageMagick/trunk

Timestamp:

07/12/08 01:01:08 (6 weeks ago)

Author:

anthony

Message:

Separate coefficient calculation from DistortImages?() (future development)

Perspective and Bilinear Distortions fallback to Affine
when less than 4 point pairs are provided.

Affine Distortion now handles 2 point pairs (SRT of a single line)
and a simple translation if only a single point pair is provided.

Location:

ImageMagick/trunk

Files:

• ChangeLog (modified) (1 diff)
• magick/distort.c (modified) (32 diffs)

ImageMagick/trunk/ChangeLog

@@ +1 @@
+2008-07-12  6.4.2-2 Anthony <anthony@griffith...>
+  * Separate coefficient calculation from DistortImages (future development)
+  * Perspective and Bilinear Distortions fallback to Affine
+    when less than 4 point pairs are provided.
+  * Affine Distortion now handles 2 point pairs (SRT of a single line)
+    and a simple translation if only a single point pair is provided.
+
 2008-06-29  6.4.2-2 Cristy  <quetzlzacatenango@image...>
   * Add log2 and round as -fx operators.

ImageMagick/trunk/magick/distort.c

@@ -86 +86 @@
 %  account in the mapping.
 %
+%
+%  The format of the DistortImage() method is:
+%
+%      Image *DistortImage(const Image *image,const DistortImageMethod method,
+%        const unsigned long number_arguments,const double *arguments,
+%        MagickBooleanType bestfit, ExceptionInfo *exception)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image to be distorted.
+%
+%    o method: the method of image distortion.
+%
+%    o number_arguments: the number of arguments given.
+%
+%    o arguments: the arguments for this distortion method.
+%
+%    o bestfit: Attempt to 'bestfit' the size of the resulting image.
+%         This also forces the resulting image to be a 'layered' virtual
+%         canvas image.  Can be overridden using 'distort:viewport' setting.
+%
+%    o exception: Return any errors or warnings in this structure
+%
+% Specific Distortion notes...
+%
 %  ArcDistortion will always ignore source image offset, and always 'bestfit'
 %  the destination image with the top left corner offset relative to the polar
@@ -93 +118 @@
 %  of image distortion.
 %
-%  The format of the DistortImage() method is:
-%
-%      Image *DistortImage(const Image *image,const DistortImageMethod method,
-%        const unsigned long number_arguments,const double *arguments,
-%        MagickBooleanType bestfit, ExceptionInfo *exception)
-%
-%  A description of each parameter follows:
-%
-%    o image: the image to be distorted.
-%
-%    o method: the method of image distortion.
-%
-%    o number_arguments: the number of arguments given.
-%
-%    o arguments: the arguments for this distortion method.
-%
-%    o bestfit: Attempt to resize destination to fit distorted source.
-%
-%    o exception: Return any errors or warnings in this structure
+%  Affine, Perspective, and Bilinear, will do least squares fitting of the
+%  distrotion when more than the minimum number of control point pairs are
+%  provided.
+%
+%  Perspective, and Bilinear, will fall back to a Affine distortion when less
+%  that 4 control point pairs are provided. While Affine distortions will let
+%  you use any number of control point pairs, that is Zero pairs is a No-Op
+%  (viewport only) distrotion, one pair is a translation and two pairs of
+%  control points will do a scale-rotate-translate, without any shearing.
 %
 */
@@ -192 +207 @@
 
 
-MagickExport Image *DistortImage(Image *image,const DistortImageMethod method,
-  const unsigned long number_arguments,const double *arguments,
-  MagickBooleanType bestfit,ExceptionInfo *exception)
+static double *DistortCoefficents(Image *image, DistortImageMethod *method,
+     const unsigned long number_arguments, const double *arguments,
+     MagickBooleanType bestfit, ExceptionInfo *exception)
 {
-#define DistortImageTag  "Distort/Image"
-
-  const char
-    *artifact;
+  /* Determine correct coefficients for the given input arguments.
+     The given distortion method may be modified, to use a simplier
+     distortion method to better handle the requested arguments.
+
+     The returned array of coefficients, will need to be freed by the caller.
+     by calling    RelinquishMagickMemory()
+
+     If NULL is returned an exception has occured.
+
+     WARNING: this routine will change, to allow for generating coefficients
+     for a FUTURE 2d gradient generation.
+  */
 
   double
-    coefficients[9],
-    validity;
-
-  Image
-    *distort_image;
-
-  long
-    j,
-    y;
+    *coefficients;
+
+  register int
+    i;
 
   MagickBooleanType
     status;
 
-  PointInfo
-    point;          /* point to sample (center of filtered resample of area) */
-
-  RectangleInfo
-    geometry;
-
-  register long
-    i,
-    x;
-
-  ResampleFilter
-    **resample_filter;
-
-  ViewInfo
-    **distort_view;
-
-  assert(image != (Image *) NULL);
-  assert(image->signature == MagickSignature);
-  if (image->debug != MagickFalse)
-    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  assert(exception != (ExceptionInfo *) NULL);
-  assert(exception->signature == MagickSignature);
-  (void) ResetMagickMemory(coefficients,0,sizeof(coefficients));
-
-  /*
-    Convert input arguments into mapping coefficents for distortion
+  /* currently the array is always 9 doubles, that may change in the future */
+  coefficients = AcquireQuantumMemory(9,sizeof(*coefficients));
+  if (coefficients == (double *) NULL) {
+    (void) ThrowMagickException(exception,GetMagickModule(),
+                  ResourceLimitError,"MemoryAllocationFailed",
+                  "%s", "DistortCoefficients");
+    return((double *) NULL);
+  }
+  /* zero out coeffiecents array */
+  for (i=0; i < 9L; i++)
+    coefficients[i] = 0.0;
+
+  /* If not enough control point pairs are found for specific distortions
+    fall back to Affine distortion (allowing 0 to 3 point pairs)
   */
-  switch (method)
+  if ( number_arguments < 16 &&
+       (  *method == BilinearDistortion
+       || *method == PerspectiveDistortion
+       ) )
+    *method = AffineDistortion;
+
+  switch (*method)
   {
     case AffineDistortion:
     {
-      double
-        **matrix,
-        **vectors,
-        terms[3];
-
-      /* Affine Distortion
-            u =   c0*x + c2*y + c4
-            v =   c1*x + c3*y + c5
-
-         Input Arguments are pairs of distorted coodinates (minimum 3 sets)
-         Which will be used to generate the coefficients of the above.
-            u1,v1, x1,y1,  u2,v2, x2,y2,  u3,v3, x3,y3, ...
+      /* Input Arguments are any number of pairs of distorted coodinates
+         Which will be used to generate coefficients for Distortion.
+            u1,v1, x1,y1,  u2,v2, x2,y2,  u3,v3, x3,y3,  u4,v4, x4,y4 ...
       */
-      if (number_arguments < 12) {
-        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-                     "InvalidArgument","%s : '%s'","distort Affine",
-                     "Minimum of 3 pairs of coords needed (12 numbers)");
-        return((Image *) NULL);
-      }
       if (number_arguments%4 != 0) {
         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-                     "InvalidArgument","%s : '%s'","distort Affine",
+                     "InvalidArgument","%s : '%s'","distort Map",
                      "Requires pairs of coords (sX,sY,dX,dY)");
-        return((Image *) NULL);
-      }
-      matrix = AcquireMagickMatrix(3UL,3UL);
-      vectors = AcquireMagickMatrix(2UL,3UL);
-      if (matrix == (double **) NULL || vectors == (double **) NULL)
-      {
+        RelinquishMagickMemory(coefficients);
+        return((double *) NULL);
+      }
+      /* handle special cases of not enough arguments */
+      if ( number_arguments == 0 ) {
+        /* null distortion -- use No-Op Affine Distortion
+            coefficients: sx, rx, ry, sy, tx, ty
+        */
+        *method = AffineDistortion;
+        coefficients[0] = 1.0;
+        coefficients[3] = 1.0;
+      }
+      else if ( number_arguments == 4 ) {
+        /* if only 1 pair of points - Affine Translation only
+            coefficients: sx, rx, ry, sy, tx, ty
+        */
+        *method = AffineDistortion;
+        coefficients[0] = 1.0;
+        coefficients[3] = 1.0;
+        coefficients[4] = arguments[0] - arguments[2];
+        coefficients[5] = arguments[1] - arguments[3];
+      }
+      else {
+        /*
+           Affine Distortion
+              u =  c0*x + c2*y + c4
+              v =  c1*x + c3*y + c5
+        */
+        double
+          **matrix,
+          **vectors,
+          terms[3];
+
+        *method = AffineDistortion;
+        matrix = AcquireMagickMatrix(3UL,3UL);
+        vectors = AcquireMagickMatrix(2UL,3UL);
+        if (matrix == (double **) NULL || vectors == (double **) NULL)
+        {
+          matrix = RelinquishMagickMatrix(matrix, 3UL);
+          vectors = RelinquishMagickMatrix(vectors, 2UL);
+          coefficients = RelinquishMagickMemory(coefficients);
+          (void) ThrowMagickException(exception,GetMagickModule(),
+                  ResourceLimitError,"MemoryAllocationFailed",
+                  "%s", "DistortCoefficients");
+          return((double *) NULL);
+        }
+        /* Add given control point pairs for least squares solving */
+        for (i=0; i < (long) number_arguments; i+=4) {
+          terms[0] = arguments[i+2];  /* x */
+          terms[1] = arguments[i+3];  /* y */
+          terms[2] = 1;               /* 1 */
+          LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i]),3UL,2UL);
+        }
+        if ( number_arguments == 8 ) {
+          /* Only two pairs were given, but we need 3 to solve the affine.
+             Fake a 3rd control point pair by p1 around p0 by 90 degrees
+             That is given    u0,v0 x0,y0   u1,v1 x1,y1
+             Then   u2 = u0 - (v1-v0)   v2 = v0 + (u1-u0)
+                    x2 = x0 - (y1-y0)   y2 = y0 + (x1-x0)
+          */
+          double
+            uv2[2];
+
+          uv2[0] = arguments[0] - arguments[5] + arguments[1];   /* u2 */
+          uv2[1] = arguments[1] + arguments[4] + arguments[0];   /* v2 */
+          terms[0] = arguments[2] - arguments[7] + arguments[3]; /* x2 */
+          terms[1] = arguments[3] + arguments[6] + arguments[2]; /* y2 */
+          terms[2] = 1;                                          /* 1 */
+          LeastSquaresAddTerms(matrix,vectors,terms,uv2,3UL,2UL);
+        }
+        status=GaussJordanElimination(matrix,vectors,3UL,2UL);
+        if ( status == MagickTrue ) {
+          /* Affine is a bit weird in its ordering of coefficients.
+            This ordering is for backward compatibility with the older
+            affine functions of IM.
+
+            coefficients order: sx, rx, ry, sy, tx, ty
+            Where    u = sx*x + ry*y + tx
+                     v = rx*x + sy*y + ty
+          */
+          coefficients[0] = vectors[0][0]; /* sx */
+          coefficients[1] = vectors[1][0]; /* rx */
+          coefficients[2] = vectors[0][1]; /* ry */
+          coefficients[3] = vectors[1][1]; /* sy */
+          coefficients[4] = vectors[0][2]; /* tx */
+          coefficients[5] = vectors[1][2]; /* ty */
+        }
         matrix = RelinquishMagickMatrix(matrix, 3UL);
         vectors = RelinquishMagickMatrix(vectors, 2UL);
-        ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
-        return((Image *) NULL);
-      }
-      /* Add control points for least squares solving */
-      for (i=0; i < (long) number_arguments; i+=4) {
-        terms[0] = arguments[i+2];  /* x */
-        terms[1] = arguments[i+3];  /* y */
-        terms[2] = 1;               /* 1 */
-        LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i]),3UL,2UL);
-      }
-      status=GaussJordanElimination(matrix,vectors,3UL,2UL);
-      matrix = RelinquishMagickMatrix(matrix, 3UL);
-      if ( status == MagickTrue ) {
-        /* Affine has a bit of a weird ordering of coefficents.
-           This ordering is for backward compatibility with the older
-           affine functions of IM. (See AffineProjection next)
-        */
-        coefficients[0] = vectors[0][0]; /* u = ... */
-        coefficients[2] = vectors[0][1];
-        coefficients[4] = vectors[0][2];
-        coefficients[1] = vectors[1][0]; /* v = ... */
-        coefficients[3] = vectors[1][1];
-        coefficients[5] = vectors[1][2];
-      }
-      vectors = RelinquishMagickMatrix(vectors, 2UL);
-      if ( status == MagickFalse ) {
-        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-                     "InvalidArgument","%s : '%s'","distort Affine",
-                     "Unsolvable Matrix");
-        return((Image *) NULL);
-      }
-      break;
+        if ( status == MagickFalse ) {
+          coefficients = RelinquishMagickMemory(coefficients);
+          (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
+                      "InvalidArgument","%s : '%s'","distort Affine",
+                      "Unsolvable Matrix");
+          return((double *) NULL);
+        }
+      }
+      if ( GetImageArtifact(image,"distort:verbose") != (const char *) NULL ) {
+        double *inverse = AcquireQuantumMemory(8,sizeof(coefficients[i]));
+        if (inverse == (double *) NULL) {
+          coefficients = RelinquishMagickMemory(coefficients);
+          (void) ThrowMagickException(exception,GetMagickModule(),
+                  ResourceLimitError,"MemoryAllocationFailed",
+                  "%s", "DistortCoefficients");
+          return((double *) NULL);
+        }
+        InvertAffineCoefficients(coefficients, inverse);
+        fprintf(stderr, "Affine ReverseMap\n");
+        for (i=0; i<6; i++)
+          fprintf(stderr, "  %lf\n", coefficients[i]);
+        fprintf(stderr, "Affine Projection\n");
+        for (i=0; i<6; i++)
+          fprintf(stderr, "  %lf\n", inverse[i]);
+        inverse = RelinquishMagickMemory(inverse);
+      }
+      return(coefficients);
     }
     case AffineProjectionDistortion:
@@ -314 +387 @@
       /*
         Arguments: Affine Matrix (forward mapping)
-           sx, rx, ry, sy, tx, ty
+        Arguments  sx, rx, ry, sy, tx, ty
+        Where      u = sx*x + ry*y + tx
+                   v = rx*x + sy*y + ty
       */
       if (number_arguments != 6) {
+        coefficients = RelinquishMagickMemory(coefficients);
         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
                      "InvalidArgument","%s : '%s'","distort AffineProjection",
                      "Needs 6 numbers");
-        return((Image *) NULL);
+        return((double *) NULL);
       }
       InvertAffineCoefficients(arguments, coefficients);
-      break;
-    }
-    case BilinearDistortion:
-    {
-      double
-        **matrix,
-        *vectors[2],
-        terms[4];
-
-      /* Bilinear Distortion (reversed)
-
-            u = c0*x + c1*y + c2*x*y + c3;
-
-            v = c4*x + c5*y + c6*x*y + c7;
-
-         Input Arguments are pairs of distorted coodinates (minimum 4 pairs)
-         Which will be used to generate the coefficients of the above.
-            u1,v1, x1,y1,  u2,v2, x2,y2,  u3,v3, x3,y3,  u4,v4, x4,y4 ...
-      */
-      if (number_arguments < 16) {
-        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-                     "InvalidArgument","%s : '%s'","distort Bilinear",
-                     "Minimum of 4 pairs of coords needed (16 numbers)");
-        return((Image *) NULL);
-      }
-      if (number_arguments%4 != 0) {
-        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-                     "InvalidArgument","%s : '%s'","distort Bilinear",
-                     "Requires pairs of coords (sX,sY,dX,dY)");
-        return((Image *) NULL);
-      }
-      matrix = AcquireMagickMatrix(4UL,4UL);
-      if (matrix == (double **) NULL) {
-        matrix = RelinquishMagickMatrix(matrix, 4UL);
-        ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
-        return((Image *) NULL);
-      }
-      /* fake a 2x3 vectors matrix from coefficients array */
-      for (i=0; i < 8L; i++)
-        coefficients[i] = 0.0;
-      vectors[0] = &(coefficients[0]);
-      vectors[1] = &(coefficients[4]);
-      /* Add control points for least squares solving */
-      for (i=0; i < (long) number_arguments; i+=4) {
-        terms[0] = arguments[i+2];     /*  x  */
-        terms[1] = arguments[i+3];     /*  y  */
-        terms[2] = terms[0]*terms[1];  /* x*y */
-        terms[3] = 1;                  /*  1  */
-        LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i]),4UL,2UL);
-      }
-      status=GaussJordanElimination(matrix,vectors,4UL,2UL);
-      matrix = RelinquishMagickMatrix(matrix, 4UL);
-      if ( status == MagickFalse ) {
-        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-                     "InvalidArgument","%s : '%s'","distort Bilinear",
-                     "Unsolvable Matrix");
-        return((Image *) NULL);
-      }
-      break;
+      if ( GetImageArtifact(image,"distort:verbose") != (const char *) NULL ) {
+        fprintf(stderr, "Affine ReverseMap\n");
+        for (i=0; i<6; i++)
+          fprintf(stderr, "  %lf\n", coefficients[i]);
+        fprintf(stderr, "AffineProjection\n");
+        for (i=0; i<6; i++)
+          fprintf(stderr, "  %lf\n", arguments[i]);
+      }
+      return(coefficients);
     }
     case PerspectiveDistortion:
-    {
+    { /*
+         Perspective Distortion (a ratio of affine distortions)
+
+                p()     c0*x + c1*y + c2
+            u = ----- = ------------------
+                r()     c6*x + c7*y + 1
+
+                q()     c3*x + c4*y + c5
+            v = ----- = ------------------
+                r()     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.
+
+         WARNING: Perspective Equations are linked (not separatable) as c6
+         and c7 are shared by both equations.  All 8 coefficients need to be
+         determined simultaneously.
+      */
       double
         **matrix,
@@ -390 +435 @@
         terms[8];
 
-      /* Perspective Distortion (a ratio of affine distortions)
-
-                 p()     c0*x + c1*y + c2
-            u = ----- = ------------------
-                 r()     c6*x + c7*y + 1
-
-                 q()     c3*x + c4*y + c5
-            v = ----- = ------------------
-                 r()     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.
-
-         WARNING: Equations are linked, and not separable.
-
-         Input Arguments are pairs of distorted coodinates (minimum 4 pairs)
-         Which will be used to generate the coefficients of the above.
-            u1,v1, x1,y1,  u2,v2, x2,y2,  u3,v3, x3,y3,  u4,v4, x4,y4 ...
-      */
-      if (number_arguments < 16) {
-        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-                     "InvalidArgument","%s : '%s'","distort Bilinear",
-                     "Minimum of 4 pairs of coords needed (16 numbers)");
-        return((Image *) NULL);
-      }
-      if (number_arguments%4 != 0) {
-        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-                     "InvalidArgument","%s : '%s'","distort Bilinear",
-                     "Requires pairs of coords (sX,sY,dX,dY)");
-        return((Image *) NULL);
-      }
+      *method = PerspectiveDistortion;
+      /* fake 1x8 vectors matrix directly from coefficients array */
+      vectors[0] = &(coefficients[0]);
+      /* 8x8 least-squares matrix (zeroed) */
       matrix = AcquireMagickMatrix(8UL,8UL);
       if (matrix == (double **) NULL) {
-        matrix = RelinquishMagickMatrix(matrix, 8UL);
-        ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
-        return((Image *) NULL);
-      }
-      /* fake a 1x8 vectors matrix from coefficients array */
-      for (i=0; i < 8L; i++)
-        coefficients[i] = 0.0;
-      vectors[0] = &(coefficients[0]);
+        (void) ThrowMagickException(exception,GetMagickModule(),
+                  ResourceLimitError,"MemoryAllocationFailed",
+                  "%s", "DistortCoefficients");
+        return((double *) NULL);
+      }
       /* Add control points for least squares solving */
       for (i=0; i < (long) number_arguments; i+=4) {
@@ -443 +457 @@
         terms[7]=-terms[1]*arguments[i];  /* 1/(c7*y) */
         LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i]),
-             8UL,1UL);
+            8UL,1UL);
 
         terms[0]=0.0;
@@ -454 +468 @@
         terms[7]=-terms[4]*arguments[i+1];  /* 1/(c7*y) */
         LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i+1]),
-             8UL,1UL);
+            8UL,1UL);
       }
       status=GaussJordanElimination(matrix,vectors,8UL,1UL);
       matrix = RelinquishMagickMatrix(matrix, 8UL);
       if ( status == MagickFalse ) {
+        coefficients = RelinquishMagickMemory(coefficients);
         (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
-                     "InvalidArgument","%s : '%s'","distort Prespective",
-                     "Unsolvable Matrix");
-        return((Image *) NULL);
+                    "InvalidArgument","%s : '%s'","distort Prespective",
+                    "Unsolvable Matrix");
+        return((double *) NULL);
       }
       /*
@@ -471 +486 @@
       */
       coefficients[8] = coefficients[6]*arguments[2]
-                           + coefficients[7]*arguments[3] + 1.0;
+                          + coefficients[7]*arguments[3] + 1.0;
       coefficients[8] = (coefficients[8] < 0.0) ? -1.0 : +1.0;
-      break;
+
+      if ( GetImageArtifact(image,"distort:verbose") != (const char *) NULL ) {
+        double *inverse = AcquireQuantumMemory(8,sizeof(coefficients[i]));
+        if (inverse == (double *) NULL) {
+          coefficients = RelinquishMagickMemory(coefficients);
+          (void) ThrowMagickException(exception,GetMagickModule(),
+                  ResourceLimitError,"MemoryAllocationFailed",
+                  "%s", "DistortCoefficients");
+          return((double *) NULL);
+        }
+        InvertPerspectiveCoefficients(coefficients, inverse);
+        fprintf(stderr, "Perspective ReverseMap\n");
+        for (i=0; i<8; i++)
+          fprintf(stderr, "  %lf\n", coefficients[i]);
+        fprintf(stderr, "PerspectiveProjection\n");
+        for (i=0; i<8; i++)
+          fprintf(stderr, "  %lf\n", inverse[i]);
+        inverse = RelinquishMagickMemory(inverse);
+      }
+      return(coefficients);
     }
     case PerspectiveProjectionDistortion:
@@ -484 +518 @@
                      "InvalidArgument","%s : '%s'",
                      "distort PerspectiveProjection", "Needs 8 numbers");
-        return((Image *) NULL);
+        return((double *) NULL);
       }
       InvertPerspectiveCoefficients(arguments, coefficients);
+      if ( GetImageArtifact(image,"distort:verbose") != (const char *) NULL ) {
+        fprintf(stderr, "Perspective ReverseMap\n");
+        for (i=0; i<8; i++)
+          fprintf(stderr, "  %lf\n", coefficients[i]);
+        fprintf(stderr, "PerspectiveProjection\n");
+        for (i=0; i<8; i++)
+          fprintf(stderr, "  %lf\n", arguments[i]);
+      }
       /*
         Calculate 9'th coefficient! The ground-sky determination.
@@ -492 +534 @@
         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 get the sign of denominator of that.
+        c2,c5 in the distorted image, so set the sign of denominator of that.
       */
       coefficients[8] = coefficients[6]*arguments[2]
                            + coefficients[7]*arguments[5] + 1.0;
       coefficients[8] = (coefficients[8] < 0.0) ? -1.0 : +1.0;
-      break;
+      return(coefficients);
+    }
+    case BilinearDistortion:
+    {
+      double
+        **matrix,
+        *vectors[2],
+        terms[4];
+
+      /* Bilinear Distortion (currently reversed -- this needs to be fixed)
+            u = c0*x + c1*y + c2*x*y + c3;
+            v = c4*x + c5*y + c6*x*y + c7;
+
+         Input Arguments are pairs of distorted coodinates (minimum 4 pairs)
+         Which will be used to generate the coefficients of the above.
+            u1,v1, x1,y1,  u2,v2, x2,y2,  u3,v3, x3,y3,  u4,v4, x4,y4 ...
+      */
+      if (number_arguments%4 != 0) {
+        coefficients = RelinquishMagickMemory(coefficients);
+        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
+                     "InvalidArgument","%s : '%s'","distort Bilinear",
+                     "Requires pairs of coords (sX,sY,dX,dY)");
+        return((double *) NULL);
+      }
+      matrix = AcquireMagickMatrix(4UL,4UL);
+      if (matrix == (double **) NULL) {
+        coefficients = RelinquishMagickMemory(coefficients);
+        (void) ThrowMagickException(exception,GetMagickModule(),
+                  ResourceLimitError,"MemoryAllocationFailed",
+                  "%s", "DistortCoefficients");
+        return((double *) NULL);
+      }
+      /* fake a 2x3 vectors matrix from coefficients array */
+      vectors[0] = &(coefficients[0]);
+      vectors[1] = &(coefficients[4]);
+      /* Add control points for least squares solving */
+      for (i=0; i < (long) number_arguments; i+=4) {
+        terms[0] = arguments[i+2];     /*  x  */
+        terms[1] = arguments[i+3];     /*  y  */
+        terms[2] = terms[0]*terms[1];  /* x*y */
+        terms[3] = 1;                  /*  1  */
+        LeastSquaresAddTerms(matrix,vectors,terms,&(arguments[i]),4UL,2UL);
+      }
+      status=GaussJordanElimination(matrix,vectors,4UL,2UL);
+      matrix = RelinquishMagickMatrix(matrix, 4UL);
+      if ( status == MagickFalse ) {
+        coefficients = RelinquishMagickMemory(coefficients);
+        (void) ThrowMagickException(exception,GetMagickModule(),OptionError,
+                     "InvalidArgument","%s : '%s'","distort Bilinear",
+                     "Unsolvable Matrix");
+        return((double *) NULL);
+      }
+      return(coefficients);
     }
     case ScaleRotateTranslateDistortion:
@@ -530 +624 @@
       switch ( number_arguments ) {
       case 0: