granicus.if.org Git - imagemagick/blob - MagickCore/feature.c

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %               FFFFF  EEEEE   AAA   TTTTT  U   U  RRRR   EEEEE               %
   7 %               F      E      A   A    T    U   U  R   R  E                   %
   8 %               FFF    EEE    AAAAA    T    U   U  RRRR   EEE                 %
   9 %               F      E      A   A    T    U   U  R R    E                   %
  10 %               F      EEEEE  A   A    T     UUU   R  R   EEEEE               %
  11 %                                                                             %
  12 %                                                                             %
  13 %                      MagickCore Image Feature Methods                       %
  14 %                                                                             %
  15 %                              Software Design                                %
  16 %                                   Cristy                                    %
  17 %                                 July 1992                                   %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "MagickCore/studio.h"
  44 #include "MagickCore/animate.h"
  45 #include "MagickCore/artifact.h"
  46 #include "MagickCore/blob.h"
  47 #include "MagickCore/blob-private.h"
  48 #include "MagickCore/cache.h"
  49 #include "MagickCore/cache-private.h"
  50 #include "MagickCore/cache-view.h"
  51 #include "MagickCore/client.h"
  52 #include "MagickCore/color.h"
  53 #include "MagickCore/color-private.h"
  54 #include "MagickCore/colorspace.h"
  55 #include "MagickCore/colorspace-private.h"
  56 #include "MagickCore/composite.h"
  57 #include "MagickCore/composite-private.h"
  58 #include "MagickCore/compress.h"
  59 #include "MagickCore/constitute.h"
  60 #include "MagickCore/display.h"
  61 #include "MagickCore/draw.h"
  62 #include "MagickCore/enhance.h"
  63 #include "MagickCore/exception.h"
  64 #include "MagickCore/exception-private.h"
  65 #include "MagickCore/feature.h"
  66 #include "MagickCore/gem.h"
  67 #include "MagickCore/geometry.h"
  68 #include "MagickCore/list.h"
  69 #include "MagickCore/image-private.h"
  70 #include "MagickCore/magic.h"
  71 #include "MagickCore/magick.h"
  72 #include "MagickCore/matrix.h"
  73 #include "MagickCore/memory_.h"
  74 #include "MagickCore/module.h"
  75 #include "MagickCore/monitor.h"
  76 #include "MagickCore/monitor-private.h"
  77 #include "MagickCore/morphology-private.h"
  78 #include "MagickCore/option.h"
  79 #include "MagickCore/paint.h"
  80 #include "MagickCore/pixel-accessor.h"
  81 #include "MagickCore/profile.h"
  82 #include "MagickCore/property.h"
  83 #include "MagickCore/quantize.h"
  84 #include "MagickCore/quantum-private.h"
  85 #include "MagickCore/random_.h"
  86 #include "MagickCore/resource_.h"
  87 #include "MagickCore/segment.h"
  88 #include "MagickCore/semaphore.h"
  89 #include "MagickCore/signature-private.h"
  90 #include "MagickCore/string_.h"
  91 #include "MagickCore/thread-private.h"
  92 #include "MagickCore/timer.h"
  93 #include "MagickCore/utility.h"
  94 #include "MagickCore/version.h"
  95 \f
  96 /*
  97 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  98 %                                                                             %
  99 %                                                                             %
 100 %                                                                             %
 101 %     C a n n y E d g e I m a g e                                             %
 102 %                                                                             %
 103 %                                                                             %
 104 %                                                                             %
 105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 106 %
 107 %  CannyEdgeImage() uses a multi-stage algorithm to detect a wide range of
 108 %  edges in images.
 109 %
 110 %  The format of the CannyEdgeImage method is:
 111 %
 112 %      Image *CannyEdgeImage(const Image *image,const double radius,
 113 %        const double sigma,const double lower_percent,
 114 %        const double upper_percent,ExceptionInfo *exception)
 115 %
 116 %  A description of each parameter follows:
 117 %
 118 %    o image: the image.
 119 %
 120 %    o radius: the radius of the gaussian smoothing filter.
 121 %
 122 %    o sigma: the sigma of the gaussian smoothing filter.
 123 %
 124 %    o lower_precent: percentage of edge pixels in the lower threshold.
 125 %
 126 %    o upper_percent: percentage of edge pixels in the upper threshold.
 127 %
 128 %    o exception: return any errors or warnings in this structure.
 129 %
 130 */
 131
 132 typedef struct _CannyInfo
 133 {
 134   double
 135     magnitude,
 136     intensity;
 137
 138   int
 139     orientation;
 140
 141   ssize_t
 142     x,
 143     y;
 144 } CannyInfo;
 145
 146 static inline MagickBooleanType IsAuthenticPixel(const Image *image,
 147   const ssize_t x,const ssize_t y)
 148 {
 149   if ((x < 0) || (x >= (ssize_t) image->columns))
 150     return(MagickFalse);
 151   if ((y < 0) || (y >= (ssize_t) image->rows))
 152     return(MagickFalse);
 153   return(MagickTrue);
 154 }
 155
 156 static MagickBooleanType TraceEdges(Image *edge_image,CacheView *edge_view,
 157   MatrixInfo *canny_cache,const ssize_t x,const ssize_t y,
 158   const double lower_threshold,ExceptionInfo *exception)
 159 {
 160   CannyInfo
 161     edge,
 162     pixel;
 163
 164   MagickBooleanType
 165     status;
 166
 167   register Quantum
 168     *q;
 169
 170   register ssize_t
 171     i;
 172
 173   q=GetCacheViewAuthenticPixels(edge_view,x,y,1,1,exception);
 174   if (q == (Quantum *) NULL)
 175     return(MagickFalse);
 176   *q=QuantumRange;
 177   status=SyncCacheViewAuthenticPixels(edge_view,exception);
 178   if (status == MagickFalse)
 179     return(MagickFalse);;
 180   if (GetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
 181     return(MagickFalse);
 182   edge.x=x;
 183   edge.y=y;
 184   if (SetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
 185     return(MagickFalse);
 186   for (i=1; i != 0; )
 187   {
 188     ssize_t
 189       v;
 190
 191     i--;
 192     status=GetMatrixElement(canny_cache,i,0,&edge);
 193     if (status == MagickFalse)
 194       return(MagickFalse);
 195     for (v=(-1); v <= 1; v++)
 196     {
 197       ssize_t
 198         u;
 199
 200       for (u=(-1); u <= 1; u++)
 201       {
 202         if ((u == 0) && (v == 0))
 203           continue;
 204         if (IsAuthenticPixel(edge_image,edge.x+u,edge.y+v) == MagickFalse)
 205           continue;
 206         /*
 207           Not an edge if gradient value is below the lower threshold.
 208         */
 209         q=GetCacheViewAuthenticPixels(edge_view,edge.x+u,edge.y+v,1,1,
 210           exception);
 211         if (q == (Quantum *) NULL)
 212           return(MagickFalse);
 213         status=GetMatrixElement(canny_cache,edge.x+u,edge.y+v,&pixel);
 214         if (status == MagickFalse)
 215           return(MagickFalse);
 216         if ((GetPixelIntensity(edge_image,q) == 0.0) &&
 217             (pixel.intensity >= lower_threshold))
 218           {
 219             *q=QuantumRange;
 220             status=SyncCacheViewAuthenticPixels(edge_view,exception);
 221             if (status == MagickFalse)
 222               return(MagickFalse);
 223             edge.x+=u;
 224             edge.y+=v;
 225             status=SetMatrixElement(canny_cache,i,0,&edge);
 226             if (status == MagickFalse)
 227               return(MagickFalse);
 228             i++;
 229           }
 230       }
 231     }
 232   }
 233   return(MagickTrue);
 234 }
 235
 236 MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
 237   const double sigma,const double lower_percent,const double upper_percent,
 238   ExceptionInfo *exception)
 239 {
 240 #define CannyEdgeImageTag  "CannyEdge/Image"
 241
 242   CacheView
 243     *edge_view;
 244
 245   CannyInfo
 246     pixel;
 247
 248   char
 249     geometry[MaxTextExtent];
 250
 251   double
 252     lower_threshold,
 253     max,
 254     min,
 255     upper_threshold;
 256
 257   Image
 258     *edge_image;
 259
 260   KernelInfo
 261     *kernel_info;
 262
 263   MagickBooleanType
 264     status;
 265
 266   MagickOffsetType
 267     progress;
 268
 269   MatrixInfo
 270     *canny_cache;
 271
 272   ssize_t
 273     y;
 274
 275   assert(image != (const Image *) NULL);
 276   assert(image->signature == MagickSignature);
 277   if (image->debug != MagickFalse)
 278     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 279   assert(exception != (ExceptionInfo *) NULL);
 280   assert(exception->signature == MagickSignature);
 281   /*
 282     Filter out noise.
 283   */
 284   (void) FormatLocaleString(geometry,MaxTextExtent,
 285     "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
 286   kernel_info=AcquireKernelInfo(geometry);
 287   if (kernel_info == (KernelInfo *) NULL)
 288     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 289   edge_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info,
 290     UndefinedCompositeOp,0.0,exception);
 291   kernel_info=DestroyKernelInfo(kernel_info);
 292   if (edge_image == (Image *) NULL)
 293     return((Image *) NULL);
 294   if (SetImageColorspace(edge_image,GRAYColorspace,exception) == MagickFalse)
 295     {
 296       edge_image=DestroyImage(edge_image);
 297       return((Image *) NULL);
 298     }
 299   /*
 300     Find the intensity gradient of the image.
 301   */
 302   canny_cache=AcquireMatrixInfo(edge_image->columns,edge_image->rows,
 303     sizeof(CannyInfo),exception);
 304   if (canny_cache == (MatrixInfo *) NULL)
 305     {
 306       edge_image=DestroyImage(edge_image);
 307       return((Image *) NULL);
 308     }
 309   status=MagickTrue;
 310   edge_view=AcquireVirtualCacheView(edge_image,exception);
 311 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 312   #pragma omp parallel for schedule(static,4) shared(status) \
 313     magick_threads(edge_image,edge_image,edge_image->rows,1)
 314 #endif
 315   for (y=0; y < (ssize_t) edge_image->rows; y++)
 316   {
 317     register const Quantum
 318       *restrict p;
 319
 320     register ssize_t
 321       x;
 322
 323     if (status == MagickFalse)
 324       continue;
 325     p=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns+1,2,
 326       exception);
 327     if (p == (const Quantum *) NULL)
 328       {
 329         status=MagickFalse;
 330         continue;
 331       }
 332     for (x=0; x < (ssize_t) edge_image->columns; x++)
 333     {
 334       CannyInfo
 335         pixel;
 336
 337       double
 338         dx,
 339         dy;
 340
 341       register const Quantum
 342         *restrict kernel_pixels;
 343
 344       ssize_t
 345         v;
 346
 347       static double
 348         Gx[2][2] =
 349         {
 350           { -1.0,  +1.0 },
 351           { -1.0,  +1.0 }
 352         },
 353         Gy[2][2] =
 354         {
 355           { +1.0, +1.0 },
 356           { -1.0, -1.0 }
 357         };
 358
 359       (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
 360       dx=0.0;
 361       dy=0.0;
 362       kernel_pixels=p;
 363       for (v=0; v < 2; v++)
 364       {
 365         ssize_t
 366           u;
 367
 368         for (u=0; u < 2; u++)
 369         {
 370           double
 371             intensity;
 372
 373           intensity=GetPixelIntensity(edge_image,kernel_pixels+u);
 374           dx+=0.5*Gx[v][u]*intensity;
 375           dy+=0.5*Gy[v][u]*intensity;
 376         }
 377         kernel_pixels+=edge_image->columns+1;
 378       }
 379       pixel.magnitude=hypot(dx,dy);
 380       pixel.orientation=0;
 381       if (fabs(dx) > MagickEpsilon)
 382         {
 383           double
 384             slope;
 385
 386           slope=dy/dx;
 387           if (slope < 0.0)
 388             {
 389               if (slope < -2.41421356237)
 390                 pixel.orientation=0;
 391               else
 392                 if (slope < -0.414213562373)
 393                   pixel.orientation=1;
 394                 else
 395                   pixel.orientation=2;
 396             }
 397           else
 398             {
 399               if (slope > 2.41421356237)
 400                 pixel.orientation=0;
 401               else
 402                 if (slope > 0.414213562373)
 403                   pixel.orientation=3;
 404                 else
 405                   pixel.orientation=2;
 406             }
 407         }
 408       if (SetMatrixElement(canny_cache,x,y,&pixel) == MagickFalse)
 409         continue;
 410       p+=GetPixelChannels(edge_image);
 411     }
 412   }
 413   edge_view=DestroyCacheView(edge_view);
 414   /*
 415     Non-maxima suppression, remove pixels that are not considered to be part
 416     of an edge.
 417   */
 418   progress=0;
 419   (void) GetMatrixElement(canny_cache,0,0,&pixel);
 420   max=pixel.intensity;
 421   min=pixel.intensity;
 422   edge_view=AcquireAuthenticCacheView(edge_image,exception);
 423 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 424   #pragma omp parallel for schedule(static,4) shared(status) \
 425     magick_threads(edge_image,edge_image,edge_image->rows,1)
 426 #endif
 427   for (y=0; y < (ssize_t) edge_image->rows; y++)
 428   {
 429     register Quantum
 430       *restrict q;
 431
 432     register ssize_t
 433       x;
 434
 435     if (status == MagickFalse)
 436       continue;
 437     q=GetCacheViewAuthenticPixels(edge_view,0,y,edge_image->columns,1,
 438       exception);
 439     if (q == (Quantum *) NULL)
 440       {
 441         status=MagickFalse;
 442         continue;
 443       }
 444     for (x=0; x < (ssize_t) edge_image->columns; x++)
 445     {
 446       CannyInfo
 447         alpha_pixel,
 448         beta_pixel,
 449         pixel;
 450
 451       (void) GetMatrixElement(canny_cache,x,y,&pixel);
 452       switch (pixel.orientation)
 453       {
 454         case 0:
 455         default:
 456         {
 457           /*
 458             0 degrees, north and south.
 459           */
 460           (void) GetMatrixElement(canny_cache,x,y-1,&alpha_pixel);
 461           (void) GetMatrixElement(canny_cache,x,y+1,&beta_pixel);
 462           break;
 463         }
 464         case 1:
 465         {
 466           /*
 467             45 degrees, northwest and southeast.
 468           */
 469           (void) GetMatrixElement(canny_cache,x-1,y-1,&alpha_pixel);
 470           (void) GetMatrixElement(canny_cache,x+1,y+1,&beta_pixel);
 471           break;
 472         }
 473         case 2:
 474         {
 475           /*
 476             90 degrees, east and west.
 477           */
 478           (void) GetMatrixElement(canny_cache,x-1,y,&alpha_pixel);
 479           (void) GetMatrixElement(canny_cache,x+1,y,&beta_pixel);
 480           break;
 481         }
 482         case 3:
 483         {
 484           /*
 485             135 degrees, northeast and southwest.
 486           */
 487           (void) GetMatrixElement(canny_cache,x+1,y-1,&beta_pixel);
 488           (void) GetMatrixElement(canny_cache,x-1,y+1,&alpha_pixel);
 489           break;
 490         }
 491       }
 492       pixel.intensity=pixel.magnitude;
 493       if ((pixel.magnitude < alpha_pixel.magnitude) ||
 494           (pixel.magnitude < beta_pixel.magnitude))
 495         pixel.intensity=0;
 496       (void) SetMatrixElement(canny_cache,x,y,&pixel);
 497 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 498       #pragma omp critical (MagickCore_CannyEdgeImage)
 499 #endif
 500       {
 501         if (pixel.intensity < min)
 502           min=pixel.intensity;
 503         if (pixel.intensity > max)
 504           max=pixel.intensity;
 505       }
 506       *q=0;
 507       q+=GetPixelChannels(edge_image);
 508     }
 509     if (SyncCacheViewAuthenticPixels(edge_view,exception) == MagickFalse)
 510       status=MagickFalse;
 511   }
 512   edge_view=DestroyCacheView(edge_view);
 513   /*
 514     Estimate hysteresis threshold.
 515   */
 516   lower_threshold=lower_percent*(max-min)+min;
 517   upper_threshold=upper_percent*(max-min)+min;
 518   /*
 519     Hysteresis threshold.
 520   */
 521   edge_view=AcquireAuthenticCacheView(edge_image,exception);
 522   for (y=0; y < (ssize_t) edge_image->rows; y++)
 523   {
 524     register ssize_t
 525       x;
 526
 527     if (status == MagickFalse)
 528       continue;
 529     for (x=0; x < (ssize_t) edge_image->columns; x++)
 530     {
 531       CannyInfo
 532         pixel;
 533
 534       register const Quantum
 535         *restrict p;
 536
 537       /*
 538         Edge if pixel gradient higher than upper threshold.
 539       */
 540       p=GetCacheViewVirtualPixels(edge_view,x,y,1,1,exception);
 541       if (p == (const Quantum *) NULL)
 542         continue;
 543       status=GetMatrixElement(canny_cache,x,y,&pixel);
 544       if (status == MagickFalse)
 545         continue;
 546       if ((GetPixelIntensity(edge_image,p) == 0.0) &&
 547           (pixel.intensity >= upper_threshold))
 548         status=TraceEdges(edge_image,edge_view,canny_cache,x,y,lower_threshold,
 549           exception);
 550     }
 551     if (image->progress_monitor != (MagickProgressMonitor) NULL)
 552       {
 553         MagickBooleanType
 554           proceed;
 555
 556 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 557         #pragma omp critical (MagickCore_CannyEdgeImage)
 558 #endif
 559         proceed=SetImageProgress(image,CannyEdgeImageTag,progress++,
 560           image->rows);
 561         if (proceed == MagickFalse)
 562           status=MagickFalse;
 563       }
 564   }
 565   edge_view=DestroyCacheView(edge_view);
 566   /*
 567     Free resources.
 568   */
 569   canny_cache=DestroyMatrixInfo(canny_cache);
 570   return(edge_image);
 571 }
 572 \f
 573 /*
 574 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 575 %                                                                             %
 576 %                                                                             %
 577 %                                                                             %
 578 %   G e t I m a g e F e a t u r e s                                           %
 579 %                                                                             %
 580 %                                                                             %
 581 %                                                                             %
 582 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 583 %
 584 %  GetImageFeatures() returns features for each channel in the image in
 585 %  each of four directions (horizontal, vertical, left and right diagonals)
 586 %  for the specified distance.  The features include the angular second
 587 %  moment, contrast, correlation, sum of squares: variance, inverse difference
 588 %  moment, sum average, sum varience, sum entropy, entropy, difference variance,%  difference entropy, information measures of correlation 1, information
 589 %  measures of correlation 2, and maximum correlation coefficient.  You can
 590 %  access the red channel contrast, for example, like this:
 591 %
 592 %      channel_features=GetImageFeatures(image,1,exception);
 593 %      contrast=channel_features[RedPixelChannel].contrast[0];
 594 %
 595 %  Use MagickRelinquishMemory() to free the features buffer.
 596 %
 597 %  The format of the GetImageFeatures method is:
 598 %
 599 %      ChannelFeatures *GetImageFeatures(const Image *image,
 600 %        const size_t distance,ExceptionInfo *exception)
 601 %
 602 %  A description of each parameter follows:
 603 %
 604 %    o image: the image.
 605 %
 606 %    o distance: the distance.
 607 %
 608 %    o exception: return any errors or warnings in this structure.
 609 %
 610 */
 611
 612 static inline ssize_t MagickAbsoluteValue(const ssize_t x)
 613 {
 614   if (x < 0)
 615     return(-x);
 616   return(x);
 617 }
 618
 619 static inline double MagickLog10(const double x)
 620 {
 621 #define Log10Epsilon  (1.0e-11)
 622
 623  if (fabs(x) < Log10Epsilon)
 624    return(log10(Log10Epsilon));
 625  return(log10(fabs(x)));
 626 }
 627
 628 MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
 629   const size_t distance,ExceptionInfo *exception)
 630 {
 631   typedef struct _ChannelStatistics
 632   {
 633     PixelInfo
 634       direction[4];  /* horizontal, vertical, left and right diagonals */
 635   } ChannelStatistics;
 636
 637   CacheView
 638     *image_view;
 639
 640   ChannelFeatures
 641     *channel_features;
 642
 643   ChannelStatistics
 644     **cooccurrence,
 645     correlation,
 646     *density_x,
 647     *density_xy,
 648     *density_y,
 649     entropy_x,
 650     entropy_xy,
 651     entropy_xy1,
 652     entropy_xy2,
 653     entropy_y,
 654     mean,
 655     **Q,
 656     *sum,
 657     sum_squares,
 658     variance;
 659
 660   PixelPacket
 661     gray,
 662     *grays;
 663
 664   MagickBooleanType
 665     status;
 666
 667   register ssize_t
 668     i;
 669
 670   size_t
 671     length;
 672
 673   ssize_t
 674     y;
 675
 676   unsigned int
 677     number_grays;
 678
 679   assert(image != (Image *) NULL);
 680   assert(image->signature == MagickSignature);
 681   if (image->debug != MagickFalse)
 682     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 683   if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
 684     return((ChannelFeatures *) NULL);
 685   length=CompositeChannels+1UL;
 686   channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
 687     sizeof(*channel_features));
 688   if (channel_features == (ChannelFeatures *) NULL)
 689     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 690   (void) ResetMagickMemory(channel_features,0,length*
 691     sizeof(*channel_features));
 692   /*
 693     Form grays.
 694   */
 695   grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays));
 696   if (grays == (PixelPacket *) NULL)
 697     {
 698       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 699         channel_features);
 700       (void) ThrowMagickException(exception,GetMagickModule(),
 701         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 702       return(channel_features);
 703     }
 704   for (i=0; i <= (ssize_t) MaxMap; i++)
 705   {
 706     grays[i].red=(~0U);
 707     grays[i].green=(~0U);
 708     grays[i].blue=(~0U);
 709     grays[i].alpha=(~0U);
 710     grays[i].black=(~0U);
 711   }
 712   status=MagickTrue;
 713   image_view=AcquireVirtualCacheView(image,exception);
 714 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 715   #pragma omp parallel for schedule(static,4) shared(status) \
 716     magick_threads(image,image,image->rows,1)
 717 #endif
 718   for (y=0; y < (ssize_t) image->rows; y++)
 719   {
 720     register const Quantum
 721       *restrict p;
 722
 723     register ssize_t
 724       x;
 725
 726     if (status == MagickFalse)
 727       continue;
 728     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 729     if (p == (const Quantum *) NULL)
 730       {
 731         status=MagickFalse;
 732         continue;
 733       }
 734     for (x=0; x < (ssize_t) image->columns; x++)
 735     {
 736       grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
 737         ScaleQuantumToMap(GetPixelRed(image,p));
 738       grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
 739         ScaleQuantumToMap(GetPixelGreen(image,p));
 740       grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
 741         ScaleQuantumToMap(GetPixelBlue(image,p));
 742       if (image->colorspace == CMYKColorspace)
 743         grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
 744           ScaleQuantumToMap(GetPixelBlack(image,p));
 745       if (image->alpha_trait == BlendPixelTrait)
 746         grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
 747           ScaleQuantumToMap(GetPixelAlpha(image,p));
 748       p+=GetPixelChannels(image);
 749     }
 750   }
 751   image_view=DestroyCacheView(image_view);
 752   if (status == MagickFalse)
 753     {
 754       grays=(PixelPacket *) RelinquishMagickMemory(grays);
 755       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 756         channel_features);
 757       return(channel_features);
 758     }
 759   (void) ResetMagickMemory(&gray,0,sizeof(gray));
 760   for (i=0; i <= (ssize_t) MaxMap; i++)
 761   {
 762     if (grays[i].red != ~0U)
 763       grays[gray.red++].red=grays[i].red;
 764     if (grays[i].green != ~0U)
 765       grays[gray.green++].green=grays[i].green;
 766     if (grays[i].blue != ~0U)
 767       grays[gray.blue++].blue=grays[i].blue;
 768     if (image->colorspace == CMYKColorspace)
 769       if (grays[i].black != ~0U)
 770         grays[gray.black++].black=grays[i].black;
 771     if (image->alpha_trait == BlendPixelTrait)
 772       if (grays[i].alpha != ~0U)
 773         grays[gray.alpha++].alpha=grays[i].alpha;
 774   }
 775   /*
 776     Allocate spatial dependence matrix.
 777   */
 778   number_grays=gray.red;
 779   if (gray.green > number_grays)
 780     number_grays=gray.green;
 781   if (gray.blue > number_grays)
 782     number_grays=gray.blue;
 783   if (image->colorspace == CMYKColorspace)
 784     if (gray.black > number_grays)
 785       number_grays=gray.black;
 786   if (image->alpha_trait == BlendPixelTrait)
 787     if (gray.alpha > number_grays)
 788       number_grays=gray.alpha;
 789   cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays,
 790     sizeof(*cooccurrence));
 791   density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
 792     sizeof(*density_x));
 793   density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
 794     sizeof(*density_xy));
 795   density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
 796     sizeof(*density_y));
 797   Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q));
 798   sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum));
 799   if ((cooccurrence == (ChannelStatistics **) NULL) ||
 800       (density_x == (ChannelStatistics *) NULL) ||
 801       (density_xy == (ChannelStatistics *) NULL) ||
 802       (density_y == (ChannelStatistics *) NULL) ||
 803       (Q == (ChannelStatistics **) NULL) ||
 804       (sum == (ChannelStatistics *) NULL))
 805     {
 806       if (Q != (ChannelStatistics **) NULL)
 807         {
 808           for (i=0; i < (ssize_t) number_grays; i++)
 809             Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
 810           Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
 811         }
 812       if (sum != (ChannelStatistics *) NULL)
 813         sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
 814       if (density_y != (ChannelStatistics *) NULL)
 815         density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
 816       if (density_xy != (ChannelStatistics *) NULL)
 817         density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
 818       if (density_x != (ChannelStatistics *) NULL)
 819         density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
 820       if (cooccurrence != (ChannelStatistics **) NULL)
 821         {
 822           for (i=0; i < (ssize_t) number_grays; i++)
 823             cooccurrence[i]=(ChannelStatistics *)
 824               RelinquishMagickMemory(cooccurrence[i]);
 825           cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
 826             cooccurrence);
 827         }
 828       grays=(PixelPacket *) RelinquishMagickMemory(grays);
 829       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 830         channel_features);
 831       (void) ThrowMagickException(exception,GetMagickModule(),
 832         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 833       return(channel_features);
 834     }
 835   (void) ResetMagickMemory(&correlation,0,sizeof(correlation));
 836   (void) ResetMagickMemory(density_x,0,2*(number_grays+1)*sizeof(*density_x));
 837   (void) ResetMagickMemory(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
 838   (void) ResetMagickMemory(density_y,0,2*(number_grays+1)*sizeof(*density_y));
 839   (void) ResetMagickMemory(&mean,0,sizeof(mean));
 840   (void) ResetMagickMemory(sum,0,number_grays*sizeof(*sum));
 841   (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
 842   (void) ResetMagickMemory(density_xy,0,2*number_grays*sizeof(*density_xy));
 843   (void) ResetMagickMemory(&entropy_x,0,sizeof(entropy_x));
 844   (void) ResetMagickMemory(&entropy_xy,0,sizeof(entropy_xy));
 845   (void) ResetMagickMemory(&entropy_xy1,0,sizeof(entropy_xy1));
 846   (void) ResetMagickMemory(&entropy_xy2,0,sizeof(entropy_xy2));
 847   (void) ResetMagickMemory(&entropy_y,0,sizeof(entropy_y));
 848   (void) ResetMagickMemory(&variance,0,sizeof(variance));
 849   for (i=0; i < (ssize_t) number_grays; i++)
 850   {
 851     cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
 852       sizeof(**cooccurrence));
 853     Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
 854     if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
 855         (Q[i] == (ChannelStatistics *) NULL))
 856       break;
 857     (void) ResetMagickMemory(cooccurrence[i],0,number_grays*
 858       sizeof(**cooccurrence));
 859     (void) ResetMagickMemory(Q[i],0,number_grays*sizeof(**Q));
 860   }
 861   if (i < (ssize_t) number_grays)
 862     {
 863       for (i--; i >= 0; i--)
 864       {
 865         if (Q[i] != (ChannelStatistics *) NULL)
 866           Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
 867         if (cooccurrence[i] != (ChannelStatistics *) NULL)
 868           cooccurrence[i]=(ChannelStatistics *)
 869             RelinquishMagickMemory(cooccurrence[i]);
 870       }
 871       Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
 872       cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
 873       sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
 874       density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
 875       density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
 876       density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
 877       grays=(PixelPacket *) RelinquishMagickMemory(grays);
 878       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 879         channel_features);
 880       (void) ThrowMagickException(exception,GetMagickModule(),
 881         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 882       return(channel_features);
 883     }
 884   /*
 885     Initialize spatial dependence matrix.
 886   */
 887   status=MagickTrue;
 888   image_view=AcquireVirtualCacheView(image,exception);
 889   for (y=0; y < (ssize_t) image->rows; y++)
 890   {
 891     register const Quantum
 892       *restrict p;
 893
 894     register ssize_t
 895       x;
 896
 897     ssize_t
 898       i,
 899       offset,
 900       u,
 901       v;
 902
 903     if (status == MagickFalse)
 904       continue;
 905     p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,y,image->columns+
 906       2*distance,distance+2,exception);
 907     if (p == (const Quantum *) NULL)
 908       {
 909         status=MagickFalse;
 910         continue;
 911       }
 912     p+=distance*GetPixelChannels(image);;
 913     for (x=0; x < (ssize_t) image->columns; x++)
 914     {
 915       for (i=0; i < 4; i++)
 916       {
 917         switch (i)
 918         {
 919           case 0:
 920           default:
 921           {
 922             /*
 923               Horizontal adjacency.
 924             */
 925             offset=(ssize_t) distance;
 926             break;
 927           }
 928           case 1:
 929           {
 930             /*
 931               Vertical adjacency.
 932             */
 933             offset=(ssize_t) (image->columns+2*distance);
 934             break;
 935           }
 936           case 2:
 937           {
 938             /*
 939               Right diagonal adjacency.
 940             */
 941             offset=(ssize_t) ((image->columns+2*distance)-distance);
 942             break;
 943           }
 944           case 3:
 945           {
 946             /*
 947               Left diagonal adjacency.
 948             */
 949             offset=(ssize_t) ((image->columns+2*distance)+distance);
 950             break;
 951           }
 952         }
 953         u=0;
 954         v=0;
 955         while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p)))
 956           u++;
 957         while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image))))
 958           v++;
 959         cooccurrence[u][v].direction[i].red++;
 960         cooccurrence[v][u].direction[i].red++;
 961         u=0;
 962         v=0;
 963         while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p)))
 964           u++;
 965         while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image))))
 966           v++;
 967         cooccurrence[u][v].direction[i].green++;
 968         cooccurrence[v][u].direction[i].green++;
 969         u=0;
 970         v=0;
 971         while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p)))
 972           u++;
 973         while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image))))
 974           v++;
 975         cooccurrence[u][v].direction[i].blue++;
 976         cooccurrence[v][u].direction[i].blue++;
 977         if (image->colorspace == CMYKColorspace)
 978           {
 979             u=0;
 980             v=0;
 981             while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p)))
 982               u++;
 983             while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image))))
 984               v++;
 985             cooccurrence[u][v].direction[i].black++;
 986             cooccurrence[v][u].direction[i].black++;
 987           }
 988         if (image->alpha_trait == BlendPixelTrait)
 989           {
 990             u=0;
 991             v=0;
 992             while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p)))
 993               u++;
 994             while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image))))
 995               v++;
 996             cooccurrence[u][v].direction[i].alpha++;
 997             cooccurrence[v][u].direction[i].alpha++;
 998           }
 999       }
