morphology.c

Go to the documentation of this file.

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%    M   M    OOO    RRRR   PPPP   H   H   OOO   L       OOO    GGGG  Y   Y   %
%    MM MM   O   O   R   R  P   P  H   H  O   O  L      O   O  G       Y Y    %
%    M M M   O   O   RRRR   PPPP   HHHHH  O   O  L      O   O  G GGG    Y     %
%    M   M   O   O   R R    P      H   H  O   O  L      O   O  G   G    Y     %
%    M   M    OOO    R  R   P      H   H   OOO   LLLLL   OOO    GGG     Y     %
%                                                                             %
%                                                                             %
%                        MagickCore Morphology Methods                        %
%                                                                             %
%                              Software Design                                %
%                              Anthony Thyssen                                %
%                               January 2010                                  %
%                                                                             %
%                                                                             %
%  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
%  dedicated to making software imaging solutions freely available.           %
%                                                                             %
%  You may not use this file except in compliance with the License.  You may  %
%  obtain a copy of the License at                                            %
%                                                                             %
%    http://www.imagemagick.org/script/license.php                            %
%                                                                             %
%  Unless required by applicable law or agreed to in writing, software        %
%  distributed under the License is distributed on an "AS IS" BASIS,          %
%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
%  See the License for the specific language governing permissions and        %
%  limitations under the License.                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Morpology is the the application of various kernels, of any size and even
% shape, to a image in various ways (typically binary, but not always).
%
% Convolution (weighted sum or average) is just one specific type of
% morphology. Just one that is very common for image bluring and sharpening
% effects.  Not only 2D Gaussian blurring, but also 2-pass 1D Blurring.
%
% This module provides not only a general morphology function, and the ability
% to apply more advanced or iterative morphologies, but also functions for the
% generation of many different types of kernel arrays from user supplied
% arguments. Prehaps even the generation of a kernel from a small image.
*/

/*
  Include declarations.
*/
#include "MagickCore/studio.h"
#include "MagickCore/artifact.h"
#include "MagickCore/cache-view.h"
#include "MagickCore/color-private.h"
#include "MagickCore/enhance.h"
#include "MagickCore/exception.h"
#include "MagickCore/exception-private.h"
#include "MagickCore/gem.h"
#include "MagickCore/gem-private.h"
#include "MagickCore/hashmap.h"
#include "MagickCore/image.h"
#include "MagickCore/image-private.h"
#include "MagickCore/list.h"
#include "MagickCore/magick.h"
#include "MagickCore/memory_.h"
#include "MagickCore/monitor-private.h"
#include "MagickCore/morphology.h"
#include "MagickCore/morphology-private.h"
#include "MagickCore/option.h"
#include "MagickCore/pixel-accessor.h"
#include "MagickCore/prepress.h"
#include "MagickCore/quantize.h"
#include "MagickCore/registry.h"
#include "MagickCore/semaphore.h"
#include "MagickCore/splay-tree.h"
#include "MagickCore/statistic.h"
#include "MagickCore/string_.h"
#include "MagickCore/string-private.h"
#include "MagickCore/token.h"
#include "MagickCore/utility.h"
#include "MagickCore/utility-private.h"


/*
** 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))

/*
  Other global definitions used by module.
*/
static inline double MagickMin(const double x,const double y)
{
  return( x < y ? x : y);
}
static inline double MagickMax(const double x,const double y)
{
  return( x > y ? x : y);
}
#define Minimize(assign,value) assign=MagickMin(assign,value)
#define Maximize(assign,value) assign=MagickMax(assign,value)

/* Currently these are only internal to this module */
static void
  CalcKernelMetaData(KernelInfo *),
  ExpandMirrorKernelInfo(KernelInfo *),
  ExpandRotateKernelInfo(KernelInfo *, const double),
  RotateKernelInfo(KernelInfo *, double);


