]> granicus.if.org Git - imagemagick/blobdiff - MagickCore/feature.c
Update web pages
[imagemagick] / MagickCore / feature.c
index 8ae6005128298b1530833acb6d461a3db8ee2364..4f5b23700ff17c2576f867290122dd673030220e 100644 (file)
@@ -17,7 +17,7 @@
 %                                 July 1992                                   %
 %                                                                             %
 %                                                                             %
-%  Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization      %
+%  Copyright 1999-2015 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  %
@@ -48,6 +48,7 @@
 #include "MagickCore/cache.h"
 #include "MagickCore/cache-private.h"
 #include "MagickCore/cache-view.h"
+#include "MagickCore/channel.h"
 #include "MagickCore/client.h"
 #include "MagickCore/color.h"
 #include "MagickCore/color-private.h"
 %  CannyEdgeImage() uses a multi-stage algorithm to detect a wide range of
 %  edges in images.
 %
-%  The format of the EdgeImage method is:
+%  The format of the CannyEdgeImage method is:
 %
 %      Image *CannyEdgeImage(const Image *image,const double radius,
 %        const double sigma,const double lower_percent,
@@ -237,6 +238,8 @@ MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
   const double sigma,const double lower_percent,const double upper_percent,
   ExceptionInfo *exception)
 {
+#define CannyEdgeImageTag  "CannyEdge/Image"
+
   CacheView
     *edge_view;
 
@@ -244,7 +247,7 @@ MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
     pixel;
 
   char
-    geometry[MaxTextExtent];
+    geometry[MagickPathExtent];
 
   double
     lower_threshold,
@@ -261,6 +264,9 @@ MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
   MagickBooleanType
     status;
 
+  MagickOffsetType
+    progress;
+
   MatrixInfo
     *canny_cache;
 
@@ -268,17 +274,17 @@ MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
     y;
 
   assert(image != (const Image *) NULL);
-  assert(image->signature == MagickSignature);
+  assert(image->signature == MagickCoreSignature);
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
   assert(exception != (ExceptionInfo *) NULL);
-  assert(exception->signature == MagickSignature);
+  assert(exception->signature == MagickCoreSignature);
   /*
     Filter out noise.
   */
-  (void) FormatLocaleString(geometry,MaxTextExtent,
+  (void) FormatLocaleString(geometry,MagickPathExtent,
     "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
-  kernel_info=AcquireKernelInfo(geometry);
+  kernel_info=AcquireKernelInfo(geometry,exception);
   if (kernel_info == (KernelInfo *) NULL)
     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
   edge_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info,
@@ -291,6 +297,7 @@ MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
       edge_image=DestroyImage(edge_image);
       return((Image *) NULL);
     }
+  (void) SetImageAlphaChannel(edge_image,OffAlphaChannel,exception);
   /*
     Find the intensity gradient of the image.
   */
@@ -410,6 +417,7 @@ MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
     Non-maxima suppression, remove pixels that are not considered to be part
     of an edge.
   */
+  progress=0;
   (void) GetMatrixElement(canny_cache,0,0,&pixel);
   max=pixel.intensity;
   min=pixel.intensity;
@@ -542,6 +550,19 @@ MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
         status=TraceEdges(edge_image,edge_view,canny_cache,x,y,lower_threshold,
           exception);
     }
+    if (image->progress_monitor != (MagickProgressMonitor) NULL)
+      {
+        MagickBooleanType
+          proceed;
+
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+        #pragma omp critical (MagickCore_CannyEdgeImage)
+#endif
+        proceed=SetImageProgress(image,CannyEdgeImageTag,progress++,
+          image->rows);
+        if (proceed == MagickFalse)
+          status=MagickFalse;
+      }
   }
   edge_view=DestroyCacheView(edge_view);
   /*
@@ -556,112 +577,139 @@ MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
 %                                                                             %
 %                                                                             %
 %                                                                             %
-%     H o u g h L i n e s I m a g e                                           %
+%   G e t I m a g e F e a t u r e s                                           %
 %                                                                             %
 %                                                                             %
 %                                                                             %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %
-%  HoughLinesImage() identifies lines in the image.
+%  GetImageFeatures() returns features for each channel in the image in
+%  each of four directions (horizontal, vertical, left and right diagonals)
+%  for the specified distance.  The features include the angular second
+%  moment, contrast, correlation, sum of squares: variance, inverse difference
+%  moment, sum average, sum varience, sum entropy, entropy, difference variance,%  difference entropy, information measures of correlation 1, information
+%  measures of correlation 2, and maximum correlation coefficient.  You can
+%  access the red channel contrast, for example, like this:
+%
+%      channel_features=GetImageFeatures(image,1,exception);
+%      contrast=channel_features[RedPixelChannel].contrast[0];
 %
-%  The format of the HoughLinesImage method is:
+%  Use MagickRelinquishMemory() to free the features buffer.
 %
-%      Image *HoughLinesImage(const Image *image,const size_t width,
-%        const size_t height,const size_t threshold,ExceptionInfo *exception)
+%  The format of the GetImageFeatures method is:
+%
+%      ChannelFeatures *GetImageFeatures(const Image *image,
+%        const size_t distance,ExceptionInfo *exception)
 %
 %  A description of each parameter follows:
 %
 %    o image: the image.
 %
-%    o width, height: find line pairs as local maxima in this neighborhood.
-%
-%    o threshold: the line count threshold.
+%    o distance: the distance.
 %
 %    o exception: return any errors or warnings in this structure.
 %
 */
 
-static inline double MagickRound(double x)
+static inline double MagickLog10(const double x)
 {
-  /*
-    Round the fraction to nearest integer.
-  */
-  if ((x-floor(x)) < (ceil(x)-x))
-    return(floor(x));
-  return(ceil(x));
+#define Log10Epsilon  (1.0e-11)
+
+ if (fabs(x) < Log10Epsilon)
+   return(log10(Log10Epsilon));
+ return(log10(fabs(x)));
 }
 
-MagickExport Image *HoughLinesImage(const Image *image,const size_t width,
-  const size_t height,const size_t threshold,ExceptionInfo *exception)
+MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
+  const size_t distance,ExceptionInfo *exception)
 {
+  typedef struct _ChannelStatistics
+  {
+    PixelInfo
+      direction[4];  /* horizontal, vertical, left and right diagonals */
+  } ChannelStatistics;
+
   CacheView
     *image_view;
 
-  char
-    message[MaxTextExtent],
-    path[MaxTextExtent];
-
-  const char
-    *artifact;
-
-  double
-    hough_height;
-
-  Image
-    *lines_image = NULL;
+  ChannelFeatures
+    *channel_features;
 
-  ImageInfo
-    *image_info;
+  ChannelStatistics
+    **cooccurrence,
+    correlation,
+    *density_x,
+    *density_xy,
+    *density_y,
+    entropy_x,
+    entropy_xy,
+    entropy_xy1,
+    entropy_xy2,
+    entropy_y,
+    mean,
+    **Q,
+    *sum,
+    sum_squares,
+    variance;
 
-  int
-    file;
+  PixelPacket
+    gray,
+    *grays;
 
   MagickBooleanType
     status;
 
-  MatrixInfo
-    *accumulator;
+  register ssize_t
+    i;
 
-  PointInfo
-    center;
+  size_t
+    length;
 
-  register ssize_t
+  ssize_t
     y;
 
-  size_t
-    accumulator_height,
-    accumulator_width,
-    line_count;
+  unsigned int
+    number_grays;
 
-  /*
-    Create the accumulator.
-  */
-  assert(image != (const Image *) NULL);
-  assert(image->signature == MagickSignature);
+  assert(image != (Image *) NULL);
+  assert(image->signature == MagickCoreSignature);
   if (image->debug != MagickFalse)
     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  assert(exception != (ExceptionInfo *) NULL);
-  assert(exception->signature == MagickSignature);
-  accumulator_width=180;
-  hough_height=((sqrt(2.0)*(double) (image->rows > image->columns ?
-    image->rows : image->columns))/2.0);
-  accumulator_height=(size_t) (2.0*hough_height);
-  accumulator=AcquireMatrixInfo(accumulator_width,accumulator_height,
-    sizeof(double),exception);
-  if (accumulator == (MatrixInfo *) NULL)
-    ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
-  if (NullMatrix(accumulator) == MagickFalse)
-    {
-      accumulator=DestroyMatrixInfo(accumulator);
-      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
-    }
+  if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
+    return((ChannelFeatures *) NULL);
+  length=CompositeChannels+1UL;
+  channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
+    sizeof(*channel_features));
+  if (channel_features == (ChannelFeatures *) NULL)
+    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
+  (void) ResetMagickMemory(channel_features,0,length*
+    sizeof(*channel_features));
   /*
-    Populate the accumulator.
+    Form grays.
   */
