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