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