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