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