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