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