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

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