]> granicus.if.org Git - imagemagick/blobdiff - magick/morphology.c
(no commit message)
[imagemagick] / magick / morphology.c
index ddf34053c1f7f88fd9a0c3eae2c20dc5dc3e72f3..982325a75b4a1348c5ee160a935c43b2a564bf35 100644 (file)
 #include "magick/string-private.h"
 #include "magick/token.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 +112,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 *, double),
+  ExpandMirrorKernelInfo(KernelInfo *),
+  ExpandRotateKernelInfo(KernelInfo *, const double),
   RotateKernelInfo(KernelInfo *, double);
 \f
 
@@ -144,7 +148,7 @@ static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
 %  anywhere within that array of values.
 %
 %  Previously IM was restricted to a square of odd size using the exact
-%  center as origin, this is no longer the case, and any rectangular kernel
+%  center as origin, this is no ssize_ter the case, and any rectangular kernel
 %  with any value being declared the origin. This in turn allows the use of
 %  highly asymmetrical kernels.
 %
@@ -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
@@ -186,9 +186,16 @@ static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
 %
 %     " 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
@@ -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 == ',')
@@ -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,19 +362,24 @@ 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
+  ssize_t
     type;
 
   const char
@@ -416,9 +428,9 @@ static KernelInfo *ParseKernelName(const char *kernel_string)
       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);
+        args.xi = (double)(((ssize_t)args.rho-1)/2);
       if ( (flags & YValue) == 0 )
-        args.psi = (double)(((long)args.sigma-1)/2);
+        args.psi = (double)(((ssize_t)args.sigma-1)/2);
       break;
     case SquareKernel:
     case DiamondKernel:
@@ -434,7 +446,7 @@ static KernelInfo *ParseKernelName(const char *kernel_string)
         args.xi = 1.0;
       break;
     case ChebyshevKernel:
-    case ManhattenKernel:
+    case ManhattanKernel:
     case EuclideanKernel:
       if ( (flags & HeightValue) == 0 )           /* no distance scale */
         args.sigma = 100.0;                       /* default distance scaling */
@@ -447,7 +459,19 @@ static KernelInfo *ParseKernelName(const char *kernel_string)
       break;
   }
 
-  return(AcquireKernelBuiltIn((KernelInfoType)type, &args));
+  kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args);
+
+  /* global expand to rotated kernel list - only for single kernels */
+  if ( kernel->next == (KernelInfo *) NULL ) {
+    if ( (flags & AreaValue) != 0 )         /* '@' symbol in kernel args */
+      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);
 }
 
 MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
@@ -463,7 +487,7 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
   const char
     *p;
 
-  unsigned long
+  size_t
     kernel_number;
 
   p = kernel_string;
@@ -483,7 +507,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);
@@ -557,19 +582,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
@@ -583,15 +608,6 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %       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,
@@ -627,39 +643,144 @@ 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
-%    FreiChen:{angle}
-%      Frei-Chen 'Edge' convolution kernel (3x3)
-%             -1,     0,    1
-%           -sqrt(2), 0,  sqrt(2)
-%             -1,     0,    1
+%          | -1, 0, 1 |
+%          | -2, 0,-2 |
+%          | -1, 0, 1 |
+%
+%    Sobel:{type},{angle}
+%      Type 0:  default un-nomalized version shown above.
+%
+%      Type 1:  As default but pre-normalized
+%          | 1, 0, -1 |
+%          | 2, 0, -2 |  / 4
+%          | 1, 0, -1 |
+%
+%      Type 2:  Diagonal version with same normalization as 1
+%          | 1, 0, -1 |
+%          | 2, 0, -2 |  / 4
+%          | 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 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     |
+%
+%      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 10: All 9 of the following pre-weighted kernels...
+%
+%        Type 11: |   1,     0,   -1     |
+%                 | sqrt(2), 0, -sqrt(2) | / 2*sqrt(2)
+%                 |   1,     0,   -1     |
+%
+%        Type 12: | 1, sqrt(2), 1 |
+%                 | 0,   0,     0 | / 2*sqrt(2)
+%                 | 1, sqrt(2), 1 |
+%
+%        Type 13: | sqrt(2), -1,    0     |
+%                 |  -1,      0,    1     | / 2*sqrt(2)
+%                 |   0,      1, -sqrt(2) |
+%
+%        Type 14: |    0,     1, -sqrt(2) |
+%                 |   -1,     0,     1    | / 2*sqrt(2)
+%                 | sqrt(2), -1,     0    |
+%
+%        Type 15: | 0, -1, 0 |
+%                 | 1,  0, 1 | / 2
+%                 | 0, -1, 0 |
+%
+%        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
 %
@@ -729,17 +850,35 @@ 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
-%    LineEnds
+%       Find 90 degree corners of a binary shape
+%    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
+%    Skeleton:type
+%       Traditional skeleton generating kernels.
+%         Type 1: Tradional Skeleton kernel (4 connected skeleton)
+%         Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
+%         Type 3: Experimental Variation to try to present left-right symmetry
+%         Type 4: Experimental Variation to preserve left-right symmetry
 %
 %  Distance Measuring Kernels
 %
@@ -759,8 +898,8 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %       '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
+%    Manhattan:[{radius}][x{scale}[%!]]
+%       Manhattan 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
@@ -781,7 +920,7 @@ MagickExport KernelInfo *AcquireKernelInfo(const char *kernel_string)
 %       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
+%       "Manhattan" 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
@@ -795,10 +934,10 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
   KernelInfo
     *kernel;
 
-  register long
+  register ssize_t
     i;
 
-  register long
+  register ssize_t
     u,
     v;
 
@@ -808,28 +947,32 @@ 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:
+    case TestKernel:
       break;
-    case LaplacianKernel:  /* Named Descrete Convolution Kernels */
+    case UnityKernel:      /* Named Descrete Convolution Kernels */
+    case LaplacianKernel:
     case SobelKernel:
     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 LineEndsKernel:
     case LineJunctionsKernel:
-    case ThinningKernel:
+    case RidgesKernel:
     case ConvexHullKernel:
     case SkeletonKernel:
-      /* A pre-generated kernel is not needed */
-      break;
+      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 GaussianKernel:
-    case DOGKernel:
+    case DoGKernel:
+    case LoGKernel:
     case BlurKernel:
-    case DOBKernel:
     case CometKernel:
     case DiamondKernel:
     case SquareKernel:
@@ -840,10 +983,11 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
     case RingKernel:
     case PeaksKernel:
     case ChebyshevKernel:
-    case ManhattenKernel:
+    case ManhattanKernel:
     case EuclideanKernel:
-#endif
+#else
     default:
+#endif
       /* Generate the base Kernel Structure */
       kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
       if (kernel == (KernelInfo *) NULL)
@@ -860,21 +1004,21 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
   switch(type) {
     /* Convolution Kernels */
     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)
@@ -887,13 +1031,13 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
          * 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++)
+                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 */
@@ -903,27 +1047,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++)
+                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++)
+                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;
                     }
@@ -954,20 +1098,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,
@@ -991,63 +1131,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);
-            = 1.0/(MagickSQ2PI*sigma );
+            alpha = 1.0/(2.0*sigma*sigma);
+            beta= 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
@@ -1066,7 +1179,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:
@@ -1077,7 +1190,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;
@@ -1100,7 +1213,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 */
@@ -1111,12 +1224,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);
@@ -1141,8 +1254,7 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
 
     /* Convolution Kernels - Well Known Constants */
     case LaplacianKernel:
-      {
-        switch ( (int) args->rho ) {
+      { switch ( (int) args->rho ) {
           case 0:
           default: /* laplacian square filter -- default */
             kernel=ParseKernelArray("3: -1,-1,-1  -1,8,-1  -1,-1,-1");
@@ -1164,14 +1276,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)
@@ -1180,28 +1292,52 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         break;
       }
     case SobelKernel:
-      {
-        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 */
+#if 0
+      { /* Sobel with optional 'sub-types' */
+        switch ( (int) args->rho ) {
+          default:
+          case 0:
+            kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            break;
+          case 1:
+            kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            ScaleKernelInfo(kernel, 0.25, NoValue);
+            break;
+          case 2:
+            kernel=ParseKernelArray("3: 1,2,0  2,0,-2  0,-2,-1");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            ScaleKernelInfo(kernel, 0.25, NoValue);
+            break;
+        }
+        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;
       }
-    case FreiChenKernel:
-      {
-        kernel=ParseKernelArray("3: -1,0,1  -2,0,2  -1,0,1");
+#else
+      { /* Simple Sobel Kernel */
+        kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
         if (kernel == (KernelInfo *) NULL)
           return(kernel);
-        kernel->values[3] = -MagickSQ2;
-        kernel->values[5] = +MagickSQ2;
-        CalcKernelMetaData(kernel);     /* recalculate meta-data */
-        RotateKernelInfo(kernel, args->rho);  /* Rotate by angle */
+        kernel->type = type;
+        RotateKernelInfo(kernel, args->rho);
         break;
       }
+#endif
     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;
@@ -1210,7 +1346,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;
@@ -1219,7 +1355,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;
@@ -1228,21 +1364,137 @@ 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;
         RotateKernelInfo(kernel, args->rho);
         break;
       }
+    case FreiChenKernel:
+      /* 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 0:
+            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);     /* recalculate meta-data */
+            break;
+          case 2:
+            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, 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);     /* recalculate meta-data */
+            ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
+            break;
+          case 12:
+            kernel=ParseKernelArray("3: 1,2,1  0,0,0  1,2,1");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            kernel->values[1] = +MagickSQ2;
+            kernel->values[7] = +MagickSQ2;
+            CalcKernelMetaData(kernel);
+            ScaleKernelInfo(kernel, 1.0/2.0*MagickSQ2, NoValue);
+            break;
+          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);
+            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, 1.0/2.0*MagickSQ2, NoValue);
+            break;
+          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 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 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 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 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;
+        }
+        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 */
     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));
@@ -1250,9 +1502,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;
@@ -1268,21 +1520,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,
@@ -1291,7 +1543,7 @@ 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 */
@@ -1300,12 +1552,14 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
       }
     case DiskKernel:
       {
-        long limit = (long)(args->rho*args->rho);
-        if (args->rho < 0.1)             /* default radius approx 3.5 */
+        ssize_t
+         limit = (ssize_t)(args->rho*args->rho);
+
+        if (args->rho < 0.4)           /* 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->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));
@@ -1313,8 +1567,8 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
           return(DestroyKernelInfo(kernel));
 
         /* 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++)
+        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
@@ -1327,8 +1581,8 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         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->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));
@@ -1336,8 +1590,8 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
           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++)
+        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);
@@ -1348,8 +1602,8 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         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->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));
@@ -1357,8 +1611,8 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
           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++)
+        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);
@@ -1368,38 +1622,38 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
     case RingKernel:
     case PeaksKernel:
       {
-        long
+        ssize_t
           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 = ((size_t)args->sigma)*2+1;
+            limit1 = (ssize_t)(args->rho*args->rho);
+            limit2 = (ssize_t)(args->sigma*args->sigma);
           }
         else
           {
-            kernel->width = ((unsigned long)args->rho)*2+1;
-            limit1 = (long)args->sigma*args->sigma;
-            limit2 = (long)args->rho*args->rho;
+            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 = (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)
           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;
+        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
@@ -1420,7 +1674,7 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         if (kernel == (KernelInfo *) NULL)
           return(kernel);
         kernel->type = type;
-        ExpandKernelInfo(kernel, 90.0); /* Create a list of 4 rotated kernels */
+        ExpandMirrorKernelInfo(kernel); /* mirror expansion of other kernels */
         break;
       }
     case CornersKernel:
@@ -1429,77 +1683,276 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
         if (kernel == (KernelInfo *) NULL)
           return(kernel);
         kernel->type = type;
-        ExpandKernelInfo(kernel, 90.0); /* Create a list of 4 rotated kernels */
+        ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
         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;