/* Quick function to find last kernel in a kernel list */
static inline KernelInfo *LastKernelInfo(KernelInfo *kernel)
{
  while (kernel->next != (KernelInfo *) NULL)
    kernel = kernel->next;
  return(kernel);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     A c q u i r e K e r n e l I n f o                                       %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AcquireKernelInfo() takes the given string (generally supplied by the
%  user) and converts it into a Morphology/Convolution Kernel.  This allows
%  users to specify a kernel from a number of pre-defined kernels, or to fully
%  specify their own kernel for a specific Convolution or Morphology
%  Operation.
%
%  The kernel so generated can be any rectangular array of floating point
%  values (doubles) with the 'control point' or 'pixel being affected'
%  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
%  with any value being declared the origin. This in turn allows the use of
%  highly asymmetrical kernels.
%
%  The floating point values in the kernel can also include a special value
%  known as 'nan' or 'not a number' to indicate that this value is not part
%  of the kernel array. This allows you to shaped the kernel within its
%  rectangular area. That is 'nan' values provide a 'mask' for the kernel
%  shape.  However at least one non-nan value must be provided for correct
%  working of a kernel.
%
%  The returned kernel should be freed using the DestroyKernelInfo() when you
%  are finished with it.  Do not free this memory yourself.
%
%  Input kernel defintion strings can consist of any of three types.
%
%    "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 ..."
%         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.
%
%    "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
%         square kernel, 25 for a 5x5 square kernel, 49 for 7x7, etc.
%         Values can be space or comma separated.  This is not recommended.
%
%  You can define a 'list of kernels' which can be used by some morphology
%  operators A list is defined as a semi-colon separated list kernels.
%
%     " kernel ; kernel ; kernel ; "
%
%  Any extra ';' characters, at 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
%  numbers generating a odd-sized square kernel has been given.
%
%  The format of the AcquireKernal method is:
%
%      KernelInfo *AcquireKernelInfo(const char *kernel_string)
%
%  A description of each parameter follows:
%
%    o kernel_string: the Morphology/Convolution kernel wanted.
%
*/

/* This was separated so that it could be used as a separate
** array input handling function, such as for -color-matrix
*/
static KernelInfo *ParseKernelArray(const char *kernel_string)
{
  KernelInfo
    *kernel;

  char
    token[MaxTextExtent];

  const char
    *p,
    *end;

  register ssize_t
    i;

  double
    nan = sqrt((double)-1.0);  /* Special Value : Not A Number */

  MagickStatusType
    flags;

  GeometryInfo
    args;

  kernel=(KernelInfo *) AcquireQuantumMemory(1,sizeof(*kernel));
  if (kernel == (KernelInfo *)NULL)
    return(kernel);
  (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
  kernel->minimum = kernel->maximum = kernel->angle = 0.0;
  kernel->negative_range = kernel->positive_range = 0.0;
  kernel->type = UserDefinedKernel;
  kernel->next = (KernelInfo *) NULL;
  kernel->signature = MagickSignature;
  if (kernel_string == (const char *) NULL)
    return(kernel);

  /* find end of this specific kernel definition string */
  end = strchr(kernel_string, ';');
  if ( end == (char *) NULL )
    end = strchr(kernel_string, '\0');

  /* clear flags - for Expanding kernel lists thorugh rotations */
   flags = NoValue;

  /* Has a ':' in argument - New user kernel specification */
  p = strchr(kernel_string, ':');
  if ( p != (char *) NULL && p < end)
    {
      /* ParseGeometry() needs the geometry separated! -- Arrgghh */
      memcpy(token, kernel_string, (size_t) (p-kernel_string));
      token[p-kernel_string] = '\0';
      SetGeometryInfo(&args);
      flags = ParseGeometry(token, &args);

      /* Size handling and checks of geometry settings */
      if ( (flags & WidthValue) == 0 ) /* if no width then */
        args.rho = args.sigma;         /* then  width = height */
      if ( args.rho < 1.0 )            /* if width too small */
         args.rho = 1.0;               /* then  width = 1 */
      if ( args.sigma < 1.0 )          /* if height too small */
        args.sigma = args.rho;         /* then  height = width */
      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) ? (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 ':' */
    }
  else
    { /* ELSE - Old old specification, forming odd-square kernel */
      /* count up number of values given */
      p=(const char *) kernel_string;
      while ((isspace((int) ((unsigned char) *p)) != 0) || (*p == '\''))
        p++;  /* ignore "'" chars for convolve filter usage - Cristy */
      for (i=0; p < end; i++)
      {
        GetMagickToken(p,&p,token);
        if (*token == ',')
          GetMagickToken(p,&p,token);
      }
      /* set the size of the kernel - old sized square */
      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 */
    }

  /* Read in the kernel values from rest of input string argument */
  kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
    kernel->height*sizeof(*kernel->values));
  if (kernel->values == (MagickRealType *) NULL)
    return(DestroyKernelInfo(kernel));
  kernel->minimum = +MagickHuge;
  kernel->maximum = -MagickHuge;
  kernel->negative_range = kernel->positive_range = 0.0;
  for (i=0; (i < (ssize_t) (kernel->width*kernel->height)) && (p < end); i++)
  {
    GetMagickToken(p,&p,token);
    if (*token == ',')
      GetMagickToken(p,&p,token);
    if (    LocaleCompare("nan",token) == 0
        || LocaleCompare("-",token) == 0 ) {
      kernel->values[i] = nan; /* do not include this value in kernel */
    }
    else {
      kernel->values[i]=StringToDouble(token,(char **) NULL);
      ( kernel->values[i] < 0)
          ?  ( kernel->negative_range += kernel->values[i] )
          :  ( kernel->positive_range += kernel->values[i] );
      Minimize(kernel->minimum, kernel->values[i]);
      Maximize(kernel->maximum, kernel->values[i]);
    }
  }

  /* sanity check -- no more values in kernel definition */
  GetMagickToken(p,&p,token);
  if ( *token != '\0' && *token != ';' && *token != '\'' )
    return(DestroyKernelInfo(kernel));

#if 0
  /* this was the old method of handling a incomplete kernel */
  if ( i < (ssize_t) (kernel->width*kernel->height) ) {
    Minimize(kernel->minimum, kernel->values[i]);
    Maximize(kernel->maximum, kernel->values[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 < (ssize_t) (kernel->width*kernel->height) )
    return(DestroyKernelInfo(kernel));
#endif

  /* check that we recieved at least one real (non-nan) value! */
  if ( kernel->minimum == MagickHuge )
    return(DestroyKernelInfo(kernel));

  if ( (flags & AreaValue) != 0 )         /* '@' symbol in kernel size */
    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)
{
  char
    token[MaxTextExtent];

  const char
    *p,
    *end;

  GeometryInfo
    args;

  KernelInfo
    *kernel;

  MagickStatusType
    flags;

  ssize_t
    type;

  /* Parse special 'named' kernel */
  GetMagickToken(kernel_string,&p,token);
  type=ParseCommandOption(MagickKernelOptions,MagickFalse,token);
  if ( type < 0 || type == UserDefinedKernel )
    return((KernelInfo *)NULL);  /* not a valid named kernel */

  while (((isspace((int) ((unsigned char) *p)) != 0) ||
          (*p == ',') || (*p == ':' )) && (*p != '\0') && (*p != ';'))
    p++;

  end = strchr(p, ';'); /* end of this kernel defintion */
  if ( end == (char *) NULL )
    end = strchr(p, '\0');

  /* ParseGeometry() needs the geometry separated! -- Arrgghh */
  memcpy(token, p, (size_t) (end-p));
  token[end-p] = '\0';
  SetGeometryInfo(&args);
  flags = ParseGeometry(token, &args);

#if 0
  /* For Debugging Geometry Input */
  (void) FormatLocaleFile(stderr, "Geometry = 0x%04X : %lg x %lg %+lg %+lg\n",
    flags, args.rho, args.sigma, args.xi, args.psi );
#endif

  /* special handling of missing values in input string */
  switch( type ) {
    /* Shape Kernel Defaults */
    case UnityKernel:
      if ( (flags & WidthValue) == 0 )
        args.rho = 1.0;    /* Default scale = 1.0, zero is valid */
      break;
    case SquareKernel:
    case DiamondKernel:
    case OctagonKernel:
    case DiskKernel:
    case PlusKernel:
    case CrossKernel:
      if ( (flags & HeightValue) == 0 )
        args.sigma = 1.0;    /* Default scale = 1.0, zero is valid */
      break;
    case RingKernel:
      if ( (flags & XValue) == 0 )
        args.xi = 1.0;       /* Default scale = 1.0, zero is valid */
      break;
    case RectangleKernel:    /* Rectangle - set size defaults */
      if ( (flags & WidthValue) == 0 ) /* if no width then */
        args.rho = args.sigma;         /* then  width = height */
      if ( args.rho < 1.0 )            /* if width too small */
          args.rho = 3;                /* then  width = 3 */
      if ( args.sigma < 1.0 )          /* if height too small */
        args.sigma = args.rho;         /* then  height = width */
      if ( (flags & XValue) == 0 )     /* center offset if not defined */
        args.xi = (double)(((ssize_t)args.rho-1)/2);
      if ( (flags & YValue) == 0 )
        args.psi = (double)(((ssize_t)args.sigma-1)/2);
      break;
    /* Distance Kernel Defaults */
    case ChebyshevKernel:
    case ManhattanKernel:
    case OctagonalKernel:
    case EuclideanKernel:
      if ( (flags & HeightValue) == 0 )           /* no distance scale */
        args.sigma = 100.0;                       /* default distance scaling */
      else if ( (flags & AspectValue ) != 0 )     /* '!' flag */
        args.sigma = QuantumRange/(args.sigma+1); /* maximum pixel distance */
      else if ( (flags & PercentValue ) != 0 )    /* '%' flag */
        args.sigma *= QuantumRange/100.0;         /* percentage of color range */
      break;
    default:
      break;
  }

  kernel = AcquireKernelBuiltIn((KernelInfoType)type, &args);
  if ( kernel == (KernelInfo *) NULL )
    return(kernel);

  /* global expand to rotated kernel list - only for single kernels */
  if ( kernel->next == (KernelInfo *) NULL ) {
    if ( (flags & AreaValue) != 0 )         /* '@' symbol in kernel args */
      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)
{

  KernelInfo
    *kernel,
    *new_kernel;

  char
    token[MaxTextExtent];

  const char
    *p;

  size_t
    kernel_number;

  if (kernel_string == (const char *) NULL)
    return(ParseKernelArray(kernel_string));
  p = kernel_string;
  kernel = NULL;
  kernel_number = 0;

  while ( GetMagickToken(p,NULL,token),  *token != '\0' ) {

    /* ignore extra or multiple ';' kernel separators */
    if ( *token != ';' ) {

      /* tokens starting with alpha is a Named kernel */
      if (isalpha((int) *token) != 0)
        new_kernel = ParseKernelName(p);
      else /* otherwise a user defined kernel array */
        new_kernel = ParseKernelArray(p);

      /* Error handling -- this is not proper error handling! */
      if ( new_kernel == (KernelInfo *) NULL ) {
        (void) FormatLocaleFile(stderr, "Failed to parse kernel number #%.20g\n",
          (double) kernel_number);
        if ( kernel != (KernelInfo *) NULL )
          kernel=DestroyKernelInfo(kernel);
        return((KernelInfo *) NULL);
      }

      /* initialise or append the kernel list */
      if ( kernel == (KernelInfo *) NULL )
        kernel = new_kernel;
      else
        LastKernelInfo(kernel)->next = new_kernel;
    }

    /* look for the next kernel in list */
    p = strchr(p, ';');
    if ( p == (char *) NULL )
      break;
    p++;

  }
  return(kernel);
}


/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     A c q u i r e K e r n e l B u i l t I n                                 %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AcquireKernelBuiltIn() returned one of the 'named' built-in types of
%  kernels used for special purposes such as gaussian blurring, skeleton
%  pruning, and edge distance determination.
%
%  They take a KernelType, and a set of geometry style arguments, which were
%  typically decoded from a user supplied string, or from a more complex
%  Morphology Method that was requested.
%
%  The format of the AcquireKernalBuiltIn method is:
%
%      KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
%           const GeometryInfo args)
%
%  A description of each parameter follows:
%
%    o type: the pre-defined type of kernel wanted
%
%    o args: arguments defining or modifying the kernel
%
%  Convolution Kernels
%
%    Unity
%       The a No-Op or Scaling single element kernel.
%
%    Gaussian:{radius},{sigma}
%       Generate a two-dimensional gaussian kernel, as used by -gaussian.
%       The sigma for the curve is required.  The resulting kernel is
%       normalized,
%
%       If 'sigma' is zero, you get a single pixel on a field of zeros.
%
%       NOTE: that the 'radius' is optional, but if provided can limit (clip)
%       the final size of the resulting kernel to a square 2*radius+1 in size.
%       The radius should be at least 2 times that of the sigma value, or
%       sever clipping and aliasing may result.  If not given or set to 0 the
%       radius will be determined so as to produce the best minimal error
%       result, which is usally much larger than is normally needed.
%
%    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.
%
%    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
%       kernel is clipped to a width of 2*radius+1.  Kernel can be rotated
%       by a 90 degree angle.
%
%       If 'sigma' is zero, you get a single pixel on a field of zeros.
%
%       Note that two convolutions with two "Blur" kernels perpendicular to
%       each other, is equivalent to a far larger "Gaussian" kernel with the
%       same sigma value, However it is much faster to apply. This is how the
%       "-blur" operator actually works.
%
%    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,
%       Adding two such blurs in opposite directions produces a Blur Kernel.
%       Angle can be rotated in multiples of 90 degrees.
%
%       Note that the first argument is the width of the kernel and not the
%       radius of the kernel.
%
%    # Still to be implemented...
%    #
%    # Filter2D
%    # Filter1D
%    #    Set kernel values using a resize filter, and given scale (sigma)
%    #    Cylindrical or Linear.   Is this possible with an image?
%    #
%
%  Named Constant Convolution Kernels
%
%  All these are unscaled, zero-summing kernels by default. As such for
%  non-HDRI version of ImageMagick some form of normalization, user scaling,
%  and biasing the results is recommended, to prevent the resulting image
%  being 'clipped'.
%
%  The 3x3 kernels (most of these) can be circularly rotated in multiples of
%  45 degrees to generate the 8 angled varients of each of the kernels.
%
%    Laplacian:{type}
%      Discrete Lapacian Kernels, (without normalization)
%        Type 0 :  3x3 with center:8 surounded by -1  (8 neighbourhood)
%        Type 1 :  3x3 with center:4 edge:-1 corner:0 (4 neighbourhood)
%        Type 2 :  3x3 with center:4 edge:1 corner:-2
%        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)
%
%    Sobel:{angle}
%      Sobel 'Edge' convolution kernel (3x3)
%          | -1, 0, 1 |
%          | -2, 0,-2 |
%          | -1, 0, 1 |
%
%    Roberts:{angle}
%      Roberts convolution kernel (3x3)
%          |  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 |
%
%    Compass:{angle}
%      Prewitt's "Compass" convolution kernel (3x3)
%          | -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 |
%
%    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
%
%    Diamond:[{radius}[,{scale}]]
%       Generate a diamond shaped kernel with given radius to the points.
%       Kernel size will again be radius*2+1 square and defaults to radius 1,
%       generating a 3x3 kernel that is slightly larger than a square.
%
%    Square:[{radius}[,{scale}]]
%       Generate a square shaped kernel of size radius*2+1, and defaulting
%       to a 3x3 (radius 1).
%
%    Octagon:[{radius}[,{scale}]]
%       Generate octagonal shaped kernel of given radius and constant scale.
%       Default radius is 3 producing a 7x7 kernel. A radius of 1 will result
%       in "Diamond" kernel.
%
%    Disk:[{radius}[,{scale}]]
%       Generate a binary disk, thresholded at the radius given, the radius
%       may be a float-point value. Final Kernel size is floor(radius)*2+1
%       square. A radius of 5.3 is the default.
%
%       NOTE: That a low radii Disk kernels produce the same results as
%       many of the previously defined kernels, but differ greatly at larger
%       radii.  Here is a table of equivalences...
%          "Disk:1"    => "Diamond", "Octagon:1", or "Cross:1"
%          "Disk:1.5"  => "Square"
%          "Disk:2"    => "Diamond:2"
%          "Disk:2.5"  => "Octagon"
%          "Disk:2.9"  => "Square:2"
%          "Disk:3.5"  => "Octagon:3"
%          "Disk:4.5"  => "Octagon:4"
%          "Disk:5.4"  => "Octagon:5"
%          "Disk:6.4"  => "Octagon:6"
%       All other Disk shapes are unique to this kernel, but because a "Disk"
%       is more circular when using a larger radius, using a larger radius is
%       preferred over iterating the morphological operation.
%
%    Rectangle:{geometry}
%       Simply generate a rectangle of 1's with the size given. You can also
%       specify the location of the 'control point', otherwise the closest
%       pixel to the center of the rectangle is selected.
%
%       Properly centered and odd sized rectangles work the best.
%
%  Symbol Dilation Kernels
%
%    These kernel is not a good general morphological kernel, but is used
%    more for highlighting and marking any single pixels in an image using,
%    a "Dilate" method as appropriate.
%
%    For the same reasons iterating these kernels does not produce the
%    same result as using a larger radius for the symbol.
%
%    Plus:[{radius}[,{scale}]]
%    Cross:[{radius}[,{scale}]]
%       Generate a kernel in the shape of a 'plus' or a 'cross' with
%       a each arm the length of the given radius (default 2).
%
%       NOTE: "plus:1" is equivalent to a "Diamond" kernel.
%
%    Ring:{radius1},{radius2}[,{scale}]
%       A ring of the values given that falls between the two radii.
%       Defaults to a ring of approximataly 3 radius in a 7x7 kernel.
%       This is the 'edge' pixels of the default "Disk" kernel,
%       More specifically, "Ring" -> "Ring:2.5,3.5,1.0"
%
%  Hit and Miss Kernels
%
%    Peak:radius1,radius2
%       Find any peak larger than the pixels the fall between the two radii.
%       The default ring of pixels is as per "Ring".
%    Edges
%       Find flat orthogonal edges of a binary shape
%    Corners
%       Find 90 degree corners of a binary shape
%    Diagonals:type
%       A special kernel to thin the 'outside' of diagonals
%    LineEnds:type
%       Find end points of lines (for pruning a skeletion)
%       Two types of lines ends (default to both) can be searched for
%         Type 0: All line ends
%         Type 1: single kernel for 4-conneected line ends
%         Type 2: single kernel for simple line ends
%    LineJunctions
%       Find three line junctions (within a skeletion)
%         Type 0: all line junctions
%         Type 1: Y Junction kernel
%         Type 2: Diagonal T Junction kernel
%         Type 3: Orthogonal T Junction kernel
%         Type 4: Diagonal X Junction kernel
%         Type 5: Orthogonal + Junction kernel
%    Ridges:type
%       Find single pixel ridges or thin lines
%         Type 1: Fine single pixel thick lines and ridges
%         Type 2: Find two pixel thick lines and ridges
%    ConvexHull
%       Octagonal Thickening Kernel, to generate convex hulls of 45 degrees
%    Skeleton:type
%       Traditional skeleton generating kernels.
%         Type 1: Tradional Skeleton kernel (4 connected skeleton)
%         Type 2: HIPR2 Skeleton kernel (8 connected skeleton)
%         Type 3: Thinning skeleton based on a ressearch paper by
%                 Dan S. Bloomberg (Default Type)
%    ThinSE:type
%       A huge variety of Thinning Kernels designed to preserve conectivity.
%       many other kernel sets use these kernels as source definitions.
%       Type numbers are 41-49, 81-89, 481, and 482 which are based on
%       the super and sub notations used in the source research paper.
%
%  Distance Measuring Kernels
%
%    Different types of distance measuring methods, which are used with the
%    a 'Distance' morphology method for generating a gradient based on
%    distance from an edge of a binary shape, though there is a technique
%    for handling a anti-aliased shape.
%
%    See the 'Distance' Morphological Method, for information of how it is
%    applied.
%
%    Chebyshev:[{radius}][x{scale}[%!]]
%       Chebyshev Distance (also known as Tchebychev or Chessboard distance)
%       is a value of one to any neighbour, orthogonal or diagonal. One why
%       of thinking of it is the number of squares a 'King' or 'Queen' in
%       chess needs to traverse reach any other position on a chess board.
%       It results in a 'square' like distance function, but one where
%       diagonals are given a value that is closer than expected.
%
%    Manhattan:[{radius}][x{scale}[%!]]
%       Manhattan Distance (also known as Rectilinear, City Block, or the Taxi
%       Cab distance metric), it is the distance needed when you can only
%       travel in horizontal or vertical directions only.  It is the
%       distance a 'Rook' in chess would have to travel, and results in a
%       diamond like distances, where diagonals are further than expected.
%
%    Octagonal:[{radius}][x{scale}[%!]]
%       An interleving of Manhatten and Chebyshev metrics producing an
%       increasing octagonally shaped distance.  Distances matches those of
%       the "Octagon" shaped kernel of the same radius.  The minimum radius
%       and default is 2, producing a 5x5 kernel.
%
%    Euclidean:[{radius}][x{scale}[%!]]
%       Euclidean distance is the 'direct' or 'as the crow flys' distance.
%       However by default the kernel size only has a radius of 1, which
%       limits the distance to 'Knight' like moves, with only orthogonal and
%       diagonal measurements being correct.  As such for the default kernel
%       you will get octagonal like distance function.
%
%       However using a larger radius such as "Euclidean:4" you will get a
%       much smoother distance gradient from the edge of the shape. Especially
%       if the image is pre-processed to include any anti-aliasing pixels.
%       Of course a larger kernel is slower to use, and not always needed.
%
%    The first three Distance Measuring Kernels will only generate distances
%    of exact multiples of {scale} in binary images. As such you can use a
%    scale of 1 without loosing any information.  However you also need some
%    scaling when handling non-binary anti-aliased shapes.
%
%    The "Euclidean" Distance Kernel however does generate a non-integer
%    fractional results, and as such scaling is vital even for binary shapes.
%
*/

MagickExport KernelInfo *AcquireKernelBuiltIn(const KernelInfoType type,
   const GeometryInfo *args)
{
  KernelInfo
    *kernel;

  register ssize_t
    i;

  register ssize_t
    u,
    v;

  double
    nan = sqrt((double)-1.0);  /* Special Value : Not A Number */

  /* Generate a new empty kernel if needed */
  kernel=(KernelInfo *) NULL;
  switch(type) {
    case UndefinedKernel:    /* These should not call this function */
    case UserDefinedKernel:
      assert("Should not call this function" != (char *)NULL);
      break;
    case LaplacianKernel:   /* Named Descrete Convolution Kernels */
    case SobelKernel:       /* these are defined using other kernels */
    case RobertsKernel:
    case PrewittKernel:
    case CompassKernel:
    case KirschKernel:
    case FreiChenKernel:
    case EdgesKernel:       /* Hit and Miss kernels */
    case CornersKernel:
    case DiagonalsKernel:
    case LineEndsKernel:
    case LineJunctionsKernel:
    case RidgesKernel:
    case ConvexHullKernel:
    case SkeletonKernel:
    case ThinSEKernel:
      break;               /* A pre-generated kernel is not needed */
#if 0
    /* set to 1 to do a compile-time check that we haven't missed anything */
    case UnityKernel:
    case GaussianKernel:
    case DoGKernel:
    case LoGKernel:
    case BlurKernel:
    case CometKernel:
    case DiamondKernel:
    case SquareKernel:
    case RectangleKernel:
    case OctagonKernel:
    case DiskKernel:
    case PlusKernel:
    case CrossKernel:
    case RingKernel:
    case PeaksKernel:
    case ChebyshevKernel:
    case ManhattanKernel:
    case OctangonalKernel:
    case EuclideanKernel:
#else
    default:
#endif
      /* Generate the base Kernel Structure */
      kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
      if (kernel == (KernelInfo *) NULL)
        return(kernel);
      (void) ResetMagickMemory(kernel,0,sizeof(*kernel));
      kernel->minimum = kernel->maximum = kernel->angle = 0.0;
      kernel->negative_range = kernel->positive_range = 0.0;
      kernel->type = type;
      kernel->next = (KernelInfo *) NULL;
      kernel->signature = MagickSignature;
      break;
  }

  switch(type) {
    /*
      Convolution Kernels
    */
    case UnityKernel:
      {
        kernel->height = kernel->width = (size_t) 1;
        kernel->x = kernel->y = (ssize_t) 0;
        kernel->values=(MagickRealType *) AcquireAlignedMemory(1,
          sizeof(*kernel->values));
        if (kernel->values == (MagickRealType *) NULL)
          return(DestroyKernelInfo(kernel));
        kernel->maximum = kernel->values[0] = args->rho;
        break;
      }
      break;
    case GaussianKernel:
    case DoGKernel:
    case LoGKernel:
      { double
          sigma = fabs(args->sigma),
          sigma2 = fabs(args->xi),
          A, B, R;

        if ( args->rho >= 1.0 )
          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 = (ssize_t) (kernel->width-1)/2;
        kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
          kernel->height*sizeof(*kernel->values));
        if (kernel->values == (MagickRealType *) NULL)
          return(DestroyKernelInfo(kernel));

        /* WARNING: The following generates a 'sampled gaussian' kernel.
         * What we really want is a 'discrete gaussian' kernel.
         *
         * How to do this is I don't know, but appears to be basied on the
         * Error Function 'erf()' (intergral of a gaussian)
         */

        if ( type == GaussianKernel || type == DoGKernel )
          { /* Calculate a Gaussian,  OR positive half of a DoG */
            if ( sigma > MagickEpsilon )
              { A = 1.0/(2.0*sigma*sigma);  /* simplify loop expressions */
                B = (double) (1.0/(Magick2PI*sigma*sigma));
                for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
                  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
                      kernel->values[i] = exp(-((double)(u*u+v*v))*A)*B;
              }
            else /* limiting case - a unity (normalized Dirac) kernel */
              { (void) ResetMagickMemory(kernel->values,0, (size_t)
                            kernel->width*kernel->height*sizeof(double));
                kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
              }
          }

        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 = (double) (1.0/(Magick2PI*sigma*sigma));
                for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
                  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
                    kernel->values[i] -= exp(-((double)(u*u+v*v))*A)*B;
              }
            else /* limiting case - a unity (normalized Dirac) kernel */
              kernel->values[kernel->x+kernel->y*kernel->width] -= 1.0;
          }

        if ( type == LoGKernel )
          { /* Calculate a Laplacian of a Gaussian - Or Mexician Hat */
            if ( sigma > MagickEpsilon )
              { A = 1.0/(2.0*sigma*sigma);  /* simplify loop expressions */
                B = (double) (1.0/(MagickPI*sigma*sigma*sigma*sigma));
                for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
                  for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
                    { R = ((double)(u*u+v*v))*A;
                      kernel->values[i] = (1-R)*exp(-R)*B;
                    }
              }
            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;
              }
          }

        /* Note the above kernels may have been 'clipped' by a user defined
        ** radius, producing a smaller (darker) kernel.  Also for very small
        ** sigma's (> 0.1) the central value becomes larger than one, and thus
        ** producing a very bright kernel.
        **
        ** Normalization will still be needed.
        */

        /* Normalize the 2D Gaussian Kernel
        **
        ** NB: a CorrelateNormalize performs a normal Normalize if
        ** there are no negative values.
        */
        CalcKernelMetaData(kernel);  /* the other kernel meta-data */
        ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);

        break;
      }
    case BlurKernel:
      { double
          sigma = fabs(args->sigma),
          alpha, beta;

        if ( args->rho >= 1.0 )
          kernel->width = (size_t)args->rho*2+1;
        else
          kernel->width = GetOptimalKernelWidth1D(args->rho,sigma);
        kernel->height = 1;
        kernel->x = (ssize_t) (kernel->width-1)/2;
        kernel->y = 0;
        kernel->negative_range = kernel->positive_range = 0.0;
        kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
          kernel->height*sizeof(*kernel->values));
        if (kernel->values == (MagickRealType *) NULL)
          return(DestroyKernelInfo(kernel));

