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

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