2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 %
13 % MagickCore Image Feature Methods %
20 % Copyright 1999-2014 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
26 % http://www.imagemagick.org/script/license.php %
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. %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43 #include "MagickCore/studio.h"
44 #include "MagickCore/property.h"
45 #include "MagickCore/animate.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/client.h"
52 #include "MagickCore/color.h"
53 #include "MagickCore/color-private.h"
54 #include "MagickCore/colorspace.h"
55 #include "MagickCore/colorspace-private.h"
56 #include "MagickCore/composite.h"
57 #include "MagickCore/composite-private.h"
58 #include "MagickCore/compress.h"
59 #include "MagickCore/constitute.h"
60 #include "MagickCore/display.h"
61 #include "MagickCore/draw.h"
62 #include "MagickCore/enhance.h"
63 #include "MagickCore/exception.h"
64 #include "MagickCore/exception-private.h"
65 #include "MagickCore/feature.h"
66 #include "MagickCore/gem.h"
67 #include "MagickCore/geometry.h"
68 #include "MagickCore/list.h"
69 #include "MagickCore/image-private.h"
70 #include "MagickCore/magic.h"
71 #include "MagickCore/magick.h"
72 #include "MagickCore/matrix.h"
73 #include "MagickCore/memory_.h"
74 #include "MagickCore/module.h"
75 #include "MagickCore/monitor.h"
76 #include "MagickCore/monitor-private.h"
77 #include "MagickCore/morphology-private.h"
78 #include "MagickCore/option.h"
79 #include "MagickCore/paint.h"
80 #include "MagickCore/pixel-accessor.h"
81 #include "MagickCore/profile.h"
82 #include "MagickCore/quantize.h"
83 #include "MagickCore/quantum-private.h"
84 #include "MagickCore/random_.h"
85 #include "MagickCore/resource_.h"
86 #include "MagickCore/segment.h"
87 #include "MagickCore/semaphore.h"
88 #include "MagickCore/signature-private.h"
89 #include "MagickCore/string_.h"
90 #include "MagickCore/thread-private.h"
91 #include "MagickCore/timer.h"
92 #include "MagickCore/utility.h"
93 #include "MagickCore/version.h"
96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
100 % C a n n y E d g e I m a g e %
104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
106 % CannyEdgeImage() uses a multi-stage algorithm to detect a wide range of
109 % The format of the EdgeImage method is:
111 % Image *CannyEdgeImage(const Image *image,const double radius,
112 % const double sigma,const double lower_percent,
113 % const double upper_percent,ExceptionInfo *exception)
115 % A description of each parameter follows:
117 % o image: the image.
119 % o radius: the radius of the gaussian smoothing filter.
121 % o sigma: the sigma of the gaussian smoothing filter.
123 % o lower_precent: percentage of edge pixels in the lower threshold.
125 % o upper_percent: percentage of edge pixels in the upper threshold.
127 % o exception: return any errors or warnings in this structure.
131 typedef struct _CannyInfo
145 static inline MagickBooleanType IsAuthenticPixel(const Image *image,
146 const ssize_t x,const ssize_t y)
148 if ((x < 0) || (x >= (ssize_t) image->columns))
150 if ((y < 0) || (y >= (ssize_t) image->rows))
155 static MagickBooleanType TraceEdges(Image *edge_image,CacheView *edge_view,
156 MatrixInfo *canny_cache,const ssize_t x,const ssize_t y,
157 const double lower_threshold,ExceptionInfo *exception)
172 q=GetCacheViewAuthenticPixels(edge_view,x,y,1,1,exception);
173 if (q == (Quantum *) NULL)
176 status=SyncCacheViewAuthenticPixels(edge_view,exception);
177 if (status == MagickFalse)
178 return(MagickFalse);;
179 if (GetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
183 if (SetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
191 status=GetMatrixElement(canny_cache,i,0,&edge);
192 if (status == MagickFalse)
194 for (v=(-1); v <= 1; v++)
199 for (u=(-1); u <= 1; u++)
201 if ((u == 0) && (v == 0))
203 if (IsAuthenticPixel(edge_image,edge.x+u,edge.y+v) == MagickFalse)
206 Not an edge if gradient value is below the lower threshold.
208 q=GetCacheViewAuthenticPixels(edge_view,edge.x+u,edge.y+v,1,1,
210 if (q == (Quantum *) NULL)
212 status=GetMatrixElement(canny_cache,edge.x+u,edge.y+v,&pixel);
213 if (status == MagickFalse)
215 if ((GetPixelIntensity(edge_image,q) == 0.0) &&
216 (pixel.intensity >= lower_threshold))
219 status=SyncCacheViewAuthenticPixels(edge_view,exception);
220 if (status == MagickFalse)
224 status=SetMatrixElement(canny_cache,i,0,&edge);
225 if (status == MagickFalse)
235 MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
236 const double sigma,const double lower_percent,const double upper_percent,
237 ExceptionInfo *exception)
246 geometry[MaxTextExtent];
269 assert(image != (const Image *) NULL);
270 assert(image->signature == MagickSignature);
271 if (image->debug != MagickFalse)
272 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
273 assert(exception != (ExceptionInfo *) NULL);
274 assert(exception->signature == MagickSignature);
278 (void) FormatLocaleString(geometry,MaxTextExtent,
279 "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
280 kernel_info=AcquireKernelInfo(geometry);
281 if (kernel_info == (KernelInfo *) NULL)
282 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
283 edge_image=MorphologyApply(image,ConvolveMorphology,1,kernel_info,
284 UndefinedCompositeOp,0.0,exception);
285 kernel_info=DestroyKernelInfo(kernel_info);
286 if (edge_image == (Image *) NULL)
287 return((Image *) NULL);
288 if (SetImageColorspace(edge_image,GRAYColorspace,exception) == MagickFalse)
290 edge_image=DestroyImage(edge_image);
291 return((Image *) NULL);
294 Find the intensity gradient of the image.
296 canny_cache=AcquireMatrixInfo(edge_image->columns,edge_image->rows,
297 sizeof(CannyInfo),exception);
298 if (canny_cache == (MatrixInfo *) NULL)
300 edge_image=DestroyImage(edge_image);
301 return((Image *) NULL);
304 edge_view=AcquireVirtualCacheView(edge_image,exception);
305 #if defined(MAGICKCORE_OPENMP_SUPPORT)
306 #pragma omp parallel for schedule(static,4) shared(status) \
307 magick_threads(edge_image,edge_image,edge_image->rows,1)
309 for (y=0; y < (ssize_t) edge_image->rows; y++)
311 register const Quantum
317 if (status == MagickFalse)
319 p=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns+1,2,
321 if (p == (const Quantum *) NULL)
326 for (x=0; x < (ssize_t) edge_image->columns; x++)
335 register const Quantum
336 *restrict kernel_pixels;
353 (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
357 for (v=0; v < 2; v++)
362 for (u=0; u < 2; u++)
367 intensity=GetPixelIntensity(edge_image,kernel_pixels+u);
368 dx+=0.5*Gx[v][u]*intensity;
369 dy+=0.5*Gy[v][u]*intensity;
371 kernel_pixels+=edge_image->columns+1;
373 pixel.magnitude=hypot(dx,dy);
375 if (fabs(dx) > MagickEpsilon)
383 if (slope < -2.41421356237)
386 if (slope < -0.414213562373)
393 if (slope > 2.41421356237)
396 if (slope > 0.414213562373)
402 if (SetMatrixElement(canny_cache,x,y,&pixel) == MagickFalse)
404 p+=GetPixelChannels(edge_image);
407 edge_view=DestroyCacheView(edge_view);
409 Non-maxima suppression, remove pixels that are not considered to be part
412 (void) GetMatrixElement(canny_cache,0,0,&pixel);
415 edge_view=AcquireAuthenticCacheView(edge_image,exception);
416 #if defined(MAGICKCORE_OPENMP_SUPPORT)
417 #pragma omp parallel for schedule(static,4) shared(status) \
418 magick_threads(edge_image,edge_image,edge_image->rows,1)
420 for (y=0; y < (ssize_t) edge_image->rows; y++)
428 if (status == MagickFalse)
430 q=GetCacheViewAuthenticPixels(edge_view,0,y,edge_image->columns,1,
432 if (q == (Quantum *) NULL)
437 for (x=0; x < (ssize_t) edge_image->columns; x++)
444 (void) GetMatrixElement(canny_cache,x,y,&pixel);
445 switch (pixel.orientation)
451 0 degrees, north and south.
453 (void) GetMatrixElement(canny_cache,x,y-1,&alpha_pixel);
454 (void) GetMatrixElement(canny_cache,x,y+1,&beta_pixel);
460 45 degrees, northwest and southeast.
462 (void) GetMatrixElement(canny_cache,x-1,y-1,&alpha_pixel);
463 (void) GetMatrixElement(canny_cache,x+1,y+1,&beta_pixel);
469 90 degrees, east and west.
471 (void) GetMatrixElement(canny_cache,x-1,y,&alpha_pixel);
472 (void) GetMatrixElement(canny_cache,x+1,y,&beta_pixel);
478 135 degrees, northeast and southwest.
480 (void) GetMatrixElement(canny_cache,x+1,y-1,&beta_pixel);
481 (void) GetMatrixElement(canny_cache,x-1,y+1,&alpha_pixel);
485 pixel.intensity=pixel.magnitude;
486 if ((pixel.magnitude < alpha_pixel.magnitude) ||
487 (pixel.magnitude < beta_pixel.magnitude))
489 (void) SetMatrixElement(canny_cache,x,y,&pixel);
490 #if defined(MAGICKCORE_OPENMP_SUPPORT)
491 #pragma omp critical (MagickCore_CannyEdgeImage)
494 if (pixel.intensity < min)
496 if (pixel.intensity > max)
500 q+=GetPixelChannels(edge_image);
502 if (SyncCacheViewAuthenticPixels(edge_view,exception) == MagickFalse)
505 edge_view=DestroyCacheView(edge_view);
507 Estimate hysteresis threshold.
509 lower_threshold=lower_percent*(max-min)+min;
510 upper_threshold=upper_percent*(max-min)+min;
512 Hysteresis threshold.
514 edge_view=AcquireAuthenticCacheView(edge_image,exception);
515 for (y=0; y < (ssize_t) edge_image->rows; y++)
520 if (status == MagickFalse)
522 for (x=0; x < (ssize_t) edge_image->columns; x++)
527 register const Quantum
531 Edge if pixel gradient higher than upper threshold.
533 p=GetCacheViewVirtualPixels(edge_view,x,y,1,1,exception);
534 if (p == (const Quantum *) NULL)
536 status=GetMatrixElement(canny_cache,x,y,&pixel);
537 if (status == MagickFalse)
539 if ((GetPixelIntensity(edge_image,p) == 0.0) &&
540 (pixel.intensity >= upper_threshold))
541 status=TraceEdges(edge_image,edge_view,canny_cache,x,y,lower_threshold,
545 edge_view=DestroyCacheView(edge_view);
549 canny_cache=DestroyMatrixInfo(canny_cache);
554 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
558 % H o u g h L i n e s I m a g e %
562 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
564 % HoughLinesImage() identifies lines in an image.
566 % The format of the HoughLinesImage method is:
568 % Image *HoughLinesImage(const Image *image,const size_t width,
569 % const size_t height,const size_t threshold,ExceptionInfo *exception)
571 % A description of each parameter follows:
573 % o image: the image.
575 % o width, height: find line pairs as local maxima in this neighborhood.
577 % o threshold: the line count threshold.
579 % o exception: return any errors or warnings in this structure.
582 MagickExport Image *HoughLinesImage(const Image *image,const size_t width,
583 const size_t height,const size_t threshold,ExceptionInfo *exception)
585 return((Image *) NULL);
589 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
593 % G e t I m a g e F e a t u r e s %
597 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
599 % GetImageFeatures() returns features for each channel in the image in
600 % each of four directions (horizontal, vertical, left and right diagonals)
601 % for the specified distance. The features include the angular second
602 % moment, contrast, correlation, sum of squares: variance, inverse difference
603 % moment, sum average, sum varience, sum entropy, entropy, difference variance,% difference entropy, information measures of correlation 1, information
604 % measures of correlation 2, and maximum correlation coefficient. You can
605 % access the red channel contrast, for example, like this:
607 % channel_features=GetImageFeatures(image,1,exception);
608 % contrast=channel_features[RedPixelChannel].contrast[0];
610 % Use MagickRelinquishMemory() to free the features buffer.
612 % The format of the GetImageFeatures method is:
614 % ChannelFeatures *GetImageFeatures(const Image *image,
615 % const size_t distance,ExceptionInfo *exception)
617 % A description of each parameter follows:
619 % o image: the image.
621 % o distance: the distance.
623 % o exception: return any errors or warnings in this structure.
627 static inline ssize_t MagickAbsoluteValue(const ssize_t x)
634 static inline double MagickLog10(const double x)
636 #define Log10Epsilon (1.0e-11)
638 if (fabs(x) < Log10Epsilon)
639 return(log10(Log10Epsilon));
640 return(log10(fabs(x)));
643 MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
644 const size_t distance,ExceptionInfo *exception)
646 typedef struct _ChannelStatistics
649 direction[4]; /* horizontal, vertical, left and right diagonals */
694 assert(image != (Image *) NULL);
695 assert(image->signature == MagickSignature);
696 if (image->debug != MagickFalse)
697 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
698 if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
699 return((ChannelFeatures *) NULL);
700 length=CompositeChannels+1UL;
701 channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
702 sizeof(*channel_features));
703 if (channel_features == (ChannelFeatures *) NULL)
704 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
705 (void) ResetMagickMemory(channel_features,0,length*
706 sizeof(*channel_features));
710 grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays));
711 if (grays == (PixelPacket *) NULL)
713 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
715 (void) ThrowMagickException(exception,GetMagickModule(),
716 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
717 return(channel_features);
719 for (i=0; i <= (ssize_t) MaxMap; i++)
722 grays[i].green=(~0U);
724 grays[i].alpha=(~0U);
725 grays[i].black=(~0U);
728 image_view=AcquireVirtualCacheView(image,exception);
729 #if defined(MAGICKCORE_OPENMP_SUPPORT)
730 #pragma omp parallel for schedule(static,4) shared(status) \
731 magick_threads(image,image,image->rows,1)
733 for (y=0; y < (ssize_t) image->rows; y++)
735 register const Quantum
741 if (status == MagickFalse)
743 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
744 if (p == (const Quantum *) NULL)
749 for (x=0; x < (ssize_t) image->columns; x++)
751 grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
752 ScaleQuantumToMap(GetPixelRed(image,p));
753 grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
754 ScaleQuantumToMap(GetPixelGreen(image,p));
755 grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
756 ScaleQuantumToMap(GetPixelBlue(image,p));
757 if (image->colorspace == CMYKColorspace)
758 grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
759 ScaleQuantumToMap(GetPixelBlack(image,p));
760 if (image->alpha_trait == BlendPixelTrait)
761 grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
762 ScaleQuantumToMap(GetPixelAlpha(image,p));
763 p+=GetPixelChannels(image);
766 image_view=DestroyCacheView(image_view);
767 if (status == MagickFalse)
769 grays=(PixelPacket *) RelinquishMagickMemory(grays);
770 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
772 return(channel_features);
774 (void) ResetMagickMemory(&gray,0,sizeof(gray));
775 for (i=0; i <= (ssize_t) MaxMap; i++)
777 if (grays[i].red != ~0U)
778 grays[gray.red++].red=grays[i].red;
779 if (grays[i].green != ~0U)
780 grays[gray.green++].green=grays[i].green;
781 if (grays[i].blue != ~0U)
782 grays[gray.blue++].blue=grays[i].blue;
783 if (image->colorspace == CMYKColorspace)
784 if (grays[i].black != ~0U)
785 grays[gray.black++].black=grays[i].black;
786 if (image->alpha_trait == BlendPixelTrait)
787 if (grays[i].alpha != ~0U)
788 grays[gray.alpha++].alpha=grays[i].alpha;
791 Allocate spatial dependence matrix.
793 number_grays=gray.red;
794 if (gray.green > number_grays)
795 number_grays=gray.green;
796 if (gray.blue > number_grays)
797 number_grays=gray.blue;
798 if (image->colorspace == CMYKColorspace)
799 if (gray.black > number_grays)
800 number_grays=gray.black;
801 if (image->alpha_trait == BlendPixelTrait)
802 if (gray.alpha > number_grays)
803 number_grays=gray.alpha;
804 cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays,
805 sizeof(*cooccurrence));
806 density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
808 density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
809 sizeof(*density_xy));
810 density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
812 Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q));
813 sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum));
814 if ((cooccurrence == (ChannelStatistics **) NULL) ||
815 (density_x == (ChannelStatistics *) NULL) ||
816 (density_xy == (ChannelStatistics *) NULL) ||
817 (density_y == (ChannelStatistics *) NULL) ||
818 (Q == (ChannelStatistics **) NULL) ||
819 (sum == (ChannelStatistics *) NULL))
821 if (Q != (ChannelStatistics **) NULL)
823 for (i=0; i < (ssize_t) number_grays; i++)
824 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
825 Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
827 if (sum != (ChannelStatistics *) NULL)
828 sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
829 if (density_y != (ChannelStatistics *) NULL)
830 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
831 if (density_xy != (ChannelStatistics *) NULL)
832 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
833 if (density_x != (ChannelStatistics *) NULL)
834 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
835 if (cooccurrence != (ChannelStatistics **) NULL)
837 for (i=0; i < (ssize_t) number_grays; i++)
838 cooccurrence[i]=(ChannelStatistics *)
839 RelinquishMagickMemory(cooccurrence[i]);
840 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
843 grays=(PixelPacket *) RelinquishMagickMemory(grays);
844 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
846 (void) ThrowMagickException(exception,GetMagickModule(),
847 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
848 return(channel_features);
850 (void) ResetMagickMemory(&correlation,0,sizeof(correlation));
851 (void) ResetMagickMemory(density_x,0,2*(number_grays+1)*sizeof(*density_x));
852 (void) ResetMagickMemory(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
853 (void) ResetMagickMemory(density_y,0,2*(number_grays+1)*sizeof(*density_y));
854 (void) ResetMagickMemory(&mean,0,sizeof(mean));
855 (void) ResetMagickMemory(sum,0,number_grays*sizeof(*sum));
856 (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
857 (void) ResetMagickMemory(density_xy,0,2*number_grays*sizeof(*density_xy));
858 (void) ResetMagickMemory(&entropy_x,0,sizeof(entropy_x));
859 (void) ResetMagickMemory(&entropy_xy,0,sizeof(entropy_xy));
860 (void) ResetMagickMemory(&entropy_xy1,0,sizeof(entropy_xy1));
861 (void) ResetMagickMemory(&entropy_xy2,0,sizeof(entropy_xy2));
862 (void) ResetMagickMemory(&entropy_y,0,sizeof(entropy_y));
863 (void) ResetMagickMemory(&variance,0,sizeof(variance));
864 for (i=0; i < (ssize_t) number_grays; i++)
866 cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
867 sizeof(**cooccurrence));
868 Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
869 if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
870 (Q[i] == (ChannelStatistics *) NULL))
872 (void) ResetMagickMemory(cooccurrence[i],0,number_grays*
873 sizeof(**cooccurrence));
874 (void) ResetMagickMemory(Q[i],0,number_grays*sizeof(**Q));
876 if (i < (ssize_t) number_grays)
878 for (i--; i >= 0; i--)
880 if (Q[i] != (ChannelStatistics *) NULL)
881 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
882 if (cooccurrence[i] != (ChannelStatistics *) NULL)
883 cooccurrence[i]=(ChannelStatistics *)
884 RelinquishMagickMemory(cooccurrence[i]);
886 Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
887 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
888 sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
889 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
890 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
891 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
892 grays=(PixelPacket *) RelinquishMagickMemory(grays);
893 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
895 (void) ThrowMagickException(exception,GetMagickModule(),
896 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
897 return(channel_features);
900 Initialize spatial dependence matrix.
903 image_view=AcquireVirtualCacheView(image,exception);
904 for (y=0; y < (ssize_t) image->rows; y++)
906 register const Quantum
918 if (status == MagickFalse)
920 p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,y,image->columns+
921 2*distance,distance+2,exception);
922 if (p == (const Quantum *) NULL)
927 p+=distance*GetPixelChannels(image);;
928 for (x=0; x < (ssize_t) image->columns; x++)
930 for (i=0; i < 4; i++)
938 Horizontal adjacency.
940 offset=(ssize_t) distance;
948 offset=(ssize_t) (image->columns+2*distance);
954 Right diagonal adjacency.
956 offset=(ssize_t) ((image->columns+2*distance)-distance);
962 Left diagonal adjacency.
964 offset=(ssize_t) ((image->columns+2*distance)+distance);
970 while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p)))
972 while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image))))
974 cooccurrence[u][v].direction[i].red++;
975 cooccurrence[v][u].direction[i].red++;
978 while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p)))
980 while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image))))
982 cooccurrence[u][v].direction[i].green++;
983 cooccurrence[v][u].direction[i].green++;
986 while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p)))
988 while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image))))
990 cooccurrence[u][v].direction[i].blue++;
991 cooccurrence[v][u].direction[i].blue++;
992 if (image->colorspace == CMYKColorspace)
996 while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p)))
998 while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image))))
1000 cooccurrence[u][v].direction[i].black++;
1001 cooccurrence[v][u].direction[i].black++;
1003 if (image->alpha_trait == BlendPixelTrait)
1007 while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p)))
1009 while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image))))
1011 cooccurrence[u][v].direction[i].alpha++;
1012 cooccurrence[v][u].direction[i].alpha++;
1015 p+=GetPixelChannels(image);
1018 grays=(PixelPacket *) RelinquishMagickMemory(grays);
1019 image_view=DestroyCacheView(image_view);
1020 if (status == MagickFalse)
1022 for (i=0; i < (ssize_t) number_grays; i++)
1023 cooccurrence[i]=(ChannelStatistics *)
1024 RelinquishMagickMemory(cooccurrence[i]);
1025 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1026 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1028 (void) ThrowMagickException(exception,GetMagickModule(),
1029 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1030 return(channel_features);
1033 Normalize spatial dependence matrix.
1035 for (i=0; i < 4; i++)
1049 Horizontal adjacency.
1051 normalize=2.0*image->rows*(image->columns-distance);
1059 normalize=2.0*(image->rows-distance)*image->columns;
1065 Right diagonal adjacency.
1067 normalize=2.0*(image->rows-distance)*(image->columns-distance);
1073 Left diagonal adjacency.
1075 normalize=2.0*(image->rows-distance)*(image->columns-distance);
1079 normalize=PerceptibleReciprocal(normalize);
1080 for (y=0; y < (ssize_t) number_grays; y++)
1085 for (x=0; x < (ssize_t) number_grays; x++)
1087 cooccurrence[x][y].direction[i].red*=normalize;
1088 cooccurrence[x][y].direction[i].green*=normalize;
1089 cooccurrence[x][y].direction[i].blue*=normalize;
1090 if (image->colorspace == CMYKColorspace)
1091 cooccurrence[x][y].direction[i].black*=normalize;
1092 if (image->alpha_trait == BlendPixelTrait)
1093 cooccurrence[x][y].direction[i].alpha*=normalize;
1098 Compute texture features.
1100 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1101 #pragma omp parallel for schedule(static,4) shared(status) \
1102 magick_threads(image,image,number_grays,1)
1104 for (i=0; i < 4; i++)
1109 for (y=0; y < (ssize_t) number_grays; y++)
1114 for (x=0; x < (ssize_t) number_grays; x++)
1117 Angular second moment: measure of homogeneity of the image.
1119 channel_features[RedPixelChannel].angular_second_moment[i]+=
1120 cooccurrence[x][y].direction[i].red*
1121 cooccurrence[x][y].direction[i].red;
1122 channel_features[GreenPixelChannel].angular_second_moment[i]+=
1123 cooccurrence[x][y].direction[i].green*
1124 cooccurrence[x][y].direction[i].green;
1125 channel_features[BluePixelChannel].angular_second_moment[i]+=
1126 cooccurrence[x][y].direction[i].blue*
1127 cooccurrence[x][y].direction[i].blue;
1128 if (image->colorspace == CMYKColorspace)
1129 channel_features[BlackPixelChannel].angular_second_moment[i]+=
1130 cooccurrence[x][y].direction[i].black*
1131 cooccurrence[x][y].direction[i].black;
1132 if (image->alpha_trait == BlendPixelTrait)
1133 channel_features[AlphaPixelChannel].angular_second_moment[i]+=
1134 cooccurrence[x][y].direction[i].alpha*
1135 cooccurrence[x][y].direction[i].alpha;
1137 Correlation: measure of linear-dependencies in the image.
1139 sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1140 sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1141 sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1142 if (image->colorspace == CMYKColorspace)
1143 sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black;
1144 if (image->alpha_trait == BlendPixelTrait)
1145 sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha;
1146 correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
1147 correlation.direction[i].green+=x*y*
1148 cooccurrence[x][y].direction[i].green;
1149 correlation.direction[i].blue+=x*y*
1150 cooccurrence[x][y].direction[i].blue;
1151 if (image->colorspace == CMYKColorspace)
1152 correlation.direction[i].black+=x*y*
1153 cooccurrence[x][y].direction[i].black;
1154 if (image->alpha_trait == BlendPixelTrait)
1155 correlation.direction[i].alpha+=x*y*
1156 cooccurrence[x][y].direction[i].alpha;
1158 Inverse Difference Moment.
1160 channel_features[RedPixelChannel].inverse_difference_moment[i]+=
1161 cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
1162 channel_features[GreenPixelChannel].inverse_difference_moment[i]+=
1163 cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
1164 channel_features[BluePixelChannel].inverse_difference_moment[i]+=
1165 cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
1166 if (image->colorspace == CMYKColorspace)
1167 channel_features[BlackPixelChannel].inverse_difference_moment[i]+=
1168 cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1);
1169 if (image->alpha_trait == BlendPixelTrait)
1170 channel_features[AlphaPixelChannel].inverse_difference_moment[i]+=
1171 cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
1175 density_xy[y+x+2].direction[i].red+=
1176 cooccurrence[x][y].direction[i].red;
1177 density_xy[y+x+2].direction[i].green+=
1178 cooccurrence[x][y].direction[i].green;
1179 density_xy[y+x+2].direction[i].blue+=
1180 cooccurrence[x][y].direction[i].blue;
1181 if (image->colorspace == CMYKColorspace)
1182 density_xy[y+x+2].direction[i].black+=
1183 cooccurrence[x][y].direction[i].black;
1184 if (image->alpha_trait == BlendPixelTrait)
1185 density_xy[y+x+2].direction[i].alpha+=
1186 cooccurrence[x][y].direction[i].alpha;
1190 channel_features[RedPixelChannel].entropy[i]-=
1191 cooccurrence[x][y].direction[i].red*
1192 MagickLog10(cooccurrence[x][y].direction[i].red);
1193 channel_features[GreenPixelChannel].entropy[i]-=
1194 cooccurrence[x][y].direction[i].green*
1195 MagickLog10(cooccurrence[x][y].direction[i].green);
1196 channel_features[BluePixelChannel].entropy[i]-=
1197 cooccurrence[x][y].direction[i].blue*
1198 MagickLog10(cooccurrence[x][y].direction[i].blue);
1199 if (image->colorspace == CMYKColorspace)
1200 channel_features[BlackPixelChannel].entropy[i]-=
1201 cooccurrence[x][y].direction[i].black*
1202 MagickLog10(cooccurrence[x][y].direction[i].black);
1203 if (image->alpha_trait == BlendPixelTrait)
1204 channel_features[AlphaPixelChannel].entropy[i]-=
1205 cooccurrence[x][y].direction[i].alpha*
1206 MagickLog10(cooccurrence[x][y].direction[i].alpha);
1208 Information Measures of Correlation.
1210 density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
1211 density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
1212 density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1213 if (image->alpha_trait == BlendPixelTrait)
1214 density_x[x].direction[i].alpha+=
1215 cooccurrence[x][y].direction[i].alpha;
1216 if (image->colorspace == CMYKColorspace)
1217 density_x[x].direction[i].black+=
1218 cooccurrence[x][y].direction[i].black;
1219 density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1220 density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1221 density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1222 if (image->colorspace == CMYKColorspace)
1223 density_y[y].direction[i].black+=
1224 cooccurrence[x][y].direction[i].black;
1225 if (image->alpha_trait == BlendPixelTrait)
1226 density_y[y].direction[i].alpha+=
1227 cooccurrence[x][y].direction[i].alpha;
1229 mean.direction[i].red+=y*sum[y].direction[i].red;
1230 sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
1231 mean.direction[i].green+=y*sum[y].direction[i].green;
1232 sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
1233 mean.direction[i].blue+=y*sum[y].direction[i].blue;
1234 sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
1235 if (image->colorspace == CMYKColorspace)
1237 mean.direction[i].black+=y*sum[y].direction[i].black;
1238 sum_squares.direction[i].black+=y*y*sum[y].direction[i].black;
1240 if (image->alpha_trait == BlendPixelTrait)
1242 mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
1243 sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
1247 Correlation: measure of linear-dependencies in the image.
1249 channel_features[RedPixelChannel].correlation[i]=
1250 (correlation.direction[i].red-mean.direction[i].red*
1251 mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
1252 (mean.direction[i].red*mean.direction[i].red))*sqrt(
1253 sum_squares.direction[i].red-(mean.direction[i].red*
1254 mean.direction[i].red)));
1255 channel_features[GreenPixelChannel].correlation[i]=
1256 (correlation.direction[i].green-mean.direction[i].green*
1257 mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
1258 (mean.direction[i].green*mean.direction[i].green))*sqrt(
1259 sum_squares.direction[i].green-(mean.direction[i].green*
1260 mean.direction[i].green)));
1261 channel_features[BluePixelChannel].correlation[i]=
1262 (correlation.direction[i].blue-mean.direction[i].blue*
1263 mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
1264 (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
1265 sum_squares.direction[i].blue-(mean.direction[i].blue*
1266 mean.direction[i].blue)));
1267 if (image->colorspace == CMYKColorspace)
1268 channel_features[BlackPixelChannel].correlation[i]=
1269 (correlation.direction[i].black-mean.direction[i].black*
1270 mean.direction[i].black)/(sqrt(sum_squares.direction[i].black-
1271 (mean.direction[i].black*mean.direction[i].black))*sqrt(
1272 sum_squares.direction[i].black-(mean.direction[i].black*
1273 mean.direction[i].black)));
1274 if (image->alpha_trait == BlendPixelTrait)
1275 channel_features[AlphaPixelChannel].correlation[i]=
1276 (correlation.direction[i].alpha-mean.direction[i].alpha*
1277 mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
1278 (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt(
1279 sum_squares.direction[i].alpha-(mean.direction[i].alpha*
1280 mean.direction[i].alpha)));
1283 Compute more texture features.
1285 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1286 #pragma omp parallel for schedule(static,4) shared(status) \
1287 magick_threads(image,image,number_grays,1)
1289 for (i=0; i < 4; i++)
1294 for (x=2; x < (ssize_t) (2*number_grays); x++)
1299 channel_features[RedPixelChannel].sum_average[i]+=
1300 x*density_xy[x].direction[i].red;
1301 channel_features[GreenPixelChannel].sum_average[i]+=
1302 x*density_xy[x].direction[i].green;
1303 channel_features[BluePixelChannel].sum_average[i]+=
1304 x*density_xy[x].direction[i].blue;
1305 if (image->colorspace == CMYKColorspace)
1306 channel_features[BlackPixelChannel].sum_average[i]+=
1307 x*density_xy[x].direction[i].black;
1308 if (image->alpha_trait == BlendPixelTrait)
1309 channel_features[AlphaPixelChannel].sum_average[i]+=
1310 x*density_xy[x].direction[i].alpha;
1314 channel_features[RedPixelChannel].sum_entropy[i]-=
1315 density_xy[x].direction[i].red*
1316 MagickLog10(density_xy[x].direction[i].red);
1317 channel_features[GreenPixelChannel].sum_entropy[i]-=
1318 density_xy[x].direction[i].green*
1319 MagickLog10(density_xy[x].direction[i].green);
1320 channel_features[BluePixelChannel].sum_entropy[i]-=
1321 density_xy[x].direction[i].blue*
1322 MagickLog10(density_xy[x].direction[i].blue);
1323 if (image->colorspace == CMYKColorspace)
1324 channel_features[BlackPixelChannel].sum_entropy[i]-=
1325 density_xy[x].direction[i].black*
1326 MagickLog10(density_xy[x].direction[i].black);
1327 if (image->alpha_trait == BlendPixelTrait)
1328 channel_features[AlphaPixelChannel].sum_entropy[i]-=
1329 density_xy[x].direction[i].alpha*
1330 MagickLog10(density_xy[x].direction[i].alpha);
1334 channel_features[RedPixelChannel].sum_variance[i]+=
1335 (x-channel_features[RedPixelChannel].sum_entropy[i])*
1336 (x-channel_features[RedPixelChannel].sum_entropy[i])*
1337 density_xy[x].direction[i].red;
1338 channel_features[GreenPixelChannel].sum_variance[i]+=
1339 (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1340 (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1341 density_xy[x].direction[i].green;
1342 channel_features[BluePixelChannel].sum_variance[i]+=
1343 (x-channel_features[BluePixelChannel].sum_entropy[i])*
1344 (x-channel_features[BluePixelChannel].sum_entropy[i])*
1345 density_xy[x].direction[i].blue;
1346 if (image->colorspace == CMYKColorspace)
1347 channel_features[BlackPixelChannel].sum_variance[i]+=
1348 (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1349 (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1350 density_xy[x].direction[i].black;
1351 if (image->alpha_trait == BlendPixelTrait)
1352 channel_features[AlphaPixelChannel].sum_variance[i]+=
1353 (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1354 (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1355 density_xy[x].direction[i].alpha;
1359 Compute more texture features.
1361 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1362 #pragma omp parallel for schedule(static,4) shared(status) \
1363 magick_threads(image,image,number_grays,1)
1365 for (i=0; i < 4; i++)
1370 for (y=0; y < (ssize_t) number_grays; y++)
1375 for (x=0; x < (ssize_t) number_grays; x++)
1378 Sum of Squares: Variance
1380 variance.direction[i].red+=(y-mean.direction[i].red+1)*
1381 (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
1382 variance.direction[i].green+=(y-mean.direction[i].green+1)*
1383 (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
1384 variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
1385 (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
1386 if (image->colorspace == CMYKColorspace)
1387 variance.direction[i].black+=(y-mean.direction[i].black+1)*
1388 (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black;
1389 if (image->alpha_trait == BlendPixelTrait)
1390 variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
1391 (y-mean.direction[i].alpha+1)*
1392 cooccurrence[x][y].direction[i].alpha;
1394 Sum average / Difference Variance.
1396 density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
1397 cooccurrence[x][y].direction[i].red;
1398 density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
1399 cooccurrence[x][y].direction[i].green;
1400 density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
1401 cooccurrence[x][y].direction[i].blue;
1402 if (image->colorspace == CMYKColorspace)
1403 density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
1404 cooccurrence[x][y].direction[i].black;
1405 if (image->alpha_trait == BlendPixelTrait)
1406 density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
1407 cooccurrence[x][y].direction[i].alpha;
1409 Information Measures of Correlation.
1411 entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
1412 MagickLog10(cooccurrence[x][y].direction[i].red);
1413 entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
1414 MagickLog10(cooccurrence[x][y].direction[i].green);
1415 entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
1416 MagickLog10(cooccurrence[x][y].direction[i].blue);
1417 if (image->colorspace == CMYKColorspace)
1418 entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black*
1419 MagickLog10(cooccurrence[x][y].direction[i].black);
1420 if (image->alpha_trait == BlendPixelTrait)
1421 entropy_xy.direction[i].alpha-=
1422 cooccurrence[x][y].direction[i].alpha*MagickLog10(
1423 cooccurrence[x][y].direction[i].alpha);
1424 entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
1425 MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red));
1426 entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
1427 MagickLog10(density_x[x].direction[i].green*
1428 density_y[y].direction[i].green));
1429 entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
1430 MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue));
1431 if (image->colorspace == CMYKColorspace)
1432 entropy_xy1.direction[i].black-=(
1433 cooccurrence[x][y].direction[i].black*MagickLog10(
1434 density_x[x].direction[i].black*density_y[y].direction[i].black));
1435 if (image->alpha_trait == BlendPixelTrait)
1436 entropy_xy1.direction[i].alpha-=(
1437 cooccurrence[x][y].direction[i].alpha*MagickLog10(
1438 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1439 entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
1440 density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red*
1441 density_y[y].direction[i].red));
1442 entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
1443 density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green*
1444 density_y[y].direction[i].green));
1445 entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
1446 density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue*
1447 density_y[y].direction[i].blue));
1448 if (image->colorspace == CMYKColorspace)
1449 entropy_xy2.direction[i].black-=(density_x[x].direction[i].black*
1450 density_y[y].direction[i].black*MagickLog10(
1451 density_x[x].direction[i].black*density_y[y].direction[i].black));
1452 if (image->alpha_trait == BlendPixelTrait)
1453 entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha*
1454 density_y[y].direction[i].alpha*MagickLog10(
1455 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1458 channel_features[RedPixelChannel].variance_sum_of_squares[i]=
1459 variance.direction[i].red;
1460 channel_features[GreenPixelChannel].variance_sum_of_squares[i]=
1461 variance.direction[i].green;
1462 channel_features[BluePixelChannel].variance_sum_of_squares[i]=
1463 variance.direction[i].blue;
1464 if (image->colorspace == CMYKColorspace)
1465 channel_features[BlackPixelChannel].variance_sum_of_squares[i]=
1466 variance.direction[i].black;
1467 if (image->alpha_trait == BlendPixelTrait)
1468 channel_features[AlphaPixelChannel].variance_sum_of_squares[i]=
1469 variance.direction[i].alpha;
1472 Compute more texture features.
1474 (void) ResetMagickMemory(&variance,0,sizeof(variance));
1475 (void) ResetMagickMemory(&sum_squares,0,sizeof(sum_squares));
1476 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1477 #pragma omp parallel for schedule(static,4) shared(status) \
1478 magick_threads(image,image,number_grays,1)
1480 for (i=0; i < 4; i++)
1485 for (x=0; x < (ssize_t) number_grays; x++)
1488 Difference variance.
1490 variance.direction[i].red+=density_xy[x].direction[i].red;
1491 variance.direction[i].green+=density_xy[x].direction[i].green;
1492 variance.direction[i].blue+=density_xy[x].direction[i].blue;
1493 if (image->colorspace == CMYKColorspace)
1494 variance.direction[i].black+=density_xy[x].direction[i].black;
1495 if (image->alpha_trait == BlendPixelTrait)
1496 variance.direction[i].alpha+=density_xy[x].direction[i].alpha;
1497 sum_squares.direction[i].red+=density_xy[x].direction[i].red*
1498 density_xy[x].direction[i].red;
1499 sum_squares.direction[i].green+=density_xy[x].direction[i].green*
1500 density_xy[x].direction[i].green;
1501 sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
1502 density_xy[x].direction[i].blue;
1503 if (image->colorspace == CMYKColorspace)
1504 sum_squares.direction[i].black+=density_xy[x].direction[i].black*
1505 density_xy[x].direction[i].black;
1506 if (image->alpha_trait == BlendPixelTrait)
1507 sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
1508 density_xy[x].direction[i].alpha;
1512 channel_features[RedPixelChannel].difference_entropy[i]-=
1513 density_xy[x].direction[i].red*
1514 MagickLog10(density_xy[x].direction[i].red);
1515 channel_features[GreenPixelChannel].difference_entropy[i]-=
1516 density_xy[x].direction[i].green*
1517 MagickLog10(density_xy[x].direction[i].green);
1518 channel_features[BluePixelChannel].difference_entropy[i]-=
1519 density_xy[x].direction[i].blue*
1520 MagickLog10(density_xy[x].direction[i].blue);
1521 if (image->colorspace == CMYKColorspace)
1522 channel_features[BlackPixelChannel].difference_entropy[i]-=
1523 density_xy[x].direction[i].black*
1524 MagickLog10(density_xy[x].direction[i].black);
1525 if (image->alpha_trait == BlendPixelTrait)
1526 channel_features[AlphaPixelChannel].difference_entropy[i]-=
1527 density_xy[x].direction[i].alpha*
1528 MagickLog10(density_xy[x].direction[i].alpha);
1530 Information Measures of Correlation.
1532 entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1533 MagickLog10(density_x[x].direction[i].red));
1534 entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1535 MagickLog10(density_x[x].direction[i].green));
1536 entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1537 MagickLog10(density_x[x].direction[i].blue));
1538 if (image->colorspace == CMYKColorspace)
1539 entropy_x.direction[i].black-=(density_x[x].direction[i].black*
1540 MagickLog10(density_x[x].direction[i].black));
1541 if (image->alpha_trait == BlendPixelTrait)
1542 entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha*
1543 MagickLog10(density_x[x].direction[i].alpha));
1544 entropy_y.direction[i].red-=(density_y[x].direction[i].red*
1545 MagickLog10(density_y[x].direction[i].red));
1546 entropy_y.direction[i].green-=(density_y[x].direction[i].green*
1547 MagickLog10(density_y[x].direction[i].green));
1548 entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
1549 MagickLog10(density_y[x].direction[i].blue));
1550 if (image->colorspace == CMYKColorspace)
1551 entropy_y.direction[i].black-=(density_y[x].direction[i].black*
1552 MagickLog10(density_y[x].direction[i].black));
1553 if (image->alpha_trait == BlendPixelTrait)
1554 entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
1555 MagickLog10(density_y[x].direction[i].alpha));
1558 Difference variance.
1560 channel_features[RedPixelChannel].difference_variance[i]=
1561 (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1562 (variance.direction[i].red*variance.direction[i].red))/
1563 ((double) number_grays*number_grays*number_grays*number_grays);
1564 channel_features[GreenPixelChannel].difference_variance[i]=
1565 (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1566 (variance.direction[i].green*variance.direction[i].green))/
1567 ((double) number_grays*number_grays*number_grays*number_grays);
1568 channel_features[BluePixelChannel].difference_variance[i]=
1569 (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1570 (variance.direction[i].blue*variance.direction[i].blue))/
1571 ((double) number_grays*number_grays*number_grays*number_grays);
1572 if (image->colorspace == CMYKColorspace)
1573 channel_features[BlackPixelChannel].difference_variance[i]=
1574 (((double) number_grays*number_grays*sum_squares.direction[i].black)-
1575 (variance.direction[i].black*variance.direction[i].black))/
1576 ((double) number_grays*number_grays*number_grays*number_grays);
1577 if (image->alpha_trait == BlendPixelTrait)
1578 channel_features[AlphaPixelChannel].difference_variance[i]=
1579 (((double) number_grays*number_grays*sum_squares.direction[i].alpha)-
1580 (variance.direction[i].alpha*variance.direction[i].alpha))/
1581 ((double) number_grays*number_grays*number_grays*number_grays);
1583 Information Measures of Correlation.
1585 channel_features[RedPixelChannel].measure_of_correlation_1[i]=
1586 (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1587 (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1588 entropy_x.direction[i].red : entropy_y.direction[i].red);
1589 channel_features[GreenPixelChannel].measure_of_correlation_1[i]=
1590 (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1591 (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1592 entropy_x.direction[i].green : entropy_y.direction[i].green);
1593 channel_features[BluePixelChannel].measure_of_correlation_1[i]=
1594 (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1595 (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1596 entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1597 if (image->colorspace == CMYKColorspace)
1598 channel_features[BlackPixelChannel].measure_of_correlation_1[i]=
1599 (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/
1600 (entropy_x.direction[i].black > entropy_y.direction[i].black ?
1601 entropy_x.direction[i].black : entropy_y.direction[i].black);
1602 if (image->alpha_trait == BlendPixelTrait)
1603 channel_features[AlphaPixelChannel].measure_of_correlation_1[i]=
1604 (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/
1605 (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ?
1606 entropy_x.direction[i].alpha : entropy_y.direction[i].alpha);
1607 channel_features[RedPixelChannel].measure_of_correlation_2[i]=
1608 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].red-
1609 entropy_xy.direction[i].red)))));
1610 channel_features[GreenPixelChannel].measure_of_correlation_2[i]=
1611 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].green-
1612 entropy_xy.direction[i].green)))));
1613 channel_features[BluePixelChannel].measure_of_correlation_2[i]=
1614 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].blue-
1615 entropy_xy.direction[i].blue)))));
1616 if (image->colorspace == CMYKColorspace)
1617 channel_features[BlackPixelChannel].measure_of_correlation_2[i]=
1618 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].black-
1619 entropy_xy.direction[i].black)))));
1620 if (image->alpha_trait == BlendPixelTrait)
1621 channel_features[AlphaPixelChannel].measure_of_correlation_2[i]=
1622 (sqrt(fabs(1.0-exp(-2.0*(entropy_xy2.direction[i].alpha-
1623 entropy_xy.direction[i].alpha)))));
1626 Compute more texture features.
1628 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1629 #pragma omp parallel for schedule(static,4) shared(status) \
1630 magick_threads(image,image,number_grays,1)
1632 for (i=0; i < 4; i++)
1637 for (z=0; z < (ssize_t) number_grays; z++)
1645 (void) ResetMagickMemory(&pixel,0,sizeof(pixel));
1646 for (y=0; y < (ssize_t) number_grays; y++)
1651 for (x=0; x < (ssize_t) number_grays; x++)
1654 Contrast: amount of local variations present in an image.
1656 if (((y-x) == z) || ((x-y) == z))
1658 pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1659 pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1660 pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1661 if (image->colorspace == CMYKColorspace)
1662 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
1663 if (image->alpha_trait == BlendPixelTrait)
1664 pixel.direction[i].alpha+=
1665 cooccurrence[x][y].direction[i].alpha;
1668 Maximum Correlation Coefficient.
1670 Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1671 cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1672 density_y[x].direction[i].red;
1673 Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1674 cooccurrence[y][x].direction[i].green/
1675 density_x[z].direction[i].green/density_y[x].direction[i].red;
1676 Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1677 cooccurrence[y][x].direction[i].blue/density_x[z].direction[i].blue/
1678 density_y[x].direction[i].blue;
1679 if (image->colorspace == CMYKColorspace)
1680 Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
1681 cooccurrence[y][x].direction[i].black/
1682 density_x[z].direction[i].black/density_y[x].direction[i].black;
1683 if (image->alpha_trait == BlendPixelTrait)
1684 Q[z][y].direction[i].alpha+=
1685 cooccurrence[z][x].direction[i].alpha*
1686 cooccurrence[y][x].direction[i].alpha/
1687 density_x[z].direction[i].alpha/
1688 density_y[x].direction[i].alpha;
1691 channel_features[RedPixelChannel].contrast[i]+=z*z*
1692 pixel.direction[i].red;
1693 channel_features[GreenPixelChannel].contrast[i]+=z*z*
1694 pixel.direction[i].green;
1695 channel_features[BluePixelChannel].contrast[i]+=z*z*
1696 pixel.direction[i].blue;
1697 if (image->colorspace == CMYKColorspace)
1698 channel_features[BlackPixelChannel].contrast[i]+=z*z*
1699 pixel.direction[i].black;
1700 if (image->alpha_trait == BlendPixelTrait)
1701 channel_features[AlphaPixelChannel].contrast[i]+=z*z*
1702 pixel.direction[i].alpha;
1705 Maximum Correlation Coefficient.
1706 Future: return second largest eigenvalue of Q.
1708 channel_features[RedPixelChannel].maximum_correlation_coefficient[i]=
1709 sqrt((double) -1.0);
1710 channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]=
1711 sqrt((double) -1.0);
1712 channel_features[BluePixelChannel].maximum_correlation_coefficient[i]=
1713 sqrt((double) -1.0);
1714 if (image->colorspace == CMYKColorspace)
1715 channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]=
1716 sqrt((double) -1.0);
1717 if (image->alpha_trait == BlendPixelTrait)
1718 channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]=
1719 sqrt((double) -1.0);
1722 Relinquish resources.
1724 sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1725 for (i=0; i < (ssize_t) number_grays; i++)
1726 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1727 Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1728 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1729 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1730 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1731 for (i=0; i < (ssize_t) number_grays; i++)
1732 cooccurrence[i]=(ChannelStatistics *)
1733 RelinquishMagickMemory(cooccurrence[i]);
1734 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1735 return(channel_features);