#if 1
#define KernelRank 3
        /* Formula derived from GetBlurKernel() in "effect.c" (plus bug fix).
        ** It generates a gaussian 3 times the width, and compresses it into
        ** the expected range.  This produces a closer normalization of the
        ** resulting kernel, especially for very low sigma values.
        ** As such while wierd it is prefered.
        **
        ** I am told this method originally came from Photoshop.
        **
        ** A properly normalized curve is generated (apart from edge clipping)
        ** even though we later normalize the result (for edge clipping)
        ** to allow the correct generation of a "Difference of Blurs".
        */

        /* initialize */
        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 */
            alpha = 1.0/(2.0*sigma*sigma);
            beta= (double) (1.0/(MagickSQ2PI*sigma ));
            for ( u=-v; u <= v; u++) {
              kernel->values[(u+v)/KernelRank] +=
                              exp(-((double)(u*u))*alpha)*beta;
            }
          }
        else /* special case - generate a unity kernel */
          kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
#else
        /* Direct calculation without curve averaging */

        /* Calculate a Positive Gaussian */
        if ( sigma > MagickEpsilon )
          { 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;
          }
#endif
        /* Note the above kernel may have been 'clipped' by a user defined
        ** radius, producing a smaller (darker) kernel.  Also for very small
        ** sigma's (> 0.1) the central value becomes larger than one, and thus
        ** producing a very bright kernel.
        **
        ** Normalization will still be needed.
        */

        /* Normalize the 1D Gaussian Kernel
        **
        ** NB: a CorrelateNormalize performs a normal Normalize if
        ** there are no negative values.
        */
        CalcKernelMetaData(kernel);  /* the other kernel meta-data */
        ScaleKernelInfo(kernel, 1.0, CorrelateNormalizeValue);

        /* rotate the 1D kernel by given angle */
        RotateKernelInfo(kernel, args->xi );
        break;
      }
    case CometKernel:
      { double
          sigma = fabs(args->sigma),
          A;

        if ( args->rho < 1.0 )
          kernel->width = (GetOptimalKernelWidth1D(args->rho,sigma)-1)/2+1;
        else
          kernel->width = (size_t)args->rho;
        kernel->x = kernel->y = 0;
        kernel->height = 1;
        kernel->negative_range = kernel->positive_range = 0.0;
        kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
          kernel->height*sizeof(*kernel->values));
        if (kernel->values == (MagickRealType *) NULL)
          return(DestroyKernelInfo(kernel));

        /* A comet blur is half a 1D gaussian curve, so that the object is
        ** blurred in one direction only.  This may not be quite the right
        ** curve to use so may change in the future. The function must be
        ** normalised after generation, which also resolves any clipping.
        **
        ** As we are normalizing and not subtracting gaussians,
        ** there is no need for a divisor in the gaussian formula
        **
        ** It is less comples
        */
        if ( sigma > MagickEpsilon )
          {
#if 1
#define KernelRank 3
            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 */
            A = 1.0/(2.0*sigma*sigma);
            /* B = 1.0/(MagickSQ2PI*sigma); */
            for ( u=0; u < v; u++) {
              kernel->values[u/KernelRank] +=
                  exp(-((double)(u*u))*A);
              /*  exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
            }
            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 < (ssize_t) kernel->width; i++)
              kernel->positive_range +=
                kernel->values[i] =
                  exp(-((double)(i*i))*A);
                /* exp(-((double)(i*i))/2.0*sigma*sigma)/(MagickSQ2PI*sigma); */
#endif
          }
        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;
            kernel->positive_range = 1.0;
          }

        kernel->minimum = 0.0;
        kernel->maximum = kernel->values[0];
        kernel->negative_range = 0.0;

        ScaleKernelInfo(kernel, 1.0, NormalizeValue); /* Normalize */
        RotateKernelInfo(kernel, args->xi); /* Rotate by angle */
        break;
      }

    /*
      Convolution Kernels - Well Known Named Constant Kernels
    */
    case LaplacianKernel:
      { switch ( (int) args->rho ) {
          case 0:
          default: /* laplacian square filter -- default */
            kernel=ParseKernelArray("3: -1,-1,-1  -1,8,-1  -1,-1,-1");
            break;
          case 1:  /* laplacian diamond filter */
            kernel=ParseKernelArray("3: 0,-1,0  -1,4,-1  0,-1,0");
            break;
          case 2:
            kernel=ParseKernelArray("3: -2,1,-2  1,4,1  -2,1,-2");
            break;
          case 3:
            kernel=ParseKernelArray("3: 1,-2,1  -2,4,-2  1,-2,1");
            break;
          case 5:   /* a 5x5 laplacian */
            kernel=ParseKernelArray(
              "5: -4,-1,0,-1,-4  -1,2,3,2,-1  0,3,4,3,0  -1,2,3,2,-1  -4,-1,0,-1,-4");
            break;
          case 7:   /* a 7x7 laplacian */
            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) */
            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) */
            /* 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");
            break;
        }
        if (kernel == (KernelInfo *) NULL)
          return(kernel);
        kernel->type = type;
        break;
      }
    case SobelKernel:
      { /* Simple Sobel Kernel */
        kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
        if (kernel == (KernelInfo *) NULL)
          return(kernel);
        kernel->type = type;
        RotateKernelInfo(kernel, args->rho);
        break;
      }
    case RobertsKernel:
      {
        kernel=ParseKernelArray("3: 0,0,0  1,-1,0  0,0,0");
        if (kernel == (KernelInfo *) NULL)
          return(kernel);
        kernel->type = type;
        RotateKernelInfo(kernel, args->rho);
        break;
      }
    case PrewittKernel:
      {
        kernel=ParseKernelArray("3: 1,0,-1  1,0,-1  1,0,-1");
        if (kernel == (KernelInfo *) NULL)
          return(kernel);
        kernel->type = type;
        RotateKernelInfo(kernel, args->rho);
        break;
      }
    case CompassKernel:
      {
        kernel=ParseKernelArray("3: 1,1,-1  1,-2,-1  1,1,-1");
        if (kernel == (KernelInfo *) NULL)
          return(kernel);
        kernel->type = type;
        RotateKernelInfo(kernel, args->rho);
        break;
      }
    case KirschKernel:
      {
        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, (double) (1.0/2.0*MagickSQ2), NoValue);
            break;
          case 10:
            kernel=AcquireKernelInfo("FreiChen:11;FreiChen:12;FreiChen:13;FreiChen:14;FreiChen:15;FreiChen:16;FreiChen:17;FreiChen:18;FreiChen:19");
            if (kernel == (KernelInfo *) NULL)
              return(kernel);
            break;
          case 1:
          case 11:
            kernel=ParseKernelArray("3: 1,0,-1  2,0,-2  1,0,-1");
            if (kernel == (KernelInfo *) NULL)
              return(kernel);
            kernel->type = type;
            kernel->values[3] = +MagickSQ2;
            kernel->values[5] = -MagickSQ2;
            CalcKernelMetaData(kernel);     /* recalculate meta-data */
            ScaleKernelInfo(kernel, (double) (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, (double) (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, (double) (1.0/2.0*MagickSQ2), NoValue);
            break;
          case 14:
            kernel=ParseKernelArray("3: 0,1,-2  -1,0,1  2,-1,0");
            if (kernel == (KernelInfo *) NULL)
              return(kernel);
            kernel->type = type;
            kernel->values[2] = -MagickSQ2;
            kernel->values[6] = +MagickSQ2;
            CalcKernelMetaData(kernel);
            ScaleKernelInfo(kernel, (double) (1.0/2.0*MagickSQ2), NoValue);
            break;
          case 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 or Shaped Kernels
    */
    case DiamondKernel:
      {
        if (args->rho < 1.0)
          kernel->width = kernel->height = 3;  /* default radius = 1 */
        else
          kernel->width = kernel->height = ((size_t)args->rho)*2+1;
        kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;

        kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
          kernel->height*sizeof(*kernel->values));
        if (kernel->values == (MagickRealType *) NULL)
          return(DestroyKernelInfo(kernel));

        /* set all kernel values within diamond area to scale given */
        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;
        kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
        break;
      }
    case SquareKernel:
    case RectangleKernel:
      { double
          scale;
        if ( type == SquareKernel )
          {
            if (args->rho < 1.0)
              kernel->width = kernel->height = 3;  /* default radius = 1 */
            else
              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 = (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 = (ssize_t) args->xi;
            kernel->y = (ssize_t) args->psi;
            scale = 1.0;
          }
        kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
          kernel->height*sizeof(*kernel->values));
        if (kernel->values == (MagickRealType *) NULL)
          return(DestroyKernelInfo(kernel));

        /* set all kernel values to scale given */
        u=(ssize_t) (kernel->width*kernel->height);
        for ( i=0; i < u; i++)
            kernel->values[i] = scale;
        kernel->minimum = kernel->maximum = scale;   /* a flat shape */
        kernel->positive_range = scale*u;
        break;
      }
      case OctagonKernel:
        {
          if (args->rho < 1.0)
            kernel->width = kernel->height = 5;  /* default radius = 2 */
          else
            kernel->width = kernel->height = ((size_t)args->rho)*2+1;
          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;

          kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
            kernel->height*sizeof(*kernel->values));
          if (kernel->values == (MagickRealType *) NULL)
            return(DestroyKernelInfo(kernel));

          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
              if ( (labs((long) u)+labs((long) v)) <=
                        ((long)kernel->x + (long)(kernel->x/2)) )
                kernel->positive_range += kernel->values[i] = args->sigma;
              else
                kernel->values[i] = nan;
          kernel->minimum = kernel->maximum = args->sigma;  /* a flat shape */
          break;
        }
      case DiskKernel:
        {
          ssize_t
            limit = (ssize_t)(args->rho*args->rho);

          if (args->rho < 0.4)           /* default radius approx 4.3 */
            kernel->width = kernel->height = 9L, limit = 18L;
          else
            kernel->width = kernel->height = (size_t)fabs(args->rho)*2+1;
          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;

          kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
            kernel->height*sizeof(*kernel->values));
          if (kernel->values == (MagickRealType *) NULL)
            return(DestroyKernelInfo(kernel));

          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
              if ((u*u+v*v) <= limit)
                kernel->positive_range += kernel->values[i] = args->sigma;
              else
                kernel->values[i] = nan;
          kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
          break;
        }
      case PlusKernel:
        {
          if (args->rho < 1.0)
            kernel->width = kernel->height = 5;  /* default radius 2 */
          else
            kernel->width = kernel->height = ((size_t)args->rho)*2+1;
          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;

          kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
            kernel->height*sizeof(*kernel->values));
          if (kernel->values == (MagickRealType *) NULL)
            return(DestroyKernelInfo(kernel));

          /* set all kernel values along axises to given scale */
          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
              kernel->values[i] = (u == 0 || v == 0) ? args->sigma : nan;
          kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
          kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
          break;
        }
      case CrossKernel:
        {
          if (args->rho < 1.0)
            kernel->width = kernel->height = 5;  /* default radius 2 */
          else
            kernel->width = kernel->height = ((size_t)args->rho)*2+1;
          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;

          kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
            kernel->height*sizeof(*kernel->values));
          if (kernel->values == (MagickRealType *) NULL)
            return(DestroyKernelInfo(kernel));

          /* set all kernel values along axises to given scale */
          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
              kernel->values[i] = (u == v || u == -v) ? args->sigma : nan;
          kernel->minimum = kernel->maximum = args->sigma;   /* a flat shape */
          kernel->positive_range = args->sigma*(kernel->width*2.0 - 1.0);
          break;
        }
      /*
        HitAndMiss Kernels
      */
      case RingKernel:
      case PeaksKernel:
        {
          ssize_t
            limit1,
            limit2,
            scale;

          if (args->rho < args->sigma)
            {
              kernel->width = ((size_t)args->sigma)*2+1;
              limit1 = (ssize_t)(args->rho*args->rho);
              limit2 = (ssize_t)(args->sigma*args->sigma);
            }
          else
            {
              kernel->width = ((size_t)args->rho)*2+1;
              limit1 = (ssize_t)(args->sigma*args->sigma);
              limit2 = (ssize_t)(args->rho*args->rho);
            }
          if ( limit2 <= 0 )
            kernel->width = 7L, limit1 = 7L, limit2 = 11L;

          kernel->height = kernel->width;
          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;
          kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
            kernel->height*sizeof(*kernel->values));
          if (kernel->values == (MagickRealType *) NULL)
            return(DestroyKernelInfo(kernel));

          /* set a ring of points of 'scale' ( 0.0 for PeaksKernel ) */
          scale = (ssize_t) (( type == PeaksKernel) ? 0.0 : args->xi);
          for ( i=0, v= -kernel->y; v <= (ssize_t)kernel->y; v++)
            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
              { ssize_t radius=u*u+v*v;
                if (limit1 < radius && radius <= limit2)
                  kernel->positive_range += kernel->values[i] = (double) scale;
                else
                  kernel->values[i] = nan;
              }
          kernel->minimum = kernel->maximum = (double) scale;
          if ( type == PeaksKernel ) {
            /* set the central point in the middle */
            kernel->values[kernel->x+kernel->y*kernel->width] = 1.0;
            kernel->positive_range = 1.0;
            kernel->maximum = 1.0;
          }
          break;
        }
      case EdgesKernel:
        {
          kernel=AcquireKernelInfo("ThinSE:482");
          if (kernel == (KernelInfo *) NULL)
            return(kernel);
          kernel->type = type;
          ExpandMirrorKernelInfo(kernel); /* mirror expansion of kernels */
          break;
        }
      case CornersKernel:
        {
          kernel=AcquireKernelInfo("ThinSE:87");
          if (kernel == (KernelInfo *) NULL)
            return(kernel);
          kernel->type = type;
          ExpandRotateKernelInfo(kernel, 90.0); /* Expand 90 degree rotations */
          break;
        }
      case DiagonalsKernel:
        {
          switch ( (int) args->rho ) {
            case 0:
            default:
              { KernelInfo
                  *new_kernel;
                kernel=ParseKernelArray("3: 0,0,0  0,-,1  1,1,-");
                if (kernel == (KernelInfo *) NULL)
                  return(kernel);
                kernel->type = type;
                new_kernel=ParseKernelArray("3: 0,0,1  0,-,1  0,1,-");
                if (new_kernel == (KernelInfo *) NULL)
                  return(DestroyKernelInfo(kernel));
                new_kernel->type = type;
                LastKernelInfo(kernel)->next = new_kernel;
                ExpandMirrorKernelInfo(kernel);
                return(kernel);
              }
            case 1:
              kernel=ParseKernelArray("3: 0,0,0  0,-,1  1,1,-");
              break;
            case 2:
              kernel=ParseKernelArray("3: 0,0,1  0,-,1  0,1,-");
              break;
          }
          if (kernel == (KernelInfo *) NULL)
            return(kernel);
          kernel->type = type;
          RotateKernelInfo(kernel, args->sigma);
          break;
        }
      case LineEndsKernel:
        { /* Kernels for finding the end of thin lines */
          switch ( (int) args->rho ) {
            case 0:
            default:
              /* set of kernels to find all end of lines */
              return(AcquireKernelInfo("LineEnds:1>;LineEnds:2>"));
            case 1:
              /* kernel for 4-connected line ends - no rotation */
              kernel=ParseKernelArray("3: 0,0,-  0,1,1  0,0,-");
              break;
          case 2:
              /* kernel to add for 8-connected lines - no rotation */
              kernel=ParseKernelArray("3: 0,0,0  0,1,0  0,0,1");
              break;
          case 3:
              /* kernel to add for orthogonal line ends - does not find corners */
              kernel=ParseKernelArray("3: 0,0,0  0,1,1  0,0,0");
              break;
          case 4:
              /* traditional line end - fails on last T end */
              kernel=ParseKernelArray("3: 0,0,0  0,1,-  0,0,-");
              break;
          }
          if (kernel == (KernelInfo *) NULL)
            return(kernel);
          kernel->type = type;
          RotateKernelInfo(kernel, args->sigma);
          break;
        }
      case LineJunctionsKernel:
        { /* kernels for finding the junctions of multiple lines */
          switch ( (int) args->rho ) {
            case 0:
            default:
              /* set of kernels to find all line junctions */
              return(AcquireKernelInfo("LineJunctions:1@;LineJunctions:2>"));
            case 1:
              /* Y Junction */
              kernel=ParseKernelArray("3: 1,-,1  -,1,-  -,1,-");
              break;
            case 2:
              /* Diagonal T Junctions */
              kernel=ParseKernelArray("3: 1,-,-  -,1,-  1,-,1");
              break;
            case 3:
              /* Orthogonal T Junctions */
              kernel=ParseKernelArray("3: -,-,-  1,1,1  -,1,-");
              break;
            case 4:
              /* Diagonal X Junctions */
              kernel=ParseKernelArray("3: 1,-,1  -,1,-  1,-,1");
              break;
            case 5:
              /* Orthogonal X Junctions - minimal diamond kernel */
              kernel=ParseKernelArray("3: -,1,-  1,1,1  -,1,-");
              break;
          }
          if (kernel == (KernelInfo *) NULL)
            return(kernel);
          kernel->type = type;
          RotateKernelInfo(kernel, args->sigma);
          break;
        }
      case RidgesKernel:
        { /* Ridges - Ridge finding kernels */
          KernelInfo
            *new_kernel;
          switch ( (int) args->rho ) {
            case 1:
            default:
              kernel=ParseKernelArray("3x1:0,1,0");
              if (kernel == (KernelInfo *) NULL)
                return(kernel);
              kernel->type = type;
              ExpandRotateKernelInfo(kernel, 90.0); /* 2 rotated kernels (symmetrical) */
              break;
            case 2:
              kernel=ParseKernelArray("4x1:0,1,1,0");
              if (kernel == (KernelInfo *) NULL)
                return(kernel);
              kernel->type = type;
              ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotated kernels */

              /* Kernels to find a stepped 'thick' line, 4 rotates + mirrors */
              /* Unfortunatally we can not yet rotate a non-square kernel */
              /* But then we can't flip a non-symetrical kernel either */
              new_kernel=ParseKernelArray("4x3+1+1:0,1,1,- -,1,1,- -,1,1,0");
              if (new_kernel == (KernelInfo *) NULL)
                return(DestroyKernelInfo(kernel));
              new_kernel->type = type;
              LastKernelInfo(kernel)->next = new_kernel;
              new_kernel=ParseKernelArray("4x3+2+1:0,1,1,- -,1,1,- -,1,1,0");
              if (new_kernel == (KernelInfo *) NULL)
                return(DestroyKernelInfo(kernel));
              new_kernel->type = type;
              LastKernelInfo(kernel)->next = new_kernel;
              new_kernel=ParseKernelArray("4x3+1+1:-,1,1,0 -,1,1,- 0,1,1,-");
              if (new_kernel == (KernelInfo *) NULL)
                return(DestroyKernelInfo(kernel));
              new_kernel->type = type;
              LastKernelInfo(kernel)->next = new_kernel;
              new_kernel=ParseKernelArray("4x3+2+1:-,1,1,0 -,1,1,- 0,1,1,-");
              if (new_kernel == (KernelInfo *) NULL)
                return(DestroyKernelInfo(kernel));
              new_kernel->type = type;
              LastKernelInfo(kernel)->next = new_kernel;
              new_kernel=ParseKernelArray("3x4+1+1:0,-,- 1,1,1 1,1,1 -,-,0");
              if (new_kernel == (KernelInfo *) NULL)
                return(DestroyKernelInfo(kernel));
              new_kernel->type = type;
              LastKernelInfo(kernel)->next = new_kernel;
              new_kernel=ParseKernelArray("3x4+1+2:0,-,- 1,1,1 1,1,1 -,-,0");
              if (new_kernel == (KernelInfo *) NULL)
                return(DestroyKernelInfo(kernel));
              new_kernel->type = type;
              LastKernelInfo(kernel)->next = new_kernel;
              new_kernel=ParseKernelArray("3x4+1+1:-,-,0 1,1,1 1,1,1 0,-,-");
              if (new_kernel == (KernelInfo *) NULL)
                return(DestroyKernelInfo(kernel));
              new_kernel->type = type;
              LastKernelInfo(kernel)->next = new_kernel;
              new_kernel=ParseKernelArray("3x4+1+2:-,-,0 1,1,1 1,1,1 0,-,-");
              if (new_kernel == (KernelInfo *) NULL)
                return(DestroyKernelInfo(kernel));
              new_kernel->type = type;
              LastKernelInfo(kernel)->next = new_kernel;
              break;
          }
          break;
        }
      case ConvexHullKernel:
        {
          KernelInfo
            *new_kernel;
          /* first set of 8 kernels */
          kernel=ParseKernelArray("3: 1,1,-  1,0,-  1,-,0");
          if (kernel == (KernelInfo *) NULL)
            return(kernel);
          kernel->type = type;
          ExpandRotateKernelInfo(kernel, 90.0);
          /* append the mirror versions too - no flip function yet */
          new_kernel=ParseKernelArray("3: 1,1,1  1,0,-  -,-,0");
          if (new_kernel == (KernelInfo *) NULL)
            return(DestroyKernelInfo(kernel));
          new_kernel->type = type;
          ExpandRotateKernelInfo(new_kernel, 90.0);
          LastKernelInfo(kernel)->next = new_kernel;
          break;
        }
      case SkeletonKernel:
        {
          switch ( (int) args->rho ) {
            case 1:
            default:
              /* Traditional Skeleton...
              ** A cyclically rotated single kernel
              */
              kernel=AcquireKernelInfo("ThinSE:482");
              if (kernel == (KernelInfo *) NULL)
                return(kernel);
              kernel->type = type;
              ExpandRotateKernelInfo(kernel, 45.0); /* 8 rotations */
              break;
            case 2:
              /* HIPR Variation of the cyclic skeleton
              ** Corners of the traditional method made more forgiving,
              ** but the retain the same cyclic order.
              */
              kernel=AcquireKernelInfo("ThinSE:482; ThinSE:87x90;");
              if (kernel == (KernelInfo *) NULL)
                return(kernel);
              if (kernel->next == (KernelInfo *) NULL)
                return(DestroyKernelInfo(kernel));
              kernel->type = type;
              kernel->next->type = type;
              ExpandRotateKernelInfo(kernel, 90.0); /* 4 rotations of the 2 kernels */
              break;
            case 3:
              /* Dan Bloomberg Skeleton, from his paper on 3x3 thinning SE's
              ** "Connectivity-Preserving Morphological Image Thransformations"
              ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
              **   http://www.leptonica.com/papers/conn.pdf
              */
              kernel=AcquireKernelInfo(
                            "ThinSE:41; ThinSE:42; ThinSE:43");
              if (kernel == (KernelInfo *) NULL)
                return(kernel);
              kernel->type = type;
              kernel->next->type = type;
              kernel->next->next->type = type;
              ExpandMirrorKernelInfo(kernel); /* 12 kernels total */
              break;
           }
          break;
        }
      case ThinSEKernel:
        { /* Special kernels for general thinning, while preserving connections
          ** "Connectivity-Preserving Morphological Image Thransformations"
          ** by Dan S. Bloomberg, available on Leptonica, Selected Papers,
          **   http://www.leptonica.com/papers/conn.pdf
          ** And
          **   http://tpgit.github.com/Leptonica/ccthin_8c_source.html
          **
          ** Note kernels do not specify the origin pixel, allowing them
          ** to be used for both thickening and thinning operations.
          */
          switch ( (int) args->rho ) {
            /* SE for 4-connected thinning */
            case 41: /* SE_4_1 */
              kernel=ParseKernelArray("3: -,-,1  0,-,1  -,-,1");
              break;
            case 42: /* SE_4_2 */
              kernel=ParseKernelArray("3: -,-,1  0,-,1  -,0,-");
              break;
            case 43: /* SE_4_3 */
              kernel=ParseKernelArray("3: -,0,-  0,-,1  -,-,1");
              break;
            case 44: /* SE_4_4 */
              kernel=ParseKernelArray("3: -,0,-  0,-,1  -,0,-");
              break;
            case 45: /* SE_4_5 */
              kernel=ParseKernelArray("3: -,0,1  0,-,1  -,0,-");
              break;
            case 46: /* SE_4_6 */
              kernel=ParseKernelArray("3: -,0,-  0,-,1  -,0,1");
              break;
            case 47: /* SE_4_7 */
              kernel=ParseKernelArray("3: -,1,1  0,-,1  -,0,-");
              break;
            case 48: /* SE_4_8 */
              kernel=ParseKernelArray("3: -,-,1  0,-,1  0,-,1");
              break;
            case 49: /* SE_4_9 */
              kernel=ParseKernelArray("3: 0,-,1  0,-,1  -,-,1");
              break;
            /* SE for 8-connected thinning - negatives of the above */
            case 81: /* SE_8_0 */
              kernel=ParseKernelArray("3: -,1,-  0,-,1  -,1,-");
              break;
            case 82: /* SE_8_2 */
              kernel=ParseKernelArray("3: -,1,-  0,-,1  0,-,-");
              break;
            case 83: /* SE_8_3 */
              kernel=ParseKernelArray("3: 0,-,-  0,-,1  -,1,-");
              break;
            case 84: /* SE_8_4 */
              kernel=ParseKernelArray("3: 0,-,-  0,-,1  0,-,-");
              break;
            case 85: /* SE_8_5 */
              kernel=ParseKernelArray("3: 0,-,1  0,-,1  0,-,-");
              break;
            case 86: /* SE_8_6 */
              kernel=ParseKernelArray("3: 0,-,-  0,-,1  0,-,1");
              break;
            case 87: /* SE_8_7 */
              kernel=ParseKernelArray("3: -,1,-  0,-,1  0,0,-");
              break;
            case 88: /* SE_8_8 */
              kernel=ParseKernelArray("3: -,1,-  0,-,1  0,1,-");
              break;
            case 89: /* SE_8_9 */
              kernel=ParseKernelArray("3: 0,1,-  0,-,1  -,1,-");
              break;
            /* Special combined SE kernels */
            case 423: /* SE_4_2 , SE_4_3 Combined Kernel */
              kernel=ParseKernelArray("3: -,-,1  0,-,-  -,0,-");
              break;
            case 823: /* SE_8_2 , SE_8_3 Combined Kernel */
              kernel=ParseKernelArray("3: -,1,-  -,-,1  0,-,-");
              break;
            case 481: /* SE_48_1 - General Connected Corner Kernel */
              kernel=ParseKernelArray("3: -,1,1  0,-,1  0,0,-");
              break;
            default:
            case 482: /* SE_48_2 - General Edge Kernel */
              kernel=ParseKernelArray("3: 0,-,1  0,-,1  0,-,1");
              break;
          }
          if (kernel == (KernelInfo *) NULL)
            return(kernel);
          kernel->type = type;
          RotateKernelInfo(kernel, args->sigma);
          break;
        }
      /*
        Distance Measuring Kernels
      */
      case ChebyshevKernel:
        {
          if (args->rho < 1.0)
            kernel->width = kernel->height = 3;  /* default radius = 1 */
          else
            kernel->width = kernel->height = ((size_t)args->rho)*2+1;
          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;

          kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
            kernel->height*sizeof(*kernel->values));
          if (kernel->values == (MagickRealType *) NULL)
            return(DestroyKernelInfo(kernel));

          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
              kernel->positive_range += ( kernel->values[i] =
                  args->sigma*MagickMax(fabs((double)u),fabs((double)v)) );
          kernel->maximum = kernel->values[0];
          break;
        }
      case ManhattanKernel:
        {
          if (args->rho < 1.0)
            kernel->width = kernel->height = 3;  /* default radius = 1 */
          else
            kernel->width = kernel->height = ((size_t)args->rho)*2+1;
          kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;

          kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
            kernel->height*sizeof(*kernel->values));
          if (kernel->values == (MagickRealType *) NULL)
            return(DestroyKernelInfo(kernel));

          for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
            for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
              kernel->positive_range += ( kernel->values[i] =
                  args->sigma*(labs((long) u)+labs((long) v)) );
          kernel->maximum = kernel->values[0];
          break;
        }
      case OctagonalKernel:
      {
        if (args->rho < 2.0)
          kernel->width = kernel->height = 5;  /* default/minimum radius = 2 */
        else
          kernel->width = kernel->height = ((size_t)args->rho)*2+1;
        kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;

        kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
          kernel->height*sizeof(*kernel->values));
        if (kernel->values == (MagickRealType *) NULL)
          return(DestroyKernelInfo(kernel));

        for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
          for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
            {
              double
                r1 = MagickMax(fabs((double)u),fabs((double)v)),
                r2 = floor((double)(labs((long)u)+labs((long)v)+1)/1.5);
              kernel->positive_range += kernel->values[i] =
                        args->sigma*MagickMax(r1,r2);
            }
        kernel->maximum = kernel->values[0];
        break;
      }
    case EuclideanKernel:
      {
        if (args->rho < 1.0)
          kernel->width = kernel->height = 3;  /* default radius = 1 */
        else
          kernel->width = kernel->height = ((size_t)args->rho)*2+1;
        kernel->x = kernel->y = (ssize_t) (kernel->width-1)/2;

        kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
          kernel->height*sizeof(*kernel->values));
        if (kernel->values == (MagickRealType *) NULL)
          return(DestroyKernelInfo(kernel));

        for ( i=0, v=-kernel->y; v <= (ssize_t)kernel->y; v++)
          for ( u=-kernel->x; u <= (ssize_t)kernel->x; u++, i++)
            kernel->positive_range += ( kernel->values[i] =
                 args->sigma*sqrt((double)(u*u+v*v)) );
        kernel->maximum = kernel->values[0];
        break;
      }
    default:
      {
        /* No-Op Kernel - Basically just a single pixel on its own */
        kernel=ParseKernelArray("1:1");
        if (kernel == (KernelInfo *) NULL)
          return(kernel);
        kernel->type = UndefinedKernel;
        break;
      }
      break;
  }
  return(kernel);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     C l o n e K e r n e l I n f o                                           %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  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 DestroyKernelInfo() when no longer needed.
