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