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