%
%  The format of the CloneKernelInfo method is:
%
%      KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
%
%  A description of each parameter follows:
%
%    o kernel: the Morphology/Convolution kernel to be cloned
%
*/
MagickExport KernelInfo *CloneKernelInfo(const KernelInfo *kernel)
{
  register ssize_t
    i;

  KernelInfo
    *new_kernel;

  assert(kernel != (KernelInfo *) NULL);
  new_kernel=(KernelInfo *) AcquireMagickMemory(sizeof(*kernel));
  if (new_kernel == (KernelInfo *) NULL)
    return(new_kernel);
  *new_kernel=(*kernel); /* copy values in structure */

  /* replace the values with a copy of the values */
  new_kernel->values=(MagickRealType *) AcquireAlignedMemory(kernel->width,
    kernel->height*sizeof(*kernel->values));
  if (new_kernel->values == (MagickRealType *) NULL)
    return(DestroyKernelInfo(new_kernel));
  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 */
  if ( kernel->next != (KernelInfo *) NULL ) {
    new_kernel->next = CloneKernelInfo(kernel->next);
    if ( new_kernel->next == (KernelInfo *) NULL )
      return(DestroyKernelInfo(new_kernel));
  }

  return(new_kernel);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     D e s t r o y K e r n e l I n f o                                       %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  DestroyKernelInfo() frees the memory used by a Convolution/Morphology
%  kernel.
%
%  The format of the DestroyKernelInfo method is:
%
%      KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
%
%  A description of each parameter follows:
%
%    o kernel: the Morphology/Convolution kernel to be destroyed
%
*/
MagickExport KernelInfo *DestroyKernelInfo(KernelInfo *kernel)
{
  assert(kernel != (KernelInfo *) NULL);
  if ( kernel->next != (KernelInfo *) NULL )
    kernel->next=DestroyKernelInfo(kernel->next);
  kernel->values=(MagickRealType *) RelinquishAlignedMemory(kernel->values);
  kernel=(KernelInfo *) RelinquishMagickMemory(kernel);
  return(kernel);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+     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;
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+     E x p a n d R o t a t e K e r n e l I n f o                             %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ExpandRotateKernelInfo() takes a kernel list, and expands it by rotating
%  incrementally by the angle given, until the kernel repeats.
%
%  WARNING: 45 degree rotations only works for 3x3 kernels.
%  While 90 degree roatations only works for linear and square kernels
%
%  The format of the ExpandRotateKernelInfo method is:
%
%      void ExpandRotateKernelInfo(KernelInfo *kernel, double angle)
%
%  A description of each parameter follows:
%
%    o kernel: the Morphology/Convolution kernel
%
%    o angle: angle to rotate in degrees
%
% 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.
*/

/* 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 equivalence */
    if ( IsNan(kernel1->values[i]) && !IsNan(kernel2->values[i]) )
      return MagickFalse;
    if ( IsNan(kernel2->values[i]) && !IsNan(kernel1->values[i]) )
      return MagickFalse;
    /* Test actual values are equivalent */
    if ( fabs(kernel1->values[i] - kernel2->values[i]) > MagickEpsilon )
      return MagickFalse;
  }

  return MagickTrue;
}

