MagickCore  6.9.13-47
Convert, Edit, Or Compose Bitmap Images
feature.c
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % FFFFF EEEEE AAA TTTTT U U RRRR EEEEE %
7 % F E A A T U U R R E %
8 % FFF EEE AAAAA T U U RRRR EEE %
9 % F E A A T U U R R E %
10 % F EEEEE A A T UUU R R EEEEE %
11 % %
12 % %
13 % MagickCore Image Feature Methods %
14 % %
15 % Software Design %
16 % Cristy %
17 % July 1992 %
18 % %
19 % %
20 % Copyright 1999 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
22 % %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
25 % %
26 % https://imagemagick.org/license/ %
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 
170  PixelPacket
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  {
437  PixelPacket
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  {
626  DoublePixelPacket
627  direction[4]; /* horizontal, vertical, left and right diagonals */
628  } ChannelStatistics;
629 
630  CacheView
631  *image_view;
632 
633  ChannelFeatures
634  *channel_features;
635 
636  ChannelStatistics
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 
653  LongPixelPacket
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  Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
805  if (sum != (ChannelStatistics *) NULL)
806  sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
807  if (density_y != (ChannelStatistics *) NULL)
808  density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
809  if (density_xy != (ChannelStatistics *) NULL)
810  density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
811  if (density_x != (ChannelStatistics *) NULL)
812  density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
813  if (cooccurrence != (ChannelStatistics **) NULL)
814  cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
815  cooccurrence);
816  grays=(LongPixelPacket *) RelinquishMagickMemory(grays);
817  channel_features=(ChannelFeatures *) RelinquishMagickMemory(
818  channel_features);
819  (void) ThrowMagickException(exception,GetMagickModule(),
820  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
821  return(channel_features);
822  }
823  (void) memset(&correlation,0,sizeof(correlation));
824  (void) memset(density_x,0,2*(number_grays+1)*sizeof(*density_x));
825  (void) memset(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
826  (void) memset(density_y,0,2*(number_grays+1)*sizeof(*density_y));
827  (void) memset(&mean,0,sizeof(mean));
828  (void) memset(sum,0,number_grays*sizeof(*sum));
829  (void) memset(&sum_squares,0,sizeof(sum_squares));
830  (void) memset(density_xy,0,2*number_grays*sizeof(*density_xy));
831  (void) memset(&entropy_x,0,sizeof(entropy_x));
832  (void) memset(&entropy_xy,0,sizeof(entropy_xy));
833  (void) memset(&entropy_xy1,0,sizeof(entropy_xy1));
834  (void) memset(&entropy_xy2,0,sizeof(entropy_xy2));
835  (void) memset(&entropy_y,0,sizeof(entropy_y));
836  (void) memset(&variance,0,sizeof(variance));
837  for (i=0; i < (ssize_t) number_grays; i++)
838  {
839  cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
840  sizeof(**cooccurrence));
841  Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
842  if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
843  (Q[i] == (ChannelStatistics *) NULL))
844  break;
845  (void) memset(cooccurrence[i],0,number_grays*
846  sizeof(**cooccurrence));
847  (void) memset(Q[i],0,number_grays*sizeof(**Q));
848  }
849  if (i < (ssize_t) number_grays)
850  {
851  for (i--; i >= 0; i--)
852  {
853  if (Q[i] != (ChannelStatistics *) NULL)
854  Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
855  if (cooccurrence[i] != (ChannelStatistics *) NULL)
856  cooccurrence[i]=(ChannelStatistics *)
857  RelinquishMagickMemory(cooccurrence[i]);
858  }
859  Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
860  cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
861  sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
862  density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
863  density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
864  density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
865  grays=(LongPixelPacket *) RelinquishMagickMemory(grays);
866  channel_features=(ChannelFeatures *) RelinquishMagickMemory(
867  channel_features);
868  (void) ThrowMagickException(exception,GetMagickModule(),
869  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
870  return(channel_features);
871  }
872  /*
873  Initialize spatial dependence matrix.
874  */
875  status=MagickTrue;
876  image_view=AcquireVirtualCacheView(image,exception);
877  for (y=0; y < (ssize_t) image->rows; y++)
878  {
879  const IndexPacket
880  *magick_restrict indexes;
881 
882  const PixelPacket
883  *magick_restrict p;
884 
885  ssize_t
886  x;
887 
888  ssize_t
889  i,
890  offset,
891  u,
892  v;
893 
894  if (status == MagickFalse)
895  continue;
896  p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,y,image->columns+
897  2*distance,distance+2,exception);
898  if (p == (const PixelPacket *) NULL)
899  {
900  status=MagickFalse;
901  continue;
902  }
903  indexes=GetCacheViewVirtualIndexQueue(image_view);
904  p+=(ptrdiff_t) distance;
905  indexes+=distance;
906  for (x=0; x < (ssize_t) image->columns; x++)
907  {
908  for (i=0; i < 4; i++)
909  {
910  switch (i)
911  {
912  case 0:
913  default:
914  {
915  /*
916  Horizontal adjacency.
917  */
918  offset=(ssize_t) distance;
919  break;
920  }
921  case 1:
922  {
923  /*
924  Vertical adjacency.
925  */
926  offset=(ssize_t) (image->columns+2*distance);
927  break;
928  }
929  case 2:
930  {
931  /*
932  Right diagonal adjacency.
933  */
934  offset=(ssize_t) ((image->columns+2*distance)-distance);
935  break;
936  }
937  case 3:
938  {
939  /*
940  Left diagonal adjacency.
941  */
942  offset=(ssize_t) ((image->columns+2*distance)+distance);
943  break;
944  }
945  }
946  u=0;
947  v=0;
948  while (grays[u].red != ScaleQuantumToMap(GetPixelRed(p)))
949  u++;
950  while (grays[v].red != ScaleQuantumToMap(GetPixelRed(p+offset)))
951  v++;
952  cooccurrence[u][v].direction[i].red++;
953  cooccurrence[v][u].direction[i].red++;
954  u=0;
955  v=0;
956  while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(p)))
957  u++;
958  while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(p+offset)))
959  v++;
960  cooccurrence[u][v].direction[i].green++;
961  cooccurrence[v][u].direction[i].green++;
962  u=0;
963  v=0;
964  while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(p)))
965  u++;
966  while (grays[v].blue != ScaleQuantumToMap((p+offset)->blue))
967  v++;
968  cooccurrence[u][v].direction[i].blue++;
969  cooccurrence[v][u].direction[i].blue++;
970  if (image->colorspace == CMYKColorspace)
971  {
972  u=0;
973  v=0;
974  while (grays[u].index != ScaleQuantumToMap(GetPixelIndex(indexes+x)))
975  u++;
976  while (grays[v].index != ScaleQuantumToMap(GetPixelIndex(indexes+x+offset)))
977  v++;
978  cooccurrence[u][v].direction[i].index++;
979  cooccurrence[v][u].direction[i].index++;
980  }
981  if (image->matte != MagickFalse)
982  {
983  u=0;
984  v=0;
985  while (grays[u].opacity != ScaleQuantumToMap(GetPixelOpacity(p)))
986  u++;
987  while (grays[v].opacity != ScaleQuantumToMap((p+offset)->opacity))
988  v++;
989  cooccurrence[u][v].direction[i].opacity++;
990  cooccurrence[v][u].direction[i].opacity++;
991  }
992  }
993  p++;
994  }
995  }
996  grays=(LongPixelPacket *) RelinquishMagickMemory(grays);
997  image_view=DestroyCacheView(image_view);
998  if (status == MagickFalse)
999  {
1000  for (i=0; i < (ssize_t) number_grays; i++)
1001  cooccurrence[i]=(ChannelStatistics *)
1002  RelinquishMagickMemory(cooccurrence[i]);
1003  cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1004  channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1005  channel_features);
1006  (void) ThrowMagickException(exception,GetMagickModule(),
1007  ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1008  return(channel_features);
1009  }
1010  /*
1011  Normalize spatial dependence matrix.
1012  */
1013  for (i=0; i < 4; i++)
1014  {
1015  double
1016  normalize;
1017 
1018  ssize_t
1019  y;
1020 
1021  switch (i)
1022  {
1023  case 0:
1024  default:
1025  {
1026  /*
1027  Horizontal adjacency.
1028  */
1029  normalize=2.0*image->rows*(image->columns-distance);
1030  break;
1031  }
1032  case 1:
1033  {
1034  /*
1035  Vertical adjacency.
1036  */
1037  normalize=2.0*(image->rows-distance)*image->columns;
1038  break;
1039  }
1040  case 2:
1041  {
1042  /*
1043  Right diagonal adjacency.
1044  */
1045  normalize=2.0*(image->rows-distance)*(image->columns-distance);
1046  break;
1047  }
1048  case 3:
1049  {
1050  /*
1051  Left diagonal adjacency.
1052  */
1053  normalize=2.0*(image->rows-distance)*(image->columns-distance);
1054  break;
1055  }
1056  }
1057  normalize=MagickSafeReciprocal(normalize);
1058  for (y=0; y < (ssize_t) number_grays; y++)
1059  {
1060  ssize_t
1061  x;
1062 
1063  for (x=0; x < (ssize_t) number_grays; x++)
1064  {
1065  cooccurrence[x][y].direction[i].red*=normalize;
1066  cooccurrence[x][y].direction[i].green*=normalize;
1067  cooccurrence[x][y].direction[i].blue*=normalize;
1068  if (image->colorspace == CMYKColorspace)
1069  cooccurrence[x][y].direction[i].index*=normalize;
1070  if (image->matte != MagickFalse)
1071  cooccurrence[x][y].direction[i].opacity*=normalize;
1072  }
1073  }
1074  }
1075  /*
1076  Compute texture features.
1077  */
1078 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1079  #pragma omp parallel for schedule(static) shared(status) \
1080  magick_number_threads(image,image,number_grays,1)
1081 #endif
1082  for (i=0; i < 4; i++)
1083  {
1084  ssize_t
1085  y;
1086 
1087  for (y=0; y < (ssize_t) number_grays; y++)
1088  {
1089  ssize_t
1090  x;
1091 
1092  for (x=0; x < (ssize_t) number_grays; x++)
1093  {
1094  /*
1095  Angular second moment: measure of homogeneity of the image.
1096  */
1097  channel_features[RedChannel].angular_second_moment[i]+=
1098  cooccurrence[x][y].direction[i].red*
1099  cooccurrence[x][y].direction[i].red;
1100  channel_features[GreenChannel].angular_second_moment[i]+=
1101  cooccurrence[x][y].direction[i].green*
1102  cooccurrence[x][y].direction[i].green;
1103  channel_features[BlueChannel].angular_second_moment[i]+=
1104  cooccurrence[x][y].direction[i].blue*
1105  cooccurrence[x][y].direction[i].blue;
1106  if (image->colorspace == CMYKColorspace)
1107  channel_features[BlackChannel].angular_second_moment[i]+=
1108  cooccurrence[x][y].direction[i].index*
1109  cooccurrence[x][y].direction[i].index;
1110  if (image->matte != MagickFalse)
1111  channel_features[OpacityChannel].angular_second_moment[i]+=
1112  cooccurrence[x][y].direction[i].opacity*
1113  cooccurrence[x][y].direction[i].opacity;
1114  /*
1115  Correlation: measure of linear-dependencies in the image.
1116  */
1117  sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1118  sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1119  sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1120  if (image->colorspace == CMYKColorspace)
1121  sum[y].direction[i].index+=cooccurrence[x][y].direction[i].index;
1122  if (image->matte != MagickFalse)
1123  sum[y].direction[i].opacity+=cooccurrence[x][y].direction[i].opacity;
1124  correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
1125  correlation.direction[i].green+=x*y*
1126  cooccurrence[x][y].direction[i].green;
1127  correlation.direction[i].blue+=x*y*
1128  cooccurrence[x][y].direction[i].blue;
1129  if (image->colorspace == CMYKColorspace)
1130  correlation.direction[i].index+=x*y*
1131  cooccurrence[x][y].direction[i].index;
1132  if (image->matte != MagickFalse)
1133  correlation.direction[i].opacity+=x*y*
1134  cooccurrence[x][y].direction[i].opacity;
1135  /*
1136  Inverse Difference Moment.
1137  */
1138  channel_features[RedChannel].inverse_difference_moment[i]+=
1139  cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
1140  channel_features[GreenChannel].inverse_difference_moment[i]+=
1141  cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
1142  channel_features[BlueChannel].inverse_difference_moment[i]+=
1143  cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
1144  if (image->colorspace == CMYKColorspace)
1145  channel_features[IndexChannel].inverse_difference_moment[i]+=
1146  cooccurrence[x][y].direction[i].index/((y-x)*(y-x)+1);
1147  if (image->matte != MagickFalse)
1148  channel_features[OpacityChannel].inverse_difference_moment[i]+=
1149  cooccurrence[x][y].direction[i].opacity/((y-x)*(y-x)+1);
1150  /*
1151  Sum average.
1152  */
1153  density_xy[y+x+2].direction[i].red+=
1154  cooccurrence[x][y].direction[i].red;
1155  density_xy[y+x+2].direction[i].green+=
1156  cooccurrence[x][y].direction[i].green;
1157  density_xy[y+x+2].direction[i].blue+=
1158  cooccurrence[x][y].direction[i].blue;
1159  if (image->colorspace == CMYKColorspace)
1160  density_xy[y+x+2].direction[i].index+=
1161  cooccurrence[x][y].direction[i].index;
1162  if (image->matte != MagickFalse)
1163  density_xy[y+x+2].direction[i].opacity+=
1164  cooccurrence[x][y].direction[i].opacity;
1165  /*
1166  Entropy.
1167  */
1168  channel_features[RedChannel].entropy[i]-=
1169  cooccurrence[x][y].direction[i].red*
1170  log2(cooccurrence[x][y].direction[i].red);
1171  channel_features[GreenChannel].entropy[i]-=
1172  cooccurrence[x][y].direction[i].green*
1173  log2(cooccurrence[x][y].direction[i].green);
1174  channel_features[BlueChannel].entropy[i]-=
1175  cooccurrence[x][y].direction[i].blue*
1176  log2(cooccurrence[x][y].direction[i].blue);
1177  if (image->colorspace == CMYKColorspace)
1178  channel_features[IndexChannel].entropy[i]-=
1179  cooccurrence[x][y].direction[i].index*
1180  log2(cooccurrence[x][y].direction[i].index);
1181  if (image->matte != MagickFalse)
1182  channel_features[OpacityChannel].entropy[i]-=
1183  cooccurrence[x][y].direction[i].opacity*
1184  log2(cooccurrence[x][y].direction[i].opacity);
1185  /*
1186  Information Measures of Correlation.
1187  */
1188  density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
1189  density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
1190  density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1191  if (image->colorspace == CMYKColorspace)
1192  density_x[x].direction[i].index+=
1193  cooccurrence[x][y].direction[i].index;
1194  if (image->matte != MagickFalse)
1195  density_x[x].direction[i].opacity+=
1196  cooccurrence[x][y].direction[i].opacity;
1197  density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1198  density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1199  density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1200  if (image->colorspace == CMYKColorspace)
1201  density_y[y].direction[i].index+=
1202  cooccurrence[x][y].direction[i].index;
1203  if (image->matte != MagickFalse)
1204  density_y[y].direction[i].opacity+=
1205  cooccurrence[x][y].direction[i].opacity;
1206  }
1207  mean.direction[i].red+=y*sum[y].direction[i].red;
1208  sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
1209  mean.direction[i].green+=y*sum[y].direction[i].green;
1210  sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
1211  mean.direction[i].blue+=y*sum[y].direction[i].blue;
1212  sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
1213  if (image->colorspace == CMYKColorspace)
1214  {
1215  mean.direction[i].index+=y*sum[y].direction[i].index;
1216  sum_squares.direction[i].index+=y*y*sum[y].direction[i].index;
1217  }
1218  if (image->matte != MagickFalse)
1219  {
1220  mean.direction[i].opacity+=y*sum[y].direction[i].opacity;
1221  sum_squares.direction[i].opacity+=y*y*sum[y].direction[i].opacity;
1222  }
1223  }
1224  /*
1225  Correlation: measure of linear-dependencies in the image.
1226  */
1227  channel_features[RedChannel].correlation[i]=
1228  (correlation.direction[i].red-mean.direction[i].red*
1229  mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
1230  (mean.direction[i].red*mean.direction[i].red))*sqrt(
1231  sum_squares.direction[i].red-(mean.direction[i].red*
1232  mean.direction[i].red)));
1233  channel_features[GreenChannel].correlation[i]=
1234  (correlation.direction[i].green-mean.direction[i].green*
1235  mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
1236  (mean.direction[i].green*mean.direction[i].green))*sqrt(
1237  sum_squares.direction[i].green-(mean.direction[i].green*
1238  mean.direction[i].green)));
1239  channel_features[BlueChannel].correlation[i]=
1240  (correlation.direction[i].blue-mean.direction[i].blue*
1241  mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
1242  (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
1243  sum_squares.direction[i].blue-(mean.direction[i].blue*
1244  mean.direction[i].blue)));
1245  if (image->colorspace == CMYKColorspace)
1246  channel_features[IndexChannel].correlation[i]=
1247  (correlation.direction[i].index-mean.direction[i].index*
1248  mean.direction[i].index)/(sqrt(sum_squares.direction[i].index-
1249  (mean.direction[i].index*mean.direction[i].index))*sqrt(
1250  sum_squares.direction[i].index-(mean.direction[i].index*
1251  mean.direction[i].index)));
1252  if (image->matte != MagickFalse)
1253  channel_features[OpacityChannel].correlation[i]=
1254  (correlation.direction[i].opacity-mean.direction[i].opacity*
1255  mean.direction[i].opacity)/(sqrt(sum_squares.direction[i].opacity-
1256  (mean.direction[i].opacity*mean.direction[i].opacity))*sqrt(
1257  sum_squares.direction[i].opacity-(mean.direction[i].opacity*
1258  mean.direction[i].opacity)));
1259  }
1260  /*
1261  Compute more texture features.
1262  */
1263 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1264  #pragma omp parallel for schedule(static) shared(status) \
1265  magick_number_threads(image,image,number_grays,1)
1266 #endif
1267  for (i=0; i < 4; i++)
1268  {
1269  ssize_t
1270  x;
1271 
1272  for (x=2; x < (ssize_t) (2*number_grays); x++)
1273  {
1274  /*
1275  Sum average.
1276  */
1277  channel_features[RedChannel].sum_average[i]+=
1278  x*density_xy[x].direction[i].red;
1279  channel_features[GreenChannel].sum_average[i]+=
1280  x*density_xy[x].direction[i].green;
1281  channel_features[BlueChannel].sum_average[i]+=
1282  x*density_xy[x].direction[i].blue;
1283  if (image->colorspace == CMYKColorspace)
1284  channel_features[IndexChannel].sum_average[i]+=
1285  x*density_xy[x].direction[i].index;
1286  if (image->matte != MagickFalse)
1287  channel_features[OpacityChannel].sum_average[i]+=
1288  x*density_xy[x].direction[i].opacity;
1289  /*
1290  Sum entropy.
1291  */
1292  channel_features[RedChannel].sum_entropy[i]-=
1293  density_xy[x].direction[i].red*
1294  log2(density_xy[x].direction[i].red);
1295  channel_features[GreenChannel].sum_entropy[i]-=
1296  density_xy[x].direction[i].green*
1297  log2(density_xy[x].direction[i].green);
1298  channel_features[BlueChannel].sum_entropy[i]-=
1299  density_xy[x].direction[i].blue*
1300  log2(density_xy[x].direction[i].blue);
1301  if (image->colorspace == CMYKColorspace)
1302  channel_features[IndexChannel].sum_entropy[i]-=
1303  density_xy[x].direction[i].index*
1304  log2(density_xy[x].direction[i].index);
1305  if (image->matte != MagickFalse)
1306  channel_features[OpacityChannel].sum_entropy[i]-=
1307  density_xy[x].direction[i].opacity*
1308  log2(density_xy[x].direction[i].opacity);
1309  /*
1310  Sum variance.
1311  */
1312  channel_features[RedChannel].sum_variance[i]+=
1313  (x-channel_features[RedChannel].sum_entropy[i])*
1314  (x-channel_features[RedChannel].sum_entropy[i])*
1315  density_xy[x].direction[i].red;
1316  channel_features[GreenChannel].sum_variance[i]+=
1317  (x-channel_features[GreenChannel].sum_entropy[i])*
1318  (x-channel_features[GreenChannel].sum_entropy[i])*
1319  density_xy[x].direction[i].green;
1320  channel_features[BlueChannel].sum_variance[i]+=
1321  (x-channel_features[BlueChannel].sum_entropy[i])*
1322  (x-channel_features[BlueChannel].sum_entropy[i])*
1323  density_xy[x].direction[i].blue;
1324  if (image->colorspace == CMYKColorspace)
1325  channel_features[IndexChannel].sum_variance[i]+=
1326  (x-channel_features[IndexChannel].sum_entropy[i])*
1327  (x-channel_features[IndexChannel].sum_entropy[i])*
1328  density_xy[x].direction[i].index;
1329  if (image->matte != MagickFalse)
1330  channel_features[OpacityChannel].sum_variance[i]+=
1331  (x-channel_features[OpacityChannel].sum_entropy[i])*
1332  (x-channel_features[OpacityChannel].sum_entropy[i])*
1333  density_xy[x].direction[i].opacity;
1334  }
1335  }
1336  /*
1337  Compute more texture features.
1338  */
1339 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1340  #pragma omp parallel for schedule(static) shared(status) \
1341  magick_number_threads(image,image,number_grays,1)
1342 #endif
1343  for (i=0; i < 4; i++)
1344  {
1345  ssize_t
1346  y;
1347 
1348  for (y=0; y < (ssize_t) number_grays; y++)
1349  {
1350  ssize_t
1351  x;
1352 
1353  for (x=0; x < (ssize_t) number_grays; x++)
1354  {
1355  /*
1356  Sum of Squares: Variance
1357  */
1358  variance.direction[i].red+=(y-mean.direction[i].red+1)*
1359  (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
1360  variance.direction[i].green+=(y-mean.direction[i].green+1)*
1361  (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
1362  variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
1363  (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
1364  if (image->colorspace == CMYKColorspace)
1365  variance.direction[i].index+=(y-mean.direction[i].index+1)*
1366  (y-mean.direction[i].index+1)*cooccurrence[x][y].direction[i].index;
1367  if (image->matte != MagickFalse)
1368  variance.direction[i].opacity+=(y-mean.direction[i].opacity+1)*
1369  (y-mean.direction[i].opacity+1)*
1370  cooccurrence[x][y].direction[i].opacity;
1371  /*
1372  Sum average / Difference Variance.
1373  */
1374  density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
1375  cooccurrence[x][y].direction[i].red;
1376  density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
1377  cooccurrence[x][y].direction[i].green;
1378  density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
1379  cooccurrence[x][y].direction[i].blue;
1380  if (image->colorspace == CMYKColorspace)
1381  density_xy[MagickAbsoluteValue(y-x)].direction[i].index+=
1382  cooccurrence[x][y].direction[i].index;
1383  if (image->matte != MagickFalse)
1384  density_xy[MagickAbsoluteValue(y-x)].direction[i].opacity+=
1385  cooccurrence[x][y].direction[i].opacity;
1386  /*
1387  Information Measures of Correlation.
1388  */
1389  entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
1390  log2(cooccurrence[x][y].direction[i].red);
1391  entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
1392  log2(cooccurrence[x][y].direction[i].green);
1393  entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
1394  log2(cooccurrence[x][y].direction[i].blue);
1395  if (image->colorspace == CMYKColorspace)
1396  entropy_xy.direction[i].index-=cooccurrence[x][y].direction[i].index*
1397  log2(cooccurrence[x][y].direction[i].index);
1398  if (image->matte != MagickFalse)
1399  entropy_xy.direction[i].opacity-=
1400  cooccurrence[x][y].direction[i].opacity*log2(
1401  cooccurrence[x][y].direction[i].opacity);
1402  entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
1403  log2(density_x[x].direction[i].red*
1404  density_y[y].direction[i].red));
1405  entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
1406  log2(density_x[x].direction[i].green*
1407  density_y[y].direction[i].green));
1408  entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
1409  log2(density_x[x].direction[i].blue*
1410  density_y[y].direction[i].blue));
1411  if (image->colorspace == CMYKColorspace)
1412  entropy_xy1.direction[i].index-=(
1413  cooccurrence[x][y].direction[i].index*log2(
1414  density_x[x].direction[i].index*density_y[y].direction[i].index));
1415  if (image->matte != MagickFalse)
1416  entropy_xy1.direction[i].opacity-=(
1417  cooccurrence[x][y].direction[i].opacity*log2(
1418  density_x[x].direction[i].opacity*
1419  density_y[y].direction[i].opacity));
1420  entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
1421  density_y[y].direction[i].red*log2(
1422  density_x[x].direction[i].red*density_y[y].direction[i].red));
1423  entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
1424  density_y[y].direction[i].green*log2(
1425  density_x[x].direction[i].green*density_y[y].direction[i].green));
1426  entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
1427  density_y[y].direction[i].blue*log2(
1428  density_x[x].direction[i].blue*density_y[y].direction[i].blue));
1429  if (image->colorspace == CMYKColorspace)
1430  entropy_xy2.direction[i].index-=(density_x[x].direction[i].index*
1431  density_y[y].direction[i].index*log2(
1432  density_x[x].direction[i].index*density_y[y].direction[i].index));
1433  if (image->matte != MagickFalse)
1434  entropy_xy2.direction[i].opacity-=(density_x[x].direction[i].opacity*
1435  density_y[y].direction[i].opacity*log2(
1436  density_x[x].direction[i].opacity*
1437  density_y[y].direction[i].opacity));
1438  }
1439  }
1440  channel_features[RedChannel].variance_sum_of_squares[i]=
1441  variance.direction[i].red;
1442  channel_features[GreenChannel].variance_sum_of_squares[i]=
1443  variance.direction[i].green;
1444  channel_features[BlueChannel].variance_sum_of_squares[i]=
1445  variance.direction[i].blue;
1446  if (image->colorspace == CMYKColorspace)
1447  channel_features[RedChannel].variance_sum_of_squares[i]=
1448  variance.direction[i].index;
1449  if (image->matte != MagickFalse)
1450  channel_features[RedChannel].variance_sum_of_squares[i]=
1451  variance.direction[i].opacity;
1452  }
1453  /*
1454  Compute more texture features.
1455  */
1456  (void) memset(&variance,0,sizeof(variance));
1457  (void) memset(&sum_squares,0,sizeof(sum_squares));
1458 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1459  #pragma omp parallel for schedule(static) shared(status) \
1460  magick_number_threads(image,image,number_grays,1)
1461 #endif
1462  for (i=0; i < 4; i++)
1463  {
1464  ssize_t
1465  x;
1466 
1467  for (x=0; x < (ssize_t) number_grays; x++)
1468  {
1469  /*
1470  Difference variance.
1471  */
1472  variance.direction[i].red+=density_xy[x].direction[i].red;
1473  variance.direction[i].green+=density_xy[x].direction[i].green;
1474  variance.direction[i].blue+=density_xy[x].direction[i].blue;
1475  if (image->colorspace == CMYKColorspace)
1476  variance.direction[i].index+=density_xy[x].direction[i].index;
1477  if (image->matte != MagickFalse)
1478  variance.direction[i].opacity+=density_xy[x].direction[i].opacity;
1479  sum_squares.direction[i].red+=density_xy[x].direction[i].red*
1480  density_xy[x].direction[i].red;
1481  sum_squares.direction[i].green+=density_xy[x].direction[i].green*
1482  density_xy[x].direction[i].green;
1483  sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
1484  density_xy[x].direction[i].blue;
1485  if (image->colorspace == CMYKColorspace)
1486  sum_squares.direction[i].index+=density_xy[x].direction[i].index*
1487  density_xy[x].direction[i].index;
1488  if (image->matte != MagickFalse)
1489  sum_squares.direction[i].opacity+=density_xy[x].direction[i].opacity*
1490  density_xy[x].direction[i].opacity;
1491  /*
1492  Difference entropy.
1493  */
1494  channel_features[RedChannel].difference_entropy[i]-=
1495  density_xy[x].direction[i].red*
1496  log2(density_xy[x].direction[i].red);
1497  channel_features[GreenChannel].difference_entropy[i]-=
1498  density_xy[x].direction[i].green*
1499  log2(density_xy[x].direction[i].green);
1500  channel_features[BlueChannel].difference_entropy[i]-=
1501  density_xy[x].direction[i].blue*
1502  log2(density_xy[x].direction[i].blue);
1503  if (image->colorspace == CMYKColorspace)
1504  channel_features[IndexChannel].difference_entropy[i]-=
1505  density_xy[x].direction[i].index*
1506  log2(density_xy[x].direction[i].index);
1507  if (image->matte != MagickFalse)
1508  channel_features[OpacityChannel].difference_entropy[i]-=
1509  density_xy[x].direction[i].opacity*
1510  log2(density_xy[x].direction[i].opacity);
1511  /*
1512  Information Measures of Correlation.
1513  */
1514  entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1515  log2(density_x[x].direction[i].red));
1516  entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1517  log2(density_x[x].direction[i].green));
1518  entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1519  log2(density_x[x].direction[i].blue));
1520  if (image->colorspace == CMYKColorspace)
1521  entropy_x.direction[i].index-=(density_x[x].direction[i].index*
1522  log2(density_x[x].direction[i].index));
1523  if (image->matte != MagickFalse)
1524  entropy_x.direction[i].opacity-=(density_x[x].direction[i].opacity*
1525  log2(density_x[x].direction[i].opacity));
1526  entropy_y.direction[i].red-=(density_y[x].direction[i].red*
1527  log2(density_y[x].direction[i].red));
1528  entropy_y.direction[i].green-=(density_y[x].direction[i].green*
1529  log2(density_y[x].direction[i].green));
1530  entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
1531  log2(density_y[x].direction[i].blue));
1532  if (image->colorspace == CMYKColorspace)
1533  entropy_y.direction[i].index-=(density_y[x].direction[i].index*
1534  log2(density_y[x].direction[i].index));
1535  if (image->matte != MagickFalse)
1536  entropy_y.direction[i].opacity-=(density_y[x].direction[i].opacity*
1537  log2(density_y[x].direction[i].opacity));
1538  }
1539  /*
1540  Difference variance.
1541  */
1542  channel_features[RedChannel].difference_variance[i]=
1543  (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1544  (variance.direction[i].red*variance.direction[i].red))/
1545  ((double) number_grays*number_grays*number_grays*number_grays);
1546  channel_features[GreenChannel].difference_variance[i]=
1547  (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1548  (variance.direction[i].green*variance.direction[i].green))/
1549  ((double) number_grays*number_grays*number_grays*number_grays);
1550  channel_features[BlueChannel].difference_variance[i]=
1551  (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1552  (variance.direction[i].blue*variance.direction[i].blue))/
1553  ((double) number_grays*number_grays*number_grays*number_grays);
1554  if (image->matte != MagickFalse)
1555  channel_features[OpacityChannel].difference_variance[i]=
1556  (((double) number_grays*number_grays*sum_squares.direction[i].opacity)-
1557  (variance.direction[i].opacity*variance.direction[i].opacity))/
1558  ((double) number_grays*number_grays*number_grays*number_grays);
1559  if (image->colorspace == CMYKColorspace)
1560  channel_features[IndexChannel].difference_variance[i]=
1561  (((double) number_grays*number_grays*sum_squares.direction[i].index)-
1562  (variance.direction[i].index*variance.direction[i].index))/
1563  ((double) number_grays*number_grays*number_grays*number_grays);
1564  /*
1565  Information Measures of Correlation.
1566  */
1567  channel_features[RedChannel].measure_of_correlation_1[i]=
1568  (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1569  (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1570  entropy_x.direction[i].red : entropy_y.direction[i].red);
1571  channel_features[GreenChannel].measure_of_correlation_1[i]=
1572  (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1573  (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1574  entropy_x.direction[i].green : entropy_y.direction[i].green);
1575  channel_features[BlueChannel].measure_of_correlation_1[i]=
1576  (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1577  (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1578  entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1579  if (image->colorspace == CMYKColorspace)
1580  channel_features[IndexChannel].measure_of_correlation_1[i]=
1581  (entropy_xy.direction[i].index-entropy_xy1.direction[i].index)/
1582  (entropy_x.direction[i].index > entropy_y.direction[i].index ?
1583  entropy_x.direction[i].index : entropy_y.direction[i].index);
1584  if (image->matte != MagickFalse)
1585  channel_features[OpacityChannel].measure_of_correlation_1[i]=
1586  (entropy_xy.direction[i].opacity-entropy_xy1.direction[i].opacity)/
1587  (entropy_x.direction[i].opacity > entropy_y.direction[i].opacity ?
1588  entropy_x.direction[i].opacity : entropy_y.direction[i].opacity);
1589  channel_features[RedChannel].measure_of_correlation_2[i]=
1590  (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].red-
1591  entropy_xy.direction[i].red)))));
1592  channel_features[GreenChannel].measure_of_correlation_2[i]=
1593  (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].green-
1594  entropy_xy.direction[i].green)))));
1595  channel_features[BlueChannel].measure_of_correlation_2[i]=
1596  (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].blue-
1597  entropy_xy.direction[i].blue)))));
1598  if (image->colorspace == CMYKColorspace)
1599  channel_features[IndexChannel].measure_of_correlation_2[i]=
1600  (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].index-
1601  entropy_xy.direction[i].index)))));
1602  if (image->matte != MagickFalse)
1603  channel_features[OpacityChannel].measure_of_correlation_2[i]=
1604  (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].opacity-
1605  entropy_xy.direction[i].opacity)))));
1606  }
1607  /*
1608  Compute more texture features.
1609  */
1610 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1611  #pragma omp parallel for schedule(static) shared(status) \
1612  magick_number_threads(image,image,number_grays,1)
1613 #endif
1614  for (i=0; i < 4; i++)
1615  {
1616  ssize_t
1617  z;
1618 
1619  for (z=0; z < (ssize_t) number_grays; z++)
1620  {
1621  ssize_t
1622  y;
1623 
1624  ChannelStatistics
1625  pixel;
1626 
1627  (void) memset(&pixel,0,sizeof(pixel));
1628  for (y=0; y < (ssize_t) number_grays; y++)
1629  {
1630  ssize_t
1631  x;
1632 
1633  for (x=0; x < (ssize_t) number_grays; x++)
1634  {
1635  /*
1636  Contrast: amount of local variations present in an image.
1637  */
1638  if (((y-x) == z) || ((x-y) == z))
1639  {
1640  pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1641  pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1642  pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1643  if (image->colorspace == CMYKColorspace)
1644  pixel.direction[i].index+=cooccurrence[x][y].direction[i].index;
1645  if (image->matte != MagickFalse)
1646  pixel.direction[i].opacity+=
1647  cooccurrence[x][y].direction[i].opacity;
1648  }
1649  /*
1650  Maximum Correlation Coefficient.
1651  */
1652  if ((fabs(density_x[z].direction[i].red) > MagickEpsilon) &&
1653  (fabs(density_y[x].direction[i].red) > MagickEpsilon))
1654  Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1655  cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1656  density_y[x].direction[i].red;
1657  if ((fabs(density_x[z].direction[i].green) > MagickEpsilon) &&
1658  (fabs(density_y[x].direction[i].red) > MagickEpsilon))
1659  Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1660  cooccurrence[y][x].direction[i].green/
1661  density_x[z].direction[i].green/density_y[x].direction[i].red;
1662  if ((fabs(density_x[z].direction[i].blue) > MagickEpsilon) &&
1663  (fabs(density_y[x].direction[i].blue) > MagickEpsilon))
1664  Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1665  cooccurrence[y][x].direction[i].blue/
1666  density_x[z].direction[i].blue/density_y[x].direction[i].blue;
1667  if (image->colorspace == CMYKColorspace)
1668  if ((fabs(density_x[z].direction[i].index) > MagickEpsilon) &&
1669  (fabs(density_y[x].direction[i].index) > MagickEpsilon))
1670  Q[z][y].direction[i].index+=cooccurrence[z][x].direction[i].index*
1671  cooccurrence[y][x].direction[i].index/
1672  density_x[z].direction[i].index/density_y[x].direction[i].index;
1673  if (image->matte != MagickFalse)
1674  if ((fabs(density_x[z].direction[i].opacity) > MagickEpsilon) &&
1675  (fabs(density_y[x].direction[i].opacity) > MagickEpsilon))
1676  Q[z][y].direction[i].opacity+=
1677  cooccurrence[z][x].direction[i].opacity*
1678  cooccurrence[y][x].direction[i].opacity/
1679  density_x[z].direction[i].opacity/
1680  density_y[x].direction[i].opacity;
1681  }
1682  }
1683  channel_features[RedChannel].contrast[i]+=z*z*pixel.direction[i].red;
1684  channel_features[GreenChannel].contrast[i]+=z*z*pixel.direction[i].green;
1685  channel_features[BlueChannel].contrast[i]+=z*z*pixel.direction[i].blue;
1686  if (image->colorspace == CMYKColorspace)
1687  channel_features[BlackChannel].contrast[i]+=z*z*
1688  pixel.direction[i].index;
1689  if (image->matte != MagickFalse)
1690  channel_features[OpacityChannel].contrast[i]+=z*z*
1691  pixel.direction[i].opacity;
1692  }
1693  /*
1694  Maximum Correlation Coefficient.
1695  Future: return second largest eigenvalue of Q.
1696  */
1697  channel_features[RedChannel].maximum_correlation_coefficient[i]=
1698  sqrt(-1.0);
1699  channel_features[GreenChannel].maximum_correlation_coefficient[i]=
1700  sqrt(-1.0);
1701  channel_features[BlueChannel].maximum_correlation_coefficient[i]=
1702  sqrt(-1.0);
1703  if (image->colorspace == CMYKColorspace)
1704  channel_features[IndexChannel].maximum_correlation_coefficient[i]=
1705  sqrt(-1.0);
1706  if (image->matte != MagickFalse)
1707  channel_features[OpacityChannel].maximum_correlation_coefficient[i]=
1708  sqrt(-1.0);
1709  }
1710  /*
1711  Relinquish resources.
1712  */
1713  sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1714  for (i=0; i < (ssize_t) number_grays; i++)
1715  Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1716  Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1717  density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1718  density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1719  density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1720  for (i=0; i < (ssize_t) number_grays; i++)
1721  cooccurrence[i]=(ChannelStatistics *)
1722  RelinquishMagickMemory(cooccurrence[i]);
1723  cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1724  return(channel_features);
1725 }
1726 ␌
1727 /*
1728 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1729 % %
1730 % %
1731 % %
1732 % H o u g h L i n e I m a g e %
1733 % %
1734 % %
1735 % %
1736 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1737 %
1738 % Use HoughLineImage() in conjunction with any binary edge extracted image (we
1739 % recommand Canny) to identify lines in the image. The algorithm accumulates
1740 % counts for every white pixel for every possible orientation (for angles from
1741 % 0 to 179 in 1 degree increments) and distance from the center of the image to
1742 % the corner (in 1 px increments) and stores the counts in an accumulator
1743 % 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
1744 % of the peaks to slope and intercept in the normal x,y input image space. Use
1745 % the slope/intercepts to find the endpoints clipped to the bounds of the
1746 % image. The lines are then drawn. The counts are a measure of the length of
1747 % the lines.
1748 %
1749 % The format of the HoughLineImage method is:
1750 %
1751 % Image *HoughLineImage(const Image *image,const size_t width,
1752 % const size_t height,const size_t threshold,ExceptionInfo *exception)
1753 %
1754 % A description of each parameter follows:
1755 %
1756 % o image: the image.
1757 %
1758 % o width, height: find line pairs as local maxima in this neighborhood.
1759 %
1760 % o threshold: the line count threshold.
1761 %
1762 % o exception: return any errors or warnings in this structure.
1763 %
1764 */
1765 
1766 static inline double MagickRound(double x)
1767 {
1768  /*
1769  Round the fraction to nearest integer.
1770  */
1771  if ((x-floor(x)) < (ceil(x)-x))
1772  return(floor(x));
1773  return(ceil(x));
1774 }
1775 
1776 static Image *RenderHoughLines(const ImageInfo *image_info,const size_t columns,
1777  const size_t rows,ExceptionInfo *exception)
1778 {
1779 #define BoundingBox "viewbox"
1780 
1781  DrawInfo
1782  *draw_info;
1783 
1784  Image
1785  *image;
1786 
1787  MagickBooleanType
1788  status;
1789 
1790  /*
1791  Open image.
1792  */
1793  image=AcquireImage(image_info);
1794  status=OpenBlob(image_info,image,ReadBinaryBlobMode,exception);
1795  if (status == MagickFalse)
1796  {
1797  image=DestroyImageList(image);
1798  return((Image *) NULL);
1799  }
1800  image->columns=columns;
1801  image->rows=rows;
1802  draw_info=CloneDrawInfo(image_info,(DrawInfo *) NULL);
1803  draw_info->affine.sx=image->x_resolution == 0.0 ? 1.0 : image->x_resolution/
1804  DefaultResolution;
1805  draw_info->affine.sy=image->y_resolution == 0.0 ? 1.0 : image->y_resolution/
1806  DefaultResolution;
1807  image->columns=(size_t) (draw_info->affine.sx*image->columns);
1808  image->rows=(size_t) (draw_info->affine.sy*image->rows);
1809  status=SetImageExtent(image,image->columns,image->rows);
1810  if (status == MagickFalse)
1811  return(DestroyImageList(image));
1812  if (SetImageBackgroundColor(image) == MagickFalse)
1813  {
1814  image=DestroyImageList(image);
1815  return((Image *) NULL);
1816  }
1817  /*
1818  Render drawing.
1819  */
1820  if (GetBlobStreamData(image) == (unsigned char *) NULL)
1821  draw_info->primitive=FileToString(image->filename,~0UL,exception);
1822  else
1823  {
1824  draw_info->primitive=(char *) AcquireQuantumMemory(1,(size_t)
1825  GetBlobSize(image)+1);
1826  if (draw_info->primitive != (char *) NULL)
1827  {
1828  (void) memcpy(draw_info->primitive,GetBlobStreamData(image),
1829  (size_t) GetBlobSize(image));
1830  draw_info->primitive[GetBlobSize(image)]='\0';
1831  }
1832  }
1833  (void) DrawImage(image,draw_info);
1834  draw_info=DestroyDrawInfo(draw_info);
1835  if (CloseBlob(image) == MagickFalse)
1836  image=DestroyImageList(image);
1837  return(GetFirstImageInList(image));
1838 }
1839 
1840 MagickExport Image *HoughLineImage(const Image *image,const size_t width,
1841  const size_t height,const size_t threshold,ExceptionInfo *exception)
1842 {
1843 #define HoughLineImageTag "HoughLine/Image"
1844 
1845  CacheView
1846  *image_view;
1847 
1848  char
1849  message[MaxTextExtent],
1850  path[MaxTextExtent];
1851 
1852  const char
1853  *artifact;
1854 
1855  double
1856  hough_height;
1857 
1858  Image
1859  *lines_image = NULL;
1860 
1861  ImageInfo
1862  *image_info;
1863 
1864  int
1865  file;
1866 
1867  MagickBooleanType
1868  status;
1869 
1870  MagickOffsetType
1871  progress;
1872 
1873  MatrixInfo
1874  *accumulator;
1875 
1876  PointInfo
1877  center;
1878 
1879  ssize_t
1880  y;
1881 
1882  size_t
1883  accumulator_height,
1884  accumulator_width,
1885  line_count;
1886 
1887  /*
1888  Create the accumulator.
1889  */
1890  assert(image != (const Image *) NULL);
1891  assert(image->signature == MagickCoreSignature);
1892  assert(exception != (ExceptionInfo *) NULL);
1893  assert(exception->signature == MagickCoreSignature);
1894  if (IsEventLogging() != MagickFalse)
1895  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1896  accumulator_width=180;
1897  hough_height=((sqrt(2.0)*(double) (image->rows > image->columns ?
1898  image->rows : image->columns))/2.0);
1899  accumulator_height=(size_t) (2.0*hough_height);
1900  accumulator=AcquireMatrixInfo(accumulator_width,accumulator_height,
1901  sizeof(double),exception);
1902  if (accumulator == (MatrixInfo *) NULL)
1903  ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1904  if (NullMatrix(accumulator) == MagickFalse)
1905  {
1906  accumulator=DestroyMatrixInfo(accumulator);
1907  ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1908  }
1909  /*
1910  Populate the accumulator.
1911  */
1912  status=MagickTrue;
1913  progress=0;
1914  center.x=(double) image->columns/2.0;
1915  center.y=(double) image->rows/2.0;
1916  image_view=AcquireVirtualCacheView(image,exception);
1917  for (y=0; y < (ssize_t) image->rows; y++)
1918  {
1919  const PixelPacket
1920  *magick_restrict p;
1921 
1922  ssize_t
1923  x;
1924 
1925  if (status == MagickFalse)
1926  continue;
1927  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1928  if (p == (PixelPacket *) NULL)
1929  {
1930  status=MagickFalse;
1931  continue;
1932  }
1933  for (x=0; x < (ssize_t) image->columns; x++)
1934  {
1935  if (GetPixelIntensity(image,p) > ((MagickRealType) QuantumRange/2.0))
1936  {
1937  ssize_t
1938  i;
1939 
1940  for (i=0; i < 180; i++)
1941  {
1942  double
1943  count,
1944  radius;
1945 
1946  radius=(((double) x-center.x)*cos(DegreesToRadians((double) i)))+
1947  (((double) y-center.y)*sin(DegreesToRadians((double) i)));
1948  (void) GetMatrixElement(accumulator,i,(ssize_t)
1949  MagickRound(radius+hough_height),&count);
1950  count++;
1951  (void) SetMatrixElement(accumulator,i,(ssize_t)
1952  MagickRound(radius+hough_height),&count);
1953  }
1954  }
1955  p++;
1956  }
1957  if (image->progress_monitor != (MagickProgressMonitor) NULL)
1958  {
1959  MagickBooleanType
1960  proceed;
1961 
1962 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1963  #pragma omp atomic
1964 #endif
1965  progress++;
1966  proceed=SetImageProgress(image,HoughLineImageTag,progress,image->rows);
1967  if (proceed == MagickFalse)
1968  status=MagickFalse;
1969  }
1970  }
1971  image_view=DestroyCacheView(image_view);
1972  if (status == MagickFalse)
1973  {
1974  accumulator=DestroyMatrixInfo(accumulator);
1975  return((Image *) NULL);
1976  }
1977  /*
1978  Generate line segments from accumulator.
1979  */
1980  file=AcquireUniqueFileResource(path);
1981  if (file == -1)
1982  {
1983  accumulator=DestroyMatrixInfo(accumulator);
1984  return((Image *) NULL);
1985  }
1986  (void) FormatLocaleString(message,MaxTextExtent,
1987  "# Hough line transform: %.20gx%.20g%+.20g\n",(double) width,
1988  (double) height,(double) threshold);
1989  if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1990  status=MagickFalse;
1991  (void) FormatLocaleString(message,MaxTextExtent,"viewbox 0 0 %.20g %.20g\n",
1992  (double) image->columns,(double) image->rows);
1993  if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1994  status=MagickFalse;
1995  (void) FormatLocaleString(message,MaxTextExtent,
1996  "# x1,y1 x2,y2 # count angle distance\n");
1997  if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1998  status=MagickFalse;
1999  line_count=image->columns > image->rows ? image->columns/4 : image->rows/4;
2000  if (threshold != 0)
2001  line_count=threshold;
2002  for (y=0; y < (ssize_t) accumulator_height; y++)
2003  {
2004  ssize_t
2005  x;
2006 
2007  for (x=0; x < (ssize_t) accumulator_width; x++)
2008  {
2009  double
2010  count;
2011 
2012  (void) GetMatrixElement(accumulator,x,y,&count);
2013  if (count >= (double) line_count)
2014  {
2015  double
2016  maxima;
2017 
2018  SegmentInfo
2019  line;
2020 
2021  ssize_t
2022  v;
2023 
2024  /*
2025  Is point a local maxima?
2026  */
2027  maxima=count;
2028  for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
2029  {
2030  ssize_t
2031  u;
2032 
2033  for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
2034  {
2035  if ((u != 0) || (v !=0))
2036  {
2037  (void) GetMatrixElement(accumulator,x+u,y+v,&count);
2038  if (count > maxima)
2039  {
2040  maxima=count;
2041  break;
2042  }
2043  }
2044  }
2045  if (u < (ssize_t) (width/2))
2046  break;
2047  }
2048  (void) GetMatrixElement(accumulator,x,y,&count);
2049  if (maxima > count)
2050  continue;
2051  if ((x >= 45) && (x <= 135))
2052  {
2053  /*
2054  y = (r-x cos(t))/sin(t)
2055  */
2056  line.x1=0.0;
2057  line.y1=((double) (y-(accumulator_height/2.0))-((line.x1-
2058  (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
2059  sin(DegreesToRadians((double) x))+(image->rows/2.0);
2060  line.x2=(double) image->columns;
2061  line.y2=((double) (y-(accumulator_height/2.0))-((line.x2-
2062  (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
2063  sin(DegreesToRadians((double) x))+(image->rows/2.0);
2064  }
2065  else
2066  {
2067  /*
2068  x = (r-y cos(t))/sin(t)
2069  */
2070  line.y1=0.0;
2071  line.x1=((double) (y-(accumulator_height/2.0))-((line.y1-
2072  (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
2073  cos(DegreesToRadians((double) x))+(image->columns/2.0);
2074  line.y2=(double) image->rows;
2075  line.x2=((double) (y-(accumulator_height/2.0))-((line.y2-
2076  (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
2077  cos(DegreesToRadians((double) x))+(image->columns/2.0);
2078  }
2079  (void) FormatLocaleString(message,MaxTextExtent,
2080  "line %g,%g %g,%g # %g %g %g\n",line.x1,line.y1,line.x2,line.y2,
2081  maxima,(double) x,(double) y);
2082  if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
2083  status=MagickFalse;
2084  }
2085  }
2086  }
2087  (void) close(file);
2088  /*
2089  Render lines to image canvas.
2090  */
2091  image_info=AcquireImageInfo();
2092  image_info->background_color=image->background_color;
2093  (void) FormatLocaleString(image_info->filename,MaxTextExtent,"%s",path);
2094  artifact=GetImageArtifact(image,"background");
2095  if (artifact != (const char *) NULL)
2096  (void) SetImageOption(image_info,"background",artifact);
2097  artifact=GetImageArtifact(image,"fill");
2098  if (artifact != (const char *) NULL)
2099  (void) SetImageOption(image_info,"fill",artifact);
2100  artifact=GetImageArtifact(image,"stroke");
2101  if (artifact != (const char *) NULL)
2102  (void) SetImageOption(image_info,"stroke",artifact);
2103  artifact=GetImageArtifact(image,"strokewidth");
2104  if (artifact != (const char *) NULL)
2105  (void) SetImageOption(image_info,"strokewidth",artifact);
2106  lines_image=RenderHoughLines(image_info,image->columns,image->rows,exception);
2107  artifact=GetImageArtifact(image,"hough-lines:accumulator");
2108  if ((lines_image != (Image *) NULL) &&
2109  (IsMagickTrue(artifact) != MagickFalse))
2110  {
2111  Image
2112  *accumulator_image;
2113 
2114  accumulator_image=MatrixToImage(accumulator,exception);
2115  if (accumulator_image != (Image *) NULL)
2116  AppendImageToList(&lines_image,accumulator_image);
2117  }
2118  /*
2119  Free resources.
2120  */
2121  accumulator=DestroyMatrixInfo(accumulator);
2122  image_info=DestroyImageInfo(image_info);
2123  (void) RelinquishUniqueFileResource(path);
2124  return(GetFirstImageInList(lines_image));
2125 }
2126 ␌
2127 /*
2128 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2129 % %
2130 % %
2131 % %
2132 % M e a n S h i f t I m a g e %
2133 % %
2134 % %
2135 % %
2136 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2137 %
2138 % MeanShiftImage() delineate arbitrarily shaped clusters in the image. For
2139 % each pixel, it visits all the pixels in the neighborhood specified by
2140 % the window centered at the pixel and excludes those that are outside the
2141 % radius=(window-1)/2 surrounding the pixel. From those pixels, it finds those
2142 % that are within the specified color distance from the current mean, and
2143 % computes a new x,y centroid from those coordinates and a new mean. This new
2144 % x,y centroid is used as the center for a new window. This process iterates
2145 % until it converges and the final mean is replaces the (original window
2146 % center) pixel value. It repeats this process for the next pixel, etc.,
2147 % until it processes all pixels in the image. Results are typically better with
2148 % colorspaces other than sRGB. We recommend YIQ, YUV or YCbCr.
2149 %
2150 % The format of the MeanShiftImage method is:
2151 %
2152 % Image *MeanShiftImage(const Image *image,const size_t width,
2153 % const size_t height,const double color_distance,
2154 % ExceptionInfo *exception)
2155 %
2156 % A description of each parameter follows:
2157 %
2158 % o image: the image.
2159 %
2160 % o width, height: find pixels in this neighborhood.
2161 %
2162 % o color_distance: the color distance.
2163 %
2164 % o exception: return any errors or warnings in this structure.
2165 %
2166 */
2167 MagickExport Image *MeanShiftImage(const Image *image,const size_t width,
2168  const size_t height,const double color_distance,ExceptionInfo *exception)
2169 {
2170 #define MaxMeanShiftIterations 100
2171 #define MeanShiftImageTag "MeanShift/Image"
2172 
2173  CacheView
2174  *image_view,
2175  *mean_view,
2176  *pixel_view;
2177 
2178  Image
2179  *mean_image;
2180 
2181  MagickBooleanType
2182  status;
2183 
2184  MagickOffsetType
2185  progress;
2186 
2187  ssize_t
2188  y;
2189 
2190  assert(image != (const Image *) NULL);
2191  assert(image->signature == MagickCoreSignature);
2192  assert(exception != (ExceptionInfo *) NULL);
2193  assert(exception->signature == MagickCoreSignature);
2194  if (IsEventLogging() != MagickFalse)
2195  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2196  mean_image=CloneImage(image,0,0,MagickTrue,exception);
2197  if (mean_image == (Image *) NULL)
2198  return((Image *) NULL);
2199  if (SetImageStorageClass(mean_image,DirectClass) == MagickFalse)
2200  {
2201  InheritException(exception,&mean_image->exception);
2202  mean_image=DestroyImage(mean_image);
2203  return((Image *) NULL);
2204  }
2205  status=MagickTrue;
2206  progress=0;
2207  image_view=AcquireVirtualCacheView(image,exception);
2208  pixel_view=AcquireVirtualCacheView(image,exception);
2209  mean_view=AcquireAuthenticCacheView(mean_image,exception);
2210 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2211  #pragma omp parallel for schedule(static) shared(status,progress) \
2212  magick_number_threads(mean_image,mean_image,mean_image->rows,1)
2213 #endif
2214  for (y=0; y < (ssize_t) mean_image->rows; y++)
2215  {
2216  const IndexPacket
2217  *magick_restrict indexes;
2218 
2219  const PixelPacket
2220  *magick_restrict p;
2221 
2222  PixelPacket
2223  *magick_restrict q;
2224 
2225  ssize_t
2226  x;
2227 
2228  if (status == MagickFalse)
2229  continue;
2230  p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2231  q=GetCacheViewAuthenticPixels(mean_view,0,y,mean_image->columns,1,
2232  exception);
2233  if ((p == (const PixelPacket *) NULL) || (q == (PixelPacket *) NULL))
2234  {
2235  status=MagickFalse;
2236  continue;
2237  }
2238  indexes=GetCacheViewVirtualIndexQueue(image_view);
2239  for (x=0; x < (ssize_t) mean_image->columns; x++)
2240  {
2241  MagickPixelPacket
2242  mean_pixel,
2243  previous_pixel;
2244 
2245  PointInfo
2246  mean_location,
2247  previous_location;
2248 
2249  ssize_t
2250  i;
2251 
2252  GetMagickPixelPacket(image,&mean_pixel);
2253  SetMagickPixelPacket(image,p,indexes+x,&mean_pixel);
2254  mean_location.x=(double) x;
2255  mean_location.y=(double) y;
2256  for (i=0; i < MaxMeanShiftIterations; i++)
2257  {
2258  double
2259  distance,
2260  gamma = 1.0;
2261 
2262  MagickPixelPacket
2263  sum_pixel;
2264 
2265  PointInfo
2266  sum_location;
2267 
2268  ssize_t
2269  count,
2270  v;
2271 
2272  sum_location.x=0.0;
2273  sum_location.y=0.0;
2274  GetMagickPixelPacket(image,&sum_pixel);
2275  previous_location=mean_location;
2276  previous_pixel=mean_pixel;
2277  count=0;
2278  for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
2279  {
2280  ssize_t
2281  u;
2282 
2283  for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
2284  {
2285  if ((v*v+u*u) <= (ssize_t) ((width/2)*(height/2)))
2286  {
2287  PixelPacket
2288  pixel;
2289 
2290  status=GetOneCacheViewVirtualPixel(pixel_view,(ssize_t)
2291  MagickRound(mean_location.x+u),(ssize_t) MagickRound(
2292  mean_location.y+v),&pixel,exception);
2293  distance=((MagickRealType) mean_pixel.red-(MagickRealType)
2294  pixel.red)*((MagickRealType) mean_pixel.red-(MagickRealType)
2295  pixel.red)+((MagickRealType) mean_pixel.green-
2296  (MagickRealType) pixel.green)*((MagickRealType)
2297  mean_pixel.green-(MagickRealType) pixel.green)+
2298  ((MagickRealType) mean_pixel.blue-(MagickRealType)
2299  pixel.blue)*((MagickRealType) mean_pixel.blue-
2300  (MagickRealType) pixel.blue);
2301  if (distance <= (color_distance*color_distance))
2302  {
2303  sum_location.x+=mean_location.x+u;
2304  sum_location.y+=mean_location.y+v;
2305  sum_pixel.red+=(MagickRealType) pixel.red;
2306  sum_pixel.green+=(MagickRealType) pixel.green;
2307  sum_pixel.blue+=(MagickRealType) pixel.blue;
2308  sum_pixel.opacity+=(MagickRealType) pixel.opacity;
2309  count++;
2310  }
2311  }
2312  }
2313  }
2314  if (count != 0)
2315  gamma=MagickSafeReciprocal((double) count);
2316  mean_location.x=gamma*sum_location.x;
2317  mean_location.y=gamma*sum_location.y;
2318  mean_pixel.red=gamma*sum_pixel.red;
2319  mean_pixel.green=gamma*sum_pixel.green;
2320  mean_pixel.blue=gamma*sum_pixel.blue;
2321  mean_pixel.opacity=gamma*sum_pixel.opacity;
2322  distance=(mean_location.x-previous_location.x)*
2323  (mean_location.x-previous_location.x)+
2324  (mean_location.y-previous_location.y)*
2325  (mean_location.y-previous_location.y)+
2326  255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)*
2327  255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)+
2328  255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)*
2329  255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)+
2330  255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue)*
2331  255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue);
2332  if (distance <= 3.0)
2333  break;
2334  }
2335  q->red=ClampToQuantum(mean_pixel.red);
2336  q->green=ClampToQuantum(mean_pixel.green);
2337  q->blue=ClampToQuantum(mean_pixel.blue);
2338  q->opacity=ClampToQuantum(mean_pixel.opacity);
2339  p++;
2340  q++;
2341  }
2342  if (SyncCacheViewAuthenticPixels(mean_view,exception) == MagickFalse)
2343  status=MagickFalse;
2344  if (image->progress_monitor != (MagickProgressMonitor) NULL)
2345  {
2346  MagickBooleanType
2347  proceed;
2348 
2349 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2350  #pragma omp atomic
2351 #endif
2352  progress++;
2353  proceed=SetImageProgress(image,MeanShiftImageTag,progress,image->rows);
2354  if (proceed == MagickFalse)
2355  status=MagickFalse;
2356  }
2357  }
2358  mean_view=DestroyCacheView(mean_view);
2359  pixel_view=DestroyCacheView(pixel_view);
2360  image_view=DestroyCacheView(image_view);
2361  return(mean_image);
2362 }
_CacheView
Definition: cache-view.c:66
_CannyInfo
Definition: feature.c:136
_ChannelFeatures
Definition: feature.h:29
_ChannelStatistics
Definition: statistic.h:31
_DoublePixelPacket
Definition: pixel.h:85
_DrawInfo
Definition: draw.h:212
_ExceptionInfo
Definition: exception.h:103
_ImageInfo
Definition: image.h:346
_Image
Definition: image.h:134
_LongPixelPacket
Definition: pixel.h:95
_MagickPixelPacket
Definition: pixel.h:105
_MatrixInfo
Definition: matrix.c:60
_PixelPacket
Definition: pixel.h:132
_PointInfo
Definition: draw.h:140
_SegmentInfo
Definition: image.h:87