+      { /* Kernels for finding the end of thin lines */
+        switch ( (int) args->rho ) {
+          case 0:
+          default:
+            /* set of kernels to find all end of lines */
+            kernel=AcquireKernelInfo("LineEnds:1>;LineEnds:2>");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            break;
+          case 1:
+            /* kernel for 4-connected line ends - no rotation */
+            kernel=ParseKernelArray("3: 0,0,0  0,1,0  -,1,-");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            break;
+         case 2:
+            /* kernel to add for 8-connected lines - no rotation */
+            kernel=ParseKernelArray("3: 0,0,0  0,1,0  0,0,1");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            break;
+        }
         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 */
+            kernel=AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            break;
+          case 1:
+            /* Y Junction */
+            kernel=ParseKernelArray("3: 1,-,1  -,1,-  -,1,-");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            break;
+          case 2:
+            /* Diagonal T Junctions */
+            kernel=ParseKernelArray("3: 1,-,-  -,1,-  1,-,1");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            break;
+          case 3:
+            /* Orthogonal T Junctions */
+            kernel=ParseKernelArray("3: -,-,-  1,1,1  -,1,-");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            break;
+          case 4:
+            /* Diagonal X Junctions */
+            kernel=ParseKernelArray("3: 1,-,1  -,1,-  1,-,1");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            break;
+          case 5:
+            /* Orthogonal X Junctions - minimal diamond kernel */
+            kernel=ParseKernelArray("3: -,1,-  1,1,1  -,1,-");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            break;
+        }
+        break;
+      }
+    case RidgesKernel:
+      { /* Ridges - Ridge finding kernels */
         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;
+        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 4 kernels */
+        /* first set of 8 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");
+        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;
-        ExpandKernelInfo(new_kernel, 90.0);
+        ExpandRotateKernelInfo(new_kernel, 90.0);
         LastKernelInfo(kernel)->next = new_kernel;
         break;
       }
-    case ThinningKernel:
-      { /* Thinning Kernel ??  -- filled corner and edge */
-        kernel=ParseKernelArray("3: 0,0,-  0,1,1  -,1,1");
-        if (kernel == (KernelInfo *) NULL)
-          return(kernel);
-        kernel->type = type;
-        ExpandKernelInfo(kernel, 45);
-        break;
-      }
     case SkeletonKernel:
       {
-        kernel=AcquireKernelInfo("Edges;Corners");
+        KernelInfo
+          *new_kernel;
+        switch ( (int) args->rho ) {
+          case 1:
+          default:
+            /* Traditional Skeleton...
+            ** A cyclically rotated single kernel
+            */
+            kernel=ParseKernelArray("3: 0,0,0  -,1,-  1,1,1");
+            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=ParseKernelArray("3: 0,0,0  -,1,-  1,1,1");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            new_kernel=ParseKernelArray("3: -,0,0  1,1,0  -,1,-");
+            if (new_kernel == (KernelInfo *) NULL)
+              return(new_kernel);
+            new_kernel->type = type;
+            LastKernelInfo(kernel)->next = new_kernel;
+            ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
+            break;
+          case 3:
+            /* Jittered Skeleton: do top, then bottom, then each sides */
+            /* Do top edge */
+            kernel=ParseKernelArray("3: 0,0,0  -,1,-  1,1,1");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            new_kernel=ParseKernelArray("3: 0,0,-  0,1,1  -,1,-");
+            if (new_kernel == (KernelInfo *) NULL)
+              return(new_kernel);
+            new_kernel->type = type;
+            LastKernelInfo(kernel)->next = new_kernel;
+            new_kernel=ParseKernelArray("3: -,0,0  1,1,0  -,1,-");
+            if (new_kernel == (KernelInfo *) NULL)
+              return(new_kernel);
+            new_kernel->type = type;
+            LastKernelInfo(kernel)->next = new_kernel;
+            /* Do Bottom edge */
+            new_kernel=ParseKernelArray("3: 1,1,1  -,1,-  0,0,0");
+            if (new_kernel == (KernelInfo *) NULL)
+              return(new_kernel);
+            new_kernel->type = type;
+            LastKernelInfo(kernel)->next = new_kernel;
+            new_kernel=ParseKernelArray("3: -,1,-  1,1,0  -,0,0");
+            if (new_kernel == (KernelInfo *) NULL)
+              return(new_kernel);
+            new_kernel->type = type;
+            LastKernelInfo(kernel)->next = new_kernel;
+            new_kernel=ParseKernelArray("3: -,1,-  0,1,1  0,0,-");
+            if (new_kernel == (KernelInfo *) NULL)
+              return(new_kernel);
+            new_kernel->type = type;
+            LastKernelInfo(kernel)->next = new_kernel;
+            /* Last the two sides */
+            new_kernel=ParseKernelArray("3: 0,-,1  0,1,1  0,-,1");
+            if (new_kernel == (KernelInfo *) NULL)
+              return(new_kernel);
+            new_kernel->type = type;
+            LastKernelInfo(kernel)->next = new_kernel;
+            new_kernel=ParseKernelArray("3: 1,-,0  1,1,0  1,-,0");
+            if (new_kernel == (KernelInfo *) NULL)
+              return(new_kernel);
+            new_kernel->type = type;
+            LastKernelInfo(kernel)->next = new_kernel;
+            break;
+          case 4:
+            /* Just a simple 'Edge' kernel, but with a extra two kernels
+            ** to finish off diagonal lines,  top then bottom then sides.
+            ** Works well for test case but fails for general case.
+            */
+            kernel=ParseKernelArray("3: 0,0,0  -,1,-  1,1,1");
+            if (kernel == (KernelInfo *) NULL)
+              return(kernel);
+            kernel->type = type;
+            new_kernel=ParseKernelArray("3: 0,0,0  0,1,1  1,1,-");
+            if (new_kernel == (KernelInfo *) NULL)
+              return(DestroyKernelInfo(kernel));
+            new_kernel->type = type;
+            LastKernelInfo(kernel)->next = new_kernel;
+            new_kernel=ParseKernelArray("3: 0,0,0  1,1,0  -,1,1");
+            if (new_kernel == (KernelInfo *) NULL)
+              return(DestroyKernelInfo(kernel));
+            new_kernel->type = type;
+            LastKernelInfo(kernel)->next = new_kernel;
+            ExpandMirrorKernelInfo(kernel);
+            /* Append a set of corner kernels */
+            new_kernel=ParseKernelArray("3: 0,0,-  0,1,1  -,1,-");
+            if (new_kernel == (KernelInfo *) NULL)
+              return(DestroyKernelInfo(kernel));
+            new_kernel->type = type;
+            ExpandRotateKernelInfo(new_kernel, 90.0);
+            LastKernelInfo(kernel)->next = new_kernel;
+            break;
+        }
         break;
       }
     /* Distance Measuring Kernels */
@@ -1508,38 +1961,38 @@ 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)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));
 
-        for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-          for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
+        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(u)>labs(v)) ? labs(u) : labs(v)) );
+                 args->sigma*((labs((long) u)>labs((long) v)) ? labs((long) u) : labs((long) v)) );
         kernel->maximum = kernel->values[0];
         break;
       }
-    case ManhattenKernel:
+    case ManhattanKernel:
       {
         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));
         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++)
+        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(u)+labs(v)) );
+                 args->sigma*(labs((long) u)+labs((long) v)) );
         kernel->maximum = kernel->values[0];
         break;
       }
@@ -1548,16 +2001,16 @@ 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)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));
 
-        for ( i=0, v=-kernel->y; v <= (long)kernel->y; v++)
-          for ( u=-kernel->x; u <= (long)kernel->x; u++, i++)
+        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];
@@ -1566,8 +2019,8 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
     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");
+        /* Unity or No-Op Kernel - Basically just a single pixel on its own */
+        kernel=ParseKernelArray("1:1");
         if (kernel == (KernelInfo *) NULL)
           return(kernel);
         kernel->type = ( type == UnityKernel ) ? UnityKernel : UndefinedKernel;
@@ -1592,7 +2045,7 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
 %
 %  CloneKernelInfo() creates a new clone of the given Kernel List so that its
 %  can be modified without effecting the original.  The cloned kernel should
-%  be destroyed using DestoryKernelInfo() when no longer needed.
+%  be destroyed using DestoryKernelInfo() when no ssize_ter needed.
 %
 %  The format of the CloneKernelInfo method is:
 %
@@ -1605,7 +2058,7 @@ MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
 */
 MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
 {
-  register long
+  register ssize_t
     i;
 
   KernelInfo
@@ -1622,7 +2075,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 */
@@ -1658,7 +2111,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);
@@ -1676,21 +2128,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                             %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  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                             %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  ExpandKernelInfo() takes a single kernel, and expands it into a list
-%  of kernels each incrementally rotated the angle given.
+%  ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
+%  incrementally by the angle given, until the first 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:
 %
@@ -1702,22 +2229,53 @@ MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
 % especially with regard to non-orthogonal angles, and rotation of larger
 % 2D kernels.
 */
