]> granicus.if.org Git - imagemagick/blobdiff - magick/morphology.c
(no commit message)
[imagemagick] / magick / morphology.c
index 008dea0d98aa18bc897518edf39187d4051887ee..db3175c8a39ddd783071438fe593dd6e6dbb30ad 100644 (file)
@@ -17,7 +17,7 @@
 %                               January 2010                                  %
 %                                                                             %
 %                                                                             %
-%  Copyright 1999-2010 ImageMagick Studio LLC, a non-profit organization      %
+%  Copyright 1999-2011 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  %
 #include "magick/string_.h"
 #include "magick/string-private.h"
 #include "magick/token.h"
+#include "magick/utility.h"
 \f
-/*
-  The following test is for special floating point numbers of value NaN (not
-  a number), that may be used within a Kernel Definition.  NaN's are defined
-  as part of the IEEE standard for floating point number representation.
-
-  These are used a Kernel value of NaN means that that kernel position is not
-  part of the normal convolution or morphology process, and thus allowing the
-  use of 'shaped' kernels.
 
-  Special properities two NaN's are never equal, even if they are from the
-  same variable That is the IsNaN() macro is only true if the value is NaN.
+/*
+** The following test is for special floating point numbers of value NaN (not
+** a number), that may be used within a Kernel Definition.  NaN's are defined
+** as part of the IEEE standard for floating point number representation.
+**
+** These are used as a Kernel value to mean that this kernel position is not
+** part of the kernel neighbourhood for convolution or morphology processing,
+** and thus should be ignored.  This allows the use of 'shaped' kernels.
+**
+** The special properity that two NaN's are never equal, even if they are from
+** the same variable allow you to test if a value is special NaN value.
+**
+** This macro  IsNaN() is thus is only true if the value given is NaN.
 */
 #define IsNan(a)   ((a)!=(a))
 
@@ -109,7 +113,8 @@ static inline double MagickMax(const double x,const double y)
 /* Currently these are only internal to this module */
 static void
   CalcKernelMetaData(KernelInfo *),
-  ExpandKernelInfo(KernelInfo *, const double),
+  ExpandMirrorKernelInfo(KernelInfo *),
+  ExpandRotateKernelInfo(KernelInfo *, const double),
   RotateKernelInfo(KernelInfo *, double);
 \f
 
@@ -121,7 +126,6 @@ static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
   return(kernel);
 }
 
-
 /*
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %                                                                             %
@@ -160,21 +164,17 @@ static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
 %
 %  Input kernel defintion strings can consist of any of three types.
 %
-%    "name:args"
+%    "name:args[[@><]"
 %         Select from one of the built in kernels, using the name and
 %         geometry arguments supplied.  See AcquireKernelBuiltIn()
 %
-%    "WxH[+X+Y][^@]:num, num, num ..."
+%    "WxH[+X+Y][@><]:num, num, num ..."
 %         a kernel of size W by H, with W*H floating point numbers following.
 %         the 'center' can be optionally be defined at +X+Y (such that +0+0
 %         is top left corner). If not defined the pixel in the center, for
 %         odd sizes, or to the immediate top or left of center for even sizes
 %         is automatically selected.
 %
-%         If a '^' is included the kernel expanded with 90-degree rotations,
-%         While a '@' will allow you to expand a 3x3 kernel using 45-degree
-%         circular rotates.
-%
 %    "num, num, num, num, ..."
 %         list of floating point numbers defining an 'old style' odd sized
 %         square kernel.  At least 9 values should be provided for a 3x3
@@ -182,13 +182,20 @@ static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
 %         Values can be space or comma separated.  This is not recommended.
 %
 %  You can define a 'list of kernels' which can be used by some morphology
-%  operators A list is defined as a semi-colon seperated list kernels.
+%  operators A list is defined as a semi-colon separated list kernels.
 %
 %     " kernel ; kernel ; kernel ; "
 %
-%  Any extra ';' characters (at start, end or between kernel defintions are
+%  Any extra ';' charactersat start, end or between kernel defintions are
 %  simply ignored.
 %
+%  The special flags will expand a single kernel, into a list of rotated
+%  kernels. A '@' flag will expand a 3x3 kernel into a list of 45-degree
+%  cyclic rotations, while a '>' will generate a list of 90-degree rotations.
+%  The '<' also exands using 90-degree rotates, but giving a 180-degree
+%  reflected kernel before the +/- 90-degree rotations, which can be important
+%  for Thinning operations.
+%
 %  Note that 'name' kernels will start with an alphabetic character while the
 %  new kernel specification has a ':' character in its specification string.
 %  If neither is the case, it is assumed an old style of a simple list of
@@ -219,7 +226,7 @@ static KernelInfo *ParseKernelArray(const char *kernel_string)
     *p,
     *end;
 
-  register long
+  register ssize_t
     i;
 
   double
@@ -246,7 +253,7 @@ static KernelInfo *ParseKernelArray(const char *kernel_string)
   if ( end == (char *) NULL )
     end = strchr(kernel_string, '\0');
 
-  /* clear flags - for Expanding kernal lists thorugh rotations */
+  /* clear flags - for Expanding kernel lists thorugh rotations */
    flags = NoValue;
 
   /* Has a ':' in argument - New user kernel specification */
@@ -266,18 +273,18 @@ static KernelInfo *ParseKernelArray(const char *kernel_string)
          args.rho = 1.0;               /* then  width = 1 */
       if ( args.sigma < 1.0 )          /* if height too small */
         args.sigma = args.rho;         /* then  height = width */
-      kernel->width = (unsigned long)args.rho;
-      kernel->height = (unsigned long)args.sigma;
+      kernel->width = (size_t)args.rho;
+      kernel->height = (size_t)args.sigma;
 
       /* Offset Handling and Checks */
       if ( args.xi  < 0.0 || args.psi < 0.0 )
         return(DestroyKernelInfo(kernel));
-      kernel->x = ((flags & XValue)!=0) ? (long)args.xi
-                                               : (long) (kernel->width-1)/2;
-      kernel->y = ((flags & YValue)!=0) ? (long)args.psi
-                                               : (long) (kernel->height-1)/2;
-      if ( kernel->x >= (long) kernel->width ||
-           kernel->y >= (long) kernel->height )
+      kernel->x = ((flags & XValue)!=0) ? (ssize_t)args.xi
+                                               : (ssize_t) (kernel->width-1)/2;
+      kernel->y = ((flags & YValue)!=0) ? (ssize_t)args.psi
+                                               : (ssize_t) (kernel->height-1)/2;
+      if ( kernel->x >= (ssize_t) kernel->width ||
+           kernel->y >= (ssize_t) kernel->height )
         return(DestroyKernelInfo(kernel));
 
       p++; /* advance beyond the ':' */
@@ -295,8 +302,8 @@ static KernelInfo *ParseKernelArray(const char *kernel_string)
           GetMagickToken(p,&p,token);
       }
       /* set the size of the kernel - old sized square */
-      kernel->width = kernel->height= (unsigned long) sqrt((double) i+1.0);
-      kernel->x = kernel->y = (long) (kernel->width-1)/2;
+      kernel->width = kernel->height= (size_t) sqrt((double) i+1.0);
+      kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
       p=(const char *) kernel_string;
       while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
         p++;  /* ignore "'" chars for convolve filter usage - Cristy */
@@ -312,7 +319,7 @@ static KernelInfo *ParseKernelArray(const char *kernel_string)
   kernel->maximum = -MagickHuge;
   kernel->negative_range = kernel->positive_range = 0.0;
 
