MagickCore  7.1.0
Convert, Edit, Or Compose Bitmap Images
feature.c
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % FFFFF EEEEE AAA TTTTT U U RRRR EEEEE %
7 % F E A A T U U R R E %
8 % FFF EEE AAAAA T U U RRRR EEE %
9 % F E A A T U U R R E %
10 % F EEEEE A A T UUU R R EEEEE %
11 % %
12 % %
13 % MagickCore Image Feature Methods %
14 % %
15 % Software Design %
16 % Cristy %
17 % July 1992 %
18 % %
19 % %
20 % Copyright @ 1999 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 % https://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/animate.h"
45 #include "MagickCore/artifact.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/channel.h"
52 #include "MagickCore/client.h"
53 #include "MagickCore/color.h"
54 #include "MagickCore/color-private.h"
55 #include "MagickCore/colorspace.h"
56 #include "MagickCore/colorspace-private.h"
57 #include "MagickCore/composite.h"
58 #include "MagickCore/composite-private.h"
59 #include "MagickCore/compress.h"
60 #include "MagickCore/constitute.h"
61 #include "MagickCore/display.h"
62 #include "MagickCore/draw.h"
63 #include "MagickCore/enhance.h"
64 #include "MagickCore/exception.h"
65 #include "MagickCore/exception-private.h"
66 #include "MagickCore/feature.h"
67 #include "MagickCore/gem.h"
68 #include "MagickCore/geometry.h"
69 #include "MagickCore/list.h"
70 #include "MagickCore/image-private.h"
71 #include "MagickCore/magic.h"
72 #include "MagickCore/magick.h"
73 #include "MagickCore/matrix.h"
74 #include "MagickCore/memory_.h"
75 #include "MagickCore/module.h"
76 #include "MagickCore/monitor.h"
77 #include "MagickCore/monitor-private.h"
78 #include "MagickCore/morphology-private.h"
79 #include "MagickCore/option.h"
80 #include "MagickCore/paint.h"
81 #include "MagickCore/pixel-accessor.h"
82 #include "MagickCore/profile.h"
83 #include "MagickCore/property.h"
84 #include "MagickCore/quantize.h"
85 #include "MagickCore/quantum-private.h"
86 #include "MagickCore/random_.h"
87 #include "MagickCore/resource_.h"
88 #include "MagickCore/segment.h"
89 #include "MagickCore/semaphore.h"
90 #include "MagickCore/signature-private.h"
91 #include "MagickCore/string_.h"
92 #include "MagickCore/thread-private.h"
93 #include "MagickCore/timer.h"
94 #include "MagickCore/utility.h"
95 #include "MagickCore/version.h"
96 ␌
97 /*
98 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99 % %
100 % %
101 % %
102 % C a n n y E d g e I m a g e %
103 % %
104 % %
105 % %
106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107 %
108 % CannyEdgeImage() uses a multi-stage algorithm to detect a wide range of
109 % edges in images.
110 %
111 % The format of the CannyEdgeImage method is:
112 %
113 % Image *CannyEdgeImage(const Image *image,const double radius,
114 % const double sigma,const double lower_percent,
115 % const double upper_percent,ExceptionInfo *exception)
116 %
117 % A description of each parameter follows:
118 %
119 % o image: the image.
120 %
121 % o radius: the radius of the gaussian smoothing filter.
122 %
123 % o sigma: the sigma of the gaussian smoothing filter.
124 %
125 % o lower_percent: percentage of edge pixels in the lower threshold.
126 %
127 % o upper_percent: percentage of edge pixels in the upper threshold.
128 %
129 % o exception: return any errors or warnings in this structure.
130 %
131 */
132 
133 typedef struct _CannyInfo
134 {
135  double
136  magnitude,
137  intensity;
138 
139  int
140  orientation;
141 
142  ssize_t
143  x,
144  y;
145 } CannyInfo;
146 
147 static inline MagickBooleanType IsAuthenticPixel(const Image *image,
148  const ssize_t x,const ssize_t y)
149 {
150  if ((x < 0) || (x >= (ssize_t) image->columns))
151  return(MagickFalse);
152  if ((y < 0) || (y >= (ssize_t) image->rows))
153  return(MagickFalse);
154  return(MagickTrue);
155 }
156 
157 static MagickBooleanType TraceEdges(Image *edge_image,CacheView *edge_view,
158  MatrixInfo *canny_cache,const ssize_t x,const ssize_t y,
159  const double lower_threshold,ExceptionInfo *exception)
160 {
161  CannyInfo
162  edge,
163  pixel;
164 
165  MagickBooleanType
166  status;
167 
168  Quantum
169  *q;
170 
171  ssize_t
172  i;
173 
174  q=GetCacheViewAuthenticPixels(edge_view,x,y,1,1,exception);
175  if (q == (Quantum *) NULL)
176  return(MagickFalse);
177  *q=QuantumRange;
178  status=SyncCacheViewAuthenticPixels(edge_view,exception);
179  if (status == MagickFalse)
180  return(MagickFalse);
181  if (GetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
182  return(MagickFalse);
183  edge.x=x;
184  edge.y=y;
185  if (SetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
186  return(MagickFalse);
187  for (i=1; i != 0; )
188  {
189  ssize_t
190  v;
191 
192  i--;
193  status=GetMatrixElement(canny_cache,i,0,&edge);
194  if (status == MagickFalse)
195  return(MagickFalse);
196  for (v=(-1); v <= 1; v++)
197  {
198  ssize_t
199  u;
200 
201  for (u=(-1); u <= 1; u++)
202  {
203  if ((u == 0) && (v == 0))
204  continue;
205  if (IsAuthenticPixel(edge_image,edge.x+u,edge.y+v) == MagickFalse)
206  continue;
207  /*
208  Not an edge if gradient value is below the lower threshold.
209  */
210  q=GetCacheViewAuthenticPixels(edge_view,edge.x+u,edge.y+v,1,1,
211  exception);
212  if (q == (Quantum *) NULL)
213  return(MagickFalse);
214  status=GetMatrixElement(canny_cache,edge.x+u,edge.y+v,&pixel);
215  if (status == MagickFalse)
216  return(MagickFalse);
217  if ((GetPixelIntensity(edge_image,q) == 0.0) &&
218  (pixel.intensity >= lower_threshold))
219  {
220  *q=QuantumRange;
221  status=SyncCacheViewAuthenticPixels(edge_view,exception);
222  if (status == MagickFalse)
223  return(MagickFalse);
224  edge.x+=u;
225  edge.y+=v;
226  status=SetMatrixElement(canny_cache,i,0,&edge);
227  if (status == MagickFalse)
228  return(MagickFalse);
229  i++;
230  }
231  }
232  }
233  }
234  return(MagickTrue);
235 }
236 
237 MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
238  const double sigma,const double lower_percent,const double upper_percent,
239  ExceptionInfo *exception)
240 {
241 #define CannyEdgeImageTag "CannyEdge/Image"
242 
243  CacheView
244  *edge_view;
245 
246  CannyInfo
247  element;
248 
249  char
250  geometry[MagickPathExtent];
251 
252  double
253  lower_threshold,
254  max,
255  min,
256  upper_threshold;
257 
258  Image
259  *edge_image;
260 
261  KernelInfo
262  *kernel_info;
263 
264  MagickBooleanType
265  status;
266 
267  MagickOffsetType
268  progress;
269 
270  MatrixInfo
271  *canny_cache;
272 
273  ssize_t
274  y;
275 
276  assert(image != (const Image *) NULL);
277  assert(image->signature == MagickCoreSignature);
278  assert(exception != (ExceptionInfo *) NULL);
279  assert(exception->signature == MagickCoreSignature);
280  if (IsEventLogging() != MagickFalse)
281  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
282  /*
283  Filter out noise.
284  */
285  (void) FormatLocaleString(geometry,MagickPathExtent,
286  "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
287  kernel_info=AcquireKernelInfo(geometry,exception);
288  if (kernel_info == (KernelInfo *) NULL)
289  ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
290  edge_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,exception);
291  kernel_info=DestroyKernelInfo(kernel_info);
292  if (edge_image == (Image *) NULL)
293  return((Image *) NULL);
294  if (TransformImageColorspace(edge_image,GRAYColorspace,exception) == MagickFalse)
295  {
296  edge_image=DestroyImage(edge_image);
297  return((Image *) NULL);
298  }
299  (void) SetImageAlphaChannel(edge_image,OffAlphaChannel,exception);
300  /*
301  Find the intensity gradient of the image.
302  */
303  canny_cache=AcquireMatrixInfo(edge_image->columns,edge_image->rows,
304  sizeof(CannyInfo),exception);
305  if (canny_cache == (MatrixInfo *) NULL)
306  {
307  edge_image=DestroyImage(edge_image);
308  return((Image *) NULL);
309  }
310  status=MagickTrue;
311  edge_view=AcquireVirtualCacheView(edge_image,exception);
312 #if defined(MAGICKCORE_OPENMP_SUPPORT)
313  #pragma omp parallel for schedule(static) shared(status) \
314  magick_number_threads(edge_image,edge_image,edge_image->rows,1)
315 #endif
316  for (y=0; y < (ssize_t) edge_image->rows; y++)
317  {
318  const Quantum
319  *magick_restrict p;
320 
321  ssize_t
322  x;
323 
324  if (status == MagickFalse)
325  continue;
326  p=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns+1,2,
327  exception);
328  if (p == (const Quantum *) NULL)
329  {
330  status=MagickFalse;
331  continue;
332  }
333  for (x=0; x < (ssize_t) edge_image->columns; x++)
334  {
335  CannyInfo
336  pixel;
337 
338  double
339  dx,
340  dy;
341 
342  const Quantum
343  *magick_restrict kernel_pixels;
344 
345  ssize_t
346  v;
347 
348  static double
349  Gx[2][2] =
350  {
351  { -1.0, +1.0 },
352  { -1.0, +1.0 }
353  },
354  Gy[2][2] =
355  {
356  { +1.0, +1.0 },
357  { -1.0, -1.0 }
358  };
359 
360  (void) memset(&pixel,0,sizeof(pixel));
361  dx=0.0;
362  dy=0.0;
363  kernel_pixels=p;
364  for (v=0; v < 2; v++)
365  {
366  ssize_t
367  u;
368 
369  for (u=0; u < 2; u++)
370  {
371  double
372  intensity;
373 
374  intensity=GetPixelIntensity(edge_image,kernel_pixels+u);
375  dx+=0.5*Gx[v][u]*intensity;
376  dy+=0.5*Gy[v][u]*intensity;
377  }
378  kernel_pixels+=edge_image->columns+1;
379  }
380  pixel.magnitude=hypot(dx,dy);
381  pixel.orientation=0;
382  if (fabs(dx) > MagickEpsilon)
383  {
384  double
385  slope;
386 
387  slope=dy/dx;
388  if (slope < 0.0)
389  {
390  if (slope < -2.41421356237)
391  pixel.orientation=0;
392  else
393  if (slope < -0.414213562373)
394  pixel.orientation=1;
395  else
396  pixel.orientation=2;
397  }
398  else
399  {
400  if (slope > 2.41421356237)
401  pixel.orientation=0;
402  else
403  if (slope > 0.414213562373)
404  pixel.orientation=3;
405  else
406  pixel.orientation=2;
407  }
408  }
409  if (SetMatrixElement(canny_cache,x,y,&pixel) == MagickFalse)
410  continue;
411  p+=GetPixelChannels(edge_image);
412  }
413  }
414  edge_view=DestroyCacheView(edge_view);
415  /*
416  Non-maxima suppression, remove pixels that are not considered to be part
417  of an edge.
418  */
419  progress=0;
420  (void) GetMatrixElement(canny_cache,0,0,&element);
421  max=element.intensity;
422  min=element.intensity;
423  edge_view=AcquireAuthenticCacheView(edge_image,exception);
424 #if defined(MAGICKCORE_OPENMP_SUPPORT)
425  #pragma omp parallel for schedule(static) shared(status) \
426  magick_number_threads(edge_image,edge_image,edge_image->rows,1)
427 #endif
428  for (y=0; y < (ssize_t) edge_image->rows; y++)
429  {
430  Quantum
431  *magick_restrict q;
432 
433  ssize_t
434  x;
435 
436  if (status == MagickFalse)
437  continue;
438  q=GetCacheViewAuthenticPixels(edge_view,0,y,edge_image->columns,1,
439  exception);
440  if (q == (Quantum *) NULL)
441  {
442  status=MagickFalse;
443  continue;
444  }
445  for (x=0; x < (ssize_t) edge_image->columns; x++)
446  {
447  CannyInfo
448  alpha_pixel,
449  beta_pixel,
450  pixel;
451 
452  (void) GetMatrixElement(canny_cache,x,y,&pixel);
453  switch (pixel.orientation)
454  {
455  case 0:
456  default:
457  {
458  /*
459  0 degrees, north and south.
460  */
461  (void) GetMatrixElement(canny_cache,x,y-1,&alpha_pixel);
462  (void) GetMatrixElement(canny_cache,x,y+1,&beta_pixel);
463  break;
464  }
465  case 1:
466  {
467  /*
468  45 degrees, northwest and southeast.
469  */
470  (void) GetMatrixElement(canny_cache,x-1,y-1,&alpha_pixel);
471  (void) GetMatrixElement(canny_cache,x+1,y+1,&beta_pixel);
472  break;
473  }
474  case 2:
475  {
476  /*
477  90 degrees, east and west.
478  */
479  (void) GetMatrixElement(canny_cache,x-1,y,&alpha_pixel);
480  (void) GetMatrixElement(canny_cache,x+1,y,&beta_pixel);
481  break;
482  }
483  case 3:
484  {
485  /*
486  135 degrees, northeast and southwest.
487  */
488  (void) GetMatrixElement(canny_cache,x+1,y-1,&beta_pixel);
489  (void) GetMatrixElement(canny_cache,x-1,y+1,&alpha_pixel);
490  break;
491  }
492  }
493  pixel.intensity=pixel.magnitude;
494  if ((pixel.magnitude < alpha_pixel.magnitude) ||
495  (pixel.magnitude < beta_pixel.magnitude))
496  pixel.intensity=0;
497  (void) SetMatrixElement(canny_cache,x,y,&pixel);
498 #if defined(MAGICKCORE_OPENMP_SUPPORT)
499  #pragma omp critical (MagickCore_CannyEdgeImage)
500 #endif
501  {
502  if (pixel.intensity < min)
503  min=pixel.intensity;
504  if (pixel.intensity > max)
505  max=pixel.intensity;
506  }
507  *q=0;
508  q+=GetPixelChannels(edge_image);
509  }
510  if (SyncCacheViewAuthenticPixels(edge_view,exception) == MagickFalse)
511  status=MagickFalse;
512  }
513  edge_view=DestroyCacheView(edge_view);
514  /*
515  Estimate hysteresis threshold.
516  */
517  lower_threshold=lower_percent*(max-min)+min;
518  upper_threshold=upper_percent*(max-min)+min;
519  /*
520  Hysteresis threshold.
521  */
522  edge_view=AcquireAuthenticCacheView(edge_image,exception);
523  for (y=0; y < (ssize_t) edge_image->rows; y++)
524  {
525  ssize_t
526  x;
527 
528  if (status == MagickFalse)
529  continue;
530  for (x=0; x < (ssize_t) edge_image->columns; x++)
531  {
532  CannyInfo
533  pixel;
534 
535  const Quantum
536  *magick_restrict p;
537 
538  /*
539  Edge if pixel gradient higher than upper threshold.
540  */
541  p=GetCacheViewVirtualPixels(edge_view,x,y,1,1,exception);
542  if (p == (const Quantum *) NULL)
543  continue;
544  status=GetMatrixElement(canny_cache,x,y,&pixel);
545  if (status == MagickFalse)
546  continue;
547  if ((GetPixelIntensity(edge_image,p) == 0.0) &&
548  (pixel.intensity >= upper_threshold))
549  status=TraceEdges(edge_image,edge_view,canny_cache,x,y,lower_threshold,
550  exception);
551  }
552  if (image->progress_monitor != (MagickProgressMonitor) NULL)
553  {
554  MagickBooleanType
555  proceed;
556 
557 #if defined(MAGICKCORE_OPENMP_SUPPORT)
558  #pragma omp atomic
559 #endif
560  progress++;
561  proceed=SetImageProgress(image,CannyEdgeImageTag,progress,image->rows);
562  if (proceed == MagickFalse)
563  status=MagickFalse;
564  }
565  }
566  edge_view=DestroyCacheView(edge_view);
567  /*
568  Free resources.
569  */
570  canny_cache=DestroyMatrixInfo(canny_cache);
571  return(edge_image);
572 }
573 ␌
574 /*
575 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
576 % %
577 % %
578 % %
579 % G e t I m a g e F e a t u r e s %
580 % %
581 % %
582 % %
583 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
584 %
585 % GetImageFeatures() returns features for each channel in the image in
586 % each of four directions (horizontal, vertical, left and right diagonals)
587 % for the specified distance. The features include the angular second
588 % moment, contrast, correlation, sum of squares: variance, inverse difference
589 % moment, sum average, sum varience, sum entropy, entropy, difference variance,
590 % difference entropy, information measures of correlation 1, information
591 % measures of correlation 2, and maximum correlation coefficient. You can
592 % access the red channel contrast, for example, like this:
593 %
594 % channel_features=GetImageFeatures(image,1,exception);
595 % contrast=channel_features[RedPixelChannel].contrast[0];
596 %
597 % Use MagickRelinquishMemory() to free the features buffer.
598 %
599 % The format of the GetImageFeatures method is:
600 %
601 % ChannelFeatures *GetImageFeatures(const Image *image,
602 % const size_t distance,ExceptionInfo *exception)
603 %
604 % A description of each parameter follows:
605 %
606 % o image: the image.
607 %
608 % o distance: the distance.
609 %
610 % o exception: return any errors or warnings in this structure.
611 %
612 */
613 
614 static inline double MagickLog10(const double x)
615 {
616  if (fabs(x) < MagickEpsilon)
617  return(-INFINITY);
618  return(log10(fabs(x)));
619 }
620 
621 MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
622  const size_t distance,ExceptionInfo *exception)
623 {
624  typedef struct _ChannelStatistics
625  {
626  PixelInfo
627  direction[4]; /* horizontal, vertical, left and right diagonals */
629 
630  CacheView
631  *image_view;
632 
634  *channel_features;
635 
637  **cooccurrence,
638  correlation,
639  *density_x,
640  *density_xy,
641  *density_y,
642  entropy_x,
643  entropy_xy,
644  entropy_xy1,
645  entropy_xy2,
646  entropy_y,
647  mean,
648  **Q,
649  *sum,
650  sum_squares,
651  variance;
652 
654  gray,
655  *grays;
656 
657  MagickBooleanType
658  status;
659 
660  ssize_t
661  i,
662  r;
663 
664  size_t
665  length;
666 
667  unsigned int
668  number_grays;
669 
670  assert(image != (Image *) NULL);
671  assert(image->signature == MagickCoreSignature);
672  if (IsEventLogging() != MagickFalse)
673  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
674  if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
675  return((ChannelFeatures *) NULL);
676  length=MaxPixelChannels+1UL;
677  channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
678  sizeof(*channel_features));
679  if (channel_features == (ChannelFeatures *) NULL)
680  ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
681  (void) memset(channel_features,0,length*
682  sizeof(*channel_features));
683  /*
684  Form grays.
685  */
686  grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays));
687  if (grays == (PixelPacket *) NULL)
688  {
689  channel_features=(ChannelFeatures *) RelinquishMagickMemory(
690  channel_features);
691  (void) ThrowMagickException(exception,GetMagickModule(),
692  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
693  return(channel_features);
694  }
695  for (i=0; i <= (ssize_t) MaxMap; i++)
696  {
697  grays[i].red=(~0U);
698  grays[i].green=(~0U);
699  grays[i].blue=(~0U);
700  grays[i].alpha=(~0U);
701  grays[i].black=(~0U);
702  }
703  status=MagickTrue;
704  image_view=AcquireVirtualCacheView(image,exception);
705 #if defined(MAGICKCORE_OPENMP_SUPPORT)
706  #pragma omp parallel for schedule(static) shared(status) \
707  magick_number_threads(image,image,image->rows,1)
708 #endif
709  for (r=0; r < (ssize_t) image->rows; r++)
710  {
711  const Quantum
712  *magick_restrict p;
713 
714  ssize_t
715  x;
716 
717  if (status == MagickFalse)
718  continue;
719  p=GetCacheViewVirtualPixels(image_view,0,r,image->columns,1,exception);
720  if (p == (const Quantum *) NULL)
721  {
722  status=MagickFalse;
723  continue;
724  }
725  for (x=0; x < (ssize_t) image->columns; x++)
726  {
727  grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
728  ScaleQuantumToMap(GetPixelRed(image,p));
729  grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
730  ScaleQuantumToMap(GetPixelGreen(image,p));
731  grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
732  ScaleQuantumToMap(GetPixelBlue(image,p));
733  if (image->colorspace == CMYKColorspace)
734  grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
735  ScaleQuantumToMap(GetPixelBlack(image,p));
736  if (image->alpha_trait != UndefinedPixelTrait)
737  grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
738  ScaleQuantumToMap(GetPixelAlpha(image,p));
739  p+=GetPixelChannels(image);
740  }
741  }
742  image_view=DestroyCacheView(image_view);
743  if (status == MagickFalse)
744  {
745  grays=(PixelPacket *) RelinquishMagickMemory(grays);
746  channel_features=(ChannelFeatures *) RelinquishMagickMemory(
747  channel_features);
748  return(channel_features);
749  }
750  (void) memset(&gray,0,sizeof(gray));
751  for (i=0; i <= (ssize_t) MaxMap; i++)
752  {
753  if (grays[i].red != ~0U)
754  grays[gray.red++].red=grays[i].red;
755  if (grays[i].green != ~0U)
756  grays[gray.green++].green=grays[i].green;
757  if (grays[i].blue != ~0U)
758  grays[gray.blue++].blue=grays[i].blue;
759  if (image->colorspace == CMYKColorspace)
760  if (grays[i].black != ~0U)
761  grays[gray.black++].black=grays[i].black;
762  if (image->alpha_trait != UndefinedPixelTrait)
763  if (grays[i].alpha != ~0U)
764  grays[gray.alpha++].alpha=grays[i].alpha;
765  }
766  /*
767  Allocate spatial dependence matrix.
768  */
769  number_grays=gray.red;
770  if (gray.green > number_grays)
771  number_grays=gray.green;
772  if (gray.blue > number_grays)
773  number_grays=gray.blue;
774  if (image->colorspace == CMYKColorspace)
775  if (gray.black > number_grays)
776  number_grays=gray.black;
777  if (image->alpha_trait != UndefinedPixelTrait)
778  if (gray.alpha > number_grays)
779  number_grays=gray.alpha;
780  cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays,
781  sizeof(*cooccurrence));
782  density_x=(ChannelStatistics *) AcquireQuantumMemory(number_grays+1,
783  2*sizeof(*density_x));
784  density_xy=(ChannelStatistics *) AcquireQuantumMemory(number_grays+1,
785  2*sizeof(*density_xy));
786  density_y=(ChannelStatistics *) AcquireQuantumMemory(number_grays+1,
787  2*sizeof(*density_y));
788  Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q));
789  sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum));
790  if ((cooccurrence == (ChannelStatistics **) NULL) ||
791  (density_x == (ChannelStatistics *) NULL) ||
792  (density_xy == (ChannelStatistics *) NULL) ||
793  (density_y == (ChannelStatistics *) NULL) ||
794  (Q == (ChannelStatistics **) NULL) ||
795  (sum == (ChannelStatistics *) NULL))
796  {
797  if (Q != (ChannelStatistics **) NULL)
798  {
799  for (i=0; i < (ssize_t) number_grays; i++)
800  Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
801  Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
802  }
803  if (sum != (ChannelStatistics *) NULL)
804  sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
805  if (density_y != (ChannelStatistics *) NULL)
806  density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
807  if (density_xy != (ChannelStatistics *) NULL)
808  density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
809  if (density_x != (ChannelStatistics *) NULL)
810  density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
811  if (cooccurrence != (ChannelStatistics **) NULL)
812  {
813  for (i=0; i < (ssize_t) number_grays; i++)
814  cooccurrence[i]=(ChannelStatistics *)
815  RelinquishMagickMemory(cooccurrence[i]);
816  cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
817  cooccurrence);
818  }
819  grays=(PixelPacket *) RelinquishMagickMemory(grays);
820  channel_features=(ChannelFeatures *) RelinquishMagickMemory(
821  channel_features);
822  (void) ThrowMagickException(exception,GetMagickModule(),
823  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
824  return(channel_features);
825  }
826  (void) memset(&correlation,0,sizeof(correlation));
827  (void) memset(density_x,0,2*(number_grays+1)*sizeof(*density_x));
828  (void) memset(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
829  (void) memset(density_y,0,2*(number_grays+1)*sizeof(*density_y));
830  (void) memset(&mean,0,sizeof(mean));
831  (void) memset(sum,0,number_grays*sizeof(*sum));
832  (void) memset(&sum_squares,0,sizeof(sum_squares));
833  (void) memset(density_xy,0,2*number_grays*sizeof(*density_xy));
834  (void) memset(&entropy_x,0,sizeof(entropy_x));
835  (void) memset(&entropy_xy,0,sizeof(entropy_xy));
836  (void) memset(&entropy_xy1,0,sizeof(entropy_xy1));
837  (void) memset(&entropy_xy2,0,sizeof(entropy_xy2));
838  (void) memset(&entropy_y,0,sizeof(entropy_y));
839  (void) memset(&variance,0,sizeof(variance));
840  for (i=0; i < (ssize_t) number_grays; i++)
841  {
842  cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
843  sizeof(**cooccurrence));
844  Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
845  if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
846  (Q[i] == (ChannelStatistics *) NULL))
847  break;
848  (void) memset(cooccurrence[i],0,number_grays*
849  sizeof(**cooccurrence));
850  (void) memset(Q[i],0,number_grays*sizeof(**Q));
851  }
852  if (i < (ssize_t) number_grays)
853  {
854  for (i--; i >= 0; i--)
855  {
856  if (Q[i] != (ChannelStatistics *) NULL)
857  Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
858  if (cooccurrence[i] != (ChannelStatistics *) NULL)
859  cooccurrence[i]=(ChannelStatistics *)
860  RelinquishMagickMemory(cooccurrence[i]);
861  }
862  Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
863  cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
864  sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
865  density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
866  density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
867  density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
868  grays=(PixelPacket *) RelinquishMagickMemory(grays);
869  channel_features=(ChannelFeatures *) RelinquishMagickMemory(
870  channel_features);
871  (void) ThrowMagickException(exception,GetMagickModule(),
872  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
873  return(channel_features);
874  }
875  /*
876  Initialize spatial dependence matrix.
877  */
878  status=MagickTrue;
879  image_view=AcquireVirtualCacheView(image,exception);
880  for (r=0; r < (ssize_t) image->rows; r++)
881  {
882  const Quantum
883  *magick_restrict p;
884 
885  ssize_t
886  x;
887 
888  ssize_t
889  offset,
890  u,
891  v;
892 
893  if (status == MagickFalse)
894  continue;
895  p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,r,image->columns+
896  2*distance,distance+2,exception);
897  if (p == (const Quantum *) NULL)
898  {
899  status=MagickFalse;
900  continue;
901  }
902  p+=distance*GetPixelChannels(image);;
903  for (x=0; x < (ssize_t) image->columns; x++)
904  {
905  for (i=0; i < 4; i++)
906  {
907  switch (i)
908  {
909  case 0:
910  default:
911  {
912  /*
913  Horizontal adjacency.
914  */
915  offset=(ssize_t) distance;
916  break;
917  }
918  case 1:
919  {
920  /*
921  Vertical adjacency.
922  */
923  offset=(ssize_t) (image->columns+2*distance);
924  break;
925  }
926  case 2:
927  {
928  /*
929  Right diagonal adjacency.
930  */
931  offset=(ssize_t) ((image->columns+2*distance)-distance);
932  break;
933  }
934  case 3:
935  {
936  /*
937  Left diagonal adjacency.
938  */
939  offset=(ssize_t) ((image->columns+2*distance)+distance);
940  break;
941  }
942  }
943  u=0;
944  v=0;
945  while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p)))
946  u++;
947  while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image))))
948  v++;
949  cooccurrence[u][v].direction[i].red++;
950  cooccurrence[v][u].direction[i].red++;
951  u=0;
952  v=0;
953  while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p)))
954  u++;
955  while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image))))
956  v++;
957  cooccurrence[u][v].direction[i].green++;
958  cooccurrence[v][u].direction[i].green++;
959  u=0;
960  v=0;
961  while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p)))
962  u++;
963  while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image))))
964  v++;
965  cooccurrence[u][v].direction[i].blue++;
966  cooccurrence[v][u].direction[i].blue++;
967  if (image->colorspace == CMYKColorspace)
968  {
969  u=0;
970  v=0;
971  while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p)))
972  u++;
973  while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image))))
974  v++;
975  cooccurrence[u][v].direction[i].black++;
976  cooccurrence[v][u].direction[i].black++;
977  }
978  if (image->alpha_trait != UndefinedPixelTrait)
979  {
980  u=0;
981  v=0;
982  while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p)))
983  u++;
984  while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image))))
985  v++;
986  cooccurrence[u][v].direction[i].alpha++;
987  cooccurrence[v][u].direction[i].alpha++;
988  }
989  }
990  p+=GetPixelChannels(image);
991  }
992  }
993  grays=(PixelPacket *) RelinquishMagickMemory(grays);
994  image_view=DestroyCacheView(image_view);
995  if (status == MagickFalse)
996  {
997  for (i=0; i < (ssize_t) number_grays; i++)
998  cooccurrence[i]=(ChannelStatistics *)
999  RelinquishMagickMemory(cooccurrence[i]);
1000  cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1001  channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1002  channel_features);
1003  (void) ThrowMagickException(exception,GetMagickModule(),
1004  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1005  return(channel_features);
1006  }
1007  /*
1008  Normalize spatial dependence matrix.
1009  */
1010  for (i=0; i < 4; i++)
1011  {
1012  double
1013  normalize;
1014 
1015  ssize_t
1016  y;
1017 
1018  switch (i)
1019  {
1020  case 0:
1021  default:
1022  {
1023  /*
1024  Horizontal adjacency.
1025  */
1026  normalize=2.0*image->rows*(image->columns-distance);
1027  break;
1028  }
1029  case 1:
1030  {
1031  /*
1032  Vertical adjacency.
1033  */
1034  normalize=2.0*(image->rows-distance)*image->columns;
1035  break;
1036  }
1037  case 2:
1038  {
1039  /*
1040  Right diagonal adjacency.
1041  */
1042  normalize=2.0*(image->rows-distance)*(image->columns-distance);
1043  break;
1044  }
1045  case 3:
1046  {
1047  /*
1048  Left diagonal adjacency.
1049  */
1050  normalize=2.0*(image->rows-distance)*(image->columns-distance);
1051  break;
1052  }
1053  }
1054  normalize=PerceptibleReciprocal(normalize);
1055  for (y=0; y < (ssize_t) number_grays; y++)
1056  {
1057  ssize_t
1058  x;
1059 
1060  for (x=0; x < (ssize_t) number_grays; x++)
1061  {
1062  cooccurrence[x][y].direction[i].red*=normalize;
1063  cooccurrence[x][y].direction[i].green*=normalize;
1064  cooccurrence[x][y].direction[i].blue*=normalize;
1065  if (image->colorspace == CMYKColorspace)
1066  cooccurrence[x][y].direction[i].black*=normalize;
1067  if (image->alpha_trait != UndefinedPixelTrait)
1068  cooccurrence[x][y].direction[i].alpha*=normalize;
1069  }
1070  }
1071  }
1072  /*
1073  Compute texture features.
1074  */
1075 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1076  #pragma omp parallel for schedule(static) shared(status) \
1077  magick_number_threads(image,image,number_grays,1)
1078 #endif
1079  for (i=0; i < 4; i++)
1080  {
1081  ssize_t
1082  y;
1083 
1084  for (y=0; y < (ssize_t) number_grays; y++)
1085  {
1086  ssize_t
1087  x;
1088 
1089  for (x=0; x < (ssize_t) number_grays; x++)
1090  {
1091  /*
1092  Angular second moment: measure of homogeneity of the image.
1093  */
1094  channel_features[RedPixelChannel].angular_second_moment[i]+=
1095  cooccurrence[x][y].direction[i].red*
1096  cooccurrence[x][y].direction[i].red;
1097  channel_features[GreenPixelChannel].angular_second_moment[i]+=
1098  cooccurrence[x][y].direction[i].green*
1099  cooccurrence[x][y].direction[i].green;
1100  channel_features[BluePixelChannel].angular_second_moment[i]+=
1101  cooccurrence[x][y].direction[i].blue*
1102  cooccurrence[x][y].direction[i].blue;
1103  if (image->colorspace == CMYKColorspace)
1104  channel_features[BlackPixelChannel].angular_second_moment[i]+=
1105  cooccurrence[x][y].direction[i].black*
1106  cooccurrence[x][y].direction[i].black;
1107  if (image->alpha_trait != UndefinedPixelTrait)
1108  channel_features[AlphaPixelChannel].angular_second_moment[i]+=
1109  cooccurrence[x][y].direction[i].alpha*
1110  cooccurrence[x][y].direction[i].alpha;
1111  /*
1112  Correlation: measure of linear-dependencies in the image.
1113  */
1114  sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1115  sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1116  sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1117  if (image->colorspace == CMYKColorspace)
1118  sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black;
1119  if (image->alpha_trait != UndefinedPixelTrait)
1120  sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha;
1121  correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
1122  correlation.direction[i].green+=x*y*
1123  cooccurrence[x][y].direction[i].green;
1124  correlation.direction[i].blue+=x*y*
1125  cooccurrence[x][y].direction[i].blue;
1126  if (image->colorspace == CMYKColorspace)
1127  correlation.direction[i].black+=x*y*
1128  cooccurrence[x][y].direction[i].black;
1129  if (image->alpha_trait != UndefinedPixelTrait)
1130  correlation.direction[i].alpha+=x*y*
1131  cooccurrence[x][y].direction[i].alpha;
1132  /*
1133  Inverse Difference Moment.
1134  */
1135  channel_features[RedPixelChannel].inverse_difference_moment[i]+=
1136  cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
1137  channel_features[GreenPixelChannel].inverse_difference_moment[i]+=
1138  cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
1139  channel_features[BluePixelChannel].inverse_difference_moment[i]+=
1140  cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
1141  if (image->colorspace == CMYKColorspace)
1142  channel_features[BlackPixelChannel].inverse_difference_moment[i]+=
1143  cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1);
1144  if (image->alpha_trait != UndefinedPixelTrait)
1145  channel_features[AlphaPixelChannel].inverse_difference_moment[i]+=
1146  cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
1147  /*
1148  Sum average.
1149  */
1150  density_xy[y+x+2].direction[i].red+=
1151  cooccurrence[x][y].direction[i].red;
1152  density_xy[y+x+2].direction[i].green+=
1153  cooccurrence[x][y].direction[i].green;
1154  density_xy[y+x+2].direction[i].blue+=
1155  cooccurrence[x][y].direction[i].blue;
1156  if (image->colorspace == CMYKColorspace)
1157  density_xy[y+x+2].direction[i].black+=
1158  cooccurrence[x][y].direction[i].black;
1159  if (image->alpha_trait != UndefinedPixelTrait)
1160  density_xy[y+x+2].direction[i].alpha+=
1161  cooccurrence[x][y].direction[i].alpha;
1162  /*
1163  Entropy.
1164  */
1165  channel_features[RedPixelChannel].entropy[i]-=
1166  cooccurrence[x][y].direction[i].red*
1167  MagickLog10(cooccurrence[x][y].direction[i].red);
1168  channel_features[GreenPixelChannel].entropy[i]-=
1169  cooccurrence[x][y].direction[i].green*
1170  MagickLog10(cooccurrence[x][y].direction[i].green);
1171  channel_features[BluePixelChannel].entropy[i]-=
1172  cooccurrence[x][y].direction[i].blue*
1173  MagickLog10(cooccurrence[x][y].direction[i].blue);
1174  if (image->colorspace == CMYKColorspace)
1175  channel_features[BlackPixelChannel].entropy[i]-=
1176  cooccurrence[x][y].direction[i].black*
1177  MagickLog10(cooccurrence[x][y].direction[i].black);
1178  if (image->alpha_trait != UndefinedPixelTrait)
1179  channel_features[AlphaPixelChannel].entropy[i]-=
1180  cooccurrence[x][y].direction[i].alpha*
1181  MagickLog10(cooccurrence[x][y].direction[i].alpha);
1182  /*
1183  Information Measures of Correlation.
1184  */
1185  density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
1186  density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
1187  density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1188  if (image->alpha_trait != UndefinedPixelTrait)
1189  density_x[x].direction[i].alpha+=
1190  cooccurrence[x][y].direction[i].alpha;
1191  if (image->colorspace == CMYKColorspace)
1192  density_x[x].direction[i].black+=
1193  cooccurrence[x][y].direction[i].black;
1194  density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1195  density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1196  density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1197  if (image->colorspace == CMYKColorspace)
1198  density_y[y].direction[i].black+=
1199  cooccurrence[x][y].direction[i].black;
1200  if (image->alpha_trait != UndefinedPixelTrait)
1201  density_y[y].direction[i].alpha+=
1202  cooccurrence[x][y].direction[i].alpha;
1203  }
1204  mean.direction[i].red+=y*sum[y].direction[i].red;
1205  sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
1206  mean.direction[i].green+=y*sum[y].direction[i].green;
1207  sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
1208  mean.direction[i].blue+=y*sum[y].direction[i].blue;
1209  sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
1210  if (image->colorspace == CMYKColorspace)
1211  {
1212  mean.direction[i].black+=y*sum[y].direction[i].black;
1213  sum_squares.direction[i].black+=y*y*sum[y].direction[i].black;
1214  }
1215  if (image->alpha_trait != UndefinedPixelTrait)
1216  {
1217  mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
1218  sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
1219  }
1220  }
1221  /*
1222  Correlation: measure of linear-dependencies in the image.
1223  */
1224  channel_features[RedPixelChannel].correlation[i]=
1225  (correlation.direction[i].red-mean.direction[i].red*
1226  mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
1227  (mean.direction[i].red*mean.direction[i].red))*sqrt(
1228  sum_squares.direction[i].red-(mean.direction[i].red*
1229  mean.direction[i].red)));
1230  channel_features[GreenPixelChannel].correlation[i]=
1231  (correlation.direction[i].green-mean.direction[i].green*
1232  mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
1233  (mean.direction[i].green*mean.direction[i].green))*sqrt(
1234  sum_squares.direction[i].green-(mean.direction[i].green*
1235  mean.direction[i].green)));
1236  channel_features[BluePixelChannel].correlation[i]=
1237  (correlation.direction[i].blue-mean.direction[i].blue*
1238  mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
1239  (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
1240  sum_squares.direction[i].blue-(mean.direction[i].blue*
1241  mean.direction[i].blue)));
1242  if (image->colorspace == CMYKColorspace)
1243  channel_features[BlackPixelChannel].correlation[i]=
1244  (correlation.direction[i].black-mean.direction[i].black*
1245  mean.direction[i].black)/(sqrt(sum_squares.direction[i].black-
1246  (mean.direction[i].black*mean.direction[i].black))*sqrt(
1247  sum_squares.direction[i].black-(mean.direction[i].black*
1248  mean.direction[i].black)));
1249  if (image->alpha_trait != UndefinedPixelTrait)
1250  channel_features[AlphaPixelChannel].correlation[i]=
1251  (correlation.direction[i].alpha-mean.direction[i].alpha*
1252  mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
1253  (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt(
1254  sum_squares.direction[i].alpha-(mean.direction[i].alpha*
1255  mean.direction[i].alpha)));
1256  }
1257  /*
1258  Compute more texture features.
1259  */
1260 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1261  #pragma omp parallel for schedule(static) shared(status) \
1262  magick_number_threads(image,image,number_grays,1)
1263 #endif
1264  for (i=0; i < 4; i++)
1265  {
1266  ssize_t
1267  x;
1268 
1269  for (x=2; x < (ssize_t) (2*number_grays); x++)
1270  {
1271  /*
1272  Sum average.
1273  */
1274  channel_features[RedPixelChannel].sum_average[i]+=
1275  x*density_xy[x].direction[i].red;
1276  channel_features[GreenPixelChannel].sum_average[i]+=
1277  x*density_xy[x].direction[i].green;
1278  channel_features[BluePixelChannel].sum_average[i]+=
1279  x*density_xy[x].direction[i].blue;
1280  if (image->colorspace == CMYKColorspace)
1281  channel_features[BlackPixelChannel].sum_average[i]+=
1282  x*density_xy[x].direction[i].black;
1283  if (image->alpha_trait != UndefinedPixelTrait)
1284  channel_features[AlphaPixelChannel].sum_average[i]+=
1285  x*density_xy[x].direction[i].alpha;
1286  /*
1287  Sum entropy.
1288  */
1289  channel_features[RedPixelChannel].sum_entropy[i]-=
1290  density_xy[x].direction[i].red*
1291  MagickLog10(density_xy[x].direction[i].red);
1292  channel_features[GreenPixelChannel].sum_entropy[i]-=
1293  density_xy[x].direction[i].green*
1294  MagickLog10(density_xy[x].direction[i].green);
1295  channel_features[BluePixelChannel].sum_entropy[i]-=
1296  density_xy[x].direction[i].blue*
1297  MagickLog10(density_xy[x].direction[i].blue);
1298  if (image->colorspace == CMYKColorspace)
1299  channel_features[BlackPixelChannel].sum_entropy[i]-=
1300  density_xy[x].direction[i].black*
1301  MagickLog10(density_xy[x].direction[i].black);
1302  if (image->alpha_trait != UndefinedPixelTrait)
1303  channel_features[AlphaPixelChannel].sum_entropy[i]-=
1304  density_xy[x].direction[i].alpha*
1305  MagickLog10(density_xy[x].direction[i].alpha);
1306  /*
1307  Sum variance.
1308  */
1309  channel_features[RedPixelChannel].sum_variance[i]+=
1310  (x-channel_features[RedPixelChannel].sum_entropy[i])*
1311  (x-channel_features[RedPixelChannel].sum_entropy[i])*
1312  density_xy[x].direction[i].red;
1313  channel_features[GreenPixelChannel].sum_variance[i]+=
1314  (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1315  (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1316  density_xy[x].direction[i].green;
1317  channel_features[BluePixelChannel].sum_variance[i]+=
1318  (x-channel_features[BluePixelChannel].sum_entropy[i])*
1319  (x-channel_features[BluePixelChannel].sum_entropy[i])*
1320  density_xy[x].direction[i].blue;
1321  if (image->colorspace == CMYKColorspace)
1322  channel_features[BlackPixelChannel].sum_variance[i]+=
1323  (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1324  (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1325  density_xy[x].direction[i].black;
1326  if (image->alpha_trait != UndefinedPixelTrait)
1327  channel_features[AlphaPixelChannel].sum_variance[i]+=
1328  (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1329  (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1330  density_xy[x].direction[i].alpha;
1331  }
1332  }
1333  /*
1334  Compute more texture features.
1335  */
1336 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1337  #pragma omp parallel for schedule(static) shared(status) \
1338  magick_number_threads(image,image,number_grays,1)
1339 #endif
1340  for (i=0; i < 4; i++)
1341  {
1342  ssize_t
1343  y;
1344 
1345  for (y=0; y < (ssize_t) number_grays; y++)
1346  {
1347  ssize_t
1348  x;
1349 
1350  for (x=0; x < (ssize_t) number_grays; x++)
1351  {
1352  /*
1353  Sum of Squares: Variance
1354  */
1355  variance.direction[i].red+=(y-mean.direction[i].red+1)*
1356  (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
1357  variance.direction[i].green+=(y-mean.direction[i].green+1)*
1358  (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
1359  variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
1360  (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
1361  if (image->colorspace == CMYKColorspace)
1362  variance.direction[i].black+=(y-mean.direction[i].black+1)*
1363  (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black;
1364  if (image->alpha_trait != UndefinedPixelTrait)
1365  variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
1366  (y-mean.direction[i].alpha+1)*
1367  cooccurrence[x][y].direction[i].alpha;
1368  /*
1369  Sum average / Difference Variance.
1370  */
1371  density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
1372  cooccurrence[x][y].direction[i].red;
1373  density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
1374  cooccurrence[x][y].direction[i].green;
1375  density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
1376  cooccurrence[x][y].direction[i].blue;
1377  if (image->colorspace == CMYKColorspace)
1378  density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
1379  cooccurrence[x][y].direction[i].black;
1380  if (image->alpha_trait != UndefinedPixelTrait)
1381  density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
1382  cooccurrence[x][y].direction[i].alpha;
1383  /*
1384  Information Measures of Correlation.
1385  */
1386  entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
1387  MagickLog10(cooccurrence[x][y].direction[i].red);
1388  entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
1389  MagickLog10(cooccurrence[x][y].direction[i].green);
1390  entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
1391  MagickLog10(cooccurrence[x][y].direction[i].blue);
1392  if (image->colorspace == CMYKColorspace)
1393  entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black*
1394  MagickLog10(cooccurrence[x][y].direction[i].black);
1395  if (image->alpha_trait != UndefinedPixelTrait)
1396  entropy_xy.direction[i].alpha-=
1397  cooccurrence[x][y].direction[i].alpha*MagickLog10(
1398  cooccurrence[x][y].direction[i].alpha);
1399  entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
1400  MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red));
1401  entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
1402  MagickLog10(density_x[x].direction[i].green*
1403  density_y[y].direction[i].green));
1404  entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
1405  MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue));
1406  if (image->colorspace == CMYKColorspace)
1407  entropy_xy1.direction[i].black-=(
1408  cooccurrence[x][y].direction[i].black*MagickLog10(
1409  density_x[x].direction[i].black*density_y[y].direction[i].black));
1410  if (image->alpha_trait != UndefinedPixelTrait)
1411  entropy_xy1.direction[i].alpha-=(
1412  cooccurrence[x][y].direction[i].alpha*MagickLog10(
1413  density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1414  entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
1415  density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red*
1416  density_y[y].direction[i].red));
1417  entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
1418  density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green*
1419  density_y[y].direction[i].green));
1420  entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
1421  density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue*
1422  density_y[y].direction[i].blue));
1423  if (image->colorspace == CMYKColorspace)
1424  entropy_xy2.direction[i].black-=(density_x[x].direction[i].black*
1425  density_y[y].direction[i].black*MagickLog10(
1426  density_x[x].direction[i].black*density_y[y].direction[i].black));
1427  if (image->alpha_trait != UndefinedPixelTrait)
1428  entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha*
1429  density_y[y].direction[i].alpha*MagickLog10(
1430  density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1431  }
1432  }
1433  channel_features[RedPixelChannel].variance_sum_of_squares[i]=
1434  variance.direction[i].red;
1435  channel_features[GreenPixelChannel].variance_sum_of_squares[i]=
1436  variance.direction[i].green;
1437  channel_features[BluePixelChannel].variance_sum_of_squares[i]=
1438  variance.direction[i].blue;
1439  if (image->colorspace == CMYKColorspace)
1440  channel_features[BlackPixelChannel].variance_sum_of_squares[i]=
1441  variance.direction[i].black;
1442  if (image->alpha_trait != UndefinedPixelTrait)
1443  channel_features[AlphaPixelChannel].variance_sum_of_squares[i]=
1444  variance.direction[i].alpha;
1445  }
1446  /*
1447  Compute more texture features.
1448  */
1449  (void) memset(&variance,0,sizeof(variance));
1450  (void) memset(&sum_squares,0,sizeof(sum_squares));
1451 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1452  #pragma omp parallel for schedule(static) shared(status) \
1453  magick_number_threads(image,image,number_grays,1)
1454 #endif
1455  for (i=0; i < 4; i++)
1456  {
1457  ssize_t
1458  x;
1459 
1460  for (x=0; x < (ssize_t) number_grays; x++)
1461  {
1462  /*
1463  Difference variance.
1464  */
1465  variance.direction[i].red+=density_xy[x].direction[i].red;
1466  variance.direction[i].green+=density_xy[x].direction[i].green;
1467  variance.direction[i].blue+=density_xy[x].direction[i].blue;
1468  if (image->colorspace == CMYKColorspace)
1469  variance.direction[i].black+=density_xy[x].direction[i].black;
1470  if (image->alpha_trait != UndefinedPixelTrait)
1471  variance.direction[i].alpha+=density_xy[x].direction[i].alpha;
1472  sum_squares.direction[i].red+=density_xy[x].direction[i].red*
1473  density_xy[x].direction[i].red;
1474  sum_squares.direction[i].green+=density_xy[x].direction[i].green*
1475  density_xy[x].direction[i].green;
1476  sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
1477  density_xy[x].direction[i].blue;
1478  if (image->colorspace == CMYKColorspace)
1479  sum_squares.direction[i].black+=density_xy[x].direction[i].black*
1480  density_xy[x].direction[i].black;
1481  if (image->alpha_trait != UndefinedPixelTrait)
1482  sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
1483  density_xy[x].direction[i].alpha;
1484  /*
1485  Difference entropy.
1486  */
1487  channel_features[RedPixelChannel].difference_entropy[i]-=
1488  density_xy[x].direction[i].red*
1489  MagickLog10(density_xy[x].direction[i].red);
1490  channel_features[GreenPixelChannel].difference_entropy[i]-=
1491  density_xy[x].direction[i].green*
1492  MagickLog10(density_xy[x].direction[i].green);
1493  channel_features[BluePixelChannel].difference_entropy[i]-=
1494  density_xy[x].direction[i].blue*
1495  MagickLog10(density_xy[x].direction[i].blue);
1496  if (image->colorspace == CMYKColorspace)
1497  channel_features[BlackPixelChannel].difference_entropy[i]-=
1498  density_xy[x].direction[i].black*
1499  MagickLog10(density_xy[x].direction[i].black);
1500  if (image->alpha_trait != UndefinedPixelTrait)
1501  channel_features[AlphaPixelChannel].difference_entropy[i]-=
1502  density_xy[x].direction[i].alpha*
1503  MagickLog10(density_xy[x].direction[i].alpha);
1504  /*
1505  Information Measures of Correlation.
1506  */
1507  entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1508  MagickLog10(density_x[x].direction[i].red));
1509  entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1510  MagickLog10(density_x[x].direction[i].green));
1511  entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1512  MagickLog10(density_x[x].direction[i].blue));
1513  if (image->colorspace == CMYKColorspace)
1514  entropy_x.direction[i].black-=(density_x[x].direction[i].black*
1515  MagickLog10(density_x[x].direction[i].black));
1516  if (image->alpha_trait != UndefinedPixelTrait)
1517  entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha*
1518  MagickLog10(density_x[x].direction[i].alpha));
1519  entropy_y.direction[i].red-=(density_y[x].direction[i].red*
1520  MagickLog10(density_y[x].direction[i].red));
1521  entropy_y.direction[i].green-=(density_y[x].direction[i].green*
1522  MagickLog10(density_y[x].direction[i].green));
1523  entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
1524  MagickLog10(density_y[x].direction[i].blue));
1525  if (image->colorspace == CMYKColorspace)
1526  entropy_y.direction[i].black-=(density_y[x].direction[i].black*
1527  MagickLog10(density_y[x].direction[i].black));
1528  if (image->alpha_trait != UndefinedPixelTrait)
1529  entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
1530  MagickLog10(density_y[x].direction[i].alpha));
1531  }
1532  /*
1533  Difference variance.
1534  */
1535  channel_features[RedPixelChannel].difference_variance[i]=
1536  (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1537  (variance.direction[i].red*variance.direction[i].red))/
1538  ((double) number_grays*number_grays*number_grays*number_grays);
1539  channel_features[GreenPixelChannel].difference_variance[i]=
1540  (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1541  (variance.direction[i].green*variance.direction[i].green))/
1542  ((double) number_grays*number_grays*number_grays*number_grays);
1543  channel_features[BluePixelChannel].difference_variance[i]=
1544  (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1545  (variance.direction[i].blue*variance.direction[i].blue))/
1546  ((double) number_grays*number_grays*number_grays*number_grays);
1547  if (image->colorspace == CMYKColorspace)
1548  channel_features[BlackPixelChannel].difference_variance[i]=
1549  (((double) number_grays*number_grays*sum_squares.direction[i].black)-
1550  (variance.direction[i].black*variance.direction[i].black))/
1551  ((double) number_grays*number_grays*number_grays*number_grays);
1552  if (image->alpha_trait != UndefinedPixelTrait)
1553  channel_features[AlphaPixelChannel].difference_variance[i]=
1554  (((double) number_grays*number_grays*sum_squares.direction[i].alpha)-
1555  (variance.direction[i].alpha*variance.direction[i].alpha))/
1556  ((double) number_grays*number_grays*number_grays*number_grays);
1557  /*
1558  Information Measures of Correlation.
1559  */
1560  channel_features[RedPixelChannel].measure_of_correlation_1[i]=
1561  (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1562  (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1563  entropy_x.direction[i].red : entropy_y.direction[i].red);
1564  channel_features[GreenPixelChannel].measure_of_correlation_1[i]=
1565  (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1566  (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1567  entropy_x.direction[i].green : entropy_y.direction[i].green);
1568  channel_features[BluePixelChannel].measure_of_correlation_1[i]=
1569  (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1570  (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1571  entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1572  if (image->colorspace == CMYKColorspace)
1573  channel_features[BlackPixelChannel].measure_of_correlation_1[i]=
1574  (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/
1575  (entropy_x.direction[i].black > entropy_y.direction[i].black ?
1576  entropy_x.direction[i].black : entropy_y.direction[i].black);
1577  if (image->alpha_trait != UndefinedPixelTrait)
1578  channel_features[AlphaPixelChannel].measure_of_correlation_1[i]=
1579  (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/
1580  (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ?
1581  entropy_x.direction[i].alpha : entropy_y.direction[i].alpha);
1582  channel_features[RedPixelChannel].measure_of_correlation_2[i]=
1583  (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].red-
1584  entropy_xy.direction[i].red)))));
1585  channel_features[GreenPixelChannel].measure_of_correlation_2[i]=
1586  (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].green-
1587  entropy_xy.direction[i].green)))));
1588  channel_features[BluePixelChannel].measure_of_correlation_2[i]=
1589  (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].blue-
1590  entropy_xy.direction[i].blue)))));
1591  if (image->colorspace == CMYKColorspace)
1592  channel_features[BlackPixelChannel].measure_of_correlation_2[i]=
1593  (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].black-
1594  entropy_xy.direction[i].black)))));
1595  if (image->alpha_trait != UndefinedPixelTrait)
1596  channel_features[AlphaPixelChannel].measure_of_correlation_2[i]=
1597  (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].alpha-
1598  entropy_xy.direction[i].alpha)))));
1599  }
1600  /*
1601  Compute more texture features.
1602  */
1603 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1604  #pragma omp parallel for schedule(static) shared(status) \
1605  magick_number_threads(image,image,number_grays,1)
1606 #endif
1607  for (i=0; i < 4; i++)
1608  {
1609  ssize_t
1610  z;
1611 
1612  for (z=0; z < (ssize_t) number_grays; z++)
1613  {
1614  ssize_t
1615  y;
1616 
1618  pixel;
1619 
1620  (void) memset(&pixel,0,sizeof(pixel));
1621  for (y=0; y < (ssize_t) number_grays; y++)
1622  {
1623  ssize_t
1624  x;
1625 
1626  for (x=0; x < (ssize_t) number_grays; x++)
1627  {
1628  /*
1629  Contrast: amount of local variations present in an image.
1630  */
1631  if (((y-x) == z) || ((x-y) == z))
1632  {
1633  pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1634  pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1635  pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1636  if (image->colorspace == CMYKColorspace)
1637  pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
1638  if (image->alpha_trait != UndefinedPixelTrait)
1639  pixel.direction[i].alpha+=
1640  cooccurrence[x][y].direction[i].alpha;
1641  }
1642  /*
1643  Maximum Correlation Coefficient.
1644  */
1645  if ((fabs(density_x[z].direction[i].red) > MagickEpsilon) &&
1646  (fabs(density_y[x].direction[i].red) > MagickEpsilon))
1647  Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1648  cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1649  density_y[x].direction[i].red;
1650  if ((fabs(density_x[z].direction[i].green) > MagickEpsilon) &&
1651  (fabs(density_y[x].direction[i].red) > MagickEpsilon))
1652  Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1653  cooccurrence[y][x].direction[i].green/
1654  density_x[z].direction[i].green/density_y[x].direction[i].red;
1655  if ((fabs(density_x[z].direction[i].blue) > MagickEpsilon) &&
1656  (fabs(density_y[x].direction[i].blue) > MagickEpsilon))
1657  Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1658  cooccurrence[y][x].direction[i].blue/
1659  density_x[z].direction[i].blue/density_y[x].direction[i].blue;
1660  if (image->colorspace == CMYKColorspace)
1661  if ((fabs(density_x[z].direction[i].black) > MagickEpsilon) &&
1662  (fabs(density_y[x].direction[i].black) > MagickEpsilon))
1663  Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
1664  cooccurrence[y][x].direction[i].black/
1665  density_x[z].direction[i].black/density_y[x].direction[i].black;
1666  if (image->alpha_trait != UndefinedPixelTrait)
1667  if ((fabs(density_x[z].direction[i].alpha) > MagickEpsilon) &&
1668  (fabs(density_y[x].direction[i].alpha) > MagickEpsilon))
1669  Q[z][y].direction[i].alpha+=
1670  cooccurrence[z][x].direction[i].alpha*
1671  cooccurrence[y][x].direction[i].alpha/
1672  density_x[z].direction[i].alpha/
1673  density_y[x].direction[i].alpha;
1674  }
1675  }
1676  channel_features[RedPixelChannel].contrast[i]+=z*z*
1677  pixel.direction[i].red;
1678  channel_features[GreenPixelChannel].contrast[i]+=z*z*
1679  pixel.direction[i].green;
1680  channel_features[BluePixelChannel].contrast[i]+=z*z*
1681  pixel.direction[i].blue;
1682  if (image->colorspace == CMYKColorspace)
1683  channel_features[BlackPixelChannel].contrast[i]+=z*z*
1684  pixel.direction[i].black;
1685  if (image->alpha_trait != UndefinedPixelTrait)
1686  channel_features[AlphaPixelChannel].contrast[i]+=z*z*
1687  pixel.direction[i].alpha;
1688  }
1689  /*
1690  Maximum Correlation Coefficient.
1691  Future: return second largest eigenvalue of Q.
1692  */
1693  channel_features[RedPixelChannel].maximum_correlation_coefficient[i]=
1694  sqrt((double) -1.0);
1695  channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]=
1696  sqrt((double) -1.0);
1697  channel_features[BluePixelChannel].maximum_correlation_coefficient[i]=
1698  sqrt((double) -1.0);
1699  if (image->colorspace == CMYKColorspace)
1700  channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]=
1701  sqrt((double) -1.0);
1702  if (image->alpha_trait != UndefinedPixelTrait)
1703  channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]=
1704  sqrt((double) -1.0);
1705  }
1706  /*
1707  Relinquish resources.
1708  */
1709  sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1710  for (i=0; i < (ssize_t) number_grays; i++)
1711  Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1712  Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1713  density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1714  density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1715  density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1716  for (i=0; i < (ssize_t) number_grays; i++)
1717  cooccurrence[i]=(ChannelStatistics *)
1718  RelinquishMagickMemory(cooccurrence[i]);
1719  cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1720  return(channel_features);
1721 }
1722 ␌
1723 /*
1724 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1725 % %
1726 % %
1727 % %
1728 % H o u g h L i n e I m a g e %
1729 % %
1730 % %
1731 % %
1732 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1733 %
1734 % HoughLineImage() can be used in conjunction with any binary edge extracted
1735 % image (we recommend Canny) to identify lines in the image. The algorithm
1736 % accumulates counts for every white pixel for every possible orientation (for
1737 % angles from 0 to 179 in 1 degree increments) and distance from the center of
1738 % the image to the corner (in 1 px increments) and stores the counts in an
1739 % accumulator matrix of angle vs distance. The size of the accumulator is
1740 % 180x(diagonal/2). Next it searches this space for peaks in counts and
1741 % converts the locations of the peaks to slope and intercept in the normal
1742 % x,y input image space. Use the slope/intercepts to find the endpoints
1743 % clipped to the bounds of the image. The lines are then drawn. The counts
1744 % are a measure of the length of the lines.
1745 %
1746 % The format of the HoughLineImage method is:
1747 %
1748 % Image *HoughLineImage(const Image *image,const size_t width,
1749 % const size_t height,const size_t threshold,ExceptionInfo *exception)
1750 %
1751 % A description of each parameter follows:
1752 %
1753 % o image: the image.
1754 %
1755 % o width, height: find line pairs as local maxima in this neighborhood.
1756 %
1757 % o threshold: the line count threshold.
1758 %
1759 % o exception: return any errors or warnings in this structure.
1760 %
1761 */
1762 
1763 static inline double MagickRound(double x)
1764 {
1765  /*
1766  Round the fraction to nearest integer.
1767  */
1768  if ((x-floor(x)) < (ceil(x)-x))
1769  return(floor(x));
1770  return(ceil(x));
1771 }
1772 
1773 static Image *RenderHoughLines(const ImageInfo *image_info,const size_t columns,
1774  const size_t rows,ExceptionInfo *exception)
1775 {
1776 #define BoundingBox "viewbox"
1777 
1778  DrawInfo
1779  *draw_info;
1780 
1781  Image
1782  *image;
1783 
1784  MagickBooleanType
1785  status;
1786 
1787  /*
1788  Open image.
1789  */
1790  image=AcquireImage(image_info,exception);
1791  status=OpenBlob(image_info,image,ReadBinaryBlobMode,exception);
1792  if (status == MagickFalse)
1793  {
1794  image=DestroyImageList(image);
1795  return((Image *) NULL);
1796  }
1797  image->columns=columns;
1798  image->rows=rows;
1799  draw_info=CloneDrawInfo(image_info,(DrawInfo *) NULL);
1800  draw_info->affine.sx=image->resolution.x == 0.0 ? 1.0 : image->resolution.x/
1801  DefaultResolution;
1802  draw_info->affine.sy=image->resolution.y == 0.0 ? 1.0 : image->resolution.y/
1803  DefaultResolution;
1804  image->columns=(size_t) (draw_info->affine.sx*image->columns);
1805  image->rows=(size_t) (draw_info->affine.sy*image->rows);
1806  status=SetImageExtent(image,image->columns,image->rows,exception);
1807  if (status == MagickFalse)
1808  return(DestroyImageList(image));
1809  if (SetImageBackgroundColor(image,exception) == MagickFalse)
1810  {
1811  image=DestroyImageList(image);
1812  return((Image *) NULL);
1813  }
1814  /*
1815  Render drawing.
1816  */
1817  if (GetBlobStreamData(image) == (unsigned char *) NULL)
1818  draw_info->primitive=FileToString(image->filename,~0UL,exception);
1819  else
1820  {
1821  draw_info->primitive=(char *) AcquireQuantumMemory(1,(size_t)
1822  GetBlobSize(image)+1);
1823  if (draw_info->primitive != (char *) NULL)
1824  {
1825  (void) memcpy(draw_info->primitive,GetBlobStreamData(image),
1826  (size_t) GetBlobSize(image));
1827  draw_info->primitive[GetBlobSize(image)]='\0';
1828  }
1829  }
1830  (void) DrawImage(image,draw_info,exception);
1831  draw_info=DestroyDrawInfo(draw_info);
1832  (void) CloseBlob(image);
1833  return(GetFirstImageInList(image));
1834 }
1835 
1836 MagickExport Image *HoughLineImage(const Image *image,const size_t width,
1837  const size_t height,const size_t threshold,ExceptionInfo *exception)
1838 {
1839 #define HoughLineImageTag "HoughLine/Image"
1840 
1841  CacheView
1842  *image_view;
1843 
1844  char
1845  message[MagickPathExtent],
1846  path[MagickPathExtent];
1847 
1848  const char
1849  *artifact;
1850 
1851  double
1852  hough_height;
1853 
1854  Image
1855  *lines_image = NULL;
1856 
1857  ImageInfo
1858  *image_info;
1859 
1860  int
1861  file;
1862 
1863  MagickBooleanType
1864  status;
1865 
1866  MagickOffsetType
1867  progress;
1868 
1869  MatrixInfo
1870  *accumulator;
1871 
1872  PointInfo
1873  center;
1874 
1875  ssize_t
1876  y;
1877 
1878  size_t
1879  accumulator_height,
1880  accumulator_width,
1881  line_count;
1882 
1883  /*
1884  Create the accumulator.
1885  */
1886  assert(image != (const Image *) NULL);
1887  assert(image->signature == MagickCoreSignature);
1888  assert(exception != (ExceptionInfo *) NULL);
1889  assert(exception->signature == MagickCoreSignature);
1890  if (IsEventLogging() != MagickFalse)
1891  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1892  accumulator_width=180;
1893  hough_height=((sqrt(2.0)*(double) (image->rows > image->columns ?
1894  image->rows : image->columns))/2.0);
1895  accumulator_height=(size_t) (2.0*hough_height);
1896  accumulator=AcquireMatrixInfo(accumulator_width,accumulator_height,
1897  sizeof(double),exception);
1898  if (accumulator == (MatrixInfo *) NULL)
1899  ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1900  if (NullMatrix(accumulator) == MagickFalse)
1901  {
1902  accumulator=DestroyMatrixInfo(accumulator);
1903  ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1904  }
1905  /*
1906  Populate the accumulator.
1907  */
1908  status=MagickTrue;
1909  progress=0;
1910  center.x=(double) image->columns/2.0;
1911  center.y=(double) image->rows/2.0;
1912  image_view=AcquireVirtualCacheView(image,exception);
1913  for (y=0; y < (ssize_t) image->rows; y++)
1914  {
1915  const Quantum
1916  *magick_restrict p;
1917 
1918  ssize_t
1919  x;
1920 
1921  if (status == MagickFalse)
1922  continue;
1923  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1924  if (p == (Quantum *) NULL)
1925  {
1926  status=MagickFalse;
1927  continue;
1928  }
1929  for (x=0; x < (ssize_t) image->columns; x++)
1930  {
1931  if (GetPixelIntensity(image,p) > (QuantumRange/2.0))
1932  {
1933  ssize_t
1934  i;
1935 
1936  for (i=0; i < 180; i++)
1937  {
1938  double
1939  count,
1940  radius;
1941 
1942  radius=(((double) x-center.x)*cos(DegreesToRadians((double) i)))+
1943  (((double) y-center.y)*sin(DegreesToRadians((double) i)));
1944  (void) GetMatrixElement(accumulator,i,(ssize_t)
1945  MagickRound(radius+hough_height),&count);
1946  count++;
1947  (void) SetMatrixElement(accumulator,i,(ssize_t)
1948  MagickRound(radius+hough_height),&count);
1949  }
1950  }
1951  p+=GetPixelChannels(image);
1952  }
1953  if (image->progress_monitor != (MagickProgressMonitor) NULL)
1954  {
1955  MagickBooleanType
1956  proceed;
1957 
1958 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1959  #pragma omp atomic
1960 #endif
1961  progress++;
1962  proceed=SetImageProgress(image,CannyEdgeImageTag,progress,image->rows);
1963  if (proceed == MagickFalse)
1964  status=MagickFalse;
1965  }
1966  }
1967  image_view=DestroyCacheView(image_view);
1968  if (status == MagickFalse)
1969  {
1970  accumulator=DestroyMatrixInfo(accumulator);
1971  return((Image *) NULL);
1972  }
1973  /*
1974  Generate line segments from accumulator.
1975  */
1976  file=AcquireUniqueFileResource(path);
1977  if (file == -1)
1978  {
1979  accumulator=DestroyMatrixInfo(accumulator);
1980  return((Image *) NULL);
1981  }
1982  (void) FormatLocaleString(message,MagickPathExtent,
1983  "# Hough line transform: %.20gx%.20g%+.20g\n",(double) width,
1984  (double) height,(double) threshold);
1985  if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1986  status=MagickFalse;
1987  (void) FormatLocaleString(message,MagickPathExtent,
1988  "viewbox 0 0 %.20g %.20g\n",(double) image->columns,(double) image->rows);
1989  if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1990  status=MagickFalse;
1991  (void) FormatLocaleString(message,MagickPathExtent,
1992  "# x1,y1 x2,y2 # count angle distance\n");
1993  if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1994  status=MagickFalse;
1995  line_count=image->columns > image->rows ? image->columns/4 : image->rows/4;
1996  if (threshold != 0)
1997  line_count=threshold;
1998  for (y=0; y < (ssize_t) accumulator_height; y++)
1999  {
2000  ssize_t
2001  x;
2002 
2003  for (x=0; x < (ssize_t) accumulator_width; x++)
2004  {
2005  double
2006  count;
2007 
2008  (void) GetMatrixElement(accumulator,x,y,&count);
2009  if (count >= (double) line_count)
2010  {
2011  double
2012  maxima;
2013 
2014  SegmentInfo
2015  line;
2016 
2017  ssize_t
2018  v;
2019 
2020  /*
2021  Is point a local maxima?
2022  */
2023  maxima=count;
2024  for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
2025  {
2026  ssize_t
2027  u;
2028 
2029  for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
2030  {
2031  if ((u != 0) || (v !=0))
2032  {
2033  (void) GetMatrixElement(accumulator,x+u,y+v,&count);
2034  if (count > maxima)
2035  {
2036  maxima=count;
2037  break;
2038  }
2039  }
2040  }
2041  if (u < (ssize_t) (width/2))
2042  break;
2043  }
2044  (void) GetMatrixElement(accumulator,x,y,&count);
2045  if (maxima > count)
2046  continue;
2047  if ((x >= 45) && (x <= 135))
2048  {
2049  /*
2050  y = (r-x cos(t))/sin(t)
2051  */
2052  line.x1=0.0;
2053  line.y1=((double) (y-(accumulator_height/2.0))-((line.x1-
2054  (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
2055  sin(DegreesToRadians((double) x))+(image->rows/2.0);
2056  line.x2=(double) image->columns;
2057  line.y2=((double) (y-(accumulator_height/2.0))-((line.x2-
2058  (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
2059  sin(DegreesToRadians((double) x))+(image->rows/2.0);
2060  }
2061  else
2062  {
2063  /*
2064  x = (r-y cos(t))/sin(t)
2065  */
2066  line.y1=0.0;
2067  line.x1=((double) (y-(accumulator_height/2.0))-((line.y1-
2068  (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
2069  cos(DegreesToRadians((double) x))+(image->columns/2.0);
2070  line.y2=(double) image->rows;
2071  line.x2=((double) (y-(accumulator_height/2.0))-((line.y2-
2072  (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
2073  cos(DegreesToRadians((double) x))+(image->columns/2.0);
2074  }
2075  (void) FormatLocaleString(message,MagickPathExtent,
2076  "line %g,%g %g,%g # %g %g %g\n",line.x1,line.y1,line.x2,line.y2,
2077  maxima,(double) x,(double) y);
2078  if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
2079  status=MagickFalse;
2080  }
2081  }
2082  }
2083  (void) close(file);
2084  /*
2085  Render lines to image canvas.
2086  */
2087  image_info=AcquireImageInfo();
2088  image_info->background_color=image->background_color;
2089  (void) FormatLocaleString(image_info->filename,MagickPathExtent,"%s",path);
2090  artifact=GetImageArtifact(image,"background");
2091  if (artifact != (const char *) NULL)
2092  (void) SetImageOption(image_info,"background",artifact);
2093  artifact=GetImageArtifact(image,"fill");
2094  if (artifact != (const char *) NULL)
2095  (void) SetImageOption(image_info,"fill",artifact);
2096  artifact=GetImageArtifact(image,"stroke");
2097  if (artifact != (const char *) NULL)
2098  (void) SetImageOption(image_info,"stroke",artifact);
2099  artifact=GetImageArtifact(image,"strokewidth");
2100  if (artifact != (const char *) NULL)
2101  (void) SetImageOption(image_info,"strokewidth",artifact);
2102  lines_image=RenderHoughLines(image_info,image->columns,image->rows,exception);
2103  artifact=GetImageArtifact(image,"hough-lines:accumulator");
2104  if ((lines_image != (Image *) NULL) &&
2105  (IsStringTrue(artifact) != MagickFalse))
2106  {
2107  Image
2108  *accumulator_image;
2109 
2110  accumulator_image=MatrixToImage(accumulator,exception);
2111  if (accumulator_image != (Image *) NULL)
2112  AppendImageToList(&lines_image,accumulator_image);
2113  }
2114  /*
2115  Free resources.
2116  */
2117  accumulator=DestroyMatrixInfo(accumulator);
2118  image_info=DestroyImageInfo(image_info);
2119  (void) RelinquishUniqueFileResource(path);
2120  return(GetFirstImageInList(lines_image));
2121 }
2122 ␌
2123 /*
2124 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2125 % %
2126 % %
2127 % %
2128 % M e a n S h i f t I m a g e %
2129 % %
2130 % %
2131 % %
2132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2133 %
2134 % MeanShiftImage() delineate arbitrarily shaped clusters in the image. For
2135 % each pixel, it visits all the pixels in the neighborhood specified by
2136 % the window centered at the pixel and excludes those that are outside the
2137 % radius=(window-1)/2 surrounding the pixel. From those pixels, it finds those
2138 % that are within the specified color distance from the current mean, and
2139 % computes a new x,y centroid from those coordinates and a new mean. This new
2140 % x,y centroid is used as the center for a new window. This process iterates
2141 % until it converges and the final mean is replaces the (original window
2142 % center) pixel value. It repeats this process for the next pixel, etc.,
2143 % until it processes all pixels in the image. Results are typically better with
2144 % colorspaces other than sRGB. We recommend YIQ, YUV or YCbCr.
2145 %
2146 % The format of the MeanShiftImage method is:
2147 %
2148 % Image *MeanShiftImage(const Image *image,const size_t width,
2149 % const size_t height,const double color_distance,
2150 % ExceptionInfo *exception)
2151 %
2152 % A description of each parameter follows:
2153 %
2154 % o image: the image.
2155 %
2156 % o width, height: find pixels in this neighborhood.
2157 %
2158 % o color_distance: the color distance.
2159 %
2160 % o exception: return any errors or warnings in this structure.
2161 %
2162 */
2163 MagickExport Image *MeanShiftImage(const Image *image,const size_t width,
2164  const size_t height,const double color_distance,ExceptionInfo *exception)
2165 {
2166 #define MaxMeanShiftIterations 100
2167 #define MeanShiftImageTag "MeanShift/Image"
2168 
2169  CacheView
2170  *image_view,
2171  *mean_view,
2172  *pixel_view;
2173 
2174  Image
2175  *mean_image;
2176 
2177  MagickBooleanType
2178  status;
2179 
2180  MagickOffsetType
2181  progress;
2182 
2183  ssize_t
2184  y;
2185 
2186  assert(image != (const Image *) NULL);
2187  assert(image->signature == MagickCoreSignature);
2188  assert(exception != (ExceptionInfo *) NULL);
2189  assert(exception->signature == MagickCoreSignature);
2190  if (IsEventLogging() != MagickFalse)
2191  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2192  mean_image=CloneImage(image,0,0,MagickTrue,exception);
2193  if (mean_image == (Image *) NULL)
2194  return((Image *) NULL);
2195  if (SetImageStorageClass(mean_image,DirectClass,exception) == MagickFalse)
2196  {
2197  mean_image=DestroyImage(mean_image);
2198  return((Image *) NULL);
2199  }
2200  status=MagickTrue;
2201  progress=0;
2202  image_view=AcquireVirtualCacheView(image,exception);
2203  pixel_view=AcquireVirtualCacheView(image,exception);
2204  mean_view=AcquireAuthenticCacheView(mean_image,exception);
2205 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2206  #pragma omp parallel for schedule(static) shared(status,progress) \
2207  magick_number_threads(mean_image,mean_image,mean_image->rows,1)
2208 #endif
2209  for (y=0; y < (ssize_t) mean_image->rows; y++)
2210  {
2211  const Quantum
2212  *magick_restrict p;
2213 
2214  Quantum
2215  *magick_restrict q;
2216 
2217  ssize_t
2218  x;
2219 
2220  if (status == MagickFalse)
2221  continue;
2222  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2223  q=GetCacheViewAuthenticPixels(mean_view,0,y,mean_image->columns,1,
2224  exception);
2225  if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2226  {
2227  status=MagickFalse;
2228  continue;
2229  }
2230  for (x=0; x < (ssize_t) mean_image->columns; x++)
2231  {
2232  PixelInfo
2233  mean_pixel,
2234  previous_pixel;
2235 
2236  PointInfo
2237  mean_location,
2238  previous_location;
2239 
2240  ssize_t
2241  i;
2242 
2243  GetPixelInfo(image,&mean_pixel);
2244  GetPixelInfoPixel(image,p,&mean_pixel);
2245  mean_location.x=(double) x;
2246  mean_location.y=(double) y;
2247  for (i=0; i < MaxMeanShiftIterations; i++)
2248  {
2249  double
2250  distance,
2251  gamma;
2252 
2253  PixelInfo
2254  sum_pixel;
2255 
2256  PointInfo
2257  sum_location;
2258 
2259  ssize_t
2260  count,
2261  v;
2262 
2263  sum_location.x=0.0;
2264  sum_location.y=0.0;
2265  GetPixelInfo(image,&sum_pixel);
2266  previous_location=mean_location;
2267  previous_pixel=mean_pixel;
2268  count=0;
2269  for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
2270  {
2271  ssize_t
2272  u;
2273 
2274  for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
2275  {
2276  if ((v*v+u*u) <= (ssize_t) ((width/2)*(height/2)))
2277  {
2278  PixelInfo
2279  pixel;
2280 
2281  status=GetOneCacheViewVirtualPixelInfo(pixel_view,(ssize_t)
2282  MagickRound(mean_location.x+u),(ssize_t) MagickRound(
2283  mean_location.y+v),&pixel,exception);
2284  distance=(mean_pixel.red-pixel.red)*(mean_pixel.red-pixel.red)+
2285  (mean_pixel.green-pixel.green)*(mean_pixel.green-pixel.green)+
2286  (mean_pixel.blue-pixel.blue)*(mean_pixel.blue-pixel.blue);
2287  if (distance <= (color_distance*color_distance))
2288  {
2289  sum_location.x+=mean_location.x+u;
2290  sum_location.y+=mean_location.y+v;
2291  sum_pixel.red+=pixel.red;
2292  sum_pixel.green+=pixel.green;
2293  sum_pixel.blue+=pixel.blue;
2294  sum_pixel.alpha+=pixel.alpha;
2295  count++;
2296  }
2297  }
2298  }
2299  }
2300  gamma=PerceptibleReciprocal(count);
2301  mean_location.x=gamma*sum_location.x;
2302  mean_location.y=gamma*sum_location.y;
2303  mean_pixel.red=gamma*sum_pixel.red;
2304  mean_pixel.green=gamma*sum_pixel.green;
2305  mean_pixel.blue=gamma*sum_pixel.blue;
2306  mean_pixel.alpha=gamma*sum_pixel.alpha;
2307  distance=(mean_location.x-previous_location.x)*
2308  (mean_location.x-previous_location.x)+
2309  (mean_location.y-previous_location.y)*
2310  (mean_location.y-previous_location.y)+
2311  255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)*
2312  255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)+
2313  255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)*
2314  255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)+
2315  255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue)*
2316  255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue);
2317  if (distance <= 3.0)
2318  break;
2319  }
2320  SetPixelRed(mean_image,ClampToQuantum(mean_pixel.red),q);
2321  SetPixelGreen(mean_image,ClampToQuantum(mean_pixel.green),q);
2322  SetPixelBlue(mean_image,ClampToQuantum(mean_pixel.blue),q);
2323  SetPixelAlpha(mean_image,ClampToQuantum(mean_pixel.alpha),q);
2324  p+=GetPixelChannels(image);
2325  q+=GetPixelChannels(mean_image);
2326  }
2327  if (SyncCacheViewAuthenticPixels(mean_view,exception) == MagickFalse)
2328  status=MagickFalse;
2329  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2330  {
2331  MagickBooleanType
2332  proceed;
2333 
2334 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2335  #pragma omp atomic
2336 #endif
2337  progress++;
2338  proceed=SetImageProgress(image,MeanShiftImageTag,progress,image->rows);
2339  if (proceed == MagickFalse)
2340  status=MagickFalse;
2341  }
2342  }
2343  mean_view=DestroyCacheView(mean_view);
2344  pixel_view=DestroyCacheView(pixel_view);
2345  image_view=DestroyCacheView(image_view);
2346  return(mean_image);
2347 }
Definition: image.h:152