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