-  for (i=0; (i < (long) (kernel->width*kernel->height)) && (p < end); i++)
+  for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
   {
     GetMagickToken(p,&p,token);
     if (*token == ',')
@@ -322,7 +329,7 @@ static KernelInfo *ParseKernelArray(const char *kernel_string)
       kernel->values[i] = nan; /* do not include this value in kernel */
     }
     else {
-      kernel->values[i] = StringToDouble(token);
+      kernel->values[i] = StringToDouble(token,(char **) NULL);
       ( kernel->values[i] < 0)
           ?  ( kernel->negative_range += kernel->values[i] )
           :  ( kernel->positive_range += kernel->values[i] );
@@ -338,15 +345,15 @@ static KernelInfo *ParseKernelArray(const char *kernel_string)
 
 #if 0
   /* this was the old method of handling a incomplete kernel */
-  if ( i < (long) (kernel->width*kernel->height) ) {
+  if ( i < (ssize_t) (kernel->width*kernel->height) ) {
     Minimize(kernel->minimum, kernel->values[i]);
     Maximize(kernel->maximum, kernel->values[i]);
-    for ( ; i < (long) (kernel->width*kernel->height); i++)
+    for ( ; i < (ssize_t) (kernel->width*kernel->height); i++)
       kernel->values[i]=0.0;
   }
 #else
   /* Number of values for kernel was not enough - Report Error */
-  if ( i < (long) (kernel->width*kernel->height) )
+  if ( i < (ssize_t) (kernel->width*kernel->height) )
     return(DestroyKernelInfo(kernel));
 #endif
 
@@ -355,37 +362,39 @@ static KernelInfo *ParseKernelArray(const char *kernel_string)
     return(DestroyKernelInfo(kernel));
 
   if ( (flags & AreaValue) != 0 )         /* '@' symbol in kernel size */
-    ExpandKernelInfo(kernel, 45.0);
-  else if ( (flags & MinimumValue) != 0 ) /* '^' symbol in kernel size */
-    ExpandKernelInfo(kernel, 90.0);
+    ExpandRotateKernelInfo(kernel, 45.0); /* cyclic rotate 3x3 kernels */
+  else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
+    ExpandRotateKernelInfo(kernel, 90.0); /* 90 degree rotate of kernel */
+  else if ( (flags & LessValue) != 0 )    /* '<' symbol in kernel args */
+    ExpandMirrorKernelInfo(kernel);       /* 90 degree mirror rotate */
 
   return(kernel);
 }
 
 static KernelInfo *ParseKernelName(const char *kernel_string)
 {
-  KernelInfo
-    *kernel;
-
   char
     token[MaxTextExtent];
 
-  long
-    type;
-
   const char
     *p,
     *end;
 
+  GeometryInfo
+    args;
+
+  KernelInfo
+    *kernel;
+
   MagickStatusType
     flags;
 
-  GeometryInfo
-    args;
+  ssize_t
+    type;
 
   /* Parse special 'named' kernel */
   GetMagickToken(kernel_string,&p,token);
-  type=ParseMagickOption(MagickKernelOptions,MagickFalse,token);
+  type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
   if ( type < 0 || type == UserDefinedKernel )
     return((KernelInfo *)NULL);  /* not a valid named kernel */
 
@@ -411,33 +420,40 @@ static KernelInfo *ParseKernelName(const char *kernel_string)
 
   /* special handling of missing values in input string */
   switch( type ) {
-    case RectangleKernel:
-      if ( (flags & WidthValue) == 0 ) /* if no width then */
-        args.rho = args.sigma;         /* then  width = height */
-      if ( args.rho < 1.0 )            /* if width too small */
-          args.rho = 3;                 /* then  width = 3 */
-      if ( args.sigma < 1.0 )          /* if height too small */
-        args.sigma = args.rho;         /* then  height = width */
-      if ( (flags & XValue) == 0 )     /* center offset if not defined */
-        args.xi = (double)(((long)args.rho-1)/2);
-      if ( (flags & YValue) == 0 )
-        args.psi = (double)(((long)args.sigma-1)/2);
+    /* Shape Kernel Defaults */
+    case UnityKernel:
+      if ( (flags & WidthValue) == 0 )
+        args.rho = 1.0;    /* Default scale = 1.0, zero is valid */
       break;
     case SquareKernel:
     case DiamondKernel:
+    case OctagonKernel:
     case DiskKernel:
     case PlusKernel:
     case CrossKernel:
-      /* If no scale given (a 0 scale is valid! - set it to 1.0 */
       if ( (flags & HeightValue) == 0 )
-        args.sigma = 1.0;
+        args.sigma = 1.0;    /* Default scale = 1.0, zero is valid */
       break;
     case RingKernel:
       if ( (flags & XValue) == 0 )
-        args.xi = 1.0;
+        args.xi = 1.0;       /* Default scale = 1.0, zero is valid */
+      break;
+    case RectangleKernel:    /* Rectangle - set size defaults */
+      if ( (flags & WidthValue) == 0 ) /* if no width then */
+        args.rho = args.sigma;         /* then  width = height */
+      if ( args.rho < 1.0 )            /* if width too small */
+          args.rho = 3;                /* then  width = 3 */
+      if ( args.sigma < 1.0 )          /* if height too small */
+        args.sigma = args.rho;         /* then  height = width */
+      if ( (flags & XValue) == 0 )     /* center offset if not defined */
+        args.xi = (double)(((ssize_t)args.rho-1)/2);
+      if ( (flags & YValue) == 0 )
+        args.psi = (double)(((ssize_t)args.sigma-1)/2);
       break;
+    /* Distance Kernel Defaults */
     case ChebyshevKernel:
-    case ManhattenKernel:
+    case ManhattanKernel:
+    case OctagonalKernel:
     case EuclideanKernel:
       if ( (flags & HeightValue) == 0 )           /* no distance scale */
         args.sigma = 100.0;                       /* default distance scaling */
@@ -451,13 +467,17 @@ static KernelInfo *ParseKernelName(const char *kernel_string)
   }
 
   kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args);
+  if ( kernel == (KernelInfo *) NULL )
+    return(kernel);
 
   /* global expand to rotated kernel list - only for single kernels */
   if ( kernel->next == (KernelInfo *) NULL ) {
     if ( (flags & AreaValue) != 0 )         /* '@' symbol in kernel args */
-      ExpandKernelInfo(kernel, 45.0);
-    else if ( (flags & MinimumValue) != 0 ) /* '^' symbol in kernel args */
-      ExpandKernelInfo(kernel, 90.0);
+      ExpandRotateKernelInfo(kernel, 45.0);
+    else if ( (flags & GreaterValue) != 0 ) /* '>' symbol in kernel args */
+      ExpandRotateKernelInfo(kernel, 90.0);
+    else if ( (flags & LessValue) != 0 )    /* '<' symbol in kernel args */
+      ExpandMirrorKernelInfo(kernel);
   }
 
   return(kernel);
@@ -476,7 +496,7 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
   const char
     *p;
 
-  unsigned long
+  size_t
     kernel_number;
 
   p = kernel_string;
@@ -485,7 +505,7 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 
   while ( GetMagickToken(p,NULL,token),  *token != '\0' ) {
 
-    /* ignore extra or multiple ';' kernel seperators */
+    /* ignore extra or multiple ';' kernel separators */
     if ( *token != ';' ) {
 
       /* tokens starting with alpha is a Named kernel */
@@ -496,7 +516,8 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 
       /* Error handling -- this is not proper error handling! */
       if ( new_kernel == (KernelInfo *) NULL ) {
-        fprintf(stderr, "Failed to parse kernel number #%lu\n", kernel_number);
+        fprintf(stderr, "Failed to parse kernel number #%.20g\n",(double)
+          kernel_number);
         if ( kernel != (KernelInfo *) NULL )
           kernel=DestroyKernelInfo(kernel);
         return((KernelInfo *) NULL);
@@ -553,8 +574,7 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %  Convolution Kernels
 %
 %    Unity
-%       the No-Op kernel, also requivelent to  Gaussian of sigma zero.
-%       Basically a 3x3 kernel of a 1 surrounded by zeros.
+%       The a No-Op or Scaling single element kernel.
 %
 %    Gaussian:{radius},{sigma}
 %       Generate a two-dimentional gaussian kernel, as used by -gaussian.
@@ -570,19 +590,19 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %       radius will be determined so as to produce the best minimal error
 %       result, which is usally much larger than is normally needed.
 %
-%    DOG:{radius},{sigma1},{sigma2}
+%    LoG:{radius},{sigma}
+%        "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
+%        The supposed ideal edge detection, zero-summing kernel.
+%
+%        An alturnative to this kernel is to use a "DoG" with a sigma ratio of
+%        approx 1.6 (according to wikipedia).
+%
+%    DoG:{radius},{sigma1},{sigma2}
 %        "Difference of Gaussians" Kernel.
 %        As "Gaussian" but with a gaussian produced by 'sigma2' subtracted
 %        from the gaussian produced by 'sigma1'. Typically sigma2 > sigma1.
 %        The result is a zero-summing kernel.
 %
-%    LOG:{radius},{sigma}
-%        "Laplacian of a Gaussian" or "Mexician Hat" Kernel.
-%        The supposed ideal edge detection, zero-summing kernel.
-%
-%        An alturnative to this kernel is to use a "DOG" with a sigma ratio of
-%        approx 1.6, which can also be applied as a 2 pass "DOB" (see below).
-%
 %    Blur:{radius},{sigma}[,{angle}]
 %       Generates a 1 dimensional or linear gaussian blur, at the angle given
 %       (current restricted to orthogonal angles).  If a 'radius' is given the
@@ -592,19 +612,10 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %       If 'sigma' is zero, you get a single pixel on a field of zeros.
 %
 %       Note that two convolutions with two "Blur" kernels perpendicular to
-%       each other, is equivelent to a far larger "Gaussian" kernel with the
+%       each other, is equivalent to a far larger "Gaussian" kernel with the
 %       same sigma value, However it is much faster to apply. This is how the
 %       "-blur" operator actually works.
 %
-%    DOB:{radius},{sigma1},{sigma2}[,{angle}]
-%        "Difference of Blurs" Kernel.
-%        As "Blur" but with the 1D gaussian produced by 'sigma2' subtracted
-%        from thethe 1D gaussian produced by 'sigma1'.
-%        The result is a  zero-summing kernel.
-%
-%        This can be used to generate a faster "DOG" convolution, in the same
-%        way "Blur" can.
-%
 %    Comet:{width},{sigma},{angle}
 %       Blur in one direction only, much like how a bright object leaves
 %       a comet like trail.  The Kernel is actually half a gaussian curve,
@@ -640,82 +651,131 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %        Type 3 :  3x3 with center:4 edge:-2 corner:1
 %        Type 5 :  5x5 laplacian
 %        Type 7 :  7x7 laplacian
-%        Type 15 : 5x5 LOG (sigma approx 1.4)
-%        Type 19 : 9x9 LOG (sigma approx 1.4)
+%        Type 15 : 5x5 LoG (sigma approx 1.4)
+%        Type 19 : 9x9 LoG (sigma approx 1.4)
 %
 %    Sobel:{angle}
 %      Sobel 'Edge' convolution kernel (3x3)
-%           -1, 0, 1
-%           -2, 0,-2
-%           -1, 0, 1
+%          | -1, 0, 1 |
+%          | -2, 0,-2 |
+%          | -1, 0, 1 |
 %
 %    Roberts:{angle}
 %      Roberts convolution kernel (3x3)
-%            0, 0, 0
-%           -1, 1, 0
-%            0, 0, 0
+%          |  0, 0, 0 |
+%          | -1, 1, 0 |
+%          |  0, 0, 0 |
+%
 %    Prewitt:{angle}
 %      Prewitt Edge convolution kernel (3x3)
-%           -1, 0, 1
-%           -1, 0, 1
-%           -1, 0, 1
+%          | -1, 0, 1 |
+%          | -1, 0, 1 |
+%          | -1, 0, 1 |
+%
 %    Compass:{angle}
 %      Prewitt's "Compass" convolution kernel (3x3)
-%           -1, 1, 1
-%           -1,-2, 1
-%           -1, 1, 1
+%          | -1, 1, 1 |
+%          | -1,-2, 1 |
+%          | -1, 1, 1 |
+%
 %    Kirsch:{angle}
 %      Kirsch's "Compass" convolution kernel (3x3)
-%           -3,-3, 5
-%           -3, 0, 5
-%           -3,-3, 5
+%          | -3,-3, 5 |
+%          | -3, 0, 5 |
+%          | -3,-3, 5 |
+%
+%    FreiChen:{angle}
+%      Frei-Chen Edge Detector is based on a kernel that is similar to
+%      the Sobel Kernel, but is designed to be isotropic. That is it takes
+%      into account the distance of the diagonal in the kernel.
+%
+%          |   1,     0,   -1     |
+%          | sqrt(2), 0, -sqrt(2) |
+%          |   1,     0,   -1     |
 %
 %    FreiChen:{type},{angle}
-%      Frei-Chen Edge Detector is a set of 9 unique convolution kernels that
-%      are specially weighted.  They should not be normalized. After applying
-%      each to the original image, the results is then added together.  The
-%      square root of the resulting image is the cosine of the edge, and the
-%      direction of the feature detection.
 %
-%        Type 1: |  1,   sqrt(2),  1 |
-%                |  0,     0,      0 | / 2*sqrt(2)
-%                | -1,  -sqrt(2), -1 |
+%      Frei-Chen Pre-weighted kernels...
+%
+%        Type 0:  default un-nomalized version shown above.
+%
+%        Type 1: Orthogonal Kernel (same as type 11 below)
+%          |   1,     0,   -1     |
+%          | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
+%          |   1,     0,   -1     |
+%
+%        Type 2: Diagonal form of Kernel...
+%          |   1,     sqrt(2),    0     |
+%          | sqrt(2),   0,     -sqrt(2) | / 2*sqrt(2)
+%          |   0,    -sqrt(2)    -1     |
 %
-%        Type 2: |   1,     0,   1     |
-%                | sqrt(2), 0, sqrt(2) | / 2*sqrt(2)
-%                |   1,     0,   1     |
+%      However this kernel is als at the heart of the FreiChen Edge Detection
+%      Process which uses a set of 9 specially weighted kernel.  These 9
+%      kernels not be normalized, but directly applied to the image. The
+%      results is then added together, to produce the intensity of an edge in
+%      a specific direction.  The square root of the pixel value can then be
+%      taken as the cosine of the edge, and at least 2 such runs at 90 degrees
+%      from each other, both the direction and the strength of the edge can be
+%      determined.
 %
-%        Type 3: |    0,    -1, sqrt(2) |
-%                |    1,     0,   -1    | / 2*sqrt(2)
-%                | -sqrt(2), 1,    0    |
+%        Type 10: All 9 of the following pre-weighted kernels...
 %
-%        Type 4: | sqrt(2), -1,     0    |
-%                |   -1,     0,     1    | / 2*sqrt(2)
-%                |    0,     1, -sqrt(2) |
+%        Type 11: |   1,     0,   -1     |
+%                 | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
+%                 |   1,     0,   -1     |
 %
-%        Type 5: |  0, 1,  0 |
-%                | -1, 0, -1 | / 2
-%                |  0, 1,  0 |
+%        Type 12: | 1, sqrt(2), 1 |
+%                 | 0,   0,     0 | / 2*sqrt(2)
+%                 | 1, sqrt(2), 1 |
 %
-%        Type 6: | -1, 0,  1 |
-%                |  0, 0,  0 | / 2
-%                |  1, 0, -1 |
+%        Type 13: | sqrt(2), -1,    0     |
+%                 |  -1,      0,    1     | / 2*sqrt(2)
+%                 |   0,      1, -sqrt(2) |
 %
-%        Type 7: |  1, -2,  1 |
-%                | -2,  4, -2 | / 6
-%                |  1, -2,  1 |
+%        Type 14: |    0,     1, -sqrt(2) |
+%                 |   -1,     0,     1    | / 2*sqrt(2)
+%                 | sqrt(2), -1,     0    |
 %
-%        Type 8: | -2, 1, -2 |
-%                |  1, 4,  1 | / 6
-%                | -2, 1, -2 |
+%        Type 15: | 0, -1, 0 |
+%                 | 1,  0, 1 | / 2
+%                 | 0, -1, 0 |
 %
-%        Type 9: | 1, 1, 1 |
-%                | 1, 1, 1 | / 3
-%                | 1, 1, 1 |
+%        Type 16: |  1, 0, -1 |
+%                 |  0, 0,  0 | / 2
+%                 | -1, 0,  1 |
+%
+%        Type 17: |  1, -2,  1 |
+%                 | -2,  4, -2 | / 6
+%                 | -1, -2,  1 |
+%
+%        Type 18: | -2, 1, -2 |
+%                 |  1, 4,  1 | / 6
+%                 | -2, 1, -2 |
+%
+%        Type 19: | 1, 1, 1 |
+%                 | 1, 1, 1 | / 3
+%                 | 1, 1, 1 |
 %
 %      The first 4 are for edge detection, the next 4 are for line detection
 %      and the last is to add a average component to the results.
 %
+%      Using a special type of '-1' will return all 9 pre-weighted kernels
+%      as a multi-kernel list, so that you can use them directly (without
+%      normalization) with the special "-set option:morphology:compose Plus"
+%      setting to apply the full FreiChen Edge Detection Technique.
+%
+%      If 'type' is large it will be taken to be an actual rotation angle for
+%      the default FreiChen (type 0) kernel.  As such  FreiChen:45  will look
+%      like a  Sobel:45  but with 'sqrt(2)' instead of '2' values.
+%
+%      WARNING: The above was layed out as per
+%          http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf
+%      But rotated 90 degrees so direction is from left rather than the top.
+%      I have yet to find any secondary confirmation of the above. The only
+%      other source found was actual source code at
+%          http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf
+%      Neigher paper defineds the kernels in a way that looks locical or
+%      correct when taken as a whole.
 %
 %  Boolean Kernels
 %
@@ -728,10 +788,31 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %       Generate a square shaped kernel of size radius*2+1, and defaulting
 %       to a 3x3 (radius 1).
 %
-%       Note that using a larger radius for the "Square" or the "Diamond" is
-%       also equivelent to iterating the basic morphological method that many
-%       times. However iterating with the smaller radius is actually faster
-%       than using a larger kernel radius.
+%    Octagon:[{radius}[,{scale}]]
+%       Generate octagonal shaped kernel of given radius and constant scale.
+%       Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
+%       in "Diamond" kernel.
+%
+%    Disk:[{radius}[,{scale}]]
+%       Generate a binary disk, thresholded at the radius given, the radius
+%       may be a float-point value. Final Kernel size is floor(radius)*2+1
+%       square. A radius of 5.3 is the default.
+%
+%       NOTE: That a low radii Disk kernels produce the same results as
+%       many of the previously defined kernels, but differ greatly at larger
+%       radii.  Here is a table of equivalences...
+%          "Disk:1"    => "Diamond", "Octagon:1", or "Cross:1"
+%          "Disk:1.5"  => "Square"
+%          "Disk:2"    => "Diamond:2"
+%          "Disk:2.5"  => "Octagon"
+%          "Disk:2.9"  => "Square:2"
+%          "Disk:3.5"  => "Octagon:3"
+%          "Disk:4.5"  => "Octagon:4"
+%          "Disk:5.4"  => "Octagon:5"
+%          "Disk:6.4"  => "Octagon:6"
+%       All other Disk shapes are unique to this kernel, but because a "Disk"
+%       is more circular when using a larger radius, using a larger radius is
+%       preferred over iterating the morphological operation.
 %
 %    Rectangle:{geometry}
 %       Simply generate a rectangle of 1's with the size given. You can also
@@ -740,23 +821,6 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %
 %       Properly centered and odd sized rectangles work the best.
 %
-%    Disk:[{radius}[,{scale}]]
-%       Generate a binary disk of the radius given, radius may be a float.
-%       Kernel size will be ceil(radius)*2+1 square.
-%       NOTE: Here are some disk shapes of specific interest
-%          "Disk:1"    => "diamond" or "cross:1"
-%          "Disk:1.5"  => "square"
-%          "Disk:2"    => "diamond:2"
-%          "Disk:2.5"  => a general disk shape of radius 2
-%          "Disk:2.9"  => "square:2"
-%          "Disk:3.5"  => default - octagonal/disk shape of radius 3
-%          "Disk:4.2"  => roughly octagonal shape of radius 4
-%          "Disk:4.3"  => a general disk shape of radius 4
-%       After this all the kernel shape becomes more and more circular.
-%
-%       Because a "disk" is more circular when using a larger radius, using a
-%       larger radius is preferred over iterating the morphological operation.
-%
 %  Symbol Dilation Kernels
 %
 %    These kernel is not a good general morphological kernel, but is used
@@ -771,7 +835,7 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %       Generate a kernel in the shape of a 'plus' or a 'cross' with
 %       a each arm the length of the given radius (default 2).
 %
-%       NOTE: "plus:1" is equivelent to a "Diamond" kernel.
+%       NOTE: "plus:1" is equivalent to a "Diamond" kernel.
 %
 %    Ring:{radius1},{radius2}[,{scale}]
 %       A ring of the values given that falls between the two radii.
@@ -785,23 +849,42 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %       Find any peak larger than the pixels the fall between the two radii.
 %       The default ring of pixels is as per "Ring".
 %    Edges
-%       Find edges of a binary shape
+%       Find flat orthogonal edges of a binary shape
 %    Corners
-%       Find corners of a binary shape
-%    Ridges
-%       Find single pixel ridges or thin lines
-%    Ridges2
-%       Find 2 pixel thick ridges or lines
-%    Ridges3
-%       Find 2 pixel thick diagonal ridges (experimental)
-%    LineEnds
+%       Find 90 degree corners of a binary shape
+%    Diagonals:type
+%       A special kernel to thin the 'outside' of diagonals
+%    LineEnds:type
 %       Find end points of lines (for pruning a skeletion)
+%       Two types of lines ends (default to both) can be searched for
+%         Type 0: All line ends
+%         Type 1: single kernel for 4-conneected line ends
+%         Type 2: single kernel for simple line ends
 %    LineJunctions
 %       Find three line junctions (within a skeletion)
+%         Type 0: all line junctions
+%         Type 1: Y Junction kernel
+%         Type 2: Diagonal T Junction kernel
+%         Type 3: Orthogonal T Junction kernel
+%         Type 4: Diagonal X Junction kernel
+%         Type 5: Orthogonal + Junction kernel
+%    Ridges:type
+%       Find single pixel ridges or thin lines
+%         Type 1: Fine single pixel thick lines and ridges
+%         Type 2: Find two pixel thick lines and ridges
 %    ConvexHull
-%       Octagonal thicken kernel, to generate convex hulls of 45 degrees
-%    Skeleton
-%       Thinning kernel, which leaves behind a skeletion of a shape
+%       Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
+%    Skeleton:type
+%       Traditional skeleton generating kernels.
+%         Type 1: Tradional Skeleton kernel (4 connected skeleton)
+%         Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
+%         Type 3: Thinning skeleton based on a ressearch paper by
+%                 Dan S. Bloomberg (Default Type)
+%    ThinSE:type
+%       A huge variety of Thinning Kernels designed to preserve conectivity.
+%       many other kernel sets use these kernels as source definitions.
+%       Type numbers are 41-49, 81-89, 481, and 482 which are based on
+%       the super and sub notations used in the source research paper.
 %
 %  Distance Measuring Kernels
 %
@@ -814,40 +897,45 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %    applied.
 %
 %    Chebyshev:[{radius}][x{scale}[%!]]
-%       Chebyshev Distance (also known as Tchebychev Distance) is a value of
-%       one to any neighbour, orthogonal or diagonal. One why of thinking of
-%       it is the number of squares a 'King' or 'Queen' in chess needs to
-%       traverse reach any other position on a chess board.  It results in a
-%       'square' like distance function, but one where diagonals are closer
-%       than expected.
-%
-%    Manhatten:[{radius}][x{scale}[%!]]
-%       Manhatten Distance (also known as Rectilinear Distance, or the Taxi
-%       Cab metric), is the distance needed when you can only travel in
-%       orthogonal (horizontal or vertical) only.  It is the distance a 'Rook'
-%       in chess would travel. It results in a diamond like distances, where
-%       diagonals are further than expected.
+%       Chebyshev Distance (also known as Tchebychev or Chessboard distance)
+%       is a value of one to any neighbour, orthogonal or diagonal. One why
+%       of thinking of it is the number of squares a 'King' or 'Queen' in
+%       chess needs to traverse reach any other position on a chess board.
+%       It results in a 'square' like distance function, but one where
+%       diagonals are given a value that is closer than expected.
+%
+%    Manhattan:[{radius}][x{scale}[%!]]
+%       Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
+%       Cab distance metric), it is the distance needed when you can only
+%       travel in horizontal or vertical directions only.  It is the
+%       distance a 'Rook' in chess would have to travel, and results in a
+%       diamond like distances, where diagonals are further than expected.
+%
+%    Octagonal:[{radius}][x{scale}[%!]]
+%       An interleving of Manhatten and Chebyshev metrics producing an
+%       increasing octagonally shaped distance.  Distances matches those of
+%       the "Octagon" shaped kernel of the same radius.  The minimum radius
+%       and default is 2, producing a 5x5 kernel.
 %
 %    Euclidean:[{radius}][x{scale}[%!]]
-%       Euclidean Distance is the 'direct' or 'as the crow flys distance.
+%       Euclidean distance is the 'direct' or 'as the crow flys' distance.
 %       However by default the kernel size only has a radius of 1, which
 %       limits the distance to 'Knight' like moves, with only orthogonal and
 %       diagonal measurements being correct.  As such for the default kernel
-%       you will get octagonal like distance function, which is reasonally
-%       accurate.
-%
-%       However if you use a larger radius such as "Euclidean:4" you will
-%       get a much smoother distance gradient from the edge of the shape.
-%       Of course a larger kernel is slower to use, and generally not needed.
-%
-%       To allow the use of fractional distances that you get with diagonals
-%       the actual distance is scaled by a fixed value which the user can
-%       provide.  This is not actually nessary for either ""Chebyshev" or
-%       "Manhatten" distance kernels, but is done for all three distance
-%       kernels.  If no scale is provided it is set to a value of 100,
-%       allowing for a maximum distance measurement of 655 pixels using a Q16
-%       version of IM, from any edge.  However for small images this can
-%       result in quite a dark gradient.
+%       you will get octagonal like distance function.
+%
+%       However using a larger radius such as "Euclidean:4" you will get a
+%       much smoother distance gradient from the edge of the shape. Especially
+%       if the image is pre-processed to include any anti-aliasing pixels.
+%       Of course a larger kernel is slower to use, and not always needed.
+%
+%    The first three Distance Measuring Kernels will only generate distances
+%    of exact multiples of {scale} in binary images. As such you can use a
+%    scale of 1 without loosing any information.  However you also need some
+%    scaling when handling non-binary anti-aliased shapes.
+%
+%    The "Euclidean" Distance Kernel however does generate a non-integer
+%    fractional results, and as such scaling is vital even for binary shapes.
 %
 */
 
@@ -857,10 +945,10 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
   KernelInfo
     *kernel;
 
-  register long
+  register ssize_t
     i;
 
-  register long
+  register ssize_t
     u,
     v;
 
@@ -870,41 +958,51 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
   /* Generate a new empty kernel if needed */
   kernel=(KernelInfo *) NULL;
   switch(type) {
-    case UndefinedKernel:      /* These should not be used here */
+    case UndefinedKernel:    /* These should not call this function */
     case UserDefinedKernel:
+      assert("Should not call this function" != (char *)NULL);
       break;
-    case LaplacianKernel:  /* Named Descrete Convolution Kernels */
-    case SobelKernel:
+    case LaplacianKernel:   /* Named Descrete Convolution Kernels */
+    case SobelKernel:       /* these are defined using other kernels */
     case RobertsKernel:
     case PrewittKernel:
     case CompassKernel:
     case KirschKernel:
-    case CornersKernel:    /* Hit and Miss kernels */
+    case FreiChenKernel:
+    case EdgesKernel:       /* Hit and Miss kernels */
+    case CornersKernel:
+    case DiagonalsKernel:
     case LineEndsKernel:
     case LineJunctionsKernel:
+    case RidgesKernel:
     case ConvexHullKernel:
     case SkeletonKernel:
-      /* A pre-generated kernel is not needed */
-      break;
+    case ThinSEKernel:
+      break;               /* A pre-generated kernel is not needed */
 #if 0
+    /* set to 1 to do a compile-time check that we haven't missed anything */
+    case UnityKernel:
     case GaussianKernel:
-    case DOGKernel:
+    case DoGKernel:
+    case LoGKernel:
     case BlurKernel:
-    case DOBKernel:
     case CometKernel:
     case DiamondKernel:
     case SquareKernel:
     case RectangleKernel:
+    case OctagonKernel:
     case DiskKernel:
     case PlusKernel:
     case CrossKernel:
     case RingKernel:
     case PeaksKernel:
     case ChebyshevKernel:
-    case ManhattenKernel:
+    case ManhattanKernel:
+    case OctangonalKernel:
     case EuclideanKernel:
-#endif
+#else
     default:
+#endif
       /* Generate the base Kernel Structure */
       kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
       if (kernel == (KernelInfo *) NULL)
@@ -919,23 +1017,36 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
   }
 
   switch(type) {
-    /* Convolution Kernels */
+    /*
+      Convolution Kernels
+    */
+    case UnityKernel:
+      {
+        kernel->height = kernel->width = (size_t) 1;
+        kernel->x = kernel->y = (ssize_t) 0;
+        kernel->values=(double *) AcquireQuantumMemory(1,sizeof(double));
+        if (kernel->values == (double *) NULL)
+          return(DestroyKernelInfo(kernel));
+        kernel->maximum = kernel->values[0] = args->rho;
+        break;
+      }
+      break;
     case GaussianKernel:
-    case DOGKernel:
-    case LOGKernel:
+    case DoGKernel:
+    case LoGKernel:
       { double
           sigma = fabs(args->sigma),
           sigma2 = fabs(args->xi),
           A, B, R;
 
         if ( args->rho >= 1.0 )
-          kernel->width = (unsigned long)args->rho*2+1;
-        else if ( (type != DOGKernel) || (sigma >= sigma2) )
+          kernel->width = (size_t)args->rho*2+1;
+        else if ( (type != DoGKernel) || (sigma >= sigma2) )
           kernel->width = GetOptimalKernelWidth2D(args->rho,sigma);
         else
           kernel->width = GetOptimalKernelWidth2D(args->rho,sigma2);
         kernel->height = kernel->width;
-        kernel->x = kernel->y = (long) (kernel->width-1)/2;
+        kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
                               kernel->height*sizeof(double));
         if (kernel->values == (double *) NULL)
@@ -944,17 +1055,17 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         /* WARNING: The following generates a 'sampled gaussian' kernel.
          * What we really want is a 'discrete gaussian' kernel.
          *
-         * How to do this is currently not known, but appears to be
-         * basied on the Error Function 'erf()' (intergral of a gaussian)
+         * How to do this is I don't know, but appears to be basied on the
+         * Error Function 'erf()' (intergral of a gaussian)
          */
 
-        if ( type == GaussianKernel || type == DOGKernel )
-          { /* Calculate a Gaussian,  OR positive half of a DOG */
+        if ( type == GaussianKernel || type == DoGKernel )
+          { /* Calculate a Gaussian,  OR positive half of a DoG */
             if ( sigma > MagickEpsilon )
               { A = 1.0/(2.0*sigma*sigma);  /* simplify loop expressions */
-                B = 1.0/(Magick2PI*sigma*sigma);
-                for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-                  for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
+                B = (double) (1.0/(Magick2PI*sigma*sigma));
+                for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+                  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
                       kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
               }
             else /* limiting case - a unity (normalized Dirac) kernel */
@@ -964,27 +1075,27 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
               }
           }
 
-        if ( type == DOGKernel )
+        if ( type == DoGKernel )
           { /* Subtract a Negative Gaussian for "Difference of Gaussian" */
             if ( sigma2 > MagickEpsilon )
               { sigma = sigma2;                /* simplify loop expressions */
                 A = 1.0/(2.0*sigma*sigma);
-                B = 1.0/(Magick2PI*sigma*sigma);
-                for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-                  for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
+                B = (double) (1.0/(Magick2PI*sigma*sigma));
+                for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+                  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
                     kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
               }
             else /* limiting case - a unity (normalized Dirac) kernel */
               kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
           }
 
-        if ( type == LOGKernel )
+        if ( type == LoGKernel )
           { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
             if ( sigma > MagickEpsilon )
               { A = 1.0/(2.0*sigma*sigma);  /* simplify loop expressions */
-                B = 1.0/(MagickPI*sigma*sigma*sigma*sigma);
-                for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-                  for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
+                B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
+                for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+                  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
                     { R = ((double)(u*u+v*v))*A;
                       kernel->values[i] = (1-R)*exp(-R)*B;
                     }
@@ -1015,20 +1126,16 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         break;
       }
     case BlurKernel:
-    case DOBKernel:
       { double
           sigma = fabs(args->sigma),
-          sigma2 = fabs(args->xi),
-          A, B;
+          alpha, beta;
 
         if ( args->rho >= 1.0 )
-          kernel->width = (unsigned long)args->rho*2+1;
-        else if ( (type == BlurKernel) || (sigma >= sigma2) )
-          kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
+          kernel->width = (size_t)args->rho*2+1;
         else
-          kernel->width = GetOptimalKernelWidth1D(args->rho,sigma2);
+          kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
         kernel->height = 1;
-        kernel->x = (long) (kernel->width-1)/2;
+        kernel->x = (ssize_t) (kernel->width-1)/2;
         kernel->y = 0;
         kernel->negative_range = kernel->positive_range = 0.0;
         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
@@ -1052,63 +1159,36 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         */
 
         /* initialize */
-        v = (long) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
+        v = (ssize_t) (kernel->width*KernelRank-1)/2; /* start/end points to fit range */
         (void) ResetMagickMemory(kernel->values,0, (size_t)
                      kernel->width*kernel->height*sizeof(double));
         /* Calculate a Positive 1D Gaussian */
         if ( sigma > MagickEpsilon )
           { sigma *= KernelRank;               /* simplify loop expressions */
-            A = 1.0/(2.0*sigma*sigma);
-            B = 1.0/(MagickSQ2PI*sigma );
+            alpha = 1.0/(2.0*sigma*sigma);
+            beta= (double) (1.0/(MagickSQ2PI*sigma ));
             for ( u=-v; u <= v; u++) {
-              kernel->values[(u+v)/KernelRank] += exp(-((double)(u*u))*A)*B;
+              kernel->values[(u+v)/KernelRank] +=
+                              exp(-((double)(u*u))*alpha)*beta;
             }
           }
         else /* special case - generate a unity kernel */
           kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
-
-        /* Subtract a Second 1D Gaussian for "Difference of Blur" */
-        if ( type == DOBKernel )
-          {
-            if ( sigma2 > MagickEpsilon )
-              { sigma = sigma2*KernelRank;      /* simplify loop expressions */
-                A = 1.0/(2.0*sigma*sigma);
-                B = 1.0/(MagickSQ2PI*sigma);
-                for ( u=-v; u <= v; u++)
-                  kernel->values[(u+v)/KernelRank] -= exp(-((double)(u*u))*A)*B;
-              }
-            else /* limiting case - a unity (normalized Dirac) kernel */
-              kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
-          }
 #else
         /* Direct calculation without curve averaging */
 
         /* Calculate a Positive Gaussian */
         if ( sigma > MagickEpsilon )
-          { A = 1.0/(2.0*sigma*sigma);     /* simplify loop expressions */
-            B = 1.0/(MagickSQ2PI*sigma);
-            for ( i=0, u=-kernel->x; u <= (long)kernel->x; u++, i++)
-              kernel->values[i] = exp(-((double)(u*u))*A)*B;
+          { alpha = 1.0/(2.0*sigma*sigma);    /* simplify loop expressions */
+            beta = 1.0/(MagickSQ2PI*sigma);
+            for ( i=0, u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
+              kernel->values[i] = exp(-((double)(u*u))*alpha)*beta;
           }
         else /* special case - generate a unity kernel */
           { (void) ResetMagickMemory(kernel->values,0, (size_t)
                          kernel->width*kernel->height*sizeof(double));
             kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
           }
-
-        /* Subtract a Second 1D Gaussian for "Difference of Blur" */
-        if ( type == DOBKernel )
-          {
-            if ( sigma2 > MagickEpsilon )
-              { sigma = sigma2;                /* simplify loop expressions */
-                A = 1.0/(2.0*sigma*sigma);
-                B = 1.0/(MagickSQ2PI*sigma);
-                for ( i=0, u=-kernel->x; u <= (long)kernel->x; u++, i++)
-                  kernel->values[i] -= exp(-((double)(u*u))*A)*B;
-              }
-            else /* limiting case - a unity (normalized Dirac) kernel */
-              kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
-          }
 #endif
         /* Note the above kernel may have been 'clipped' by a user defined
         ** radius, producing a smaller (darker) kernel.  Also for very small
@@ -1127,7 +1207,7 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);
 
         /* rotate the 1D kernel by given angle */
-        RotateKernelInfo(kernel, (type == BlurKernel) ? args->xi : args->psi );
+        RotateKernelInfo(kernel, args->xi );
         break;
       }
     case CometKernel:
@@ -1138,7 +1218,7 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         if ( args->rho < 1.0 )
           kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
         else
-          kernel->width = (unsigned long)args->rho;
+          kernel->width = (size_t)args->rho;
         kernel->x = kernel->y = 0;
         kernel->height = 1;
         kernel->negative_range = kernel->positive_range = 0.0;
@@ -1161,7 +1241,7 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
           {
 #if 1
 #define KernelRank 3
-            v = (long) kernel->width*KernelRank; /* start/end points */
+            v = (ssize_t) kernel->width*KernelRank; /* start/end points */
             (void) ResetMagickMemory(kernel->values,0, (size_t)
                           kernel->width*sizeof(double));
             sigma *= KernelRank;            /* simplify the loop expression */
@@ -1172,12 +1252,12 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
                   exp(-((double)(u*u))*A);
               /*  exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
             }
-            for (i=0; i < (long) kernel->width; i++)
+            for (i=0; i < (ssize_t) kernel->width; i++)
               kernel->positive_range += kernel->values[i];
 #else
             A = 1.0/(2.0*sigma*sigma);     /* simplify the loop expression */
             /* B = 1.0/(MagickSQ2PI*sigma); */
-            for ( i=0; i < (long) kernel->width; i++)
+            for ( i=0; i < (ssize_t) kernel->width; i++)
               kernel->positive_range +=
                 kernel->values[i] =
                   exp(-((double)(i*i))*A);
@@ -1200,7 +1280,9 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         break;
       }
 
-    /* Convolution Kernels - Well Known Constants */
+    /*
+      Convolution Kernels - Well Known Named Constant Kernels
+    */
     case LaplacianKernel:
       { switch ( (int) args->rho ) {
           case 0:
@@ -1224,14 +1306,14 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
             kernel=ParseKernelArray(
               "7:-10,-5,-2,-1,-2,-5,-10 -5,0,3,4,3,0,-5 -2,3,6,7,6,3,-2 -1,4,7,8,7,4,-1 -2,3,6,7,6,3,-2 -5,0,3,4,3,0,-5 -10,-5,-2,-1,-2,-5,-10" );
             break;
-          case 15:  /* a 5x5 LOG (sigma approx 1.4) */
+          case 15:  /* a 5x5 LoG (sigma approx 1.4) */
             kernel=ParseKernelArray(
               "5: 0,0,-1,0,0  0,-1,-2,-1,0  -1,-2,16,-2,-1  0,-1,-2,-1,0  0,0,-1,0,0");
             break;
-          case 19:  /* a 9x9 LOG (sigma approx 1.4) */
+          case 19:  /* a 9x9 LoG (sigma approx 1.4) */
             /* http://www.cscjournals.org/csc/manuscript/Journals/IJIP/volume3/Issue1/IJIP-15.pdf */
             kernel=ParseKernelArray(
-              "9: 0,-1,-1,-2,-2,-2,-1,-1,0  -1,-2,-4,-5,-5,-5,-4,-2,-1  -1,-4,-5,-3,-0,-3,-5,-4,-1  -2,-5,-3,@12,@24,@12,-3,-5,-2  -2,-5,-0,@24,@40,@24,-0,-5,-2  -2,-5,-3,@12,@24,@12,-3,-5,-2  -1,-4,-5,-3,-0,-3,-5,-4,-1  -1,-2,-4,-5,-5,-5,-4,-2,-1  0,-1,-1,-2,-2,-2,-1,-1,0");
+              "9: 0,-1,-1,-2,-2,-2,-1,-1,0  -1,-2,-4,-5,-5,-5,-4,-2,-1  -1,-4,-5,-3,-0,-3,-5,-4,-1  -2,-5,-3,12,24,12,-3,-5,-2  -2,-5,-0,24,40,24,-0,-5,-2  -2,-5,-3,12,24,12,-3,-5,-2  -1,-4,-5,-3,-0,-3,-5,-4,-1  -1,-2,-4,-5,-5,-5,-4,-2,-1  0,-1,-1,-2,-2,-2,-1,-1,0");
             break;
         }
         if (kernel == (KernelInfo *) NULL)
@@ -1240,17 +1322,17 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         break;
       }
     case SobelKernel:
-      {
-        kernel=ParseKernelArray("3: -1,0,1  -2,0,2  -1,0,1");
+      { /* Simple Sobel Kernel */
+        kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
         if (kernel == (KernelInfo *) NULL)
           return(kernel);
         kernel->type = type;
-        RotateKernelInfo(kernel, args->rho); /* Rotate by angle */
+        RotateKernelInfo(kernel, args->rho);
         break;
       }
     case RobertsKernel:
       {
-        kernel=ParseKernelArray("3: 0,0,0  -1,1,0  0,0,0");
+        kernel=ParseKernelArray("3: 0,0,0  1,-1,0  0,0,0");
         if (kernel == (KernelInfo *) NULL)
           return(kernel);
         kernel->type = type;
@@ -1259,7 +1341,7 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
       }
     case PrewittKernel:
       {
-        kernel=ParseKernelArray("3: -1,1,1  0,0,0  -1,1,1");
+        kernel=ParseKernelArray("3: 1,0,-1  1,0,-1  1,0,-1");
         if (kernel == (KernelInfo *) NULL)
           return(kernel);
         kernel->type = type;
@@ -1268,7 +1350,7 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
       }
     case CompassKernel:
       {
-        kernel=ParseKernelArray("3: -1,1,1  -1,-2,1  -1,1,1");
+        kernel=ParseKernelArray("3: 1,1,-1  1,-2,-1  1,1,-1");
         if (kernel == (KernelInfo *) NULL)
           return(kernel);
         kernel->type = type;
@@ -1277,7 +1359,7 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
       }
     case KirschKernel:
       {
-        kernel=ParseKernelArray("3: -3,-3,5  -3,0,5  -3,-3,5");
+        kernel=ParseKernelArray("3: 5,-3,-3  5,0,-3  5,-3,-3");
         if (kernel == (KernelInfo *) NULL)
           return(kernel);
         kernel->type = type;
@@ -1285,89 +1367,131 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         break;
       }
     case FreiChenKernel:
-      /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf */
-      /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf */
+      /* Direction is set to be left to right positive */
+      /* http://www.math.tau.ac.il/~turkel/notes/edge_detectors.pdf -- RIGHT? */
+      /* http://ltswww.epfl.ch/~courstiv/exos_labos/sol3.pdf -- WRONG? */
       { switch ( (int) args->rho ) {
           default:
-          case 1:
-            kernel=ParseKernelArray("3: 1,2,1  0,0,0  -1,2,-1");
+          case 0:
+            kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
             if (kernel == (KernelInfo *) NULL)
               return(kernel);
-            kernel->values[1] = +MagickSQ2;
-            kernel->values[7] = -MagickSQ2;
+            kernel->type = type;
+            kernel->values[3] = +MagickSQ2;
+            kernel->values[5] = -MagickSQ2;
             CalcKernelMetaData(kernel);     /* recalculate meta-data */
-            ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
             break;
           case 2:
-            kernel=ParseKernelArray("3: 1,0,1  2,0,2  1,0,1");
+            kernel=ParseKernelArray("3: 1,2,0  2,0,-2  0,-2,-1");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            kernel->values[1] = kernel->values[3] = +MagickSQ2;
+            kernel->values[5] = kernel->values[7] = -MagickSQ2;
+            CalcKernelMetaData(kernel);     /* recalculate meta-data */
+            ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
+            break;
+          case 10:
+            kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            break;
+          case 1:
+          case 11:
+            kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
             if (kernel == (KernelInfo *) NULL)
               return(kernel);
+            kernel->type = type;
             kernel->values[3] = +MagickSQ2;
-            kernel->values[5] = +MagickSQ2;
-            CalcKernelMetaData(kernel);
-            ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
+            kernel->values[5] = -MagickSQ2;
+            CalcKernelMetaData(kernel);     /* recalculate meta-data */
+            ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
             break;
-          case 3:
-            kernel=ParseKernelArray("3: 0,-1,2  1,0,-1  -2,1,0");
+          case 12:
+            kernel=ParseKernelArray("3: 1,2,1  0,0,0  1,2,1");
             if (kernel == (KernelInfo *) NULL)
               return(kernel);
-            kernel->values[2] = +MagickSQ2;
-            kernel->values[6] = -MagickSQ2;
+            kernel->type = type;
+            kernel->values[1] = +MagickSQ2;
+            kernel->values[7] = +MagickSQ2;
             CalcKernelMetaData(kernel);
-            ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
+            ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
             break;
-          case 4:
+          case 13:
             kernel=ParseKernelArray("3: 2,-1,0  -1,0,1  0,1,-2");
             if (kernel == (KernelInfo *) NULL)
               return(kernel);
+            kernel->type = type;
             kernel->values[0] = +MagickSQ2;
             kernel->values[8] = -MagickSQ2;
             CalcKernelMetaData(kernel);
-            ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
+            ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
+            break;
+          case 14:
+            kernel=ParseKernelArray("3: 0,1,-2  -1,0,1  2,-1,0");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            kernel->values[2] = -MagickSQ2;
+            kernel->values[6] = +MagickSQ2;
+            CalcKernelMetaData(kernel);
+            ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
             break;
-          case 5:
-            kernel=ParseKernelArray("3: 0,1,0  -1,0,-1  0,1,0");
+          case 15:
+            kernel=ParseKernelArray("3: 0,-1,0  1,0,1  0,-1,0");
             if (kernel == (KernelInfo *) NULL)
               return(kernel);
+            kernel->type = type;
             ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
             break;
-          case 6:
-            kernel=ParseKernelArray("3: -1,0,1  0,0,0  1,0,-1");
+          case 16:
+            kernel=ParseKernelArray("3: 1,0,-1  0,0,0  -1,0,1");
             if (kernel == (KernelInfo *) NULL)
               return(kernel);
+            kernel->type = type;
             ScaleKernelInfo(kernel, 1.0/2.0, NoValue);
             break;
-          case 7:
-            kernel=ParseKernelArray("3: 1,-2,1  -2,4,-2  1,-2,1");
+          case 17:
+            kernel=ParseKernelArray("3: 1,-2,1  -2,4,-2  -1,-2,1");
             if (kernel == (KernelInfo *) NULL)
               return(kernel);
+            kernel->type = type;
             ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
             break;
-          case 8:
+          case 18:
             kernel=ParseKernelArray("3: -2,1,-2  1,4,1  -2,1,-2");
             if (kernel == (KernelInfo *) NULL)
               return(kernel);
+            kernel->type = type;
             ScaleKernelInfo(kernel, 1.0/6.0, NoValue);
             break;
-          case 9:
-            kernel=ParseKernelName("3: 1,1,1  1,1,1  1,1,1");
+          case 19:
+            kernel=ParseKernelArray("3: 1,1,1  1,1,1  1,1,1");
             if (kernel == (KernelInfo *) NULL)
               return(kernel);
+            kernel->type = type;
             ScaleKernelInfo(kernel, 1.0/3.0, NoValue);
             break;
         }
-        RotateKernelInfo(kernel, args->sigma);  /* Rotate by angle */
+        if ( fabs(args->sigma) > MagickEpsilon )
+          /* Rotate by correctly supplied 'angle' */
+          RotateKernelInfo(kernel, args->sigma);
+        else if ( args->rho > 30.0 || args->rho < -30.0 )
+          /* Rotate by out of bounds 'type' */
+          RotateKernelInfo(kernel, args->rho);
         break;
       }
 
-    /* Boolean Kernels */
+    /*
+      Boolean or Shaped Kernels
+    */
     case DiamondKernel:
       {
         if (args->rho < 1.0)
           kernel->width = kernel->height = 3;  /* default radius = 1 */
         else
-          kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
-        kernel->x = kernel->y = (long) (kernel->width-1)/2;
+          kernel->width = kernel->height = ((size_t)args->rho)*2+1;
+        kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
 
         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
                               kernel->height*sizeof(double));
@@ -1375,9 +1499,9 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
           return(DestroyKernelInfo(kernel));
 
         /* set all kernel values within diamond area to scale given */
-        for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-          for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
-            if ((labs(u)+labs(v)) <= (long)kernel->x)
+        for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+          for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
+            if ( (labs((long) u)+labs((long) v)) <= (long) kernel->x)
               kernel->positive_range += kernel->values[i] = args->sigma;
             else
               kernel->values[i] = nan;
@@ -1393,21 +1517,21 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
             if (args->rho < 1.0)
               kernel->width = kernel->height = 3;  /* default radius = 1 */
             else
-              kernel->width = kernel->height = (unsigned long) (2*args->rho+1);
-            kernel->x = kernel->y = (long) (kernel->width-1)/2;
+              kernel->width = kernel->height = (size_t) (2*args->rho+1);
+            kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
             scale = args->sigma;
           }
         else {
             /* NOTE: user defaults set in "AcquireKernelInfo()" */
             if ( args->rho < 1.0 || args->sigma < 1.0 )
               return(DestroyKernelInfo(kernel));    /* invalid args given */
-            kernel->width = (unsigned long)args->rho;
-            kernel->height = (unsigned long)args->sigma;
+            kernel->width = (size_t)args->rho;
+            kernel->height = (size_t)args->sigma;
             if ( args->xi  < 0.0 || args->xi  > (double)kernel->width ||
                  args->psi < 0.0 || args->psi > (double)kernel->height )
               return(DestroyKernelInfo(kernel));    /* invalid args given */
-            kernel->x = (long) args->xi;
-            kernel->y = (long) args->psi;
+            kernel->x = (ssize_t) args->xi;
+            kernel->y = (ssize_t) args->psi;
             scale = 1.0;
           }
         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
@@ -1416,384 +1540,589 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
           return(DestroyKernelInfo(kernel));
 
         /* set all kernel values to scale given */
-        u=(long) kernel->width*kernel->height;
+        u=(ssize_t) (kernel->width*kernel->height);
         for ( i=0; i < u; i++)
             kernel->values[i] = scale;
         kernel->minimum = kernel->maximum = scale;   /* a flat shape */
         kernel->positive_range = scale*u;
         break;
       }
-    case DiskKernel:
-      {
-        long limit = (long)(args->rho*args->rho);
-        if (args->rho < 0.1)             /* default radius approx 3.5 */
-          kernel->width = kernel->height = 7L, limit = 10L;
-        else
-           kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
-        kernel->x = kernel->y = (long) (kernel->width-1)/2;
-
-        kernel->values=(double *) AcquireQuantumMemory(kernel->width,
-                              kernel->height*sizeof(double));
-        if (kernel->values == (double *) NULL)
-          return(DestroyKernelInfo(kernel));
+      case OctagonKernel:
+        {
+          if (args->rho < 1.0)
+            kernel->width = kernel->height = 5;  /* default radius = 2 */
+          else
+            kernel->width = kernel->height = ((size_t)args->rho)*2+1;
+          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
+
+          kernel->values=(double *) AcquireQuantumMemory(kernel->width,
+                                kernel->height*sizeof(double));
+          if (kernel->values == (double *) NULL)
+            return(DestroyKernelInfo(kernel));
+
+          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
+              if ( (labs((long) u)+labs((long) v)) <=
+                        ((long)kernel->x + (long)(kernel->x/2)) )
+                kernel->positive_range += kernel->values[i] = args->sigma;
+              else
+                kernel->values[i] = nan;
+          kernel->minimum = kernel->maximum = args->sigma;  /* a flat shape */
+          break;
+        }
+      case DiskKernel:
+        {
+          ssize_t
+            limit = (ssize_t)(args->rho*args->rho);
 
-        /* set all kernel values within disk area to scale given */
-        for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-          for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
-            if ((u*u+v*v) <= limit)
-              kernel->positive_range += kernel->values[i] = args->sigma;
-            else
-              kernel->values[i] = nan;
-        kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
-        break;
-      }
-    case PlusKernel:
+          if (args->rho < 0.4)           /* default radius approx 4.3 */
+            kernel->width = kernel->height = 9L, limit = 18L;
+          else
+            kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
+          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
+
+          kernel->values=(double *) AcquireQuantumMemory(kernel->width,
+                                kernel->height*sizeof(double));
+          if (kernel->values == (double *) NULL)
+            return(DestroyKernelInfo(kernel));
+
+          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
+              if ((u*u+v*v) <= limit)
+                kernel->positive_range += kernel->values[i] = args->sigma;
+              else
+                kernel->values[i] = nan;
+          kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
+          break;
+        }
+      case PlusKernel:
+        {
+          if (args->rho < 1.0)
+            kernel->width = kernel->height = 5;  /* default radius 2 */
+          else
+            kernel->width = kernel->height = ((size_t)args->rho)*2+1;
+          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
+
+          kernel->values=(double *) AcquireQuantumMemory(kernel->width,
+                                kernel->height*sizeof(double));
+          if (kernel->values == (double *) NULL)
+            return(DestroyKernelInfo(kernel));
+
+          /* set all kernel values along axises to given scale */
+          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
+              kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
+          kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
+          kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
+          break;
+        }
+      case CrossKernel:
+        {
+          if (args->rho < 1.0)
+            kernel->width = kernel->height = 5;  /* default radius 2 */
+          else
+            kernel->width = kernel->height = ((size_t)args->rho)*2+1;
+          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
+
+          kernel->values=(double *) AcquireQuantumMemory(kernel->width,
+                                kernel->height*sizeof(double));
+          if (kernel->values == (double *) NULL)
+            return(DestroyKernelInfo(kernel));
+
+          /* set all kernel values along axises to given scale */
+          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
+              kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
+          kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
+          kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
+          break;
+        }
+      /*
+        HitAndMiss Kernels
+      */
+      case RingKernel:
+      case PeaksKernel:
+        {
+          ssize_t
+            limit1,
+            limit2,
+            scale;
+
+          if (args->rho < args->sigma)
+            {
+              kernel->width = ((size_t)args->sigma)*2+1;
+              limit1 = (ssize_t)(args->rho*args->rho);
+              limit2 = (ssize_t)(args->sigma*args->sigma);
+            }
+          else
+            {
+              kernel->width = ((size_t)args->rho)*2+1;
+              limit1 = (ssize_t)(args->sigma*args->sigma);
+              limit2 = (ssize_t)(args->rho*args->rho);
+            }
+          if ( limit2 <= 0 )
+            kernel->width = 7L, limit1 = 7L, limit2 = 11L;
+
+          kernel->height = kernel->width;
+          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
+          kernel->values=(double *) AcquireQuantumMemory(kernel->width,
+                                kernel->height*sizeof(double));
+          if (kernel->values == (double *) NULL)
+            return(DestroyKernelInfo(kernel));
+
+          /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
+          scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
+          for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
+            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
+              { ssize_t radius=u*u+v*v;
+                if (limit1 < radius && radius <= limit2)
+                  kernel->positive_range += kernel->values[i] = (double) scale;
+                else
+                  kernel->values[i] = nan;
+              }
+          kernel->minimum = kernel->maximum = (double) scale;
+          if ( type == PeaksKernel ) {
+            /* set the central point in the middle */
+            kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
+            kernel->positive_range = 1.0;
+            kernel->maximum = 1.0;
+          }
+          break;
+        }
+      case EdgesKernel:
+        {
+          kernel=AcquireKernelInfo("ThinSE:482");
+          if (kernel == (KernelInfo *) NULL)
+            return(kernel);
+          kernel->type = type;
+          ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
+          break;
+        }
+      case CornersKernel:
+        {
+          kernel=AcquireKernelInfo("ThinSE:87");
+          if (kernel == (KernelInfo *) NULL)
+            return(kernel);
+          kernel->type = type;
+          ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
+          break;
+        }
+      case DiagonalsKernel:
+        {
+          switch ( (int) args->rho ) {
+            case 0:
+            default:
+              { KernelInfo
+                  *new_kernel;
+                kernel=ParseKernelArray("3: 0,0,0  0,-,1  1,1,-");
+                if (kernel == (KernelInfo *) NULL)
+                  return(kernel);
+                kernel->type = type;
+                new_kernel=ParseKernelArray("3: 0,0,1  0,-,1  0,1,-");
+                if (new_kernel == (KernelInfo *) NULL)
+                  return(DestroyKernelInfo(kernel));
+                new_kernel->type = type;
+                LastKernelInfo(kernel)->next = new_kernel;
+                ExpandMirrorKernelInfo(kernel);
+                return(kernel);
+              }
+            case 1:
+              kernel=ParseKernelArray("3: 0,0,0  0,-,1  1,1,-");
+              break;
+            case 2:
+              kernel=ParseKernelArray("3: 0,0,1  0,-,1  0,1,-");
+              break;
+          }
+          if (kernel == (KernelInfo *) NULL)
+            return(kernel);
+          kernel->type = type;
+          RotateKernelInfo(kernel, args->sigma);
+          break;
+        }
+      case LineEndsKernel:
+        { /* Kernels for finding the end of thin lines */
+          switch ( (int) args->rho ) {
+            case 0:
+            default:
+              /* set of kernels to find all end of lines */
+              return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>"));
+            case 1:
+              /* kernel for 4-connected line ends - no rotation */
+              kernel=ParseKernelArray("3: 0,0,-  0,1,1  0,0,-");
+              break;
+          case 2:
+              /* kernel to add for 8-connected lines - no rotation */
+              kernel=ParseKernelArray("3: 0,0,0  0,1,0  0,0,1");
+              break;
+          case 3:
+              /* kernel to add for orthogonal line ends - does not find corners */
+              kernel=ParseKernelArray("3: 0,0,0  0,1,1  0,0,0");
+              break;
+          case 4:
+              /* traditional line end - fails on last T end */
+              kernel=ParseKernelArray("3: 0,0,0  0,1,-  0,0,-");
+              break;
+          }
+          if (kernel == (KernelInfo *) NULL)
+            return(kernel);
+          kernel->type = type;
+          RotateKernelInfo(kernel, args->sigma);
+          break;
+        }
+      case LineJunctionsKernel:
+        { /* kernels for finding the junctions of multiple lines */
+          switch ( (int) args->rho ) {
+            case 0:
+            default:
+              /* set of kernels to find all line junctions */
+              return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>"));
+            case 1:
+              /* Y Junction */
+              kernel=ParseKernelArray("3: 1,-,1  -,1,-  -,1,-");
+              break;
+            case 2:
+              /* Diagonal T Junctions */
+              kernel=ParseKernelArray("3: 1,-,-  -,1,-  1,-,1");
+              break;
+            case 3:
+              /* Orthogonal T Junctions */
+              kernel=ParseKernelArray("3: -,-,-  1,1,1  -,1,-");
+              break;
+            case 4:
+              /* Diagonal X Junctions */
+              kernel=ParseKernelArray("3: 1,-,1  -,1,-  1,-,1");
+              break;
+            case 5:
+              /* Orthogonal X Junctions - minimal diamond kernel */
+              kernel=ParseKernelArray("3: -,1,-  1,1,1  -,1,-");
+              break;
+          }
+          if (kernel == (KernelInfo *) NULL)
+            return(kernel);
+          kernel->type = type;
+          RotateKernelInfo(kernel, args->sigma);
+          break;
+        }
+      case RidgesKernel:
+        { /* Ridges - Ridge finding kernels */
+          KernelInfo
+            *new_kernel;
+          switch ( (int) args->rho ) {
+            case 1:
+            default:
+              kernel=ParseKernelArray("3x1:0,1,0");
+              if (kernel == (KernelInfo *) NULL)
+                return(kernel);
+              kernel->type = type;
+              ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
+              break;
+            case 2:
+              kernel=ParseKernelArray("4x1:0,1,1,0");
+              if (kernel == (KernelInfo *) NULL)
+                return(kernel);
+              kernel->type = type;
+              ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */
+
+              /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
+              /* Unfortunatally we can not yet rotate a non-square kernel */
+              /* But then we can't flip a non-symetrical kernel either */
+              new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
+              if (new_kernel == (KernelInfo *) NULL)
+                return(DestroyKernelInfo(kernel));
+              new_kernel->type = type;
+              LastKernelInfo(kernel)->next = new_kernel;
+              new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
+              if (new_kernel == (KernelInfo *) NULL)
+                return(DestroyKernelInfo(kernel));
+              new_kernel->type = type;
+              LastKernelInfo(kernel)->next = new_kernel;
+              new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
+              if (new_kernel == (KernelInfo *) NULL)
+                return(DestroyKernelInfo(kernel));
+              new_kernel->type = type;
+              LastKernelInfo(kernel)->next = new_kernel;
+              new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
+              if (new_kernel == (KernelInfo *) NULL)
+                return(DestroyKernelInfo(kernel));
+              new_kernel->type = type;
+              LastKernelInfo(kernel)->next = new_kernel;
+              new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
+              if (new_kernel == (KernelInfo *) NULL)
+                return(DestroyKernelInfo(kernel));
+              new_kernel->type = type;
+              LastKernelInfo(kernel)->next = new_kernel;
+              new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
+              if (new_kernel == (KernelInfo *) NULL)
+                return(DestroyKernelInfo(kernel));
+              new_kernel->type = type;
+              LastKernelInfo(kernel)->next = new_kernel;
+              new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
+              if (new_kernel == (KernelInfo *) NULL)
+                return(DestroyKernelInfo(kernel));
+              new_kernel->type = type;
+              LastKernelInfo(kernel)->next = new_kernel;
+              new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
+              if (new_kernel == (KernelInfo *) NULL)
+                return(DestroyKernelInfo(kernel));
+              new_kernel->type = type;
+              LastKernelInfo(kernel)->next = new_kernel;
+              break;
+          }
+          break;
+        }
+      case ConvexHullKernel:
+        {
+          KernelInfo
+            *new_kernel;
+          /* first set of 8 kernels */
+          kernel=ParseKernelArray("3: 1,1,-  1,0,-  1,-,0");
+          if (kernel == (KernelInfo *) NULL)
+            return(kernel);
+          kernel->type = type;
+          ExpandRotateKernelInfo(kernel, 90.0);
+          /* append the mirror versions too - no flip function yet */
+          new_kernel=ParseKernelArray("3: 1,1,1  1,0,-  -,-,0");
+          if (new_kernel == (KernelInfo *) NULL)
+            return(DestroyKernelInfo(kernel));
+          new_kernel->type = type;
+          ExpandRotateKernelInfo(new_kernel, 90.0);
+          LastKernelInfo(kernel)->next = new_kernel;
+          break;
+        }
+      case SkeletonKernel:
+        {
+          switch ( (int) args->rho ) {
+            case 1:
+            default:
+              /* Traditional Skeleton...
+              ** A cyclically rotated single kernel
+              */
+              kernel=AcquireKernelInfo("ThinSE:482");
+              if (kernel == (KernelInfo *) NULL)
+                return(kernel);
+              kernel->type = type;
+              ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
+              break;
+            case 2:
+              /* HIPR Variation of the cyclic skeleton
+              ** Corners of the traditional method made more forgiving,
+              ** but the retain the same cyclic order.
+              */
+              kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;");
+              if (kernel == (KernelInfo *) NULL)
+                return(kernel);
+              if (kernel->next == (KernelInfo *) NULL)
+                return(DestroyKernelInfo(kernel));
+              kernel->type = type;
+              kernel->next->type = type;
+              ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
+              break;
+            case 3:
+              /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
+              ** "Connectivity-Preserving Morphological Image Thransformations"
+              ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
+              **   http://www.leptonica.com/papers/conn.pdf
+              */
+              kernel=AcquireKernelInfo(
+                            "ThinSE:41; ThinSE:42; ThinSE:43");
+              if (kernel == (KernelInfo *) NULL)
+                return(kernel);
+              kernel->type = type;
+              kernel->next->type = type;
+              kernel->next->next->type = type;
+              ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
+              break;
+           }
+          break;
+        }
+      case ThinSEKernel:
+        { /* Special kernels for general thinning, while preserving connections
+          ** "Connectivity-Preserving Morphological Image Thransformations"
+          ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
+          **   http://www.leptonica.com/papers/conn.pdf
+          ** And
+          **   http://tpgit.github.com/Leptonica/ccthin_8c_source.html
+          **
+          ** Note kernels do not specify the origin pixel, allowing them
+          ** to be used for both thickening and thinning operations.
+          */
+          switch ( (int) args->rho ) {
+            /* SE for 4-connected thinning */
+            case 41: /* SE_4_1 */
+              kernel=ParseKernelArray("3: -,-,1  0,-,1  -,-,1");
+              break;
+            case 42: /* SE_4_2 */
+              kernel=ParseKernelArray("3: -,-,1  0,-,1  -,0,-");
+              break;
+            case 43: /* SE_4_3 */
+              kernel=ParseKernelArray("3: -,0,-  0,-,1  -,-,1");
+              break;
+            case 44: /* SE_4_4 */
+              kernel=ParseKernelArray("3: -,0,-  0,-,1  -,0,-");
+              break;
+            case 45: /* SE_4_5 */
+              kernel=ParseKernelArray("3: -,0,1  0,-,1  -,0,-");
+              break;
+            case 46: /* SE_4_6 */
+              kernel=ParseKernelArray("3: -,0,-  0,-,1  -,0,1");
+              break;
+            case 47: /* SE_4_7 */
+              kernel=ParseKernelArray("3: -,1,1  0,-,1  -,0,-");
+              break;
+            case 48: /* SE_4_8 */
+              kernel=ParseKernelArray("3: -,-,1  0,-,1  0,-,1");
+              break;
+            case 49: /* SE_4_9 */
+              kernel=ParseKernelArray("3: 0,-,1  0,-,1  -,-,1");
+              break;
+            /* SE for 8-connected thinning - negatives of the above */
+            case 81: /* SE_8_0 */
+              kernel=ParseKernelArray("3: -,1,-  0,-,1  -,1,-");
+              break;
+            case 82: /* SE_8_2 */
+              kernel=ParseKernelArray("3: -,1,-  0,-,1  0,-,-");
+              break;
+            case 83: /* SE_8_3 */
+              kernel=ParseKernelArray("3: 0,-,-  0,-,1  -,1,-");
+              break;
+            case 84: /* SE_8_4 */
+              kernel=ParseKernelArray("3: 0,-,-  0,-,1  0,-,-");
+              break;
+            case 85: /* SE_8_5 */
+              kernel=ParseKernelArray("3: 0,-,1  0,-,1  0,-,-");
+              break;
+            case 86: /* SE_8_6 */
+              kernel=ParseKernelArray("3: 0,-,-  0,-,1  0,-,1");
+              break;
+            case 87: /* SE_8_7 */
+              kernel=ParseKernelArray("3: -,1,-  0,-,1  0,0,-");
+              break;
+            case 88: /* SE_8_8 */
+              kernel=ParseKernelArray("3: -,1,-  0,-,1  0,1,-");
+              break;
+            case 89: /* SE_8_9 */
+              kernel=ParseKernelArray("3: 0,1,-  0,-,1  -,1,-");
+              break;
+            /* Special combined SE kernels */
+            case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
+              kernel=ParseKernelArray("3: -,-,1  0,-,-  -,0,-");
+              break;
+            case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
+              kernel=ParseKernelArray("3: -,1,-  -,-,1  0,-,-");
+              break;
+            case 481: /* SE_48_1 - General Connected Corner Kernel */
+              kernel=ParseKernelArray("3: -,1,1  0,-,1  0,0,-");
+              break;
+            default:
+            case 482: /* SE_48_2 - General Edge Kernel */
+              kernel=ParseKernelArray("3: 0,-,1  0,-,1  0,-,1");
+              break;
+          }
+          if (kernel == (KernelInfo *) NULL)
+            return(kernel);
+          kernel->type = type;
+          RotateKernelInfo(kernel, args->sigma);
+          break;
+        }
+      /*
+        Distance Measuring Kernels
+      */
+      case ChebyshevKernel:
+        {
+          if (args->rho < 1.0)
+            kernel->width = kernel->height = 3;  /* default radius = 1 */
+          else
+            kernel->width = kernel->height = ((size_t)args->rho)*2+1;
+          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
+
+          kernel->values=(double *) AcquireQuantumMemory(kernel->width,
+                                kernel->height*sizeof(double));
+          if (kernel->values == (double *) NULL)
+            return(DestroyKernelInfo(kernel));
+
+          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
+              kernel->positive_range += ( kernel->values[i] =
+                  args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
+          kernel->maximum = kernel->values[0];
+          break;
+        }
+      case ManhattanKernel:
+        {
+          if (args->rho < 1.0)
+            kernel->width = kernel->height = 3;  /* default radius = 1 */
+          else
+            kernel->width = kernel->height = ((size_t)args->rho)*2+1;
+          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
+
+          kernel->values=(double *) AcquireQuantumMemory(kernel->width,
+                                kernel->height*sizeof(double));
+          if (kernel->values == (double *) NULL)
+            return(DestroyKernelInfo(kernel));
+
+          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
+              kernel->positive_range += ( kernel->values[i] =
+                  args->sigma*(labs((long) u)+labs((long) v)) );
+          kernel->maximum = kernel->values[0];
+          break;
+        }
+      case OctagonalKernel:
       {
-        if (args->rho < 1.0)
-          kernel->width = kernel->height = 5;  /* default radius 2 */
+        if (args->rho < 2.0)
+          kernel->width = kernel->height = 5;  /* default/minimum radius = 2 */
         else
-           kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
-        kernel->x = kernel->y = (long) (kernel->width-1)/2;
+          kernel->width = kernel->height = ((size_t)args->rho)*2+1;
+        kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
 
         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
                               kernel->height*sizeof(double));
         if (kernel->values == (double *) NULL)
           return(DestroyKernelInfo(kernel));
 
-        /* set all kernel values along axises to given scale */
-        for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-          for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
-            kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
-        kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
-        kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
+        for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+          for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
+            {
+              double
+                r1 = MagickMax(fabs((double)u),fabs((double)v)),
+                r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
+              kernel->positive_range += kernel->values[i] =
+                        args->sigma*MagickMax(r1,r2);
+            }
+        kernel->maximum = kernel->values[0];
         break;
       }
-    case CrossKernel:
+    case EuclideanKernel:
       {
         if (args->rho < 1.0)
-          kernel->width = kernel->height = 5;  /* default radius 2 */
-        else
-           kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
-        kernel->x = kernel->y = (long) (kernel->width-1)/2;
-
-        kernel->values=(double *) AcquireQuantumMemory(kernel->width,
-                              kernel->height*sizeof(double));
-        if (kernel->values == (double *) NULL)
-          return(DestroyKernelInfo(kernel));
-
-        /* set all kernel values along axises to given scale */
-        for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-          for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
-            kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
-        kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
-        kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
-        break;
-      }
-    /* HitAndMiss Kernels */
-    case RingKernel:
-    case PeaksKernel:
-      {
-        long
-          limit1,
-          limit2,
-          scale;
-
-        if (args->rho < args->sigma)
-          {
-            kernel->width = ((unsigned long)args->sigma)*2+1;
-            limit1 = (long)args->rho*args->rho;
-            limit2 = (long)args->sigma*args->sigma;
-          }
+          kernel->width = kernel->height = 3;  /* default radius = 1 */
         else
-          {
-            kernel->width = ((unsigned long)args->rho)*2+1;
-            limit1 = (long)args->sigma*args->sigma;
-            limit2 = (long)args->rho*args->rho;
-          }
-        if ( limit2 <= 0 )
-          kernel->width = 7L, limit1 = 7L, limit2 = 11L;
+          kernel->width = kernel->height = ((size_t)args->rho)*2+1;
+        kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
 
-        kernel->height = kernel->width;
-        kernel->x = kernel->y = (long) (kernel->width-1)/2;
         kernel->values=(double *) AcquireQuantumMemory(kernel->width,
                               kernel->height*sizeof(double));
         if (kernel->values == (double *) NULL)
           return(DestroyKernelInfo(kernel));
 
-        /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
-        scale = (long) (( type == PeaksKernel) ? 0.0 : args->xi);
-        for ( i=0, v= -kernel->y; v <= (long)kernel->y; v++)
-          for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
-            { long radius=u*u+v*v;
-              if (limit1 < radius && radius <= limit2)
-                kernel->positive_range += kernel->values[i] = (double) scale;
-              else
-                kernel->values[i] = nan;
-            }
-        kernel->minimum = kernel->minimum = (double) scale;
-        if ( type == PeaksKernel ) {
-          /* set the central point in the middle */
-          kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
-          kernel->positive_range = 1.0;
-          kernel->maximum = 1.0;
-        }
-        break;
-      }
-    case EdgesKernel:
-      {
-        kernel=ParseKernelArray("3: 0,0,0  -,1,-  1,1,1");
-        if (kernel == (KernelInfo *) NULL)
-          return(kernel);
-        kernel->type = type;
-        ExpandKernelInfo(kernel, 90.0); /* Create a list of 4 rotated kernels */
+        for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
+          for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
+            kernel->positive_range += ( kernel->values[i] =
+                 args->sigma*sqrt((double)(u*u+v*v)) );
+        kernel->maximum = kernel->values[0];
         break;
       }
-    case CornersKernel:
+    default:
       {
-        kernel=ParseKernelArray("3: 0,0,-  0,1,1  -,1,-");
+        /* No-Op Kernel - Basically just a single pixel on its own */
+        kernel=ParseKernelArray("1:1");
         if (kernel == (KernelInfo *) NULL)
           return(kernel);
-        kernel->type = type;
-        ExpandKernelInfo(kernel, 90.0); /* Create a list of 4 rotated kernels */
-        break;
-      }
-    case RidgesKernel:
-      {
-        kernel=ParseKernelArray("3: 0,1,0 ");
-        if (kernel == (KernelInfo *) NULL)
-          return(kernel);
-        kernel->type = type;
-        ExpandKernelInfo(kernel, 45.0); /* 4 rotated kernels (symmetrical) */
-        break;
-      }
-    case Ridges2Kernel:
-      {
-        KernelInfo
-          *new_kernel;
-        kernel=ParseKernelArray("4x1^:0,1,1,0");
-        if (kernel == (KernelInfo *) NULL)
-          return(kernel);
-        kernel->type = type;
-        ExpandKernelInfo(kernel, 90.0); /* 4 rotated kernels */
-#if 0
-        /* 2 pixel diagonaly thick - 4 rotates - not needed? */
-        new_kernel=ParseKernelArray("4x4^:0,-,-,- -,1,-,- -,-,1,- -,-,-,0'");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        ExpandKernelInfo(new_kernel, 90.0);  /* 4 rotated kernels */
-        LastKernelInfo(kernel)->next = new_kernel;
-#endif
-        /* kernels to find a stepped 'thick' line - 4 rotates * mirror */
-        /* Unfortunatally we can not yet rotate a non-square kernel */
-        /* But then we can't flip a non-symetrical kernel either */
-        new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        LastKernelInfo(kernel)->next = new_kernel;
-        new_kernel=ParseKernelArray("4x3+2+1^:0,1,1,- -,1,1,- -,1,1,0");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        LastKernelInfo(kernel)->next = new_kernel;
-        new_kernel=ParseKernelArray("4x3+1+1^:-,1,1,0 -,1,1,- 0,1,1,-");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        LastKernelInfo(kernel)->next = new_kernel;
-        new_kernel=ParseKernelArray("4x3+2+1^:-,1,1,0 -,1,1,- 0,1,1,-");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        LastKernelInfo(kernel)->next = new_kernel;
-        new_kernel=ParseKernelArray("3x4+1+1^:0,-,- 1,1,1 1,1,1 -,-,0");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        LastKernelInfo(kernel)->next = new_kernel;
-        new_kernel=ParseKernelArray("3x4+1+2^:0,-,- 1,1,1 1,1,1 -,-,0");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        LastKernelInfo(kernel)->next = new_kernel;
-        new_kernel=ParseKernelArray("3x4+1+1^:-,-,0 1,1,1 1,1,1 0,-,-");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        LastKernelInfo(kernel)->next = new_kernel;
-        new_kernel=ParseKernelArray("3x4+1+2^:-,-,0 1,1,1 1,1,1 0,-,-");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        LastKernelInfo(kernel)->next = new_kernel;
-        break;
-      }
-    case LineEndsKernel:
-      {
-        KernelInfo
-          *new_kernel;
-        kernel=ParseKernelArray("3: 0,0,0  0,1,0  -,1,-");
-        if (kernel == (KernelInfo *) NULL)
-          return(kernel);
-        kernel->type = type;
-        ExpandKernelInfo(kernel, 90.0);
-        /* append second set of 4 kernels */
-        new_kernel=ParseKernelArray("3: 0,0,0  0,1,0  0,0,1");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        ExpandKernelInfo(new_kernel, 90.0);
-        LastKernelInfo(kernel)->next = new_kernel;
-        break;
-      }
-    case LineJunctionsKernel:
-      {
-        KernelInfo
-          *new_kernel;
-        /* first set of 4 kernels */
-        kernel=ParseKernelArray("3: -,1,-  -,1,-  1,-,1");
-        if (kernel == (KernelInfo *) NULL)
-          return(kernel);
-        kernel->type = type;
-        ExpandKernelInfo(kernel, 45.0);
-        /* append second set of 4 kernels */
-        new_kernel=ParseKernelArray("3: 1,-,-  -,1,-  1,-,1");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        ExpandKernelInfo(new_kernel, 90.0);
-        LastKernelInfo(kernel)->next = new_kernel;
-        break;
-      }
-    case ConvexHullKernel:
-      {
-        KernelInfo
-          *new_kernel;
-        /* first set of 4 kernels */
-        kernel=ParseKernelArray("3: 1,1,-  1,0,-  1,-,0");
-        if (kernel == (KernelInfo *) NULL)
-          return(kernel);
-        kernel->type = type;
-        ExpandKernelInfo(kernel, 90.0);
-        /* append second set of 4 kernels */
-        new_kernel=ParseKernelArray("3: 1,1,1  1,0,0  -,-,0");
-        if (new_kernel == (KernelInfo *) NULL)
-          return(DestroyKernelInfo(kernel));
-        new_kernel->type = type;
-        ExpandKernelInfo(new_kernel, 90.0);
-        LastKernelInfo(kernel)->next = new_kernel;
-        break;
-      }
-    case SkeletonKernel:
-      { /* what is the best form for skeletonization by thinning? */
-#if 0
-#  if 0
-        kernel=AcquireKernelInfo("Corners;Edges");
-#  else
-        kernel=AcquireKernelInfo("Edges;Corners");
-#  endif
-#else
-        kernel=ParseKernelArray("3: 0,0,-  0,1,1  -,1,1");
-        if (kernel == (KernelInfo *) NULL)
-          return(kernel);
-        kernel->type = type;
-        ExpandKernelInfo(kernel, 45);
-        break;
-#endif
-        break;
-      }
-    case MatKernel: /* experimental - MAT from a Distance Gradient */
-      {
-        KernelInfo
-          *new_kernel;
-        /* Ridge Kernel but without the diagonal */
-        kernel=ParseKernelArray("3x1: 0,1,0");
-        if (kernel == (KernelInfo *) NULL)
-          return(kernel);
-        kernel->type = RidgesKernel;
-        ExpandKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
-        /* Plus the 2 pixel ridges kernel - no diagonal */
-        new_kernel=AcquireKernelBuiltIn(Ridges2Kernel,args);
-        if (new_kernel == (KernelInfo *) NULL)
-          return(kernel);
-        LastKernelInfo(kernel)->next = new_kernel;
-        break;
-      }
-    /* Distance Measuring Kernels */
-    case ChebyshevKernel:
-      {
-        if (args->rho < 1.0)
-          kernel->width = kernel->height = 3;  /* default radius = 1 */
-        else
-          kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
-        kernel->x = kernel->y = (long) (kernel->width-1)/2;
-
-        kernel->values=(double *) AcquireQuantumMemory(kernel->width,
-                              kernel->height*sizeof(double));
-        if (kernel->values == (double *) NULL)
-          return(DestroyKernelInfo(kernel));
-
-        for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-          for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
-            kernel->positive_range += ( kernel->values[i] =
-                 args->sigma*((labs(u)>labs(v)) ? labs(u) : labs(v)) );
-        kernel->maximum = kernel->values[0];
-        break;
-      }
-    case ManhattenKernel:
-      {
-        if (args->rho < 1.0)
-          kernel->width = kernel->height = 3;  /* default radius = 1 */
-        else
-           kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
-        kernel->x = kernel->y = (long) (kernel->width-1)/2;
-
-        kernel->values=(double *) AcquireQuantumMemory(kernel->width,
-                              kernel->height*sizeof(double));
-        if (kernel->values == (double *) NULL)
-          return(DestroyKernelInfo(kernel));
-
-        for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-          for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
-            kernel->positive_range += ( kernel->values[i] =
-                 args->sigma*(labs(u)+labs(v)) );
-        kernel->maximum = kernel->values[0];
-        break;
-      }
-    case EuclideanKernel:
-      {
-        if (args->rho < 1.0)
-          kernel->width = kernel->height = 3;  /* default radius = 1 */
-        else
-           kernel->width = kernel->height = ((unsigned long)args->rho)*2+1;
-        kernel->x = kernel->y = (long) (kernel->width-1)/2;
-
-        kernel->values=(double *) AcquireQuantumMemory(kernel->width,
-                              kernel->height*sizeof(double));
-        if (kernel->values == (double *) NULL)
-          return(DestroyKernelInfo(kernel));
-
-        for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-          for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
-            kernel->positive_range += ( kernel->values[i] =
-                 args->sigma*sqrt((double)(u*u+v*v)) );
-        kernel->maximum = kernel->values[0];
-        break;
-      }
-    case UnityKernel:
-    default:
-      {
-        /* Unity or No-Op Kernel - 3x3 with 1 in center */
-        kernel=ParseKernelArray("3:0,0,0,0,1,0,0,0,0");
-        if (kernel == (KernelInfo *) NULL)
-          return(kernel);
-        kernel->type = ( type == UnityKernel ) ? UnityKernel : UndefinedKernel;
+        kernel->type = UndefinedKernel;
         break;
       }
       break;
   }
-
   return(kernel);
 }
 \f
@@ -1823,7 +2152,7 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
 */
 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
 {
-  register long
+  register ssize_t
     i;
 
   KernelInfo
@@ -1840,7 +2169,7 @@ MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
     kernel->height*sizeof(double));
   if (new_kernel->values == (double *) NULL)
     return(DestroyKernelInfo(new_kernel));
-  for (i=0; i < (long) (kernel->width*kernel->height); i++)
+  for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
     new_kernel->values[i]=kernel->values[i];
 
   /* Also clone the next kernel in the kernel list */
@@ -1876,7 +2205,6 @@ MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
 %    o kernel: the Morphology/Convolution kernel to be destroyed
 %
 */
-
 MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
 {
   assert(kernel != (KernelInfo *) NULL);
@@ -1894,21 +2222,96 @@ MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%     E x p a n d K e r n e l I n f o                                         %
++     E x p a n d M i r r o r K e r n e l I n f o                             %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  ExpandKernelInfo() takes a single kernel, and expands it into a list
-%  of kernels each incrementally rotated the angle given.
+%  ExpandMirrorKernelInfo() takes a single kernel, and expands it into a
+%  sequence of 90-degree rotated kernels but providing a reflected 180
+%  rotatation, before the -/+ 90-degree rotations.
+%
+%  This special rotation order produces a better, more symetrical thinning of
+%  objects.
+%
+%  The format of the ExpandMirrorKernelInfo method is:
+%
+%      void ExpandMirrorKernelInfo(KernelInfo *kernel)
+%
+%  A description of each parameter follows:
+%
+%    o kernel: the Morphology/Convolution kernel
+%
+% This function is only internel to this module, as it is not finalized,
+% especially with regard to non-orthogonal angles, and rotation of larger
+% 2D kernels.
+*/
+
+#if 0
+static void FlopKernelInfo(KernelInfo *kernel)
+    { /* Do a Flop by reversing each row. */
+      size_t
+        y;
+      register ssize_t
+        x,r;
+      register double
+        *k,t;
+
+      for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
+        for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
+          t=k[x],  k[x]=k[r],  k[r]=t;
+
+      kernel->x = kernel->width - kernel->x - 1;
+      angle = fmod(angle+180.0, 360.0);
+    }
+#endif
+
+static void ExpandMirrorKernelInfo(KernelInfo *kernel)
+{
+  KernelInfo
+    *clone,
+    *last;
+
+  last = kernel;
+
+  clone = CloneKernelInfo(last);
+  RotateKernelInfo(clone, 180);   /* flip */
+  LastKernelInfo(last)->next = clone;
+  last = clone;
+
+  clone = CloneKernelInfo(last);
+  RotateKernelInfo(clone, 90);   /* transpose */
+  LastKernelInfo(last)->next = clone;
+  last = clone;
+
+  clone = CloneKernelInfo(last);
+  RotateKernelInfo(clone, 180);  /* flop */
+  LastKernelInfo(last)->next = clone;
+
+  return;
+}
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
++     E x p a n d R o t a t e K e r n e l I n f o                             %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
+%  incrementally by the angle given, until the kernel repeats.
 %
 %  WARNING: 45 degree rotations only works for 3x3 kernels.
 %  While 90 degree roatations only works for linear and square kernels
 %
-%  The format of the RotateKernelInfo method is:
+%  The format of the ExpandRotateKernelInfo method is:
 %
-%      void ExpandKernelInfo(KernelInfo *kernel, double angle)
+%      void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
 %
 %  A description of each parameter follows:
 %
@@ -1925,7 +2328,7 @@ MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
 static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
      const KernelInfo *kernel2)
 {
-  register unsigned long
+  register size_t
     i;
 
   /* check size and origin location */
@@ -1937,12 +2340,12 @@ static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
 
   /* check actual kernel values */
   for (i=0; i < (kernel1->width*kernel1->height); i++) {
-    /* Test for Nan equivelence */
+    /* Test for Nan equivalence */
     if ( IsNan(kernel1->values[i]) && !IsNan(kernel2->values[i]) )
       return MagickFalse;
     if ( IsNan(kernel2->values[i]) && !IsNan(kernel1->values[i]) )
       return MagickFalse;
-    /* Test actual values are equivelent */
+    /* Test actual values are equivalent */
     if ( fabs(kernel1->values[i] - kernel2->values[i]) > MagickEpsilon )
       return MagickFalse;
   }
@@ -1950,7 +2353,7 @@ static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
   return MagickTrue;
 }
 
-static void ExpandKernelInfo(KernelInfo *kernel, const double angle)
+static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
 {
   KernelInfo
     *clone,
@@ -1962,10 +2365,10 @@ static void ExpandKernelInfo(KernelInfo *kernel, const double angle)
     RotateKernelInfo(clone, angle);
     if ( SameKernelInfo(kernel, clone) == MagickTrue )
       break;
-    last->next = clone;
+    LastKernelInfo(last)->next = clone;
     last = clone;
   }
-  clone = DestroyKernelInfo(clone); /* This was the same as the first - junk */
+  clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
   return;
 }
 \f
@@ -2007,7 +2410,7 @@ static void ExpandKernelInfo(KernelInfo *kernel, const double angle)
 */
 static void CalcKernelMetaData(KernelInfo *kernel)
 {
-  register unsigned long
+  register size_t
     i;
 
   kernel->minimum = kernel->maximum = 0.0;
@@ -2040,28 +2443,36 @@ static void CalcKernelMetaData(KernelInfo *kernel)
 %  MorphologyApply() applies a morphological method, multiple times using
 %  a list of multiple kernels.
 %
-%  It is basically equivelent to as MorphologyImageChannel() (see below) but
-%  without user controls, that that function extracts and applies to kernels
-%  and morphology methods.
+%  It is basically equivalent to as MorphologyImageChannel() (see below) but
+%  without any user controls.  This allows internel programs to use this
+%  function, to actually perform a specific task without posible interference
+%  by any API user supplied settings.
+%
+%  It is MorphologyImageChannel() task to extract any such user controls, and
+%  pass them to this function for processing.
 %
 %  More specifically kernels are not normalized/scaled/blended by the
-%  'convolve:scale' Image Artifact (-set setting), and the convolve bias
-%  (-bias setting or image->bias) is passed directly to this function,
-%  and not extracted from an image.
+%  'convolve:scale' Image Artifact (setting), nor is the convolve bias
+%  (-bias setting or image->bias) loooked at, but must be supplied from the
+%  function arguments.
 %
 %  The format of the MorphologyApply method is:
 %
 %      Image *MorphologyApply(const Image *image,MorphologyMethod method,
-%        const long iterations,const KernelInfo *kernel,
-%        const CompositeMethod compose, const double bias,
-%        ExceptionInfo *exception)
+%        const ChannelType channel, const ssize_t iterations,
+%        const KernelInfo *kernel, const CompositeMethod compose,
+%        const double bias, ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
-%    o image: the image.
+%    o image: the source image
 %
 %    o method: the morphology method to be applied.
 %
+%    o channel: the channels to which the operations are applied
+%               The channel 'sync' flag determines if 'alpha weighting' is
+%               applied for convolution style operations.
+%
 %    o iterations: apply the operation this many times (or no change).
 %                  A value of -1 means loop until no change found.
 %                  How this is applied may depend on the morphology method.
@@ -2070,13 +2481,11 @@ static void CalcKernelMetaData(KernelInfo *kernel)
 %    o channel: the channel type.
 %
 %    o kernel: An array of double representing the morphology kernel.
-%              Warning: kernel may be normalized for the Convolve method.
 %
 %    o compose: How to handle or merge multi-kernel results.
-%               If 'Undefined' use default of the Morphology method.
-%               If 'No' force image to be re-iterated by each kernel.
-%               Otherwise merge the results using the mathematical compose
-%               method given.
+%          If 'UndefinedCompositeOp' use default for the Morphology method.
+%          If 'NoCompositeOp' force image to be re-iterated by each kernel.
+%          Otherwise merge the results using the compose method given.
 %
 %    o bias: Convolution Output Bias.
 %
@@ -2084,28 +2493,42 @@ static void CalcKernelMetaData(KernelInfo *kernel)
 %
 */
 
-
 /* Apply a Morphology Primative to an image using the given kernel.
-** Two pre-created images must be provided, no image is created.
-** Returning the number of pixels that changed.
+** Two pre-created images must be provided, and no image is created.
+** It returns the number of pixels that changed between the images
+** for result convergence determination.
 */
-static unsigned long MorphologyPrimitive(const Image *image, Image
-     *result_image, const MorphologyMethod method, const ChannelType channel,
+static ssize_t MorphologyPrimitive(const Image *image, Image *result_image,
+     const MorphologyMethod method, const ChannelType channel,
      const KernelInfo *kernel,const double bias,ExceptionInfo *exception)
 {
 #define MorphologyTag  "Morphology/Image"
 
-  long
-    progress,
-    y, offx, offy,
+  CacheView
+    *p_view,
+    *q_view;
+
+  ssize_t
+    y, offx, offy;
+
+  size_t
+    virt_width,
     changed;
 
   MagickBooleanType
     status;
 
-  CacheView
-    *p_view,
-    *q_view;
+  MagickOffsetType
+    progress;
+
+  assert(image != (Image *) NULL);
+  assert(image->signature == MagickSignature);
+  assert(result_image != (Image *) NULL);
+  assert(result_image->signature == MagickSignature);
+  assert(kernel != (KernelInfo *) NULL);
+  assert(kernel->signature == MagickSignature);
+  assert(exception != (ExceptionInfo *) NULL);
+  assert(exception->signature == MagickSignature);
 
   status=MagickTrue;
   changed=0;
@@ -2113,6 +2536,7 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
 
   p_view=AcquireCacheView(image);
   q_view=AcquireCacheView(result_image);
+  virt_width=image->columns+kernel->width-1;
 
   /* Some methods (including convolve) needs use a reflected kernel.
    * Adjust 'origin' offsets to loop though kernel as a reflection.
@@ -2123,31 +2547,229 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
     case ConvolveMorphology:
     case DilateMorphology:
     case DilateIntensityMorphology:
-    case DistanceMorphology:
+    /*case DistanceMorphology:*/
       /* kernel needs to used with reflection about origin */
-      offx = (long) kernel->width-offx-1;
-      offy = (long) kernel->height-offy-1;
+      offx = (ssize_t) kernel->width-offx-1;
+      offy = (ssize_t) kernel->height-offy-1;
       break;
     case ErodeMorphology:
     case ErodeIntensityMorphology:
     case HitAndMissMorphology:
     case ThinningMorphology:
     case ThickenMorphology:
-      /* kernel is user as is, without reflection */
+      /* kernel is used as is, without reflection */
       break;
     default:
       assert("Not a Primitive Morphology Method" != (char *) NULL);
       break;
   }
 
+
+  if ( method == ConvolveMorphology && kernel->width == 1 )
+  { /* Special handling (for speed) of vertical (blur) kernels.
+    ** This performs its handling in columns rather than in rows.
+    ** This is only done for convolve as it is the only method that
+    ** generates very large 1-D vertical kernels (such as a 'BlurKernel')
+    **
+    ** Timing tests (on single CPU laptop)
+    ** Using a vertical 1-d Blue with normal row-by-row (below)
+    **   time convert logo: -morphology Convolve Blur:0x10+90 null:
+    **      0.807u
+    ** Using this column method
+    **   time convert logo: -morphology Convolve Blur:0x10+90 null:
+    **      0.620u
+    **
+    ** Anthony Thyssen, 14 June 2010
+    */
+    register ssize_t
+      x;
+
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+#pragma omp parallel for schedule(dynamic,4) shared(progress,status)
+#endif
+    for (x=0; x < (ssize_t) image->columns; x++)
+    {
+      register const PixelPacket
+        *restrict p;
+
+      register const IndexPacket
+        *restrict p_indexes;
+
+      register PixelPacket
+        *restrict q;
+
+      register IndexPacket
+        *restrict q_indexes;
+
+      register ssize_t
+        y;
+
+      ssize_t
+        r;
+
+      if (status == MagickFalse)
+        continue;
+      p=GetCacheViewVirtualPixels(p_view, x,  -offy,1,
+          image->rows+kernel->height-1, exception);
+      q=GetCacheViewAuthenticPixels(q_view,x,0,1,result_image->rows,exception);
+      if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
+        {
+          status=MagickFalse;
+          continue;
+        }
+      p_indexes=GetCacheViewVirtualIndexQueue(p_view);
+      q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
+
+      /* offset to origin in 'p'. while 'q' points to it directly */
+      r = offy;
+
+      for (y=0; y < (ssize_t) image->rows; y++)
+      {
+        register ssize_t
+          v;
+
+        register const double
+          *restrict k;
+
+        register const PixelPacket
+          *restrict k_pixels;
+
+        register const IndexPacket
+          *restrict k_indexes;
+
+        MagickPixelPacket
+          result;
+
+        /* Copy input image to the output image for unused channels
+        * This removes need for 'cloning' a new image every iteration
+        */
+        *q = p[r];
+        if (image->colorspace == CMYKColorspace)
+          SetIndexPixelComponent(q_indexes+y,GetIndexPixelComponent(
+            p_indexes+r));
+
+        /* Set the bias of the weighted average output */
+        result.red     =
+        result.green   =
+        result.blue    =
+        result.opacity =
+        result.index   = bias;
+
+
+        /* Weighted Average of pixels using reflected kernel
+        **
+        ** NOTE for correct working of this operation for asymetrical
+        ** kernels, the kernel needs to be applied in its reflected form.
+        ** That is its values needs to be reversed.
+        */
+        k = &kernel->values[ kernel->height-1 ];
+        k_pixels = p;
+        k_indexes = p_indexes;
+        if ( ((channel & SyncChannels) == 0 ) ||
+                             (image->matte == MagickFalse) )
+          { /* No 'Sync' involved.
+            ** Convolution is simple greyscale channel operation
+            */
+            for (v=0; v < (ssize_t) kernel->height; v++) {
+              if ( IsNan(*k) ) continue;
+              result.red     += (*k)*GetRedPixelComponent(k_pixels);
+              result.green   += (*k)*GetGreenPixelComponent(k_pixels);
+              result.blue    += (*k)*GetBluePixelComponent(k_pixels);
+              result.opacity += (*k)*GetOpacityPixelComponent(k_pixels);
+              if ( image->colorspace == CMYKColorspace)
+                result.index += (*k)*(*k_indexes);
+              k--;
+              k_pixels++;
+              k_indexes++;
+            }
+            if ((channel & RedChannel) != 0)
+              SetRedPixelComponent(q,ClampToQuantum(result.red));
+            if ((channel & GreenChannel) != 0)
+              SetGreenPixelComponent(q,ClampToQuantum(result.green));
+            if ((channel & BlueChannel) != 0)
+              SetBluePixelComponent(q,ClampToQuantum(result.blue));
+            if ((channel & OpacityChannel) != 0
+                && image->matte == MagickTrue )
+              SetOpacityPixelComponent(q,ClampToQuantum(result.opacity));
+            if ((channel & IndexChannel) != 0
+                && image->colorspace == CMYKColorspace)
+              SetIndexPixelComponent(q_indexes+x,ClampToQuantum(result.index));
+          }
+        else
+          { /* Channel 'Sync' Flag, and Alpha Channel enabled.
+            ** Weight the color channels with Alpha Channel so that
+            ** transparent pixels are not part of the results.
+            */
+            MagickRealType
+              alpha,  /* alpha weighting of colors : kernel*alpha  */
+              gamma;  /* divisor, sum of color weighting values */
+
+            gamma=0.0;
+            for (v=0; v < (ssize_t) kernel->height; v++) {
+              if ( IsNan(*k) ) continue;
+              alpha=(*k)*(QuantumScale*(QuantumRange-GetOpacityPixelComponent(k_pixels)));
+              gamma += alpha;
+              result.red     += alpha*GetRedPixelComponent(k_pixels);
+              result.green   += alpha*GetGreenPixelComponent(k_pixels);
+              result.blue    += alpha*GetBluePixelComponent(k_pixels);
+              result.opacity += (*k)*GetOpacityPixelComponent(k_pixels);
+              if ( image->colorspace == CMYKColorspace)
+                result.index += alpha*(*k_indexes);
+              k--;
+              k_pixels++;
+              k_indexes++;
+            }
+            /* Sync'ed channels, all channels are modified */
+            gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
+            SetRedPixelComponent(q,ClampToQuantum(gamma*result.red));
+            SetGreenPixelComponent(q,ClampToQuantum(gamma*result.green));
+            SetBluePixelComponent(q,ClampToQuantum(gamma*result.blue));
+            SetOpacityPixelComponent(q,ClampToQuantum(result.opacity));
+            if (image->colorspace == CMYKColorspace)
+              SetIndexPixelComponent(q_indexes+x,ClampToQuantum(gamma*
+                result.index));
+          }
+
+        /* Count up changed pixels */
+        if (   ( p[r].red != GetRedPixelComponent(q))
+            || ( p[r].green != GetGreenPixelComponent(q))
+            || ( p[r].blue != GetBluePixelComponent(q))
+            || ( p[r].opacity != GetOpacityPixelComponent(q))
+            || ( image->colorspace == CMYKColorspace &&
+                GetIndexPixelComponent(p_indexes+r) != GetIndexPixelComponent(q_indexes+x) ) )
+          changed++;  /* The pixel was changed in some way! */
+        p++;
+        q++;
+      } /* y */
+      if ( SyncCacheViewAuthenticPixels(q_view,exception) == MagickFalse)
+        status=MagickFalse;
+      if (image->progress_monitor != (MagickProgressMonitor) NULL)
+        {
+          MagickBooleanType
+            proceed;
+
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp critical (MagickCore_MorphologyImage)
+#endif
+          proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
+          if (proceed == MagickFalse)
+            status=MagickFalse;
+        }
+    } /* x */
+    result_image->type=image->type;
+    q_view=DestroyCacheView(q_view);
+    p_view=DestroyCacheView(p_view);
+    return(status ? (ssize_t) changed : 0);
+  }
+
+  /*
+  ** Normal handling of horizontal or rectangular kernels (row by row)
+  */
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
   #pragma omp parallel for schedule(dynamic,4) shared(progress,status)
 #endif
-  for (y=0; y < (long) image->rows; y++)
+  for (y=0; y < (ssize_t) image->rows; y++)
   {
-    MagickBooleanType
-      sync;
-
     register const PixelPacket
       *restrict p;
 
@@ -2160,16 +2782,16 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
     register IndexPacket
       *restrict q_indexes;
 
-    register long
+    register ssize_t
       x;
 
-    unsigned long
+    size_t
       r;
 
     if (status == MagickFalse)
       continue;
-    p=GetCacheViewVirtualPixels(p_view, -offx,  y-offy,
-         image->columns+kernel->width,  kernel->height,  exception);
+    p=GetCacheViewVirtualPixels(p_view, -offx, y-offy, virt_width,
+         kernel->height,  exception);
     q=GetCacheViewAuthenticPixels(q_view,0,y,result_image->columns,1,
          exception);
     if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
@@ -2179,14 +2801,16 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
       }
     p_indexes=GetCacheViewVirtualIndexQueue(p_view);
     q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
-    r = (image->columns+kernel->width)*offy+offx; /* constant */
 
-    for (x=0; x < (long) image->columns; x++)
+    /* offset to origin in 'p'. while 'q' points to it directly */
+    r = virt_width*offy + offx;
+
+    for (x=0; x < (ssize_t) image->columns; x++)
     {
-       long
+       ssize_t
         v;
 
-      register long
+      register ssize_t
         u;
 
       register const double
@@ -2203,12 +2827,12 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
         min,
         max;
 
-      /* Copy input to ouput image for unused channels
+      /* Copy input image to the output image for unused channels
        * This removes need for 'cloning' a new image every iteration
        */
       *q = p[r];
       if (image->colorspace == CMYKColorspace)
-        q_indexes[x] = p_indexes[r];
+        SetIndexPixelComponent(q_indexes+x,GetIndexPixelComponent(p_indexes+r));
 
       /* Defaults */
       min.red     =
@@ -2228,11 +2852,11 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
       result.opacity = QuantumRange - (MagickRealType) p[r].opacity;
       result.index   = 0.0;
       if ( image->colorspace == CMYKColorspace)
-         result.index   = (MagickRealType) p_indexes[r];
+         result.index   = (MagickRealType) GetIndexPixelComponent(p_indexes+r);
 
       switch (method) {
         case ConvolveMorphology:
-          /* Set the user defined bias of the weighted average output */
+          /* Set the bias of the weighted average output */
           result.red     =
           result.green   =
           result.blue    =
@@ -2268,22 +2892,53 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
             ** For more details of Correlation vs Convolution see
             **   http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
             */
-            if (((channel & SyncChannels) != 0 ) &&
-                      (image->matte == MagickTrue))
-              { /* Channel has a 'Sync' Flag, and Alpha Channel enabled.
+            k = &kernel->values[ kernel->width*kernel->height-1 ];
+            k_pixels = p;
+            k_indexes = p_indexes;
+            if ( ((channel & SyncChannels) == 0 ) ||
+                                 (image->matte == MagickFalse) )
+              { /* No 'Sync' involved.
+                ** Convolution is simple greyscale channel operation
+                */
+                for (v=0; v < (ssize_t) kernel->height; v++) {
+                  for (u=0; u < (ssize_t) kernel->width; u++, k--) {
+                    if ( IsNan(*k) ) continue;
+                    result.red     += (*k)*k_pixels[u].red;
+                    result.green   += (*k)*k_pixels[u].green;
+                    result.blue    += (*k)*k_pixels[u].blue;
+                    result.opacity += (*k)*k_pixels[u].opacity;
+                    if ( image->colorspace == CMYKColorspace)
+                      result.index += (*k)*GetIndexPixelComponent(k_indexes+u);
+                  }
+                  k_pixels += virt_width;
+                  k_indexes += virt_width;
+                }
+                if ((channel & RedChannel) != 0)
+                  SetRedPixelComponent(q,ClampToQuantum(result.red));
+                if ((channel & GreenChannel) != 0)
+                  SetGreenPixelComponent(q,ClampToQuantum(result.green));
+                if ((channel & BlueChannel) != 0)
+                  SetBluePixelComponent(q,ClampToQuantum(result.blue));
+                if ((channel & OpacityChannel) != 0
+                    && image->matte == MagickTrue )
+                  SetOpacityPixelComponent(q,ClampToQuantum(result.opacity));
+                if ((channel & IndexChannel) != 0
+                    && image->colorspace == CMYKColorspace)
+                  SetIndexPixelComponent(q_indexes+x,ClampToQuantum(
+                    result.index));
+              }
+            else
+              { /* Channel 'Sync' Flag, and Alpha Channel enabled.
                 ** Weight the color channels with Alpha Channel so that
                 ** transparent pixels are not part of the results.
                 */
                 MagickRealType
-                  alpha,  /* color channel weighting : kernel*alpha  */
-                  gamma;  /* divisor, sum of weighting values */
+                  alpha,  /* alpha weighting of colors : kernel*alpha  */
+                  gamma;  /* divisor, sum of color weighting values */
 
                 gamma=0.0;
-                k = &kernel->values[ kernel->width*kernel->height-1 ];
-                k_pixels = p;
-                k_indexes = p_indexes;
-                for (v=0; v < (long) kernel->height; v++) {
-                  for (u=0; u < (long) kernel->width; u++, k--) {
+                for (v=0; v < (ssize_t) kernel->height; v++) {
+                  for (u=0; u < (ssize_t) kernel->width; u++, k--) {
                     if ( IsNan(*k) ) continue;
                     alpha=(*k)*(QuantumScale*(QuantumRange-
                                           k_pixels[u].opacity));
@@ -2291,41 +2946,22 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
                     result.red     += alpha*k_pixels[u].red;
                     result.green   += alpha*k_pixels[u].green;
                     result.blue    += alpha*k_pixels[u].blue;
-                    result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
+                    result.opacity += (*k)*k_pixels[u].opacity;
                     if ( image->colorspace == CMYKColorspace)
-                      result.index   += alpha*k_indexes[u];
+                      result.index+=alpha*GetIndexPixelComponent(k_indexes+u);
                   }
-                  k_pixels += image->columns+kernel->width;
-                  k_indexes += image->columns+kernel->width;
+                  k_pixels += virt_width;
+                  k_indexes += virt_width;
                 }
+                /* Sync'ed channels, all channels are modified */
                 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
-                result.red *= gamma;
-                result.green *= gamma;
-                result.blue *= gamma;
-                result.opacity *= gamma;
-                result.index *= gamma;
-              }
-            else
-              {
-                /* No 'Sync' flag, or no Alpha involved.
-                ** Convolution is simple individual channel weigthed sum.
-                */
-                k = &kernel->values[ kernel->width*kernel->height-1 ];
-                k_pixels = p;
-                k_indexes = p_indexes;
-                for (v=0; v < (long) kernel->height; v++) {
-                  for (u=0; u < (long) kernel->width; u++, k--) {
-                    if ( IsNan(*k) ) continue;
-                    result.red     += (*k)*k_pixels[u].red;
-                    result.green   += (*k)*k_pixels[u].green;
-                    result.blue    += (*k)*k_pixels[u].blue;
-                    result.opacity += (*k)*(QuantumRange-k_pixels[u].opacity);
-                    if ( image->colorspace == CMYKColorspace)
-                      result.index   += (*k)*k_indexes[u];
-                  }
-                  k_pixels += image->columns+kernel->width;
-                  k_indexes += image->columns+kernel->width;
-                }
+                SetRedPixelComponent(q,ClampToQuantum(gamma*result.red));
+                SetGreenPixelComponent(q,ClampToQuantum(gamma*result.green));
+                SetBluePixelComponent(q,ClampToQuantum(gamma*result.blue));
+                SetOpacityPixelComponent(q,ClampToQuantum(result.opacity));
+                if (image->colorspace == CMYKColorspace)
+                  SetIndexPixelComponent(q_indexes+x,ClampToQuantum(gamma*
+                   result.index));
               }
             break;
 
@@ -2341,8 +2977,8 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
             k = kernel->values;
             k_pixels = p;
             k_indexes = p_indexes;
-            for (v=0; v < (long) kernel->height; v++) {
-              for (u=0; u < (long) kernel->width; u++, k++) {
+            for (v=0; v < (ssize_t) kernel->height; v++) {
+              for (u=0; u < (ssize_t) kernel->width; u++, k++) {
                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
                 Minimize(min.red,     (double) k_pixels[u].red);
                 Minimize(min.green,   (double) k_pixels[u].green);
@@ -2350,14 +2986,14 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
                 Minimize(min.opacity,
                             QuantumRange-(double) k_pixels[u].opacity);
                 if ( image->colorspace == CMYKColorspace)
-                  Minimize(min.index,   (double) k_indexes[u]);
+                  Minimize(min.index,(double) GetIndexPixelComponent(
+                    k_indexes+u));
               }
-              k_pixels += image->columns+kernel->width;
-              k_indexes += image->columns+kernel->width;
+              k_pixels += virt_width;
+              k_indexes += virt_width;
             }
             break;
 
-
         case DilateMorphology:
             /* Maximum Value within kernel neighbourhood
             **
@@ -2373,8 +3009,8 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
             k = &kernel->values[ kernel->width*kernel->height-1 ];
             k_pixels = p;
             k_indexes = p_indexes;
-            for (v=0; v < (long) kernel->height; v++) {
-              for (u=0; u < (long) kernel->width; u++, k--) {
+            for (v=0; v < (ssize_t) kernel->height; v++) {
+              for (u=0; u < (ssize_t) kernel->width; u++, k--) {
                 if ( IsNan(*k) || (*k) < 0.5 ) continue;
                 Maximize(max.red,     (double) k_pixels[u].red);
                 Maximize(max.green,   (double) k_pixels[u].green);
@@ -2382,10 +3018,11 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
                 Maximize(max.opacity,
                             QuantumRange-(double) k_pixels[u].opacity);
                 if ( image->colorspace == CMYKColorspace)
-                  Maximize(max.index,   (double) k_indexes[u]);
+                  Maximize(max.index,   (double) GetIndexPixelComponent(
+                    k_indexes+u));
               }
-              k_pixels += image->columns+kernel->width;
-              k_indexes += image->columns+kernel->width;
+              k_pixels += virt_width;
+              k_indexes += virt_width;
             }
             break;
 
@@ -2399,14 +3036,15 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
             ** neighbourhoods, 0.0 for background, and 1.0 for foreground
             ** with either Nan or 0.5 values for don't care.
             **
-            ** Note that this can produce negative results, though really
-            ** only a positive match has any real value.
+            ** Note that this will never produce a meaningless negative
+            ** result.  Such results can cause Thinning/Thicken to not work
+            ** correctly when used against a greyscale image.
             */
             k = kernel->values;
             k_pixels = p;
             k_indexes = p_indexes;
-            for (v=0; v < (long) kernel->height; v++) {
-              for (u=0; u < (long) kernel->width; u++, k++) {
+            for (v=0; v < (ssize_t) kernel->height; v++) {
+              for (u=0; u < (ssize_t) kernel->width; u++, k++) {
                 if ( IsNan(*k) ) continue;
                 if ( (*k) > 0.7 )
                 { /* minimim of foreground pixels */
@@ -2416,7 +3054,8 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
                   Minimize(min.opacity,
                               QuantumRange-(double) k_pixels[u].opacity);
                   if ( image->colorspace == CMYKColorspace)
-                    Minimize(min.index,   (double) k_indexes[u]);
+                    Minimize(min.index,(double) GetIndexPixelComponent(
+                      k_indexes+u));
                 }
                 else if ( (*k) < 0.3 )
                 { /* maximum of background pixels */
@@ -2426,13 +3065,14 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
                   Maximize(max.opacity,
                               QuantumRange-(double) k_pixels[u].opacity);
                   if ( image->colorspace == CMYKColorspace)
-                    Maximize(max.index,   (double) k_indexes[u]);
+                    Maximize(max.index,   (double) GetIndexPixelComponent(
+                      k_indexes+u));
                 }
               }
-              k_pixels += image->columns+kernel->width;
-              k_indexes += image->columns+kernel->width;
+              k_pixels += virt_width;
+              k_indexes += virt_width;
             }
-            /* Pattern Match  only if min fg larger than min bg pixels */
+            /* Pattern Match if difference is positive */
             min.red     -= max.red;     Maximize( min.red,     0.0 );
             min.green   -= max.green;   Maximize( min.green,   0.0 );
             min.blue    -= max.blue;    Maximize( min.blue,    0.0 );
@@ -2444,7 +3084,7 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
             /* Select Pixel with Minimum Intensity within kernel neighbourhood
             **
             ** WARNING: the intensity test fails for CMYK and does not
-            ** take into account the moderating effect of teh alpha channel
+            ** take into account the moderating effect of the alpha channel
             ** on the intensity.
             **
             ** NOTE that the kernel is not reflected for this operation!
@@ -2452,185 +3092,649 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
             k = kernel->values;
             k_pixels = p;
             k_indexes = p_indexes;
-            for (v=0; v < (long) kernel->height; v++) {
-              for (u=0; u < (long) kernel->width; u++, k++) {
-                if ( IsNan(*k) || (*k) < 0.5 ) continue;
-                if ( result.red == 0.0 ||
-                     PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
-                  /* copy the whole pixel - no channel selection */
-                  *q = k_pixels[u];
-                  if ( result.red > 0.0 ) changed++;
-                  result.red = 1.0;
-                }
+            for (v=0; v < (ssize_t) kernel->height; v++) {
+              for (u=0; u < (ssize_t) kernel->width; u++, k++) {
+                if ( IsNan(*k) || (*k) < 0.5 ) continue;
+                if ( result.red == 0.0 ||
+                     PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
+                  /* copy the whole pixel - no channel selection */
+                  *q = k_pixels[u];
+                  if ( result.red > 0.0 ) changed++;
+                  result.red = 1.0;
+                }
+              }
+              k_pixels += virt_width;
+              k_indexes += virt_width;
+            }
+            break;
+
+        case DilateIntensityMorphology:
+            /* Select Pixel with Maximum Intensity within kernel neighbourhood
+            **
+            ** WARNING: the intensity test fails for CMYK and does not
+            ** take into account the moderating effect of the alpha channel
+            ** on the intensity (yet).
+            **
+            ** NOTE for correct working of this operation for asymetrical
+            ** kernels, the kernel needs to be applied in its reflected form.
+            ** That is its values needs to be reversed.
+            */
+            k = &kernel->values[ kernel->width*kernel->height-1 ];
+            k_pixels = p;
+            k_indexes = p_indexes;
+            for (v=0; v < (ssize_t) kernel->height; v++) {
+              for (u=0; u < (ssize_t) kernel->width; u++, k--) {
+                if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
+                if ( result.red == 0.0 ||
+                     PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
+                  /* copy the whole pixel - no channel selection */
+                  *q = k_pixels[u];
+                  if ( result.red > 0.0 ) changed++;
+                  result.red = 1.0;
+                }
+              }
+              k_pixels += virt_width;
+              k_indexes += virt_width;
+            }
+            break;
+#if 0
+  This code has been obsoleted by the MorphologyPrimitiveDirect() function.
+  However it is still (almost) correct coding for Grayscale Morphology.
+  That is...
+
+  GrayErode    is equivalent but with kernel values subtracted from pixels
+               without the kernel rotation
+  GreyDilate   is equivalent but using Maximum() instead of Minimum()
+               useing kernel rotation
+
+        case DistanceMorphology:
+            /* Add kernel Value and select the minimum value found.
+            ** The result is a iterative distance from edge of image shape.
+            **
+            ** All Distance Kernels are symetrical, but that may not always
+            ** be the case. For example how about a distance from left edges?
+            ** To work correctly with asymetrical kernels the reflected kernel
+            ** needs to be applied.
+            */
+            k = &kernel->values[ kernel->width*kernel->height-1 ];
+            k_pixels = p;
+            k_indexes = p_indexes;
+            for (v=0; v < (ssize_t) kernel->height; v++) {
+              for (u=0; u < (ssize_t) kernel->width; u++, k--) {
+                if ( IsNan(*k) ) continue;
+                Minimize(result.red,     (*k)+k_pixels[u].red);
+                Minimize(result.green,   (*k)+k_pixels[u].green);
+                Minimize(result.blue,    (*k)+k_pixels[u].blue);
+                Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
+                if ( image->colorspace == CMYKColorspace)
+                  Minimize(result.index,(*k)+GetIndexPixelComponent(
+                    k_indexes+u));
+              }
+              k_pixels += virt_width;
+              k_indexes += virt_width;
+            }
+            break;
+#endif
+        case UndefinedMorphology:
+        default:
+            break; /* Do nothing */
+      }
+      /* Final mathematics of results (combine with original image?)
+      **
+      ** NOTE: Difference Morphology operators Edge* and *Hat could also
+      ** be done here but works better with iteration as a image difference
+      ** in the controling function (below).  Thicken and Thinning however
+      ** should be done here so thay can be iterated correctly.
+      */
+      switch ( method ) {
+        case HitAndMissMorphology:
+        case ErodeMorphology:
+          result = min;    /* minimum of neighbourhood */
+          break;
+        case DilateMorphology:
+          result = max;    /* maximum of neighbourhood */
+          break;
+        case ThinningMorphology:
+          /* subtract pattern match from original */
+          result.red     -= min.red;
+          result.green   -= min.green;
+          result.blue    -= min.blue;
+          result.opacity -= min.opacity;
+          result.index   -= min.index;
+          break;
+        case ThickenMorphology:
+          /* Add the pattern matchs to the original */
+          result.red     += min.red;
+          result.green   += min.green;
+          result.blue    += min.blue;
+          result.opacity += min.opacity;
+          result.index   += min.index;
+          break;
+        default:
+          /* result directly calculated or assigned */
+          break;
+      }
+      /* Assign the resulting pixel values - Clamping Result */
+      switch ( method ) {
+        case UndefinedMorphology:
+        case ConvolveMorphology:
+        case DilateIntensityMorphology:
+        case ErodeIntensityMorphology:
+          break;  /* full pixel was directly assigned - not a channel method */
+        default:
+          if ((channel & RedChannel) != 0)
+            SetRedPixelComponent(q,ClampToQuantum(result.red));
+          if ((channel & GreenChannel) != 0)
+            SetGreenPixelComponent(q,ClampToQuantum(result.green));
+          if ((channel & BlueChannel) != 0)
+            SetBluePixelComponent(q,ClampToQuantum(result.blue));
+          if ((channel & OpacityChannel) != 0
+              && image->matte == MagickTrue )
+            SetOpacityPixelComponent(q,ClampToQuantum(QuantumRange-result.opacity));
+          if ((channel & IndexChannel) != 0
+              && image->colorspace == CMYKColorspace)
+            SetIndexPixelComponent(q_indexes+x,ClampToQuantum(result.index));
+          break;
+      }
+      /* Count up changed pixels */
+      if (   ( p[r].red != GetRedPixelComponent(q) )
+          || ( p[r].green != GetGreenPixelComponent(q) )
+          || ( p[r].blue != GetBluePixelComponent(q) )
+          || ( p[r].opacity != GetOpacityPixelComponent(q) )
+          || ( image->colorspace == CMYKColorspace &&
+               GetIndexPixelComponent(p_indexes+r) != GetIndexPixelComponent(q_indexes+x) ) )
+        changed++;  /* The pixel was changed in some way! */
+      p++;
+      q++;
+    } /* x */
+    if ( SyncCacheViewAuthenticPixels(q_view,exception) == MagickFalse)
+      status=MagickFalse;
+    if (image->progress_monitor != (MagickProgressMonitor) NULL)
+      {
+        MagickBooleanType
+          proceed;
+
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp critical (MagickCore_MorphologyImage)
+#endif
+        proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
+        if (proceed == MagickFalse)
+          status=MagickFalse;
+      }
+  } /* y */
+  q_view=DestroyCacheView(q_view);
+  p_view=DestroyCacheView(p_view);
+  return(status ? (ssize_t)changed : -1);
+}
+
+/* This is almost identical to the MorphologyPrimative() function above,
+** but will apply the primitive directly to the image in two passes.
+**
+** That is after each row is 'Sync'ed' into the image, the next row will
+** make use of those values as part of the calculation of the next row.
+** It then repeats, but going in the oppisite (bottom-up) direction.
+**
+** Because of this 'iterative' handling this function can not make use
+** of multi-threaded, parellel processing.
+*/
+static ssize_t MorphologyPrimitiveDirect(Image *image,
+     const MorphologyMethod method, const ChannelType channel,
+     const KernelInfo *kernel,ExceptionInfo *exception)
+{
+  CacheView
+    *auth_view,
+    *virt_view;
+
+  MagickBooleanType
+    status;
+
+  MagickOffsetType
+    progress;
+
+  ssize_t
+    y, offx, offy;
+
+  size_t
+    virt_width,
+    changed;
+
+  status=MagickTrue;
+  changed=0;
+  progress=0;
+
+  assert(image != (Image *) NULL);
+  assert(image->signature == MagickSignature);
+  assert(kernel != (KernelInfo *) NULL);
+  assert(kernel->signature == MagickSignature);
+  assert(exception != (ExceptionInfo *) NULL);
+  assert(exception->signature == MagickSignature);
+
+  /* Some methods (including convolve) needs use a reflected kernel.
+   * Adjust 'origin' offsets to loop though kernel as a reflection.
+   */
+  offx = kernel->x;
+  offy = kernel->y;
+  switch(method) {
+    case DistanceMorphology:
+    case VoronoiMorphology:
+      /* kernel needs to used with reflection about origin */
+      offx = (ssize_t) kernel->width-offx-1;
+      offy = (ssize_t) kernel->height-offy-1;
+      break;
+#if 0
+    case ?????Morphology:
+      /* kernel is used as is, without reflection */
+      break;
+#endif
+    default:
+      assert("Not a PrimativeDirect Morphology Method" != (char *) NULL);
+      break;
+  }
+
+  /* DO NOT THREAD THIS CODE! */
+  /* two views into same image (virtual, and actual) */
+  virt_view=AcquireCacheView(image);
+  auth_view=AcquireCacheView(image);
+  virt_width=image->columns+kernel->width-1;
+
+  for (y=0; y < (ssize_t) image->rows; y++)
+  {
+    register const PixelPacket
+      *restrict p;
+
+    register const IndexPacket
+      *restrict p_indexes;
+
+    register PixelPacket
+      *restrict q;
+
+    register IndexPacket
+      *restrict q_indexes;
+
+    register ssize_t
+      x;
+
+    ssize_t
+      r;
+
+    /* NOTE read virtual pixels, and authentic pixels, from the same image!
+    ** we read using virtual to get virtual pixel handling, but write back
+    ** into the same image.
+    **
+    ** Only top half of kernel is processed as we do a single pass downward
+    ** through the image iterating the distance function as we go.
+    */
+    if (status == MagickFalse)
+      break;
+    p=GetCacheViewVirtualPixels(virt_view, -offx,  y-offy, virt_width, (size_t) offy+1,
+         exception);
+    q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
+         exception);
+    if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
+      status=MagickFalse;
+    if (status == MagickFalse)
+      break;
+    p_indexes=GetCacheViewVirtualIndexQueue(virt_view);
+    q_indexes=GetCacheViewAuthenticIndexQueue(auth_view);
+
+    /* offset to origin in 'p'. while 'q' points to it directly */
+    r = (ssize_t) virt_width*offy + offx;
+
+    for (x=0; x < (ssize_t) image->columns; x++)
+    {
+      ssize_t
+        v;
+
+      register ssize_t
+        u;
+
+      register const double
+        *restrict k;
+
+      register const PixelPacket
+        *restrict k_pixels;
+
+      register const IndexPacket
+        *restrict k_indexes;
+
+      MagickPixelPacket
+        result;
+
+      /* Starting Defaults */
+      GetMagickPixelPacket(image,&result);
+      SetMagickPixelPacket(image,q,q_indexes,&result);
+      if ( method != VoronoiMorphology )
+        result.opacity = QuantumRange - result.opacity;
+
+      switch ( method ) {
+        case DistanceMorphology:
+            /* Add kernel Value and select the minimum value found. */
+            k = &kernel->values[ kernel->width*kernel->height-1 ];
+            k_pixels = p;
+            k_indexes = p_indexes;
+            for (v=0; v <= (ssize_t) offy; v++) {
+              for (u=0; u < (ssize_t) kernel->width; u++, k--) {
+                if ( IsNan(*k) ) continue;
+                Minimize(result.red,     (*k)+k_pixels[u].red);
+                Minimize(result.green,   (*k)+k_pixels[u].green);
+                Minimize(result.blue,    (*k)+k_pixels[u].blue);
+                Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
+                if ( image->colorspace == CMYKColorspace)
+                  Minimize(result.index,   (*k)+GetIndexPixelComponent(k_indexes+u));
               }
-              k_pixels += image->columns+kernel->width;
-              k_indexes += image->columns+kernel->width;
+              k_pixels += virt_width;
+              k_indexes += virt_width;
             }
+            /* repeat with the just processed pixels of this row */
+            k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
+            k_pixels = q-offx;
+            k_indexes = q_indexes-offx;
+              for (u=0; u < (ssize_t) offx; u++, k--) {
+                if ( x+u-offx < 0 ) continue;  /* off the edge! */
+                if ( IsNan(*k) ) continue;
+                Minimize(result.red,     (*k)+k_pixels[u].red);
+                Minimize(result.green,   (*k)+k_pixels[u].green);
+                Minimize(result.blue,    (*k)+k_pixels[u].blue);
+                Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
+                if ( image->colorspace == CMYKColorspace)
+                  Minimize(result.index,   (*k)+GetIndexPixelComponent(k_indexes+u));
+              }
             break;
-
-        case DilateIntensityMorphology:
-            /* Select Pixel with Maximum Intensity within kernel neighbourhood
-            **
-            ** WARNING: the intensity test fails for CMYK and does not
-            ** take into account the moderating effect of the alpha channel
-            ** on the intensity (yet).
+        case VoronoiMorphology:
+            /* Apply Distance to 'Matte' channel, coping the closest color.
             **
-            ** NOTE for correct working of this operation for asymetrical
-            ** kernels, the kernel needs to be applied in its reflected form.
-            ** That is its values needs to be reversed.
+            ** This is experimental, and realy the 'alpha' component should
+            ** be completely separate 'masking' channel.
             */
             k = &kernel->values[ kernel->width*kernel->height-1 ];
             k_pixels = p;
             k_indexes = p_indexes;
-            for (v=0; v < (long) kernel->height; v++) {
-              for (u=0; u < (long) kernel->width; u++, k--) {
-                if ( IsNan(*k) || (*k) < 0.5 ) continue; /* boolean kernel */
-                if ( result.red == 0.0 ||
-                     PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
-                  /* copy the whole pixel - no channel selection */
-                  *q = k_pixels[u];
-                  if ( result.red > 0.0 ) changed++;
-                  result.red = 1.0;
-                }
+            for (v=0; v <= (ssize_t) offy; v++) {
+              for (u=0; u < (ssize_t) kernel->width; u++, k--) {
+                if ( IsNan(*k) ) continue;
+                if( result.opacity > (*k)+k_pixels[u].opacity )
+                  {
+                    SetMagickPixelPacket(image,&k_pixels[u],&k_indexes[u],
+                      &result);
+                    result.opacity += *k;
+                  }
               }
-              k_pixels += image->columns+kernel->width;
-              k_indexes += image->columns+kernel->width;
+              k_pixels += virt_width;
+              k_indexes += virt_width;
             }
+            /* repeat with the just processed pixels of this row */
+            k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
+            k_pixels = q-offx;
+            k_indexes = q_indexes-offx;
+              for (u=0; u < (ssize_t) offx; u++, k--) {
+                if ( x+u-offx < 0 ) continue;  /* off the edge! */
+                if ( IsNan(*k) ) continue;
+                if( result.opacity > (*k)+k_pixels[u].opacity )
+                  {
+                    SetMagickPixelPacket(image,&k_pixels[u],&k_indexes[u],
+                      &result);
+                    result.opacity += *k;
+                  }
+              }
             break;
+        default:
+          /* result directly calculated or assigned */
+          break;
+      }
+      /* Assign the resulting pixel values - Clamping Result */
+      switch ( method ) {
+        case VoronoiMorphology:
+          SetPixelPacket(image,&result,q,q_indexes);
+          break;
+        default:
+          if ((channel & RedChannel) != 0)
+            SetRedPixelComponent(q,ClampToQuantum(result.red));
+          if ((channel & GreenChannel) != 0)
+            SetGreenPixelComponent(q,ClampToQuantum(result.green));
+          if ((channel & BlueChannel) != 0)
+            SetBluePixelComponent(q,ClampToQuantum(result.blue));
+          if ((channel & OpacityChannel) != 0 && image->matte == MagickTrue )
+            SetOpacityPixelComponent(q,ClampToQuantum(QuantumRange-result.opacity));
+          if ((channel & IndexChannel) != 0
+              && image->colorspace == CMYKColorspace)
+            SetIndexPixelComponent(q_indexes+x,ClampToQuantum(result.index));
+          break;
+      }
+      /* Count up changed pixels */
+      if (   ( p[r].red != GetRedPixelComponent(q) )
+          || ( p[r].green != GetGreenPixelComponent(q) )
+          || ( p[r].blue != GetBluePixelComponent(q) )
+          || ( p[r].opacity != GetOpacityPixelComponent(q) )
+          || ( image->colorspace == CMYKColorspace &&
+               GetIndexPixelComponent(p_indexes+r) != GetIndexPixelComponent(q_indexes+x) ) )
+        changed++;  /* The pixel was changed in some way! */
+
+      p++; /* increment pixel buffers */
+      q++;
+    } /* x */
+
+    if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
+      status=MagickFalse;
+    if (image->progress_monitor != (MagickProgressMonitor) NULL)
+      if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
+                == MagickFalse )
+        status=MagickFalse;
+
+  } /* y */
+
+  /* Do the reversed pass through the image */
+  for (y=(ssize_t)image->rows-1; y >= 0; y--)
+  {
+    register const PixelPacket
+      *restrict p;
+
+    register const IndexPacket
+      *restrict p_indexes;
+
+    register PixelPacket
+      *restrict q;
+
+    register IndexPacket
+      *restrict q_indexes;
+
+    register ssize_t
+      x;
+
+    ssize_t
+      r;
+
+    if (status == MagickFalse)
+      break;
+    /* NOTE read virtual pixels, and authentic pixels, from the same image!
+    ** we read using virtual to get virtual pixel handling, but write back
+    ** into the same image.
+    **
+    ** Only the bottom half of the kernel will be processes as we
+    ** up the image.
+    */
+    p=GetCacheViewVirtualPixels(virt_view, -offx, y, virt_width, (size_t) kernel->y+1,
+         exception);
+    q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
+         exception);
+    if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
+      status=MagickFalse;
+    if (status == MagickFalse)
+      break;
+    p_indexes=GetCacheViewVirtualIndexQueue(virt_view);
+    q_indexes=GetCacheViewAuthenticIndexQueue(auth_view);
+
+    /* adjust positions to end of row */
+    p += image->columns-1;
+    q += image->columns-1;
+
+    /* offset to origin in 'p'. while 'q' points to it directly */
+    r = offx;
+
+    for (x=(ssize_t)image->columns-1; x >= 0; x--)
+    {
+      ssize_t
+        v;
+
+      register ssize_t
+        u;
+
+      register const double
+        *restrict k;
+
+      register const PixelPacket
+        *restrict k_pixels;
+
+      register const IndexPacket
+        *restrict k_indexes;
+
+      MagickPixelPacket
+        result;
 
+      /* Default - previously modified pixel */
+      GetMagickPixelPacket(image,&result);
+      SetMagickPixelPacket(image,q,q_indexes,&result);
+      if ( method != VoronoiMorphology )
+        result.opacity = QuantumRange - result.opacity;
 
+      switch ( method ) {
         case DistanceMorphology:
-            /* Add kernel Value and select the minimum value found.
-            ** The result is a iterative distance from edge of image shape.
-            **
-            ** All Distance Kernels are symetrical, but that may not always
-            ** be the case. For example how about a distance from left edges?
-            ** To work correctly with asymetrical kernels the reflected kernel
-            ** needs to be applied.
-            **
-            ** Actually this is really a GreyErode with a negative kernel!
-            **
-            */
-            k = &kernel->values[ kernel->width*kernel->height-1 ];
+            /* Add kernel Value and select the minimum value found. */
+            k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
             k_pixels = p;
             k_indexes = p_indexes;
-            for (v=0; v < (long) kernel->height; v++) {
-              for (u=0; u < (long) kernel->width; u++, k--) {
+            for (v=offy; v < (ssize_t) kernel->height; v++) {
+              for (u=0; u < (ssize_t) kernel->width; u++, k--) {
+                if ( IsNan(*k) ) continue;
+                Minimize(result.red,     (*k)+k_pixels[u].red);
+                Minimize(result.green,   (*k)+k_pixels[u].green);
+                Minimize(result.blue,    (*k)+k_pixels[u].blue);
+                Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
+                if ( image->colorspace == CMYKColorspace)
+                  Minimize(result.index,(*k)+GetIndexPixelComponent(k_indexes+u));
+              }
+              k_pixels += virt_width;
+              k_indexes += virt_width;
+            }
+            /* repeat with the just processed pixels of this row */
+            k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
+            k_pixels = q-offx;
+            k_indexes = q_indexes-offx;
+              for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
+                if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
                 if ( IsNan(*k) ) continue;
                 Minimize(result.red,     (*k)+k_pixels[u].red);
                 Minimize(result.green,   (*k)+k_pixels[u].green);
                 Minimize(result.blue,    (*k)+k_pixels[u].blue);
                 Minimize(result.opacity, (*k)+QuantumRange-k_pixels[u].opacity);
                 if ( image->colorspace == CMYKColorspace)
-                  Minimize(result.index,   (*k)+k_indexes[u]);
+                  Minimize(result.index,   (*k)+GetIndexPixelComponent(k_indexes+u));
+              }
+            break;
+        case VoronoiMorphology:
+            /* Apply Distance to 'Matte' channel, coping the closest color.
+            **
+            ** This is experimental, and realy the 'alpha' component should
+            ** be completely separate 'masking' channel.
+            */
+            k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
+            k_pixels = p;
+            k_indexes = p_indexes;
+            for (v=offy; v < (ssize_t) kernel->height; v++) {
+              for (u=0; u < (ssize_t) kernel->width; u++, k--) {
+                if ( IsNan(*k) ) continue;
+                if( result.opacity > (*k)+k_pixels[u].opacity )
+                  {
+                    SetMagickPixelPacket(image,&k_pixels[u],&k_indexes[u],
+                      &result);
+                    result.opacity += *k;
+                  }
               }
-              k_pixels += image->columns+kernel->width;
-              k_indexes += image->columns+kernel->width;
+              k_pixels += virt_width;
+              k_indexes += virt_width;
             }
+            /* repeat with the just processed pixels of this row */
+            k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
+            k_pixels = q-offx;
+            k_indexes = q_indexes-offx;
+              for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
+                if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
+                if ( IsNan(*k) ) continue;
+                if( result.opacity > (*k)+k_pixels[u].opacity )
+                  {
+                    SetMagickPixelPacket(image,&k_pixels[u],&k_indexes[u],
+                      &result);
+                    result.opacity += *k;
+                  }
+              }
             break;
-
-        case UndefinedMorphology:
-        default:
-            break; /* Do nothing */
-      }
-      /* Final mathematics of results (combine with original image?)
-      **
-      ** NOTE: Difference Morphology operators Edge* and *Hat could also
-      ** be done here but works better with iteration as a image difference
-      ** in the controling function (below).  Thicken and Thinning however
-      ** should be done here so thay can be iterated correctly.
-      */
-      switch ( method ) {
-        case HitAndMissMorphology:
-        case ErodeMorphology:
-          result = min;    /* minimum of neighbourhood */
-          break;
-        case DilateMorphology:
-          result = max;    /* maximum of neighbourhood */
-          break;
-        case ThinningMorphology:
-          /* subtract pattern match from original */
-          result.red     -= min.red;
-          result.green   -= min.green;
-          result.blue    -= min.blue;
-          result.opacity -= min.opacity;
-          result.index   -= min.index;
-          break;
-        case ThickenMorphology:
-          /* Union with original image (maximize) - or should this be + */
-          Maximize( result.red,     min.red );
-          Maximize( result.green,   min.green );
-          Maximize( result.blue,    min.blue );
-          Maximize( result.opacity, min.opacity );
-          Maximize( result.index,   min.index );
-          break;
         default:
           /* result directly calculated or assigned */
           break;
       }
       /* Assign the resulting pixel values - Clamping Result */
       switch ( method ) {
-        case UndefinedMorphology:
-        case DilateIntensityMorphology:
-        case ErodeIntensityMorphology:
-          break;  /* full pixel was directly assigned - not a channel method */
+        case VoronoiMorphology:
+          SetPixelPacket(image,&result,q,q_indexes);
+          break;
         default:
           if ((channel & RedChannel) != 0)
-            q->red = ClampToQuantum(result.red);
+            SetRedPixelComponent(q,ClampToQuantum(result.red));
           if ((channel & GreenChannel) != 0)
-            q->green = ClampToQuantum(result.green);
+            SetGreenPixelComponent(q,ClampToQuantum(result.green));
           if ((channel & BlueChannel) != 0)
-            q->blue = ClampToQuantum(result.blue);
-          if ((channel & OpacityChannel) != 0
-              && image->matte == MagickTrue )
-            q->opacity = ClampToQuantum(QuantumRange-result.opacity);
+            SetBluePixelComponent(q,ClampToQuantum(result.blue));
+          if ((channel & OpacityChannel) != 0 && image->matte == MagickTrue )
+            SetOpacityPixelComponent(q,ClampToQuantum(QuantumRange-result.opacity));
           if ((channel & IndexChannel) != 0
               && image->colorspace == CMYKColorspace)
-            q_indexes[x] = ClampToQuantum(result.index);
+            SetIndexPixelComponent(q_indexes+x,ClampToQuantum(result.index));
           break;
       }
       /* Count up changed pixels */
-      if (   ( p[r].red != q->red )
-          || ( p[r].green != q->green )
-          || ( p[r].blue != q->blue )
-          || ( p[r].opacity != q->opacity )
+      if (   ( p[r].red != GetRedPixelComponent(q) )
+          || ( p[r].green != GetGreenPixelComponent(q) )
+          || ( p[r].blue != GetBluePixelComponent(q) )
+          || ( p[r].opacity != GetOpacityPixelComponent(q) )
           || ( image->colorspace == CMYKColorspace &&
-                  p_indexes[r] != q_indexes[x] ) )
-        changed++;  /* The pixel had some value changed! */
-      p++;
-      q++;
+               GetIndexPixelComponent(p_indexes+r) != GetIndexPixelComponent(q_indexes+x) ) )
+        changed++;  /* The pixel was changed in some way! */
+
+      p--; /* go backward through pixel buffers */
+      q--;
     } /* x */
-    sync=SyncCacheViewAuthenticPixels(q_view,exception);
-    if (sync == MagickFalse)
+    if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
       status=MagickFalse;
     if (image->progress_monitor != (MagickProgressMonitor) NULL)
-      {
-        MagickBooleanType
-          proceed;
+      if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
+                == MagickFalse )
+        status=MagickFalse;
 
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp critical (MagickCore_MorphologyImage)
-#endif
-        proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
-        if (proceed == MagickFalse)
-          status=MagickFalse;
-      }
   } /* y */
-  result_image->type=image->type;
-  q_view=DestroyCacheView(q_view);
-  p_view=DestroyCacheView(p_view);
-  return(status ? (unsigned long) changed : 0);
-}
 
+  auth_view=DestroyCacheView(auth_view);
+  virt_view=DestroyCacheView(virt_view);
+  return(status ? (ssize_t) changed : -1);
+}
 
+/* Apply a Morphology by calling theabove low level primitive application
+** functions.  This function handles any iteration loops, composition or
+** re-iteration of results, and compound morphology methods that is based
+** on multiple low-level (staged) morphology methods.
+**
+** Basically this provides the complex grue between the requested morphology
+** method and raw low-level implementation (above).
+*/
 MagickExport Image *MorphologyApply(const Image *image, const ChannelType
-     channel,const MorphologyMethod method, const long iterations,
+     channel,const MorphologyMethod method, const ssize_t iterations,
      const KernelInfo *kernel, const CompositeOperator compose,
      const double bias, ExceptionInfo *exception)
 {
+  CompositeOperator
+    curr_compose;
+
   Image
     *curr_image,    /* Image we are working with or iterating */
-    *work_image,    /* secondary image for primative iteration */
+    *work_image,    /* secondary image for primitive iteration */
     *save_image,    /* saved image - for 'edge' method only */
     *rslt_image;    /* resultant image - after multi-kernel handling */
 
@@ -2641,27 +3745,30 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
     *this_kernel;      /* the kernel being applied */
 
   MorphologyMethod
-    primative;      /* the current morphology primative being applied */
+    primitive;      /* the current morphology primitive being applied */
 
   CompositeOperator
     rslt_compose;   /* multi-kernel compose method for results to use */
 
   MagickBooleanType
+    special,        /* do we use a direct modify function? */
     verbose;        /* verbose output of results */
 
-  unsigned long
-    method_loop,    /* Loop 1: number of compound method iterations */
+  size_t
+    method_loop,    /* Loop 1: number of compound method iterations (norm 1) */
     method_limit,   /*         maximum number of compound method iterations */
     kernel_number,  /* Loop 2: the kernel number being applied */
-    stage_loop,     /* Loop 3: primative loop for compound morphology */
-    stage_limit,    /*         how many primatives in this compound */
-    kernel_loop,    /* Loop 4: iterate the kernel (basic morphology) */
+    stage_loop,     /* Loop 3: primitive loop for compound morphology */
+    stage_limit,    /*         how many primitives are in this compound */
+    kernel_loop,    /* Loop 4: iterate the kernel over image */
     kernel_limit,   /*         number of times to iterate kernel */
-    count,          /* total count of primative steps applied */
-    changed,        /* number pixels changed by last primative operation */
+    count,          /* total count of primitive steps applied */
     kernel_changed, /* total count of changed using iterated kernel */
     method_changed; /* total count of changed over method iteration */
 
+  ssize_t
+    changed;        /* number pixels changed by last primitive operation */
+
   char
     v_info[80];
 
@@ -2672,35 +3779,37 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
   assert(exception != (ExceptionInfo *) NULL);
   assert(exception->signature == MagickSignature);
 
-  count = 0;      /* number of low-level morphology primatives performed */
+  count = 0;      /* number of low-level morphology primitives performed */
   if ( iterations == 0 )
     return((Image *)NULL);   /* null operation - nothing to do! */
 
-  kernel_limit = (unsigned long) iterations;
+  kernel_limit = (size_t) iterations;
   if ( iterations < 0 )  /* negative interations = infinite (well alomst) */
-     kernel_limit = image->columns > image->rows ? image->columns : image->rows;
+     kernel_limit = image->columns>image->rows ? image->columns : image->rows;
 
-  verbose = ( GetImageArtifact(image,"verbose") != (const char *) NULL ) ?
-    MagickTrue : MagickFalse;
+  verbose = IsMagickTrue(GetImageArtifact(image,"verbose"));
 
   /* initialise for cleanup */
   curr_image = (Image *) image;
+  curr_compose = image->compose;
+  (void) curr_compose;
   work_image = save_image = rslt_image = (Image *) NULL;
   reflected_kernel = (KernelInfo *) NULL;
 
   /* Initialize specific methods
    * + which loop should use the given iteratations
-   * + how many primatives make up the compound morphology
+   * + how many primitives make up the compound morphology
    * + multi-kernel compose method to use (by default)
    */
   method_limit = 1;       /* just do method once, unless otherwise set */
-  stage_limit = 1;        /* assume method is not a compount */
+  stage_limit = 1;        /* assume method is not a compound */
+  special = MagickFalse;   /* assume it is NOT a direct modify primitive */
   rslt_compose = compose; /* and we are composing multi-kernels as given */
   switch( method ) {
-    case SmoothMorphology:  /* 4 primative compound morphology */
+    case SmoothMorphology:  /* 4 primitive compound morphology */
       stage_limit = 4;
       break;
-    case OpenMorphology:    /* 2 primative compound morphology */
+    case OpenMorphology:    /* 2 primitive compound morphology */
     case OpenIntensityMorphology:
     case TopHatMorphology:
     case CloseMorphology:
@@ -2710,27 +3819,63 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
       stage_limit = 2;
       break;
     case HitAndMissMorphology:
-      kernel_limit = 1;          /* no method or kernel iteration */
       rslt_compose = LightenCompositeOp;  /* Union of multi-kernel results */
-      break;
+      /* FALL THUR */
     case ThinningMorphology:
     case ThickenMorphology:
-    case DistanceMorphology:
-      method_limit = kernel_limit;  /* iterate method with each kernel */
+      method_limit = kernel_limit;  /* iterate the whole method */
       kernel_limit = 1;             /* do not do kernel iteration  */
-      rslt_compose = NoCompositeOp; /* Re-iterate with multiple kernels */
+      break;
+    case DistanceMorphology:
+    case VoronoiMorphology:
+      special = MagickTrue;
       break;
     default:
       break;
   }
 
+  /* Apply special methods with special requirments
+  ** For example, single run only, or post-processing requirements
+  */
+  if ( special == MagickTrue )
+    {
+      rslt_image=CloneImage(image,0,0,MagickTrue,exception);
+      if (rslt_image == (Image *) NULL)
+        goto error_cleanup;
+      if (SetImageStorageClass(rslt_image,DirectClass) == MagickFalse)
+        {
+          InheritException(exception,&rslt_image->exception);
+          goto error_cleanup;
+        }
+
+      changed = MorphologyPrimitiveDirect(rslt_image, method,
+                      channel, kernel, exception);
+
+      if ( verbose == MagickTrue )
+        (void) fprintf(stderr, "%s:%.20g.%.20g #%.20g => Changed %.20g\n",
+            CommandOptionToMnemonic(MagickMorphologyOptions, method),
+            1.0,0.0,1.0, (double) changed);
+
+      if ( changed < 0 )
+        goto error_cleanup;
+
+      if ( method == VoronoiMorphology ) {
+        /* Preserve the alpha channel of input image - but turned off */
+        (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel);
+        (void) CompositeImageChannel(rslt_image, DefaultChannels,
+          CopyOpacityCompositeOp, image, 0, 0);
+        (void) SetImageAlphaChannel(rslt_image, DeactivateAlphaChannel);
+      }
+      goto exit_cleanup;
+    }
+
   /* Handle user (caller) specified multi-kernel composition method */
   if ( compose != UndefinedCompositeOp )
     rslt_compose = compose;  /* override default composition for method */
   if ( rslt_compose == UndefinedCompositeOp )
     rslt_compose = NoCompositeOp; /* still not defined! Then re-iterate */
 
-  /* Some methods require a reflected kernel to use with primatives.
+  /* Some methods require a reflected kernel to use with primitives.
    * Create the reflected kernel for those methods. */
   switch ( method ) {
     case CorrelateMorphology:
@@ -2747,6 +3892,9 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
       break;
   }
 
+  /* Loops around more primitive morpholgy methods
+  **  erose, dilate, open, close, smooth, edge, etc...
+  */
   /* Loop 1:  iterate the compound method */
   method_loop = 0;
   method_changed = 1;
@@ -2756,7 +3904,9 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
 
     /* Loop 2:  iterate over each kernel in a multi-kernel list */
     norm_kernel = (KernelInfo *) kernel;
+    this_kernel = (KernelInfo *) kernel;
     rflt_kernel = reflected_kernel;
+
     kernel_number = 0;
     while ( norm_kernel != NULL ) {
 
@@ -2765,65 +3915,66 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
       while ( stage_loop < stage_limit ) {
         stage_loop++;   /* The stage of the compound morphology */
 
-        /* Select primative morphology for this stage of compound method */
+        /* Select primitive morphology for this stage of compound method */
         this_kernel = norm_kernel; /* default use unreflected kernel */
-        primative = method;        /* Assume method is a primative */
+        primitive = method;        /* Assume method is a primitive */
         switch( method ) {
           case ErodeMorphology:      /* just erode */
           case EdgeInMorphology:     /* erode and image difference */
-            primative = ErodeMorphology;
+            primitive = ErodeMorphology;
             break;
           case DilateMorphology:     /* just dilate */
           case EdgeOutMorphology:    /* dilate and image difference */
-            primative = DilateMorphology;
+            primitive = DilateMorphology;
             break;
           case OpenMorphology:       /* erode then dialate */
           case TopHatMorphology:     /* open and image difference */
-            primative = ErodeMorphology;
+            primitive = ErodeMorphology;
             if ( stage_loop == 2 )
-              primative = DilateMorphology;
+              primitive = DilateMorphology;
             break;
           case OpenIntensityMorphology:
-            primative = ErodeIntensityMorphology;
+            primitive = ErodeIntensityMorphology;
             if ( stage_loop == 2 )
-              primative = DilateIntensityMorphology;
+              primitive = DilateIntensityMorphology;
+            break;
           case CloseMorphology:      /* dilate, then erode */
           case BottomHatMorphology:  /* close and image difference */
             this_kernel = rflt_kernel; /* use the reflected kernel */
-            primative = DilateMorphology;
+            primitive = DilateMorphology;
             if ( stage_loop == 2 )
-              primative = ErodeMorphology;
+              primitive = ErodeMorphology;
             break;
           case CloseIntensityMorphology:
             this_kernel = rflt_kernel; /* use the reflected kernel */
-            primative = DilateIntensityMorphology;
+            primitive = DilateIntensityMorphology;
             if ( stage_loop == 2 )
-              primative = ErodeIntensityMorphology;
+              primitive = ErodeIntensityMorphology;
             break;
           case SmoothMorphology:         /* open, close */
             switch ( stage_loop ) {
               case 1: /* start an open method, which starts with Erode */
-                primative = ErodeMorphology;
+                primitive = ErodeMorphology;
                 break;
               case 2:  /* now Dilate the Erode */
-                primative = DilateMorphology;
+                primitive = DilateMorphology;
                 break;
               case 3:  /* Reflect kernel a close */
                 this_kernel = rflt_kernel; /* use the reflected kernel */
-                primative = DilateMorphology;
+                primitive = DilateMorphology;
                 break;
               case 4:  /* Finish the Close */
                 this_kernel = rflt_kernel; /* use the reflected kernel */
-                primative = ErodeMorphology;
+                primitive = ErodeMorphology;
                 break;
             }
             break;
           case EdgeMorphology:        /* dilate and erode difference */
-            primative = DilateMorphology;
+            primitive = DilateMorphology;
             if ( stage_loop == 2 ) {
               save_image = curr_image;      /* save the image difference */
               curr_image = (Image *) image;
-              primative = ErodeMorphology;
+              primitive = ErodeMorphology;
             }
             break;
           case CorrelateMorphology:
@@ -2837,34 +3988,35 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
             ** default).
             */
             this_kernel = rflt_kernel; /* use the reflected kernel */
-            primative = ConvolveMorphology;
+            primitive = ConvolveMorphology;
             break;
           default:
             break;
         }
+        assert( this_kernel != (KernelInfo *) NULL );
 
         /* Extra information for debugging compound operations */
         if ( verbose == MagickTrue ) {
           if ( stage_limit > 1 )
-            (void) FormatMagickString(v_info, MaxTextExtent, "%s:%lu.%lu -> ",
-                 MagickOptionToMnemonic(MagickMorphologyOptions, method),
-                 method_loop, stage_loop );
-          else if ( primative != method )
-            (void) FormatMagickString(v_info, MaxTextExtent, "%s:%lu -> ",
-                 MagickOptionToMnemonic(MagickMorphologyOptions, method),
-                 method_loop );
+            (void) FormatMagickString(v_info,MaxTextExtent,"%s:%.20g.%.20g -> ",
+             CommandOptionToMnemonic(MagickMorphologyOptions,method),(double)
+             method_loop,(double) stage_loop);
+          else if ( primitive != method )
+            (void) FormatMagickString(v_info, MaxTextExtent, "%s:%.20g -> ",
+              CommandOptionToMnemonic(MagickMorphologyOptions, method),(double)
+              method_loop);
           else
             v_info[0] = '\0';
         }
 
-        /* Loop 4: Iterate the kernel with primative */
+        /* Loop 4: Iterate the kernel with primitive */
         kernel_loop = 0;
         kernel_changed = 0;
         changed = 1;
         while ( kernel_loop < kernel_limit && changed > 0 ) {
           kernel_loop++;     /* the iteration of this kernel */
 
-          /* Create a destination image, if not yet defined */
+          /* Create a clone as the destination image, if not yet defined */
           if ( work_image == (Image *) NULL )
             {
               work_image=CloneImage(image,0,0,MagickTrue,exception);
@@ -2875,23 +4027,28 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
                   InheritException(exception,&work_image->exception);
                   goto error_cleanup;
                 }
+              /* work_image->type=image->type; ??? */
             }
 
           /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
           count++;
-          changed = MorphologyPrimitive(curr_image, work_image, primative,
-                        channel, this_kernel, bias, exception);
-          kernel_changed += changed;
-          method_changed += changed;
+          changed = MorphologyPrimitive(curr_image, work_image, primitive,
+                       channel, this_kernel, bias, exception);
 
           if ( verbose == MagickTrue ) {
             if ( kernel_loop > 1 )
               fprintf(stderr, "\n"); /* add end-of-line from previous */
-            fprintf(stderr, "%s%s%s:%lu.%lu #%lu => Changed %lu", v_info,
-                MagickOptionToMnemonic(MagickMorphologyOptions, primative),
-                 ( this_kernel == rflt_kernel ) ? "*" : "",
-               method_loop+kernel_loop-1, kernel_number, count, changed);
+            (void) fprintf(stderr, "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
+              v_info,CommandOptionToMnemonic(MagickMorphologyOptions,
+              primitive),(this_kernel == rflt_kernel ) ? "*" : "",
+              (double) (method_loop+kernel_loop-1),(double) kernel_number,
+              (double) count,(double) changed);
           }
+          if ( changed < 0 )
+            goto error_cleanup;
+          kernel_changed += changed;
+          method_changed += changed;
+
           /* prepare next loop */
           { Image *tmp = work_image;   /* swap images for iteration */
             work_image = curr_image;
@@ -2900,10 +4057,10 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
           if ( work_image == image )
             work_image = (Image *) NULL; /* replace input 'image' */
 
-        } /* End Loop 4: Iterate the kernel with primative */
+        } /* End Loop 4: Iterate the kernel with primitive */
 
-        if ( verbose == MagickTrue && kernel_changed != changed )
-          fprintf(stderr, "   Total %lu", kernel_changed);
+        if ( verbose == MagickTrue && kernel_changed != (size_t)changed )
+          fprintf(stderr, "   Total %.20g",(double) kernel_changed);
         if ( verbose == MagickTrue && stage_loop < stage_limit )
           fprintf(stderr, "\n"); /* add end-of-line before looping */
 
@@ -2931,7 +4088,7 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
         case BottomHatMorphology:
           if ( verbose == MagickTrue )
             fprintf(stderr, "\n%s: Difference with original image",
-                 MagickOptionToMnemonic(MagickMorphologyOptions, method) );
+                 CommandOptionToMnemonic(MagickMorphologyOptions, method) );
           (void) CompositeImageChannel(curr_image,
                   (ChannelType) (channel & ~SyncChannels),
                   DifferenceCompositeOp, image, 0, 0);
@@ -2939,7 +4096,7 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
         case EdgeMorphology:
           if ( verbose == MagickTrue )
             fprintf(stderr, "\n%s: Difference of Dilate and Erode",
-                 MagickOptionToMnemonic(MagickMorphologyOptions, method) );
+                 CommandOptionToMnemonic(MagickMorphologyOptions, method) );
           (void) CompositeImageChannel(curr_image,
                   (ChannelType) (channel & ~SyncChannels),
                   DifferenceCompositeOp, save_image, 0, 0);
@@ -2968,19 +4125,20 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
           curr_image = (Image *) image;  /* continue with original image */
         }
       else
-        { /* add the new 'current' result to the composition
+        { /* Add the new 'current' result to the composition
           **
           ** The removal of any 'Sync' channel flag in the Image Compositon
           ** below ensures the methematical compose method is applied in a
           ** purely mathematical way, and only to the selected channels.
-          ** Turn off SVG composition 'alpha blending'.
+          ** IE: Turn off SVG composition 'alpha blending'.
           */
           if ( verbose == MagickTrue )
             fprintf(stderr, " (compose \"%s\")",
-                 MagickOptionToMnemonic(MagickComposeOptions, rslt_compose) );
+                 CommandOptionToMnemonic(MagickComposeOptions, rslt_compose) );
           (void) CompositeImageChannel(rslt_image,
                (ChannelType) (channel & ~SyncChannels), rslt_compose,
                curr_image, 0, 0);
+          curr_image = DestroyImage(curr_image);
           curr_image = (Image *) image;  /* continue with original image */
         }
       if ( verbose == MagickTrue )