-static void ExpandKernelInfo(KernelInfo *kernel, double angle)
+
+/* Internal Routine - Return true if two kernels are the same */
+static MagickBooleanType SameKernelInfo(const KernelInfo *kernel1,
+     const KernelInfo *kernel2)
+{
+  register size_t
+    i;
+
+  /* check size and origin location */
+  if (    kernel1->width != kernel2->width
+       || kernel1->height != kernel2->height
+       || kernel1->x != kernel2->x
+       || kernel1->y != kernel2->y )
+    return MagickFalse;
+
+  /* check actual kernel values */
+  for (i=0; i < (kernel1->width*kernel1->height); i++) {
+    /* Test for Nan equivelence */
+    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 */
+    if ( fabs(kernel1->values[i] - kernel2->values[i]) > MagickEpsilon )
+      return MagickFalse;
+  }
+
+  return MagickTrue;
+}
+
+static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
 {
   KernelInfo
     *clone,
     *last;
 
-  double
-    a;
-
   last = kernel;
-  for (a=angle; a<355.0; a+=angle) {
+  while(1) {
     clone = CloneKernelInfo(last);
     RotateKernelInfo(clone, angle);
-    last->next = clone;
+    if ( SameKernelInfo(kernel, clone) == MagickTrue )
+      break;
+    LastKernelInfo(last)->next = clone;
     last = clone;
   }
+  clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
+  return;
 }
 \f
 /*
@@ -1758,7 +2316,7 @@ static void ExpandKernelInfo(KernelInfo *kernel, double angle)
 */
 static void CalcKernelMetaData(KernelInfo *kernel)
 {
-  register unsigned long
+  register size_t
     i;
 
   kernel->minimum = kernel->maximum = 0.0;
@@ -1792,18 +2350,23 @@ static void CalcKernelMetaData(KernelInfo *kernel)
 %  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.
+%  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 MorphologyImage method is:
+%  The format of the MorphologyApply method is:
 %
 %      Image *MorphologyApply(const Image *image,MorphologyMethod method,
-%        const long iterations,const KernelInfo *kernel,const double bias,
+%        const ssize_t iterations,const KernelInfo *kernel,
+%        const CompositeMethod compose, const double bias,
 %        ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
@@ -1822,7 +2385,13 @@ static void CalcKernelMetaData(KernelInfo *kernel)
 %    o kernel: An array of double representing the morphology kernel.
 %              Warning: kernel may be normalized for the Convolve method.
 %
-%    o bias: Convolution Bias to use.
+%    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.
+%
+%    o bias: Convolution Output Bias.
 %
 %    o exception: return any errors or warnings in this structure.
 %
@@ -1833,23 +2402,34 @@ static void CalcKernelMetaData(KernelInfo *kernel)
 ** Two pre-created images must be provided, no image is created.
 ** Returning the number of pixels that changed.
 */
-static unsigned long MorphologyPrimitive(const Image *image, Image
+static size_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,
+  CacheView
+    *p_view,
+    *q_view;
+
+  ssize_t
     y, offx, offy,
     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;
@@ -1869,15 +2449,15 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
     case DilateIntensityMorphology:
     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);
@@ -1887,7 +2467,7 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
 #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;
@@ -1904,10 +2484,10 @@ 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)
@@ -1925,12 +2505,12 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
     q_indexes=GetCacheViewAuthenticIndexQueue(q_view);
     r = (image->columns+kernel->width)*offy+offx; /* constant */
 
-    for (x=0; x < (long) image->columns; x++)
+    for (x=0; x < (ssize_t) image->columns; x++)
     {
-       long
+       ssize_t
         v;
 
-      register long
+      register ssize_t
         u;
 
       register const double
@@ -2026,8 +2606,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) ) continue;
                     alpha=(*k)*(QuantumScale*(QuantumRange-
                                           k_pixels[u].opacity));
@@ -2057,8 +2637,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) ) continue;
                     result.red     += (*k)*k_pixels[u].red;
                     result.green   += (*k)*k_pixels[u].green;
@@ -2085,8 +2665,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);
@@ -2117,8 +2697,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);
@@ -2143,14 +2723,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 */
@@ -2176,7 +2757,7 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
               k_pixels += image->columns+kernel->width;
               k_indexes += image->columns+kernel->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 );
@@ -2196,8 +2777,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;
                 if ( result.red == 0.0 ||
                      PixelIntensity(&(k_pixels[u])) < PixelIntensity(q) ) {
@@ -2226,8 +2807,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; /* boolean kernel */
                 if ( result.red == 0.0 ||
                      PixelIntensity(&(k_pixels[u])) > PixelIntensity(q) ) {
@@ -2258,8 +2839,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) ) continue;
                 Minimize(result.red,     (*k)+k_pixels[u].red);
                 Minimize(result.green,   (*k)+k_pixels[u].green);
@@ -2301,12 +2882,12 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
           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 );
+          /* 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 */
@@ -2363,41 +2944,51 @@ static unsigned long MorphologyPrimitive(const Image *image, Image
   result_image->type=image->type;
   q_view=DestroyCacheView(q_view);
   p_view=DestroyCacheView(p_view);
-  return(status ? (unsigned long) changed : 0);
+  return(status ? (size_t) changed : 0);
 }
 
 
 MagickExport Image *MorphologyApply(const Image *image, const ChannelType
-     channel,const MorphologyMethod method, const long iterations,
-     const KernelInfo *kernel,const double bias,ExceptionInfo *exception)
+     channel,const MorphologyMethod method, const ssize_t iterations,
+     const KernelInfo *kernel, const CompositeOperator compose,
+     const double bias, ExceptionInfo *exception)
 {
   Image
-    *curr_image,   /* Image we are working with */
-    *work_image,   /* secondary working image */
-    *save_image;   /* save image for later use */
+    *curr_image,    /* Image we are working with or iterating */
+    *work_image,    /* secondary image for primative iteration */
+    *save_image,    /* saved image - for 'edge' method only */
+    *rslt_image;    /* resultant image - after multi-kernel handling */
 
   KernelInfo
-    *curr_kernel,  /* current kernel list to apply */
-    *this_kernel;  /* current individual kernel to apply */
+    *reflected_kernel, /* A reflected copy of the kernel (if needed) */
+    *norm_kernel,      /* the current normal un-reflected kernel */
+    *rflt_kernel,      /* the current reflected kernel (if needed) */
+    *this_kernel;      /* the kernel being applied */
 
   MorphologyMethod
-    primitive;     /* the current morphology primitive being applied */
+    primative;      /* the current morphology primative being applied */
+
+  CompositeOperator
+    rslt_compose;   /* multi-kernel compose method for results to use */
 
   MagickBooleanType
-    verbose;           /* verbose output of results */
+    verbose;        /* verbose output of results */
+
+  size_t
+    method_loop,    /* Loop 1: number of compound method iterations */
+    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) */
+    kernel_limit,   /*         number of times to iterate kernel */
+    count,          /* total count of primative steps applied */
+    changed,        /* number pixels changed by last primative operation */
+    kernel_changed, /* total count of changed using iterated kernel */
+    method_changed; /* total count of changed over method iteration */
 
-  CompositeOperator
-    kernel_compose;  /* Handling the result of multiple kernels*/
-
-  unsigned long
-    count,         /* count of primitive steps applied */
-    loop,          /* number of times though kernel list (iterations) */
-    loop_limit,    /* finish looping after this many times */
-    stage,         /* stage number for compound morphology */
-    changed,       /* number pixels changed by one primitive operation */
-    loop_changed,  /* changes made over loop though of kernels */
-    total_changed, /* total count of all changes to image */
-    kernel_number; /* kernel number being applied */
+  char
+    v_info[80];
 
   assert(image != (Image *) NULL);
   assert(image->signature == MagickSignature);
@@ -2406,144 +2997,198 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
   assert(exception != (ExceptionInfo *) NULL);
   assert(exception->signature == MagickSignature);
 
-  loop_limit = (unsigned long) iterations;
-  if ( iterations < 0 )
-     loop_limit = image->columns > image->rows ? image->columns : image->rows;
+  count = 0;      /* number of low-level morphology primatives performed */
   if ( iterations == 0 )
-    return((Image *)NULL); /* null operation - nothing to do! */
+    return((Image *)NULL);   /* null operation - nothing to do! */
 
-  /* kernel must be valid at this point
-   * (except maybe for posible future morphology methods like "Prune"
-   */
-  assert(kernel != (KernelInfo *)NULL);
+  kernel_limit = (size_t) iterations;
+  if ( iterations < 0 )  /* negative interations = infinite (well alomst) */
+     kernel_limit = image->columns > image->rows ? image->columns : image->rows;
 
   verbose = ( GetImageArtifact(image,"verbose") != (const char *) NULL ) ?
     MagickTrue : MagickFalse;
 
   /* initialise for cleanup */