static void ExpandRotateKernelInfo(KernelInfo *kernel, const double angle)
{
  KernelInfo
    *clone,
    *last;

  last = kernel;
  while(1) {
    clone = CloneKernelInfo(last);
    RotateKernelInfo(clone, angle);
    if ( SameKernelInfo(kernel, clone) == MagickTrue )
      break;
    LastKernelInfo(last)->next = clone;
    last = clone;
  }
  clone = DestroyKernelInfo(clone); /* kernel has repeated - junk the clone */
  return;
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+     C a l c M e t a K e r n a l I n f o                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  CalcKernelMetaData() recalculate the KernelInfo meta-data of this kernel only,
%  using the kernel values.  This should only ne used if it is not possible to
%  calculate that meta-data in some easier way.
%
%  It is important that the meta-data is correct before ScaleKernelInfo() is
%  used to perform kernel normalization.
%
%  The format of the CalcKernelMetaData method is:
%
%      void CalcKernelMetaData(KernelInfo *kernel, const double scale )
%
%  A description of each parameter follows:
%
%    o kernel: the Morphology/Convolution kernel to modify
%
%  WARNING: Minimum and Maximum values are assumed to include zero, even if
%  zero is not part of the kernel (as in Gaussian Derived kernels). This
%  however is not true for flat-shaped morphological kernels.
%
%  WARNING: Only the specific kernel pointed to is modified, not a list of
%  multiple kernels.
%
% This is an internal function and not expected to be useful outside this
% module.  This could change however.
*/
static void CalcKernelMetaData(KernelInfo *kernel)
{
  register size_t
    i;

  kernel->minimum = kernel->maximum = 0.0;
  kernel->negative_range = kernel->positive_range = 0.0;
  for (i=0; i < (kernel->width*kernel->height); i++)
    {
      if ( fabs(kernel->values[i]) < MagickEpsilon )
        kernel->values[i] = 0.0;
      ( kernel->values[i] < 0)
          ?  ( kernel->negative_range += kernel->values[i] )
          :  ( kernel->positive_range += kernel->values[i] );
      Minimize(kernel->minimum, kernel->values[i]);
      Maximize(kernel->maximum, kernel->values[i]);
    }

  return;
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     M o r p h o l o g y A p p l y                                           %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  MorphologyApply() applies a morphological method, multiple times using
%  a list of multiple kernels.
%
%  It is basically equivalent to as MorphologyImage() (see below) but
%  without any user controls.  This allows internel programs to use this
%  function, to actually perform a specific task without possible interference
%  by any API user supplied settings.
%
%  It is MorphologyImage() 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 (setting), nor is the convolve bias
%  (-bias setting or image->bias) loooked at, but must be supplied from the
%  function arguments.
%
%  The format of the MorphologyApply method is:
%
%      Image *MorphologyApply(const Image *image,MorphologyMethod method,
%        const ssize_t iterations,const KernelInfo *kernel,
%        const CompositeMethod compose,const double bias,
%        ExceptionInfo *exception)
%
%  A description of each parameter follows:
%
%    o image: the source image
%
%    o method: the morphology method to be applied.
%
%    o iterations: apply the operation this many times (or no change).
%                  A value of -1 means loop until no change found.
%                  How this is applied may depend on the morphology method.
%                  Typically this is a value of 1.
%
%    o channel: the channel type.
%
%    o kernel: An array of double representing the morphology kernel.
%
%    o compose: How to handle or merge multi-kernel results.
%          If 'UndefinedCompositeOp' use default for the Morphology method.
%          If 'NoCompositeOp' force image to be re-iterated by each kernel.
%          Otherwise merge the results using the compose method given.
%
%    o bias: Convolution Output Bias.
%
%    o exception: return any errors or warnings in this structure.
%
*/

/* Apply a Morphology Primative to an image using the given kernel.
** Two pre-created images must be provided, and no image is created.
** It returns the number of pixels that changed between the images
** for result convergence determination.
*/
static ssize_t MorphologyPrimitive(const Image *image,Image *morphology_image,
  const MorphologyMethod method,const KernelInfo *kernel,const double bias,
  ExceptionInfo *exception)
{
#define MorphologyTag  "Morphology/Image"

  CacheView
    *image_view,
    *morphology_view;

  ssize_t
    y, offx, offy;

  size_t
    virt_width,
    changed;

  MagickBooleanType
    status;

  MagickOffsetType
    progress;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  assert(morphology_image != (Image *) NULL);
  assert(morphology_image->signature == MagickSignature);
  assert(kernel != (KernelInfo *) NULL);
  assert(kernel->signature == MagickSignature);
  assert(exception != (ExceptionInfo *) NULL);
  assert(exception->signature == MagickSignature);

  status=MagickTrue;
  changed=0;
  progress=0;

  image_view=AcquireCacheView(image);
  morphology_view=AcquireCacheView(morphology_image);
  virt_width=image->columns+kernel->width-1;

  /* Some methods (including convolve) needs use a reflected kernel.
   * Adjust 'origin' offsets to loop though kernel as a reflection.
   */
  offx = kernel->x;
  offy = kernel->y;
  switch(method) {
    case ConvolveMorphology:
    case DilateMorphology:
    case DilateIntensityMorphology:
    /*case DistanceMorphology:*/
      /* kernel needs to used with reflection about origin */
      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 used as is, without reflection */
      break;
    default:
      assert("Not a Primitive Morphology Method" != (char *) NULL);
      break;
  }

  if ( method == ConvolveMorphology && kernel->width == 1 )
  { /* Special handling (for speed) of vertical (blur) kernels.
    ** This performs its handling in columns rather than in rows.
    ** This is only done for convolve as it is the only method that
    ** generates very large 1-D vertical kernels (such as a 'BlurKernel')
    **
    ** Timing tests (on single CPU laptop)
    ** Using a vertical 1-d Blue with normal row-by-row (below)
    **   time convert logo: -morphology Convolve Blur:0x10+90 null:
    **      0.807u
    ** Using this column method
    **   time convert logo: -morphology Convolve Blur:0x10+90 null:
    **      0.620u
    **
    ** Anthony Thyssen, 14 June 2010
    */
    register ssize_t
      x;

#if defined(MAGICKCORE_OPENMP_SUPPORT)
#pragma omp parallel for schedule(static,4) shared(progress,status)
#endif
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      register const Quantum
        *restrict p;

      register Quantum
        *restrict q;

      register ssize_t
        y;

      ssize_t
        r;

      if (status == MagickFalse)
        continue;
      p=GetCacheViewVirtualPixels(image_view,x,-offy,1,image->rows+
        kernel->height-1,exception);
      q=GetCacheViewAuthenticPixels(morphology_view,x,0,1,
        morphology_image->rows,exception);
      if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
        {
          status=MagickFalse;
          continue;
        }
      /* offset to origin in 'p'. while 'q' points to it directly */
      r = offy;

      for (y=0; y < (ssize_t) image->rows; y++)
      {
        PixelInfo
          result;

        register ssize_t
          v;

        register const MagickRealType
          *restrict k;

        register const Quantum
          *restrict k_pixels;

        /* Copy input image to the output image for unused channels
        * This removes need for 'cloning' a new image every iteration
        */
        SetPixelRed(morphology_image,GetPixelRed(image,p+r*
          GetPixelChannels(image)),q);
        SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
          GetPixelChannels(image)),q);
        SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
          GetPixelChannels(image)),q);
        if (image->colorspace == CMYKColorspace)
          SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
            GetPixelChannels(image)),q);

        /* Set the bias of the weighted average output */
        result.red     =
        result.green   =
        result.blue    =
        result.alpha =
        result.black   = bias;


        /* Weighted Average of pixels using reflected kernel
        **
        ** NOTE for correct working of this operation for asymetrical
        ** kernels, the kernel needs to be applied in its reflected form.
        ** That is its values needs to be reversed.
        */
        k = &kernel->values[ kernel->height-1 ];
        k_pixels = p;
        if ( (image->channel_mask != DefaultChannels) || (image->matte == MagickFalse) )
          { /* No 'Sync' involved.
            ** Convolution is simple greyscale channel operation
            */
            for (v=0; v < (ssize_t) kernel->height; v++) {
              if ( IsNan(*k) ) continue;
              result.red     += (*k)*GetPixelRed(image,k_pixels);
              result.green   += (*k)*GetPixelGreen(image,k_pixels);
              result.blue    += (*k)*GetPixelBlue(image,k_pixels);
              if (image->colorspace == CMYKColorspace)
                result.black+=(*k)*GetPixelBlack(image,k_pixels);
              result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
              k--;
              k_pixels+=GetPixelChannels(image);
            }
            if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
              SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
            if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
              SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
            if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
              SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
            if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
                (image->colorspace == CMYKColorspace))
              SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
            if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
                (image->matte == MagickTrue))
              SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
          }
        else
          { /* Channel 'Sync' Flag, and Alpha Channel enabled.
            ** Weight the color channels with Alpha Channel so that
            ** transparent pixels are not part of the results.
            */
            MagickRealType
              alpha,  /* alpha weighting of colors : kernel*alpha  */
              gamma;  /* divisor, sum of color weighting values */

            gamma=0.0;
            for (v=0; v < (ssize_t) kernel->height; v++) {
              if ( IsNan(*k) ) continue;
              alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels));
              gamma += alpha;
              result.red     += alpha*GetPixelRed(image,k_pixels);
              result.green   += alpha*GetPixelGreen(image,k_pixels);
              result.blue    += alpha*GetPixelBlue(image,k_pixels);
              if (image->colorspace == CMYKColorspace)
                result.black += alpha*GetPixelBlack(image,k_pixels);
              result.alpha += (*k)*GetPixelAlpha(image,k_pixels);
              k--;
              k_pixels+=GetPixelChannels(image);
            }
            /* Sync'ed channels, all channels are modified */
            gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
            SetPixelRed(morphology_image,
              ClampToQuantum(gamma*result.red),q);
            SetPixelGreen(morphology_image,
              ClampToQuantum(gamma*result.green),q);
            SetPixelBlue(morphology_image,
              ClampToQuantum(gamma*result.blue),q);
            if (image->colorspace == CMYKColorspace)
              SetPixelBlack(morphology_image,
                ClampToQuantum(gamma*result.black),q);
            SetPixelAlpha(morphology_image,
              ClampToQuantum(result.alpha),q);
          }

        /* Count up changed pixels */
        if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q))
            || (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q))
            || (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q))
            || (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q))
            || ((image->colorspace == CMYKColorspace) &&
                (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
          changed++;  /* The pixel was changed in some way! */
        p+=GetPixelChannels(image);
        q+=GetPixelChannels(morphology_image);
      } /* y */
      if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
        status=MagickFalse;
      if (image->progress_monitor != (MagickProgressMonitor) NULL)
        {
          MagickBooleanType
            proceed;

#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp critical (MagickCore_MorphologyImage)
#endif
          proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
          if (proceed == MagickFalse)
            status=MagickFalse;
        }
    } /* x */
    morphology_image->type=image->type;
    morphology_view=DestroyCacheView(morphology_view);
    image_view=DestroyCacheView(image_view);
    return(status ? (ssize_t) changed : 0);
  }

  /*
  ** Normal handling of horizontal or rectangular kernels (row by row)
  */
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static,4) shared(progress,status)
#endif
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    register const Quantum
      *restrict p;

    register Quantum
      *restrict q;

    register ssize_t
      x;

    size_t
      r;

    if (status == MagickFalse)
      continue;
    p=GetCacheViewVirtualPixels(image_view, -offx, y-offy, virt_width,
      kernel->height,  exception);
    q=GetCacheViewAuthenticPixels(morphology_view,0,y,
      morphology_image->columns,1,exception);
    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
      {
        status=MagickFalse;
        continue;
      }
    /* offset to origin in 'p'. while 'q' points to it directly */
    r = virt_width*offy + offx;

    for (x=0; x < (ssize_t) image->columns; x++)
    {
       ssize_t
        v;

      register ssize_t
        u;

      register const MagickRealType
        *restrict k;

      register const Quantum
        *restrict k_pixels;

      PixelInfo
        result,
        min,
        max;

      /* Copy input image to the output image for unused channels
       * This removes need for 'cloning' a new image every iteration
       */
      SetPixelRed(morphology_image,GetPixelRed(image,p+r*
        GetPixelChannels(image)),q);
      SetPixelGreen(morphology_image,GetPixelGreen(image,p+r*
        GetPixelChannels(image)),q);
      SetPixelBlue(morphology_image,GetPixelBlue(image,p+r*
        GetPixelChannels(image)),q);
      if (image->colorspace == CMYKColorspace)
        SetPixelBlack(morphology_image,GetPixelBlack(image,p+r*
          GetPixelChannels(image)),q);

      /* Defaults */
      min.red     =
      min.green   =
      min.blue    =
      min.alpha =
      min.black   = (MagickRealType) QuantumRange;
      max.red     =
      max.green   =
      max.blue    =
      max.alpha =
      max.black   = (MagickRealType) 0;
      /* default result is the original pixel value */
      result.red     = (MagickRealType) GetPixelRed(image,p+r*GetPixelChannels(image));
      result.green   = (MagickRealType) GetPixelGreen(image,p+r*GetPixelChannels(image));
      result.blue    = (MagickRealType) GetPixelBlue(image,p+r*GetPixelChannels(image));
      result.black   = 0.0;
      if (image->colorspace == CMYKColorspace)
        result.black = (MagickRealType) GetPixelBlack(image,p+r*GetPixelChannels(image));
      result.alpha=(MagickRealType) GetPixelAlpha(image,p+r*GetPixelChannels(image));

      switch (method) {
        case ConvolveMorphology:
          /* Set the bias of the weighted average output */
          result.red     =
          result.green   =
          result.blue    =
          result.alpha =
          result.black   = bias;
          break;
        case DilateIntensityMorphology:
        case ErodeIntensityMorphology:
          /* use a boolean flag indicating when first match found */
          result.red = 0.0;  /* result is not used otherwise */
          break;
        default:
          break;
      }

      switch ( method ) {
        case ConvolveMorphology:
            /* Weighted Average of pixels using reflected kernel
            **
            ** NOTE for correct working of this operation for asymetrical
            ** kernels, the kernel needs to be applied in its reflected form.
            ** That is its values needs to be reversed.
            **
            ** Correlation is actually the same as this but without reflecting
            ** the kernel, and thus 'lower-level' that Convolution.  However
            ** as Convolution is the more common method used, and it does not
            ** really cost us much in terms of processing to use a reflected
            ** kernel, so it is Convolution that is implemented.
            **
            ** Correlation will have its kernel reflected before calling
            ** this function to do a Convolve.
            **
            ** For more details of Correlation vs Convolution see
            **   http://www.cs.umd.edu/~djacobs/CMSC426/Convolution.pdf
            */
            k = &kernel->values[ kernel->width*kernel->height-1 ];
            k_pixels = p;
            if ( (image->channel_mask != DefaultChannels) ||
                 (image->matte == MagickFalse) )
              { /* No 'Sync' involved.
                ** Convolution is simple greyscale channel operation
                */
                for (v=0; v < (ssize_t) kernel->height; v++) {
                  for (u=0; u < (ssize_t) kernel->width; u++, k--) {
                    if ( IsNan(*k) ) continue;
                    result.red     += (*k)*
                      GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
                    result.green   += (*k)*
                      GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
                    result.blue    += (*k)*
                      GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
                    if (image->colorspace == CMYKColorspace)
                      result.black += (*k)*
                        GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
                    result.alpha += (*k)*
                      GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
                  }
                  k_pixels += virt_width*GetPixelChannels(image);
                }
                if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
                  SetPixelRed(morphology_image,ClampToQuantum(result.red),
                    q);
                if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
                  SetPixelGreen(morphology_image,ClampToQuantum(result.green),
                    q);
                if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
                  SetPixelBlue(morphology_image,ClampToQuantum(result.blue),
                    q);
                if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
                    (image->colorspace == CMYKColorspace))
                  SetPixelBlack(morphology_image,ClampToQuantum(result.black),
                    q);
                if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
                    (image->matte == MagickTrue))
                  SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),
                    q);
              }
            else
              { /* Channel 'Sync' Flag, and Alpha Channel enabled.
                ** Weight the color channels with Alpha Channel so that
                ** transparent pixels are not part of the results.
                */
                MagickRealType
                  alpha,  /* alpha weighting of colors : kernel*alpha  */
                  gamma;  /* divisor, sum of color weighting values */

                gamma=0.0;
                for (v=0; v < (ssize_t) kernel->height; v++) {
                  for (u=0; u < (ssize_t) kernel->width; u++, k--) {
                    if ( IsNan(*k) ) continue;
                    alpha=(*k)*(QuantumScale*GetPixelAlpha(image,k_pixels+u*
                      GetPixelChannels(image)));
                    gamma += alpha;
                    result.red     += alpha*
                      GetPixelRed(image,k_pixels+u*GetPixelChannels(image));
                    result.green   += alpha*
                      GetPixelGreen(image,k_pixels+u*GetPixelChannels(image));
                    result.blue    += alpha*
                      GetPixelBlue(image,k_pixels+u*GetPixelChannels(image));
                    if (image->colorspace == CMYKColorspace)
                      result.black+=alpha*
                        GetPixelBlack(image,k_pixels+u*GetPixelChannels(image));
                    result.alpha += (*k)*
                      GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image));
                  }
                  k_pixels += virt_width*GetPixelChannels(image);
                }
                /* Sync'ed channels, all channels are modified */
                gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
                SetPixelRed(morphology_image,
                  ClampToQuantum(gamma*result.red),q);
                SetPixelGreen(morphology_image,
                  ClampToQuantum(gamma*result.green),q);
                SetPixelBlue(morphology_image,
                  ClampToQuantum(gamma*result.blue),q);
                if (image->colorspace == CMYKColorspace)
                  SetPixelBlack(morphology_image,
                    ClampToQuantum(gamma*result.black),q);
                SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
              }
            break;

        case ErodeMorphology:
            /* Minimum Value within kernel neighbourhood
            **
            ** NOTE that the kernel is not reflected for this operation!
            **
            ** NOTE: in normal Greyscale Morphology, the kernel value should
            ** be added to the real value, this is currently not done, due to
            ** the nature of the boolean kernels being used.
            */
            k = kernel->values;
            k_pixels = p;
            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)
                  GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(min.green,   (double) 
                  GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(min.blue,    (double)
                  GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
                if (image->colorspace == CMYKColorspace)
                  Minimize(min.black,(double)
                    GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(min.alpha,(double)
                  GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
              }
              k_pixels += virt_width*GetPixelChannels(image);
            }
            break;

        case DilateMorphology:
            /* Maximum Value within kernel neighbourhood
            **
            ** NOTE for correct working of this operation for asymetrical
            ** kernels, the kernel needs to be applied in its reflected form.
            ** That is its values needs to be reversed.
            **
            ** NOTE: in normal Greyscale Morphology, the kernel value should
            ** be added to the real value, this is currently not done, due to
            ** the nature of the boolean kernels being used.
            **
            */
            k = &kernel->values[ kernel->width*kernel->height-1 ];
            k_pixels = p;
            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)
                  GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
                Maximize(max.green,   (double) 
                  GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
                Maximize(max.blue,    (double) 
                  GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
                if (image->colorspace == CMYKColorspace)
                  Maximize(max.black,   (double)
                    GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
                Maximize(max.alpha,(double) 
                  GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
              }
              k_pixels += virt_width*GetPixelChannels(image);
            }
            break;

        case HitAndMissMorphology:
        case ThinningMorphology:
        case ThickenMorphology:
            /* Minimum of Foreground Pixel minus Maxumum of Background Pixels
            **
            ** NOTE that the kernel is not reflected for this operation,
            ** and consists of both foreground and background pixel
            ** neighbourhoods, 0.0 for background, and 1.0 for foreground
            ** with either Nan or 0.5 values for don't care.
            **
            ** 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;
            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 */
                  Minimize(min.red,     (double)
                    GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
                  Minimize(min.green,   (double)
                    GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
                  Minimize(min.blue,    (double)
                    GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
                  if ( image->colorspace == CMYKColorspace)
                    Minimize(min.black,(double)
                      GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
                  Minimize(min.alpha,(double)
                    GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
                }
                else if ( (*k) < 0.3 )
                { /* maximum of background pixels */
                  Maximize(max.red,     (double)
                    GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
                  Maximize(max.green,   (double)
                    GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
                  Maximize(max.blue,    (double)
                    GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
                  if (image->colorspace == CMYKColorspace)
                    Maximize(max.black,   (double)
                      GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
                  Maximize(max.alpha,(double)
                    GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
                }
              }
              k_pixels += virt_width*GetPixelChannels(image);
            }
            /* 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 );
            min.black   -= max.black;   Maximize( min.black,   0.0 );
            min.alpha -= max.alpha; Maximize( min.alpha, 0.0 );
            break;

        case ErodeIntensityMorphology:
            /* Select Pixel with Minimum Intensity within kernel neighbourhood
            **
            ** WARNING: the intensity test fails for CMYK and does not
            ** take into account the moderating effect of the alpha channel
            ** on the intensity.
            **
            ** NOTE that the kernel is not reflected for this operation!
            */
            k = kernel->values;
            k_pixels = p;
            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 ||
                     GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) < GetPixelIntensity(morphology_image,q) ) {
                  /* copy the whole pixel - no channel selection */
                  SetPixelRed(morphology_image,GetPixelRed(image,
                    k_pixels+u*GetPixelChannels(image)),q);
                  SetPixelGreen(morphology_image,GetPixelGreen(image,
                    k_pixels+u*GetPixelChannels(image)),q);
                  SetPixelBlue(morphology_image,GetPixelBlue(image,
                    k_pixels+u*GetPixelChannels(image)),q);
                  SetPixelAlpha(morphology_image,GetPixelAlpha(image,
                    k_pixels+u*GetPixelChannels(image)),q);
                  if ( result.red > 0.0 ) changed++;
                  result.red = 1.0;
                }
              }
              k_pixels += virt_width*GetPixelChannels(image);
            }
            break;

        case DilateIntensityMorphology:
            /* Select Pixel with Maximum Intensity within kernel neighbourhood
            **
            ** WARNING: the intensity test fails for CMYK and does not
            ** take into account the moderating effect of the alpha channel
            ** on the intensity (yet).
            **
            ** NOTE for correct working of this operation for asymetrical
            ** kernels, the kernel needs to be applied in its reflected form.
            ** That is its values needs to be reversed.
            */
            k = &kernel->values[ kernel->width*kernel->height-1 ];
            k_pixels = p;
            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 ||
                     GetPixelIntensity(image,k_pixels+u*GetPixelChannels(image)) > GetPixelIntensity(morphology_image,q) ) {
                  /* copy the whole pixel - no channel selection */
                  SetPixelRed(morphology_image,GetPixelRed(image,
                    k_pixels+u*GetPixelChannels(image)),q);
                  SetPixelGreen(morphology_image,GetPixelGreen(image,
                    k_pixels+u*GetPixelChannels(image)),q);
                  SetPixelBlue(morphology_image,GetPixelBlue(image,
                    k_pixels+u*GetPixelChannels(image)),q);
                  SetPixelAlpha(morphology_image,GetPixelAlpha(image,
                    k_pixels+u*GetPixelChannels(image)),q);
                  if ( result.red > 0.0 ) changed++;
                  result.red = 1.0;
                }
              }
              k_pixels += virt_width*GetPixelChannels(image);
            }
            break;
#if 0
  This code has been obsoleted by the MorphologyPrimitiveDirect() function.
  However it is still (almost) correct coding for Grayscale Morphology.
  That is...

  GrayErode    is equivalent but with kernel values subtracted from pixels
               without the kernel rotation
  GreyDilate   is equivalent but using Maximum() instead of Minimum()
               using kernel rotation

  It has thus been preserved for future implementation of those methods.

        case DistanceMorphology:
            /* Add kernel Value and select the minimum value found.
            ** The result is a iterative distance from edge of image shape.
            **
            ** All Distance Kernels are symetrical, but that may not always
            ** be the case. For example how about a distance from left edges?
            ** To work correctly with asymetrical kernels the reflected kernel
            ** needs to be applied.
            */
            k = &kernel->values[ kernel->width*kernel->height-1 ];
            k_pixels = p;
            for (v=0; v < (ssize_t) kernel->height; v++) {
              for (u=0; u < (ssize_t) kernel->width; u++, k--) {
                if ( IsNan(*k) ) continue;
                Minimize(result.red,     (*k)+k_pixels[u].red);
                Minimize(result.green,   (*k)+k_pixels[u].green);
                Minimize(result.blue,    (*k)+k_pixels[u].blue);
                Minimize(result.alpha, (*k)+k_pixels[u].alpha);
                if ( image->colorspace == CMYKColorspace)
                  Minimize(result.black,(*k)+GetPixelBlack(p_image,k_indexes+u));
              }
              k_pixels += virt_width*GetPixelChannels(image);
            }
            break;
#endif
        case UndefinedMorphology:
        default:
            break; /* Do nothing */
      }
      /* Final mathematics of results (combine with original image?)
      **
      ** NOTE: Difference Morphology operators Edge* and *Hat could also
      ** be done here but works better with iteration as a image difference
      ** in the controling function (below).  Thicken and Thinning however
      ** should be done here so thay can be iterated correctly.
      */
      switch ( method ) {
        case HitAndMissMorphology:
        case ErodeMorphology:
          result = min;    /* minimum of neighbourhood */
          break;
        case DilateMorphology:
          result = max;    /* maximum of neighbourhood */
          break;
        case ThinningMorphology:
          /* subtract pattern match from original */
          result.red     -= min.red;
          result.green   -= min.green;
          result.blue    -= min.blue;
          result.black   -= min.black;
          result.alpha -= min.alpha;
          break;
        case ThickenMorphology:
          /* Add the pattern matchs to the original */
          result.red     += min.red;
          result.green   += min.green;
          result.blue    += min.blue;
          result.black   += min.black;
          result.alpha += min.alpha;
          break;
        default:
          /* result directly calculated or assigned */
          break;
      }
      /* Assign the resulting pixel values - Clamping Result */
      switch ( method ) {
        case UndefinedMorphology:
        case ConvolveMorphology:
        case DilateIntensityMorphology:
        case ErodeIntensityMorphology:
          break;  /* full pixel was directly assigned - not a channel method */
        default:
          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
            SetPixelRed(morphology_image,ClampToQuantum(result.red),q);
          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
            SetPixelGreen(morphology_image,ClampToQuantum(result.green),q);
          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
            SetPixelBlue(morphology_image,ClampToQuantum(result.blue),q);
          if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
              (image->colorspace == CMYKColorspace))
            SetPixelBlack(morphology_image,ClampToQuantum(result.black),q);
          if (((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0) &&
              (image->matte == MagickTrue))
            SetPixelAlpha(morphology_image,ClampToQuantum(result.alpha),q);
          break;
      }
      /* Count up changed pixels */
      if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(morphology_image,q)) ||
          (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(morphology_image,q)) ||
          (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(morphology_image,q)) ||
          (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(morphology_image,q)) ||
          ((image->colorspace == CMYKColorspace) &&
           (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(morphology_image,q))))
        changed++;  /* The pixel was changed in some way! */
      p+=GetPixelChannels(image);
      q+=GetPixelChannels(morphology_image);
    } /* x */
    if ( SyncCacheViewAuthenticPixels(morphology_view,exception) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp critical (MagickCore_MorphologyImage)
#endif
        proceed=SetImageProgress(image,MorphologyTag,progress++,image->rows);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
  } /* y */
  morphology_view=DestroyCacheView(morphology_view);
  image_view=DestroyCacheView(image_view);
  return(status ? (ssize_t)changed : -1);
}

