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