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,exception);
 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   MagickBooleanType
 175     status;
 176
 177   register ssize_t
 178     i;
 179
 180   size_t
 181     length;
 182
 183   ssize_t
 184     y,
 185     z;
 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     register const IndexPacket
 409       *restrict indexes;
 410
 411     register const PixelPacket
 412       *restrict p;
 413
 414     register ssize_t
 415       x;
 416
 417     ssize_t
 418       offset,
 419       u,
 420       v;
 421
 422     if (status == MagickFalse)
 423       continue;
 424     p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,y,image->columns+
 425       2*distance,distance+1,exception);
 426     if (p == (const PixelPacket *) NULL)
 427       {
 428         status=MagickFalse;
 429         continue;
 430       }
 431     indexes=GetCacheViewVirtualIndexQueue(image_view);
 432     p+=distance;
 433     indexes+=distance;
 434     for (x=0; x < (ssize_t) image->columns; x++)
 435     {
 436       for (i=0; i < 4; i++)
 437       {
 438         switch (i)
 439         {
 440           case 0:
 441           default:
 442           {
 443             /*
 444               Horizontal adjacency.
 445             */
 446             offset=(ssize_t) distance;
 447             break;
 448           }
 449           case 1:
 450           {
 451             /*
 452               Vertical adjacency.
 453             */
 454             offset=(ssize_t) (image->columns+2*distance);
 455             break;
 456           }
 457           case 2:
 458           {
 459             /*
 460               Right diagonal adjacency.
 461             */
 462             offset=(ssize_t) ((image->columns+2*distance)-distance);
 463             break;
 464           }
 465           case 3:
 466           {
 467             /*
 468               Left diagonal adjacency.
 469             */
 470             offset=(ssize_t) ((image->columns+2*distance)+distance);
 471             break;
 472           }
 473         }
 474         u=0;
 475         v=0;
 476         while (grays[u].red != ScaleQuantumToMap(p->red))
 477           u++;
 478         while (grays[v].red != ScaleQuantumToMap((p+offset)->red))
 479           v++;
 480         cooccurrence[u][v].direction[i].red++;
 481         cooccurrence[v][u].direction[i].red++;
 482         u=0;
 483         v=0;
 484         while (grays[u].green != ScaleQuantumToMap(p->green))
 485           u++;
 486         while (grays[v].green != ScaleQuantumToMap((p+offset)->green))
 487           v++;
 488         cooccurrence[u][v].direction[i].green++;
 489         cooccurrence[v][u].direction[i].green++;
 490         u=0;
 491         v=0;
 492         while (grays[u].blue != ScaleQuantumToMap(p->blue))
 493           u++;
 494         while (grays[v].blue != ScaleQuantumToMap((p+offset)->blue))
 495           v++;
 496         cooccurrence[u][v].direction[i].blue++;
 497         cooccurrence[v][u].direction[i].blue++;
 498         if (image->matte != MagickFalse)
 499           {
 500             u=0;
 501             v=0;
 502             while (grays[u].opacity != ScaleQuantumToMap(p->opacity))
 503               u++;
 504             while (grays[v].opacity != ScaleQuantumToMap((p+offset)->opacity))
 505               v++;
 506             cooccurrence[u][v].direction[i].opacity++;
 507             cooccurrence[v][u].direction[i].opacity++;
 508           }
 509         if (image->colorspace == CMYKColorspace)
 510           {
 511             u=0;
 512             v=0;
 513             while (grays[u].index != ScaleQuantumToMap(indexes[x]))
 514               u++;
 515             while (grays[v].index != ScaleQuantumToMap(indexes[x+offset]))
 516               v++;
 517             cooccurrence[u][v].direction[i].index++;
 518             cooccurrence[v][u].direction[i].index++;
 519           }
 520       }
 521       p++;
 522     }
 523   }
 524   grays=(LongPixelPacket *) RelinquishMagickMemory(grays);
 525   image_view=DestroyCacheView(image_view);
 526   if (status == MagickFalse)
 527     {
 528       for (i=0; i < (ssize_t) number_grays; i++)
 529         cooccurrence[i]=(ChannelStatistics *)
 530           RelinquishMagickMemory(cooccurrence[i]);
 531       cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
 532       channel_features=(ChannelFeatures *) RelinquishMagickMemory(
 533         channel_features);
 534       (void) ThrowMagickException(exception,GetMagickModule(),
 535         ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
 536       return(channel_features);
 537     }
 538   /*
 539     Normalize spatial dependence matrix.
 540   */
 541 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 542   #pragma omp parallel for schedule(dynamic,4) shared(status)
 543 #endif
 544   for (i=0; i < 4; i++)
 545   {
 546     double
 547       normalize;
 548
 549     switch (i)
 550     {
 551       case 0:
 552       default:
 553       {
 554         /*
 555           Horizontal adjacency.
 556         */
 557         normalize=2.0*image->rows*(image->columns-distance);
 558         break;
 559       }
 560       case 1:
 561       {
 562         /*
 563           Vertical adjacency.
 564         */
 565         normalize=2.0*(image->rows-distance)*image->columns;
 566         break;
 567       }
 568       case 2:
 569       {
 570         /*
 571           Right diagonal adjacency.
 572         */
 573         normalize=2.0*(image->rows-distance)*(image->columns-distance);
 574         break;
 575       }
 576       case 3:
 577       {
 578         /*
 579           Left diagonal adjacency.
 580         */
 581         normalize=2.0*(image->rows-distance)*(image->columns-distance);
 582         break;
 583       }
 584     }
 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->matte != MagickFalse)
 596           cooccurrence[x][y].direction[i].opacity/=normalize;
 597         if (image->colorspace == CMYKColorspace)
 598           cooccurrence[x][y].direction[i].index/=normalize;
 599       }
 600     }
 601   }
 602   /*
 603     Compute texture features.
 604   */
 605 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 606   #pragma omp parallel for schedule(dynamic,4) shared(status)
 607 #endif
 608   for (i=0; i < 4; i++)
 609   {
 610     register ssize_t
 611       y;
 612
 613     for (y=0; y < (ssize_t) number_grays; y++)
 614     {
 615       register ssize_t
 616         x;
 617
 618       for (x=0; x < (ssize_t) number_grays; x++)
 619       {
 620         /*
 621           Angular second moment:  measure of homogeneity of the image.
 622         */
 623         channel_features[RedChannel].angular_second_moment[i]+=
 624           cooccurrence[x][y].direction[i].red*
 625           cooccurrence[x][y].direction[i].red;
 626         channel_features[GreenChannel].angular_second_moment[i]+=
 627           cooccurrence[x][y].direction[i].green*
 628           cooccurrence[x][y].direction[i].green;
 629         channel_features[BlueChannel].angular_second_moment[i]+=
 630           cooccurrence[x][y].direction[i].blue*
 631           cooccurrence[x][y].direction[i].blue;
 632         if (image->matte != MagickFalse)
 633           channel_features[OpacityChannel].angular_second_moment[i]+=
 634             cooccurrence[x][y].direction[i].opacity*
 635             cooccurrence[x][y].direction[i].opacity;
 636         if (image->colorspace == CMYKColorspace)
 637           channel_features[BlackChannel].angular_second_moment[i]+=
 638             cooccurrence[x][y].direction[i].index*
 639             cooccurrence[x][y].direction[i].index;
 640         /*
 641           Correlation: measure of linear-dependencies in the image.
 642         */
 643         sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
 644         sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
 645         sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
 646         if (image->matte != MagickFalse)
 647           sum[y].direction[i].opacity+=cooccurrence[x][y].direction[i].opacity;
 648         if (image->colorspace == CMYKColorspace)
 649           sum[y].direction[i].index+=cooccurrence[x][y].direction[i].index;
 650         correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
 651         correlation.direction[i].green+=x*y*
 652           cooccurrence[x][y].direction[i].green;
 653         correlation.direction[i].blue+=x*y*
 654           cooccurrence[x][y].direction[i].blue;
 655         if (image->matte != MagickFalse)
 656           correlation.direction[i].opacity+=x*y*
 657             cooccurrence[x][y].direction[i].opacity;
 658         if (image->colorspace == CMYKColorspace)
 659           correlation.direction[i].index+=x*y*
 660             cooccurrence[x][y].direction[i].index;
 661         /*
 662           Inverse Difference Moment.
 663         */
 664         channel_features[RedChannel].inverse_difference_moment[i]+=
 665           cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
 666         channel_features[GreenChannel].inverse_difference_moment[i]+=
 667           cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
 668         channel_features[BlueChannel].inverse_difference_moment[i]+=
 669           cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
 670         if (image->matte != MagickFalse)
 671           channel_features[OpacityChannel].inverse_difference_moment[i]+=
 672             cooccurrence[x][y].direction[i].opacity/((y-x)*(y-x)+1);
 673         if (image->colorspace == CMYKColorspace)
 674           channel_features[IndexChannel].inverse_difference_moment[i]+=
 675             cooccurrence[x][y].direction[i].index/((y-x)*(y-x)+1);
 676         /*
 677           Sum average.
 678         */
 679         density_xy[y+x+2].direction[i].red+=
 680           cooccurrence[x][y].direction[i].red;
 681         density_xy[y+x+2].direction[i].green+=
 682           cooccurrence[x][y].direction[i].green;
 683         density_xy[y+x+2].direction[i].blue+=
 684           cooccurrence[x][y].direction[i].blue;
 685         if (image->matte != MagickFalse)
 686           density_xy[y+x+2].direction[i].opacity+=
 687             cooccurrence[x][y].direction[i].opacity;
 688         if (image->colorspace == CMYKColorspace)
 689           density_xy[y+x+2].direction[i].index+=
 690             cooccurrence[x][y].direction[i].index;
 691         /*
 692           Entropy.
 693         */
 694         channel_features[RedChannel].entropy[i]-=
 695           cooccurrence[x][y].direction[i].red*
 696           log10(cooccurrence[x][y].direction[i].red+MagickEpsilon);
 697         channel_features[GreenChannel].entropy[i]-=
 698           cooccurrence[x][y].direction[i].green*
 699           log10(cooccurrence[x][y].direction[i].green+MagickEpsilon);
 700         channel_features[BlueChannel].entropy[i]-=
 701           cooccurrence[x][y].direction[i].blue*
 702           log10(cooccurrence[x][y].direction[i].blue+MagickEpsilon);
 703         if (image->matte != MagickFalse)
 704           channel_features[OpacityChannel].entropy[i]-=
 705             cooccurrence[x][y].direction[i].opacity*
 706             log10(cooccurrence[x][y].direction[i].opacity+MagickEpsilon);
 707         if (image->colorspace == CMYKColorspace)
 708           channel_features[IndexChannel].entropy[i]-=
 709             cooccurrence[x][y].direction[i].index*
 710             log10(cooccurrence[x][y].direction[i].index+MagickEpsilon);
 711         /*
 712           Information Measures of Correlation.
 713         */
 714         density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
 715         density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
 716         density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
 717         if (image->matte != MagickFalse)
 718           density_x[x].direction[i].opacity+=
 719             cooccurrence[x][y].direction[i].opacity;
 720         if (image->colorspace == CMYKColorspace)
 721           density_x[x].direction[i].index+=
 722             cooccurrence[x][y].direction[i].index;
 723         density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
 724         density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
 725         density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
 726         if (image->matte != MagickFalse)
 727           density_y[y].direction[i].opacity+=
 728             cooccurrence[x][y].direction[i].opacity;
 729         if (image->colorspace == CMYKColorspace)
 730           density_y[y].direction[i].index+=
 731             cooccurrence[x][y].direction[i].index;
 732       }
 733       mean.direction[i].red+=y*sum[y].direction[i].red;
 734       sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
 735       mean.direction[i].green+=y*sum[y].direction[i].green;
 736       sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
 737       mean.direction[i].blue+=y*sum[y].direction[i].blue;
 738       sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
 739       if (image->matte != MagickFalse)
 740         {
 741           mean.direction[i].opacity+=y*sum[y].direction[i].opacity;
 742           sum_squares.direction[i].opacity+=y*y*sum[y].direction[i].opacity;
 743         }
 744       if (image->colorspace == CMYKColorspace)
 745         {
 746           mean.direction[i].index+=y*sum[y].direction[i].index;
 747           sum_squares.direction[i].index+=y*y*sum[y].direction[i].index;
 748         }
 749     }
 750     /*
 751       Correlation: measure of linear-dependencies in the image.
 752     */
 753     channel_features[RedChannel].correlation[i]=
 754       (correlation.direction[i].red-mean.direction[i].red*
 755       mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
 756       (mean.direction[i].red*mean.direction[i].red))*sqrt(
 757       sum_squares.direction[i].red-(mean.direction[i].red*
 758       mean.direction[i].red)));
 759     channel_features[GreenChannel].correlation[i]=
 760       (correlation.direction[i].green-mean.direction[i].green*
 761       mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
 762       (mean.direction[i].green*mean.direction[i].green))*sqrt(
 763       sum_squares.direction[i].green-(mean.direction[i].green*
 764       mean.direction[i].green)));
 765     channel_features[BlueChannel].correlation[i]=
 766       (correlation.direction[i].blue-mean.direction[i].blue*
 767       mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
 768       (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
 769       sum_squares.direction[i].blue-(mean.direction[i].blue*
 770       mean.direction[i].blue)));
 771     if (image->matte != MagickFalse)
 772       channel_features[OpacityChannel].correlation[i]=
 773         (correlation.direction[i].opacity-mean.direction[i].opacity*
 774         mean.direction[i].opacity)/(sqrt(sum_squares.direction[i].opacity-
 775         (mean.direction[i].opacity*mean.direction[i].opacity))*sqrt(
 776         sum_squares.direction[i].opacity-(mean.direction[i].opacity*
 777         mean.direction[i].opacity)));
 778     if (image->colorspace == CMYKColorspace)
 779       channel_features[IndexChannel].correlation[i]=
 780         (correlation.direction[i].index-mean.direction[i].index*
 781         mean.direction[i].index)/(sqrt(sum_squares.direction[i].index-
 782         (mean.direction[i].index*mean.direction[i].index))*sqrt(
 783         sum_squares.direction[i].index-(mean.direction[i].index*
 784         mean.direction[i].index)));
 785   }
 786   /*
 787     Compute more texture features.
 788   */
 789 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 790   #pragma omp parallel for schedule(dynamic,4) shared(status)
 791 #endif
 792   for (i=0; i < 4; i++)
 793   {
 794     register ssize_t
 795       x;
 796
 797     for (x=2; x < (ssize_t) (2*number_grays); x++)
 798     {
 799       /*
 800         Sum average.
 801       */
 802       channel_features[RedChannel].sum_average[i]+=
 803         x*density_xy[x].direction[i].red;
 804       channel_features[GreenChannel].sum_average[i]+=
 805         x*density_xy[x].direction[i].green;
 806       channel_features[BlueChannel].sum_average[i]+=
 807         x*density_xy[x].direction[i].blue;
 808       if (image->matte != MagickFalse)
 809         channel_features[OpacityChannel].sum_average[i]+=
 810           x*density_xy[x].direction[i].opacity;
 811       if (image->colorspace == CMYKColorspace)
 812         channel_features[IndexChannel].sum_average[i]+=
 813           x*density_xy[x].direction[i].index;
 814       /*
 815         Sum entropy.
 816       */
 817       channel_features[RedChannel].sum_entropy[i]-=
 818         density_xy[x].direction[i].red*
 819         log10(density_xy[x].direction[i].red+MagickEpsilon);
 820       channel_features[GreenChannel].sum_entropy[i]-=
 821         density_xy[x].direction[i].green*
 822         log10(density_xy[x].direction[i].green+MagickEpsilon);
 823       channel_features[BlueChannel].sum_entropy[i]-=
 824         density_xy[x].direction[i].blue*
 825         log10(density_xy[x].direction[i].blue+MagickEpsilon);
 826       if (image->matte != MagickFalse)
 827         channel_features[OpacityChannel].sum_entropy[i]-=
 828           density_xy[x].direction[i].opacity*
 829           log10(density_xy[x].direction[i].opacity+MagickEpsilon);
 830       if (image->colorspace == CMYKColorspace)
 831         channel_features[IndexChannel].sum_entropy[i]-=
 832           density_xy[x].direction[i].index*
 833           log10(density_xy[x].direction[i].index+MagickEpsilon);
 834       /*
 835         Sum variance.
 836       */
 837       channel_features[RedChannel].sum_variance[i]+=
 838         (x-channel_features[RedChannel].sum_entropy[i])*
 839         (x-channel_features[RedChannel].sum_entropy[i])*
 840         density_xy[x].direction[i].red;
 841       channel_features[GreenChannel].sum_variance[i]+=
 842         (x-channel_features[GreenChannel].sum_entropy[i])*
 843         (x-channel_features[GreenChannel].sum_entropy[i])*
 844         density_xy[x].direction[i].green;
 845       channel_features[BlueChannel].sum_variance[i]+=
 846         (x-channel_features[BlueChannel].sum_entropy[i])*
 847         (x-channel_features[BlueChannel].sum_entropy[i])*
 848         density_xy[x].direction[i].blue;
 849       if (image->matte != MagickFalse)
 850         channel_features[OpacityChannel].sum_variance[i]+=
 851           (x-channel_features[OpacityChannel].sum_entropy[i])*
 852           (x-channel_features[OpacityChannel].sum_entropy[i])*
 853           density_xy[x].direction[i].opacity;
 854       if (image->colorspace == CMYKColorspace)
 855         channel_features[IndexChannel].sum_variance[i]+=
 856           (x-channel_features[IndexChannel].sum_entropy[i])*
 857           (x-channel_features[IndexChannel].sum_entropy[i])*
 858           density_xy[x].direction[i].index;
 859     }
 860   }
 861   /*
 862     Compute more texture features.
 863   */
 864 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 865   #pragma omp parallel for schedule(dynamic,4) shared(status)
 866 #endif
 867   for (i=0; i < 4; i++)
 868   {
 869     register ssize_t
 870       y;
 871
 872     for (y=0; y < (ssize_t) number_grays; y++)
 873     {
 874       register ssize_t
 875         x;
 876
 877       for (x=0; x < (ssize_t) number_grays; x++)
 878       {
 879         /*
 880           Sum of Squares: Variance
 881         */
 882         variance.direction[i].red+=(y-mean.direction[i].red+1)*
 883           (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
 884         variance.direction[i].green+=(y-mean.direction[i].green+1)*
 885           (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
 886         variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
 887           (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
 888         if (image->matte != MagickFalse)
 889           variance.direction[i].opacity+=(y-mean.direction[i].opacity+1)*
 890             (y-mean.direction[i].opacity+1)*
 891             cooccurrence[x][y].direction[i].opacity;
 892         if (image->colorspace == CMYKColorspace)
 893           variance.direction[i].index+=(y-mean.direction[i].index+1)*
 894             (y-mean.direction[i].index+1)*cooccurrence[x][y].direction[i].index;
 895         /*
 896           Sum average / Difference Variance.
 897         */
 898         density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
 899           cooccurrence[x][y].direction[i].red;
 900         density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
 901           cooccurrence[x][y].direction[i].green;
 902         density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
 903           cooccurrence[x][y].direction[i].blue;
 904         if (image->matte != MagickFalse)
 905           density_xy[MagickAbsoluteValue(y-x)].direction[i].opacity+=
 906             cooccurrence[x][y].direction[i].opacity;
 907         if (image->colorspace == CMYKColorspace)
 908           density_xy[MagickAbsoluteValue(y-x)].direction[i].index+=
 909             cooccurrence[x][y].direction[i].index;
 910         /*
 911           Information Measures of Correlation.
 912         */
 913         entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
 914           log10(cooccurrence[x][y].direction[i].red+MagickEpsilon);
 915         entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
 916           log10(cooccurrence[x][y].direction[i].green+MagickEpsilon);
 917         entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
 918           log10(cooccurrence[x][y].direction[i].blue+MagickEpsilon);
 919         if (image->matte != MagickFalse)
 920           entropy_xy.direction[i].opacity-=
 921             cooccurrence[x][y].direction[i].opacity*log10(
 922             cooccurrence[x][y].direction[i].opacity+MagickEpsilon);
 923         if (image->colorspace == CMYKColorspace)
 924           entropy_xy.direction[i].index-=cooccurrence[x][y].direction[i].index*
 925             log10(cooccurrence[x][y].direction[i].index+MagickEpsilon);
 926         entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
 927           log10(density_x[x].direction[i].red*density_y[y].direction[i].red+
 928           MagickEpsilon));
 929         entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
 930           log10(density_x[x].direction[i].green*density_y[y].direction[i].green+
 931           MagickEpsilon));
 932         entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
 933           log10(density_x[x].direction[i].blue*density_y[y].direction[i].blue+
 934           MagickEpsilon));
 935         if (image->matte != MagickFalse)
 936           entropy_xy1.direction[i].opacity-=(
 937             cooccurrence[x][y].direction[i].opacity*log10(
 938             density_x[x].direction[i].opacity*density_y[y].direction[i].opacity+
 939             MagickEpsilon));
 940         if (image->colorspace == CMYKColorspace)
 941           entropy_xy1.direction[i].index-=(
 942             cooccurrence[x][y].direction[i].index*log10(
 943             density_x[x].direction[i].index*density_y[y].direction[i].index+
 944             MagickEpsilon));
 945         entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
 946           density_y[y].direction[i].red*log10(density_x[x].direction[i].red*
 947           density_y[y].direction[i].red+MagickEpsilon));
 948         entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
 949           density_y[y].direction[i].green*log10(density_x[x].direction[i].green*
 950           density_y[y].direction[i].green+MagickEpsilon));
 951         entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
 952           density_y[y].direction[i].blue*log10(density_x[x].direction[i].blue*
 953           density_y[y].direction[i].blue+MagickEpsilon));
 954         if (image->matte != MagickFalse)
 955           entropy_xy2.direction[i].opacity-=(density_x[x].direction[i].opacity*
 956             density_y[y].direction[i].opacity*log10(
 957             density_x[x].direction[i].opacity*density_y[y].direction[i].opacity+
 958             MagickEpsilon));
 959         if (image->colorspace == CMYKColorspace)
 960           entropy_xy2.direction[i].index-=(density_x[x].direction[i].index*
 961             density_y[y].direction[i].index*log10(
 962             density_x[x].direction[i].index*density_y[y].direction[i].index+
 963             MagickEpsilon));
 964       }
 965     }
 966     channel_features[RedChannel].variance_sum_of_squares[i]=
 967       variance.direction[i].red;
 968     channel_features[GreenChannel].variance_sum_of_squares[i]=
 969       variance.direction[i].green;
 970     channel_features[BlueChannel].variance_sum_of_squares[i]=
 971       variance.direction[i].blue;
 972     if (image->matte != MagickFalse)
 973       channel_features[RedChannel].variance_sum_of_squares[i]=
 974         variance.direction[i].opacity;
 975     if (image->colorspace == CMYKColorspace)
 976       channel_features[RedChannel].variance_sum_of_squares[i]=
 977         variance.direction[i].index;
 978   }
 979   /*
 980     Compute more texture features.
 981   */
 982   (void) ResetMagickMemory(&variance,0,sizeof(variance));
 983   (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
 984 #if defined(MAGICKCORE_OPENMP_SUPPORT)
 985   #pragma omp parallel for schedule(dynamic,4) shared(status)
 986 #endif
 987   for (i=0; i < 4; i++)
 988   {
 989     register ssize_t
 990       x;
 991
 992     for (x=0; x < (ssize_t) number_grays; x++)
 993     {
 994       /*
 995         Difference variance.
 996       */
 997       variance.direction[i].red+=density_xy[x].direction[i].red;
 998       variance.direction[i].green+=density_xy[x].direction[i].green;
 999       variance.direction[i].blue+=density_xy[x].direction[i].blue;
1000       if (image->matte != MagickFalse)
1001         variance.direction[i].opacity+=density_xy[x].direction[i].opacity;
1002       if (image->colorspace == CMYKColorspace)
1003         variance.direction[i].index+=density_xy[x].direction[i].index;
1004       sum_squares.direction[i].red+=density_xy[x].direction[i].red*
1005         density_xy[x].direction[i].red;
1006       sum_squares.direction[i].green+=density_xy[x].direction[i].green*
1007         density_xy[x].direction[i].green;
1008       sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
1009         density_xy[x].direction[i].blue;
1010       if (image->matte != MagickFalse)
1011         sum_squares.direction[i].opacity+=density_xy[x].direction[i].opacity*
1012           density_xy[x].direction[i].opacity;
1013       if (image->colorspace == CMYKColorspace)
1014         sum_squares.direction[i].index+=density_xy[x].direction[i].index*
1015           density_xy[x].direction[i].index;
1016       /*
1017         Difference entropy.
1018       */
1019       channel_features[RedChannel].difference_entropy[i]-=
1020         density_xy[x].direction[i].red*
1021         log10(density_xy[x].direction[i].red+MagickEpsilon);
1022       channel_features[GreenChannel].difference_entropy[i]-=
1023         density_xy[x].direction[i].green*
1024         log10(density_xy[x].direction[i].green+MagickEpsilon);
1025       channel_features[BlueChannel].difference_entropy[i]-=
1026         density_xy[x].direction[i].blue*
1027         log10(density_xy[x].direction[i].blue+MagickEpsilon);
1028       if (image->matte != MagickFalse)
1029         channel_features[OpacityChannel].difference_entropy[i]-=
1030           density_xy[x].direction[i].opacity*
1031           log10(density_xy[x].direction[i].opacity+MagickEpsilon);
1032       if (image->colorspace == CMYKColorspace)
1033         channel_features[IndexChannel].difference_entropy[i]-=
1034           density_xy[x].direction[i].index*
1035           log10(density_xy[x].direction[i].index+MagickEpsilon);
1036       /*
1037         Information Measures of Correlation.
1038       */
1039       entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1040         log10(density_x[x].direction[i].red+MagickEpsilon));
1041       entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1042         log10(density_x[x].direction[i].green+MagickEpsilon));
1043       entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1044         log10(density_x[x].direction[i].blue+MagickEpsilon));
1045       if (image->matte != MagickFalse)
1046         entropy_x.direction[i].opacity-=(density_x[x].direction[i].opacity*
1047           log10(density_x[x].direction[i].opacity+MagickEpsilon));
1048       if (image->colorspace == CMYKColorspace)
1049         entropy_x.direction[i].index-=(density_x[x].direction[i].index*
1050           log10(density_x[x].direction[i].index+MagickEpsilon));
1051       entropy_y.direction[i].red-=(density_y[y].direction[i].red*
1052         log10(density_y[y].direction[i].red+MagickEpsilon));
1053       entropy_y.direction[i].green-=(density_y[y].direction[i].green*
1054         log10(density_y[y].direction[i].green+MagickEpsilon));
1055       entropy_y.direction[i].blue-=(density_y[y].direction[i].blue*
1056         log10(density_y[y].direction[i].blue+MagickEpsilon));
1057       if (image->matte != MagickFalse)
1058         entropy_y.direction[i].opacity-=(density_y[y].direction[i].opacity*
1059           log10(density_y[y].direction[i].opacity+MagickEpsilon));
1060       if (image->colorspace == CMYKColorspace)
1061         entropy_y.direction[i].index-=(density_y[y].direction[i].index*
1062           log10(density_y[y].direction[i].index+MagickEpsilon));
1063     }
1064     /*
1065       Difference variance.
1066     */
1067     channel_features[RedChannel].difference_variance[i]=
1068       (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1069       (variance.direction[i].red*variance.direction[i].red))/
1070       ((double) number_grays*number_grays*number_grays*number_grays);
1071     channel_features[GreenChannel].difference_variance[i]=
1072       (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1073       (variance.direction[i].green*variance.direction[i].green))/
1074       ((double) number_grays*number_grays*number_grays*number_grays);
1075     channel_features[BlueChannel].difference_variance[i]=
1076       (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1077       (variance.direction[i].blue*variance.direction[i].blue))/
1078       ((double) number_grays*number_grays*number_grays*number_grays);
1079     if (image->matte != MagickFalse)
1080       channel_features[OpacityChannel].difference_variance[i]=
1081         (((double) number_grays*number_grays*sum_squares.direction[i].opacity)-
1082         (variance.direction[i].opacity*variance.direction[i].opacity))/
1083         ((double) number_grays*number_grays*number_grays*number_grays);
1084     if (image->colorspace == CMYKColorspace)
1085       channel_features[IndexChannel].difference_variance[i]=
1086         (((double) number_grays*number_grays*sum_squares.direction[i].index)-
1087         (variance.direction[i].index*variance.direction[i].index))/
1088         ((double) number_grays*number_grays*number_grays*number_grays);
1089     /*
1090       Information Measures of Correlation.
1091     */
1092     channel_features[RedChannel].measure_of_correlation_1[i]=
1093       (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1094       (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1095        entropy_x.direction[i].red : entropy_y.direction[i].red);
1096     channel_features[GreenChannel].measure_of_correlation_1[i]=
1097       (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1098       (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1099        entropy_x.direction[i].green : entropy_y.direction[i].green);
1100     channel_features[BlueChannel].measure_of_correlation_1[i]=
1101       (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1102       (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1103        entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1104     if (image->matte != MagickFalse)
1105       channel_features[OpacityChannel].measure_of_correlation_1[i]=
1106         (entropy_xy.direction[i].opacity-entropy_xy1.direction[i].opacity)/
1107         (entropy_x.direction[i].opacity > entropy_y.direction[i].opacity ?
1108          entropy_x.direction[i].opacity : entropy_y.direction[i].opacity);
1109     if (image->colorspace == CMYKColorspace)
1110       channel_features[IndexChannel].measure_of_correlation_1[i]=
1111         (entropy_xy.direction[i].index-entropy_xy1.direction[i].index)/
1112         (entropy_x.direction[i].index > entropy_y.direction[i].index ?
1113          entropy_x.direction[i].index : entropy_y.direction[i].index);
1114     channel_features[RedChannel].measure_of_correlation_2[i]=
1115       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].red-
1116       entropy_xy.direction[i].red)))));
1117     channel_features[GreenChannel].measure_of_correlation_2[i]=
1118       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].green-
1119       entropy_xy.direction[i].green)))));
1120     channel_features[BlueChannel].measure_of_correlation_2[i]=
1121       (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].blue-
1122       entropy_xy.direction[i].blue)))));
1123     if (image->matte != MagickFalse)
1124       channel_features[OpacityChannel].measure_of_correlation_2[i]=
1125         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].opacity-
1126         entropy_xy.direction[i].opacity)))));
1127     if (image->colorspace == CMYKColorspace)
1128       channel_features[IndexChannel].measure_of_correlation_2[i]=
1129         (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].index-
1130         entropy_xy.direction[i].index)))));
1131   }
1132   /*
1133     Compute more texture features.
1134   */
1135 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1136   #pragma omp parallel for schedule(dynamic,4) shared(status)
1137 #endif
1138   for (i=0; i < 4; i++)
1139   {
1140     for (z=0; z < (ssize_t) number_grays; z++)
1141     {
1142       register ssize_t
1143         y;
1144
1145       ChannelStatistics
1146         pixel;
1147
1148       (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
1149       for (y=0; y < (ssize_t) number_grays; y++)
1150       {
1151         register ssize_t
1152           x;
1153
1154         for (x=0; x < (ssize_t) number_grays; x++)
1155         {
1156           /*
1157             Contrast:  amount of local variations present in an image.
1158           */
1159           if (((y-x) == z) || ((x-y) == z))
1160             {
1161               pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1162               pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1163               pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1164               if (image->matte != MagickFalse)
1165                 pixel.direction[i].opacity+=
1166                   cooccurrence[x][y].direction[i].opacity;
1167               if (image->colorspace == CMYKColorspace)
1168                 pixel.direction[i].index+=cooccurrence[x][y].direction[i].index;
1169             }
1170           /*
1171             Maximum Correlation Coefficient.
1172           */
1173           Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1174             cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1175             density_y[x].direction[i].red;
1176           Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1177             cooccurrence[y][x].direction[i].green/
1178             density_x[z].direction[i].green/density_y[x].direction[i].red;
1179           Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1180             cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/
1181             density_y[x].direction[i].blue;
1182           if (image->matte != MagickFalse)
1183             Q[z][y].direction[i].opacity+=
1184               cooccurrence[z][x].direction[i].opacity*
1185               cooccurrence[y][x].direction[i].opacity/
1186               density_x[z].direction[i].opacity/
1187               density_y[x].direction[i].opacity;
1188           if (image->colorspace == CMYKColorspace)
1189             Q[z][y].direction[i].index+=cooccurrence[z][x].direction[i].index*
1190               cooccurrence[y][x].direction[i].index/
1191               density_x[z].direction[i].index/density_y[x].direction[i].index;
1192         }
1193       }
1194       channel_features[RedChannel].contrast[i]+=z*z*pixel.direction[i].red;
1195       channel_features[GreenChannel].contrast[i]+=z*z*pixel.direction[i].green;
1196       channel_features[BlueChannel].contrast[i]+=z*z*pixel.direction[i].blue;
1197       if (image->matte != MagickFalse)
1198         channel_features[OpacityChannel].contrast[i]+=z*z*
1199           pixel.direction[i].opacity;
1200       if (image->colorspace == CMYKColorspace)
1201         channel_features[BlackChannel].contrast[i]+=z*z*
1202           pixel.direction[i].index;
1203     }
1204     /*
1205       Maximum Correlation Coefficient.
1206       Future: return second largest eigenvalue of Q.
1207     */
1208     channel_features[RedChannel].maximum_correlation_coefficient[i]=
1209       sqrt((double) -1.0);
1210     channel_features[GreenChannel].maximum_correlation_coefficient[i]=
1211       sqrt((double) -1.0);
1212     channel_features[BlueChannel].maximum_correlation_coefficient[i]=
1213       sqrt((double) -1.0);
1214     if (image->matte != MagickFalse)
1215       channel_features[OpacityChannel].maximum_correlation_coefficient[i]=
1216         sqrt((double) -1.0);
1217     if (image->colorspace == CMYKColorspace)
1218       channel_features[IndexChannel].maximum_correlation_coefficient[i]=
1219         sqrt((double) -1.0);
1220   }
1221   /*
1222     Relinquish resources.
1223   */
1224   sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1225   for (i=0; i < (ssize_t) number_grays; i++)
1226     Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1227   Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1228   density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1229   density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1230   density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1231   for (i=0; i < (ssize_t) number_grays; i++)
1232     cooccurrence[i]=(ChannelStatistics *)
1233       RelinquishMagickMemory(cooccurrence[i]);
1234   cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1235   return(channel_features);
1236 }