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