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/property.h"
  45 #include "MagickCore/animate.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/quantize.h"
  83 #include "MagickCore/quantum-private.h"
  84 #include "MagickCore/random_.h"
  85 #include "MagickCore/resource_.h"
  86 #include "MagickCore/segment.h"
  87 #include "MagickCore/semaphore.h"
  88 #include "MagickCore/signature-private.h"
  89 #include "MagickCore/string_.h"
  90 #include "MagickCore/thread-private.h"
  91 #include "MagickCore/timer.h"
  92 #include "MagickCore/utility.h"
  93 #include "MagickCore/version.h"
  94 \f
  95 /*
  96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  97 %                                                                             %
  98 %                                                                             %
  99 %                                                                             %
 100 %     C a n n y E d g e I m a g e                                             %
 101 %                                                                             %
 102 %                                                                             %
 103 %                                                                             %
 104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 105 %
 106 %  CannyEdgeImage() uses a multi-stage algorithm to detect a wide range of
 107 %  edges in images.
 108 %
 109 %  The format of the EdgeImage method is:
 110 %
 111 %      Image *CannyEdgeImage(const Image *image,const double radius,
 112 %        const double sigma,const double lower_percent,
 113 %        const double upper_percent,ExceptionInfo *exception)
 114 %
 115 %  A description of each parameter follows:
 116 %
 117 %    o image: the image.
 118 %
 119 %    o radius: the radius of the gaussian smoothing filter.
 120 %
 121 %    o sigma: the sigma of the gaussian smoothing filter.
 122 %
 123 %    o lower_precent: percentage of edge pixels in the lower threshold.
 124 %
 125 %    o upper_percent: percentage of edge pixels in the upper threshold.
 126 %
 127 %    o exception: return any errors or warnings in this structure.
 128 %
 129 */
 130
 131 typedef struct _CannyInfo
 132 {
 133   double
 134     magnitude,
 135     intensity;
 136
 137   int
 138     orientation;
 139
 140   ssize_t
 141     x,
 142     y;
 143 } CannyInfo;
 144
 145 static inline MagickBooleanType IsAuthenticPixel(const Image *image,
 146   const ssize_t x,const ssize_t y)
 147 {
 148   if ((x < 0) || (x >= (ssize_t) image->columns))
 149     return(MagickFalse);
 150   if ((y < 0) || (y >= (ssize_t) image->rows))
 151     return(MagickFalse);
 152   return(MagickTrue);
 153 }
 154
 155 static MagickBooleanType TraceEdges(Image *edge_image,CacheView *edge_view,
 156   MatrixInfo *canny_cache,const ssize_t x,const ssize_t y,
 157   const double lower_threshold,ExceptionInfo *exception)
 158 {
 159   CannyInfo
 160     edge,
 161     pixel;
 162
 163   MagickBooleanType
 164     status;
 165
 166   register Quantum
 167     *q;
 168
 169   register ssize_t
 170     i;
 171
 172   q=GetCacheViewAuthenticPixels(edge_view,x,y,1,1,exception);
 173   if (q == (Quantum *) NULL)
 174     return(MagickFalse);
 175   *q=QuantumRange;
 176   status=SyncCacheViewAuthenticPixels(edge_view,exception);
 177   if (status == MagickFalse)
 178     return(MagickFalse);;
 179   if (GetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
 180     return(MagickFalse);
 181   edge.x=x;
 182   edge.y=y;
 183   if (SetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
 184     return(MagickFalse);
 185   for (i=1; i != 0; )
 186   {
 187     ssize_t
 188       v;
 189
 190     i--;
 191     status=GetMatrixElement(canny_cache,i,0,&edge);
 192     if (status == MagickFalse)
 193       return(MagickFalse);
 194     for (v=(-1); v <= 1; v++)
 195     {
 196       ssize_t
 197         u;
 198
 199       for (u=(-1); u <= 1; u++)
 200       {
 201         if ((u == 0) && (v == 0))
 202           continue;
 203         if (IsAuthenticPixel(edge_image,edge.x+u,edge.y+v) == MagickFalse)
 204           continue;
 205         /*
 206           Not an edge if gradient value is below the lower threshold.
 207         */
 208         q=GetCacheViewAuthenticPixels(edge_view,edge.x+u,edge.y+v,1,1,
 209           exception);
 210         if (q == (Quantum *) NULL)
 211           return(MagickFalse);
 212         status=GetMatrixElement(canny_cache,edge.x+u,edge.y+v,&pixel);
 213         if (status == MagickFalse)
 214           return(MagickFalse);
 215         if ((GetPixelIntensity(edge_image,q) == 0.0) &&
 216             (pixel.intensity >= lower_threshold))
 217           {
 218             *q=QuantumRange;
 219             status=SyncCacheViewAuthenticPixels(edge_view,exception);
 220             if (status == MagickFalse)
 221               return(MagickFalse);
 222             edge.x+=u;
 223             edge.y+=v;
 224             status=SetMatrixElement(canny_cache,i,0,&edge);
 225             if (status == MagickFalse)
 226               return(MagickFalse);
 227             i++;
 228           }
 229       }
 230     }
 231   }
 232   return(MagickTrue);
 233 }
 234
 235 MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
 236   const double sigma,const double lower_percent,const double upper_percent,
 237   ExceptionInfo *exception)
 238 {
 239   CacheView
 240     *edge_view;
 241
 242   CannyInfo
 243     pixel;
 244
 245   char
 246     geometry[MaxTextExtent];
 247
 248   double
 249     lower_threshold,
 250     max,
 251     min,
 252     upper_threshold;
 253
 254   Image
 255     *edge_image;
 256
 257   KernelInfo
 258     *kernel_info;
 259
 260   MagickBooleanType
 261     status;
 262
 263   MatrixInfo
 264     *canny_cache;
 265
 266   ssize_t
 267     y;
 268
 269   assert(image != (const Image *) NULL);
 270   assert(image->signature == MagickSignature);
 271   if (image->debug != MagickFalse)
 272     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 273   assert(exception != (ExceptionInfo *) NULL);
 274   assert(exception->signature == MagickSignature);
 275   /*
 276     Filter out noise.
 277   */
 278   (void) FormatLocaleString(geometry,MaxTextExtent,
 279     "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
 280   kernel_info=AcquireKernelInfo(geometry);
 281   if (kernel_info == (KernelInfo *) NULL)
 282     ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
 283   edge_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info,
 284     UndefinedCompositeOp,0.0,exception);
 285   kernel_info=DestroyKernelInfo(kernel_info);
 286   if (edge_image == (Image *) NULL)
 287     return((Image *) NULL);
 288   if (SetImageColorspace(edge_image,GRAYColorspace,exception) == MagickFalse)
 289     {
 290       edge_image=DestroyImage(edge_image);
 291       return((Image *) NULL);
 292     }
 293   /*
 294     Find the intensity gradient of the image.
 295   */
 296   canny_cache=AcquireMatrixInfo(edge_image->columns,edge_image->rows,
 297     sizeof(CannyInfo),exception);
 298   if (canny_cache == (MatrixInfo *) NULL)
 299     {
 300       edge_image=DestroyImage(edge_image);
 301       return((Image *) NULL);
 302     }
 303   status=MagickTrue;
 304   edge_view=AcquireVirtualCacheView(edge_image,exception);
 305 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 306   #pragma omp parallel for schedule(static,4) shared(status) \
 307     magick_threads(edge_image,edge_image,edge_image->rows,1)
 308 #endif
 309   for (y=0; y < (ssize_t) edge_image->rows; y++)
 310   {
 311     register const Quantum
 312       *restrict p;
 313
 314     register ssize_t
 315       x;
 316
 317     if (status == MagickFalse)
 318       continue;
 319     p=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns+1,2,
 320       exception);
 321     if (p == (const Quantum *) NULL)
 322       {
 323         status=MagickFalse;
 324         continue;
 325       }
 326     for (x=0; x < (ssize_t) edge_image->columns; x++)
 327     {
 328       CannyInfo
 329         pixel;
 330
 331       double
 332         dx,
 333         dy;
 334
 335       register const Quantum
 336         *restrict kernel_pixels;
 337
 338       ssize_t
 339         v;
 340
 341       static double
 342         Gx[2][2] =
 343         {
 344           { -1.0,  +1.0 },
 345           { -1.0,  +1.0 }
 346         },
 347         Gy[2][2] =
 348         {
 349           { +1.0, +1.0 },
 350           { -1.0, -1.0 }
 351         };
 352
 353       (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
 354       dx=0.0;
 355       dy=0.0;
 356       kernel_pixels=p;
 357       for (v=0; v < 2; v++)
 358       {
 359         ssize_t
 360           u;
 361
 362         for (u=0; u < 2; u++)
 363         {
 364           double
 365             intensity;
 366
 367           intensity=GetPixelIntensity(edge_image,kernel_pixels+u);
 368           dx+=0.5*Gx[v][u]*intensity;
 369           dy+=0.5*Gy[v][u]*intensity;
 370         }
 371         kernel_pixels+=edge_image->columns+1;
 372       }
 373       pixel.magnitude=hypot(dx,dy);
 374       pixel.orientation=0;
 375       if (fabs(dx) > MagickEpsilon)
 376         {
 377           double
 378             slope;
 379
 380           slope=dy/dx;
 381           if (slope < 0.0)
 382             {
 383               if (slope < -2.41421356237)
 384                 pixel.orientation=0;
 385               else
 386                 if (slope < -0.414213562373)
 387                   pixel.orientation=1;
 388                 else
 389                   pixel.orientation=2;
 390             }
 391           else
 392             {
 393               if (slope > 2.41421356237)
 394                 pixel.orientation=0;
 395               else
 396                 if (slope > 0.414213562373)
 397                   pixel.orientation=3;
 398                 else
 399                   pixel.orientation=2;
 400             }
 401         }
 402       if (SetMatrixElement(canny_cache,x,y,&pixel) == MagickFalse)
 403         continue;
 404       p+=GetPixelChannels(edge_image);
 405     }
 406   }
 407   edge_view=DestroyCacheView(edge_view);
 408   /*
 409     Non-maxima suppression, remove pixels that are not considered to be part
 410     of an edge.
 411   */
 412   (void) GetMatrixElement(canny_cache,0,0,&pixel);
 413   max=pixel.intensity;
 414   min=pixel.intensity;
 415   edge_view=AcquireAuthenticCacheView(edge_image,exception);
 416 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 417   #pragma omp parallel for schedule(static,4) shared(status) \
 418     magick_threads(edge_image,edge_image,edge_image->rows,1)
 419 #endif
 420   for (y=0; y < (ssize_t) edge_image->rows; y++)
 421   {
 422     register Quantum
 423       *restrict q;
 424
 425     register ssize_t
 426       x;
 427
 428     if (status == MagickFalse)
 429       continue;
 430     q=GetCacheViewAuthenticPixels(edge_view,0,y,edge_image->columns,1,
 431       exception);
 432     if (q == (Quantum *) NULL)
 433       {
 434         status=MagickFalse;
 435         continue;
 436       }
 437     for (x=0; x < (ssize_t) edge_image->columns; x++)
 438     {
 439       CannyInfo
 440         alpha_pixel,
 441         beta_pixel,
 442         pixel;
 443
 444       (void) GetMatrixElement(canny_cache,x,y,&pixel);
 445       switch (pixel.orientation)
 446       {
 447         case 0:
 448         {
 449           /*
 450             0 degrees, north and south.
 451           */
 452           (void) GetMatrixElement(canny_cache,x,y-1,&alpha_pixel);
 453           (void) GetMatrixElement(canny_cache,x,y+1,&beta_pixel);
 454           break;
 455         }
 456         case 1:
 457         {
 458           /*
 459             45 degrees, northwest and southeast.
 460           */
 461           (void) GetMatrixElement(canny_cache,x-1,y-1,&alpha_pixel);
 462           (void) GetMatrixElement(canny_cache,x+1,y+1,&beta_pixel);
 463           break;
 464         }
 465         case 2:
 466         {
 467           /*
 468             90 degrees, east and west.
 469           */
 470           (void) GetMatrixElement(canny_cache,x-1,y,&alpha_pixel);
 471           (void) GetMatrixElement(canny_cache,x+1,y,&beta_pixel);
 472           break;
 473         }
 474         case 3:
 475         {
 476           /*
 477             135 degrees, northeast and southwest.
 478           */
 479           (void) GetMatrixElement(canny_cache,x+1,y-1,&beta_pixel);
 480           (void) GetMatrixElement(canny_cache,x-1,y+1,&alpha_pixel);
 481           break;
 482         }
 483       }
 484       pixel.intensity=pixel.magnitude;
 485       if ((pixel.magnitude < alpha_pixel.magnitude) ||
 486           (pixel.magnitude < beta_pixel.magnitude))
 487         pixel.intensity=0;
 488       (void) SetMatrixElement(canny_cache,x,y,&pixel);
 489 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 490       #pragma omp critical (MagickCore_CannyEdgeImage)
 491 #endif
 492       {
 493         if (pixel.intensity < min)
 494           min=pixel.intensity;
 495         if (pixel.intensity > max)
 496           max=pixel.intensity;
 497       }
 498       *q=0;
 499       q+=GetPixelChannels(edge_image);
 500     }
 501     if (SyncCacheViewAuthenticPixels(edge_view,exception) == MagickFalse)
 502       status=MagickFalse;
 503   }
 504   edge_view=DestroyCacheView(edge_view);
 505   /*
 506     Estimate hysteresis threshold.
 507   */
 508   lower_threshold=lower_percent*(max-min)+min;
 509   upper_threshold=upper_percent*(max-min)+min;
 510   /*
 511     Hysteresis threshold.
 512   */
 513   edge_view=AcquireAuthenticCacheView(edge_image,exception);
 514   for (y=0; y < (ssize_t) edge_image->rows; y++)
 515   {
 516     register ssize_t
 517       x;
 518
 519     if (status == MagickFalse)
 520       continue;
 521     for (x=0; x < (ssize_t) edge_image->columns; x++)
 522     {
 523       CannyInfo
 524         pixel;
 525
 526       register const Quantum
 527         *restrict p;
 528
 529       /*
 530         Edge if pixel gradient higher than upper threshold.
 531       */
 532       p=GetCacheViewVirtualPixels(edge_view,x,y,1,1,exception);
 533       if (p == (const Quantum *) NULL)
 534         continue;
 535       status=GetMatrixElement(canny_cache,x,y,&pixel);
 536       if (status == MagickFalse)
 537         continue;
 538       if ((GetPixelIntensity(edge_image,p) == 0.0) &&
 539           (pixel.intensity >= upper_threshold))
 540         status=TraceEdges(edge_image,edge_view,canny_cache,x,y,lower_threshold,
 541           exception);
 542     }
 543   }
 544   edge_view=DestroyCacheView(edge_view);
 545   /*
 546     Free resources.
 547   */
 548   canny_cache=DestroyMatrixInfo(canny_cache);
 549   return(edge_image);
 550 }
 551 \f
 552 /*
 553 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 554 %                                                                             %
 555 %                                                                             %
 556 %                                                                             %
 557 %     H o u g h T r a n s f o r m I m a g e                                   %
 558 %                                                                             %
 559 %                                                                             %
 560 %                                                                             %
 561 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 562 %
 563 %  HoughTransformImage() detects lines in an image.
 564 %
 565 %  The format of the HoughTransformImage method is:
 566 %
 567 %      Image *HoughTransformImage(const Image *image,const size_t width,
 568 %        const size_t height,const size_t threshold,ExceptionInfo *exception)
 569 %
 570 %  A description of each parameter follows:
 571 %
 572 %    o image: the image.
 573 %
 574 %    o width, height: find line pairs as local maxima in this neighborhood.
 575 %
 576 %    o threshold: the line count threshold.
 577 %
 578 %    o exception: return any errors or warnings in this structure.
 579 %
 580 */
 581 MagickExport Image *HoughTransformImage(const Image *image,const size_t width,
 582   const size_t height,const size_t threshold,ExceptionInfo *exception)
 583 {
 584   return((Image *) NULL);
 585 }
 586
 587 \f
 588 /*
 589 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 590 %                                                                             %
 591 %                                                                             %
 592 %                                                                             %
 593 %   G e t I m a g e F e a t u r e s                                           %
 594 %                                                                             %
 595 %                                                                             %
 596 %                                                                             %
 597 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 598 %
 599 %  GetImageFeatures() returns features for each channel in the image in
 600 %  each of four directions (horizontal, vertical, left and right diagonals)
 601 %  for the specified distance.  The features include the angular second
 602 %  moment, contrast, correlation, sum of squares: variance, inverse difference
 603 %  moment, sum average, sum varience, sum entropy, entropy, difference variance,%  difference entropy, information measures of correlation 1, information
 604 %  measures of correlation 2, and maximum correlation coefficient.  You can
 605 %  access the red channel contrast, for example, like this:
 606 %
 607 %      channel_features=GetImageFeatures(image,1,exception);
 608 %      contrast=channel_features[RedPixelChannel].contrast[0];
 609 %
 610 %  Use MagickRelinquishMemory() to free the features buffer.
 611 %
 612 %  The format of the GetImageFeatures method is:
 613 %
 614 %      ChannelFeatures *GetImageFeatures(const Image *image,
 615 %        const size_t distance,ExceptionInfo *exception)
 616 %
 617 %  A description of each parameter follows:
 618 %
 619 %    o image: the image.
 620 %
 621 %    o distance: the distance.
 622 %
 623 %    o exception: return any errors or warnings in this structure.
 624 %
 625 */
 626
 627 static inline ssize_t MagickAbsoluteValue(const ssize_t x)
 628 {
 629   if (x < 0)
 630     return(-x);
 631   return(x);
 632 }
 633
 634 static inline double MagickLog10(const double x)
 635 {
 636 #define Log10Epsilon  (1.0e-11)
 637
 638  if (fabs(x) < Log10Epsilon)
 639    return(log10(Log10Epsilon));
 640  return(log10(fabs(x)));
 641 }
 642
 643 MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
 644   const size_t distance,ExceptionInfo *exception)
 645 {
 646   typedef struct _ChannelStatistics
 647   {
 648     PixelInfo
 649       direction[4];  /* horizontal, vertical, left and right diagonals */
 650   } ChannelStatistics;
 651
 652   CacheView
 653     *image_view;
 654
 655   ChannelFeatures
 656     *channel_features;
 657
 658   ChannelStatistics
 659     **cooccurrence,
 660     correlation,
 661     *density_x,
 662     *density_xy,
 663     *density_y,
 664     entropy_x,
 665     entropy_xy,
 666     entropy_xy1,
 667     entropy_xy2,
 668     entropy_y,
 669     mean,
 670     **Q,
 671     *sum,
 672     sum_squares,
 673     variance;
 674
 675   PixelPacket
 676     gray,
 677     *grays;
 678
 679   MagickBooleanType
 680     status;
 681
 682   register ssize_t
 683     i;
 684
 685   size_t
 686     length;
 687
 688   ssize_t
 689     y;
 690
 691   unsigned int
 692     number_grays;
 693
 694   assert(image != (Image *) NULL);
 695   assert(image->signature == MagickSignature);
 696   if (image->debug != MagickFalse)
 697     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 698   if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
 699     return((ChannelFeatures *) NULL);
 700   length=CompositeChannels+1UL;
 701   channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
 702     sizeof(*channel_features));
 703   if (channel_features == (ChannelFeatures *) NULL)
 704     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 705   (void) ResetMagickMemory(channel_features,0,length*
 706     sizeof(*channel_features));
 707   /*
 708     Form grays.
 709   */
 710   grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays));
 711   if (grays == (PixelPacket *) NULL)
 712     {
 713       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 714         channel_features);
 715       (void) ThrowMagickException(exception,GetMagickModule(),
 716         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 717       return(channel_features);
 718     }
 719   for (i=0; i <= (ssize_t) MaxMap; i++)
 720   {
 721     grays[i].red=(~0U);
 722     grays[i].green=(~0U);
 723     grays[i].blue=(~0U);
 724     grays[i].alpha=(~0U);
 725     grays[i].black=(~0U);
 726   }
 727   status=MagickTrue;
 728   image_view=AcquireVirtualCacheView(image,exception);
 729 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 730   #pragma omp parallel for schedule(static,4) shared(status) \
 731     magick_threads(image,image,image->rows,1)
 732 #endif
 733   for (y=0; y < (ssize_t) image->rows; y++)
 734   {
 735     register const Quantum
 736       *restrict p;
 737
 738     register ssize_t
 739       x;
 740
 741     if (status == MagickFalse)
 742       continue;
 743     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 744     if (p == (const Quantum *) NULL)
 745       {
 746         status=MagickFalse;
 747         continue;
 748       }
 749     for (x=0; x < (ssize_t) image->columns; x++)
 750     {
 751       grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
 752         ScaleQuantumToMap(GetPixelRed(image,p));
 753       grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
 754         ScaleQuantumToMap(GetPixelGreen(image,p));
 755       grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
 756         ScaleQuantumToMap(GetPixelBlue(image,p));
 757       if (image->colorspace == CMYKColorspace)
 758         grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
 759           ScaleQuantumToMap(GetPixelBlack(image,p));
 760       if (image->alpha_trait == BlendPixelTrait)
 761         grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
 762           ScaleQuantumToMap(GetPixelAlpha(image,p));
 763       p+=GetPixelChannels(image);
 764     }
 765   }
 766   image_view=DestroyCacheView(image_view);
 767   if (status == MagickFalse)
 768     {
 769       grays=(PixelPacket *) RelinquishMagickMemory(grays);
 770       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 771         channel_features);
 772       return(channel_features);
 773     }
 774   (void) ResetMagickMemory(&gray,0,sizeof(gray));
 775   for (i=0; i <= (ssize_t) MaxMap; i++)
 776   {
 777     if (grays[i].red != ~0U)
 778       grays[gray.red++].red=grays[i].red;
 779     if (grays[i].green != ~0U)
 780       grays[gray.green++].green=grays[i].green;
 781     if (grays[i].blue != ~0U)
 782       grays[gray.blue++].blue=grays[i].blue;
 783     if (image->colorspace == CMYKColorspace)
 784       if (grays[i].black != ~0U)
 785         grays[gray.black++].black=grays[i].black;
 786     if (image->alpha_trait == BlendPixelTrait)
 787       if (grays[i].alpha != ~0U)
 788         grays[gray.alpha++].alpha=grays[i].alpha;
 789   }
 790   /*
 791     Allocate spatial dependence matrix.
 792   */
 793   number_grays=gray.red;
 794   if (gray.green > number_grays)
 795     number_grays=gray.green;
 796   if (gray.blue > number_grays)
 797     number_grays=gray.blue;
 798   if (image->colorspace == CMYKColorspace)
 799     if (gray.black > number_grays)
 800       number_grays=gray.black;
 801   if (image->alpha_trait == BlendPixelTrait)
 802     if (gray.alpha > number_grays)
 803       number_grays=gray.alpha;
 804   cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays,
 805     sizeof(*cooccurrence));
 806   density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
 807     sizeof(*density_x));
 808   density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
 809     sizeof(*density_xy));
 810   density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
 811     sizeof(*density_y));
 812   Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q));
 813   sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum));
 814   if ((cooccurrence == (ChannelStatistics **) NULL) ||
 815       (density_x == (ChannelStatistics *) NULL) ||
 816       (density_xy == (ChannelStatistics *) NULL) ||
 817       (density_y == (ChannelStatistics *) NULL) ||
 818       (Q == (ChannelStatistics **) NULL) ||
 819       (sum == (ChannelStatistics *) NULL))
 820     {
 821       if (Q != (ChannelStatistics **) NULL)
 822         {
 823           for (i=0; i < (ssize_t) number_grays; i++)
 824             Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
 825           Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
 826         }
 827       if (sum != (ChannelStatistics *) NULL)
 828         sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
 829       if (density_y != (ChannelStatistics *) NULL)
 830         density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
 831       if (density_xy != (ChannelStatistics *) NULL)
 832         density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
 833       if (density_x != (ChannelStatistics *) NULL)
 834         density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
 835       if (cooccurrence != (ChannelStatistics **) NULL)
 836         {
 837           for (i=0; i < (ssize_t) number_grays; i++)
 838             cooccurrence[i]=(ChannelStatistics *)
 839               RelinquishMagickMemory(cooccurrence[i]);
 840           cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
 841             cooccurrence);
 842         }
 843       grays=(PixelPacket *) RelinquishMagickMemory(grays);
 844       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 845         channel_features);
 846       (void) ThrowMagickException(exception,GetMagickModule(),
 847         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 848       return(channel_features);
 849     }
 850   (void) ResetMagickMemory(&correlation,0,sizeof(correlation));
 851   (void) ResetMagickMemory(density_x,0,2*(number_grays+1)*sizeof(*density_x));
 852   (void) ResetMagickMemory(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
 853   (void) ResetMagickMemory(density_y,0,2*(number_grays+1)*sizeof(*density_y));
 854   (void) ResetMagickMemory(&mean,0,sizeof(mean));
 855   (void) ResetMagickMemory(sum,0,number_grays*sizeof(*sum));
 856   (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
 857   (void) ResetMagickMemory(density_xy,0,2*number_grays*sizeof(*density_xy));
 858   (void) ResetMagickMemory(&entropy_x,0,sizeof(entropy_x));
 859   (void) ResetMagickMemory(&entropy_xy,0,sizeof(entropy_xy));
 860   (void) ResetMagickMemory(&entropy_xy1,0,sizeof(entropy_xy1));
 861   (void) ResetMagickMemory(&entropy_xy2,0,sizeof(entropy_xy2));
 862   (void) ResetMagickMemory(&entropy_y,0,sizeof(entropy_y));
 863   (void) ResetMagickMemory(&variance,0,sizeof(variance));
 864   for (i=0; i < (ssize_t) number_grays; i++)
 865   {
 866     cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
 867       sizeof(**cooccurrence));
 868     Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
 869     if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
 870         (Q[i] == (ChannelStatistics *) NULL))
 871       break;
 872     (void) ResetMagickMemory(cooccurrence[i],0,number_grays*
 873       sizeof(**cooccurrence));
 874     (void) ResetMagickMemory(Q[i],0,number_grays*sizeof(**Q));
 875   }
 876   if (i < (ssize_t) number_grays)
 877     {
 878       for (i--; i >= 0; i--)
 879       {
 880         if (Q[i] != (ChannelStatistics *) NULL)
 881           Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
 882         if (cooccurrence[i] != (ChannelStatistics *) NULL)
 883           cooccurrence[i]=(ChannelStatistics *)
 884             RelinquishMagickMemory(cooccurrence[i]);
 885       }
 886       Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
 887       cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
 888       sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
 889       density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
 890       density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
 891       density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
 892       grays=(PixelPacket *) RelinquishMagickMemory(grays);
 893       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 894         channel_features);
 895       (void) ThrowMagickException(exception,GetMagickModule(),
 896         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 897       return(channel_features);
 898     }
 899   /*
 900     Initialize spatial dependence matrix.
 901   */
 902   status=MagickTrue;
 903   image_view=AcquireVirtualCacheView(image,exception);
 904   for (y=0; y < (ssize_t) image->rows; y++)
 905   {
 906     register const Quantum
 907       *restrict p;
 908
 909     register ssize_t
 910       x;
 911
 912     ssize_t
 913       i,
 914       offset,
 915       u,
 916       v;
 917
 918     if (status == MagickFalse)
 919       continue;
 920     p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,y,image->columns+
 921       2*distance,distance+2,exception);
 922     if (p == (const Quantum *) NULL)
 923       {
 924         status=MagickFalse;
 925         continue;
 926       }
 927     p+=distance*GetPixelChannels(image);;
 928     for (x=0; x < (ssize_t) image->columns; x++)
 929     {
 930       for (i=0; i < 4; i++)
 931       {
 932         switch (i)
 933         {
 934           case 0:
 935           default:
 936           {
 937             /*
 938               Horizontal adjacency.
 939             */
 940             offset=(ssize_t) distance;
 941             break;
 942           }
 943           case 1:
 944           {
 945             /*
 946               Vertical adjacency.
 947             */
 948             offset=(ssize_t) (image->columns+2*distance);
 949             break;
 950           }
 951           case 2:
 952           {
 953             /*
 954               Right diagonal adjacency.
 955             */
 956             offset=(ssize_t) ((image->columns+2*distance)-distance);
 957             break;
 958           }
 959           case 3:
 960           {
 961             /*
 962               Left diagonal adjacency.
 963             */
 964             offset=(ssize_t) ((image->columns+2*distance)+distance);
 965             break;
 966           }
 967         }
 968         u=0;
 969         v=0;
 970         while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p)))
 971           u++;
 972         while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image))))
 973           v++;
 974         cooccurrence[u][v].direction[i].red++;
 975         cooccurrence[v][u].direction[i].red++;
 976         u=0;
 977         v=0;
 978         while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p)))
 979           u++;
 980         while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image))))
 981           v++;
 982         cooccurrence[u][v].direction[i].green++;
 983         cooccurrence[v][u].direction[i].green++;
 984         u=0;
 985         v=0;
 986         while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p)))
 987           u++;
 988         while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image))))
 989           v++;
 990         cooccurrence[u][v].direction[i].blue++;
 991         cooccurrence[v][u].direction[i].blue++;
 992         if (image->colorspace == CMYKColorspace)
 993           {
 994             u=0;
 995             v=0;
 996             while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p)))
 997               u++;
 998             while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image))))
 999               v++;