+  grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays));
+  if (grays == (PixelPacket *) NULL)
+    {
+      channel_features=(ChannelFeatures *) RelinquishMagickMemory(
+        channel_features);
+      (void) ThrowMagickException(exception,GetMagickModule(),
+        ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
+      return(channel_features);
+    }
+  for (i=0; i <= (ssize_t) MaxMap; i++)
+  {
+    grays[i].red=(~0U);
+    grays[i].green=(~0U);
+    grays[i].blue=(~0U);
+    grays[i].alpha=(~0U);
+    grays[i].black=(~0U);
+  }
   status=MagickTrue;
-  center.x=(double) image->columns/2.0;
-  center.y=(double) image->rows/2.0;
   image_view=AcquireVirtualCacheView(image,exception);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(static,4) shared(status) \
+    magick_threads(image,image,image->rows,1)
+#endif
   for (y=0; y < (ssize_t) image->rows; y++)
   {
     register const Quantum
@@ -673,511 +721,179 @@ MagickExport Image *HoughLinesImage(const Image *image,const size_t width,
     if (status == MagickFalse)
       continue;
     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
-    if (p == (Quantum *) NULL)
+    if (p == (const Quantum *) NULL)
       {
         status=MagickFalse;
         continue;
       }
     for (x=0; x < (ssize_t) image->columns; x++)
     {
-      if (GetPixelIntensity(image,p) > (QuantumRange/2))
-        {
-          register ssize_t
-            i;
-
-          for (i=0; i < 180; i++)
-          {
-            double
-              count,
-              radius;
-
-            radius=(((double) x-center.x)*cos(DegreesToRadians((double) i)))+
-              (((double) y-center.y)*sin(DegreesToRadians((double) i)));
-            (void) GetMatrixElement(accumulator,i,(ssize_t)
-              MagickRound(radius+hough_height),&count);
-            count++;
-            (void) SetMatrixElement(accumulator,i,(ssize_t)
-              MagickRound(radius+hough_height),&count);
-          }
-        }
+      grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
+        ScaleQuantumToMap(GetPixelRed(image,p));
+      grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
+        ScaleQuantumToMap(GetPixelGreen(image,p));
+      grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
+        ScaleQuantumToMap(GetPixelBlue(image,p));
+      if (image->colorspace == CMYKColorspace)
+        grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
+          ScaleQuantumToMap(GetPixelBlack(image,p));
+      if (image->alpha_trait != UndefinedPixelTrait)
+        grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
+          ScaleQuantumToMap(GetPixelAlpha(image,p));
       p+=GetPixelChannels(image);
     }
   }
   image_view=DestroyCacheView(image_view);
   if (status == MagickFalse)
     {
-      accumulator=DestroyMatrixInfo(accumulator);
-      return((Image *) NULL);
+      grays=(PixelPacket *) RelinquishMagickMemory(grays);
+      channel_features=(ChannelFeatures *) RelinquishMagickMemory(
+        channel_features);
+      return(channel_features);
     }
+  (void) ResetMagickMemory(&gray,0,sizeof(gray));
+  for (i=0; i <= (ssize_t) MaxMap; i++)
+  {
+    if (grays[i].red != ~0U)
+      grays[gray.red++].red=grays[i].red;
+    if (grays[i].green != ~0U)
+      grays[gray.green++].green=grays[i].green;
+    if (grays[i].blue != ~0U)
+      grays[gray.blue++].blue=grays[i].blue;
+    if (image->colorspace == CMYKColorspace)
+      if (grays[i].black != ~0U)
+        grays[gray.black++].black=grays[i].black;
+    if (image->alpha_trait != UndefinedPixelTrait)
+      if (grays[i].alpha != ~0U)
+        grays[gray.alpha++].alpha=grays[i].alpha;
+  }
   /*
-    Generate line segments from accumulator.
+    Allocate spatial dependence matrix.
   */
-  file=AcquireUniqueFileResource(path);
-  if (file == -1)
-    {
-      accumulator=DestroyMatrixInfo(accumulator);
-      return((Image *) NULL);
-    }
-  (void) FormatLocaleString(message,MaxTextExtent,"viewbox 0 0 %.20g %.20g\n",
-    (double) image->columns,(double) image->rows);
-  (void) write(file,message,strlen(message));
-  line_count=image->columns > image->rows ? image->columns/4 : image->rows/4;
-  if (threshold != 0)
-    line_count=threshold;
-  for (y=0; y < (ssize_t) accumulator_height; y++)
-  {
-    register ssize_t
-      x;
-
-    for (x=0; x < (ssize_t) accumulator_width; x++)
+  number_grays=gray.red;
+  if (gray.green > number_grays)
+    number_grays=gray.green;
+  if (gray.blue > number_grays)
+    number_grays=gray.blue;
+  if (image->colorspace == CMYKColorspace)
+    if (gray.black > number_grays)
+      number_grays=gray.black;
+  if (image->alpha_trait != UndefinedPixelTrait)
+    if (gray.alpha > number_grays)
+      number_grays=gray.alpha;
+  cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays,
+    sizeof(*cooccurrence));
+  density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
+    sizeof(*density_x));
+  density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
+    sizeof(*density_xy));
+  density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
+    sizeof(*density_y));
+  Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q));
+  sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum));
+  if ((cooccurrence == (ChannelStatistics **) NULL) ||
+      (density_x == (ChannelStatistics *) NULL) ||
+      (density_xy == (ChannelStatistics *) NULL) ||
+      (density_y == (ChannelStatistics *) NULL) ||
+      (Q == (ChannelStatistics **) NULL) ||
+      (sum == (ChannelStatistics *) NULL))
     {
-      double
-        count;
-
-      (void) GetMatrixElement(accumulator,x,y,&count);
-      if (count >= (double) line_count)
+      if (Q != (ChannelStatistics **) NULL)
         {
-          double
-            maxima;
-
-          SegmentInfo
-            line;
-
-          ssize_t
-            v;
-
-          /*
-            Is point a local maxima?
-          */
-          maxima=count;
-          for (v=((ssize_t) -(height/2)); v < ((ssize_t) (height/2)); v++)
-          {
-            ssize_t
-              u;
-
-            for (u=((ssize_t) -(width/2)); u < ((ssize_t) (width/2)); u++)
-            {
-              if ((u != 0) || (v !=0))
-                {
-                  (void) GetMatrixElement(accumulator,x+u,y+v,&count);
-                  if (count > maxima)
-                    {
-                      maxima=count;
-                      break;
-                    }
-                }
-            }
-            if (u < (ssize_t) (width/2))
-              break;
-          }
-          (void) GetMatrixElement(accumulator,x,y,&count);
-          if (maxima > count)
-            continue;
-          if ((x >= 45) && (x <= 135))
-            {
-              /*
-                y = (r-x cos(t))/sin(t)
-              */
-              line.x1=0.0;
-              line.y1=((double) (y-(accumulator_height/2.0))-((line.x1-
-                (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
-                sin(DegreesToRadians((double) x))+(image->rows/2.0);
-              line.x2=(double) image->columns;
-              line.y2=((double) (y-(accumulator_height/2.0))-((line.x2-
-                (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
-                sin(DegreesToRadians((double) x))+(image->rows/2.0);
-            }
-          else
-            {
-              /*
-                x = (r-y cos(t))/sin(t)
-              */
-              line.y1=0.0;
-              line.x1=((double) (y-(accumulator_height/2.0))-((line.y1-
-                (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
-                cos(DegreesToRadians((double) x))+(image->columns/2.0);
-              line.y2=(double) image->rows;
-              line.x2=((double) (y-(accumulator_height/2.0))-((line.y2-
-                (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
-                cos(DegreesToRadians((double) x))+(image->columns/2.0);
-            }
-          (void) FormatLocaleString(message,MaxTextExtent,
-            "line %.20g,%.20g %.20g,%.20g  # %.20g\n",
-            (double) ((ssize_t) line.x1),(double) ((ssize_t) line.y1),
-            (double) ((ssize_t) line.x2),(double) ((ssize_t) line.y2),maxima);
-          (void) write(file,message,strlen(message));
+          for (i=0; i < (ssize_t) number_grays; i++)
+            Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
+          Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
+        }
+      if (sum != (ChannelStatistics *) NULL)
+        sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
+      if (density_y != (ChannelStatistics *) NULL)
+        density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
+      if (density_xy != (ChannelStatistics *) NULL)
+        density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
+      if (density_x != (ChannelStatistics *) NULL)
+        density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
+      if (cooccurrence != (ChannelStatistics **) NULL)
+        {
+          for (i=0; i < (ssize_t) number_grays; i++)
+            cooccurrence[i]=(ChannelStatistics *)
+              RelinquishMagickMemory(cooccurrence[i]);
+          cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
+            cooccurrence);
         }
+      grays=(PixelPacket *) RelinquishMagickMemory(grays);
+      channel_features=(ChannelFeatures *) RelinquishMagickMemory(
+        channel_features);
+      (void) ThrowMagickException(exception,GetMagickModule(),
+        ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
+      return(channel_features);
     }
+  (void) ResetMagickMemory(&correlation,0,sizeof(correlation));
+  (void) ResetMagickMemory(density_x,0,2*(number_grays+1)*sizeof(*density_x));
+  (void) ResetMagickMemory(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
+  (void) ResetMagickMemory(density_y,0,2*(number_grays+1)*sizeof(*density_y));
+  (void) ResetMagickMemory(&mean,0,sizeof(mean));
+  (void) ResetMagickMemory(sum,0,number_grays*sizeof(*sum));
+  (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
+  (void) ResetMagickMemory(density_xy,0,2*number_grays*sizeof(*density_xy));
+  (void) ResetMagickMemory(&entropy_x,0,sizeof(entropy_x));
+  (void) ResetMagickMemory(&entropy_xy,0,sizeof(entropy_xy));
+  (void) ResetMagickMemory(&entropy_xy1,0,sizeof(entropy_xy1));
+  (void) ResetMagickMemory(&entropy_xy2,0,sizeof(entropy_xy2));
+  (void) ResetMagickMemory(&entropy_y,0,sizeof(entropy_y));
+  (void) ResetMagickMemory(&variance,0,sizeof(variance));
+  for (i=0; i < (ssize_t) number_grays; i++)
+  {
+    cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
+      sizeof(**cooccurrence));
+    Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
+    if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
+        (Q[i] == (ChannelStatistics *) NULL))
+      break;
+    (void) ResetMagickMemory(cooccurrence[i],0,number_grays*
+      sizeof(**cooccurrence));
+    (void) ResetMagickMemory(Q[i],0,number_grays*sizeof(**Q));
   }
-  (void) close(file);
-  /*
-    Render lines to image canvas.
-  */
-  image_info=AcquireImageInfo();
-  image_info->background_color=image->background_color;
-  (void) FormatLocaleString(image_info->filename,MaxTextExtent,"mvg:%s",path);
-  artifact=GetImageArtifact(image,"background");
-  if (artifact != (const char *) NULL)
-    (void) SetImageOption(image_info,"background",artifact);
-  artifact=GetImageArtifact(image,"fill");
-  if (artifact != (const char *) NULL)
-    (void) SetImageOption(image_info,"fill",artifact);
-  artifact=GetImageArtifact(image,"stroke");
-  if (artifact != (const char *) NULL)
-    (void) SetImageOption(image_info,"stroke",artifact);
-  artifact=GetImageArtifact(image,"strokewidth");
-  if (artifact != (const char *) NULL)
-    (void) SetImageOption(image_info,"strokewidth",artifact);
-  lines_image=ReadImage(image_info,exception);
-  artifact=GetImageArtifact(image,"hough-lines:accumulator");
-  if ((lines_image != (Image *) NULL) &&
-      (IsStringTrue(artifact) != MagickFalse))
+  if (i < (ssize_t) number_grays)
     {
-      Image
-        *accumulator_image;
-
-      accumulator_image=MatrixToImage(accumulator,exception);
-      if (accumulator_image != (Image *) NULL)
-        AppendImageToList(&lines_image,accumulator_image);
+      for (i--; i >= 0; i--)
+      {
+        if (Q[i] != (ChannelStatistics *) NULL)
+          Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
+        if (cooccurrence[i] != (ChannelStatistics *) NULL)
+          cooccurrence[i]=(ChannelStatistics *)
+            RelinquishMagickMemory(cooccurrence[i]);
+      }
+      Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
+      cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
+      sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
+      density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
+      density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
+      density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
+      grays=(PixelPacket *) RelinquishMagickMemory(grays);
+      channel_features=(ChannelFeatures *) RelinquishMagickMemory(
+        channel_features);
+      (void) ThrowMagickException(exception,GetMagickModule(),
+        ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
+      return(channel_features);
     }
   /*
-    Free resources.
+    Initialize spatial dependence matrix.
   */
-  accumulator=DestroyMatrixInfo(accumulator);
-  image_info=DestroyImageInfo(image_info);
-  (void) RelinquishUniqueFileResource(path);
-  return(GetFirstImageInList(lines_image));
-}
-\f
-/*
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%   G e t I m a g e F e a t u r e s                                           %
-%                                                                             %
-%                                                                             %
-%                                                                             %
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%  GetImageFeatures() returns features for each channel in the image in
-%  each of four directions (horizontal, vertical, left and right diagonals)
-%  for the specified distance.  The features include the angular second
-%  moment, contrast, correlation, sum of squares: variance, inverse difference
-%  moment, sum average, sum varience, sum entropy, entropy, difference variance,%  difference entropy, information measures of correlation 1, information
-%  measures of correlation 2, and maximum correlation coefficient.  You can
-%  access the red channel contrast, for example, like this:
-%
-%      channel_features=GetImageFeatures(image,1,exception);
-%      contrast=channel_features[RedPixelChannel].contrast[0];
-%
-%  Use MagickRelinquishMemory() to free the features buffer.
-%
-%  The format of the GetImageFeatures method is:
-%
-%      ChannelFeatures *GetImageFeatures(const Image *image,
-%        const size_t distance,ExceptionInfo *exception)
-%
-%  A description of each parameter follows:
-%
-%    o image: the image.
-%
-%    o distance: the distance.
-%
-%    o exception: return any errors or warnings in this structure.
-%
-*/
-
-static inline ssize_t MagickAbsoluteValue(const ssize_t x)
-{
-  if (x < 0)
-    return(-x);
-  return(x);
-}
-
-static inline double MagickLog10(const double x)
-{
-#define Log10Epsilon  (1.0e-11)
-
- if (fabs(x) < Log10Epsilon)
-   return(log10(Log10Epsilon));
- return(log10(fabs(x)));
-}
-
-MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
-  const size_t distance,ExceptionInfo *exception)
-{
-  typedef struct _ChannelStatistics
+  status=MagickTrue;
+  image_view=AcquireVirtualCacheView(image,exception);
+  for (y=0; y < (ssize_t) image->rows; y++)
   {
-    PixelInfo
-      direction[4];  /* horizontal, vertical, left and right diagonals */
-  } ChannelStatistics;
+    register const Quantum
+      *restrict p;
 
-  CacheView
-    *image_view;
+    register ssize_t
+      x;
 
-  ChannelFeatures
-    *channel_features;
-
-  ChannelStatistics
-    **cooccurrence,
-    correlation,
-    *density_x,
-    *density_xy,
-    *density_y,
-    entropy_x,
-    entropy_xy,
-    entropy_xy1,
-    entropy_xy2,
-    entropy_y,
-    mean,
-    **Q,
-    *sum,
-    sum_squares,
-    variance;
-
-  PixelPacket
-    gray,
-    *grays;
-
-  MagickBooleanType
-    status;
-
-  register ssize_t
-    i;
-
-  size_t
-    length;
-
-  ssize_t
-    y;
-
-  unsigned int
-    number_grays;
-
-  assert(image != (Image *) NULL);
-  assert(image->signature == MagickSignature);
-  if (image->debug != MagickFalse)
-    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
-  if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
-    return((ChannelFeatures *) NULL);
-  length=CompositeChannels+1UL;
-  channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
-    sizeof(*channel_features));
-  if (channel_features == (ChannelFeatures *) NULL)
-    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
-  (void) ResetMagickMemory(channel_features,0,length*
-    sizeof(*channel_features));
-  /*
-    Form grays.
-  */
-  grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays));
-  if (grays == (PixelPacket *) NULL)
-    {
-      channel_features=(ChannelFeatures *) RelinquishMagickMemory(
-        channel_features);
-      (void) ThrowMagickException(exception,GetMagickModule(),
-        ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
-      return(channel_features);
-    }
-  for (i=0; i <= (ssize_t) MaxMap; i++)
-  {
-    grays[i].red=(~0U);
-    grays[i].green=(~0U);
-    grays[i].blue=(~0U);
-    grays[i].alpha=(~0U);
-    grays[i].black=(~0U);
-  }
-  status=MagickTrue;
-  image_view=AcquireVirtualCacheView(image,exception);
-#if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,4) shared(status) \
-    magick_threads(image,image,image->rows,1)
-#endif
-  for (y=0; y < (ssize_t) image->rows; y++)
-  {
-    register const Quantum
-      *restrict p;
-
-    register ssize_t
-      x;
-
-    if (status == MagickFalse)
-      continue;
-    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
-    if (p == (const Quantum *) NULL)
-      {
-        status=MagickFalse;
-        continue;
-      }
-    for (x=0; x < (ssize_t) image->columns; x++)
-    {
-      grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
-        ScaleQuantumToMap(GetPixelRed(image,p));
-      grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
-        ScaleQuantumToMap(GetPixelGreen(image,p));
-      grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
-        ScaleQuantumToMap(GetPixelBlue(image,p));
-      if (image->colorspace == CMYKColorspace)
-        grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
-          ScaleQuantumToMap(GetPixelBlack(image,p));
-      if (image->alpha_trait == BlendPixelTrait)
-        grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
-          ScaleQuantumToMap(GetPixelAlpha(image,p));
-      p+=GetPixelChannels(image);
-    }
-  }
-  image_view=DestroyCacheView(image_view);
-  if (status == MagickFalse)
-    {
-      grays=(PixelPacket *) RelinquishMagickMemory(grays);
-      channel_features=(ChannelFeatures *) RelinquishMagickMemory(
-        channel_features);
-      return(channel_features);
-    }
-  (void) ResetMagickMemory(&gray,0,sizeof(gray));
-  for (i=0; i <= (ssize_t) MaxMap; i++)
-  {
-    if (grays[i].red != ~0U)
-      grays[gray.red++].red=grays[i].red;
-    if (grays[i].green != ~0U)
-      grays[gray.green++].green=grays[i].green;
-    if (grays[i].blue != ~0U)
-      grays[gray.blue++].blue=grays[i].blue;
-    if (image->colorspace == CMYKColorspace)
-      if (grays[i].black != ~0U)
-        grays[gray.black++].black=grays[i].black;
-    if (image->alpha_trait == BlendPixelTrait)
-      if (grays[i].alpha != ~0U)
-        grays[gray.alpha++].alpha=grays[i].alpha;
-  }
-  /*
-    Allocate spatial dependence matrix.
-  */
-  number_grays=gray.red;
-  if (gray.green > number_grays)
-    number_grays=gray.green;
-  if (gray.blue > number_grays)
-    number_grays=gray.blue;
-  if (image->colorspace == CMYKColorspace)
-    if (gray.black > number_grays)
-      number_grays=gray.black;
-  if (image->alpha_trait == BlendPixelTrait)
-    if (gray.alpha > number_grays)
-      number_grays=gray.alpha;
-  cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays,
-    sizeof(*cooccurrence));
-  density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
-    sizeof(*density_x));
-  density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
-    sizeof(*density_xy));
-  density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
-    sizeof(*density_y));
-  Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q));
-  sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum));
-  if ((cooccurrence == (ChannelStatistics **) NULL) ||
-      (density_x == (ChannelStatistics *) NULL) ||
-      (density_xy == (ChannelStatistics *) NULL) ||
-      (density_y == (ChannelStatistics *) NULL) ||
-      (Q == (ChannelStatistics **) NULL) ||
-      (sum == (ChannelStatistics *) NULL))
-    {
-      if (Q != (ChannelStatistics **) NULL)
-        {
-          for (i=0; i < (ssize_t) number_grays; i++)
-            Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
-          Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
-        }
-      if (sum != (ChannelStatistics *) NULL)
-        sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
-      if (density_y != (ChannelStatistics *) NULL)
-        density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
-      if (density_xy != (ChannelStatistics *) NULL)
-        density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
-      if (density_x != (ChannelStatistics *) NULL)
-        density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
-      if (cooccurrence != (ChannelStatistics **) NULL)
-        {
-          for (i=0; i < (ssize_t) number_grays; i++)
-            cooccurrence[i]=(ChannelStatistics *)
-              RelinquishMagickMemory(cooccurrence[i]);
-          cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
-            cooccurrence);
-        }
-      grays=(PixelPacket *) RelinquishMagickMemory(grays);
-      channel_features=(ChannelFeatures *) RelinquishMagickMemory(
-        channel_features);
-      (void) ThrowMagickException(exception,GetMagickModule(),
-        ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
-      return(channel_features);
-    }
-  (void) ResetMagickMemory(&correlation,0,sizeof(correlation));
-  (void) ResetMagickMemory(density_x,0,2*(number_grays+1)*sizeof(*density_x));
-  (void) ResetMagickMemory(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
-  (void) ResetMagickMemory(density_y,0,2*(number_grays+1)*sizeof(*density_y));
-  (void) ResetMagickMemory(&mean,0,sizeof(mean));
-  (void) ResetMagickMemory(sum,0,number_grays*sizeof(*sum));
-  (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
-  (void) ResetMagickMemory(density_xy,0,2*number_grays*sizeof(*density_xy));
-  (void) ResetMagickMemory(&entropy_x,0,sizeof(entropy_x));
-  (void) ResetMagickMemory(&entropy_xy,0,sizeof(entropy_xy));
-  (void) ResetMagickMemory(&entropy_xy1,0,sizeof(entropy_xy1));
-  (void) ResetMagickMemory(&entropy_xy2,0,sizeof(entropy_xy2));
-  (void) ResetMagickMemory(&entropy_y,0,sizeof(entropy_y));
-  (void) ResetMagickMemory(&variance,0,sizeof(variance));
-  for (i=0; i < (ssize_t) number_grays; i++)
-  {
-    cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
-      sizeof(**cooccurrence));
-    Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
-    if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
-        (Q[i] == (ChannelStatistics *) NULL))
-      break;
-    (void) ResetMagickMemory(cooccurrence[i],0,number_grays*
-      sizeof(**cooccurrence));
-    (void) ResetMagickMemory(Q[i],0,number_grays*sizeof(**Q));
-  }
-  if (i < (ssize_t) number_grays)
-    {
-      for (i--; i >= 0; i--)
-      {
-        if (Q[i] != (ChannelStatistics *) NULL)
-          Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
-        if (cooccurrence[i] != (ChannelStatistics *) NULL)
-          cooccurrence[i]=(ChannelStatistics *)
-            RelinquishMagickMemory(cooccurrence[i]);
-      }
-      Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
-      cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
-      sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
-      density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
-      density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
-      density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
-      grays=(PixelPacket *) RelinquishMagickMemory(grays);
-      channel_features=(ChannelFeatures *) RelinquishMagickMemory(
-        channel_features);
-      (void) ThrowMagickException(exception,GetMagickModule(),
-        ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
-      return(channel_features);
-    }
-  /*
-    Initialize spatial dependence matrix.
-  */
-  status=MagickTrue;
-  image_view=AcquireVirtualCacheView(image,exception);
-  for (y=0; y < (ssize_t) image->rows; y++)
-  {
-    register const Quantum
-      *restrict p;
-
-    register ssize_t
-      x;
-
-    ssize_t
-      i,
-      offset,
-      u,
-      v;
+    ssize_t
+      i,
+      offset,
+      u,
+      v;
 
     if (status == MagickFalse)
       continue;
@@ -1264,7 +980,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
             cooccurrence[u][v].direction[i].black++;
             cooccurrence[v][u].direction[i].black++;
           }
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           {
             u=0;
             v=0;
@@ -1353,7 +1069,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         cooccurrence[x][y].direction[i].blue*=normalize;
         if (image->colorspace == CMYKColorspace)
           cooccurrence[x][y].direction[i].black*=normalize;
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           cooccurrence[x][y].direction[i].alpha*=normalize;
       }
     }
@@ -1393,7 +1109,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
           channel_features[BlackPixelChannel].angular_second_moment[i]+=
             cooccurrence[x][y].direction[i].black*
             cooccurrence[x][y].direction[i].black;
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           channel_features[AlphaPixelChannel].angular_second_moment[i]+=
             cooccurrence[x][y].direction[i].alpha*
             cooccurrence[x][y].direction[i].alpha;
@@ -1405,7 +1121,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
         if (image->colorspace == CMYKColorspace)
           sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black;
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha;
         correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
         correlation.direction[i].green+=x*y*
@@ -1415,7 +1131,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         if (image->colorspace == CMYKColorspace)
           correlation.direction[i].black+=x*y*
             cooccurrence[x][y].direction[i].black;
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           correlation.direction[i].alpha+=x*y*
             cooccurrence[x][y].direction[i].alpha;
         /*
@@ -1430,7 +1146,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         if (image->colorspace == CMYKColorspace)
           channel_features[BlackPixelChannel].inverse_difference_moment[i]+=
             cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1);
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           channel_features[AlphaPixelChannel].inverse_difference_moment[i]+=
             cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
         /*
@@ -1445,7 +1161,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         if (image->colorspace == CMYKColorspace)
           density_xy[y+x+2].direction[i].black+=
             cooccurrence[x][y].direction[i].black;
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           density_xy[y+x+2].direction[i].alpha+=
             cooccurrence[x][y].direction[i].alpha;
         /*
@@ -1464,7 +1180,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
           channel_features[BlackPixelChannel].entropy[i]-=
             cooccurrence[x][y].direction[i].black*
             MagickLog10(cooccurrence[x][y].direction[i].black);
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           channel_features[AlphaPixelChannel].entropy[i]-=
             cooccurrence[x][y].direction[i].alpha*
             MagickLog10(cooccurrence[x][y].direction[i].alpha);
@@ -1474,7 +1190,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
         density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
         density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           density_x[x].direction[i].alpha+=
             cooccurrence[x][y].direction[i].alpha;
         if (image->colorspace == CMYKColorspace)
@@ -1486,7 +1202,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         if (image->colorspace == CMYKColorspace)
           density_y[y].direction[i].black+=
             cooccurrence[x][y].direction[i].black;
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           density_y[y].direction[i].alpha+=
             cooccurrence[x][y].direction[i].alpha;
       }
@@ -1501,7 +1217,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
           mean.direction[i].black+=y*sum[y].direction[i].black;
           sum_squares.direction[i].black+=y*y*sum[y].direction[i].black;
         }
-      if (image->alpha_trait == BlendPixelTrait)
+      if (image->alpha_trait != UndefinedPixelTrait)
         {
           mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
           sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
@@ -1535,7 +1251,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         (mean.direction[i].black*mean.direction[i].black))*sqrt(
         sum_squares.direction[i].black-(mean.direction[i].black*
         mean.direction[i].black)));
-    if (image->alpha_trait == BlendPixelTrait)
+    if (image->alpha_trait != UndefinedPixelTrait)
       channel_features[AlphaPixelChannel].correlation[i]=
         (correlation.direction[i].alpha-mean.direction[i].alpha*
         mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
@@ -1569,7 +1285,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
       if (image->colorspace == CMYKColorspace)
         channel_features[BlackPixelChannel].sum_average[i]+=
           x*density_xy[x].direction[i].black;
-      if (image->alpha_trait == BlendPixelTrait)
+      if (image->alpha_trait != UndefinedPixelTrait)
         channel_features[AlphaPixelChannel].sum_average[i]+=
           x*density_xy[x].direction[i].alpha;
       /*
@@ -1588,7 +1304,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         channel_features[BlackPixelChannel].sum_entropy[i]-=
           density_xy[x].direction[i].black*
           MagickLog10(density_xy[x].direction[i].black);
-      if (image->alpha_trait == BlendPixelTrait)
+      if (image->alpha_trait != UndefinedPixelTrait)
         channel_features[AlphaPixelChannel].sum_entropy[i]-=
           density_xy[x].direction[i].alpha*
           MagickLog10(density_xy[x].direction[i].alpha);
@@ -1612,7 +1328,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
           (x-channel_features[BlackPixelChannel].sum_entropy[i])*
           (x-channel_features[BlackPixelChannel].sum_entropy[i])*
           density_xy[x].direction[i].black;
-      if (image->alpha_trait == BlendPixelTrait)
+      if (image->alpha_trait != UndefinedPixelTrait)
         channel_features[AlphaPixelChannel].sum_variance[i]+=
           (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
           (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
@@ -1650,7 +1366,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         if (image->colorspace == CMYKColorspace)
           variance.direction[i].black+=(y-mean.direction[i].black+1)*
             (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black;
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
             (y-mean.direction[i].alpha+1)*
             cooccurrence[x][y].direction[i].alpha;
@@ -1666,7 +1382,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         if (image->colorspace == CMYKColorspace)
           density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
             cooccurrence[x][y].direction[i].black;
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
             cooccurrence[x][y].direction[i].alpha;
         /*
@@ -1681,7 +1397,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         if (image->colorspace == CMYKColorspace)
           entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black*
             MagickLog10(cooccurrence[x][y].direction[i].black);
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           entropy_xy.direction[i].alpha-=
             cooccurrence[x][y].direction[i].alpha*MagickLog10(
             cooccurrence[x][y].direction[i].alpha);
@@ -1696,7 +1412,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
           entropy_xy1.direction[i].black-=(
             cooccurrence[x][y].direction[i].black*MagickLog10(
             density_x[x].direction[i].black*density_y[y].direction[i].black));
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           entropy_xy1.direction[i].alpha-=(
             cooccurrence[x][y].direction[i].alpha*MagickLog10(
             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
@@ -1713,7 +1429,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
           entropy_xy2.direction[i].black-=(density_x[x].direction[i].black*
             density_y[y].direction[i].black*MagickLog10(
             density_x[x].direction[i].black*density_y[y].direction[i].black));
-        if (image->alpha_trait == BlendPixelTrait)
+        if (image->alpha_trait != UndefinedPixelTrait)
           entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha*
             density_y[y].direction[i].alpha*MagickLog10(
             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
@@ -1728,7 +1444,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
     if (image->colorspace == CMYKColorspace)
       channel_features[BlackPixelChannel].variance_sum_of_squares[i]=
         variance.direction[i].black;
-    if (image->alpha_trait == BlendPixelTrait)
+    if (image->alpha_trait != UndefinedPixelTrait)
       channel_features[AlphaPixelChannel].variance_sum_of_squares[i]=
         variance.direction[i].alpha;
   }
@@ -1756,7 +1472,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
       variance.direction[i].blue+=density_xy[x].direction[i].blue;
       if (image->colorspace == CMYKColorspace)
         variance.direction[i].black+=density_xy[x].direction[i].black;
-      if (image->alpha_trait == BlendPixelTrait)
+      if (image->alpha_trait != UndefinedPixelTrait)
         variance.direction[i].alpha+=density_xy[x].direction[i].alpha;
       sum_squares.direction[i].red+=density_xy[x].direction[i].red*
         density_xy[x].direction[i].red;
@@ -1767,7 +1483,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
       if (image->colorspace == CMYKColorspace)
         sum_squares.direction[i].black+=density_xy[x].direction[i].black*
           density_xy[x].direction[i].black;
-      if (image->alpha_trait == BlendPixelTrait)
+      if (image->alpha_trait != UndefinedPixelTrait)
         sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
           density_xy[x].direction[i].alpha;
       /*
@@ -1786,7 +1502,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         channel_features[BlackPixelChannel].difference_entropy[i]-=
           density_xy[x].direction[i].black*
           MagickLog10(density_xy[x].direction[i].black);
-      if (image->alpha_trait == BlendPixelTrait)
+      if (image->alpha_trait != UndefinedPixelTrait)
         channel_features[AlphaPixelChannel].difference_entropy[i]-=
           density_xy[x].direction[i].alpha*
           MagickLog10(density_xy[x].direction[i].alpha);
@@ -1802,7 +1518,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
       if (image->colorspace == CMYKColorspace)
         entropy_x.direction[i].black-=(density_x[x].direction[i].black*
           MagickLog10(density_x[x].direction[i].black));
-      if (image->alpha_trait == BlendPixelTrait)
+      if (image->alpha_trait != UndefinedPixelTrait)
         entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha*
           MagickLog10(density_x[x].direction[i].alpha));
       entropy_y.direction[i].red-=(density_y[x].direction[i].red*
@@ -1814,7 +1530,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
       if (image->colorspace == CMYKColorspace)
         entropy_y.direction[i].black-=(density_y[x].direction[i].black*
           MagickLog10(density_y[x].direction[i].black));
-      if (image->alpha_trait == BlendPixelTrait)
+      if (image->alpha_trait != UndefinedPixelTrait)
         entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
           MagickLog10(density_y[x].direction[i].alpha));
     }
@@ -1838,7 +1554,7 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         (((double) number_grays*number_grays*sum_squares.direction[i].black)-
         (variance.direction[i].black*variance.direction[i].black))/
         ((double) number_grays*number_grays*number_grays*number_grays);
-    if (image->alpha_trait == BlendPixelTrait)
+    if (image->alpha_trait != UndefinedPixelTrait)
       channel_features[AlphaPixelChannel].difference_variance[i]=
         (((double) number_grays*number_grays*sum_squares.direction[i].alpha)-
         (variance.direction[i].alpha*variance.direction[i].alpha))/
@@ -1863,138 +1579,695 @@ MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
         (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/
         (entropy_x.direction[i].black > entropy_y.direction[i].black ?
          entropy_x.direction[i].black : entropy_y.direction[i].black);
-    if (image->alpha_trait == BlendPixelTrait)
+    if (image->alpha_trait != UndefinedPixelTrait)
       channel_features[AlphaPixelChannel].measure_of_correlation_1[i]=
         (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/
         (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ?
          entropy_x.direction[i].alpha : entropy_y.direction[i].alpha);
     channel_features[RedPixelChannel].measure_of_correlation_2[i]=
-      (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].red-
+      (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].red-
       entropy_xy.direction[i].red)))));
     channel_features[GreenPixelChannel].measure_of_correlation_2[i]=
-      (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].green-
+      (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].green-
       entropy_xy.direction[i].green)))));
     channel_features[BluePixelChannel].measure_of_correlation_2[i]=
-      (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].blue-
+      (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].blue-
       entropy_xy.direction[i].blue)))));
     if (image->colorspace == CMYKColorspace)
       channel_features[BlackPixelChannel].measure_of_correlation_2[i]=
-        (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].black-
+        (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].black-
         entropy_xy.direction[i].black)))));
-    if (image->alpha_trait == BlendPixelTrait)
+    if (image->alpha_trait != UndefinedPixelTrait)
       channel_features[AlphaPixelChannel].measure_of_correlation_2[i]=
-        (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].alpha-
+        (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].alpha-
         entropy_xy.direction[i].alpha)))));
   }
   /*
-    Compute more texture features.
+    Compute more texture features.
+  */
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+  #pragma omp parallel for schedule(static,4) shared(status) \
+    magick_threads(image,image,number_grays,1)
+#endif
+  for (i=0; i < 4; i++)
+  {
+    ssize_t
+      z;
+
+    for (z=0; z < (ssize_t) number_grays; z++)
+    {
+      register ssize_t
+        y;
+
+      ChannelStatistics
+        pixel;
+
+      (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
+      for (y=0; y < (ssize_t) number_grays; y++)
+      {
+        register ssize_t
+          x;
+
+        for (x=0; x < (ssize_t) number_grays; x++)
+        {
+          /*
+            Contrast:  amount of local variations present in an image.
+          */
+          if (((y-x) == z) || ((x-y) == z))
+            {
+              pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
+              pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
+              pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
+              if (image->colorspace == CMYKColorspace)
+                pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
+              if (image->alpha_trait != UndefinedPixelTrait)
+                pixel.direction[i].alpha+=
+                  cooccurrence[x][y].direction[i].alpha;
+            }
+          /*
+            Maximum Correlation Coefficient.
+          */
+          Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
+            cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
+            density_y[x].direction[i].red;
+          Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
+            cooccurrence[y][x].direction[i].green/
+            density_x[z].direction[i].green/density_y[x].direction[i].red;
+          Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
+            cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/
+            density_y[x].direction[i].blue;
+          if (image->colorspace == CMYKColorspace)
+            Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
+              cooccurrence[y][x].direction[i].black/
+              density_x[z].direction[i].black/density_y[x].direction[i].black;
+          if (image->alpha_trait != UndefinedPixelTrait)
+            Q[z][y].direction[i].alpha+=
+              cooccurrence[z][x].direction[i].alpha*
+              cooccurrence[y][x].direction[i].alpha/
+              density_x[z].direction[i].alpha/
+              density_y[x].direction[i].alpha;
+        }
+      }
+      channel_features[RedPixelChannel].contrast[i]+=z*z*
+        pixel.direction[i].red;
+      channel_features[GreenPixelChannel].contrast[i]+=z*z*
+        pixel.direction[i].green;
+      channel_features[BluePixelChannel].contrast[i]+=z*z*
+        pixel.direction[i].blue;
+      if (image->colorspace == CMYKColorspace)
+        channel_features[BlackPixelChannel].contrast[i]+=z*z*
+          pixel.direction[i].black;
+      if (image->alpha_trait != UndefinedPixelTrait)
+        channel_features[AlphaPixelChannel].contrast[i]+=z*z*
+          pixel.direction[i].alpha;
+    }
+    /*
+      Maximum Correlation Coefficient.
+      Future: return second largest eigenvalue of Q.
+    */
+    channel_features[RedPixelChannel].maximum_correlation_coefficient[i]=
+      sqrt((double) -1.0);
+    channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]=
+      sqrt((double) -1.0);
+    channel_features[BluePixelChannel].maximum_correlation_coefficient[i]=
+      sqrt((double) -1.0);
+    if (image->colorspace == CMYKColorspace)
+      channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]=
+        sqrt((double) -1.0);
+    if (image->alpha_trait != UndefinedPixelTrait)
+      channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]=
+        sqrt((double) -1.0);
+  }
+  /*
+    Relinquish resources.
+  */
+  sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
+  for (i=0; i < (ssize_t) number_grays; i++)
+    Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
+  Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
+  density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
+  density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
+  density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
+  for (i=0; i < (ssize_t) number_grays; i++)
+    cooccurrence[i]=(ChannelStatistics *)
+      RelinquishMagickMemory(cooccurrence[i]);
+  cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
+  return(channel_features);
+}
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%     H o u g h L i n e I m a g e                                             %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  Use HoughLineImage() in conjunction with any binary edge extracted image (we
+%  recommand Canny) to identify lines in the image.  The algorithm accumulates
+%  counts for every white pixel for every possible orientation (for angles from
+%  0 to 179 in 1 degree increments) and distance from the center of the image to
+%  the corner (in 1 px increments) and stores the counts in an accumulator matrix
+%  of angle vs distance. The size of the accumulator is 180x(diagonal/2). Next
+%  it searches this space for peaks in counts and converts the locations of the
+%  peaks to slope and intercept in the normal x,y input image space. Use the
+%  slope/intercepts to find the endpoints clipped to the bounds of the image. The
+%  lines are then drawn. The counts are a measure of the length of the lines
+%
+%  The format of the HoughLineImage method is:
+%
+%      Image *HoughLineImage(const Image *image,const size_t width,
+%        const size_t height,const size_t threshold,ExceptionInfo *exception)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image.
+%
+%    o width, height: find line pairs as local maxima in this neighborhood.
+%
+%    o threshold: the line count threshold.
+%
+%    o exception: return any errors or warnings in this structure.
+%
+*/
+
+static inline double MagickRound(double x)
+{
+  /*
+    Round the fraction to nearest integer.
+  */
+  if ((x-floor(x)) < (ceil(x)-x))
+    return(floor(x));
+  return(ceil(x));
+}
+
+MagickExport Image *HoughLineImage(const Image *image,const size_t width,
+  const size_t height,const size_t threshold,ExceptionInfo *exception)
+{
+#define HoughLineImageTag  "HoughLine/Image"
+
+  CacheView
+    *image_view;
+
+  char
+    message[MagickPathExtent],
+    path[MagickPathExtent];
+
+  const char
+    *artifact;
+
+  double
+    hough_height;
+
+  Image
+    *lines_image = NULL;
+
+  ImageInfo
+    *image_info;
+
+  int
+    file;
+
+  MagickBooleanType
+    status;
+
+  MagickOffsetType
+    progress;
+
+  MatrixInfo
+    *accumulator;
+
+  PointInfo
+    center;
+
+  register ssize_t
+    y;
+
+  size_t
+    accumulator_height,
+    accumulator_width,
+    line_count;
+
+  /*
+    Create the accumulator.
+  */
+  assert(image != (const Image *) NULL);
+  assert(image->signature == MagickCoreSignature);
+  if (image->debug != MagickFalse)
+    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
+  assert(exception != (ExceptionInfo *) NULL);
+  assert(exception->signature == MagickCoreSignature);
+  accumulator_width=180;
+  hough_height=((sqrt(2.0)*(double) (image->rows > image->columns ?
+    image->rows : image->columns))/2.0);
+  accumulator_height=(size_t) (2.0*hough_height);
+  accumulator=AcquireMatrixInfo(accumulator_width,accumulator_height,
+    sizeof(double),exception);
+  if (accumulator == (MatrixInfo *) NULL)
+    ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
+  if (NullMatrix(accumulator) == MagickFalse)
+    {
+      accumulator=DestroyMatrixInfo(accumulator);
+      ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
+    }
+  /*
+    Populate the accumulator.
+  */
+  status=MagickTrue;
+  progress=0;
+  center.x=(double) image->columns/2.0;
+  center.y=(double) image->rows/2.0;
+  image_view=AcquireVirtualCacheView(image,exception);
+  for (y=0; y < (ssize_t) image->rows; y++)
+  {
+    register const Quantum
+      *restrict p;
+
+    register ssize_t
+      x;
+
+    if (status == MagickFalse)
+      continue;
+    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
+    if (p == (Quantum *) NULL)
+      {
+        status=MagickFalse;
+        continue;
+      }
+    for (x=0; x < (ssize_t) image->columns; x++)
+    {
+      if (GetPixelIntensity(image,p) > (QuantumRange/2.0))
+        {
+          register ssize_t
+            i;
+
+          for (i=0; i < 180; i++)
+          {
+            double
+              count,
+              radius;
+
+            radius=(((double) x-center.x)*cos(DegreesToRadians((double) i)))+
+              (((double) y-center.y)*sin(DegreesToRadians((double) i)));
+            (void) GetMatrixElement(accumulator,i,(ssize_t)
+              MagickRound(radius+hough_height),&count);
+            count++;
+            (void) SetMatrixElement(accumulator,i,(ssize_t)
+              MagickRound(radius+hough_height),&count);
+          }
+        }
+      p+=GetPixelChannels(image);
+    }
+    if (image->progress_monitor != (MagickProgressMonitor) NULL)
+      {
+        MagickBooleanType
+          proceed;
+
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+        #pragma omp critical (MagickCore_CannyEdgeImage)
+#endif
+        proceed=SetImageProgress(image,CannyEdgeImageTag,progress++,
+          image->rows);
+        if (proceed == MagickFalse)
+          status=MagickFalse;
+      }
+  }
+  image_view=DestroyCacheView(image_view);
+  if (status == MagickFalse)
+    {
+      accumulator=DestroyMatrixInfo(accumulator);
+      return((Image *) NULL);
+    }
+  /*
+    Generate line segments from accumulator.
+  */
+  file=AcquireUniqueFileResource(path);
+  if (file == -1)
+    {
+      accumulator=DestroyMatrixInfo(accumulator);
+      return((Image *) NULL);
+    }
+  (void) FormatLocaleString(message,MagickPathExtent,
+    "# Hough line transform: %.20gx%.20g%+.20g\n",(double) width,
+    (double) height,(double) threshold);
+  if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
+    status=MagickFalse;
+  (void) FormatLocaleString(message,MagickPathExtent,"viewbox 0 0 %.20g %.20g\n",
+    (double) image->columns,(double) image->rows);
+  if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
+    status=MagickFalse;
+  line_count=image->columns > image->rows ? image->columns/4 : image->rows/4;
+  if (threshold != 0)
+    line_count=threshold;
+  for (y=0; y < (ssize_t) accumulator_height; y++)
+  {
+    register ssize_t
+      x;
+
+    for (x=0; x < (ssize_t) accumulator_width; x++)
+    {
+      double
+        count;
+
+      (void) GetMatrixElement(accumulator,x,y,&count);
+      if (count >= (double) line_count)
+        {
+          double
+            maxima;
+
+          SegmentInfo
+            line;
+
+          ssize_t
+            v;
+
+          /*
+            Is point a local maxima?
+          */
+          maxima=count;
+          for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
+          {
+            ssize_t
+              u;
+
+            for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
+            {
+              if ((u != 0) || (v !=0))
+                {
+                  (void) GetMatrixElement(accumulator,x+u,y+v,&count);
+                  if (count > maxima)
+                    {
+                      maxima=count;
+                      break;
+                    }
+                }
+            }
+            if (u < (ssize_t) (width/2))
+              break;
+          }
+          (void) GetMatrixElement(accumulator,x,y,&count);
+          if (maxima > count)
+            continue;
+          if ((x >= 45) && (x <= 135))
+            {
+              /*
+                y = (r-x cos(t))/sin(t)
+              */
+              line.x1=0.0;
+              line.y1=((double) (y-(accumulator_height/2.0))-((line.x1-
+                (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
+                sin(DegreesToRadians((double) x))+(image->rows/2.0);
+              line.x2=(double) image->columns;
+              line.y2=((double) (y-(accumulator_height/2.0))-((line.x2-
+                (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
+                sin(DegreesToRadians((double) x))+(image->rows/2.0);
+            }
+          else
+            {
+              /*
+                x = (r-y cos(t))/sin(t)
+              */
+              line.y1=0.0;
+              line.x1=((double) (y-(accumulator_height/2.0))-((line.y1-
+                (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
+                cos(DegreesToRadians((double) x))+(image->columns/2.0);
+              line.y2=(double) image->rows;
+              line.x2=((double) (y-(accumulator_height/2.0))-((line.y2-
+                (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
+                cos(DegreesToRadians((double) x))+(image->columns/2.0);
+            }
+          (void) FormatLocaleString(message,MagickPathExtent,
+            "line %g,%g %g,%g  # %g\n",line.x1,line.y1,line.x2,line.y2,maxima);
+          if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
+            status=MagickFalse;
+        }
+    }
+  }
+  (void) close(file);
+  /*
+    Render lines to image canvas.
+  */
+  image_info=AcquireImageInfo();
+  image_info->background_color=image->background_color;
+  (void) FormatLocaleString(image_info->filename,MagickPathExtent,"mvg:%s",path);
+  artifact=GetImageArtifact(image,"background");
+  if (artifact != (const char *) NULL)
+    (void) SetImageOption(image_info,"background",artifact);
+  artifact=GetImageArtifact(image,"fill");
+  if (artifact != (const char *) NULL)
+    (void) SetImageOption(image_info,"fill",artifact);
+  artifact=GetImageArtifact(image,"stroke");
+  if (artifact != (const char *) NULL)
+    (void) SetImageOption(image_info,"stroke",artifact);
+  artifact=GetImageArtifact(image,"strokewidth");
+  if (artifact != (const char *) NULL)
+    (void) SetImageOption(image_info,"strokewidth",artifact);
+  lines_image=ReadImage(image_info,exception);
+  artifact=GetImageArtifact(image,"hough-lines:accumulator");
+  if ((lines_image != (Image *) NULL) &&
+      (IsStringTrue(artifact) != MagickFalse))
+    {
+      Image
+        *accumulator_image;
+
+      accumulator_image=MatrixToImage(accumulator,exception);
+      if (accumulator_image != (Image *) NULL)
+        AppendImageToList(&lines_image,accumulator_image);
+    }
+  /*
+    Free resources.
   */
+  accumulator=DestroyMatrixInfo(accumulator);
+  image_info=DestroyImageInfo(image_info);
+  (void) RelinquishUniqueFileResource(path);
+  return(GetFirstImageInList(lines_image));
+}
+\f
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%     M e a n S h i f t I m a g e                                             %
+%                                                                             %
+%                                                                             %
+%                                                                             %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%  MeanShiftImage() delineate arbitrarily shaped clusters in the image. For
+%  each pixel, it visits all the pixels in the neighborhood specified by
+%  the window centered at the pixel and excludes those that are outside the
+%  radius=(window-1)/2 surrounding the pixel. From those pixels, it finds those
+%  that are within the specified color distance from the current mean, and
+%  computes a new x,y centroid from those coordinates and a new mean. This new
+%  x,y centroid is used as the center for a new window. This process iterates
+%  until it converges and the final mean is replaces the (original window
+%  center) pixel value. It repeats this process for the next pixel, etc., 
+%  until it processes all pixels in the image. Results are typically better with
+%  colorspaces other than sRGB. We recommend YIQ, YUV or YCbCr.
+%
+%  The format of the MeanShiftImage method is:
+%
+%      Image *MeanShiftImage(const Image *image,const size_t width,
+%        const size_t height,const double color_distance,
+%        ExceptionInfo *exception)
+%
+%  A description of each parameter follows:
+%
+%    o image: the image.
+%
+%    o width, height: find pixels in this neighborhood.
+%
+%    o color_distance: the color distance.
+%
+%    o exception: return any errors or warnings in this structure.
+%
+*/
+MagickExport Image *MeanShiftImage(const Image *image,const size_t width,
+  const size_t height,const double color_distance,ExceptionInfo *exception)
+{
+#define MaxMeanShiftIterations  100
+#define MeanShiftImageTag  "MeanShift/Image"
+
+  CacheView
+    *image_view,
+    *mean_view,
+    *pixel_view;
+
+  Image
+    *mean_image;
+
+  MagickBooleanType
+    status;
+
+  MagickOffsetType
+    progress;
+
+  ssize_t
+    y;
+
+  assert(image != (const Image *) NULL);
+  assert(image->signature == MagickCoreSignature);
+  if (image->debug != MagickFalse)
+    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
+  assert(exception != (ExceptionInfo *) NULL);
+  assert(exception->signature == MagickCoreSignature);
+  mean_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
+  if (mean_image == (Image *) NULL)
+    return((Image *) NULL);
+  if (SetImageStorageClass(mean_image,DirectClass,exception) == MagickFalse)
+    {
+      mean_image=DestroyImage(mean_image);
+      return((Image *) NULL);
+    }
+  status=MagickTrue;
+  progress=0;
+  image_view=AcquireVirtualCacheView(image,exception);
+  pixel_view=AcquireVirtualCacheView(image,exception);
+  mean_view=AcquireAuthenticCacheView(mean_image,exception);
 #if defined(MAGICKCORE_OPENMP_SUPPORT)
-  #pragma omp parallel for schedule(static,4) shared(status) \
-    magick_threads(image,image,number_grays,1)
+  #pragma omp parallel for schedule(static,4) shared(status,progress) \
+    magick_threads(mean_image,mean_image,mean_image->rows,1)
 #endif
-  for (i=0; i < 4; i++)
+  for (y=0; y < (ssize_t) mean_image->rows; y++)
   {
-    ssize_t
-      z;
+    register const Quantum
+      *restrict p;
 
-    for (z=0; z < (ssize_t) number_grays; z++)
+    register Quantum
+      *restrict q;
+
+    register ssize_t
+      x;
+
+    if (status == MagickFalse)
+      continue;
+    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
+    q=GetCacheViewAuthenticPixels(mean_view,0,y,mean_image->columns,1,
+      exception);
+    if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
+      {
+        status=MagickFalse;
+        continue;
+      }
+    for (x=0; x < (ssize_t) mean_image->columns; x++)
     {
-      register ssize_t
-        y;
+      PixelInfo
+        mean_pixel,
+        previous_pixel;
 
-      ChannelStatistics
-        pixel;
+      PointInfo
+        mean_location,
+        previous_location;
 
-      (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
-      for (y=0; y < (ssize_t) number_grays; y++)
+      register ssize_t
+        i;
+
+      GetPixelInfo(image,&mean_pixel);
+      GetPixelInfoPixel(image,p,&mean_pixel);
+      mean_location.x=(double) x;
+      mean_location.y=(double) y;
+      for (i=0; i < MaxMeanShiftIterations; i++)
       {
-        register ssize_t
-          x;
+        double
+          distance,
+          gamma;
 
-        for (x=0; x < (ssize_t) number_grays; x++)
+        PixelInfo
+          sum_pixel;
+
+        PointInfo
+          sum_location;
+
+        ssize_t
+          count,
+          v;
+
+        sum_location.x=0.0;
+        sum_location.y=0.0;
+        GetPixelInfo(image,&sum_pixel);
+        previous_location=mean_location;
+        previous_pixel=mean_pixel;
+        count=0;
+        for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
         {
-          /*
-            Contrast:  amount of local variations present in an image.
-          */
-          if (((y-x) == z) || ((x-y) == z))
-            {
-              pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
-              pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
-              pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
-              if (image->colorspace == CMYKColorspace)
-                pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
-              if (image->alpha_trait == BlendPixelTrait)
-                pixel.direction[i].alpha+=
-                  cooccurrence[x][y].direction[i].alpha;
-            }
-          /*
-            Maximum Correlation Coefficient.
-          */
-          Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
-            cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
-            density_y[x].direction[i].red;
-          Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
-            cooccurrence[y][x].direction[i].green/
-            density_x[z].direction[i].green/density_y[x].direction[i].red;
-          Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
-            cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/
-            density_y[x].direction[i].blue;
-          if (image->colorspace == CMYKColorspace)
-            Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
-              cooccurrence[y][x].direction[i].black/
-              density_x[z].direction[i].black/density_y[x].direction[i].black;
-          if (image->alpha_trait == BlendPixelTrait)
-            Q[z][y].direction[i].alpha+=
-              cooccurrence[z][x].direction[i].alpha*
-              cooccurrence[y][x].direction[i].alpha/
-              density_x[z].direction[i].alpha/
-              density_y[x].direction[i].alpha;
+          ssize_t
+            u;
+
+          for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
+          {
+            if ((v*v+u*u) <= (ssize_t) ((width/2)*(height/2)))
+              {
+                PixelInfo
+                  pixel;
+
+                status=GetOneCacheViewVirtualPixelInfo(pixel_view,(ssize_t)
+                  MagickRound(mean_location.x+u),(ssize_t) MagickRound(
+                  mean_location.y+v),&pixel,exception);
+                distance=(mean_pixel.red-pixel.red)*(mean_pixel.red-pixel.red)+
+                  (mean_pixel.green-pixel.green)*(mean_pixel.green-pixel.green)+
+                  (mean_pixel.blue-pixel.blue)*(mean_pixel.blue-pixel.blue);
+                if (distance <= (color_distance*color_distance))
+                  {
+                    sum_location.x+=mean_location.x+u;
+                    sum_location.y+=mean_location.y+v;
+                    sum_pixel.red+=pixel.red;
+                    sum_pixel.green+=pixel.green;
+                    sum_pixel.blue+=pixel.blue;
+                    sum_pixel.alpha+=pixel.alpha;
+                    count++;
+                  }
+              }
+          }
         }
+        gamma=1.0/count;
+        mean_location.x=gamma*sum_location.x;
+        mean_location.y=gamma*sum_location.y;
+        mean_pixel.red=gamma*sum_pixel.red;
+        mean_pixel.green=gamma*sum_pixel.green;
+        mean_pixel.blue=gamma*sum_pixel.blue;
+        mean_pixel.alpha=gamma*sum_pixel.alpha;
+        distance=(mean_location.x-previous_location.x)*
+          (mean_location.x-previous_location.x)+
+          (mean_location.y-previous_location.y)*
+          (mean_location.y-previous_location.y)+
+          255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)*
+          255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)+
+          255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)*
+          255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)+
+          255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue)*
+          255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue);
+        if (distance <= 3.0)
+          break;
       }
-      channel_features[RedPixelChannel].contrast[i]+=z*z*
-        pixel.direction[i].red;
-      channel_features[GreenPixelChannel].contrast[i]+=z*z*
-        pixel.direction[i].green;
-      channel_features[BluePixelChannel].contrast[i]+=z*z*
-        pixel.direction[i].blue;
-      if (image->colorspace == CMYKColorspace)
-        channel_features[BlackPixelChannel].contrast[i]+=z*z*
-          pixel.direction[i].black;
-      if (image->alpha_trait == BlendPixelTrait)
-        channel_features[AlphaPixelChannel].contrast[i]+=z*z*
-          pixel.direction[i].alpha;
+      SetPixelRed(mean_image,ClampToQuantum(mean_pixel.red),q);
+      SetPixelGreen(mean_image,ClampToQuantum(mean_pixel.green),q);
+      SetPixelBlue(mean_image,ClampToQuantum(mean_pixel.blue),q);
+      SetPixelAlpha(mean_image,ClampToQuantum(mean_pixel.alpha),q);
+      p+=GetPixelChannels(image);
+      q+=GetPixelChannels(mean_image);
     }
-    /*
-      Maximum Correlation Coefficient.
-      Future: return second largest eigenvalue of Q.
-    */
-    channel_features[RedPixelChannel].maximum_correlation_coefficient[i]=
-      sqrt((double) -1.0);
-    channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]=
-      sqrt((double) -1.0);
-    channel_features[BluePixelChannel].maximum_correlation_coefficient[i]=
-      sqrt((double) -1.0);
-    if (image->colorspace == CMYKColorspace)
-      channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]=
-        sqrt((double) -1.0);
-    if (image->alpha_trait == BlendPixelTrait)
-      channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]=
-        sqrt((double) -1.0);
+    if (SyncCacheViewAuthenticPixels(mean_view,exception) == MagickFalse)
+      status=MagickFalse;
+    if (image->progress_monitor != (MagickProgressMonitor) NULL)
+      {
+        MagickBooleanType
+          proceed;
+
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+        #pragma omp critical (MagickCore_MeanShiftImage)
+#endif
+        proceed=SetImageProgress(image,MeanShiftImageTag,progress++,
+          image->rows);
+        if (proceed == MagickFalse)
+          status=MagickFalse;
+      }
   }
-  /*
-    Relinquish resources.
-  */
-  sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
-  for (i=0; i < (ssize_t) number_grays; i++)
-    Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
-  Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
-  density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
-  density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
-  density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
-  for (i=0; i < (ssize_t) number_grays; i++)
-    cooccurrence[i]=(ChannelStatistics *)
-      RelinquishMagickMemory(cooccurrence[i]);
-  cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
-  return(channel_features);
+  mean_view=DestroyCacheView(mean_view);
+  pixel_view=DestroyCacheView(pixel_view);
+  image_view=DestroyCacheView(image_view);
+  return(mean_image);
 }