/* This is almost identical to the MorphologyPrimative() function above,
** but will apply the primitive directly to the image in two passes.
**
** That is after each row is 'Sync'ed' into the image, the next row will
** make use of those values as part of the calculation of the next row.
** It then repeats, but going in the oppisite (bottom-up) direction.
**
** Because of this 'iterative' handling this function can not make use
** of multi-threaded, parellel processing.
*/
static ssize_t MorphologyPrimitiveDirect(Image *image,
  const MorphologyMethod method,const KernelInfo *kernel,
  ExceptionInfo *exception)
{
  CacheView
    *auth_view,
    *virt_view;

  MagickBooleanType
    status;

  MagickOffsetType
    progress;

  ssize_t
    y, offx, offy;

  size_t
    virt_width,
    changed;

  status=MagickTrue;
  changed=0;
  progress=0;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  assert(kernel != (KernelInfo *) NULL);
  assert(kernel->signature == MagickSignature);
  assert(exception != (ExceptionInfo *) NULL);
  assert(exception->signature == MagickSignature);

  /* Some methods (including convolve) needs use a reflected kernel.
   * Adjust 'origin' offsets to loop though kernel as a reflection.
   */
  offx = kernel->x;
  offy = kernel->y;
  switch(method) {
    case DistanceMorphology:
    case VoronoiMorphology:
      /* kernel needs to used with reflection about origin */
      offx = (ssize_t) kernel->width-offx-1;
      offy = (ssize_t) kernel->height-offy-1;
      break;
#if 0
    case ?????Morphology:
      /* kernel is used as is, without reflection */
      break;
#endif
    default:
      assert("Not a PrimativeDirect Morphology Method" != (char *) NULL);
      break;
  }

  /* DO NOT THREAD THIS CODE! */
  /* two views into same image (virtual, and actual) */
  virt_view=AcquireCacheView(image);
  auth_view=AcquireCacheView(image);
  virt_width=image->columns+kernel->width-1;

  for (y=0; y < (ssize_t) image->rows; y++)
  {
    register const Quantum
      *restrict p;

    register Quantum
      *restrict q;

    register ssize_t
      x;

    ssize_t
      r;

    /* NOTE read virtual pixels, and authentic pixels, from the same image!
    ** we read using virtual to get virtual pixel handling, but write back
    ** into the same image.
    **
    ** Only top half of kernel is processed as we do a single pass downward
    ** through the image iterating the distance function as we go.
    */
    if (status == MagickFalse)
      break;
    p=GetCacheViewVirtualPixels(virt_view,-offx,y-offy,virt_width,(size_t)
      offy+1,exception);
    q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
      exception);
    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
      status=MagickFalse;
    if (status == MagickFalse)
      break;

    /* offset to origin in 'p'. while 'q' points to it directly */
    r = (ssize_t) virt_width*offy + offx;

    for (x=0; x < (ssize_t) image->columns; x++)
    {
      ssize_t
        v;

      register ssize_t
        u;

      register const MagickRealType
        *restrict k;

      register const Quantum
        *restrict k_pixels;

      PixelInfo
        result;

      /* Starting Defaults */
      GetPixelInfo(image,&result);
      GetPixelInfoPixel(image,q,&result);
      if ( method != VoronoiMorphology )
        result.alpha = QuantumRange - result.alpha;

      switch ( method ) {
        case DistanceMorphology:
            /* Add kernel Value and select the minimum value found. */
            k = &kernel->values[ kernel->width*kernel->height-1 ];
            k_pixels = p;
            for (v=0; v <= (ssize_t) offy; v++) {
              for (u=0; u < (ssize_t) kernel->width; u++, k--) {
                if ( IsNan(*k) ) continue;
                Minimize(result.red,     (*k)+
                  GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.green,   (*k)+
                  GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.blue,    (*k)+
                  GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
                if (image->colorspace == CMYKColorspace)
                  Minimize(result.black,(*k)+
                    GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.alpha, (*k)+
                  GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
              }
              k_pixels += virt_width*GetPixelChannels(image);
            }
            /* repeat with the just processed pixels of this row */
            k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
            k_pixels = q-offx*GetPixelChannels(image);
              for (u=0; u < (ssize_t) offx; u++, k--) {
                if ( x+u-offx < 0 ) continue;  /* off the edge! */
                if ( IsNan(*k) ) continue;
                Minimize(result.red,     (*k)+
                  GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.green,   (*k)+
                  GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.blue,    (*k)+
                  GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
                if (image->colorspace == CMYKColorspace)
                  Minimize(result.black,(*k)+
                    GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.alpha,(*k)+
                  GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
              }
            break;
        case VoronoiMorphology:
            /* Apply Distance to 'Matte' channel, coping the closest color.
            **
            ** This is experimental, and realy the 'alpha' component should
            ** be completely separate 'masking' channel.
            */
            k = &kernel->values[ kernel->width*kernel->height-1 ];
            k_pixels = p;
            for (v=0; v <= (ssize_t) offy; v++) {
              for (u=0; u < (ssize_t) kernel->width; u++, k--) {
                if ( IsNan(*k) ) continue;
                if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
                  {
                    GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
                      &result);
                    result.alpha += *k;
                  }
              }
              k_pixels += virt_width*GetPixelChannels(image);
            }
            /* repeat with the just processed pixels of this row */
            k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
            k_pixels = q-offx*GetPixelChannels(image);
              for (u=0; u < (ssize_t) offx; u++, k--) {
                if ( x+u-offx < 0 ) continue;  /* off the edge! */
                if ( IsNan(*k) ) continue;
                if( result.alpha > (*k)+GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)) )
                  {
                    GetPixelInfoPixel(image,k_pixels+u*GetPixelChannels(image),
                      &result);
                    result.alpha += *k;
                  }
              }
            break;
        default:
          /* result directly calculated or assigned */
          break;
      }
      /* Assign the resulting pixel values - Clamping Result */
      switch ( method ) {
        case VoronoiMorphology:
          SetPixelInfoPixel(image,&result,q);
          break;
        default:
          if ((GetPixelRedTraits(image) & UpdatePixelTrait) != 0)
            SetPixelRed(image,ClampToQuantum(result.red),q);
          if ((GetPixelGreenTraits(image) & UpdatePixelTrait) != 0)
            SetPixelGreen(image,ClampToQuantum(result.green),q);
          if ((GetPixelBlueTraits(image) & UpdatePixelTrait) != 0)
            SetPixelBlue(image,ClampToQuantum(result.blue),q);
          if (((GetPixelBlackTraits(image) & UpdatePixelTrait) != 0) &&
              (image->colorspace == CMYKColorspace))
            SetPixelBlack(image,ClampToQuantum(result.black),q);
          if ((GetPixelAlphaTraits(image) & UpdatePixelTrait) != 0 &&
              (image->matte == MagickTrue))
            SetPixelAlpha(image,ClampToQuantum(result.alpha),q);
          break;
      }
      /* Count up changed pixels */
      if ((GetPixelRed(image,p+r*GetPixelChannels(image)) != GetPixelRed(image,q)) ||
          (GetPixelGreen(image,p+r*GetPixelChannels(image)) != GetPixelGreen(image,q)) ||
          (GetPixelBlue(image,p+r*GetPixelChannels(image)) != GetPixelBlue(image,q)) ||
          (GetPixelAlpha(image,p+r*GetPixelChannels(image)) != GetPixelAlpha(image,q)) ||
          ((image->colorspace == CMYKColorspace) &&
           (GetPixelBlack(image,p+r*GetPixelChannels(image)) != GetPixelBlack(image,q))))
        changed++;  /* The pixel was changed in some way! */

      p+=GetPixelChannels(image); /* increment pixel buffers */
      q+=GetPixelChannels(image);
    } /* x */

    if ( SyncCacheViewAuthenticPixels(auth_view,exception) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      if ( SetImageProgress(image,MorphologyTag,progress++,image->rows)
                == MagickFalse )
        status=MagickFalse;

  } /* y */

  /* Do the reversed pass through the image */
  for (y=(ssize_t)image->rows-1; y >= 0; y--)
  {
    register const Quantum
      *restrict p;

    register Quantum
      *restrict q;

    register ssize_t
      x;

    ssize_t
      r;

    if (status == MagickFalse)
      break;
    /* NOTE read virtual pixels, and authentic pixels, from the same image!
    ** we read using virtual to get virtual pixel handling, but write back
    ** into the same image.
    **
    ** Only the bottom half of the kernel will be processes as we
    ** up the image.
    */
    p=GetCacheViewVirtualPixels(virt_view,-offx,y,virt_width,(size_t)
      kernel->y+1,exception);
    q=GetCacheViewAuthenticPixels(auth_view, 0, y, image->columns, 1,
      exception);
    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
      status=MagickFalse;
    if (status == MagickFalse)
      break;

    /* adjust positions to end of row */
    p += (image->columns-1)*GetPixelChannels(image);
    q += (image->columns-1)*GetPixelChannels(image);

    /* offset to origin in 'p'. while 'q' points to it directly */
    r = offx;

    for (x=(ssize_t)image->columns-1; x >= 0; x--)
    {
      ssize_t
        v;

      register ssize_t
        u;

      register const MagickRealType
        *restrict k;

      register const Quantum
        *restrict k_pixels;

      PixelInfo
        result;

      /* Default - previously modified pixel */
      GetPixelInfo(image,&result);
      GetPixelInfoPixel(image,q,&result);
      if ( method != VoronoiMorphology )
        result.alpha = QuantumRange - result.alpha;

      switch ( method ) {
        case DistanceMorphology:
            /* Add kernel Value and select the minimum value found. */
            k = &kernel->values[ kernel->width*(kernel->y+1)-1 ];
            k_pixels = p;
            for (v=offy; v < (ssize_t) kernel->height; v++) {
              for (u=0; u < (ssize_t) kernel->width; u++, k--) {
                if ( IsNan(*k) ) continue;
                Minimize(result.red,     (*k)+
                  GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.green,   (*k)+
                  GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.blue,    (*k)+
                  GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
                if ( image->colorspace == CMYKColorspace)
                  Minimize(result.black,(*k)+
                    GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.alpha, (*k)+
                  GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
              }
              k_pixels += virt_width*GetPixelChannels(image);
            }
            /* repeat with the just processed pixels of this row */
            k = &kernel->values[ kernel->width*(kernel->y)+kernel->x-1 ];
            k_pixels = q-offx*GetPixelChannels(image);
              for (u=offx+1; u < (ssize_t) kernel->width; u++, k--) {
                if ( (x+u-offx) >= (ssize_t)image->columns ) continue;
                if ( IsNan(*k) ) continue;
                Minimize(result.red,     (*k)+
                  GetPixelRed(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.green,   (*k)+
                  GetPixelGreen(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.blue,    (*k)+
                  GetPixelBlue(image,k_pixels+u*GetPixelChannels(image)));
                if ( image->colorspace == CMYKColorspace)
                  Minimize(result.black,   (*k)+
                    GetPixelBlack(image,k_pixels+u*GetPixelChannels(image)));
                Minimize(result.alpha, (*k)+
                  GetPixelAlpha(image,k_pixels+u*GetPixelChannels(image)));
              }
            break;
        case