@@ -2999,16 +4157,14 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
 
   /* Yes goto's are bad, but it makes cleanup lot more efficient */
 error_cleanup:
-  if ( curr_image != (Image *) NULL &&
-       curr_image != rslt_image &&
-       curr_image != image )
-    curr_image = DestroyImage(curr_image);
+  if ( curr_image == rslt_image )
+    curr_image = (Image *) NULL;
   if ( rslt_image != (Image *) NULL )
     rslt_image = DestroyImage(rslt_image);
 exit_cleanup:
-  if ( curr_image != (Image *) NULL &&
-       curr_image != rslt_image &&
-       curr_image != image )
+  if ( curr_image == rslt_image || curr_image == image )
+    curr_image = (Image *) NULL;
+  if ( curr_image != (Image *) NULL )
     curr_image = DestroyImage(curr_image);
   if ( work_image != (Image *) NULL )
     work_image = DestroyImage(work_image);
@@ -3018,6 +4174,7 @@ exit_cleanup:
     reflected_kernel = DestroyKernelInfo(reflected_kernel);
   return(rslt_image);
 }
+
 \f
 /*
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -3045,10 +4202,10 @@ exit_cleanup:
 %  The format of the MorphologyImage method is:
 %
 %      Image *MorphologyImage(const Image *image,MorphologyMethod method,
-%        const long iterations,KernelInfo *kernel,ExceptionInfo *exception)
+%        const ssize_t iterations,KernelInfo *kernel,ExceptionInfo *exception)
 %
 %      Image *MorphologyImageChannel(const Image *image, const ChannelType
-%        channel,MorphologyMethod method,const long iterations,
+%        channel,MorphologyMethod method,const ssize_t iterations,
 %        KernelInfo *kernel,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
@@ -3073,11 +4230,8 @@ exit_cleanup:
 
 MagickExport Image *MorphologyImageChannel(const Image *image,
   const ChannelType channel,const MorphologyMethod method,
-  const long iterations,const KernelInfo *kernel,ExceptionInfo *exception)
+  const ssize_t iterations,const KernelInfo *kernel,ExceptionInfo *exception)
 {
-  const char
-    *artifact;
-
   KernelInfo
     *curr_kernel;
 
@@ -3095,8 +4249,10 @@ MagickExport Image *MorphologyImageChannel(const Image *image,
   curr_kernel = (KernelInfo *) kernel;
   if ( method == ConvolveMorphology ||  method == CorrelateMorphology )
     {
+      const char
+        *artifact;
       artifact = GetImageArtifact(image,"convolve:scale");
-      if ( artifact != (char *)NULL ) {
+      if ( artifact != (const char *)NULL ) {
         if ( curr_kernel == kernel )
           curr_kernel = CloneKernelInfo(kernel);
         if (curr_kernel == (KernelInfo *) NULL) {
@@ -3108,25 +4264,25 @@ MagickExport Image *MorphologyImageChannel(const Image *image,
     }
 
   /* display the (normalized) kernel via stderr */
