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