-  curr_image = (Image *) image;    /* result of morpholgy primitive */
-  work_image = (Image *) NULL;     /* secondary working image */
-  save_image = (Image *) NULL;     /* save image for some compound methods */
-  curr_kernel = (KernelInfo *)kernel; /* allow kernel list to be modified */
-
-  kernel_compose = NoCompositeOp;  /* iterated over all kernels */
-
-  /* Select initial primitive morphology to apply */
-  primitive = UndefinedMorphology;
+  curr_image = (Image *) image;
+  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
+   * + 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 */
+  rslt_compose = compose; /* and we are composing multi-kernels as given */
   switch( method ) {
-    case CorrelateMorphology:
-      /* A Correlation is actually a Convolution with a reflected kernel.
-      ** However a Convolution is a weighted sum using a reflected kernel.
-      ** It may seem stange to convert a Correlation into a Convolution
-      ** as the Correleation is the simplier method, but Convolution is
-      ** much more commonly used, and it makes sense to implement it directly
-      ** so as to avoid the need to duplicate the kernel when it is not
-      ** required (which is typically the default).
-      */
-      curr_kernel = CloneKernelInfo(kernel);
-      if (curr_kernel == (KernelInfo *) NULL)
-        goto error_cleanup;
-      RotateKernelInfo(curr_kernel,180);
-      /* FALL THRU to Convolve */
-    case ConvolveMorphology:
-      primitive = ConvolveMorphology;
-      kernel_compose = NoCompositeOp;
+    case SmoothMorphology:  /* 4 primative compound morphology */
+      stage_limit = 4;
       break;
-    case ErodeMorphology:      /* just erode */
-    case OpenMorphology:       /* erode then dialate */
-    case EdgeInMorphology:     /* erode and image difference */
-    case TopHatMorphology:     /* erode, dilate and image difference */
-    case SmoothMorphology:     /* erode, dilate, dilate, erode */
-      primitive = ErodeMorphology;
-      break;
-    case ErodeIntensityMorphology:
+    case OpenMorphology:    /* 2 primative compound morphology */
     case OpenIntensityMorphology:
-      primitive = ErodeIntensityMorphology;
-      break;
-    case DilateMorphology:     /* just dilate */
-    case EdgeOutMorphology:    /* dilate and image difference */
-    case EdgeMorphology:       /* dilate and erode difference */
-      primitive = DilateMorphology;
-      break;
-    case CloseMorphology:      /* dilate, then erode */
-    case BottomHatMorphology:  /* dilate and image difference */
-      curr_kernel = CloneKernelInfo(kernel);
-      if (curr_kernel == (KernelInfo *) NULL)
-        goto error_cleanup;
-      RotateKernelInfo(curr_kernel,180);
-      primitive = DilateMorphology;
-      break;
-    case DilateIntensityMorphology:
+    case TopHatMorphology:
+    case CloseMorphology:
     case CloseIntensityMorphology:
-      curr_kernel = CloneKernelInfo(kernel);
-      if (curr_kernel == (KernelInfo *) NULL)
-        goto error_cleanup;
-      RotateKernelInfo(curr_kernel,180);
-      primitive = DilateIntensityMorphology;
+    case BottomHatMorphology:
+    case EdgeMorphology:
+      stage_limit = 2;
       break;
     case HitAndMissMorphology:
-      primitive = HitAndMissMorphology;
-      loop_limit = 1;                       /* iterate only once */
-      kernel_compose = LightenCompositeOp;  /* Union of Hit-And-Miss */
-      break;
-    case ThinningMorphology:     /* iterative morphology */
+      rslt_compose = LightenCompositeOp;  /* Union of multi-kernel results */
+      /* FALL THUR */
+    case ThinningMorphology:
     case ThickenMorphology:
-    case DistanceMorphology:     /* Distance should never use multple kernels */
-    case UndefinedMorphology:
-      primitive = method;
+      method_limit = kernel_limit;  /* iterate the whole method */
+      kernel_limit = 1;             /* do not do kernel iteration  */
+      break;
+    default:
       break;
   }
 
-#if 0
-  { /* User override of results handling  -- Experimental */
-    const char
-       *artifact = GetImageArtifact(image,"morphology:style");
-    if ( artifact != (const char *) NULL ) {
-      if (LocaleCompare("union",artifact) == 0)
-        kernel_compose = LightenCompositeOp;
-      if (LocaleCompare("iterate",artifact) == 0)
-        kernel_compose = NoCompositeOp;
-      else
-        kernel_compose = (CompositeOperator) ParseMagickOption(
-                                 MagickComposeOptions,MagickFalse,artifact);
-      if ( kernel_compose == UndefinedCompositeOp )
-        perror("Invalid \"morphology:compose\" setting\n");
-    }
-  }
-#endif
-
-  /* Initialize compound morphology stages  */
-  count = 0;          /* number of low-level morphology primitives performed */
-  total_changed = 0;  /* total number of pixels changed thoughout */
-  stage = 1;          /* the compound morphology stage number */
+  /* 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 */
 