1000             cooccurrence[u][v].direction[i].black++;
1001             cooccurrence[v][u].direction[i].black++;
1002           }
1003         if (image->alpha_trait == BlendPixelTrait)
1004           {
1005             u=0;
1006             v=0;
1007             while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p)))
1008               u++;
1009             while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image))))
1010               v++;
1011             cooccurrence[u][v].direction[i].alpha++;
1012             cooccurrence[v][u].direction[i].alpha++;
1013           }
1014       }
1015       p+=GetPixelChannels(image);
1016     }
1017   }
1018   grays=(PixelPacket *) RelinquishMagickMemory(grays);
1019   image_view=DestroyCacheView(image_view);
1020   if (status == MagickFalse)
1021     {
1022       for (i=0; i < (ssize_t) number_grays; i++)
1023         cooccurrence[i]=(ChannelStatistics *)
1024           RelinquishMagickMemory(cooccurrence[i]);
1025       cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1026       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1027         channel_features);
1028       (void) ThrowMagickException(exception,GetMagickModule(),
1029         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1030       return(channel_features);
1031     }
1032   /*
1033     Normalize spatial dependence matrix.
1034   */
1035   for (i=0; i < 4; i++)
1036   {
1037     double
1038       normalize;
1039
1040     register ssize_t
1041       y;
1042
1043     switch (i)
1044     {
1045       case 0:
1046       default:
1047       {
1048         /*
1049           Horizontal adjacency.
1050         */
1051         normalize=2.0*image->rows*(image->columns-distance);
1052         break;
1053       }
1054       case 1:
1055       {
1056         /*
1057           Vertical adjacency.
1058         */
1059         normalize=2.0*(image->rows-distance)*image->columns;
1060         break;
1061       }
1062       case 2:
1063       {
1064         /*
1065           Right diagonal adjacency.
1066         */
1067         normalize=2.0*(image->rows-distance)*(image->columns-distance);
1068         break;
1069       }
1070       case 3:
1071       {
1072         /*
1073           Left diagonal adjacency.
1074         */
1075         normalize=2.0*(image->rows-distance)*(image->columns-distance);
1076         break;
1077       }
1078     }
1079     normalize=PerceptibleReciprocal(normalize);
1080     for (y=0; y < (ssize_t) number_grays; y++)
1081     {
1082       register ssize_t
1083         x;
1084
1085       for (x=0; x < (ssize_t) number_grays; x++)
1086       {
1087         cooccurrence[x][y].direction[i].red*=normalize;
1088         cooccurrence[x][y].direction[i].green*=normalize;
1089         cooccurrence[x][y].direction[i].blue*=normalize;
1090         if (image->colorspace == CMYKColorspace)
1091           cooccurrence[x][y].direction[i].black*=normalize;
1092         if (image->alpha_trait == BlendPixelTrait)
1093           cooccurrence[x][y].direction[i].alpha*=normalize;
1094       }
1095     }
1096   }
1097   /*
1098     Compute texture features.
1099   */
1100 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1101   #pragma omp parallel for schedule(static,4) shared(status) \
1102     magick_threads(image,image,number_grays,1)
1103 #endif
1104   for (i=0; i < 4; i++)
1105   {
1106     register ssize_t
1107       y;
1108
1109     for (y=0; y < (ssize_t) number_grays; y++)
1110     {
1111       register ssize_t
1112         x;
1113
1114       for (x=0; x < (ssize_t) number_grays; x++)
1115       {
1116         /*
1117           Angular second moment:  measure of homogeneity of the image.
1118         */
1119         channel_features[RedPixelChannel].angular_second_moment[i]+=
1120           cooccurrence[x][y].direction[i].red*
1121           cooccurrence[x][y].direction[i].red;
1122         channel_features[GreenPixelChannel].angular_second_moment[i]+=
1123           cooccurrence[x][y].direction[i].green*
1124           cooccurrence[x][y].direction[i].green;
1125         channel_features[BluePixelChannel].angular_second_moment[i]+=
1126           cooccurrence[x][y].direction[i].blue*
1127           cooccurrence[x][y].direction[i].blue;
1128         if (image->colorspace == CMYKColorspace)
1129           channel_features[BlackPixelChannel].angular_second_moment[i]+=
1130             cooccurrence[x][y].direction[i].black*
1131             cooccurrence[x][y].direction[i].black;
1132         if (image->alpha_trait == BlendPixelTrait)
1133           channel_features[AlphaPixelChannel].angular_second_moment[i]+=
1134             cooccurrence[x][y].direction[i].alpha*
1135             cooccurrence[x][y].direction[i].alpha;
1136         /*
1137           Correlation: measure of linear-dependencies in the image.
1138         */
1139         sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1140         sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1141         sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1142         if (image->colorspace == CMYKColorspace)
1143           sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black;
1144         if (image->alpha_trait == BlendPixelTrait)
1145           sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha;
1146         correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
1147         correlation.direction[i].green+=x*y*
1148           cooccurrence[x][y].direction[i].green;
1149         correlation.direction[i].blue+=x*y*
1150           cooccurrence[x][y].direction[i].blue;
1151         if (image->colorspace == CMYKColorspace)
1152           correlation.direction[i].black+=x*y*
1153             cooccurrence[x][y].direction[i].black;
1154         if (image->alpha_trait == BlendPixelTrait)
1155           correlation.direction[i].alpha+=x*y*
1156             cooccurrence[x][y].direction[i].alpha;
1157         /*
1158           Inverse Difference Moment.
1159         */
1160         channel_features[RedPixelChannel].inverse_difference_moment[i]+=
1161           cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
1162         channel_features[GreenPixelChannel].inverse_difference_moment[i]+=
1163           cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
1164         channel_features[BluePixelChannel].inverse_difference_moment[i]+=
1165           cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
1166         if (image->colorspace == CMYKColorspace)
1167           channel_features[BlackPixelChannel].inverse_difference_moment[i]+=
1168             cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1);
1169         if (image->alpha_trait == BlendPixelTrait)
1170           channel_features[AlphaPixelChannel].inverse_difference_moment[i]+=
1171             cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
1172         /*
1173           Sum average.
1174         */
1175         density_xy[y+x+2].direction[i].red+=
1176           cooccurrence[x][y].direction[i].red;
1177         density_xy[y+x+2].direction[i].green+=
1178           cooccurrence[x][y].direction[i].green;
1179         density_xy[y+x+2].direction[i].blue+=
1180           cooccurrence[x][y].direction[i].blue;
1181         if (image->colorspace == CMYKColorspace)
1182           density_xy[y+x+2].direction[i].black+=
1183             cooccurrence[x][y].direction[i].black;
1184         if (image->alpha_trait == BlendPixelTrait)
1185           density_xy[y+x+2].direction[i].alpha+=
1186             cooccurrence[x][y].direction[i].alpha;
1187         /*
1188           Entropy.
1189         */
1190         channel_features[RedPixelChannel].entropy[i]-=
1191           cooccurrence[x][y].direction[i].red*
1192           MagickLog10(cooccurrence[x][y].direction[i].red);
1193         channel_features[GreenPixelChannel].entropy[i]-=
1194           cooccurrence[x][y].direction[i].green*
1195           MagickLog10(cooccurrence[x][y].direction[i].green);
1196         channel_features[BluePixelChannel].entropy[i]-=
1197           cooccurrence[x][y].direction[i].blue*
1198           MagickLog10(cooccurrence[x][y].direction[i].blue);
1199         if (image->colorspace == CMYKColorspace)
1200           channel_features[BlackPixelChannel].entropy[i]-=
1201             cooccurrence[x][y].direction[i].black*
1202             MagickLog10(cooccurrence[x][y].direction[i].black);
1203         if (image->alpha_trait == BlendPixelTrait)
1204           channel_features[AlphaPixelChannel].entropy[i]-=
1205             cooccurrence[x][y].direction[i].alpha*
1206             MagickLog10(cooccurrence[x][y].direction[i].alpha);
1207         /*
1208           Information Measures of Correlation.
1209         */
1210         density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
1211         density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
1212         density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1213         if (image->alpha_trait == BlendPixelTrait)
1214           density_x[x].direction[i].alpha+=
1215             cooccurrence[x][y].direction[i].alpha;
1216         if (image->colorspace == CMYKColorspace)
1217           density_x[x].direction[i].black+=
1218             cooccurrence[x][y].direction[i].black;
1219         density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1220         density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1221         density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1222         if (image->colorspace == CMYKColorspace)
1223           density_y[y].direction[i].black+=
1224             cooccurrence[x][y].direction[i].black;
1225         if (image->alpha_trait == BlendPixelTrait)
1226           density_y[y].direction[i].alpha+=
1227             cooccurrence[x][y].direction[i].alpha;
1228       }
1229       mean.direction[i].red+=y*sum[y].direction[i].red;
1230       sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
1231       mean.direction[i].green+=y*sum[y].direction[i].green;
1232       sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
1233       mean.direction[i].blue+=y*sum[y].direction[i].blue;
1234       sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
1235       if (image->colorspace == CMYKColorspace)
1236         {
1237           mean.direction[i].black+=y*sum[y].direction[i].black;
1238           sum_squares.direction[i].black+=y*y*sum[y].direction[i].black;
1239         }
1240       if (image->alpha_trait == BlendPixelTrait)
1241         {
1242           mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
1243           sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
1244         }
1245     }
1246     /*
1247       Correlation: measure of linear-dependencies in the image.
1248     */
1249     channel_features[RedPixelChannel].correlation[i]=
1250       (correlation.direction[i].red-mean.direction[i].red*
1251       mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
1252       (mean.direction[i].red*mean.direction[i].red))*sqrt(
1253       sum_squares.direction[i].red-(mean.direction[i].red*
1254       mean.direction[i].red)));
1255     channel_features[GreenPixelChannel].correlation[i]=
1256       (correlation.direction[i].green-mean.direction[i].green*
1257       mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
1258       (mean.direction[i].green*mean.direction[i].green))*sqrt(
1259       sum_squares.direction[i].green-(mean.direction[i].green*
1260       mean.direction[i].green)));
1261     channel_features[BluePixelChannel].correlation[i]=
1262       (correlation.direction[i].blue-mean.direction[i].blue*
1263       mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
1264       (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
1265       sum_squares.direction[i].blue-(mean.direction[i].blue*
1266       mean.direction[i].blue)));
1267     if (image->colorspace == CMYKColorspace)
1268       channel_features[BlackPixelChannel].correlation[i]=
1269         (correlation.direction[i].black-mean.direction[i].black*
1270         mean.direction[i].black)/(sqrt(sum_squares.direction[i].black-
1271         (mean.direction[i].black*mean.direction[i].black))*sqrt(
1272         sum_squares.direction[i].black-(mean.direction[i].black*
1273         mean.direction[i].black)));
1274     if (image->alpha_trait == BlendPixelTrait)
1275       channel_features[AlphaPixelChannel].correlation[i]=
1276         (correlation.direction[i].alpha-mean.direction[i].alpha*
1277         mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
1278         (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt(
1279         sum_squares.direction[i].alpha-(mean.direction[i].alpha*
1280         mean.direction[i].alpha)));
1281   }
1282   /*
1283     Compute more texture features.
1284   */
1285 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1286   #pragma omp parallel for schedule(static,4) shared(status) \
1287     magick_threads(image,image,number_grays,1)
1288 #endif
1289   for (i=0; i < 4; i++)
1290   {
1291     register ssize_t
1292       x;
1293
1294     for (x=2; x < (ssize_t) (2*number_grays); x++)
1295     {
1296       /*
1297         Sum average.
1298       */
1299       channel_features[RedPixelChannel].sum_average[i]+=
1300         x*density_xy[x].direction[i].red;
1301       channel_features[GreenPixelChannel].sum_average[i]+=
1302         x*density_xy[x].direction[i].green;
1303       channel_features[BluePixelChannel].sum_average[i]+=
1304         x*density_xy[x].direction[i].blue;
1305       if (image->colorspace == CMYKColorspace)
1306         channel_features[BlackPixelChannel].sum_average[i]+=
1307           x*density_xy[x].direction[i].black;
1308       if (image->alpha_trait == BlendPixelTrait)
1309         channel_features[AlphaPixelChannel].sum_average[i]+=
1310           x*density_xy[x].direction[i].alpha;
1311       /*
1312         Sum entropy.
1313       */
1314       channel_features[RedPixelChannel].sum_entropy[i]-=
1315         density_xy[x].direction[i].red*
1316         MagickLog10(density_xy[x].direction[i].red);
1317       channel_features[GreenPixelChannel].sum_entropy[i]-=
1318         density_xy[x].direction[i].green*
1319         MagickLog10(density_xy[x].direction[i].green);
1320       channel_features[BluePixelChannel].sum_entropy[i]-=
1321         density_xy[x].direction[i].blue*
1322         MagickLog10(density_xy[x].direction[i].blue);
1323       if (image->colorspace == CMYKColorspace)
1324         channel_features[BlackPixelChannel].sum_entropy[i]-=
1325           density_xy[x].direction[i].black*
1326           MagickLog10(density_xy[x].direction[i].black);
1327       if (image->alpha_trait == BlendPixelTrait)
1328         channel_features[AlphaPixelChannel].sum_entropy[i]-=
1329           density_xy[x].direction[i].alpha*
1330           MagickLog10(density_xy[x].direction[i].alpha);
1331       /*
1332         Sum variance.
1333       */
1334       channel_features[RedPixelChannel].sum_variance[i]+=
1335         (x-channel_features[RedPixelChannel].sum_entropy[i])*
1336         (x-channel_features[RedPixelChannel].sum_entropy[i])*
1337         density_xy[x].direction[i].red;
1338       channel_features[GreenPixelChannel].sum_variance[i]+=
1339         (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1340         (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1341         density_xy[x].direction[i].green;
1342       channel_features[BluePixelChannel].sum_variance[i]+=
1343         (x-channel_features[BluePixelChannel].sum_entropy[i])*
1344         (x-channel_features[BluePixelChannel].sum_entropy[i])*
1345         density_xy[x].direction[i].blue;
1346       if (image->colorspace == CMYKColorspace)
1347         channel_features[BlackPixelChannel].sum_variance[i]+=
1348           (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1349           (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1350           density_xy[x].direction[i].black;
1351       if (image->alpha_trait == BlendPixelTrait)
1352         channel_features[AlphaPixelChannel].sum_variance[i]+=
1353           (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1354           (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1355           density_xy[x].direction[i].alpha;
1356     }
1357   }
1358   /*
1359     Compute more texture features.
1360   */
1361 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1362   #pragma omp parallel for schedule(static,4) shared(status) \
1363     magick_threads(image,image,number_grays,1)
1364 #endif
1365   for (i=0; i < 4; i++)
1366   {
1367     register ssize_t
1368       y;
1369
1370     for (y=0; y < (ssize_t) number_grays; y++)
1371     {
1372       register ssize_t
1373         x;
1374
1375       for (x=0; x < (ssize_t) number_grays; x++)
1376       {
1377         /*
1378           Sum of Squares: Variance
1379         */
1380         variance.direction[i].red+=(y-mean.direction[i].red+1)*
1381           (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
1382         variance.direction[i].green+=(y-mean.direction[i].green+1)*
1383           (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
1384         variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
1385           (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
1386         if (image->colorspace == CMYKColorspace)
1387           variance.direction[i].black+=(y-mean.direction[i].black+1)*
1388             (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black;
1389         if (image->alpha_trait == BlendPixelTrait)
1390           variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
1391             (y-mean.direction[i].alpha+1)*
1392             cooccurrence[x][y].direction[i].alpha;
1393         /*
1394           Sum average / Difference Variance.
1395         */
1396         density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
1397           cooccurrence[x][y].direction[i].red;
1398         density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
1399           cooccurrence[x][y].direction[i].green;
1400         density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
1401           cooccurrence[x][y].direction[i].blue;
1402         if (image->colorspace == CMYKColorspace)
1403           density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
1404             cooccurrence[x][y].direction[i].black;
1405         if (image->alpha_trait == BlendPixelTrait)
1406           density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
1407             cooccurrence[x][y].direction[i].alpha;
1408         /*
1409           Information Measures of Correlation.
1410         */
1411         entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
1412           MagickLog10(cooccurrence[x][y].direction[i].red);
1413         entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
1414           MagickLog10(cooccurrence[x][y].direction[i].green);
1415         entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
1416           MagickLog10(cooccurrence[x][y].direction[i].blue);
1417         if (image->colorspace == CMYKColorspace)
1418           entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black*
1419             MagickLog10(cooccurrence[x][y].direction[i].black);
1420         if (image->alpha_trait == BlendPixelTrait)
1421           entropy_xy.direction[i].alpha-=
1422             cooccurrence[x][y].direction[i].alpha*MagickLog10(
1423             cooccurrence[x][y].direction[i].alpha);
1424         entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
1425           MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red));
1426         entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
1427           MagickLog10(density_x[x].direction[i].green*
1428           density_y[y].direction[i].green));
1429         entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
1430           MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue));
1431         if (image->colorspace == CMYKColorspace)
1432           entropy_xy1.direction[i].black-=(
1433             cooccurrence[x][y].direction[i].black*MagickLog10(
1434             density_x[x].direction[i].black*density_y[y].direction[i].black));
1435         if (image->alpha_trait == BlendPixelTrait)
1436           entropy_xy1.direction[i].alpha-=(
1437             cooccurrence[x][y].direction[i].alpha*MagickLog10(
1438             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1439         entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
1440           density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red*
1441           density_y[y].direction[i].red));
1442         entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
1443           density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green*
1444           density_y[y].direction[i].green));
1445         entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
1446           density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue*
1447           density_y[y].direction[i].blue));
1448         if (image->colorspace == CMYKColorspace)
1449           entropy_xy2.direction[i].black-=(density_x[x].direction[i].black*
1450             density_y[y].direction[i].black*MagickLog10(
1451             density_x[x].direction[i].black*density_y[y].direction[i].black));
1452         if (image->alpha_trait == BlendPixelTrait)
1453           entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha*
1454             density_y[y].direction[i].alpha*MagickLog10(
1455             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1456       }
1457     }
1458     channel_features[RedPixelChannel].variance_sum_of_squares[i]=
1459       variance.direction[i].red;
1460     channel_features[GreenPixelChannel].variance_sum_of_squares[i]=
1461       variance.direction[i].green;
1462     channel_features[BluePixelChannel].variance_sum_of_squares[i]=
1463       variance.direction[i].blue;
1464     if (image->colorspace == CMYKColorspace)
1465       channel_features[BlackPixelChannel].variance_sum_of_squares[i]=
1466         variance.direction[i].black;
1467     if (image->alpha_trait == BlendPixelTrait)
1468       channel_features[AlphaPixelChannel].variance_sum_of_squares[i]=
1469         variance.direction[i].alpha;
1470   }
1471   /*
1472     Compute more texture features.
1473   */
1474   (void) ResetMagickMemory(&variance,0,sizeof(variance));
1475   (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
1476 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1477   #pragma omp parallel for schedule(static,4) shared(status) \
1478     magick_threads(image,image,number_grays,1)
1479 #endif
1480   for (i=0; i < 4; i++)
1481   {
1482     register ssize_t
1483       x;
1484
1485     for (x=0; x < (ssize_t) number_grays; x++)
1486     {
1487       /*
1488         Difference variance.
1489       */
1490       variance.direction[i].red+=density_xy[x].direction[i].red;
1491       variance.direction[i].green+=density_xy[x].direction[i].green;
1492       variance.direction[i].blue+=density_xy[x].direction[i].blue;
1493       if (image->colorspace == CMYKColorspace)
1494         variance.direction[i].black+=density_xy[x].direction[i].black;
1495       if (image->alpha_trait == BlendPixelTrait)
1496         variance.direction[i].alpha+=density_xy[x].direction[i].alpha;
1497       sum_squares.direction[i].red+=density_xy[x].direction[i].red*
1498         density_xy[x].direction[i].red;
1499       sum_squares.direction[i].green+=density_xy[x].direction[i].green*
1500         density_xy[x].direction[i].green;
1501       sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
1502         density_xy[x].direction[i].blue;
1503       if (image->colorspace == CMYKColorspace)
1504         sum_squares.direction[i].black+=density_xy[x].direction[i].black*
1505           density_xy[x].direction[i].black;
1506       if (image->alpha_trait == BlendPixelTrait)
1507         sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
1508           density_xy[x].direction[i].alpha;
1509       /*
1510         Difference entropy.
1511       */
1512       channel_features[RedPixelChannel].difference_entropy[i]-=
1513         density_xy[x].direction[i].red*
1514         MagickLog10(density_xy[x].direction[i].red);
1515       channel_features[GreenPixelChannel].difference_entropy[i]-=
1516         density_xy[x].direction[i].green*
1517         MagickLog10(density_xy[x].direction[i].green);
1518       channel_features[BluePixelChannel].difference_entropy[i]-=
1519         density_xy[x].direction[i].blue*
1520         MagickLog10(density_xy[x].direction[i].blue);
1521       if (image->colorspace == CMYKColorspace)
1522         channel_features[BlackPixelChannel].difference_entropy[i]-=
1523           density_xy[x].direction[i].black*
1524           MagickLog10(density_xy[x].direction[i].black);
1525       if (image->alpha_trait == BlendPixelTrait)
1526         channel_features[AlphaPixelChannel].difference_entropy[i]-=
1527           density_xy[x].direction[i].alpha*
1528           MagickLog10(density_xy[x].direction[i].alpha);
1529       /*
1530         Information Measures of Correlation.
1531       */
1532       entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1533         MagickLog10(density_x[x].direction[i].red));
1534       entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1535         MagickLog10(density_x[x].direction[i].green));
1536       entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1537         MagickLog10(density_x[x].direction[i].blue));
1538       if (image->colorspace == CMYKColorspace)
1539         entropy_x.direction[i].black-=(density_x[x].direction[i].black*
1540           MagickLog10(density_x[x].direction[i].black));
1541       if (image->alpha_trait == BlendPixelTrait)
1542         entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha*
1543           MagickLog10(density_x[x].direction[i].alpha));
1544       entropy_y.direction[i].red-=(density_y[x].direction[i].red*
1545         MagickLog10(density_y[x].direction[i].red));
1546       entropy_y.direction[i].green-=(density_y[x].direction[i].green*
1547         MagickLog10(density_y[x].direction[i].green));
1548       entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
1549         MagickLog10(density_y[x].direction[i].blue));
1550       if (image->colorspace == CMYKColorspace)
1551         entropy_y.direction[i].black-=(density_y[x].direction[i].black*
1552           MagickLog10(density_y[x].direction[i].black));
1553       if (image->alpha_trait == BlendPixelTrait)
1554         entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
1555           MagickLog10(density_y[x].direction[i].alpha));
1556     }
1557     /*
1558       Difference variance.
1559     */
1560     channel_features[RedPixelChannel].difference_variance[i]=
1561       (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1562       (variance.direction[i].red*variance.direction[i].red))/
1563       ((double) number_grays*number_grays*number_grays*number_grays);
1564     channel_features[GreenPixelChannel].difference_variance[i]=
1565       (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1566       (variance.direction[i].green*variance.direction[i].green))/
1567       ((double) number_grays*number_grays*number_grays*number_grays);
1568     channel_features[BluePixelChannel].difference_variance[i]=
1569       (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1570       (variance.direction[i].blue*variance.direction[i].blue))/
1571       ((double) number_grays*number_grays*number_grays*number_grays);
1572     if (image->colorspace == CMYKColorspace)
1573       channel_features[BlackPixelChannel].difference_variance[i]=
1574         (((double) number_grays*number_grays*sum_squares.direction[i].black)-
1575         (variance.direction[i].black*variance.direction[i].black))/
1576         ((double) number_grays*number_grays*number_grays*number_grays);
1577     if (image->alpha_trait == BlendPixelTrait)
1578       channel_features[AlphaPixelChannel].difference_variance[i]=
1579         (((double) number_grays*number_grays*sum_squares.direction[i].alpha)-
1580         (variance.direction[i].alpha*variance.direction[i].alpha))/
1581         ((double) number_grays*number_grays*number_grays*number_grays);
1582     /*
1583       Information Measures of Correlation.
1584     */
1585     channel_features[RedPixelChannel].measure_of_correlation_1[i]=
1586       (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1587       (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1588        entropy_x.direction[i].red : entropy_y.direction[i].red);
1589     channel_features[GreenPixelChannel].measure_of_correlation_1[i]=
1590       (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1591       (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1592        entropy_x.direction[i].green : entropy_y.direction[i].green);
1593     channel_features[BluePixelChannel].measure_of_correlation_1[i]=
1594       (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1595       (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1596        entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1597     if (image->colorspace == CMYKColorspace)
1598       channel_features[BlackPixelChannel].measure_of_correlation_1[i]=
1599         (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/
1600         (entropy_x.direction[i].black > entropy_y.direction[i].black ?
1601          entropy_x.direction[i].black : entropy_y.direction[i].black);
1602     if (image->alpha_trait == BlendPixelTrait)
1603       channel_features[AlphaPixelChannel].measure_of_correlation_1[i]=
1604         (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/
1605         (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ?
1606          entropy_x.direction[i].alpha : entropy_y.direction[i].alpha);
1607     channel_features[RedPixelChannel].measure_of_correlation_2[i]=
1608       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].red-
1609       entropy_xy.direction[i].red)))));
1610     channel_features[GreenPixelChannel].measure_of_correlation_2[i]=
1611       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].green-
1612       entropy_xy.direction[i].green)))));
1613     channel_features[BluePixelChannel].measure_of_correlation_2[i]=
1614       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].blue-
1615       entropy_xy.direction[i].blue)))));
1616     if (image->colorspace == CMYKColorspace)
1617       channel_features[BlackPixelChannel].measure_of_correlation_2[i]=
1618         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].black-
1619         entropy_xy.direction[i].black)))));
1620     if (image->alpha_trait == BlendPixelTrait)
1621       channel_features[AlphaPixelChannel].measure_of_correlation_2[i]=
1622         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].alpha-
1623         entropy_xy.direction[i].alpha)))));
1624   }
1625   /*
1626     Compute more texture features.
1627   */
1628 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1629   #pragma omp parallel for schedule(static,4) shared(status) \
1630     magick_threads(image,image,number_grays,1)
1631 #endif
1632   for (i=0; i < 4; i++)
1633   {
1634     ssize_t
1635       z;
1636
1637     for (z=0; z < (ssize_t) number_grays; z++)
1638     {
1639       register ssize_t
1640         y;
1641
1642       ChannelStatistics
1643         pixel;
1644
1645       (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
1646       for (y=0; y < (ssize_t) number_grays; y++)
1647       {
1648         register ssize_t
1649           x;
1650
1651         for (x=0; x < (ssize_t) number_grays; x++)
1652         {
1653           /*
1654             Contrast:  amount of local variations present in an image.
1655           */
1656           if (((y-x) == z) || ((x-y) == z))
1657             {
1658               pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1659               pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1660               pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1661               if (image->colorspace == CMYKColorspace)
1662                 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
1663               if (image->alpha_trait == BlendPixelTrait)
1664                 pixel.direction[i].alpha+=
1665                   cooccurrence[x][y].direction[i].alpha;
1666             }
1667           /*
1668             Maximum Correlation Coefficient.
1669           */
1670           Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1671             cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1672             density_y[x].direction[i].red;
1673           Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1674             cooccurrence[y][x].direction[i].green/
1675             density_x[z].direction[i].green/density_y[x].direction[i].red;
1676           Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1677             cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/
1678             density_y[x].direction[i].blue;
1679           if (image->colorspace == CMYKColorspace)
1680             Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
1681               cooccurrence[y][x].direction[i].black/
1682               density_x[z].direction[i].black/density_y[x].direction[i].black;
1683           if (image->alpha_trait == BlendPixelTrait)
1684             Q[z][y].direction[i].alpha+=
1685               cooccurrence[z][x].direction[i].alpha*
1686               cooccurrence[y][x].direction[i].alpha/
1687               density_x[z].direction[i].alpha/
1688               density_y[x].direction[i].alpha;
1689         }
1690       }
1691       channel_features[RedPixelChannel].contrast[i]+=z*z*
1692         pixel.direction[i].red;
1693       channel_features[GreenPixelChannel].contrast[i]+=z*z*
1694         pixel.direction[i].green;
1695       channel_features[BluePixelChannel].contrast[i]+=z*z*
1696         pixel.direction[i].blue;
1697       if (image->colorspace == CMYKColorspace)
1698         channel_features[BlackPixelChannel].contrast[i]+=z*z*
1699           pixel.direction[i].black;
1700       if (image->alpha_trait == BlendPixelTrait)
1701         channel_features[AlphaPixelChannel].contrast[i]+=z*z*
1702           pixel.direction[i].alpha;
1703     }
1704     /*
1705       Maximum Correlation Coefficient.
1706       Future: return second largest eigenvalue of Q.
1707     */
1708     channel_features[RedPixelChannel].maximum_correlation_coefficient[i]=
1709       sqrt((double) -1.0);
1710     channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]=
1711       sqrt((double) -1.0);
1712     channel_features[BluePixelChannel].maximum_correlation_coefficient[i]=
1713       sqrt((double) -1.0);
1714     if (image->colorspace == CMYKColorspace)
1715       channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]=
1716         sqrt((double) -1.0);
1717     if (image->alpha_trait == BlendPixelTrait)
1718       channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]=
1719         sqrt((double) -1.0);
1720   }
1721   /*
1722     Relinquish resources.
1723   */
1724   sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1725   for (i=0; i < (ssize_t) number_grays; i++)
1726     Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1727   Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1728   density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1729   density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1730   density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1731   for (i=0; i < (ssize_t) number_grays; i++)
1732     cooccurrence[i]=(ChannelStatistics *)
1733       RelinquishMagickMemory(cooccurrence[i]);
1734   cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1735   return(channel_features);
1736 }