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