-  /* compount morphology staging loop */
-  while ( 1 ) {
-
-#if 1
-    /* Extra information for debugging compound operations */
-    if ( verbose == MagickTrue && primitive != method )
-      fprintf(stderr, "Morphology %s: Stage %lu %s%s (%s)\n",
-        MagickOptionToMnemonic(MagickMorphologyOptions, method), stage,
-        MagickOptionToMnemonic(MagickMorphologyOptions, primitive),
-        ( curr_kernel == kernel) ? "" : "*",
-        ( kernel_compose == NoCompositeOp ) ? "iterate"
-          : MagickOptionToMnemonic(MagickComposeOptions, kernel_compose) );
-#endif
+  /* Some methods require a reflected kernel to use with primatives.
+   * Create the reflected kernel for those methods. */
+  switch ( method ) {
+    case CorrelateMorphology:
+    case CloseMorphology:
+    case CloseIntensityMorphology:
+    case BottomHatMorphology:
+    case SmoothMorphology:
+      reflected_kernel = CloneKernelInfo(kernel);
+      if (reflected_kernel == (KernelInfo *) NULL)
+        goto error_cleanup;
+      RotateKernelInfo(reflected_kernel,180);
+      break;
+    default:
+      break;
+  }
 
-    if ( kernel_compose == NoCompositeOp ) {
-      /******************************
-       ** Iterate over all Kernels **
-       ******************************/
-      loop = 0;
-      loop_changed = 1;
-      while ( loop < loop_limit && loop_changed > 0 ) {
-        loop++;       /* the iteration of this kernel */
+  /* Loop 1:  iterate the compound method */
+  method_loop = 0;
+  method_changed = 1;
+  while ( method_loop < method_limit && method_changed > 0 ) {
+    method_loop++;
+    method_changed = 0;
+
+    /* 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 ) {
+
+      /* Loop 3: Compound Morphology Staging - Select Primative to apply */
+      stage_loop = 0;          /* the compound morphology stage number */
+      while ( stage_loop < stage_limit ) {
+        stage_loop++;   /* The stage of the compound morphology */
+
+        /* Select primative morphology for this stage of compound method */
+        this_kernel = norm_kernel; /* default use unreflected kernel */
+        primative = method;        /* Assume method is a primative */
+        switch( method ) {
+          case ErodeMorphology:      /* just erode */
+          case EdgeInMorphology:     /* erode and image difference */
+            primative = ErodeMorphology;
+            break;
+          case DilateMorphology:     /* just dilate */
+          case EdgeOutMorphology:    /* dilate and image difference */
+            primative = DilateMorphology;
+            break;
+          case OpenMorphology:       /* erode then dialate */
+          case TopHatMorphology:     /* open and image difference */
+            primative = ErodeMorphology;
+            if ( stage_loop == 2 )
+              primative = DilateMorphology;
+            break;
+          case OpenIntensityMorphology:
+            primative = ErodeIntensityMorphology;
+            if ( stage_loop == 2 )
+              primative = DilateIntensityMorphology;
+            break;
+          case CloseMorphology:      /* dilate, then erode */
+          case BottomHatMorphology:  /* close and image difference */
+            this_kernel = rflt_kernel; /* use the reflected kernel */
+            primative = DilateMorphology;
+            if ( stage_loop == 2 )
+              primative = ErodeMorphology;
+            break;
+          case CloseIntensityMorphology:
+            this_kernel = rflt_kernel; /* use the reflected kernel */
+            primative = DilateIntensityMorphology;
+            if ( stage_loop == 2 )
+              primative = ErodeIntensityMorphology;
+            break;
+          case SmoothMorphology:         /* open, close */
+            switch ( stage_loop ) {
+              case 1: /* start an open method, which starts with Erode */
+                primative = ErodeMorphology;
+                break;
+              case 2:  /* now Dilate the Erode */
+                primative = DilateMorphology;
+                break;
+              case 3:  /* Reflect kernel a close */
+                this_kernel = rflt_kernel; /* use the reflected kernel */
+                primative = DilateMorphology;
+                break;
+              case 4:  /* Finish the Close */
+                this_kernel = rflt_kernel; /* use the reflected kernel */
+                primative = ErodeMorphology;
+                break;
+            }
+            break;
+          case EdgeMorphology:        /* dilate and erode difference */
+            primative = DilateMorphology;
+            if ( stage_loop == 2 ) {
+              save_image = curr_image;      /* save the image difference */
+              curr_image = (Image *) image;
+              primative = ErodeMorphology;
+            }
+            break;
+          case CorrelateMorphology:
+            /* A Correlation is a Convolution with a reflected kernel.
+            ** However a Convolution is a weighted sum using a reflected
+            ** kernel.  It may seem stange to convert a Correlation into a
+            ** Convolution as the Correlation is the simplier method, but
+            ** Convolution is much more commonly used, and it makes sense to
+            ** implement it directly so as to avoid the need to duplicate the
+            ** kernel when it is not required (which is typically the
+            ** default).
+            */
+            this_kernel = rflt_kernel; /* use the reflected kernel */
+            primative = 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:%.20g.%.20g -> ",
+             MagickOptionToMnemonic(MagickMorphologyOptions,method),(double)
+             method_loop,(double) stage_loop);
+          else if ( primative != method )
+            (void) FormatMagickString(v_info, MaxTextExtent, "%s:%.20g -> ",
+              MagickOptionToMnemonic(MagickMorphologyOptions, method),(double)
+              method_loop);
+          else
+            v_info[0] = '\0';
+        }
 
-        loop_changed = 0;
-        this_kernel = curr_kernel;
-        kernel_number = 0;
-        while ( this_kernel != NULL ) {
+        /* Loop 4: Iterate the kernel with primative */
+        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 */
           if ( work_image == (Image *) NULL )
@@ -2558,237 +3203,153 @@ MagickExport Image *MorphologyApply(const Image *image, const ChannelType
                 }
             }
 
-          /* morphological primitive  curr -> work */
+          /* APPLY THE MORPHOLOGICAL PRIMITIVE (curr -> work) */
           count++;
-          changed = MorphologyPrimitive(curr_image, work_image, primitive,
+          changed = MorphologyPrimitive(curr_image, work_image, primative,
                         channel, this_kernel, bias, exception);
-          loop_changed += changed;
-          total_changed += changed;
-
-          if ( verbose == MagickTrue )
-            fprintf(stderr, "Morphology %s:%lu.%lu #%lu => Changed %lu\n",
-                MagickOptionToMnemonic(MagickMorphologyOptions, primitive),
-                loop, kernel_number, count, changed);
-
+          kernel_changed += changed;
+          method_changed += changed;
+
+          if ( verbose == MagickTrue ) {
+            if ( kernel_loop > 1 )
+              fprintf(stderr, "\n"); /* add end-of-line from previous */
+            (void) fprintf(stderr, "%s%s%s:%.20g.%.20g #%.20g => Changed %.20g",
+              v_info,MagickOptionToMnemonic(MagickMorphologyOptions,
+              primative),(this_kernel == rflt_kernel ) ? "*" : "",
+              (double) (method_loop+kernel_loop-1),(double) kernel_number,
+              (double) count,(double) changed);
+          }
           /* prepare next loop */
           { Image *tmp = work_image;   /* swap images for iteration */
             work_image = curr_image;
             curr_image = tmp;
           }
           if ( work_image == image )
-            work_image = (Image *) NULL;  /* never assign image to 'work' */
-          this_kernel = this_kernel->next;  /* prepare next kernel (if any) */
-          kernel_number++;
-        }
+            work_image = (Image *) NULL; /* replace input 'image' */
 
-        if ( verbose == MagickTrue && kernel->next != NULL )
-          fprintf(stderr, "Morphology %s:%lu #%lu ===> Changed %lu   Total %lu\n",
-                MagickOptionToMnemonic(MagickMorphologyOptions, primitive),
-                loop, count, loop_changed, total_changed );
-      }
-    }
+        } /* End Loop 4: Iterate the kernel with primative */
 
-    else {
-      /*************************************
-       ** Composition of Iterated Kernels **
-       *************************************/
-      Image
-        *input_image,  /* starting point for kernels */
-        *union_image;
-      input_image = curr_image;
-      union_image = (Image *) NULL;
-
-      this_kernel = curr_kernel;
-      kernel_number = 0;
-      while ( this_kernel != NULL ) {
-
-        if( curr_image != (Image *) NULL && curr_image != input_image )
-          curr_image=DestroyImage(curr_image);
-        curr_image = input_image;  /* always start with the original image */
-
-        loop = 0;
-        changed = 1;
-        loop_changed = 0;
-        while ( loop < loop_limit && changed > 0 ) {
-          loop++;       /* the iteration of this kernel */
+        if ( verbose == MagickTrue && kernel_changed != changed )
+          fprintf(stderr, "   Total %.20g",(double) kernel_changed);
+        if ( verbose == MagickTrue && stage_loop < stage_limit )
+          fprintf(stderr, "\n"); /* add end-of-line before looping */
 
-          /* Create a destination image, if not defined */
-          if ( work_image == (Image *) NULL )
-            {
-              work_image=CloneImage(image,0,0,MagickTrue,exception);
-              if (work_image == (Image *) NULL)
-                goto error_cleanup;
-              if (SetImageStorageClass(work_image,DirectClass) == MagickFalse)
-                {
-                  InheritException(exception,&curr_image->exception);
-                  if( union_image != (Image *) NULL )
-                    union_image=DestroyImage(union_image);
-                  if( curr_image != input_image )
-                    curr_image = DestroyImage(curr_image);
-                  curr_image = (Image *) input_image;
-                  goto error_cleanup;
-                }
-            }
+#if 0
+    fprintf(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
+    fprintf(stderr, "      curr =0x%lx\n", (unsigned long)curr_image);
+    fprintf(stderr, "      work =0x%lx\n", (unsigned long)work_image);
+    fprintf(stderr, "      save =0x%lx\n", (unsigned long)save_image);
+    fprintf(stderr, "      union=0x%lx\n", (unsigned long)rslt_image);
+#endif
 
-          /* morphological primitive  curr -> work */
-          count++;
-          changed = MorphologyPrimitive(curr_image,work_image,primitive,
-                        channel, this_kernel, bias, exception);
-          loop_changed += changed;
-          total_changed += changed;
+      } /* End Loop 3: Primative (staging) Loop for Coumpound Methods */
 
+      /*  Final Post-processing for some Compound Methods
+      **
+      ** 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'.
+      */
+      switch( method ) {
+        case EdgeOutMorphology:
+        case EdgeInMorphology:
+        case TopHatMorphology:
+        case BottomHatMorphology:
           if ( verbose == MagickTrue )
-            fprintf(stderr, "Morphology %s:%lu.%lu #%lu => Changed %lu\n",
-                MagickOptionToMnemonic(MagickMorphologyOptions, primitive),
-                loop, kernel_number, count, changed);
-
-          /* prepare next loop */
-          { Image *tmp = work_image;   /* swap images for iteration */
-            work_image = curr_image;   /* curr_image is now the results */
-            curr_image = tmp;
-          }
-          if ( work_image == input_image )
-            work_image = (Image *) NULL;  /* clear work of the input_image */
-
-        } /* end kernel iteration */
-
-        /* make a union of the iterated kernel */
-        if ( union_image == (Image *) NULL)   /* start the union? */
-          union_image = curr_image, curr_image = (Image *)NULL;
-        else
-          (void) CompositeImageChannel(union_image,
-            (ChannelType) (channel & ~SyncChannels), kernel_compose,
-            curr_image, 0, 0);
-
-        this_kernel = this_kernel->next;  /* next kernel (if any) */
-        kernel_number++;
+            fprintf(stderr, "\n%s: Difference with original image",
+                 MagickOptionToMnemonic(MagickMorphologyOptions, method) );
+          (void) CompositeImageChannel(curr_image,
+                  (ChannelType) (channel & ~SyncChannels),
+                  DifferenceCompositeOp, image, 0, 0);
+          break;
+        case EdgeMorphology:
+          if ( verbose == MagickTrue )
+            fprintf(stderr, "\n%s: Difference of Dilate and Erode",
+                 MagickOptionToMnemonic(MagickMorphologyOptions, method) );
+          (void) CompositeImageChannel(curr_image,
+                  (ChannelType) (channel & ~SyncChannels),
+                  DifferenceCompositeOp, save_image, 0, 0);
+          save_image = DestroyImage(save_image); /* finished with save image */
+          break;
+        default:
+          break;
       }
 
-      if ( verbose == MagickTrue && kernel->next != NULL && loop_limit > 1 )
-        fprintf(stderr, "Morphology %s:%lu #%lu ===> Changed %lu   Total %lu\n",
-              MagickOptionToMnemonic(MagickMorphologyOptions, primitive),
-              loop, count, loop_changed, total_changed );
-
-#if 0
-fprintf(stderr, "--E-- image=0x%lx\n", (unsigned long)image);
-fprintf(stderr, "      input=0x%lx\n", (unsigned long)input_image);
-fprintf(stderr, "      union=0x%lx\n", (unsigned long)union_image);
-fprintf(stderr, "      curr =0x%lx\n", (unsigned long)curr_image);
-fprintf(stderr, "      work =0x%lx\n", (unsigned long)work_image);
-fprintf(stderr, "      save =0x%lx\n", (unsigned long)save_image);
-#endif
-
-      /* Finish up - return the union of results */
-      if( curr_image != (Image *) NULL && curr_image != input_image )
-          curr_image=DestroyImage(curr_image);
-      if( input_image != input_image )
-        input_image = DestroyImage(input_image);
-      curr_image = union_image;
-    }
-
-    /* Compound Morphology Operations
-     *   set next 'primitive' iteration, and continue
-     *   or break when all operations are complete.
-     */
-    stage++;   /* what is the next stage number to do */
-    switch( method ) {
-      case SmoothMorphology:           /* open, close */
-        switch ( stage ) {
-        /* case 1:  initialized above */
-        case 2:  /* open part 2 */
-          primitive = DilateMorphology;
-          continue;
-        case 3:  /* close part 1 */
-          curr_kernel = CloneKernelInfo(kernel);
-          if (curr_kernel == (KernelInfo *) NULL)
-            goto error_cleanup;
-          RotateKernelInfo(curr_kernel,180);
-          continue;
-        case 4:  /* close part 2 */
-          primitive = ErodeMorphology;
-          continue;
+      /* multi-kernel handling:  re-iterate, or compose results */
+      if ( kernel->next == (KernelInfo *) NULL )
+        rslt_image = curr_image;   /* just return the resulting image */
+      else if ( rslt_compose == NoCompositeOp )
+        { if ( verbose == MagickTrue ) {
+            if ( this_kernel->next != (KernelInfo *) NULL )
+              fprintf(stderr, " (re-iterate)");
+            else
+              fprintf(stderr, " (done)");
+          }
+          rslt_image = curr_image; /* return result, and re-iterate */
         }
-        break;
-      case OpenMorphology:      /* erode, dilate */
-      case TopHatMorphology:
-        primitive = DilateMorphology;
-        if ( stage <= 2 ) continue;
-        break;
-      case OpenIntensityMorphology:
-        primitive = DilateIntensityMorphology;
-        if ( stage <= 2 ) continue;
-        break;
-      case CloseMorphology:       /* dilate, erode */
-      case BottomHatMorphology:
-        primitive = ErodeMorphology;
-        if ( stage <= 2 ) continue;
-        break;
-      case CloseIntensityMorphology:
-        primitive = ErodeIntensityMorphology;
-        if ( stage <= 2 ) continue;
-        break;
-      case EdgeMorphology:        /* dilate and erode difference */
-        if (stage <= 2) {
-          save_image = curr_image;
-          curr_image = (Image *) image;
-          primitive = ErodeMorphology;
-          continue;
+      else if ( rslt_image == (Image *) NULL)
+        { if ( verbose == MagickTrue )
+            fprintf(stderr, " (save for compose)");
+          rslt_image = curr_image;
+          curr_image = (Image *) image;  /* continue with original image */
         }
-        break;
-      default:  /* Primitive Morphology is just finished! */
-        break;
-    }
+      else
+        { /* 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'.
+          **
+          ** The compose image order is specifically so that the new image can
+          ** be subtarcted 'Minus' from the collected result, to allow you to
+          ** convert a HitAndMiss methd into a Thinning method.
+          */
+          if ( verbose == MagickTrue )
+            fprintf(stderr, " (compose \"%s\")",
+                 MagickOptionToMnemonic(MagickComposeOptions, rslt_compose) );
+          (void) CompositeImageChannel(curr_image,
+               (ChannelType) (channel & ~SyncChannels), rslt_compose,
+               rslt_image, 0, 0);
+          rslt_image = DestroyImage(rslt_image);
+          rslt_image = curr_image;
+          curr_image = (Image *) image;  /* continue with original image */
+        }
+      if ( verbose == MagickTrue )
+        fprintf(stderr, "\n");
 
-    if ( verbose == MagickTrue && count > 1 )
-      fprintf(stderr, "Morphology %s: ======> Total %lu\n",
-           MagickOptionToMnemonic(MagickMorphologyOptions, method),
-           total_changed );
+      /* loop to the next kernel in a multi-kernel list */
+      norm_kernel = norm_kernel->next;
+      if ( rflt_kernel != (KernelInfo *) NULL )
+        rflt_kernel = rflt_kernel->next;
+      kernel_number++;
+    } /* End Loop 2: Loop over each kernel */
 
-    /* If we reach this point we are finished! - Break the Loop */
-    break;
-  }
-
-  /*  Final Post-processing for some Compound Methods
-  **
-  ** The removal of any 'Sync' channel flag in the Image Compositon below
-  ** ensures the compose method is applied in a purely mathematical way, only
-  ** the selected channels, without any normal 'alpha blending' normally
-  ** associated with the compose method.
-  **
-  ** Note "method" here is the 'original' morphological method, and not the
-  ** 'current' morphological method used above to generate "new_image".
-  */
-  switch( method ) {
-    case EdgeOutMorphology:
-    case EdgeInMorphology:
-    case TopHatMorphology:
-    case BottomHatMorphology:
-      (void) CompositeImageChannel(curr_image,
-               (ChannelType) (channel & ~SyncChannels), DifferenceCompositeOp,
-               image, 0, 0);
-      break;
-    case EdgeMorphology:
-      /* Difference the Eroded Image with a Dilate image */
-      (void) CompositeImageChannel(curr_image,
-               (ChannelType) (channel & ~SyncChannels), DifferenceCompositeOp,
-               save_image, 0, 0);
-      break;
-    default:
-      break;
-  }
+  } /* End Loop 1: compound method interation */
 
   goto exit_cleanup;
 
-  /* Yes goto's are bad, but in this case it makes cleanup lot more efficient */
+  /* Yes goto's are bad, but it makes cleanup lot more efficient */
 error_cleanup:
-  if ( curr_image != (Image *) NULL && curr_image != image )
-    (void) DestroyImage(curr_image);
+  if ( curr_image != (Image *) NULL &&
+       curr_image != rslt_image &&
+       curr_image != image )
+    curr_image = DestroyImage(curr_image);
+  if ( rslt_image != (Image *) NULL )
+    rslt_image = DestroyImage(rslt_image);
 exit_cleanup:
+  if ( curr_image != (Image *) NULL &&
+       curr_image != rslt_image &&
+       curr_image != image )
+    curr_image = DestroyImage(curr_image);
   if ( work_image != (Image *) NULL )
-    (void) DestroyImage(work_image);
+    work_image = DestroyImage(work_image);
   if ( save_image != (Image *) NULL )
-    (void) DestroyImage(save_image);
-  return(curr_image);
+    save_image = DestroyImage(save_image);
+  if ( reflected_kernel != (KernelInfo *) NULL )
+    reflected_kernel = DestroyKernelInfo(reflected_kernel);
+  return(rslt_image);
 }
 \f
 /*
@@ -2817,10 +3378,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:
@@ -2845,7 +3406,7 @@ 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;
@@ -2853,6 +3414,9 @@ MagickExport Image *MorphologyImageChannel(const Image *image,
   KernelInfo
     *curr_kernel;
 
+  CompositeOperator
+    compose;
+
   Image
     *morphology_image;
 
@@ -2865,7 +3429,7 @@ MagickExport Image *MorphologyImageChannel(const Image *image,
   if ( method == ConvolveMorphology ||  method == CorrelateMorphology )
     {
       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) {
@@ -2878,12 +3442,28 @@ 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)
     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.
+   * Default for 'HitAndMiss' is 'Lighten'.
+   */
+  compose = UndefinedCompositeOp;  /* use default for method */
+  artifact = GetImageArtifact(image,"morphology:compose");
+  if ( artifact != (const char *) NULL)
+    compose = (CompositeOperator) ParseMagickOption(
+                             MagickComposeOptions,MagickFalse,artifact);
+
   /* Apply the Morphology */
   morphology_image = MorphologyApply(image, channel, method, iterations,
-                         curr_kernel, image->bias, exception);
+                         curr_kernel, compose, image->bias, exception);
 
   /* Cleanup and Exit */
   if ( curr_kernel != kernel )
