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