cppcheck/test/bug-hunting/cve/CVE-2019-14981/feature.c

/*
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %                                                                             %
   %                                                                             %
   %                                                                             %
   %               FFFFF  EEEEE   AAA   TTTTT  U   U  RRRR   EEEEE               %
   %               F      E      A   A    T    U   U  R   R  E                   %
   %               FFF    EEE    AAAAA    T    U   U  RRRR   EEE                 %
   %               F      E      A   A    T    U   U  R R    E                   %
   %               F      EEEEE  A   A    T     UUU   R  R   EEEEE               %
   %                                                                             %
   %                                                                             %
   %                      MagickCore Image Feature Methods                       %
   %                                                                             %
   %                              Software Design                                %
   %                                   Cristy                                    %
   %                                 July 1992                                   %
   %                                                                             %
   %                                                                             %
   %  Copyright 1999-2019 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                                            %
   %                                                                             %
   %    https://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.                                             %
   %                                                                             %
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %
   %
   %
 */

/*
   Include declarations.
 */
#include "MagickCore/studio.h"
#include "MagickCore/animate.h"
#include "MagickCore/artifact.h"
#include "MagickCore/blob.h"
#include "MagickCore/blob-private.h"
#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"
#include "MagickCore/colorspace.h"
#include "MagickCore/colorspace-private.h"
#include "MagickCore/composite.h"
#include "MagickCore/composite-private.h"
#include "MagickCore/compress.h"
#include "MagickCore/constitute.h"
#include "MagickCore/display.h"
#include "MagickCore/draw.h"
#include "MagickCore/enhance.h"
#include "MagickCore/exception.h"
#include "MagickCore/exception-private.h"
#include "MagickCore/feature.h"
#include "MagickCore/gem.h"
#include "MagickCore/geometry.h"
#include "MagickCore/list.h"
#include "MagickCore/image-private.h"
#include "MagickCore/magic.h"
#include "MagickCore/magick.h"
#include "MagickCore/matrix.h"
#include "MagickCore/memory_.h"
#include "MagickCore/module.h"
#include "MagickCore/monitor.h"
#include "MagickCore/monitor-private.h"
#include "MagickCore/morphology-private.h"
#include "MagickCore/option.h"
#include "MagickCore/paint.h"
#include "MagickCore/pixel-accessor.h"
#include "MagickCore/profile.h"
#include "MagickCore/property.h"
#include "MagickCore/quantize.h"
#include "MagickCore/quantum-private.h"
#include "MagickCore/random_.h"
#include "MagickCore/resource_.h"
#include "MagickCore/segment.h"
#include "MagickCore/semaphore.h"
#include "MagickCore/signature-private.h"
#include "MagickCore/string_.h"
#include "MagickCore/thread-private.h"
#include "MagickCore/timer.h"
#include "MagickCore/utility.h"
#include "MagickCore/version.h"

/*
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %                                                                             %
   %                                                                             %
   %                                                                             %
   %     C a n n y E d g e I m a g e                                             %
   %                                                                             %
   %                                                                             %
   %                                                                             %
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %
   %  CannyEdgeImage() uses a multi-stage algorithm to detect a wide range of
   %  edges in images.
   %
   %  The format of the CannyEdgeImage method is:
   %
   %      Image *CannyEdgeImage(const Image *image,const double radius,
   %        const double sigma,const double lower_percent,
   %        const double upper_percent,ExceptionInfo *exception)
   %
   %  A description of each parameter follows:
   %
   %    o image: the image.
   %
   %    o radius: the radius of the gaussian smoothing filter.
   %
   %    o sigma: the sigma of the gaussian smoothing filter.
   %
   %    o lower_percent: percentage of edge pixels in the lower threshold.
   %
   %    o upper_percent: percentage of edge pixels in the upper threshold.
   %
   %    o exception: return any errors or warnings in this structure.
   %
 */

typedef struct _CannyInfo
{
    double
        magnitude,
        intensity;

    int
        orientation;

    ssize_t
        x,
        y;
} CannyInfo;

static inline MagickBooleanType IsAuthenticPixel(const Image *image,
                                                 const ssize_t x,const ssize_t y)
{
    if ((x < 0) || (x >= (ssize_t) image->columns))
        return(MagickFalse);
    if ((y < 0) || (y >= (ssize_t) image->rows))
        return(MagickFalse);
    return(MagickTrue);
}