@@ -2891,9 +3471,8 @@ MagickExport Image *MorphologyImageChannel(const Image *image,
   return(morphology_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
@@ -2957,12 +3536,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;
 
@@ -3001,7 +3581,18 @@ static void RotateKernelInfo(KernelInfo *kernel, double angle)
           kernel->values[5] = kernel->values[2];
           kernel->values[2] = kernel->values[1];
           kernel->values[1] = t;
-          /* NOT DONE - rotate a off-centered origin as well! */
+          /* rotate non-centered origin */
+          if ( kernel->x != 1 || 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 = (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);
         }
@@ -3011,14 +3602,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;
@@ -3032,20 +3623,27 @@ 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
-            i,j,x,y;
-          register MagickRealType
-            *k,t;
-          k=kernel->values;
-          for( i=0, x=kernel->width-1;  i<=x;   i++, x--)
-            for( j=0, y=kernel->height-1;  j<y;   j++, y--)
-              { t                    = k[i+j*kernel->width];
-                k[i+j*kernel->width] = k[j+x*kernel->width];
-                k[j+x*kernel->width] = k[x+y*kernel->width];
-                k[x+y*kernel->width] = k[y+i*kernel->width];
-                k[y+i*kernel->width] = t;
-              }
-          /* NOT DONE - rotate a off-centered origin as well! */
+          { register size_t
+              i,j,x,y;
+            register MagickRealType
+              *k,t;
+            k=kernel->values;
+            for( i=0, x=kernel->width-1;  i<=x;   i++, x--)
+              for( j=0, y=kernel->height-1;  j<y;   j++, y--)
+                { t                    = k[i+j*kernel->width];
+                  k[i+j*kernel->width] = k[j+x*kernel->width];
+                  k[j+x*kernel->width] = k[x+y*kernel->width];
+                  k[x+y*kernel->width] = k[y+i*kernel->width];
+                  k[y+i*kernel->width] = t;
+                }
+          }
+          /* rotate the origin - relative to center of array */
+          { 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);
         }
@@ -3059,7 +3657,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;
@@ -3068,8 +3666,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);
     }
