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