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