static MagickBooleanType TraceEdges(Image *edge_image,CacheView *edge_view,
                                    MatrixInfo *canny_cache,const ssize_t x,const ssize_t y,
                                    const double lower_threshold,ExceptionInfo *exception)
{
    CannyInfo
        edge,
        pixel;

    MagickBooleanType
        status;

    register Quantum
    *q;

    register ssize_t
        i;

    q=GetCacheViewAuthenticPixels(edge_view,x,y,1,1,exception);
    if (q == (Quantum *) NULL)
        return(MagickFalse);
    *q=QuantumRange;
    status=SyncCacheViewAuthenticPixels(edge_view,exception);
    if (status == MagickFalse)
        return(MagickFalse);
    if (GetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
        return(MagickFalse);
    edge.x=x;
    edge.y=y;
    if (SetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
        return(MagickFalse);
    for (i=1; i != 0;)
    {
        ssize_t
            v;

        i--;
        status=GetMatrixElement(canny_cache,i,0,&edge);
        if (status == MagickFalse)
            return(MagickFalse);
        for (v=(-1); v <= 1; v++)
        {
            ssize_t
                u;

            for (u=(-1); u <= 1; u++)
            {
                if ((u == 0) && (v == 0))
                    continue;
                if (IsAuthenticPixel(edge_image,edge.x+u,edge.y+v) == MagickFalse)
                    continue;
                /*
                   Not an edge if gradient value is below the lower threshold.
                 */
                q=GetCacheViewAuthenticPixels(edge_view,edge.x+u,edge.y+v,1,1,
                                              exception);
                if (q == (Quantum *) NULL)
                    return(MagickFalse);
                status=GetMatrixElement(canny_cache,edge.x+u,edge.y+v,&pixel);
                if (status == MagickFalse)
                    return(MagickFalse);
                if ((GetPixelIntensity(edge_image,q) == 0.0) &&
                    (pixel.intensity >= lower_threshold))
                {
                    *q=QuantumRange;
                    status=SyncCacheViewAuthenticPixels(edge_view,exception);
                    if (status == MagickFalse)
                        return(MagickFalse);
                    edge.x+=u;
                    edge.y+=v;
                    status=SetMatrixElement(canny_cache,i,0,&edge);
                    if (status == MagickFalse)
                        return(MagickFalse);
                    i++;
                }
            }
        }
    }
    return(MagickTrue);
}

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;

    CannyInfo
        element;

    char
        geometry[MagickPathExtent];

    double
        lower_threshold,
        max,
        min,
        upper_threshold;

    Image
    *edge_image;

    KernelInfo
    *kernel_info;

    MagickBooleanType
        status;

    MagickOffsetType
        progress;

    MatrixInfo
    *canny_cache;

    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);
    /*
       Filter out noise.
     */
    (void) FormatLocaleString(geometry,MagickPathExtent,
                              "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
    kernel_info=AcquireKernelInfo(geometry,exception);
    if (kernel_info == (KernelInfo *) NULL)
        ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
    edge_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,exception);
    kernel_info=DestroyKernelInfo(kernel_info);
    if (edge_image == (Image *) NULL)
        return((Image *) NULL);
    if (TransformImageColorspace(edge_image,GRAYColorspace,exception) == MagickFalse)
    {
        edge_image=DestroyImage(edge_image);
        return((Image *) NULL);
    }
    (void) SetImageAlphaChannel(edge_image,OffAlphaChannel,exception);
    /*
       Find the intensity gradient of the image.
     */
    canny_cache=AcquireMatrixInfo(edge_image->columns,edge_image->rows,
                                  sizeof(CannyInfo),exception);
    if (canny_cache == (MatrixInfo *) NULL)
    {
        edge_image=DestroyImage(edge_image);
        return((Image *) NULL);
    }
    status=MagickTrue;
    edge_view=AcquireVirtualCacheView(edge_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(edge_image,edge_image,edge_image->rows,1)
#endif
    for (y=0; y < (ssize_t) edge_image->rows; y++)
    {
        register const Quantum
        *magick_restrict p;

        register ssize_t
            x;

        if (status == MagickFalse)
            continue;
        p=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns+1,2,
                                    exception);
        if (p == (const Quantum *) NULL)
        {
            status=MagickFalse;
            continue;
        }
        for (x=0; x < (ssize_t) edge_image->columns; x++)
        {
            CannyInfo
                pixel;

            double
                dx,
                dy;

            register const Quantum
            *magick_restrict kernel_pixels;

            ssize_t
                v;

            static double
                Gx[2][2] =
            {
                { -1.0,  +1.0 },
                { -1.0,  +1.0 }
            },
                Gy[2][2] =
            {
                { +1.0, +1.0 },
                { -1.0, -1.0 }
            };

            (void) memset(&pixel,0,sizeof(pixel));
            dx=0.0;
            dy=0.0;
            kernel_pixels=p;
            for (v=0; v < 2; v++)
            {
                ssize_t
                    u;

                for (u=0; u < 2; u++)
                {
                    double
                        intensity;

                    intensity=GetPixelIntensity(edge_image,kernel_pixels+u);
                    dx+=0.5*Gx[v][u]*intensity;
                    dy+=0.5*Gy[v][u]*intensity;
                }
                kernel_pixels+=edge_image->columns+1;
            }
            pixel.magnitude=hypot(dx,dy);
            pixel.orientation=0;
            if (fabs(dx) > MagickEpsilon)
            {
                double
                    slope;

                slope=dy/dx;
                if (slope < 0.0)
                {
                    if (slope < -2.41421356237)
                        pixel.orientation=0;
                    else
                    if (slope < -0.414213562373)
                        pixel.orientation=1;
                    else
                        pixel.orientation=2;
                }
                else
                {
                    if (slope > 2.41421356237)
                        pixel.orientation=0;
                    else
                    if (slope > 0.414213562373)
                        pixel.orientation=3;
                    else
                        pixel.orientation=2;
                }
            }
            if (SetMatrixElement(canny_cache,x,y,&pixel) == MagickFalse)
                continue;
            p+=GetPixelChannels(edge_image);
        }
    }
    edge_view=DestroyCacheView(edge_view);
    /*
       Non-maxima suppression, remove pixels that are not considered to be part
       of an edge.
     */
    progress=0;
    (void) GetMatrixElement(canny_cache,0,0,&element);
    max=element.intensity;
    min=element.intensity;
    edge_view=AcquireAuthenticCacheView(edge_image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(edge_image,edge_image,edge_image->rows,1)
#endif
    for (y=0; y < (ssize_t) edge_image->rows; y++)
    {
        register Quantum
        *magick_restrict q;

        register ssize_t
            x;

        if (status == MagickFalse)
            continue;
        q=GetCacheViewAuthenticPixels(edge_view,0,y,edge_image->columns,1,
                                      exception);
        if (q == (Quantum *) NULL)
        {
            status=MagickFalse;
            continue;
        }
        for (x=0; x < (ssize_t) edge_image->columns; x++)
        {
            CannyInfo
                alpha_pixel,
                beta_pixel,
                pixel;

            (void) GetMatrixElement(canny_cache,x,y,&pixel);
            switch (pixel.orientation)
            {
            case 0:
            default:
            {
                /*
                   0 degrees, north and south.
                 */
                (void) GetMatrixElement(canny_cache,x,y-1,&alpha_pixel);
                (void) GetMatrixElement(canny_cache,x,y+1,&beta_pixel);
                break;
            }
            case 1:
            {
                /*
                   45 degrees, northwest and southeast.
                 */
                (void) GetMatrixElement(canny_cache,x-1,y-1,&alpha_pixel);
                (void) GetMatrixElement(canny_cache,x+1,y+1,&beta_pixel);
                break;
            }
            case 2:
            {
                /*
                   90 degrees, east and west.
                 */
                (void) GetMatrixElement(canny_cache,x-1,y,&alpha_pixel);
                (void) GetMatrixElement(canny_cache,x+1,y,&beta_pixel);
                break;
            }
            case 3:
            {
                /*
                   135 degrees, northeast and southwest.
                 */
                (void) GetMatrixElement(canny_cache,x+1,y-1,&beta_pixel);
                (void) GetMatrixElement(canny_cache,x-1,y+1,&alpha_pixel);
                break;
            }
            }
            pixel.intensity=pixel.magnitude;
            if ((pixel.magnitude < alpha_pixel.magnitude) ||
                (pixel.magnitude < beta_pixel.magnitude))
                pixel.intensity=0;
            (void) SetMatrixElement(canny_cache,x,y,&pixel);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
      #pragma omp critical (MagickCore_CannyEdgeImage)
#endif
            {
                if (pixel.intensity < min)
                    min=pixel.intensity;
                if (pixel.intensity > max)
                    max=pixel.intensity;
            }
            *q=0;
            q+=GetPixelChannels(edge_image);
        }
        if (SyncCacheViewAuthenticPixels(edge_view,exception) == MagickFalse)
            status=MagickFalse;
    }
    edge_view=DestroyCacheView(edge_view);
    /*
       Estimate hysteresis threshold.
     */
    lower_threshold=lower_percent*(max-min)+min;
    upper_threshold=upper_percent*(max-min)+min;
    /*
       Hysteresis threshold.
     */
    edge_view=AcquireAuthenticCacheView(edge_image,exception);
    for (y=0; y < (ssize_t) edge_image->rows; y++)
    {
        register ssize_t
            x;

        if (status == MagickFalse)
            continue;
        for (x=0; x < (ssize_t) edge_image->columns; x++)
        {
            CannyInfo
                pixel;

            register const Quantum
            *magick_restrict p;

            /*
               Edge if pixel gradient higher than upper threshold.
             */
            p=GetCacheViewVirtualPixels(edge_view,x,y,1,1,exception);
            if (p == (const Quantum *) NULL)
                continue;
            status=GetMatrixElement(canny_cache,x,y,&pixel);
            if (status == MagickFalse)
                continue;
            if ((GetPixelIntensity(edge_image,p) == 0.0) &&
                (pixel.intensity >= upper_threshold))
                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 atomic
#endif
            progress++;
            proceed=SetImageProgress(image,CannyEdgeImageTag,progress,image->rows);
            if (proceed == MagickFalse)
                status=MagickFalse;
        }
    }
    edge_view=DestroyCacheView(edge_view);
    /*
       Free resources.
     */
    canny_cache=DestroyMatrixInfo(canny_cache);
    return(edge_image);
}

/*
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %                                                                             %
   %                                                                             %
   %                                                                             %
   %   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 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
    {
        PixelInfo
            direction[4]; /* horizontal, vertical, left and right diagonals */
    } ChannelStatistics;

    CacheView
    *image_view;

    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,
        r;

    size_t
        length;

    unsigned int
        number_grays;

    assert(image != (Image *) NULL);
    assert(image->signature == MagickCoreSignature);
    if (image->debug != MagickFalse)
        (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
    if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
        return((ChannelFeatures *) NULL);
    length=MaxPixelChannels+1UL;
    channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
                                                              sizeof(*channel_features));
    if (channel_features == (ChannelFeatures *) NULL)
        ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
    (void) memset(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) shared(status) \
    magick_number_threads(image,image,image->rows,1)
