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

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