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