#endif
    for (r=0; r < (ssize_t) image->rows; r++)
    {
        register const Quantum
        *magick_restrict p;

        register ssize_t
            x;

        if (status == MagickFalse)
            continue;
        p=GetCacheViewVirtualPixels(image_view,0,r,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 != UndefinedPixelTrait)
                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) memset(&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;
    }
    /*
       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 != 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))
    {
        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) memset(&correlation,0,sizeof(correlation));
    (void) memset(density_x,0,2*(number_grays+1)*sizeof(*density_x));
    (void) memset(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
    (void) memset(density_y,0,2*(number_grays+1)*sizeof(*density_y));
    (void) memset(&mean,0,sizeof(mean));
    (void) memset(sum,0,number_grays*sizeof(*sum));
    (void) memset(&sum_squares,0,sizeof(sum_squares));
    (void) memset(density_xy,0,2*number_grays*sizeof(*density_xy));
    (void) memset(&entropy_x,0,sizeof(entropy_x));
    (void) memset(&entropy_xy,0,sizeof(entropy_xy));
    (void) memset(&entropy_xy1,0,sizeof(entropy_xy1));
    (void) memset(&entropy_xy2,0,sizeof(entropy_xy2));
    (void) memset(&entropy_y,0,sizeof(entropy_y));
    (void) memset(&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) memset(cooccurrence[i],0,number_grays*
                      sizeof(**cooccurrence));
        (void) memset(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 (r=0; r < (ssize_t) image->rows; r++)
    {
        register const Quantum
        *magick_restrict p;

        register ssize_t
            x;

        ssize_t
            offset,
            u,
            v;

        if (status == MagickFalse)
            continue;
        p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,r,image->columns+
                                    2*distance,distance+2,exception);
        if (p == (const Quantum *) NULL)
        {
            status=MagickFalse;
            continue;
        }
        p+=distance*GetPixelChannels(image);;
        for (x=0; x < (ssize_t) image->columns; x++)
        {
            for (i=0; i < 4; i++)
            {
                switch (i)
                {
                case 0:
                default:
                {
                    /*
                       Horizontal adjacency.
                     */
                    offset=(ssize_t) distance;
                    break;
                }
                case 1:
                {
                    /*
                       Vertical adjacency.
                     */
                    offset=(ssize_t) (image->columns+2*distance);
                    break;
                }
                case 2:
                {
                    /*
                       Right diagonal adjacency.
                     */
                    offset=(ssize_t) ((image->columns+2*distance)-distance);
                    break;
                }
                case 3:
                {
                    /*
                       Left diagonal adjacency.
                     */
                    offset=(ssize_t) ((image->columns+2*distance)+distance);
                    break;
                }
                }
                u=0;
                v=0;
                while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p)))
                    u++;
                while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image))))
                    v++;
                cooccurrence[u][v].direction[i].red++;
                cooccurrence[v][u].direction[i].red++;
                u=0;
                v=0;
                while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p)))
                    u++;
                while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image))))
                    v++;
                cooccurrence[u][v].direction[i].green++;
                cooccurrence[v][u].direction[i].green++;
                u=0;
                v=0;
                while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p)))
                    u++;
                while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image))))
                    v++;
                cooccurrence[u][v].direction[i].blue++;
                cooccurrence[v][u].direction[i].blue++;
                if (image->colorspace == CMYKColorspace)
                {
                    u=0;
                    v=0;
                    while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p)))
                        u++;
                    while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image))))
                        v++;
                    cooccurrence[u][v].direction[i].black++;
                    cooccurrence[v][u].direction[i].black++;
                }
                if (image->alpha_trait != UndefinedPixelTrait)
                {
                    u=0;
                    v=0;
                    while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p)))
                        u++;
                    while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image))))
                        v++;
                    cooccurrence[u][v].direction[i].alpha++;
                    cooccurrence[v][u].direction[i].alpha++;
                }
            }
            p+=GetPixelChannels(image);
        }
    }
    grays=(PixelPacket *) RelinquishMagickMemory(grays);
    image_view=DestroyCacheView(image_view);
    if (status == MagickFalse)
    {
        for (i=0; i < (ssize_t) number_grays; i++)
            cooccurrence[i]=(ChannelStatistics *)
                             RelinquishMagickMemory(cooccurrence[i]);
        cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
        channel_features=(ChannelFeatures *) RelinquishMagickMemory(
            channel_features);
        (void) ThrowMagickException(exception,GetMagickModule(),
                                    ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
        return(channel_features);
    }
    /*
       Normalize spatial dependence matrix.
     */
    for (i=0; i < 4; i++)
    {
        double
            normalize;

        register ssize_t
            y;

        switch (i)
        {
        case 0:
        default:
        {
            /*
               Horizontal adjacency.
             */
            normalize=2.0*image->rows*(image->columns-distance);
            break;
        }
        case 1:
        {
            /*
               Vertical adjacency.
             */
            normalize=2.0*(image->rows-distance)*image->columns;
            break;
        }
        case 2:
        {
            /*
               Right diagonal adjacency.
             */
            normalize=2.0*(image->rows-distance)*(image->columns-distance);
            break;
        }
        case 3:
        {
            /*
               Left diagonal adjacency.
             */
            normalize=2.0*(image->rows-distance)*(image->columns-distance);
            break;
        }
        }
        normalize=PerceptibleReciprocal(normalize);
        for (y=0; y < (ssize_t) number_grays; y++)
        {
            register ssize_t
                x;

            for (x=0; x < (ssize_t) number_grays; x++)
            {
                cooccurrence[x][y].direction[i].red*=normalize;
                cooccurrence[x][y].direction[i].green*=normalize;
                cooccurrence[x][y].direction[i].blue*=normalize;
                if (image->colorspace == CMYKColorspace)
                    cooccurrence[x][y].direction[i].black*=normalize;
                if (image->alpha_trait != UndefinedPixelTrait)
                    cooccurrence[x][y].direction[i].alpha*=normalize;
            }
        }
    }
    /*
       Compute texture features.
     */
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(image,image,number_grays,1)
#endif
    for (i=0; i < 4; i++)
    {
        register ssize_t
            y;

        for (y=0; y < (ssize_t) number_grays; y++)
        {
            register ssize_t
                x;

            for (x=0; x < (ssize_t) number_grays; x++)
            {
                /*
                   Angular second moment:  measure of homogeneity of the image.
                 */
                channel_features[RedPixelChannel].angular_second_moment[i]+=
                    cooccurrence[x][y].direction[i].red*
                    cooccurrence[x][y].direction[i].red;
                channel_features[GreenPixelChannel].angular_second_moment[i]+=
                    cooccurrence[x][y].direction[i].green*
                    cooccurrence[x][y].direction[i].green;
                channel_features[BluePixelChannel].angular_second_moment[i]+=
                    cooccurrence[x][y].direction[i].blue*
                    cooccurrence[x][y].direction[i].blue;
                if (image->colorspace == CMYKColorspace)
                    channel_features[BlackPixelChannel].angular_second_moment[i]+=
                        cooccurrence[x][y].direction[i].black*
                        cooccurrence[x][y].direction[i].black;
                if (image->alpha_trait != UndefinedPixelTrait)
                    channel_features[AlphaPixelChannel].angular_second_moment[i]+=
                        cooccurrence[x][y].direction[i].alpha*
                        cooccurrence[x][y].direction[i].alpha;
                /*
                   Correlation: measure of linear-dependencies in the image.
                 */
                sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
                sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
                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 != 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*
                                                 cooccurrence[x][y].direction[i].green;
                correlation.direction[i].blue+=x*y*
                                                cooccurrence[x][y].direction[i].blue;
                if (image->colorspace == CMYKColorspace)
                    correlation.direction[i].black+=x*y*
                                                     cooccurrence[x][y].direction[i].black;
                if (image->alpha_trait != UndefinedPixelTrait)
                    correlation.direction[i].alpha+=x*y*
                                                     cooccurrence[x][y].direction[i].alpha;
                /*
                   Inverse Difference Moment.
                 */
                channel_features[RedPixelChannel].inverse_difference_moment[i]+=
                    cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
                channel_features[GreenPixelChannel].inverse_difference_moment[i]+=
                    cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
                channel_features[BluePixelChannel].inverse_difference_moment[i]+=
                    cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
                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 != UndefinedPixelTrait)
                    channel_features[AlphaPixelChannel].inverse_difference_moment[i]+=
                        cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
                /*
                   Sum average.
                 */
                density_xy[y+x+2].direction[i].red+=
                    cooccurrence[x][y].direction[i].red;
                density_xy[y+x+2].direction[i].green+=
                    cooccurrence[x][y].direction[i].green;
                density_xy[y+x+2].direction[i].blue+=
                    cooccurrence[x][y].direction[i].blue;
                if (image->colorspace == CMYKColorspace)
                    density_xy[y+x+2].direction[i].black+=
                        cooccurrence[x][y].direction[i].black;
                if (image->alpha_trait != UndefinedPixelTrait)
                    density_xy[y+x+2].direction[i].alpha+=
                        cooccurrence[x][y].direction[i].alpha;
                /*
                   Entropy.
                 */
                channel_features[RedPixelChannel].entropy[i]-=
                    cooccurrence[x][y].direction[i].red*
                    MagickLog10(cooccurrence[x][y].direction[i].red);
                channel_features[GreenPixelChannel].entropy[i]-=
                    cooccurrence[x][y].direction[i].green*
                    MagickLog10(cooccurrence[x][y].direction[i].green);
                channel_features[BluePixelChannel].entropy[i]-=
                    cooccurrence[x][y].direction[i].blue*
                    MagickLog10(cooccurrence[x][y].direction[i].blue);
                if (image->colorspace == CMYKColorspace)
                    channel_features[BlackPixelChannel].entropy[i]-=
                        cooccurrence[x][y].direction[i].black*
                        MagickLog10(cooccurrence[x][y].direction[i].black);
                if (image->alpha_trait != UndefinedPixelTrait)
                    channel_features[AlphaPixelChannel].entropy[i]-=
                        cooccurrence[x][y].direction[i].alpha*
                        MagickLog10(cooccurrence[x][y].direction[i].alpha);
                /*
                   Information Measures of Correlation.
                 */
                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 != UndefinedPixelTrait)
                    density_x[x].direction[i].alpha+=
                        cooccurrence[x][y].direction[i].alpha;
                if (image->colorspace == CMYKColorspace)
                    density_x[x].direction[i].black+=
                        cooccurrence[x][y].direction[i].black;
                density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
                density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
                density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
                if (image->colorspace == CMYKColorspace)
                    density_y[y].direction[i].black+=
                        cooccurrence[x][y].direction[i].black;
                if (image->alpha_trait != UndefinedPixelTrait)
                    density_y[y].direction[i].alpha+=
                        cooccurrence[x][y].direction[i].alpha;
            }
            mean.direction[i].red+=y*sum[y].direction[i].red;
            sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
            mean.direction[i].green+=y*sum[y].direction[i].green;
            sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
            mean.direction[i].blue+=y*sum[y].direction[i].blue;
            sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
            if (image->colorspace == CMYKColorspace)
            {
                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 != UndefinedPixelTrait)
            {
                mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
                sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
            }
        }
        /*
           Correlation: measure of linear-dependencies in the image.
         */
        channel_features[RedPixelChannel].correlation[i]=
            (correlation.direction[i].red-mean.direction[i].red*
             mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
                                          (mean.direction[i].red*mean.direction[i].red))*sqrt(
                                         sum_squares.direction[i].red-(mean.direction[i].red*
                                                                       mean.direction[i].red)));
        channel_features[GreenPixelChannel].correlation[i]=
            (correlation.direction[i].green-mean.direction[i].green*
             mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
                                            (mean.direction[i].green*mean.direction[i].green))*sqrt(
                                           sum_squares.direction[i].green-(mean.direction[i].green*
                                                                           mean.direction[i].green)));
        channel_features[BluePixelChannel].correlation[i]=
            (correlation.direction[i].blue-mean.direction[i].blue*
             mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
                                           (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
                                          sum_squares.direction[i].blue-(mean.direction[i].blue*
                                                                         mean.direction[i].blue)));
        if (image->colorspace == CMYKColorspace)
            channel_features[BlackPixelChannel].correlation[i]=
                (correlation.direction[i].black-mean.direction[i].black*
                 mean.direction[i].black)/(sqrt(sum_squares.direction[i].black-
                                                (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 != UndefinedPixelTrait)
            channel_features[AlphaPixelChannel].correlation[i]=
                (correlation.direction[i].alpha-mean.direction[i].alpha*
                 mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
                                                (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt(
                                               sum_squares.direction[i].alpha-(mean.direction[i].alpha*
                                                                               mean.direction[i].alpha)));
    }
    /*
       Compute more texture features.
     */
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(image,image,number_grays,1)
#endif
    for (i=0; i < 4; i++)
    {
        register ssize_t
            x;

        for (x=2; x < (ssize_t) (2*number_grays); x++)
        {
            /*
               Sum average.
             */
            channel_features[RedPixelChannel].sum_average[i]+=
                x*density_xy[x].direction[i].red;
            channel_features[GreenPixelChannel].sum_average[i]+=
                x*density_xy[x].direction[i].green;
            channel_features[BluePixelChannel].sum_average[i]+=
                x*density_xy[x].direction[i].blue;
            if (image->colorspace == CMYKColorspace)
                channel_features[BlackPixelChannel].sum_average[i]+=
                    x*density_xy[x].direction[i].black;
            if (image->alpha_trait != UndefinedPixelTrait)
                channel_features[AlphaPixelChannel].sum_average[i]+=
                    x*density_xy[x].direction[i].alpha;
            /*
               Sum entropy.
             */
            channel_features[RedPixelChannel].sum_entropy[i]-=
                density_xy[x].direction[i].red*
                MagickLog10(density_xy[x].direction[i].red);
            channel_features[GreenPixelChannel].sum_entropy[i]-=
                density_xy[x].direction[i].green*
                MagickLog10(density_xy[x].direction[i].green);
            channel_features[BluePixelChannel].sum_entropy[i]-=
                density_xy[x].direction[i].blue*
                MagickLog10(density_xy[x].direction[i].blue);
            if (image->colorspace == CMYKColorspace)
                channel_features[BlackPixelChannel].sum_entropy[i]-=
                    density_xy[x].direction[i].black*
                    MagickLog10(density_xy[x].direction[i].black);
            if (image->alpha_trait != UndefinedPixelTrait)
                channel_features[AlphaPixelChannel].sum_entropy[i]-=
                    density_xy[x].direction[i].alpha*
                    MagickLog10(density_xy[x].direction[i].alpha);
            /*
               Sum variance.
             */
            channel_features[RedPixelChannel].sum_variance[i]+=
                (x-channel_features[RedPixelChannel].sum_entropy[i])*
                (x-channel_features[RedPixelChannel].sum_entropy[i])*
                density_xy[x].direction[i].red;
            channel_features[GreenPixelChannel].sum_variance[i]+=
                (x-channel_features[GreenPixelChannel].sum_entropy[i])*
                (x-channel_features[GreenPixelChannel].sum_entropy[i])*
                density_xy[x].direction[i].green;
            channel_features[BluePixelChannel].sum_variance[i]+=
                (x-channel_features[BluePixelChannel].sum_entropy[i])*
                (x-channel_features[BluePixelChannel].sum_entropy[i])*
                density_xy[x].direction[i].blue;
            if (image->colorspace == CMYKColorspace)
                channel_features[BlackPixelChannel].sum_variance[i]+=
                    (x-channel_features[BlackPixelChannel].sum_entropy[i])*
                    (x-channel_features[BlackPixelChannel].sum_entropy[i])*
                    density_xy[x].direction[i].black;
            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])*
                    density_xy[x].direction[i].alpha;
        }
    }
    /*
       Compute more texture features.
     */
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(image,image,number_grays,1)
#endif
    for (i=0; i < 4; i++)
    {
        register ssize_t
            y;

        for (y=0; y < (ssize_t) number_grays; y++)
        {
            register ssize_t
                x;

            for (x=0; x < (ssize_t) number_grays; x++)
            {
                /*
                   Sum of Squares: Variance
                 */
                variance.direction[i].red+=(y-mean.direction[i].red+1)*
                                            (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
                variance.direction[i].green+=(y-mean.direction[i].green+1)*
                                              (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
                variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
                                             (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
                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 != UndefinedPixelTrait)
                    variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
                                                  (y-mean.direction[i].alpha+1)*
                                                  cooccurrence[x][y].direction[i].alpha;
                /*
                   Sum average / Difference Variance.
                 */
                density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
                    cooccurrence[x][y].direction[i].red;
                density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
                    cooccurrence[x][y].direction[i].green;
                density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
                    cooccurrence[x][y].direction[i].blue;
                if (image->colorspace == CMYKColorspace)
                    density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
                        cooccurrence[x][y].direction[i].black;
                if (image->alpha_trait != UndefinedPixelTrait)
                    density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
                        cooccurrence[x][y].direction[i].alpha;
                /*
                   Information Measures of Correlation.
                 */
                entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
                                              MagickLog10(cooccurrence[x][y].direction[i].red);
                entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
                                                MagickLog10(cooccurrence[x][y].direction[i].green);
                entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
                                               MagickLog10(cooccurrence[x][y].direction[i].blue);
                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 != UndefinedPixelTrait)
                    entropy_xy.direction[i].alpha-=
                        cooccurrence[x][y].direction[i].alpha*MagickLog10(
                            cooccurrence[x][y].direction[i].alpha);
                entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
                                               MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red));
                entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
                                                 MagickLog10(density_x[x].direction[i].green*
                                                             density_y[y].direction[i].green));
                entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
                                                MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue));
                if (image->colorspace == CMYKColorspace)
                    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 != 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));
                entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
                                               density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red*
                                                                                         density_y[y].direction[i].red));
                entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
                                                 density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green*
                                                                                             density_y[y].direction[i].green));
                entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
                                                density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue*
                                                                                           density_y[y].direction[i].blue));
                if (image->colorspace == CMYKColorspace)
                    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 != 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));
            }
        }
        channel_features[RedPixelChannel].variance_sum_of_squares[i]=
            variance.direction[i].red;
        channel_features[GreenPixelChannel].variance_sum_of_squares[i]=
            variance.direction[i].green;
        channel_features[BluePixelChannel].variance_sum_of_squares[i]=
            variance.direction[i].blue;
        if (image->colorspace == CMYKColorspace)
            channel_features[BlackPixelChannel].variance_sum_of_squares[i]=
                variance.direction[i].black;
        if (image->alpha_trait != UndefinedPixelTrait)
            channel_features[AlphaPixelChannel].variance_sum_of_squares[i]=
                variance.direction[i].alpha;
    }
    /*
       Compute more texture features.
     */
    (void) memset(&variance,0,sizeof(variance));
    (void) memset(&sum_squares,0,sizeof(sum_squares));
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_threads(image,image,number_grays,1)
#endif
    for (i=0; i < 4; i++)
    {
        register ssize_t
            x;

        for (x=0; x < (ssize_t) number_grays; x++)
        {
            /*
               Difference variance.
             */
            variance.direction[i].red+=density_xy[x].direction[i].red;
            variance.direction[i].green+=density_xy[x].direction[i].green;
            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 != 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;
            sum_squares.direction[i].green+=density_xy[x].direction[i].green*
                                             density_xy[x].direction[i].green;
            sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
                                            density_xy[x].direction[i].blue;
            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 != UndefinedPixelTrait)
                sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
                                                 density_xy[x].direction[i].alpha;
            /*
               Difference entropy.
             */
            channel_features[RedPixelChannel].difference_entropy[i]-=
                density_xy[x].direction[i].red*
                MagickLog10(density_xy[x].direction[i].red);
            channel_features[GreenPixelChannel].difference_entropy[i]-=
                density_xy[x].direction[i].green*
                MagickLog10(density_xy[x].direction[i].green);
            channel_features[BluePixelChannel].difference_entropy[i]-=
                density_xy[x].direction[i].blue*
                MagickLog10(density_xy[x].direction[i].blue);
            if (image->colorspace == CMYKColorspace)
                channel_features[BlackPixelChannel].difference_entropy[i]-=
                    density_xy[x].direction[i].black*
                    MagickLog10(density_xy[x].direction[i].black);
            if (image->alpha_trait != UndefinedPixelTrait)
                channel_features[AlphaPixelChannel].difference_entropy[i]-=
                    density_xy[x].direction[i].alpha*
                    MagickLog10(density_xy[x].direction[i].alpha);
            /*
               Information Measures of Correlation.
             */
            entropy_x.direction[i].red-=(density_x[x].direction[i].red*
                                         MagickLog10(density_x[x].direction[i].red));
            entropy_x.direction[i].green-=(density_x[x].direction[i].green*
                                           MagickLog10(density_x[x].direction[i].green));
            entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
                                          MagickLog10(density_x[x].direction[i].blue));
            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 != 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*
                                         MagickLog10(density_y[x].direction[i].red));
            entropy_y.direction[i].green-=(density_y[x].direction[i].green*
                                           MagickLog10(density_y[x].direction[i].green));
            entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
                                          MagickLog10(density_y[x].direction[i].blue));
            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 != UndefinedPixelTrait)
                entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
                                               MagickLog10(density_y[x].direction[i].alpha));
        }
        /*
           Difference variance.
         */
        channel_features[RedPixelChannel].difference_variance[i]=
            (((double) number_grays*number_grays*sum_squares.direction[i].red)-
             (variance.direction[i].red*variance.direction[i].red))/
            ((double) number_grays*number_grays*number_grays*number_grays);
        channel_features[GreenPixelChannel].difference_variance[i]=
            (((double) number_grays*number_grays*sum_squares.direction[i].green)-
             (variance.direction[i].green*variance.direction[i].green))/
            ((double) number_grays*number_grays*number_grays*number_grays);
        channel_features[BluePixelChannel].difference_variance[i]=
            (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
             (variance.direction[i].blue*variance.direction[i].blue))/
            ((double) number_grays*number_grays*number_grays*number_grays);
        if (image->colorspace == CMYKColorspace)
            channel_features[BlackPixelChannel].difference_variance[i]=
                (((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 != 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))/
                ((double) number_grays*number_grays*number_grays*number_grays);
        /*
           Information Measures of Correlation.
         */
        channel_features[RedPixelChannel].measure_of_correlation_1[i]=
            (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
            (entropy_x.direction[i].red > entropy_y.direction[i].red ?
             entropy_x.direction[i].red : entropy_y.direction[i].red);
        channel_features[GreenPixelChannel].measure_of_correlation_1[i]=
            (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
            (entropy_x.direction[i].green > entropy_y.direction[i].green ?
             entropy_x.direction[i].green : entropy_y.direction[i].green);
        channel_features[BluePixelChannel].measure_of_correlation_1[i]=
            (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
            (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
             entropy_x.direction[i].blue : entropy_y.direction[i].blue);
        if (image->colorspace == CMYKColorspace)
            channel_features[BlackPixelChannel].measure_of_correlation_1[i]=
                (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 != 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*(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*(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*(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*(double) (entropy_xy2.direction[i].black-
                                                  entropy_xy.direction[i].black)))));
        if (image->alpha_trait != UndefinedPixelTrait)
            channel_features[AlphaPixelChannel].measure_of_correlation_2[i]=
                (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].alpha-
                                                  entropy_xy.direction[i].alpha)))));
    }
    /*
       Compute more texture features.
     */
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static) shared(status) \
    magick_number_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) memset(&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);
}

