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 %                                John Cristy                                  %
  17 %                                 July 1992                                   %
  18 %                                                                             %
  19 %                                                                             %
  20 %  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
  21 %  dedicated to making software imaging solutions freely available.           %
  22 %                                                                             %
  23 %  You may not use this file except in compliance with the License.  You may  %
  24 %  obtain a copy of the License at                                            %
  25 %                                                                             %
  26 %    http://www.imagemagick.org/script/license.php                            %
  27 %                                                                             %
  28 %  Unless required by applicable law or agreed to in writing, software        %
  29 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  30 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  31 %  See the License for the specific language governing permissions and        %
  32 %  limitations under the License.                                             %
  33 %                                                                             %
  34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  35 %
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "MagickCore/studio.h"
  44 #include "MagickCore/property.h"
  45 #include "MagickCore/animate.h"
  46 #include "MagickCore/blob.h"
  47 #include "MagickCore/blob-private.h"
  48 #include "MagickCore/cache.h"
  49 #include "MagickCore/cache-private.h"
  50 #include "MagickCore/cache-view.h"
  51 #include "MagickCore/client.h"
  52 #include "MagickCore/color.h"
  53 #include "MagickCore/color-private.h"
  54 #include "MagickCore/colorspace.h"
  55 #include "MagickCore/colorspace-private.h"
  56 #include "MagickCore/composite.h"
  57 #include "MagickCore/composite-private.h"
  58 #include "MagickCore/compress.h"
  59 #include "MagickCore/constitute.h"
  60 #include "MagickCore/display.h"
  61 #include "MagickCore/draw.h"
  62 #include "MagickCore/enhance.h"
  63 #include "MagickCore/exception.h"
  64 #include "MagickCore/exception-private.h"
  65 #include "MagickCore/feature.h"
  66 #include "MagickCore/gem.h"
  67 #include "MagickCore/geometry.h"
  68 #include "MagickCore/list.h"
  69 #include "MagickCore/image-private.h"
  70 #include "MagickCore/magic.h"
  71 #include "MagickCore/magick.h"
  72 #include "MagickCore/memory_.h"
  73 #include "MagickCore/module.h"
  74 #include "MagickCore/monitor.h"
  75 #include "MagickCore/monitor-private.h"
  76 #include "MagickCore/option.h"
  77 #include "MagickCore/paint.h"
  78 #include "MagickCore/pixel-accessor.h"
  79 #include "MagickCore/profile.h"
  80 #include "MagickCore/quantize.h"
  81 #include "MagickCore/quantum-private.h"
  82 #include "MagickCore/random_.h"
  83 #include "MagickCore/segment.h"
  84 #include "MagickCore/semaphore.h"
  85 #include "MagickCore/signature-private.h"
  86 #include "MagickCore/string_.h"
  87 #include "MagickCore/thread-private.h"
  88 #include "MagickCore/timer.h"
  89 #include "MagickCore/utility.h"
  90 #include "MagickCore/version.h"
  91 \f
  92 /*
  93 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  94 %                                                                             %
  95 %                                                                             %
  96 %                                                                             %
  97 %   G e t I m a g e F e a t u r e s                                           %
  98 %                                                                             %
  99 %                                                                             %
 100 %                                                                             %
 101 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 102 %
 103 %  GetImageFeatures() returns features for each channel in the image in
 104 %  each of four directions (horizontal, vertical, left and right diagonals)
 105 %  for the specified distance.  The features include the angular second
 106 %  moment, contrast, correlation, sum of squares: variance, inverse difference
 107 %  moment, sum average, sum varience, sum entropy, entropy, difference variance,%  difference entropy, information measures of correlation 1, information
 108 %  measures of correlation 2, and maximum correlation coefficient.  You can
 109 %  access the red channel contrast, for example, like this:
 110 %
 111 %      channel_features=GetImageFeatures(image,1,exception);
 112 %      contrast=channel_features[RedChannel].contrast[0];
 113 %
 114 %  Use MagickRelinquishMemory() to free the features buffer.
 115 %
 116 %  The format of the GetImageFeatures method is:
 117 %
 118 %      ChannelFeatures *GetImageFeatures(const Image *image,
 119 %        const size_t distance,ExceptionInfo *exception)
 120 %
 121 %  A description of each parameter follows:
 122 %
 123 %    o image: the image.
 124 %
 125 %    o distance: the distance.
 126 %
 127 %    o exception: return any errors or warnings in this structure.
 128 %
 129 */
 130
 131 static inline ssize_t MagickAbsoluteValue(const ssize_t x)
 132 {
 133   if (x < 0)
 134     return(-x);
 135   return(x);
 136 }
 137
 138 MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
 139   const size_t distance,ExceptionInfo *exception)
 140 {
 141   typedef struct _ChannelStatistics
 142   {
 143     PixelInfo
 144       direction[4];  /* horizontal, vertical, left and right diagonals */
 145   } ChannelStatistics;
 146
 147   CacheView
 148     *image_view;
 149
 150   ChannelFeatures
 151     *channel_features;
 152
 153   ChannelStatistics
 154     **cooccurrence,
 155     correlation,
 156     *density_x,
 157     *density_xy,
 158     *density_y,
 159     entropy_x,
 160     entropy_xy,
 161     entropy_xy1,
 162     entropy_xy2,
 163     entropy_y,
 164     mean,
 165     **Q,
 166     *sum,
 167     sum_squares,
 168     variance;
 169
 170   PixelPacket
 171     gray,
 172     *grays;
 173
 174   MagickBooleanType
 175     status;
 176
 177   register ssize_t
 178     i;
 179
 180   size_t
 181     length;
 182
 183   ssize_t
 184     y,
 185     z;
 186
 187   unsigned int
 188     number_grays;
 189
 190   assert(image != (Image *) NULL);
 191   assert(image->signature == MagickSignature);
 192   if (image->debug != MagickFalse)
 193     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 194   if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
 195     return((ChannelFeatures *) NULL);
 196   length=CompositeChannels+1UL;
 197   channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
 198     sizeof(*channel_features));
 199   if (channel_features == (ChannelFeatures *) NULL)
 200     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 201   (void) ResetMagickMemory(channel_features,0,length*
 202     sizeof(*channel_features));
 203   /*
 204     Form grays.
 205   */
 206   grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays));
 207   if (grays == (PixelPacket *) NULL)
 208     {
 209       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 210         channel_features);
 211       (void) ThrowMagickException(exception,GetMagickModule(),
 212         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 213       return(channel_features);
 214     }
 215   for (i=0; i <= (ssize_t) MaxMap; i++)
 216   {
 217     grays[i].red=(~0U);
 218     grays[i].green=(~0U);
 219     grays[i].blue=(~0U);
 220     grays[i].alpha=(~0U);
 221     grays[i].black=(~0U);
 222   }
 223   status=MagickTrue;
 224   image_view=AcquireCacheView(image);
 225 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 226   #pragma omp parallel for schedule(dynamic,4) shared(status)
 227 #endif
 228   for (y=0; y < (ssize_t) image->rows; y++)
 229   {
 230     register const Quantum
 231       *restrict p;
 232
 233     register ssize_t
 234       x;
 235
 236     if (status == MagickFalse)
 237       continue;
 238     p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
 239     if (p == (const Quantum *) NULL)
 240       {
 241         status=MagickFalse;
 242         continue;
 243       }
 244     for (x=0; x < (ssize_t) image->columns; x++)
 245     {
 246       grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
 247         ScaleQuantumToMap(GetPixelRed(image,p));
 248       grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
 249         ScaleQuantumToMap(GetPixelGreen(image,p));
 250       grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
 251         ScaleQuantumToMap(GetPixelBlue(image,p));
 252       if (image->colorspace == CMYKColorspace)
 253         grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
 254           ScaleQuantumToMap(GetPixelBlack(image,p));
 255       if (image->matte != MagickFalse)
 256         grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
 257           ScaleQuantumToMap(GetPixelAlpha(image,p));
 258       p+=GetPixelChannels(image);
 259     }
 260   }
 261   image_view=DestroyCacheView(image_view);
 262   if (status == MagickFalse)
 263     {
 264       grays=(PixelPacket *) RelinquishMagickMemory(grays);
 265       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 266         channel_features);
 267       return(channel_features);
 268     }
 269   (void) ResetMagickMemory(&gray,0,sizeof(gray));
 270   for (i=0; i <= (ssize_t) MaxMap; i++)
 271   {
 272     if (grays[i].red != ~0U)
 273       grays[gray.red++].red=grays[i].red;
 274     if (grays[i].green != ~0U)
 275       grays[gray.green++].green=grays[i].green;
 276     if (grays[i].blue != ~0U)
 277       grays[gray.blue++].blue=grays[i].blue;
 278     if (image->colorspace == CMYKColorspace)
 279       if (grays[i].black != ~0U)
 280         grays[gray.black++].black=grays[i].black;
 281     if (image->matte != MagickFalse)
 282       if (grays[i].alpha != ~0U)
 283         grays[gray.alpha++].alpha=grays[i].alpha;
 284   }
 285   /*
 286     Allocate spatial dependence matrix.
 287   */
 288   number_grays=gray.red;
 289   if (gray.green > number_grays)
 290     number_grays=gray.green;
 291   if (gray.blue > number_grays)
 292     number_grays=gray.blue;
 293   if (image->colorspace == CMYKColorspace)
 294     if (gray.black > number_grays)
 295       number_grays=gray.black;
 296   if (image->matte != MagickFalse)
 297     if (gray.alpha > number_grays)
 298       number_grays=gray.alpha;
 299   cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays,
 300     sizeof(*cooccurrence));
 301   density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
 302     sizeof(*density_x));
 303   density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
 304     sizeof(*density_xy));
 305   density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
 306     sizeof(*density_y));
 307   Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q));
 308   sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum));
 309   if ((cooccurrence == (ChannelStatistics **) NULL) ||
 310       (density_x == (ChannelStatistics *) NULL) ||
 311       (density_xy == (ChannelStatistics *) NULL) ||
 312       (density_y == (ChannelStatistics *) NULL) ||
 313       (Q == (ChannelStatistics **) NULL) ||
 314       (sum == (ChannelStatistics *) NULL))
 315     {
 316       if (Q != (ChannelStatistics **) NULL)
 317         {
 318           for (i=0; i < (ssize_t) number_grays; i++)
 319             Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
 320           Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
 321         }
 322       if (sum != (ChannelStatistics *) NULL)
 323         sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
 324       if (density_y != (ChannelStatistics *) NULL)
 325         density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
 326       if (density_xy != (ChannelStatistics *) NULL)
 327         density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
 328       if (density_x != (ChannelStatistics *) NULL)
 329         density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
 330       if (cooccurrence != (ChannelStatistics **) NULL)
 331         {
 332           for (i=0; i < (ssize_t) number_grays; i++)
 333             cooccurrence[i]=(ChannelStatistics *)
 334               RelinquishMagickMemory(cooccurrence[i]);
 335           cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
 336             cooccurrence);
 337         }
 338       grays=(PixelPacket *) RelinquishMagickMemory(grays);
 339       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 340         channel_features);
 341       (void) ThrowMagickException(exception,GetMagickModule(),
 342         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 343       return(channel_features);
 344     }
 345   (void) ResetMagickMemory(&correlation,0,sizeof(correlation));
 346   (void) ResetMagickMemory(density_x,0,2*(number_grays+1)*sizeof(*density_x));
 347   (void) ResetMagickMemory(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
 348   (void) ResetMagickMemory(density_y,0,2*(number_grays+1)*sizeof(*density_y));
 349   (void) ResetMagickMemory(&mean,0,sizeof(mean));
 350   (void) ResetMagickMemory(sum,0,number_grays*sizeof(*sum));
 351   (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
 352   (void) ResetMagickMemory(density_xy,0,2*number_grays*sizeof(*density_xy));
 353   (void) ResetMagickMemory(&entropy_x,0,sizeof(entropy_x));
 354   (void) ResetMagickMemory(&entropy_xy,0,sizeof(entropy_xy));
 355   (void) ResetMagickMemory(&entropy_xy1,0,sizeof(entropy_xy1));
 356   (void) ResetMagickMemory(&entropy_xy2,0,sizeof(entropy_xy2));
 357   (void) ResetMagickMemory(&entropy_y,0,sizeof(entropy_y));
 358   (void) ResetMagickMemory(&variance,0,sizeof(variance));
 359   for (i=0; i < (ssize_t) number_grays; i++)
 360   {
 361     cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
 362       sizeof(**cooccurrence));
 363     Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
 364     if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
 365         (Q[i] == (ChannelStatistics *) NULL))
 366       break;
 367     (void) ResetMagickMemory(cooccurrence[i],0,number_grays*
 368       sizeof(**cooccurrence));
 369     (void) ResetMagickMemory(Q[i],0,number_grays*sizeof(**Q));
 370   }
 371   if (i < (ssize_t) number_grays)
 372     {
 373       for (i--; i >= 0; i--)
 374       {
 375         if (Q[i] != (ChannelStatistics *) NULL)
 376           Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
 377         if (cooccurrence[i] != (ChannelStatistics *) NULL)
 378           cooccurrence[i]=(ChannelStatistics *)
 379             RelinquishMagickMemory(cooccurrence[i]);
 380       }
 381       Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
 382       cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
 383       sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
 384       density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
 385       density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
 386       density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
 387       grays=(PixelPacket *) RelinquishMagickMemory(grays);
 388       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 389         channel_features);
 390       (void) ThrowMagickException(exception,GetMagickModule(),
 391         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 392       return(channel_features);
 393     }
 394   /*
 395     Initialize spatial dependence matrix.
 396   */
 397   status=MagickTrue;
 398   image_view=AcquireCacheView(image);
 399   for (y=0; y < (ssize_t) image->rows; y++)
 400   {
 401     register const Quantum
 402       *restrict p;
 403
 404     register ssize_t
 405       x;
 406
 407     ssize_t
 408       i,
 409       offset,
 410       u,
 411       v;
 412
 413     if (status == MagickFalse)
 414       continue;
 415     p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,y,image->columns+
 416       2*distance,distance+1,exception);
 417     if (p == (const Quantum *) NULL)
 418       {
 419         status=MagickFalse;
 420         continue;
 421       }
 422     p+=distance*GetPixelChannels(image);;
 423     for (x=0; x < (ssize_t) image->columns; x++)
 424     {
 425       for (i=0; i < 4; i++)
 426       {
 427         switch (i)
 428         {
 429           case 0:
 430           default:
 431           {
 432             /*
 433               Horizontal adjacency.
 434             */
 435             offset=(ssize_t) distance;
 436             break;
 437           }
 438           case 1:
 439           {
 440             /*
 441               Vertical adjacency.
 442             */
 443             offset=(ssize_t) (image->columns+2*distance);
 444             break;
 445           }
 446           case 2:
 447           {
 448             /*
 449               Right diagonal adjacency.
 450             */
 451             offset=(ssize_t) ((image->columns+2*distance)-distance);
 452             break;
 453           }
 454           case 3:
 455           {
 456             /*
 457               Left diagonal adjacency.
 458             */
 459             offset=(ssize_t) ((image->columns+2*distance)+distance);
 460             break;
 461           }
 462         }
 463         u=0;
 464         v=0;
 465         while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p)))
 466           u++;
 467         while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image))))
 468           v++;
 469         cooccurrence[u][v].direction[i].red++;
 470         cooccurrence[v][u].direction[i].red++;
 471         u=0;
 472         v=0;
 473         while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p)))
 474           u++;
 475         while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image))))
 476           v++;
 477         cooccurrence[u][v].direction[i].green++;
 478         cooccurrence[v][u].direction[i].green++;
 479         u=0;
 480         v=0;
 481         while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p)))
 482           u++;
 483         while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image))))
 484           v++;
 485         cooccurrence[u][v].direction[i].blue++;
 486         cooccurrence[v][u].direction[i].blue++;
 487         if (image->colorspace == CMYKColorspace)
 488           {
 489             u=0;
 490             v=0;
 491             while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p)))
 492               u++;
 493             while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image))))
 494               v++;
 495             cooccurrence[u][v].direction[i].black++;
 496             cooccurrence[v][u].direction[i].black++;
 497           }
 498         if (image->matte != MagickFalse)
 499           {
 500             u=0;
 501             v=0;
 502             while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p)))
 503               u++;
 504             while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image))))
 505               v++;
 506             cooccurrence[u][v].direction[i].alpha++;
 507             cooccurrence[v][u].direction[i].alpha++;
 508           }
 509       }
 510       p+=GetPixelChannels(image);
 511     }
 512   }
 513   grays=(PixelPacket *) RelinquishMagickMemory(grays);
 514   image_view=DestroyCacheView(image_view);
 515   if (status == MagickFalse)
 516     {
 517       for (i=0; i < (ssize_t) number_grays; i++)
 518         cooccurrence[i]=(ChannelStatistics *)
 519           RelinquishMagickMemory(cooccurrence[i]);
 520       cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
 521       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 522         channel_features);
 523       (void) ThrowMagickException(exception,GetMagickModule(),
 524         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 525       return(channel_features);
 526     }
 527   /*
 528     Normalize spatial dependence matrix.
 529   */
 530   for (i=0; i < 4; i++)
 531   {
 532     double
 533       normalize;
 534
 535     register ssize_t
 536       y;
 537
 538     switch (i)
 539     {
 540       case 0:
 541       default:
 542       {
 543         /*
 544           Horizontal adjacency.
 545         */
 546         normalize=2.0*image->rows*(image->columns-distance);
 547         break;
 548       }
 549       case 1:
 550       {
 551         /*
 552           Vertical adjacency.
 553         */
 554         normalize=2.0*(image->rows-distance)*image->columns;
 555         break;
 556       }
 557       case 2:
 558       {
 559         /*
 560           Right diagonal adjacency.
 561         */
 562         normalize=2.0*(image->rows-distance)*(image->columns-distance);
 563         break;
 564       }
 565       case 3:
 566       {
 567         /*
 568           Left diagonal adjacency.
 569         */
 570         normalize=2.0*(image->rows-distance)*(image->columns-distance);
 571         break;
 572       }
 573     }
 574     normalize=1.0/(fabs((double) normalize) <= MagickEpsilon ? 1.0 : normalize);
 575     for (y=0; y < (ssize_t) number_grays; y++)
 576     {
 577       register ssize_t
 578         x;
 579
 580       for (x=0; x < (ssize_t) number_grays; x++)
 581       {
 582         cooccurrence[x][y].direction[i].red*=normalize;
 583         cooccurrence[x][y].direction[i].green*=normalize;
 584         cooccurrence[x][y].direction[i].blue*=normalize;
 585         if (image->colorspace == CMYKColorspace)
 586           cooccurrence[x][y].direction[i].black*=normalize;
 587         if (image->matte != MagickFalse)
 588           cooccurrence[x][y].direction[i].alpha*=normalize;
 589       }
 590     }
 591   }
 592   /*
 593     Compute texture features.
 594   */
 595 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 596   #pragma omp parallel for schedule(dynamic,4) shared(status)
 597 #endif
 598   for (i=0; i < 4; i++)
 599   {
 600     register ssize_t
 601       y;
 602
 603     for (y=0; y < (ssize_t) number_grays; y++)
 604     {
 605       register ssize_t
 606         x;
 607
 608       for (x=0; x < (ssize_t) number_grays; x++)
 609       {
 610         /*
 611           Angular second moment:  measure of homogeneity of the image.
 612         */
 613         channel_features[RedChannel].angular_second_moment[i]+=
 614           cooccurrence[x][y].direction[i].red*
 615           cooccurrence[x][y].direction[i].red;
 616         channel_features[GreenChannel].angular_second_moment[i]+=
 617           cooccurrence[x][y].direction[i].green*
 618           cooccurrence[x][y].direction[i].green;
 619         channel_features[BlueChannel].angular_second_moment[i]+=
 620           cooccurrence[x][y].direction[i].blue*
 621           cooccurrence[x][y].direction[i].blue;
 622         if (image->colorspace == CMYKColorspace)
 623           channel_features[BlackChannel].angular_second_moment[i]+=
 624             cooccurrence[x][y].direction[i].black*
 625             cooccurrence[x][y].direction[i].black;
 626         if (image->matte != MagickFalse)
 627           channel_features[OpacityChannel].angular_second_moment[i]+=
 628             cooccurrence[x][y].direction[i].alpha*
 629             cooccurrence[x][y].direction[i].alpha;
 630         /*
 631           Correlation: measure of linear-dependencies in the image.
 632         */
 633         sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
 634         sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
 635         sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
 636         if (image->colorspace == CMYKColorspace)
 637           sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black;
 638         if (image->matte != MagickFalse)
 639           sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha;
 640         correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
 641         correlation.direction[i].green+=x*y*
 642           cooccurrence[x][y].direction[i].green;
 643         correlation.direction[i].blue+=x*y*
 644           cooccurrence[x][y].direction[i].blue;
 645         if (image->colorspace == CMYKColorspace)
 646           correlation.direction[i].black+=x*y*
 647             cooccurrence[x][y].direction[i].black;
 648         if (image->matte != MagickFalse)
 649           correlation.direction[i].alpha+=x*y*
 650             cooccurrence[x][y].direction[i].alpha;
 651         /*
 652           Inverse Difference Moment.
 653         */
 654         channel_features[RedChannel].inverse_difference_moment[i]+=
 655           cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
 656         channel_features[GreenChannel].inverse_difference_moment[i]+=
 657           cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
 658         channel_features[BlueChannel].inverse_difference_moment[i]+=
 659           cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
 660         if (image->colorspace == CMYKColorspace)
 661           channel_features[BlackChannel].inverse_difference_moment[i]+=
 662             cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1);
 663         if (image->matte != MagickFalse)
 664           channel_features[OpacityChannel].inverse_difference_moment[i]+=
 665             cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
 666         /*
 667           Sum average.
 668         */
 669         density_xy[y+x+2].direction[i].red+=
 670           cooccurrence[x][y].direction[i].red;
 671         density_xy[y+x+2].direction[i].green+=
 672           cooccurrence[x][y].direction[i].green;
 673         density_xy[y+x+2].direction[i].blue+=
 674           cooccurrence[x][y].direction[i].blue;
 675         if (image->colorspace == CMYKColorspace)
 676           density_xy[y+x+2].direction[i].black+=
 677             cooccurrence[x][y].direction[i].black;
 678         if (image->matte != MagickFalse)
 679           density_xy[y+x+2].direction[i].alpha+=
 680             cooccurrence[x][y].direction[i].alpha;
 681         /*
 682           Entropy.
 683         */
 684         channel_features[RedChannel].entropy[i]-=
 685           cooccurrence[x][y].direction[i].red*
 686           log10(cooccurrence[x][y].direction[i].red+MagickEpsilon);
 687         channel_features[GreenChannel].entropy[i]-=
 688           cooccurrence[x][y].direction[i].green*
 689           log10(cooccurrence[x][y].direction[i].green+MagickEpsilon);
 690         channel_features[BlueChannel].entropy[i]-=
 691           cooccurrence[x][y].direction[i].blue*
 692           log10(cooccurrence[x][y].direction[i].blue+MagickEpsilon);
 693         if (image->colorspace == CMYKColorspace)
 694           channel_features[BlackChannel].entropy[i]-=
 695             cooccurrence[x][y].direction[i].black*
 696             log10(cooccurrence[x][y].direction[i].black+MagickEpsilon);
 697         if (image->matte != MagickFalse)
 698           channel_features[OpacityChannel].entropy[i]-=
 699             cooccurrence[x][y].direction[i].alpha*
 700             log10(cooccurrence[x][y].direction[i].alpha+MagickEpsilon);
 701         /*
 702           Information Measures of Correlation.
 703         */
 704         density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
 705         density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
 706         density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
 707         if (image->matte != MagickFalse)
 708           density_x[x].direction[i].alpha+=
 709             cooccurrence[x][y].direction[i].alpha;
 710         if (image->colorspace == CMYKColorspace)
 711           density_x[x].direction[i].black+=
 712             cooccurrence[x][y].direction[i].black;
 713         density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
 714         density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
 715         density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
 716         if (image->colorspace == CMYKColorspace)
 717           density_y[y].direction[i].black+=
 718             cooccurrence[x][y].direction[i].black;
 719         if (image->matte != MagickFalse)
 720           density_y[y].direction[i].alpha+=
 721             cooccurrence[x][y].direction[i].alpha;
 722       }
 723       mean.direction[i].red+=y*sum[y].direction[i].red;
 724       sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
 725       mean.direction[i].green+=y*sum[y].direction[i].green;
 726       sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
 727       mean.direction[i].blue+=y*sum[y].direction[i].blue;
 728       sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
 729       if (image->colorspace == CMYKColorspace)
 730         {
 731           mean.direction[i].black+=y*sum[y].direction[i].black;
 732           sum_squares.direction[i].black+=y*y*sum[y].direction[i].black;
 733         }
 734       if (image->matte != MagickFalse)
 735         {
 736           mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
 737           sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
 738         }
 739     }
 740     /*
 741       Correlation: measure of linear-dependencies in the image.
 742     */
 743     channel_features[RedChannel].correlation[i]=
 744       (correlation.direction[i].red-mean.direction[i].red*
 745       mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
 746       (mean.direction[i].red*mean.direction[i].red))*sqrt(
 747       sum_squares.direction[i].red-(mean.direction[i].red*
 748       mean.direction[i].red)));
 749     channel_features[GreenChannel].correlation[i]=
 750       (correlation.direction[i].green-mean.direction[i].green*
 751       mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
 752       (mean.direction[i].green*mean.direction[i].green))*sqrt(
 753       sum_squares.direction[i].green-(mean.direction[i].green*
 754       mean.direction[i].green)));
 755     channel_features[BlueChannel].correlation[i]=
 756       (correlation.direction[i].blue-mean.direction[i].blue*
 757       mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
 758       (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
 759       sum_squares.direction[i].blue-(mean.direction[i].blue*
 760       mean.direction[i].blue)));
 761     if (image->colorspace == CMYKColorspace)
 762       channel_features[BlackChannel].correlation[i]=
 763         (correlation.direction[i].black-mean.direction[i].black*
 764         mean.direction[i].black)/(sqrt(sum_squares.direction[i].black-
 765         (mean.direction[i].black*mean.direction[i].black))*sqrt(
 766         sum_squares.direction[i].black-(mean.direction[i].black*
 767         mean.direction[i].black)));
 768     if (image->matte != MagickFalse)
 769       channel_features[OpacityChannel].correlation[i]=
 770         (correlation.direction[i].alpha-mean.direction[i].alpha*
 771         mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
 772         (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt(
 773         sum_squares.direction[i].alpha-(mean.direction[i].alpha*
 774         mean.direction[i].alpha)));
 775   }
 776   /*
 777     Compute more texture features.
 778   */
 779 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 780   #pragma omp parallel for schedule(dynamic,4) shared(status)
 781 #endif
 782   for (i=0; i < 4; i++)
 783   {
 784     register ssize_t
 785       x;
 786
 787     for (x=2; x < (ssize_t) (2*number_grays); x++)
 788     {
 789       /*
 790         Sum average.
 791       */
 792       channel_features[RedChannel].sum_average[i]+=
 793         x*density_xy[x].direction[i].red;
 794       channel_features[GreenChannel].sum_average[i]+=
 795         x*density_xy[x].direction[i].green;
 796       channel_features[BlueChannel].sum_average[i]+=
 797         x*density_xy[x].direction[i].blue;
 798       if (image->colorspace == CMYKColorspace)
 799         channel_features[BlackChannel].sum_average[i]+=
 800           x*density_xy[x].direction[i].black;
 801       if (image->matte != MagickFalse)
 802         channel_features[OpacityChannel].sum_average[i]+=
 803           x*density_xy[x].direction[i].alpha;
 804       /*
 805         Sum entropy.
 806       */
 807       channel_features[RedChannel].sum_entropy[i]-=
 808         density_xy[x].direction[i].red*
 809         log10(density_xy[x].direction[i].red+MagickEpsilon);
 810       channel_features[GreenChannel].sum_entropy[i]-=
 811         density_xy[x].direction[i].green*
 812         log10(density_xy[x].direction[i].green+MagickEpsilon);
 813       channel_features[BlueChannel].sum_entropy[i]-=
 814         density_xy[x].direction[i].blue*
 815         log10(density_xy[x].direction[i].blue+MagickEpsilon);
 816       if (image->colorspace == CMYKColorspace)
 817         channel_features[BlackChannel].sum_entropy[i]-=
 818           density_xy[x].direction[i].black*
 819           log10(density_xy[x].direction[i].black+MagickEpsilon);
 820       if (image->matte != MagickFalse)
 821         channel_features[OpacityChannel].sum_entropy[i]-=
 822           density_xy[x].direction[i].alpha*
 823           log10(density_xy[x].direction[i].alpha+MagickEpsilon);
 824       /*
 825         Sum variance.
 826       */
 827       channel_features[RedChannel].sum_variance[i]+=
 828         (x-channel_features[RedChannel].sum_entropy[i])*
 829         (x-channel_features[RedChannel].sum_entropy[i])*
 830         density_xy[x].direction[i].red;
 831       channel_features[GreenChannel].sum_variance[i]+=
 832         (x-channel_features[GreenChannel].sum_entropy[i])*
 833         (x-channel_features[GreenChannel].sum_entropy[i])*
 834         density_xy[x].direction[i].green;
 835       channel_features[BlueChannel].sum_variance[i]+=
 836         (x-channel_features[BlueChannel].sum_entropy[i])*
 837         (x-channel_features[BlueChannel].sum_entropy[i])*
 838         density_xy[x].direction[i].blue;
 839       if (image->colorspace == CMYKColorspace)
 840         channel_features[BlackChannel].sum_variance[i]+=
 841           (x-channel_features[BlackChannel].sum_entropy[i])*
 842           (x-channel_features[BlackChannel].sum_entropy[i])*
 843           density_xy[x].direction[i].black;
 844       if (image->matte != MagickFalse)
 845         channel_features[OpacityChannel].sum_variance[i]+=
 846           (x-channel_features[OpacityChannel].sum_entropy[i])*
 847           (x-channel_features[OpacityChannel].sum_entropy[i])*
 848           density_xy[x].direction[i].alpha;
 849     }
 850   }
 851   /*
 852     Compute more texture features.
 853   */
 854 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 855   #pragma omp parallel for schedule(dynamic,4) shared(status)
 856 #endif
 857   for (i=0; i < 4; i++)
 858   {
 859     register ssize_t
 860       y;
 861
 862     for (y=0; y < (ssize_t) number_grays; y++)
 863     {
 864       register ssize_t
 865         x;
 866
 867       for (x=0; x < (ssize_t) number_grays; x++)
 868       {
 869         /*
 870           Sum of Squares: Variance
 871         */
 872         variance.direction[i].red+=(y-mean.direction[i].red+1)*
 873           (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
 874         variance.direction[i].green+=(y-mean.direction[i].green+1)*
 875           (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
 876         variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
 877           (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
 878         if (image->colorspace == CMYKColorspace)
 879           variance.direction[i].black+=(y-mean.direction[i].black+1)*
 880             (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black;
 881         if (image->matte != MagickFalse)
 882           variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
 883             (y-mean.direction[i].alpha+1)*
 884             cooccurrence[x][y].direction[i].alpha;
 885         /*
 886           Sum average / Difference Variance.
 887         */
 888         density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
 889           cooccurrence[x][y].direction[i].red;
 890         density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
 891           cooccurrence[x][y].direction[i].green;
 892         density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
 893           cooccurrence[x][y].direction[i].blue;
 894         if (image->colorspace == CMYKColorspace)
 895           density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
 896             cooccurrence[x][y].direction[i].black;
 897         if (image->matte != MagickFalse)
 898           density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
 899             cooccurrence[x][y].direction[i].alpha;
 900         /*
 901           Information Measures of Correlation.
 902         */
 903         entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
 904           log10(cooccurrence[x][y].direction[i].red+MagickEpsilon);
 905         entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
 906           log10(cooccurrence[x][y].direction[i].green+MagickEpsilon);
 907         entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
 908           log10(cooccurrence[x][y].direction[i].blue+MagickEpsilon);
 909         if (image->colorspace == CMYKColorspace)
 910           entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black*
 911             log10(cooccurrence[x][y].direction[i].black+MagickEpsilon);
 912         if (image->matte != MagickFalse)
 913           entropy_xy.direction[i].alpha-=
 914             cooccurrence[x][y].direction[i].alpha*log10(
 915             cooccurrence[x][y].direction[i].alpha+MagickEpsilon);
 916         entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
 917           log10(density_x[x].direction[i].red*density_y[y].direction[i].red+
 918           MagickEpsilon));
 919         entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
 920           log10(density_x[x].direction[i].green*density_y[y].direction[i].green+
 921           MagickEpsilon));
 922         entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
 923           log10(density_x[x].direction[i].blue*density_y[y].direction[i].blue+
 924           MagickEpsilon));
 925         if (image->colorspace == CMYKColorspace)
 926           entropy_xy1.direction[i].black-=(
 927             cooccurrence[x][y].direction[i].black*log10(
 928             density_x[x].direction[i].black*density_y[y].direction[i].black+
 929             MagickEpsilon));
 930         if (image->matte != MagickFalse)
 931           entropy_xy1.direction[i].alpha-=(
 932             cooccurrence[x][y].direction[i].alpha*log10(
 933             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha+
 934             MagickEpsilon));
 935         entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
 936           density_y[y].direction[i].red*log10(density_x[x].direction[i].red*
 937           density_y[y].direction[i].red+MagickEpsilon));
 938         entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
 939           density_y[y].direction[i].green*log10(density_x[x].direction[i].green*
 940           density_y[y].direction[i].green+MagickEpsilon));
 941         entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
 942           density_y[y].direction[i].blue*log10(density_x[x].direction[i].blue*
 943           density_y[y].direction[i].blue+MagickEpsilon));
 944         if (image->colorspace == CMYKColorspace)
 945           entropy_xy2.direction[i].black-=(density_x[x].direction[i].black*
 946             density_y[y].direction[i].black*log10(
 947             density_x[x].direction[i].black*density_y[y].direction[i].black+
 948             MagickEpsilon));
 949         if (image->matte != MagickFalse)
 950           entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha*
 951             density_y[y].direction[i].alpha*log10(
 952             density_x[x].direction[i].alpha*density_y[y].direction[i].alpha+
 953             MagickEpsilon));
 954       }
 955     }
 956     channel_features[RedChannel].variance_sum_of_squares[i]=
 957       variance.direction[i].red;
 958     channel_features[GreenChannel].variance_sum_of_squares[i]=
 959       variance.direction[i].green;
 960     channel_features[BlueChannel].variance_sum_of_squares[i]=
 961       variance.direction[i].blue;
 962     if (image->colorspace == CMYKColorspace)
 963       channel_features[RedChannel].variance_sum_of_squares[i]=
 964         variance.direction[i].black;
 965     if (image->matte != MagickFalse)
 966       channel_features[RedChannel].variance_sum_of_squares[i]=
 967         variance.direction[i].alpha;
 968   }
 969   /*
 970     Compute more texture features.
 971   */
 972   (void) ResetMagickMemory(&variance,0,sizeof(variance));
 973   (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
 974 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 975   #pragma omp parallel for schedule(dynamic,4) shared(status)
 976 #endif
 977   for (i=0; i < 4; i++)
 978   {
 979     register ssize_t
 980       x;
 981
 982     for (x=0; x < (ssize_t) number_grays; x++)
 983     {
 984       /*
 985         Difference variance.
 986       */
 987       variance.direction[i].red+=density_xy[x].direction[i].red;
 988       variance.direction[i].green+=density_xy[x].direction[i].green;
 989       variance.direction[i].blue+=density_xy[x].direction[i].blue;
 990       if (image->colorspace == CMYKColorspace)
 991         variance.direction[i].black+=density_xy[x].direction[i].black;
 992       if (image->matte != MagickFalse)
 993         variance.direction[i].alpha+=density_xy[x].direction[i].alpha;
 994       sum_squares.direction[i].red+=density_xy[x].direction[i].red*
 995         density_xy[x].direction[i].red;
 996       sum_squares.direction[i].green+=density_xy[x].direction[i].green*
 997         density_xy[x].direction[i].green;
 998       sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
 999         density_xy[x].direction[i].blue;
1000       if (image->colorspace == CMYKColorspace)
1001         sum_squares.direction[i].black+=density_xy[x].direction[i].black*
1002           density_xy[x].direction[i].black;
1003       if (image->matte != MagickFalse)
1004         sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
1005           density_xy[x].direction[i].alpha;
1006       /*
1007         Difference entropy.
1008       */
1009       channel_features[RedChannel].difference_entropy[i]-=
1010         density_xy[x].direction[i].red*
1011         log10(density_xy[x].direction[i].red+MagickEpsilon);
1012       channel_features[GreenChannel].difference_entropy[i]-=
1013         density_xy[x].direction[i].green*
1014         log10(density_xy[x].direction[i].green+MagickEpsilon);
1015       channel_features[BlueChannel].difference_entropy[i]-=
1016         density_xy[x].direction[i].blue*
1017         log10(density_xy[x].direction[i].blue+MagickEpsilon);
1018       if (image->colorspace == CMYKColorspace)
1019         channel_features[BlackChannel].difference_entropy[i]-=
1020           density_xy[x].direction[i].black*
1021           log10(density_xy[x].direction[i].black+MagickEpsilon);
1022       if (image->matte != MagickFalse)
1023         channel_features[OpacityChannel].difference_entropy[i]-=
1024           density_xy[x].direction[i].alpha*
1025           log10(density_xy[x].direction[i].alpha+MagickEpsilon);
1026       /*
1027         Information Measures of Correlation.
1028       */
1029       entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1030         log10(density_x[x].direction[i].red+MagickEpsilon));
1031       entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1032         log10(density_x[x].direction[i].green+MagickEpsilon));
1033       entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1034         log10(density_x[x].direction[i].blue+MagickEpsilon));
1035       if (image->colorspace == CMYKColorspace)
1036         entropy_x.direction[i].black-=(density_x[x].direction[i].black*
1037           log10(density_x[x].direction[i].black+MagickEpsilon));
1038       if (image->matte != MagickFalse)
1039         entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha*
1040           log10(density_x[x].direction[i].alpha+MagickEpsilon));
1041       entropy_y.direction[i].red-=(density_y[x].direction[i].red*
1042         log10(density_y[x].direction[i].red+MagickEpsilon));
1043       entropy_y.direction[i].green-=(density_y[x].direction[i].green*
1044         log10(density_y[x].direction[i].green+MagickEpsilon));
1045       entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
1046         log10(density_y[x].direction[i].blue+MagickEpsilon));
1047       if (image->colorspace == CMYKColorspace)
1048         entropy_y.direction[i].black-=(density_y[x].direction[i].black*
1049           log10(density_y[x].direction[i].black+MagickEpsilon));
1050       if (image->matte != MagickFalse)
1051         entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
1052           log10(density_y[x].direction[i].alpha+MagickEpsilon));
1053     }
1054     /*
1055       Difference variance.
1056     */
1057     channel_features[RedChannel].difference_variance[i]=
1058       (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1059       (variance.direction[i].red*variance.direction[i].red))/
1060       ((double) number_grays*number_grays*number_grays*number_grays);
1061     channel_features[GreenChannel].difference_variance[i]=
1062       (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1063       (variance.direction[i].green*variance.direction[i].green))/
1064       ((double) number_grays*number_grays*number_grays*number_grays);
1065     channel_features[BlueChannel].difference_variance[i]=
1066       (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1067       (variance.direction[i].blue*variance.direction[i].blue))/
1068       ((double) number_grays*number_grays*number_grays*number_grays);
1069     if (image->colorspace == CMYKColorspace)
1070       channel_features[BlackChannel].difference_variance[i]=
1071         (((double) number_grays*number_grays*sum_squares.direction[i].black)-
1072         (variance.direction[i].black*variance.direction[i].black))/
1073         ((double) number_grays*number_grays*number_grays*number_grays);
1074     if (image->matte != MagickFalse)
1075       channel_features[OpacityChannel].difference_variance[i]=
1076         (((double) number_grays*number_grays*sum_squares.direction[i].alpha)-
1077         (variance.direction[i].alpha*variance.direction[i].alpha))/
1078         ((double) number_grays*number_grays*number_grays*number_grays);
1079     /*
1080       Information Measures of Correlation.
1081     */
1082     channel_features[RedChannel].measure_of_correlation_1[i]=
1083       (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1084       (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1085        entropy_x.direction[i].red : entropy_y.direction[i].red);
1086     channel_features[GreenChannel].measure_of_correlation_1[i]=
1087       (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1088       (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1089        entropy_x.direction[i].green : entropy_y.direction[i].green);
1090     channel_features[BlueChannel].measure_of_correlation_1[i]=
1091       (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1092       (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1093        entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1094     if (image->colorspace == CMYKColorspace)
1095       channel_features[BlackChannel].measure_of_correlation_1[i]=
1096         (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/
1097         (entropy_x.direction[i].black > entropy_y.direction[i].black ?
1098          entropy_x.direction[i].black : entropy_y.direction[i].black);
1099     if (image->matte != MagickFalse)
1100       channel_features[OpacityChannel].measure_of_correlation_1[i]=
1101         (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/
1102         (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ?
1103          entropy_x.direction[i].alpha : entropy_y.direction[i].alpha);
1104     channel_features[RedChannel].measure_of_correlation_2[i]=
1105       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].red-
1106       entropy_xy.direction[i].red)))));
1107     channel_features[GreenChannel].measure_of_correlation_2[i]=
1108       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].green-
1109       entropy_xy.direction[i].green)))));
1110     channel_features[BlueChannel].measure_of_correlation_2[i]=
1111       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].blue-
1112       entropy_xy.direction[i].blue)))));
1113     if (image->colorspace == CMYKColorspace)
1114       channel_features[BlackChannel].measure_of_correlation_2[i]=
1115         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].black-
1116         entropy_xy.direction[i].black)))));
1117     if (image->matte != MagickFalse)
1118       channel_features[OpacityChannel].measure_of_correlation_2[i]=
1119         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].alpha-
1120         entropy_xy.direction[i].alpha)))));
1121   }
1122   /*
1123     Compute more texture features.
1124   */
1125 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1126   #pragma omp parallel for schedule(dynamic,4) shared(status)
1127 #endif
1128   for (i=0; i < 4; i++)
1129   {
1130     for (z=0; z < (ssize_t) number_grays; z++)
1131     {
1132       register ssize_t
1133         y;
1134
1135       ChannelStatistics
1136         pixel;
1137
1138       (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
1139       for (y=0; y < (ssize_t) number_grays; y++)
1140       {
1141         register ssize_t
1142           x;
1143
1144         for (x=0; x < (ssize_t) number_grays; x++)
1145         {
1146           /*
1147             Contrast:  amount of local variations present in an image.
1148           */
1149           if (((y-x) == z) || ((x-y) == z))
1150             {
1151               pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1152               pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1153               pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1154               if (image->colorspace == CMYKColorspace)
1155                 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
1156               if (image->matte != MagickFalse)
1157                 pixel.direction[i].alpha+=
1158                   cooccurrence[x][y].direction[i].alpha;
1159             }
1160           /*
1161             Maximum Correlation Coefficient.
1162           */
1163           Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1164             cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1165             density_y[x].direction[i].red;
1166           Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1167             cooccurrence[y][x].direction[i].green/
1168             density_x[z].direction[i].green/density_y[x].direction[i].red;
1169           Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1170             cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/
1171             density_y[x].direction[i].blue;
1172           if (image->colorspace == CMYKColorspace)
1173             Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
1174               cooccurrence[y][x].direction[i].black/
1175               density_x[z].direction[i].black/density_y[x].direction[i].black;
1176           if (image->matte != MagickFalse)
1177             Q[z][y].direction[i].alpha+=
1178               cooccurrence[z][x].direction[i].alpha*
1179               cooccurrence[y][x].direction[i].alpha/
1180               density_x[z].direction[i].alpha/
1181               density_y[x].direction[i].alpha;
1182         }
1183       }
1184       channel_features[RedChannel].contrast[i]+=z*z*pixel.direction[i].red;
1185       channel_features[GreenChannel].contrast[i]+=z*z*pixel.direction[i].green;
1186       channel_features[BlueChannel].contrast[i]+=z*z*pixel.direction[i].blue;
1187       if (image->colorspace == CMYKColorspace)
1188         channel_features[BlackChannel].contrast[i]+=z*z*
1189           pixel.direction[i].black;
1190       if (image->matte != MagickFalse)
1191         channel_features[OpacityChannel].contrast[i]+=z*z*
1192           pixel.direction[i].alpha;
1193     }
1194     /*
1195       Maximum Correlation Coefficient.
1196       Future: return second largest eigenvalue of Q.
1197     */
1198     channel_features[RedChannel].maximum_correlation_coefficient[i]=
1199       sqrt((double) -1.0);
1200     channel_features[GreenChannel].maximum_correlation_coefficient[i]=
1201       sqrt((double) -1.0);
1202     channel_features[BlueChannel].maximum_correlation_coefficient[i]=
1203       sqrt((double) -1.0);
1204     if (image->colorspace == CMYKColorspace)
1205       channel_features[BlackChannel].maximum_correlation_coefficient[i]=
1206         sqrt((double) -1.0);
1207     if (image->matte != MagickFalse)
1208       channel_features[OpacityChannel].maximum_correlation_coefficient[i]=
1209         sqrt((double) -1.0);
1210   }
1211   /*
1212     Relinquish resources.
1213   */
1214   sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1215   for (i=0; i < (ssize_t) number_grays; i++)
1216     Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1217   Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1218   density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1219   density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1220   density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1221   for (i=0; i < (ssize_t) number_grays; i++)
1222     cooccurrence[i]=(ChannelStatistics *)
1223       RelinquishMagickMemory(cooccurrence[i]);
1224   cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1225   return(channel_features);
1226 }