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