/*
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %                                                                             %
   %                                                                             %
   %                                                                             %
   %     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));
}

static Image *RenderHoughLines(const ImageInfo *image_info,const size_t columns,
                               const size_t rows,ExceptionInfo *exception)
{
#define BoundingBox  "viewbox"

    DrawInfo
    *draw_info;

    Image
    *image;

    MagickBooleanType
        status;

    /*
       Open image.
     */
    image=AcquireImage(image_info,exception);
    status=OpenBlob(image_info,image,ReadBinaryBlobMode,exception);
    if (status == MagickFalse)
    {
        image=DestroyImageList(image);
        return((Image *) NULL);
    }
    image->columns=columns;
    image->rows=rows;
    draw_info=CloneDrawInfo(image_info,(DrawInfo *) NULL);
    draw_info->affine.sx=image->resolution.x == 0.0 ? 1.0 : image->resolution.x/
                          DefaultResolution;
    draw_info->affine.sy=image->resolution.y == 0.0 ? 1.0 : image->resolution.y/
                          DefaultResolution;
    image->columns=(size_t) (draw_info->affine.sx*image->columns);
    image->rows=(size_t) (draw_info->affine.sy*image->rows);
    status=SetImageExtent(image,image->columns,image->rows,exception);
    if (status == MagickFalse)
        return(DestroyImageList(image));
    if (SetImageBackgroundColor(image,exception) == MagickFalse)
    {
        image=DestroyImageList(image);
        return((Image *) NULL);
    }
    /*
       Render drawing.
     */
    if (GetBlobStreamData(image) == (unsigned char *) NULL)
        draw_info->primitive=FileToString(image->filename,~0UL,exception);
    else
    {
        draw_info->primitive=(char *) AcquireMagickMemory((size_t)
                                                          GetBlobSize(image)+1);
        if (draw_info->primitive != (char *) NULL)
        {
            (void) memcpy(draw_info->primitive,GetBlobStreamData(image),
                          (size_t) GetBlobSize(image));
            draw_info->primitive[GetBlobSize(image)]='\0';
        }
    }
    (void) DrawImage(image,draw_info,exception);
    draw_info=DestroyDrawInfo(draw_info);
    (void) CloseBlob(image);
    return(GetFirstImageInList(image));
}

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
        *magick_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 atomic