@@ -3078,24 +3676,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
@@ -3241,7 +3821,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
@@ -3255,11 +3835,12 @@ MagickExport void ScaleKernelInfo(KernelInfo *kernel,
   /* Normalization of Kernel */
   pos_scale = 1.0;
   if ( (normalize_flags&NormalizeValue) != 0 ) {
-    /* non-zero and zero-summing kernels */
     if ( fabs(kernel->positive_range + kernel->negative_range) > MagickEpsilon )
+      /* non-zero-summing kernel (generally positive) */
       pos_scale = fabs(kernel->positive_range + kernel->negative_range);
     else
-      pos_scale = kernel->positive_range; /* special zero-summing kernel */
+      /* zero-summing kernel */
+      pos_scale = kernel->positive_range;
   }
   /* Force kernel into a normalized zero-summing kernel */
   if ( (normalize_flags&CorrelateNormalizeValue) != 0 ) {
@@ -3275,7 +3856,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;
 
@@ -3328,21 +3909,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) );
     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,
@@ -3358,12 +3938,12 @@ 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()+2, "nan");
+          fprintf(stderr," %*s", GetMagickPrecision()+3, "nan");
         else
-          fprintf(stderr," %*.*lg", GetMagickPrecision()+2,
+          fprintf(stderr," %*.*lg", GetMagickPrecision()+3,
               GetMagickPrecision(), k->values[i]);
       fprintf(stderr,"\n");
     }
@@ -3387,13 +3967,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:
 %
@@ -3449,7 +4024,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 */