-  artifact = GetImageArtifact(image,"showkernel");
-  if ( artifact == (const char *) NULL)
-    artifact = GetImageArtifact(image,"convolve:showkernel");
-  if ( artifact == (const char *) NULL)
-    artifact = GetImageArtifact(image,"morphology:showkernel");
-  if ( artifact != (const char *) NULL)
+  if ( IsMagickTrue(GetImageArtifact(image,"showkernel"))
+    || IsMagickTrue(GetImageArtifact(image,"convolve:showkernel"))
+    || IsMagickTrue(GetImageArtifact(image,"morphology:showkernel")) )
     ShowKernelInfo(curr_kernel);
 
-  /* override the default handling of multi-kernel morphology results
-   * if 'Undefined' use the default method
-   * if 'None' (default for 'Convolve') re-iterate previous result
-   * otherwise merge resulting images using compose method given
+  /* Override the default handling of multi-kernel morphology results
+   * If 'Undefined' use the default method
+   * If 'None' (default for 'Convolve') re-iterate previous result
+   * Otherwise merge resulting images using compose method given.
+   * Default for 'HitAndMiss' is 'Lighten'.
    */
-  compose = UndefinedCompositeOp;  /* use default for method */
-  artifact = GetImageArtifact(image,"morphology:compose");
-  if ( artifact != (const char *) NULL)
-    compose = (CompositeOperator) ParseMagickOption(
+  { const char
+      *artifact;
+    artifact = GetImageArtifact(image,"morphology:compose");
+    compose = UndefinedCompositeOp;  /* use default for method */
+    if ( artifact != (const char *) NULL)
+      compose = (CompositeOperator) ParseCommandOption(
                              MagickComposeOptions,MagickFalse,artifact);
-
+  }
   /* Apply the Morphology */
   morphology_image = MorphologyApply(image, channel, method, iterations,
                          curr_kernel, compose, image->bias, exception);
@@ -3138,7 +4294,7 @@ MagickExport Image *MorphologyImageChannel(const Image *image,
 }
 
 MagickExport Image *MorphologyImage(const Image *image, const MorphologyMethod
-  method, const long iterations,const KernelInfo *kernel, ExceptionInfo
+  method, const ssize_t iterations,const KernelInfo *kernel, ExceptionInfo
   *exception)
 {
   Image
@@ -3202,12 +4358,13 @@ static void RotateKernelInfo(KernelInfo *kernel, double angle)
   switch (kernel->type) {
     /* These built-in kernels are cylindrical kernels, rotating is useless */
     case GaussianKernel:
-    case DOGKernel:
+    case DoGKernel:
+    case LoGKernel:
     case DiskKernel:
     case PeaksKernel:
     case LaplacianKernel:
     case ChebyshevKernel:
-    case ManhattenKernel:
+    case ManhattanKernel:
     case EuclideanKernel:
       return;
 
@@ -3248,15 +4405,15 @@ static void RotateKernelInfo(KernelInfo *kernel, double angle)
           kernel->values[1] = t;
           /* rotate non-centered origin */
           if ( kernel->x != 1 || kernel->y != 1 ) {
-            long x,y;
-            x = (long) kernel->x-1;
-            y = (long) kernel->y-1;
+            ssize_t x,y;
+            x = (ssize_t) kernel->x-1;
+            y = (ssize_t) kernel->y-1;
                  if ( x == y  ) x = 0;
             else if ( x == 0  ) x = -y;
             else if ( x == -y ) y = 0;
             else if ( y == 0  ) y = x;
-            kernel->x = (unsigned long) x+1;
-            kernel->y = (unsigned long) y+1;
+            kernel->x = (ssize_t) x+1;
+            kernel->y = (ssize_t) y+1;
           }
           angle = fmod(angle+315.0, 360.0);  /* angle reduced 45 degrees */
           kernel->angle = fmod(kernel->angle+45.0, 360.0);
@@ -3267,14 +4424,14 @@ static void RotateKernelInfo(KernelInfo *kernel, double angle)
   if ( 45.0 < fmod(angle, 180.0)  && fmod(angle,180.0) <= 135.0 )
     {
       if ( kernel->width == 1 || kernel->height == 1 )
-        { /* Do a transpose of the image, which results in a 90
-          ** degree rotation of a 1 dimentional kernel
+        { /* Do a transpose of a 1 dimentional kernel,
+          ** which results in a fast 90 degree rotation of some type.
           */
-          long
+          ssize_t
             t;
-          t = (long) kernel->width;
+          t = (ssize_t) kernel->width;
           kernel->width = kernel->height;
-          kernel->height = (unsigned long) t;
+          kernel->height = (size_t) t;
           t = kernel->x;
           kernel->x = kernel->y;
           kernel->y = t;
@@ -3288,7 +4445,7 @@ static void RotateKernelInfo(KernelInfo *kernel, double angle)
         }
       else if ( kernel->width == kernel->height )
         { /* Rotate a square array of values by 90 degrees */
-          { register unsigned long
+          { register size_t
               i,j,x,y;
             register MagickRealType
               *k,t;
@@ -3303,11 +4460,11 @@ static void RotateKernelInfo(KernelInfo *kernel, double angle)
                 }
           }
           /* rotate the origin - relative to center of array */
-          { register long x,y;
-            x = (long) kernel->x*2-kernel->width+1;
-            y = (long) kernel->y*2-kernel->height+1;
-            kernel->x = (unsigned long) ( -y +kernel->width-1)/2;
-            kernel->y = (unsigned long) ( +x +kernel->height-1)/2;
+          { register ssize_t x,y;
+            x = (ssize_t) (kernel->x*2-kernel->width+1);
+            y = (ssize_t) (kernel->y*2-kernel->height+1);
+            kernel->x = (ssize_t) ( -y +(ssize_t) kernel->width-1)/2;
+            kernel->y = (ssize_t) ( +x +(ssize_t) kernel->height-1)/2;
           }
           angle = fmod(angle+270.0, 360.0);     /* angle reduced 90 degrees */
           kernel->angle = fmod(kernel->angle+90.0, 360.0);
@@ -3322,7 +4479,7 @@ static void RotateKernelInfo(KernelInfo *kernel, double angle)
        * Basically all that is needed is a reversal of the kernel data!
        * And a reflection of the origon
        */
-      unsigned long
+      size_t
         i,j;
       register double
         *k,t;
@@ -3331,8 +4488,8 @@ static void RotateKernelInfo(KernelInfo *kernel, double angle)
       for ( i=0, j=kernel->width*kernel->height-1;  i<j;  i++, j--)
         t=k[i],  k[i]=k[j],  k[j]=t;
 
-      kernel->x = (long) kernel->width  - kernel->x - 1;
-      kernel->y = (long) kernel->height - kernel->y - 1;
+      kernel->x = (ssize_t) kernel->width  - kernel->x - 1;
+      kernel->y = (ssize_t) kernel->height - kernel->y - 1;
       angle = fmod(angle-180.0, 360.0);   /* angle+180 degrees */
       kernel->angle = fmod(kernel->angle+180.0, 360.0);
     }
@@ -3341,24 +4498,6 @@ static void RotateKernelInfo(KernelInfo *kernel, double angle)
    * performed here, posibily with a linear kernel restriction.
    */
 
-#if 0
-    { /* Do a Flop by reversing each row.
-       */
-      unsigned long
-        y;
-      register long
-        x,r;
-      register double
-        *k,t;
-
-      for ( y=0, k=kernel->values; y < kernel->height; y++, k+=kernel->width)
-        for ( x=0, r=kernel->width-1; x<kernel->width/2; x++, r--)
-          t=k[x],  k[x]=k[r],  k[r]=t;
-
-      kernel->x = kernel->width - kernel->x - 1;
-      angle = fmod(angle+180.0, 360.0);
-    }
-#endif
   return;
 }
 \f
@@ -3382,10 +4521,10 @@ static void RotateKernelInfo(KernelInfo *kernel, double angle)
 %  is then passed to UnityAddKernelInfo() to add a scled unity kernel
 %  into the scaled/normalized kernel.
 %
-%  The format of the ScaleKernelInfo method is:
+%  The format of the ScaleGeometryKernelInfo method is:
 %
-%      void ScaleKernelInfo(KernelInfo *kernel, const double scaling_factor,
-%               const MagickStatusType normalize_flags )
+%      void ScaleGeometryKernelInfo(KernelInfo *kernel,
+%        const double scaling_factor,const MagickStatusType normalize_flags)
 %
 %  A description of each parameter follows:
 %
@@ -3471,7 +4610,7 @@ MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
 %  For special kernels designed for locating shapes using 'Correlate', (often
 %  only containing +1 and -1 values, representing foreground/brackground
 %  matching) a special normalization method is provided to scale the positive
-%  values seperatally to those of the negative values, so the kernel will be
+%  values separately to those of the negative values, so the kernel will be
 %  forced to become a zero-sum kernel better suited to such searches.
 %
 %  WARNING: Correct normalization of the kernel assumes that the '*_range'
@@ -3504,7 +4643,7 @@ MagickExport void ScaleGeometryKernelInfo (KernelInfo *kernel,
 MagickExport void ScaleKernelInfo(KernelInfo *kernel,
   const double scaling_factor,const GeometryFlags normalize_flags)
 {
-  register long
+  register ssize_t
     i;
 
   register double
@@ -3539,7 +4678,7 @@ MagickExport void ScaleKernelInfo(KernelInfo *kernel,
   pos_scale = scaling_factor/pos_scale;
   neg_scale = scaling_factor/neg_scale;
 
-  for (i=0; i < (long) (kernel->width*kernel->height); i++)
+  for (i=0; i < (ssize_t) (kernel->width*kernel->height); i++)
     if ( ! IsNan(kernel->values[i]) )
       kernel->values[i] *= (kernel->values[i] >= 0) ? pos_scale : neg_scale;
 
@@ -3592,21 +4731,20 @@ MagickExport void ShowKernelInfo(KernelInfo *kernel)
   KernelInfo
     *k;
 
-  unsigned long
+  size_t
     c, i, u, v;
 
   for (c=0, k=kernel;  k != (KernelInfo *) NULL;  c++, k=k->next ) {
 
     fprintf(stderr, "Kernel");
     if ( kernel->next != (KernelInfo *) NULL )
-      fprintf(stderr, " #%lu", c );
+      fprintf(stderr, " #%lu", (unsigned long) c );
     fprintf(stderr, " \"%s",
-          MagickOptionToMnemonic(MagickKernelOptions, k->type) );
+          CommandOptionToMnemonic(MagickKernelOptions, k->type) );
     if ( fabs(k->angle) > MagickEpsilon )
       fprintf(stderr, "@%lg", k->angle);
-    fprintf(stderr, "\" of size %lux%lu%+ld%+ld",
-          k->width, k->height,
-          k->x, k->y );
+    fprintf(stderr, "\" of size %lux%lu%+ld%+ld",(unsigned long) k->width,
+      (unsigned long) k->height,(long) k->x,(long) k->y);
     fprintf(stderr,
           " with values from %.*lg to %.*lg\n",
           GetMagickPrecision(), k->minimum,
@@ -3622,7 +4760,7 @@ MagickExport void ShowKernelInfo(KernelInfo *kernel)
       fprintf(stderr, " (Sum %.*lg)\n",
           GetMagickPrecision(), k->positive_range+k->negative_range);
     for (i=v=0; v < k->height; v++) {
-      fprintf(stderr, "%2lu:", v );
+      fprintf(stderr, "%2lu:", (unsigned long) v );
       for (u=0; u < k->width; u++, i++)
         if ( IsNan(k->values[i]) )
           fprintf(stderr," %*s", GetMagickPrecision()+3, "nan");
@@ -3651,13 +4789,8 @@ MagickExport void ShowKernelInfo(KernelInfo *kernel)
 %  value is usually provided by the user as a percentage value in the
 %  'convolve:scale' setting.
 %
-%  The resulting effect is to either convert a 'zero-summing' edge detection
-%  kernel (such as a "Laplacian", "DOG" or a "LOG") into a 'sharpening'
-%  kernel.
-%
-%  Alternativally by using a purely positive kernel, and using a negative
-%  post-normalizing scaling factor, you can convert a 'blurring' kernel (such
-%  as a "Gaussian") into a 'unsharp' kernel.
+%  The resulting effect is to convert the defined kernels into blended
+%  soft-blurs, unsharp kernels or into sharpening kernels.
 %
 %  The format of the UnityAdditionKernelInfo method is:
 %
@@ -3713,7 +4846,7 @@ MagickExport void UnityAddKernelInfo(KernelInfo *kernel,
 */
 MagickExport void ZeroKernelNans(KernelInfo *kernel)
 {
-  register unsigned long
+  register size_t
     i;
 
   /* do the other kernels in a multi-kernel list first */