#endif
            progress++;
            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;
    (void) FormatLocaleString(message,MagickPathExtent,
                              "# x1,y1 x2,y2 # count angle distance\n");
    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 %g %g\n",line.x1,line.y1,line.x2,line.y2,
                                          maxima,(double) x,(double) y);
                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,"%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=RenderHoughLines(image_info,image->columns,image->rows,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));
}

/*
   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   %                                                                             %
   %                                                                             %
   %                                                                             %
   %     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,0,0,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) shared(status,progress) \
    magick_number_threads(mean_image,mean_image,mean_image->rows,1)
#endif
    for (y=0; y < (ssize_t) mean_image->rows; y++)
    {
        register const Quantum
        *magick_restrict p;

        register Quantum
        *magick_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++)
        {
            PixelInfo
                mean_pixel,
                previous_pixel;

            PointInfo
                mean_location,
                previous_location;

            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++)
            {
                double
                    distance,
                    gamma;

                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++)
                {
                    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;
            }
            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);
        }
        if (SyncCacheViewAuthenticPixels(mean_view,exception) == MagickFalse)
            status=MagickFalse;
        if (image->progress_monitor != (MagickProgressMonitor) NULL)
        {
            MagickBooleanType
                proceed;

#if defined(MAGICKCORE_OPENMP_SUPPORT)
        #pragma omp atomic
#endif
            progress++;
            proceed=SetImageProgress(image,MeanShiftImageTag,progress,image->rows);
            if (proceed == MagickFalse)
                status=MagickFalse;
        }
    }
    mean_view=DestroyCacheView(mean_view);
    pixel_view=DestroyCacheView(pixel_view);
    image_view=DestroyCacheView(image_view);
    return(mean_image);
}