1000       p+=GetPixelChannels(image);
1001     }
1002   }
1003   grays=(PixelPacket *) RelinquishMagickMemory(grays);
1004   image_view=DestroyCacheView(image_view);
1005   if (status == MagickFalse)
1006     {
1007       for (i=0; i < (ssize_t) number_grays; i++)
1008         cooccurrence[i]=(ChannelStatistics *)
1009           RelinquishMagickMemory(cooccurrence[i]);
1010       cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1011       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1012         channel_features);
1013       (void) ThrowMagickException(exception,GetMagickModule(),
1014         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1015       return(channel_features);
1016     }
1017   /*
1018     Normalize spatial dependence matrix.
1019   */
1020   for (i=0; i < 4; i++)
1021   {
1022     double
1023       normalize;
1024
1025     register ssize_t
1026       y;
1027
1028     switch (i)
1029     {
1030       case 0:
1031       default:
1032       {
1033         /*
1034           Horizontal adjacency.
1035         */
1036         normalize=2.0*image->rows*(image->columns-distance);
1037         break;
1038       }
1039       case 1:
1040       {
1041         /*
1042           Vertical adjacency.
1043         */
1044         normalize=2.0*(image->rows-distance)*image->columns;
1045         break;
1046       }
1047       case 2:
1048       {
1049         /*
1050           Right diagonal adjacency.
1051         */
1052         normalize=2.0*(image->rows-distance)*(image->columns-distance);
1053         break;
1054       }
1055       case 3:
1056       {
1057         /*
1058           Left diagonal adjacency.
1059         */
1060         normalize=2.0*(image->rows-distance)*(image->columns-distance);
1061         break;
1062       }
1063     }
1064     normalize=PerceptibleReciprocal(normalize);
1065     for (y=0; y < (ssize_t) number_grays; y++)
1066     {
1067       register ssize_t
1068         x;
1069
1070       for (x=0; x < (ssize_t) number_grays; x++)
1071       {
1072         cooccurrence[x][y].direction[i].red*=normalize;
1073         cooccurrence[x][y].direction[i].green*=normalize;
1074         cooccurrence[x][y].direction[i].blue*=normalize;
1075         if (image->colorspace == CMYKColorspace)
1076           cooccurrence[x][y].direction[i].black*=normalize;
1077         if (image->alpha_trait == BlendPixelTrait)
1078           cooccurrence[x][y].direction[i].alpha*=normalize;
1079       }
1080     }
1081   }
1082   /*
1083     Compute texture features.
1084   */
1085 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1086   #pragma omp parallel for schedule(static,4) shared(status) \
1087     magick_threads(image,image,number_grays,1)
1088 #endif
1089   for (i=0; i < 4; i++)
1090   {
1091     register ssize_t
1092       y;
1093
1094     for (y=0; y < (ssize_t) number_grays; y++)
1095     {
1096       register ssize_t
1097         x;
1098
1099       for (x=0; x < (ssize_t) number_grays; x++)
1100       {
1101         /*
1102           Angular second moment:  measure of homogeneity of the image.
1103         */
1104         channel_features[RedPixelChannel].angular_second_moment[i]+=
1105           cooccurrence[x][y].direction[i].red*
1106           cooccurrence[x][y].direction[i].red;
1107         channel_features[GreenPixelChannel].angular_second_moment[i]+=
1108           cooccurrence[x][y].direction[i].green*
1109           cooccurrence[x][y].direction[i].green;
1110         channel_features[BluePixelChannel].angular_second_moment[i]+=
1111           cooccurrence[x][y].direction[i].blue*
1112           cooccurrence[x][y].direction[i].blue;
1113         if (image->colorspace == CMYKColorspace)
1114           channel_features[BlackPixelChannel].angular_second_moment[i]+=
1115             cooccurrence[x][y].direction[i].black*
1116             cooccurrence[x][y].direction[i].black;
1117         if (image->alpha_trait == BlendPixelTrait)
1118           channel_features[AlphaPixelChannel].angular_second_moment[i]+=
1119             cooccurrence[x][y].direction[i].alpha*
1120             cooccurrence[x][y].direction[i].alpha;
1121         /*
1122           Correlation: measure of linear-dependencies in the image.
1123         */
1124         sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1125         sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1126         sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1127         if (image->colorspace == CMYKColorspace)
1128           sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black;
1129         if (image->alpha_trait == BlendPixelTrait)
1130           sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha;
1131         correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
1132         correlation.direction[i].green+=x*y*
1133           cooccurrence[x][y].direction[i].green;
1134         correlation.direction[i].blue+=x*y*
1135           cooccurrence[x][y].direction[i].blue;
1136         if (image->colorspace == CMYKColorspace)
1137           correlation.direction[i].black+=x*y*
1138             cooccurrence[x][y].direction[i].black;
1139         if (image->alpha_trait == BlendPixelTrait)
1140           correlation.direction[i].alpha+=x*y*
1141             cooccurrence[x][y].direction[i].alpha;
1142         /*
1143           Inverse Difference Moment.
1144         */
1145         channel_features[RedPixelChannel].inverse_difference_moment[i]+=
1146           cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
1147         channel_features[GreenPixelChannel].inverse_difference_moment[i]+=
1148           cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
1149         channel_features[BluePixelChannel].inverse_difference_moment[i]+=
1150           cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
1151         if (image->colorspace == CMYKColorspace)
1152           channel_features[BlackPixelChannel].inverse_difference_moment[i]+=
1153             cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1);
1154         if (image->alpha_trait == BlendPixelTrait)
1155           channel_features[AlphaPixelChannel].inverse_difference_moment[i]+=
1156             cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
1157         /*
1158           Sum average.
1159         */
1160         density_xy[y+x+2].direction[i].red+=
1161           cooccurrence[x][y].direction[i].red;
1162         density_xy[y+x+2].direction[i].green+=
1163           cooccurrence[x][y].direction[i].green;
1164         density_xy[y+x+2].direction[i].blue+=
1165           cooccurrence[x][y].direction[i].blue;
1166         if (image->colorspace == CMYKColorspace)
1167           density_xy[y+x+2].direction[i].black+=
1168             cooccurrence[x][y].direction[i].black;
1169         if (image->alpha_trait == BlendPixelTrait)
1170           density_xy[y+x+2].direction[i].alpha+=
1171             cooccurrence[x][y].direction[i].alpha;
1172         /*
1173           Entropy.
1174         */
1175         channel_features[RedPixelChannel].entropy[i]-=
1176           cooccurrence[x][y].direction[i].red*
1177           MagickLog10(cooccurrence[x][y].direction[i].red);
1178         channel_features[GreenPixelChannel].entropy[i]-=
1179           cooccurrence[x][y].direction[i].green*
1180           MagickLog10(cooccurrence[x][y].direction[i].green);
1181         channel_features[BluePixelChannel].entropy[i]-=
1182           cooccurrence[x][y].direction[i].blue*
1183           MagickLog10(cooccurrence[x][y].direction[i].blue);
1184         if (image->colorspace == CMYKColorspace)
1185           channel_features[BlackPixelChannel].entropy[i]-=
1186             cooccurrence[x][y].direction[i].black*
1187             MagickLog10(cooccurrence[x][y].direction[i].black);
1188         if (image->alpha_trait == BlendPixelTrait)
1189           channel_features[AlphaPixelChannel].entropy[i]-=
1190             cooccurrence[x][y].direction[i].alpha*
1191             MagickLog10(cooccurrence[x][y].direction[i].alpha);
1192         /*
1193           Information Measures of Correlation.
1194         */
1195         density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
1196         density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
1197         density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1198         if (image->alpha_trait == BlendPixelTrait)
1199           density_x[x].direction[i].alpha+=
1200             cooccurrence[x][y].direction[i].alpha;
1201         if (image->colorspace == CMYKColorspace)
1202           density_x[x].direction[i].black+=
1203             cooccurrence[x][y].direction[i].black;
1204         density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1205         density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1206         density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1207         if (image->colorspace == CMYKColorspace)
1208           density_y[y].direction[i].black+=
1209             cooccurrence[x][y].direction[i].black;
1210         if (image->alpha_trait == BlendPixelTrait)
1211           density_y[y].direction[i].alpha+=
1212             cooccurrence[x][y].direction[i].alpha;
1213       }
1214       mean.direction[i].red+=y*sum[y].direction[i].red;
1215       sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
1216       mean.direction[i].green+=y*sum[y].direction[i].green;
1217       sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
1218       mean.direction[i].blue+=y*sum[y].direction[i].blue;
1219       sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
1220       if (image->colorspace == CMYKColorspace)
1221         {
1222           mean.direction[i].black+=y*sum[y].direction[i].black;
1223           sum_squares.direction[i].black+=y*y*sum[y].direction[i].black;
1224         }
1225       if (image->alpha_trait == BlendPixelTrait)
1226         {
1227           mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
1228           sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
1229         }
1230     }
1231     /*
1232       Correlation: measure of linear-dependencies in the image.
1233     */
1234     channel_features[RedPixelChannel].correlation[i]=
1235       (correlation.direction[i].red-mean.direction[i].red*
1236       mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
1237       (mean.direction[i].red*mean.direction[i].red))*sqrt(
1238       sum_squares.direction[i].red-(mean.direction[i].red*
1239       mean.direction[i].red)));
1240     channel_features[GreenPixelChannel].correlation[i]=
1241       (correlation.direction[i].green-mean.direction[i].green*
1242       mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
1243       (mean.direction[i].green*mean.direction[i].green))*sqrt(
1244       sum_squares.direction[i].green-(mean.direction[i].green*
1245       mean.direction[i].green)));
1246     channel_features[BluePixelChannel].correlation[i]=
1247       (correlation.direction[i].blue-mean.direction[i].blue*
1248       mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
1249       (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
1250       sum_squares.direction[i].blue-(mean.direction[i].blue*
1251       mean.direction[i].blue)));
1252     if (image->colorspace == CMYKColorspace)
1253       channel_features[BlackPixelChannel].correlation[i]=
1254         (correlation.direction[i].black-mean.direction[i].black*
1255         mean.direction[i].black)/(sqrt(sum_squares.direction[i].black-
1256         (mean.direction[i].black*mean.direction[i].black))*sqrt(
1257         sum_squares.direction[i].black-(mean.direction[i].black*
1258         mean.direction[i].black)));
1259     if (image->alpha_trait == BlendPixelTrait)
1260       channel_features[AlphaPixelChannel].correlation[i]=
1261         (correlation.direction[i].alpha-mean.direction[i].alpha*
1262         mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
1263         (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt(
1264         sum_squares.direction[i].alpha-(mean.direction[i].alpha*
1265         mean.direction[i].alpha)));
1266   }
1267   /*
1268     Compute more texture features.
1269   */
1270 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1271   #pragma omp parallel for schedule(static,4) shared(status) \
1272     magick_threads(image,image,number_grays,1)
1273 #endif
1274   for (i=0; i < 4; i++)
1275   {
1276     register ssize_t
1277       x;
1278
1279     for (x=2; x < (ssize_t) (2*number_grays); x++)
1280     {
1281       /*
1282         Sum average.
1283       */
1284       channel_features[RedPixelChannel].sum_average[i]+=
1285         x*density_xy[x].direction[i].red;
1286       channel_features[GreenPixelChannel].sum_average[i]+=
1287         x*density_xy[x].direction[i].green;
1288       channel_features[BluePixelChannel].sum_average[i]+=
1289         x*density_xy[x].direction[i].blue;
1290       if (image->colorspace == CMYKColorspace)
1291         channel_features[BlackPixelChannel].sum_average[i]+=
1292           x*density_xy[x].direction[i].black;
1293       if (image->alpha_trait == BlendPixelTrait)
1294         channel_features[AlphaPixelChannel].sum_average[i]+=
1295           x*density_xy[x].direction[i].alpha;
1296       /*
1297         Sum entropy.
1298       */
1299       channel_features[RedPixelChannel].sum_entropy[i]-=
1300         density_xy[x].direction[i].red*
1301         MagickLog10(density_xy[x].direction[i].red);
1302       channel_features[GreenPixelChannel].sum_entropy[i]-=
1303         density_xy[x].direction[i].green*
1304         MagickLog10(density_xy[x].direction[i].green);
1305       channel_features[BluePixelChannel].sum_entropy[i]-=
1306         density_xy[x].direction[i].blue*
1307         MagickLog10(density_xy[x].direction[i].blue);
1308       if (image->colorspace == CMYKColorspace)
1309         channel_features[BlackPixelChannel].sum_entropy[i]-=
1310           density_xy[x].direction[i].black*
1311           MagickLog10(density_xy[x].direction[i].black);
1312       if (image->alpha_trait == BlendPixelTrait)
1313         channel_features[AlphaPixelChannel].sum_entropy[i]-=
1314           density_xy[x].direction[i].alpha*
1315           MagickLog10(density_xy[x].direction[i].alpha);
1316       /*
1317         Sum variance.
1318       */
1319       channel_features[RedPixelChannel].sum_variance[i]+=
1320         (x-channel_features[RedPixelChannel].sum_entropy[i])*
1321         (x-channel_features[RedPixelChannel].sum_entropy[i])*
1322         density_xy[x].direction[i].red;
1323       channel_features[GreenPixelChannel].sum_variance[i]+=
1324         (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1325         (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1326         density_xy[x].direction[i].green;
1327       channel_features[BluePixelChannel].sum_variance[i]+=
1328         (x-channel_features[BluePixelChannel].sum_entropy[i])*
1329         (x-channel_features[BluePixelChannel].sum_entropy[i])*
1330         density_xy[x].direction[i].blue;
1331       if (image->colorspace == CMYKColorspace)
1332         channel_features[BlackPixelChannel].sum_variance[i]+=
1333           (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1334           (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1335           density_xy[x].direction[i].black;
1336       if (image->alpha_trait == BlendPixelTrait)
1337         channel_features[AlphaPixelChannel].sum_variance[i]+=
1338           (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1339           (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1340           density_xy[x].direction[i].alpha;
1341     }
1342   }
1343   /*
1344     Compute more texture features.
1345   */
1346 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1347   #pragma omp parallel for schedule(static,4) shared(status) \
1348     magick_threads(image,image,number_grays,1)
1349 #endif
1350   for (i=0; i < 4; i++)
1351   {
1352     register ssize_t
1353       y;
1354
1355     for (y=0; y < (ssize_t) number_grays; y++)
1356     {
1357       register ssize_t
1358         x;
1359
1360       for (x=0; x < (ssize_t) number_grays; x++)
1361       {
1362         /*
1363           Sum of Squares: Variance
1364         */
1365         variance.direction[i].red+=(y-mean.direction[i].red+1)*
1366           (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
1367         variance.direction[i].green+=(y-mean.direction[i].green+1)*
1368           (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
1369         variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
1370           (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
1371         if (image->colorspace == CMYKColorspace)
1372           variance.direction[i].black+=(y-mean.direction[i].black+1)*
1373             (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black;
1374         if (image->alpha_trait == BlendPixelTrait)
1375           variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
1376             (y-mean.direction[i].alpha+1)*
1377             cooccurrence[x][y].direction[i].alpha;
1378         /*
1379           Sum average / Difference Variance.
1380         */
1381         density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
1382           cooccurrence[x][y].direction[i].red;
1383         density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
1384           cooccurrence[x][y].direction[i].green;
1385         density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
1386           cooccurrence[x][y].direction[i].blue;
1387         if (image->colorspace == CMYKColorspace)
1388           density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
1389             cooccurrence[x][y].direction[i].black;
1390         if (image->alpha_trait == BlendPixelTrait)
1391           density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
1392             cooccurrence[x][y].direction[i].alpha;
1393         /*
1394           Information Measures of Correlation.
1395         */
1396         entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
1397           MagickLog10(cooccurrence[x][y].direction[i].red);
1398         entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
1399           MagickLog10(cooccurrence[x][y].direction[i].green);
1400         entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
1401           MagickLog10(cooccurrence[x][y].direction[i].blue);
1402         if (image->colorspace == CMYKColorspace)
1403           entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black*
1404             MagickLog10(cooccurrence[x][y].direction[i].black);
1405         if (image->alpha_trait == BlendPixelTrait)
1406           entropy_xy.direction[i].alpha-=
1407             cooccurrence[x][y].direction[i].alpha*MagickLog10(
1408             cooccurrence[x][y].direction[i].alpha);
1409         entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
1410           MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red));
1411         entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
1412           MagickLog10(density_x[x].direction[i].green*
1413           density_y[y].direction[i].green));
1414         entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
1415           MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue));
1416         if (image->colorspace == CMYKColorspace)
1417           entropy_xy1.direction[i].black-=(
1418             cooccurrence[x][y].direction[i].black*MagickLog10(
1419             density_x[x].direction[i].black*density_y[y].direction[i].black));
1420         if (image->alpha_trait == BlendPixelTrait)
1421           entropy_xy1.direction[i].alpha-=(
1422             cooccurrence[x][y].direction[i].alpha*MagickLog10(
1423             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1424         entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
1425           density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red*
1426           density_y[y].direction[i].red));
1427         entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
1428           density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green*
1429           density_y[y].direction[i].green));
1430         entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
1431           density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue*
1432           density_y[y].direction[i].blue));
1433         if (image->colorspace == CMYKColorspace)
1434           entropy_xy2.direction[i].black-=(density_x[x].direction[i].black*
1435             density_y[y].direction[i].black*MagickLog10(
1436             density_x[x].direction[i].black*density_y[y].direction[i].black));
1437         if (image->alpha_trait == BlendPixelTrait)
1438           entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha*
1439             density_y[y].direction[i].alpha*MagickLog10(
1440             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1441       }
1442     }
1443     channel_features[RedPixelChannel].variance_sum_of_squares[i]=
1444       variance.direction[i].red;
1445     channel_features[GreenPixelChannel].variance_sum_of_squares[i]=
1446       variance.direction[i].green;
1447     channel_features[BluePixelChannel].variance_sum_of_squares[i]=
1448       variance.direction[i].blue;
1449     if (image->colorspace == CMYKColorspace)
1450       channel_features[BlackPixelChannel].variance_sum_of_squares[i]=
1451         variance.direction[i].black;
1452     if (image->alpha_trait == BlendPixelTrait)
1453       channel_features[AlphaPixelChannel].variance_sum_of_squares[i]=
1454         variance.direction[i].alpha;
1455   }
1456   /*
1457     Compute more texture features.
1458   */
1459   (void) ResetMagickMemory(&variance,0,sizeof(variance));
1460   (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
1461 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1462   #pragma omp parallel for schedule(static,4) shared(status) \
1463     magick_threads(image,image,number_grays,1)
1464 #endif
1465   for (i=0; i < 4; i++)
1466   {
1467     register ssize_t
1468       x;
1469
1470     for (x=0; x < (ssize_t) number_grays; x++)
1471     {
1472       /*
1473         Difference variance.
1474       */
1475       variance.direction[i].red+=density_xy[x].direction[i].red;
1476       variance.direction[i].green+=density_xy[x].direction[i].green;
1477       variance.direction[i].blue+=density_xy[x].direction[i].blue;
1478       if (image->colorspace == CMYKColorspace)
1479         variance.direction[i].black+=density_xy[x].direction[i].black;
1480       if (image->alpha_trait == BlendPixelTrait)
1481         variance.direction[i].alpha+=density_xy[x].direction[i].alpha;
1482       sum_squares.direction[i].red+=density_xy[x].direction[i].red*
1483         density_xy[x].direction[i].red;
1484       sum_squares.direction[i].green+=density_xy[x].direction[i].green*
1485         density_xy[x].direction[i].green;
1486       sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
1487         density_xy[x].direction[i].blue;
1488       if (image->colorspace == CMYKColorspace)
1489         sum_squares.direction[i].black+=density_xy[x].direction[i].black*
1490           density_xy[x].direction[i].black;
1491       if (image->alpha_trait == BlendPixelTrait)
1492         sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
1493           density_xy[x].direction[i].alpha;
1494       /*
1495         Difference entropy.
1496       */
1497       channel_features[RedPixelChannel].difference_entropy[i]-=
1498         density_xy[x].direction[i].red*
1499         MagickLog10(density_xy[x].direction[i].red);
1500       channel_features[GreenPixelChannel].difference_entropy[i]-=
1501         density_xy[x].direction[i].green*
1502         MagickLog10(density_xy[x].direction[i].green);
1503       channel_features[BluePixelChannel].difference_entropy[i]-=
1504         density_xy[x].direction[i].blue*
1505         MagickLog10(density_xy[x].direction[i].blue);
1506       if (image->colorspace == CMYKColorspace)
1507         channel_features[BlackPixelChannel].difference_entropy[i]-=
1508           density_xy[x].direction[i].black*
1509           MagickLog10(density_xy[x].direction[i].black);
1510       if (image->alpha_trait == BlendPixelTrait)
1511         channel_features[AlphaPixelChannel].difference_entropy[i]-=
1512           density_xy[x].direction[i].alpha*
1513           MagickLog10(density_xy[x].direction[i].alpha);
1514       /*
1515         Information Measures of Correlation.
1516       */
1517       entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1518         MagickLog10(density_x[x].direction[i].red));
1519       entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1520         MagickLog10(density_x[x].direction[i].green));
1521       entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1522         MagickLog10(density_x[x].direction[i].blue));
1523       if (image->colorspace == CMYKColorspace)
1524         entropy_x.direction[i].black-=(density_x[x].direction[i].black*
1525           MagickLog10(density_x[x].direction[i].black));
1526       if (image->alpha_trait == BlendPixelTrait)
1527         entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha*
1528           MagickLog10(density_x[x].direction[i].alpha));
1529       entropy_y.direction[i].red-=(density_y[x].direction[i].red*
1530         MagickLog10(density_y[x].direction[i].red));
1531       entropy_y.direction[i].green-=(density_y[x].direction[i].green*
1532         MagickLog10(density_y[x].direction[i].green));
1533       entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
1534         MagickLog10(density_y[x].direction[i].blue));
1535       if (image->colorspace == CMYKColorspace)
1536         entropy_y.direction[i].black-=(density_y[x].direction[i].black*
1537           MagickLog10(density_y[x].direction[i].black));
1538       if (image->alpha_trait == BlendPixelTrait)
1539         entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
1540           MagickLog10(density_y[x].direction[i].alpha));
1541     }
1542     /*
1543       Difference variance.
1544     */
1545     channel_features[RedPixelChannel].difference_variance[i]=
1546       (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1547       (variance.direction[i].red*variance.direction[i].red))/
1548       ((double) number_grays*number_grays*number_grays*number_grays);
1549     channel_features[GreenPixelChannel].difference_variance[i]=
1550       (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1551       (variance.direction[i].green*variance.direction[i].green))/
1552       ((double) number_grays*number_grays*number_grays*number_grays);
1553     channel_features[BluePixelChannel].difference_variance[i]=
1554       (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1555       (variance.direction[i].blue*variance.direction[i].blue))/
1556       ((double) number_grays*number_grays*number_grays*number_grays);
1557     if (image->colorspace == CMYKColorspace)
1558       channel_features[BlackPixelChannel].difference_variance[i]=
1559         (((double) number_grays*number_grays*sum_squares.direction[i].black)-
1560         (variance.direction[i].black*variance.direction[i].black))/
1561         ((double) number_grays*number_grays*number_grays*number_grays);
1562     if (image->alpha_trait == BlendPixelTrait)
1563       channel_features[AlphaPixelChannel].difference_variance[i]=
1564         (((double) number_grays*number_grays*sum_squares.direction[i].alpha)-
1565         (variance.direction[i].alpha*variance.direction[i].alpha))/
1566         ((double) number_grays*number_grays*number_grays*number_grays);
1567     /*
1568       Information Measures of Correlation.
1569     */
1570     channel_features[RedPixelChannel].measure_of_correlation_1[i]=
1571       (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1572       (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1573        entropy_x.direction[i].red : entropy_y.direction[i].red);
1574     channel_features[GreenPixelChannel].measure_of_correlation_1[i]=
1575       (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1576       (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1577        entropy_x.direction[i].green : entropy_y.direction[i].green);
1578     channel_features[BluePixelChannel].measure_of_correlation_1[i]=
1579       (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1580       (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1581        entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1582     if (image->colorspace == CMYKColorspace)
1583       channel_features[BlackPixelChannel].measure_of_correlation_1[i]=
1584         (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/
1585         (entropy_x.direction[i].black > entropy_y.direction[i].black ?
1586          entropy_x.direction[i].black : entropy_y.direction[i].black);
1587     if (image->alpha_trait == BlendPixelTrait)
1588       channel_features[AlphaPixelChannel].measure_of_correlation_1[i]=
1589         (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/
1590         (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ?
1591          entropy_x.direction[i].alpha : entropy_y.direction[i].alpha);
1592     channel_features[RedPixelChannel].measure_of_correlation_2[i]=
1593       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].red-
1594       entropy_xy.direction[i].red)))));
1595     channel_features[GreenPixelChannel].measure_of_correlation_2[i]=
1596       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].green-
1597       entropy_xy.direction[i].green)))));
1598     channel_features[BluePixelChannel].measure_of_correlation_2[i]=
1599       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].blue-
1600       entropy_xy.direction[i].blue)))));
1601     if (image->colorspace == CMYKColorspace)
1602       channel_features[BlackPixelChannel].measure_of_correlation_2[i]=
1603         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].black-
1604         entropy_xy.direction[i].black)))));
1605     if (image->alpha_trait == BlendPixelTrait)
1606       channel_features[AlphaPixelChannel].measure_of_correlation_2[i]=
1607         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].alpha-
1608         entropy_xy.direction[i].alpha)))));
1609   }
1610   /*
1611     Compute more texture features.
1612   */
1613 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1614   #pragma omp parallel for schedule(static,4) shared(status) \
1615     magick_threads(image,image,number_grays,1)
1616 #endif
1617   for (i=0; i < 4; i++)
1618   {
1619     ssize_t
1620       z;
1621
1622     for (z=0; z < (ssize_t) number_grays; z++)
1623     {
1624       register ssize_t
1625         y;
1626
1627       ChannelStatistics
1628         pixel;
1629
1630       (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
1631       for (y=0; y < (ssize_t) number_grays; y++)
1632       {
1633         register ssize_t
1634           x;
1635
1636         for (x=0; x < (ssize_t) number_grays; x++)
1637         {
1638           /*
1639             Contrast:  amount of local variations present in an image.
1640           */
1641           if (((y-x) == z) || ((x-y) == z))
1642             {
1643               pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1644               pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1645               pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1646               if (image->colorspace == CMYKColorspace)
1647                 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
1648               if (image->alpha_trait == BlendPixelTrait)
1649                 pixel.direction[i].alpha+=
1650                   cooccurrence[x][y].direction[i].alpha;
1651             }
1652           /*
1653             Maximum Correlation Coefficient.
1654           */
1655           Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1656             cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1657             density_y[x].direction[i].red;
1658           Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1659             cooccurrence[y][x].direction[i].green/
1660             density_x[z].direction[i].green/density_y[x].direction[i].red;
1661           Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1662             cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/
1663             density_y[x].direction[i].blue;
1664           if (image->colorspace == CMYKColorspace)
1665             Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
1666               cooccurrence[y][x].direction[i].black/
1667               density_x[z].direction[i].black/density_y[x].direction[i].black;
1668           if (image->alpha_trait == BlendPixelTrait)
1669             Q[z][y].direction[i].alpha+=
1670               cooccurrence[z][x].direction[i].alpha*
1671               cooccurrence[y][x].direction[i].alpha/
1672               density_x[z].direction[i].alpha/
1673               density_y[x].direction[i].alpha;
1674         }
1675       }
1676       channel_features[RedPixelChannel].contrast[i]+=z*z*
1677         pixel.direction[i].red;
1678       channel_features[GreenPixelChannel].contrast[i]+=z*z*
1679         pixel.direction[i].green;
1680       channel_features[BluePixelChannel].contrast[i]+=z*z*
1681         pixel.direction[i].blue;
1682       if (image->colorspace == CMYKColorspace)
1683         channel_features[BlackPixelChannel].contrast[i]+=z*z*
1684           pixel.direction[i].black;
1685       if (image->alpha_trait == BlendPixelTrait)
1686         channel_features[AlphaPixelChannel].contrast[i]+=z*z*
1687           pixel.direction[i].alpha;
1688     }
1689     /*
1690       Maximum Correlation Coefficient.
1691       Future: return second largest eigenvalue of Q.
1692     */
1693     channel_features[RedPixelChannel].maximum_correlation_coefficient[i]=
1694       sqrt((double) -1.0);
1695     channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]=
1696       sqrt((double) -1.0);
1697     channel_features[BluePixelChannel].maximum_correlation_coefficient[i]=
1698       sqrt((double) -1.0);
1699     if (image->colorspace == CMYKColorspace)
1700       channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]=
1701         sqrt((double) -1.0);
1702     if (image->alpha_trait == BlendPixelTrait)
1703       channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]=
1704         sqrt((double) -1.0);
1705   }
1706   /*
1707     Relinquish resources.
1708   */
1709   sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1710   for (i=0; i < (ssize_t) number_grays; i++)
1711     Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1712   Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1713   density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1714   density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1715   density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1716   for (i=0; i < (ssize_t) number_grays; i++)
1717     cooccurrence[i]=(ChannelStatistics *)
1718       RelinquishMagickMemory(cooccurrence[i]);
1719   cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1720   return(channel_features);
1721 }
1722 \f
1723 /*
1724 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1725 %                                                                             %
1726 %                                                                             %
1727 %                                                                             %
1728 %     H o u g h L i n e I m a g e                                             %
1729 %                                                                             %
1730 %                                                                             %
1731 %                                                                             %
1732 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1733 %
1734 %  Use HoughLineImage() in conjunction with any binary edge extracted image (we
1735 %  recommand Canny) to identify lines in the image.  The algorithm accumulates
1736 %  counts for every white pixel for every possible orientation (for angles from
1737 %  0 to 179 in 1 degree increments) and distance from the center of the image to
1738 %  the corner (in 1 px increments) and stores the counts in an accumulator matrix
1739 %  of angle vs distance. The size of the accumulator is 180x(diagonal/2). Next
1740 %  it searches this space for peaks in counts and converts the locations of the
1741 %  peaks to slope and intercept in the normal x,y input image space. Use the
1742 %  slope/intercepts to find the endpoints clipped to the bounds of the image. The
1743 %  lines are then drawn. The counts are a measure of the length of the lines
1744 %
1745 %  The format of the HoughLineImage method is:
1746 %
1747 %      Image *HoughLineImage(const Image *image,const size_t width,
1748 %        const size_t height,const size_t threshold,ExceptionInfo *exception)
1749 %
1750 %  A description of each parameter follows:
1751 %
1752 %    o image: the image.
1753 %
1754 %    o width, height: find line pairs as local maxima in this neighborhood.
1755 %
1756 %    o threshold: the line count threshold.
1757 %
1758 %    o exception: return any errors or warnings in this structure.
1759 %
1760 */
1761
1762 static inline double MagickRound(double x)
1763 {
1764   /*
1765     Round the fraction to nearest integer.
1766   */
1767   if ((x-floor(x)) < (ceil(x)-x))
1768     return(floor(x));
1769   return(ceil(x));
1770 }
1771
1772 MagickExport Image *HoughLineImage(const Image *image,const size_t width,
1773   const size_t height,const size_t threshold,ExceptionInfo *exception)
1774 {
1775 #define HoughLineImageTag  "HoughLine/Image"
1776
1777   CacheView
1778     *image_view;
1779
1780   char
1781     message[MaxTextExtent],
1782     path[MaxTextExtent];
1783
1784   const char
1785     *artifact;
1786
1787   double
1788     hough_height;
1789
1790   Image
1791     *lines_image = NULL;
1792
1793   ImageInfo
1794     *image_info;
1795
1796   int
1797     file;
1798
1799   MagickBooleanType
1800     status;
1801
1802   MagickOffsetType
1803     progress;
1804
1805   MatrixInfo
1806     *accumulator;
1807
1808   PointInfo
1809     center;
1810
1811   register ssize_t
1812     y;
1813
1814   size_t
1815     accumulator_height,
1816     accumulator_width,
1817     line_count;
1818
1819   /*
1820     Create the accumulator.
1821   */
1822   assert(image != (const Image *) NULL);
1823   assert(image->signature == MagickSignature);
1824   if (image->debug != MagickFalse)
1825     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1826   assert(exception != (ExceptionInfo *) NULL);
1827   assert(exception->signature == MagickSignature);
1828   accumulator_width=180;
1829   hough_height=((sqrt(2.0)*(double) (image->rows > image->columns ?
1830     image->rows : image->columns))/2.0);
1831   accumulator_height=(size_t) (2.0*hough_height);
1832   accumulator=AcquireMatrixInfo(accumulator_width,accumulator_height,
1833     sizeof(double),exception);
1834   if (accumulator == (MatrixInfo *) NULL)
1835     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1836   if (NullMatrix(accumulator) == MagickFalse)
1837     {
1838       accumulator=DestroyMatrixInfo(accumulator);
1839       ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1840     }
1841   /*
1842     Populate the accumulator.
1843   */
1844   status=MagickTrue;
1845   progress=0;
1846   center.x=(double) image->columns/2.0;
1847   center.y=(double) image->rows/2.0;
1848   image_view=AcquireVirtualCacheView(image,exception);
1849   for (y=0; y < (ssize_t) image->rows; y++)
1850   {
1851     register const Quantum
1852       *restrict p;
1853
1854     register ssize_t
1855       x;
1856
1857     if (status == MagickFalse)
1858       continue;
1859     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1860     if (p == (Quantum *) NULL)
1861       {
1862         status=MagickFalse;
1863         continue;
1864       }
1865     for (x=0; x < (ssize_t) image->columns; x++)
1866     {
1867       if (GetPixelIntensity(image,p) > (QuantumRange/2))
1868         {
1869           register ssize_t
1870             i;
1871
1872           for (i=0; i < 180; i++)
1873           {
1874             double
1875               count,
1876               radius;
1877
1878             radius=(((double) x-center.x)*cos(DegreesToRadians((double) i)))+
1879               (((double) y-center.y)*sin(DegreesToRadians((double) i)));
1880             (void) GetMatrixElement(accumulator,i,(ssize_t)
1881               MagickRound(radius+hough_height),&count);
1882             count++;
1883             (void) SetMatrixElement(accumulator,i,(ssize_t)
1884               MagickRound(radius+hough_height),&count);
1885           }
1886         }
1887       p+=GetPixelChannels(image);
1888     }
1889     if (image->progress_monitor != (MagickProgressMonitor) NULL)
1890       {
1891         MagickBooleanType
1892           proceed;
1893
1894 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1895         #pragma omp critical (MagickCore_CannyEdgeImage)
1896 #endif
1897         proceed=SetImageProgress(image,CannyEdgeImageTag,progress++,
1898           image->rows);
1899         if (proceed == MagickFalse)
1900           status=MagickFalse;
1901       }
1902   }
1903   image_view=DestroyCacheView(image_view);
1904   if (status == MagickFalse)
1905     {
1906       accumulator=DestroyMatrixInfo(accumulator);
1907       return((Image *) NULL);
1908     }
1909   /*
1910     Generate line segments from accumulator.
1911   */
1912   file=AcquireUniqueFileResource(path);
1913   if (file == -1)
1914     {
1915       accumulator=DestroyMatrixInfo(accumulator);
1916       return((Image *) NULL);
1917     }
1918   (void) FormatLocaleString(message,MaxTextExtent,
1919     "# Hough line transform: %.20gx%.20g%+.20g\n",(double) width,
1920     (double) height,(double) threshold);
1921   if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1922     status=MagickFalse;
1923   (void) FormatLocaleString(message,MaxTextExtent,"viewbox 0 0 %.20g %.20g\n",
1924     (double) image->columns,(double) image->rows);
1925   if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1926     status=MagickFalse;
1927   line_count=image->columns > image->rows ? image->columns/4 : image->rows/4;
1928   if (threshold != 0)
1929     line_count=threshold;
1930   for (y=0; y < (ssize_t) accumulator_height; y++)
1931   {
1932     register ssize_t
1933       x;
1934
1935     for (x=0; x < (ssize_t) accumulator_width; x++)
1936     {
1937       double
1938         count;
1939
1940       (void) GetMatrixElement(accumulator,x,y,&count);
1941       if (count >= (double) line_count)
1942         {
1943           double
1944             maxima;
1945
1946           SegmentInfo
1947             line;
1948
1949           ssize_t
1950             v;
1951
1952           /*
1953             Is point a local maxima?
1954           */
1955           maxima=count;
1956           for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
1957           {
1958             ssize_t
1959               u;
1960
1961             for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
1962             {
1963               if ((u != 0) || (v !=0))
1964                 {
1965                   (void) GetMatrixElement(accumulator,x+u,y+v,&count);
1966                   if (count > maxima)
1967                     {
1968                       maxima=count;
1969                       break;
1970                     }
1971                 }
1972             }
1973             if (u < (ssize_t) (width/2))
1974               break;
1975           }
1976           (void) GetMatrixElement(accumulator,x,y,&count);
1977           if (maxima > count)
1978             continue;
1979           if ((x >= 45) && (x <= 135))
1980             {
1981               /*
1982                 y = (r-x cos(t))/sin(t)
1983               */
1984               line.x1=0.0;
1985               line.y1=((double) (y-(accumulator_height/2.0))-((line.x1-
1986                 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
1987                 sin(DegreesToRadians((double) x))+(image->rows/2.0);
1988               line.x2=(double) image->columns;
1989               line.y2=((double) (y-(accumulator_height/2.0))-((line.x2-
1990                 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
1991                 sin(DegreesToRadians((double) x))+(image->rows/2.0);
1992             }
1993           else
1994             {
1995               /*
1996                 x = (r-y cos(t))/sin(t)
1997               */
1998               line.y1=0.0;
1999               line.x1=((double) (y-(accumulator_height/2.0))-((line.y1-
2000                 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
2001                 cos(DegreesToRadians((double) x))+(image->columns/2.0);
2002               line.y2=(double) image->rows;
2003               line.x2=((double) (y-(accumulator_height/2.0))-((line.y2-
2004                 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
2005                 cos(DegreesToRadians((double) x))+(image->columns/2.0);
2006             }
2007           (void) FormatLocaleString(message,MaxTextExtent,
2008             "line %g,%g %g,%g  # %g\n",line.x1,line.y1,line.x2,line.y2,maxima);
2009           if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
2010             status=MagickFalse;
2011         }
2012     }
2013   }
2014   (void) close(file);
2015   /*
2016     Render lines to image canvas.
2017   */
2018   image_info=AcquireImageInfo();
2019   image_info->background_color=image->background_color;
2020   (void) FormatLocaleString(image_info->filename,MaxTextExtent,"mvg:%s",path);
2021   artifact=GetImageArtifact(image,"background");
2022   if (artifact != (const char *) NULL)
2023     (void) SetImageOption(image_info,"background",artifact);
2024   artifact=GetImageArtifact(image,"fill");
2025   if (artifact != (const char *) NULL)
2026     (void) SetImageOption(image_info,"fill",artifact);
2027   artifact=GetImageArtifact(image,"stroke");
2028   if (artifact != (const char *) NULL)
2029     (void) SetImageOption(image_info,"stroke",artifact);
2030   artifact=GetImageArtifact(image,"strokewidth");
2031   if (artifact != (const char *) NULL)
2032     (void) SetImageOption(image_info,"strokewidth",artifact);
2033   lines_image=ReadImage(image_info,exception);
2034   artifact=GetImageArtifact(image,"hough-lines:accumulator");
2035   if ((lines_image != (Image *) NULL) &&
2036       (IsStringTrue(artifact) != MagickFalse))
2037     {
2038       Image
2039         *accumulator_image;
2040
2041       accumulator_image=MatrixToImage(accumulator,exception);
2042       if (accumulator_image != (Image *) NULL)
2043         AppendImageToList(&lines_image,accumulator_image);
2044     }
2045   /*
2046     Free resources.
2047   */
2048   accumulator=DestroyMatrixInfo(accumulator);
2049   image_info=DestroyImageInfo(image_info);
2050   (void) RelinquishUniqueFileResource(path);
2051   return(GetFirstImageInList(lines_image));
2052 }
2053 \f
2054 /*
2055 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2056 %                                                                             %
2057 %                                                                             %
2058 %                                                                             %
2059 %     M e a n S h i f t I m a g e                                             %
2060 %                                                                             %
2061 %                                                                             %
2062 %                                                                             %
2063 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2064 %
2065 %  MeanShiftImage() delineate arbitrarily shaped clusters in the image. For
2066 %  each pixel, it visits all the pixels in the neighborhood specified by
2067 %  the window centered at the pixel and excludes those that are outside the
2068 %  radius=(window-1)/2 surrounding the pixel. From those pixels, it finds those
2069 %  that are within the specified color distance from the current mean, and
2070 %  computes a new x,y centroid from those coordinates and a new mean. This new
2071 %  x,y centroid is used as the center for a new window. This process iterates
2072 %  until it converges and the final mean is replaces the (original window
2073 %  center) pixel value. It repeats this process for the next pixel, etc.,
2074 %  until it processes all pixels in the image. Results are typically better with
2075 %  colorspaces other than sRGB. We recommend YIQ, YUV or YCbCr.
2076 %
2077 %  The format of the MeanShiftImage method is:
2078 %
2079 %      Image *MeanShiftImage(const Image *image,const size_t width,
2080 %        const size_t height,const double color_distance,
2081 %        ExceptionInfo *exception)
2082 %
2083 %  A description of each parameter follows:
2084 %
2085 %    o image: the image.
2086 %
2087 %    o width, height: find pixels in this neighborhood.
2088 %
2089 %    o color_distance: the color distance.
2090 %
2091 %    o exception: return any errors or warnings in this structure.
2092 %
2093 */
2094 MagickExport Image *MeanShiftImage(const Image *image,const size_t width,
2095   const size_t height,const double color_distance,ExceptionInfo *exception)
2096 {
2097 #define MaxMeanShiftIterations  100
2098 #define MeanShiftImageTag  "MeanShift/Image"
2099
2100   CacheView
2101     *image_view,
2102     *mean_view,
2103     *pixel_view;
2104
2105   Image
2106     *mean_image;
2107
2108   MagickBooleanType
2109     status;
2110
2111   MagickOffsetType
2112     progress;
2113
2114   ssize_t
2115     y;
2116
2117   assert(image != (const Image *) NULL);
2118   assert(image->signature == MagickSignature);
2119   if (image->debug != MagickFalse)
2120     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2121   assert(exception != (ExceptionInfo *) NULL);
2122   assert(exception->signature == MagickSignature);
2123   mean_image=CloneImage(image,image->columns,image->rows,MagickTrue,exception);
2124   if (mean_image == (Image *) NULL)
2125     return((Image *) NULL);
2126   if (SetImageStorageClass(mean_image,DirectClass,exception) == MagickFalse)
2127     {
2128       mean_image=DestroyImage(mean_image);
2129       return((Image *) NULL);
2130     }
2131   status=MagickTrue;
2132   progress=0;
2133   image_view=AcquireVirtualCacheView(image,exception);
2134   pixel_view=AcquireVirtualCacheView(image,exception);
2135   mean_view=AcquireAuthenticCacheView(mean_image,exception);
2136 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2137   #pragma omp parallel for schedule(static,4) shared(status,progress) \
2138     magick_threads(mean_image,mean_image,mean_image->rows,1)
2139 #endif
2140   for (y=0; y < (ssize_t) mean_image->rows; y++)
2141   {
2142     register const Quantum
2143       *restrict p;
2144
2145     register Quantum
2146       *restrict q;
2147
2148     register ssize_t
2149       x;
2150
2151     if (status == MagickFalse)
2152       continue;
2153     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2154     q=GetCacheViewAuthenticPixels(mean_view,0,y,mean_image->columns,1,
2155       exception);
2156     if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2157       {
2158         status=MagickFalse;
2159         continue;
2160       }
2161     for (x=0; x < (ssize_t) mean_image->columns; x++)
2162     {
2163       PixelInfo
2164         mean_pixel,
2165         previous_pixel;
2166
2167       PointInfo
2168         mean_location,
2169         previous_location;
2170
2171       register ssize_t
2172         i;
2173
2174       GetPixelInfo(image,&mean_pixel);
2175       GetPixelInfoPixel(image,p,&mean_pixel);
2176       mean_location.x=(double) x;
2177       mean_location.y=(double) y;
2178       for (i=0; i < MaxMeanShiftIterations; i++)
2179       {
2180         double
2181           distance,
2182           gamma;
2183
2184         PixelInfo
2185           sum_pixel;
2186
2187         PointInfo
2188           sum_location;
2189
2190         ssize_t
2191           count,
2192           v;
2193
2194         sum_location.x=0.0;
2195         sum_location.y=0.0;
2196         GetPixelInfo(image,&sum_pixel);
2197         previous_location=mean_location;
2198         previous_pixel=mean_pixel;
2199         count=0;
2200         for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
2201         {
2202           ssize_t
2203             u;
2204
2205           for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
2206           {
2207             if ((v*v+u*u) <= (ssize_t) ((width/2)*(height/2)))
2208               {
2209                 PixelInfo
2210                   pixel;
2211
2212                 status=GetOneCacheViewVirtualPixelInfo(pixel_view,(ssize_t)
2213                   MagickRound(mean_location.x+u),(ssize_t) MagickRound(
2214                   mean_location.y+v),&pixel,exception);
2215                 distance=(mean_pixel.red-pixel.red)*(mean_pixel.red-pixel.red)+
2216                   (mean_pixel.green-pixel.green)*(mean_pixel.green-pixel.green)+
2217                   (mean_pixel.blue-pixel.blue)*(mean_pixel.blue-pixel.blue);
2218                 if (distance <= (color_distance*color_distance))
2219                   {
2220                     sum_location.x+=mean_location.x+u;
2221                     sum_location.y+=mean_location.y+v;
2222                     sum_pixel.red+=pixel.red;
2223                     sum_pixel.green+=pixel.green;
2224                     sum_pixel.blue+=pixel.blue;
2225                     sum_pixel.alpha+=pixel.alpha;
2226                     count++;
2227                   }
2228               }
2229           }
2230         }
2231         gamma=1.0/count;
2232         mean_location.x=gamma*sum_location.x;
2233         mean_location.y=gamma*sum_location.y;
2234         mean_pixel.red=gamma*sum_pixel.red;
2235         mean_pixel.green=gamma*sum_pixel.green;
2236         mean_pixel.blue=gamma*sum_pixel.blue;
2237         mean_pixel.alpha=gamma*sum_pixel.alpha;
2238         distance=(mean_location.x-previous_location.x)*
2239           (mean_location.x-previous_location.x)+
2240           (mean_location.y-previous_location.y)*
2241           (mean_location.y-previous_location.y)+
2242           255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)*
2243           255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)+
2244           255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)*
2245           255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)+
2246           255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue)*
2247           255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue);
2248         if (distance <= 3.0)
2249           break;
2250       }
2251       SetPixelRed(mean_image,ClampToQuantum(mean_pixel.red),q);
2252       SetPixelGreen(mean_image,ClampToQuantum(mean_pixel.green),q);
2253       SetPixelBlue(mean_image,ClampToQuantum(mean_pixel.blue),q);
2254       SetPixelAlpha(mean_image,ClampToQuantum(mean_pixel.alpha),q);
2255       p+=GetPixelChannels(image);
2256       q+=GetPixelChannels(mean_image);
2257     }
2258     if (SyncCacheViewAuthenticPixels(mean_view,exception) == MagickFalse)
2259       status=MagickFalse;
2260     if (image->progress_monitor != (MagickProgressMonitor) NULL)
2261       {
2262         MagickBooleanType
2263           proceed;
2264
2265 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2266         #pragma omp critical (MagickCore_MeanShiftImage)
2267 #endif
2268         proceed=SetImageProgress(image,MeanShiftImageTag,progress++,
2269           image->rows);
2270         if (proceed == MagickFalse)
2271           status=MagickFalse;
2272       }
2273   }
2274   mean_view=DestroyCacheView(mean_view);
2275   pixel_view=DestroyCacheView(pixel_view);
2276   image_view=DestroyCacheView(image_view);
2277   return(mean_image);
2278 }