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