43#include "MagickCore/studio.h"
44#include "MagickCore/animate.h"
45#include "MagickCore/artifact.h"
46#include "MagickCore/blob.h"
47#include "MagickCore/blob-private.h"
48#include "MagickCore/cache.h"
49#include "MagickCore/cache-private.h"
50#include "MagickCore/cache-view.h"
51#include "MagickCore/channel.h"
52#include "MagickCore/client.h"
53#include "MagickCore/color.h"
54#include "MagickCore/color-private.h"
55#include "MagickCore/colorspace.h"
56#include "MagickCore/colorspace-private.h"
57#include "MagickCore/composite.h"
58#include "MagickCore/composite-private.h"
59#include "MagickCore/compress.h"
60#include "MagickCore/constitute.h"
61#include "MagickCore/display.h"
62#include "MagickCore/draw.h"
63#include "MagickCore/enhance.h"
64#include "MagickCore/exception.h"
65#include "MagickCore/exception-private.h"
66#include "MagickCore/feature.h"
67#include "MagickCore/gem.h"
68#include "MagickCore/geometry.h"
69#include "MagickCore/list.h"
70#include "MagickCore/image-private.h"
71#include "MagickCore/magic.h"
72#include "MagickCore/magick.h"
73#include "MagickCore/matrix.h"
74#include "MagickCore/memory_.h"
75#include "MagickCore/module.h"
76#include "MagickCore/monitor.h"
77#include "MagickCore/monitor-private.h"
78#include "MagickCore/morphology-private.h"
79#include "MagickCore/nt-base-private.h"
80#include "MagickCore/option.h"
81#include "MagickCore/paint.h"
82#include "MagickCore/pixel-accessor.h"
83#include "MagickCore/profile.h"
84#include "MagickCore/property.h"
85#include "MagickCore/quantize.h"
86#include "MagickCore/quantum-private.h"
87#include "MagickCore/random_.h"
88#include "MagickCore/resource_.h"
89#include "MagickCore/segment.h"
90#include "MagickCore/semaphore.h"
91#include "MagickCore/signature-private.h"
92#include "MagickCore/statistic-private.h"
93#include "MagickCore/string_.h"
94#include "MagickCore/thread-private.h"
95#include "MagickCore/timer.h"
96#include "MagickCore/utility.h"
97#include "MagickCore/utility-private.h"
98#include "MagickCore/version.h"
150static inline MagickBooleanType IsAuthenticPixel(
const Image *image,
151 const ssize_t x,
const ssize_t y)
153 if ((x < 0) || (x >= (ssize_t) image->columns))
155 if ((y < 0) || (y >= (ssize_t) image->rows))
160static MagickBooleanType TraceEdges(
Image *edge_image,
CacheView *edge_view,
161 MatrixInfo *canny_cache,
const ssize_t x,
const ssize_t y,
177 q=GetCacheViewAuthenticPixels(edge_view,x,y,1,1,exception);
178 if (q == (Quantum *) NULL)
181 status=SyncCacheViewAuthenticPixels(edge_view,exception);
182 if (status == MagickFalse)
184 if (GetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
188 if (SetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
196 status=GetMatrixElement(canny_cache,i,0,&edge);
197 if (status == MagickFalse)
199 for (v=(-1); v <= 1; v++)
204 for (u=(-1); u <= 1; u++)
206 if ((u == 0) && (v == 0))
208 if (IsAuthenticPixel(edge_image,edge.x+u,edge.y+v) == MagickFalse)
213 q=GetCacheViewAuthenticPixels(edge_view,edge.x+u,edge.y+v,1,1,
215 if (q == (Quantum *) NULL)
217 status=GetMatrixElement(canny_cache,edge.x+u,edge.y+v,&pixel);
218 if (status == MagickFalse)
220 if ((GetPixelIntensity(edge_image,q) == 0.0) &&
221 (pixel.intensity >= lower_threshold))
224 status=SyncCacheViewAuthenticPixels(edge_view,exception);
225 if (status == MagickFalse)
229 status=SetMatrixElement(canny_cache,i,0,&edge);
230 if (status == MagickFalse)
240MagickExport
Image *CannyEdgeImage(
const Image *image,
const double radius,
241 const double sigma,
const double lower_percent,
const double upper_percent,
244#define CannyEdgeImageTag "CannyEdge/Image"
253 geometry[MagickPathExtent];
279 assert(image != (
const Image *) NULL);
280 assert(image->signature == MagickCoreSignature);
282 assert(exception->signature == MagickCoreSignature);
283 if (IsEventLogging() != MagickFalse)
284 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
288 (void) FormatLocaleString(geometry,MagickPathExtent,
289 "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
290 kernel_info=AcquireKernelInfo(geometry,exception);
292 ThrowImageException(ResourceLimitError,
"MemoryAllocationFailed");
293 edge_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,exception);
294 kernel_info=DestroyKernelInfo(kernel_info);
295 if (edge_image == (
Image *) NULL)
296 return((
Image *) NULL);
297 if (TransformImageColorspace(edge_image,GRAYColorspace,exception) == MagickFalse)
299 edge_image=DestroyImage(edge_image);
300 return((
Image *) NULL);
302 (void) SetImageAlphaChannel(edge_image,OffAlphaChannel,exception);
306 canny_cache=AcquireMatrixInfo(edge_image->columns,edge_image->rows,
310 edge_image=DestroyImage(edge_image);
311 return((
Image *) NULL);
314 edge_view=AcquireVirtualCacheView(edge_image,exception);
315#if defined(MAGICKCORE_OPENMP_SUPPORT)
316 #pragma omp parallel for schedule(static) shared(status) \
317 magick_number_threads(edge_image,edge_image,edge_image->rows,1)
319 for (y=0; y < (ssize_t) edge_image->rows; y++)
327 if (status == MagickFalse)
329 p=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns+1,2,
331 if (p == (
const Quantum *) NULL)
336 for (x=0; x < (ssize_t) edge_image->columns; x++)
346 *magick_restrict kernel_pixels;
363 (void) memset(&pixel,0,
sizeof(pixel));
367 for (v=0; v < 2; v++)
372 for (u=0; u < 2; u++)
377 intensity=GetPixelIntensity(edge_image,kernel_pixels+u);
378 dx+=0.5*Gx[v][u]*intensity;
379 dy+=0.5*Gy[v][u]*intensity;
381 kernel_pixels+=edge_image->columns+1;
383 pixel.magnitude=hypot(dx,dy);
385 if (fabs(dx) > MagickEpsilon)
393 if (slope < -2.41421356237)
396 if (slope < -0.414213562373)
403 if (slope > 2.41421356237)
406 if (slope > 0.414213562373)
412 if (SetMatrixElement(canny_cache,x,y,&pixel) == MagickFalse)
414 p+=(ptrdiff_t) GetPixelChannels(edge_image);
417 edge_view=DestroyCacheView(edge_view);
423 (void) GetMatrixElement(canny_cache,0,0,&element);
424 max=element.intensity;
425 min=element.intensity;
426 edge_view=AcquireAuthenticCacheView(edge_image,exception);
427#if defined(MAGICKCORE_OPENMP_SUPPORT)
428 #pragma omp parallel for schedule(static) shared(status) \
429 magick_number_threads(edge_image,edge_image,edge_image->rows,1)
431 for (y=0; y < (ssize_t) edge_image->rows; y++)
439 if (status == MagickFalse)
441 q=GetCacheViewAuthenticPixels(edge_view,0,y,edge_image->columns,1,
443 if (q == (Quantum *) NULL)
448 for (x=0; x < (ssize_t) edge_image->columns; x++)
455 (void) GetMatrixElement(canny_cache,x,y,&pixel);
456 switch (pixel.orientation)
464 (void) GetMatrixElement(canny_cache,x,y-1,&alpha_pixel);
465 (void) GetMatrixElement(canny_cache,x,y+1,&beta_pixel);
473 (void) GetMatrixElement(canny_cache,x-1,y-1,&alpha_pixel);
474 (void) GetMatrixElement(canny_cache,x+1,y+1,&beta_pixel);
482 (void) GetMatrixElement(canny_cache,x-1,y,&alpha_pixel);
483 (void) GetMatrixElement(canny_cache,x+1,y,&beta_pixel);
491 (void) GetMatrixElement(canny_cache,x+1,y-1,&beta_pixel);
492 (void) GetMatrixElement(canny_cache,x-1,y+1,&alpha_pixel);
496 pixel.intensity=pixel.magnitude;
497 if ((pixel.magnitude < alpha_pixel.magnitude) ||
498 (pixel.magnitude < beta_pixel.magnitude))
500 (void) SetMatrixElement(canny_cache,x,y,&pixel);
501#if defined(MAGICKCORE_OPENMP_SUPPORT)
502 #pragma omp critical (MagickCore_CannyEdgeImage)
505 if (pixel.intensity < min)
507 if (pixel.intensity > max)
511 q+=(ptrdiff_t) GetPixelChannels(edge_image);
513 if (SyncCacheViewAuthenticPixels(edge_view,exception) == MagickFalse)
516 edge_view=DestroyCacheView(edge_view);
520 lower_threshold=lower_percent*(max-min)+min;
521 upper_threshold=upper_percent*(max-min)+min;
525 edge_view=AcquireAuthenticCacheView(edge_image,exception);
526 for (y=0; y < (ssize_t) edge_image->rows; y++)
531 if (status == MagickFalse)
533 for (x=0; x < (ssize_t) edge_image->columns; x++)
544 p=GetCacheViewVirtualPixels(edge_view,x,y,1,1,exception);
545 if (p == (
const Quantum *) NULL)
547 status=GetMatrixElement(canny_cache,x,y,&pixel);
548 if (status == MagickFalse)
550 if ((GetPixelIntensity(edge_image,p) == 0.0) &&
551 (pixel.intensity >= upper_threshold))
552 status=TraceEdges(edge_image,edge_view,canny_cache,x,y,lower_threshold,
555 if (image->progress_monitor != (MagickProgressMonitor) NULL)
560#if defined(MAGICKCORE_OPENMP_SUPPORT)
564 proceed=SetImageProgress(image,CannyEdgeImageTag,progress,image->rows);
565 if (proceed == MagickFalse)
569 edge_view=DestroyCacheView(edge_view);
573 canny_cache=DestroyMatrixInfo(canny_cache);
665 assert(image != (
Image *) NULL);
666 assert(image->signature == MagickCoreSignature);
667 if (IsEventLogging() != MagickFalse)
668 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
669 if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
671 length=MaxPixelChannels+1UL;
673 sizeof(*channel_features));
675 ThrowFatalException(ResourceLimitFatalError,
"MemoryAllocationFailed");
676 (void) memset(channel_features,0,length*
677 sizeof(*channel_features));
681 grays=(
PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,
sizeof(*grays));
686 (void) ThrowMagickException(exception,GetMagickModule(),
687 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
688 return(channel_features);
690 for (i=0; i <= (ssize_t) MaxMap; i++)
693 grays[i].green=(~0U);
695 grays[i].alpha=(~0U);
696 grays[i].black=(~0U);
699 image_view=AcquireVirtualCacheView(image,exception);
700#if defined(MAGICKCORE_OPENMP_SUPPORT)
701 #pragma omp parallel for schedule(static) shared(status) \
702 magick_number_threads(image,image,image->rows,1)
704 for (r=0; r < (ssize_t) image->rows; r++)
712 if (status == MagickFalse)
714 p=GetCacheViewVirtualPixels(image_view,0,r,image->columns,1,exception);
715 if (p == (
const Quantum *) NULL)
720 for (x=0; x < (ssize_t) image->columns; x++)
722 grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
723 ScaleQuantumToMap(GetPixelRed(image,p));
724 grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
725 ScaleQuantumToMap(GetPixelGreen(image,p));
726 grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
727 ScaleQuantumToMap(GetPixelBlue(image,p));
728 if (image->colorspace == CMYKColorspace)
729 grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
730 ScaleQuantumToMap(GetPixelBlack(image,p));
731 if (image->alpha_trait != UndefinedPixelTrait)
732 grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
733 ScaleQuantumToMap(GetPixelAlpha(image,p));
734 p+=(ptrdiff_t) GetPixelChannels(image);
737 image_view=DestroyCacheView(image_view);
738 if (status == MagickFalse)
740 grays=(
PixelPacket *) RelinquishMagickMemory(grays);
743 return(channel_features);
745 (void) memset(&gray,0,
sizeof(gray));
746 for (i=0; i <= (ssize_t) MaxMap; i++)
748 if (grays[i].red != ~0U)
749 grays[gray.red++].red=grays[i].red;
750 if (grays[i].green != ~0U)
751 grays[gray.green++].green=grays[i].green;
752 if (grays[i].blue != ~0U)
753 grays[gray.blue++].blue=grays[i].blue;
754 if (image->colorspace == CMYKColorspace)
755 if (grays[i].black != ~0U)
756 grays[gray.black++].black=grays[i].black;
757 if (image->alpha_trait != UndefinedPixelTrait)
758 if (grays[i].alpha != ~0U)
759 grays[gray.alpha++].alpha=grays[i].alpha;
764 number_grays=gray.red;
765 if (gray.green > number_grays)
766 number_grays=gray.green;
767 if (gray.blue > number_grays)
768 number_grays=gray.blue;
769 if (image->colorspace == CMYKColorspace)
770 if (gray.black > number_grays)
771 number_grays=gray.black;
772 if (image->alpha_trait != UndefinedPixelTrait)
773 if (gray.alpha > number_grays)
774 number_grays=gray.alpha;
776 sizeof(*cooccurrence));
778 2*
sizeof(*density_x));
780 2*
sizeof(*density_xy));
782 2*
sizeof(*density_y));
794 for (i=0; i < (ssize_t) number_grays; i++)
808 for (i=0; i < (ssize_t) number_grays; i++)
810 RelinquishMagickMemory(cooccurrence[i]);
814 grays=(
PixelPacket *) RelinquishMagickMemory(grays);
817 (void) ThrowMagickException(exception,GetMagickModule(),
818 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
819 return(channel_features);
821 (void) memset(&correlation,0,
sizeof(correlation));
822 (void) memset(density_x,0,2*(number_grays+1)*
sizeof(*density_x));
823 (void) memset(density_xy,0,2*(number_grays+1)*
sizeof(*density_xy));
824 (void) memset(density_y,0,2*(number_grays+1)*
sizeof(*density_y));
825 (void) memset(&mean,0,
sizeof(mean));
826 (void) memset(sum,0,number_grays*
sizeof(*sum));
827 (void) memset(&sum_squares,0,
sizeof(sum_squares));
828 (void) memset(density_xy,0,2*number_grays*
sizeof(*density_xy));
829 (void) memset(&entropy_x,0,
sizeof(entropy_x));
830 (void) memset(&entropy_xy,0,
sizeof(entropy_xy));
831 (void) memset(&entropy_xy1,0,
sizeof(entropy_xy1));
832 (void) memset(&entropy_xy2,0,
sizeof(entropy_xy2));
833 (void) memset(&entropy_y,0,
sizeof(entropy_y));
834 (void) memset(&variance,0,
sizeof(variance));
835 for (i=0; i < (ssize_t) number_grays; i++)
838 sizeof(**cooccurrence));
843 (void) memset(cooccurrence[i],0,number_grays*
844 sizeof(**cooccurrence));
845 (void) memset(Q[i],0,number_grays*
sizeof(**Q));
847 if (i < (ssize_t) number_grays)
849 for (i--; i >= 0; i--)
855 RelinquishMagickMemory(cooccurrence[i]);
863 grays=(
PixelPacket *) RelinquishMagickMemory(grays);
866 (void) ThrowMagickException(exception,GetMagickModule(),
867 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
868 return(channel_features);
874 image_view=AcquireVirtualCacheView(image,exception);
875 for (r=0; r < (ssize_t) image->rows; r++)
888 if (status == MagickFalse)
890 p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,r,image->columns+
891 2*distance,distance+2,exception);
892 if (p == (
const Quantum *) NULL)
897 p+=(ptrdiff_t) distance*GetPixelChannels(image);;
898 for (x=0; x < (ssize_t) image->columns; x++)
900 for (i=0; i < 4; i++)
910 offset=(ssize_t) distance;
918 offset=(ssize_t) (image->columns+2*distance);
926 offset=(ssize_t) ((image->columns+2*distance)-distance);
934 offset=(ssize_t) ((image->columns+2*distance)+distance);
940 while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p)))
942 while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*(ssize_t) GetPixelChannels(image))))
944 cooccurrence[u][v].direction[i].red++;
945 cooccurrence[v][u].direction[i].red++;
948 while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p)))
950 while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*(ssize_t) GetPixelChannels(image))))
952 cooccurrence[u][v].direction[i].green++;
953 cooccurrence[v][u].direction[i].green++;
956 while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p)))
958 while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*(ssize_t) GetPixelChannels(image))))
960 cooccurrence[u][v].direction[i].blue++;
961 cooccurrence[v][u].direction[i].blue++;
962 if (image->colorspace == CMYKColorspace)
966 while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p)))
968 while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*(ssize_t) GetPixelChannels(image))))
970 cooccurrence[u][v].direction[i].black++;
971 cooccurrence[v][u].direction[i].black++;
973 if (image->alpha_trait != UndefinedPixelTrait)
977 while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p)))
979 while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*(ssize_t) GetPixelChannels(image))))
981 cooccurrence[u][v].direction[i].alpha++;
982 cooccurrence[v][u].direction[i].alpha++;
985 p+=(ptrdiff_t) GetPixelChannels(image);
988 grays=(
PixelPacket *) RelinquishMagickMemory(grays);
989 image_view=DestroyCacheView(image_view);
990 if (status == MagickFalse)
992 for (i=0; i < (ssize_t) number_grays; i++)
994 RelinquishMagickMemory(cooccurrence[i]);
998 (void) ThrowMagickException(exception,GetMagickModule(),
999 ResourceLimitError,
"MemoryAllocationFailed",
"`%s'",image->filename);
1000 return(channel_features);
1005 for (i=0; i < 4; i++)
1021 normalize=2.0*image->rows*(image->columns-distance);
1029 normalize=2.0*(image->rows-distance)*image->columns;
1037 normalize=2.0*(image->rows-distance)*(image->columns-distance);
1045 normalize=2.0*(image->rows-distance)*(image->columns-distance);
1049 normalize=MagickSafeReciprocal(normalize);
1050 for (y=0; y < (ssize_t) number_grays; y++)
1055 for (x=0; x < (ssize_t) number_grays; x++)
1057 cooccurrence[x][y].direction[i].red*=normalize;
1058 cooccurrence[x][y].direction[i].green*=normalize;
1059 cooccurrence[x][y].direction[i].blue*=normalize;
1060 if (image->colorspace == CMYKColorspace)
1061 cooccurrence[x][y].direction[i].black*=normalize;
1062 if (image->alpha_trait != UndefinedPixelTrait)
1063 cooccurrence[x][y].direction[i].alpha*=normalize;
1070#if defined(MAGICKCORE_OPENMP_SUPPORT)
1071 #pragma omp parallel for schedule(static) shared(status) \
1072 magick_number_threads(image,image,number_grays,1)
1074 for (i=0; i < 4; i++)
1079 for (y=0; y < (ssize_t) number_grays; y++)
1084 for (x=0; x < (ssize_t) number_grays; x++)
1089 channel_features[RedPixelChannel].angular_second_moment[i]+=
1090 cooccurrence[x][y].direction[i].red*
1091 cooccurrence[x][y].direction[i].red;
1092 channel_features[GreenPixelChannel].angular_second_moment[i]+=
1093 cooccurrence[x][y].direction[i].green*
1094 cooccurrence[x][y].direction[i].green;
1095 channel_features[BluePixelChannel].angular_second_moment[i]+=
1096 cooccurrence[x][y].direction[i].blue*
1097 cooccurrence[x][y].direction[i].blue;
1098 if (image->colorspace == CMYKColorspace)
1099 channel_features[BlackPixelChannel].angular_second_moment[i]+=
1100 cooccurrence[x][y].direction[i].black*
1101 cooccurrence[x][y].direction[i].black;
1102 if (image->alpha_trait != UndefinedPixelTrait)
1103 channel_features[AlphaPixelChannel].angular_second_moment[i]+=
1104 cooccurrence[x][y].direction[i].alpha*
1105 cooccurrence[x][y].direction[i].alpha;
1109 sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1110 sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1111 sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1112 if (image->colorspace == CMYKColorspace)
1113 sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black;
1114 if (image->alpha_trait != UndefinedPixelTrait)
1115 sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha;
1116 correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
1117 correlation.direction[i].green+=x*y*
1118 cooccurrence[x][y].direction[i].green;
1119 correlation.direction[i].blue+=x*y*
1120 cooccurrence[x][y].direction[i].blue;
1121 if (image->colorspace == CMYKColorspace)
1122 correlation.direction[i].black+=x*y*
1123 cooccurrence[x][y].direction[i].black;
1124 if (image->alpha_trait != UndefinedPixelTrait)
1125 correlation.direction[i].alpha+=x*y*
1126 cooccurrence[x][y].direction[i].alpha;
1130 channel_features[RedPixelChannel].inverse_difference_moment[i]+=
1131 cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
1132 channel_features[GreenPixelChannel].inverse_difference_moment[i]+=
1133 cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
1134 channel_features[BluePixelChannel].inverse_difference_moment[i]+=
1135 cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
1136 if (image->colorspace == CMYKColorspace)
1137 channel_features[BlackPixelChannel].inverse_difference_moment[i]+=
1138 cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1);
1139 if (image->alpha_trait != UndefinedPixelTrait)
1140 channel_features[AlphaPixelChannel].inverse_difference_moment[i]+=
1141 cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
1145 density_xy[y+x+2].direction[i].red+=
1146 cooccurrence[x][y].direction[i].red;
1147 density_xy[y+x+2].direction[i].green+=
1148 cooccurrence[x][y].direction[i].green;
1149 density_xy[y+x+2].direction[i].blue+=
1150 cooccurrence[x][y].direction[i].blue;
1151 if (image->colorspace == CMYKColorspace)
1152 density_xy[y+x+2].direction[i].black+=
1153 cooccurrence[x][y].direction[i].black;
1154 if (image->alpha_trait != UndefinedPixelTrait)
1155 density_xy[y+x+2].direction[i].alpha+=
1156 cooccurrence[x][y].direction[i].alpha;
1160 channel_features[RedPixelChannel].entropy[i]-=
1161 cooccurrence[x][y].direction[i].red*
1162 log2(cooccurrence[x][y].direction[i].red);
1163 channel_features[GreenPixelChannel].entropy[i]-=
1164 cooccurrence[x][y].direction[i].green*
1165 log2(cooccurrence[x][y].direction[i].green);
1166 channel_features[BluePixelChannel].entropy[i]-=
1167 cooccurrence[x][y].direction[i].blue*
1168 log2(cooccurrence[x][y].direction[i].blue);
1169 if (image->colorspace == CMYKColorspace)
1170 channel_features[BlackPixelChannel].entropy[i]-=
1171 cooccurrence[x][y].direction[i].black*
1172 log2(cooccurrence[x][y].direction[i].black);
1173 if (image->alpha_trait != UndefinedPixelTrait)
1174 channel_features[AlphaPixelChannel].entropy[i]-=
1175 cooccurrence[x][y].direction[i].alpha*
1176 log2(cooccurrence[x][y].direction[i].alpha);
1180 density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
1181 density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
1182 density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1183 if (image->alpha_trait != UndefinedPixelTrait)
1184 density_x[x].direction[i].alpha+=
1185 cooccurrence[x][y].direction[i].alpha;
1186 if (image->colorspace == CMYKColorspace)
1187 density_x[x].direction[i].black+=
1188 cooccurrence[x][y].direction[i].black;
1189 density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1190 density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1191 density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1192 if (image->colorspace == CMYKColorspace)
1193 density_y[y].direction[i].black+=
1194 cooccurrence[x][y].direction[i].black;
1195 if (image->alpha_trait != UndefinedPixelTrait)
1196 density_y[y].direction[i].alpha+=
1197 cooccurrence[x][y].direction[i].alpha;
1199 mean.direction[i].red+=y*sum[y].direction[i].red;
1200 sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
1201 mean.direction[i].green+=y*sum[y].direction[i].green;
1202 sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
1203 mean.direction[i].blue+=y*sum[y].direction[i].blue;
1204 sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
1205 if (image->colorspace == CMYKColorspace)
1207 mean.direction[i].black+=y*sum[y].direction[i].black;
1208 sum_squares.direction[i].black+=y*y*sum[y].direction[i].black;
1210 if (image->alpha_trait != UndefinedPixelTrait)
1212 mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
1213 sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
1219 channel_features[RedPixelChannel].correlation[i]=
1220 (correlation.direction[i].red-mean.direction[i].red*
1221 mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
1222 (mean.direction[i].red*mean.direction[i].red))*sqrt(
1223 sum_squares.direction[i].red-(mean.direction[i].red*
1224 mean.direction[i].red)));
1225 channel_features[GreenPixelChannel].correlation[i]=
1226 (correlation.direction[i].green-mean.direction[i].green*
1227 mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
1228 (mean.direction[i].green*mean.direction[i].green))*sqrt(
1229 sum_squares.direction[i].green-(mean.direction[i].green*
1230 mean.direction[i].green)));
1231 channel_features[BluePixelChannel].correlation[i]=
1232 (correlation.direction[i].blue-mean.direction[i].blue*
1233 mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
1234 (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
1235 sum_squares.direction[i].blue-(mean.direction[i].blue*
1236 mean.direction[i].blue)));
1237 if (image->colorspace == CMYKColorspace)
1238 channel_features[BlackPixelChannel].correlation[i]=
1239 (correlation.direction[i].black-mean.direction[i].black*
1240 mean.direction[i].black)/(sqrt(sum_squares.direction[i].black-
1241 (mean.direction[i].black*mean.direction[i].black))*sqrt(
1242 sum_squares.direction[i].black-(mean.direction[i].black*
1243 mean.direction[i].black)));
1244 if (image->alpha_trait != UndefinedPixelTrait)
1245 channel_features[AlphaPixelChannel].correlation[i]=
1246 (correlation.direction[i].alpha-mean.direction[i].alpha*
1247 mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
1248 (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt(
1249 sum_squares.direction[i].alpha-(mean.direction[i].alpha*
1250 mean.direction[i].alpha)));
1255#if defined(MAGICKCORE_OPENMP_SUPPORT)
1256 #pragma omp parallel for schedule(static) shared(status) \
1257 magick_number_threads(image,image,number_grays,1)
1259 for (i=0; i < 4; i++)
1264 for (x=2; x < (ssize_t) (2*number_grays); x++)
1269 channel_features[RedPixelChannel].sum_average[i]+=
1270 x*density_xy[x].direction[i].red;
1271 channel_features[GreenPixelChannel].sum_average[i]+=
1272 x*density_xy[x].direction[i].green;
1273 channel_features[BluePixelChannel].sum_average[i]+=
1274 x*density_xy[x].direction[i].blue;
1275 if (image->colorspace == CMYKColorspace)
1276 channel_features[BlackPixelChannel].sum_average[i]+=
1277 x*density_xy[x].direction[i].black;
1278 if (image->alpha_trait != UndefinedPixelTrait)
1279 channel_features[AlphaPixelChannel].sum_average[i]+=
1280 x*density_xy[x].direction[i].alpha;
1284 channel_features[RedPixelChannel].sum_entropy[i]-=
1285 density_xy[x].direction[i].red*
1286 log2(density_xy[x].direction[i].red);
1287 channel_features[GreenPixelChannel].sum_entropy[i]-=
1288 density_xy[x].direction[i].green*
1289 log2(density_xy[x].direction[i].green);
1290 channel_features[BluePixelChannel].sum_entropy[i]-=
1291 density_xy[x].direction[i].blue*
1292 log2(density_xy[x].direction[i].blue);
1293 if (image->colorspace == CMYKColorspace)
1294 channel_features[BlackPixelChannel].sum_entropy[i]-=
1295 density_xy[x].direction[i].black*
1296 log2(density_xy[x].direction[i].black);
1297 if (image->alpha_trait != UndefinedPixelTrait)
1298 channel_features[AlphaPixelChannel].sum_entropy[i]-=
1299 density_xy[x].direction[i].alpha*
1300 log2(density_xy[x].direction[i].alpha);
1304 channel_features[RedPixelChannel].sum_variance[i]+=
1305 (x-channel_features[RedPixelChannel].sum_entropy[i])*
1306 (x-channel_features[RedPixelChannel].sum_entropy[i])*
1307 density_xy[x].direction[i].red;
1308 channel_features[GreenPixelChannel].sum_variance[i]+=
1309 (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1310 (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1311 density_xy[x].direction[i].green;
1312 channel_features[BluePixelChannel].sum_variance[i]+=
1313 (x-channel_features[BluePixelChannel].sum_entropy[i])*
1314 (x-channel_features[BluePixelChannel].sum_entropy[i])*
1315 density_xy[x].direction[i].blue;
1316 if (image->colorspace == CMYKColorspace)
1317 channel_features[BlackPixelChannel].sum_variance[i]+=
1318 (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1319 (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1320 density_xy[x].direction[i].black;
1321 if (image->alpha_trait != UndefinedPixelTrait)
1322 channel_features[AlphaPixelChannel].sum_variance[i]+=
1323 (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1324 (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1325 density_xy[x].direction[i].alpha;
1331#if defined(MAGICKCORE_OPENMP_SUPPORT)
1332 #pragma omp parallel for schedule(static) shared(status) \
1333 magick_number_threads(image,image,number_grays,1)
1335 for (i=0; i < 4; i++)
1340 for (y=0; y < (ssize_t) number_grays; y++)
1345 for (x=0; x < (ssize_t) number_grays; x++)
1350 variance.direction[i].red+=(y-mean.direction[i].red+1)*
1351 (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
1352 variance.direction[i].green+=(y-mean.direction[i].green+1)*
1353 (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
1354 variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
1355 (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
1356 if (image->colorspace == CMYKColorspace)
1357 variance.direction[i].black+=(y-mean.direction[i].black+1)*
1358 (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black;
1359 if (image->alpha_trait != UndefinedPixelTrait)
1360 variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
1361 (y-mean.direction[i].alpha+1)*
1362 cooccurrence[x][y].direction[i].alpha;
1366 density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
1367 cooccurrence[x][y].direction[i].red;
1368 density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
1369 cooccurrence[x][y].direction[i].green;
1370 density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
1371 cooccurrence[x][y].direction[i].blue;
1372 if (image->colorspace == CMYKColorspace)
1373 density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
1374 cooccurrence[x][y].direction[i].black;
1375 if (image->alpha_trait != UndefinedPixelTrait)
1376 density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
1377 cooccurrence[x][y].direction[i].alpha;
1381 entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
1382 log2(cooccurrence[x][y].direction[i].red);
1383 entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
1384 log2(cooccurrence[x][y].direction[i].green);
1385 entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
1386 log2(cooccurrence[x][y].direction[i].blue);
1387 if (image->colorspace == CMYKColorspace)
1388 entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black*
1389 log2(cooccurrence[x][y].direction[i].black);
1390 if (image->alpha_trait != UndefinedPixelTrait)
1391 entropy_xy.direction[i].alpha-=
1392 cooccurrence[x][y].direction[i].alpha*log2(
1393 cooccurrence[x][y].direction[i].alpha);
1394 entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
1395 log2(density_x[x].direction[i].red*density_y[y].direction[i].red));
1396 entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
1397 log2(density_x[x].direction[i].green*
1398 density_y[y].direction[i].green));
1399 entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
1400 log2(density_x[x].direction[i].blue*density_y[y].direction[i].blue));
1401 if (image->colorspace == CMYKColorspace)
1402 entropy_xy1.direction[i].black-=(
1403 cooccurrence[x][y].direction[i].black*log2(
1404 density_x[x].direction[i].black*density_y[y].direction[i].black));
1405 if (image->alpha_trait != UndefinedPixelTrait)
1406 entropy_xy1.direction[i].alpha-=(
1407 cooccurrence[x][y].direction[i].alpha*log2(
1408 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1409 entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
1410 density_y[y].direction[i].red*log2(density_x[x].direction[i].red*
1411 density_y[y].direction[i].red));
1412 entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
1413 density_y[y].direction[i].green*log2(density_x[x].direction[i].green*
1414 density_y[y].direction[i].green));
1415 entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
1416 density_y[y].direction[i].blue*log2(density_x[x].direction[i].blue*
1417 density_y[y].direction[i].blue));
1418 if (image->colorspace == CMYKColorspace)
1419 entropy_xy2.direction[i].black-=(density_x[x].direction[i].black*
1420 density_y[y].direction[i].black*log2(
1421 density_x[x].direction[i].black*density_y[y].direction[i].black));
1422 if (image->alpha_trait != UndefinedPixelTrait)
1423 entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha*
1424 density_y[y].direction[i].alpha*log2(
1425 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1428 channel_features[RedPixelChannel].variance_sum_of_squares[i]=
1429 variance.direction[i].red;
1430 channel_features[GreenPixelChannel].variance_sum_of_squares[i]=
1431 variance.direction[i].green;
1432 channel_features[BluePixelChannel].variance_sum_of_squares[i]=
1433 variance.direction[i].blue;
1434 if (image->colorspace == CMYKColorspace)
1435 channel_features[BlackPixelChannel].variance_sum_of_squares[i]=
1436 variance.direction[i].black;
1437 if (image->alpha_trait != UndefinedPixelTrait)
1438 channel_features[AlphaPixelChannel].variance_sum_of_squares[i]=
1439 variance.direction[i].alpha;
1444 (void) memset(&variance,0,
sizeof(variance));
1445 (void) memset(&sum_squares,0,
sizeof(sum_squares));
1446#if defined(MAGICKCORE_OPENMP_SUPPORT)
1447 #pragma omp parallel for schedule(static) shared(status) \
1448 magick_number_threads(image,image,number_grays,1)
1450 for (i=0; i < 4; i++)
1455 for (x=0; x < (ssize_t) number_grays; x++)
1460 variance.direction[i].red+=density_xy[x].direction[i].red;
1461 variance.direction[i].green+=density_xy[x].direction[i].green;
1462 variance.direction[i].blue+=density_xy[x].direction[i].blue;
1463 if (image->colorspace == CMYKColorspace)
1464 variance.direction[i].black+=density_xy[x].direction[i].black;
1465 if (image->alpha_trait != UndefinedPixelTrait)
1466 variance.direction[i].alpha+=density_xy[x].direction[i].alpha;
1467 sum_squares.direction[i].red+=density_xy[x].direction[i].red*
1468 density_xy[x].direction[i].red;
1469 sum_squares.direction[i].green+=density_xy[x].direction[i].green*
1470 density_xy[x].direction[i].green;
1471 sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
1472 density_xy[x].direction[i].blue;
1473 if (image->colorspace == CMYKColorspace)
1474 sum_squares.direction[i].black+=density_xy[x].direction[i].black*
1475 density_xy[x].direction[i].black;
1476 if (image->alpha_trait != UndefinedPixelTrait)
1477 sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
1478 density_xy[x].direction[i].alpha;
1482 channel_features[RedPixelChannel].difference_entropy[i]-=
1483 density_xy[x].direction[i].red*
1484 log2(density_xy[x].direction[i].red);
1485 channel_features[GreenPixelChannel].difference_entropy[i]-=
1486 density_xy[x].direction[i].green*
1487 log2(density_xy[x].direction[i].green);
1488 channel_features[BluePixelChannel].difference_entropy[i]-=
1489 density_xy[x].direction[i].blue*
1490 log2(density_xy[x].direction[i].blue);
1491 if (image->colorspace == CMYKColorspace)
1492 channel_features[BlackPixelChannel].difference_entropy[i]-=
1493 density_xy[x].direction[i].black*
1494 log2(density_xy[x].direction[i].black);
1495 if (image->alpha_trait != UndefinedPixelTrait)
1496 channel_features[AlphaPixelChannel].difference_entropy[i]-=
1497 density_xy[x].direction[i].alpha*
1498 log2(density_xy[x].direction[i].alpha);
1502 entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1503 log2(density_x[x].direction[i].red));
1504 entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1505 log2(density_x[x].direction[i].green));
1506 entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1507 log2(density_x[x].direction[i].blue));
1508 if (image->colorspace == CMYKColorspace)
1509 entropy_x.direction[i].black-=(density_x[x].direction[i].black*
1510 log2(density_x[x].direction[i].black));
1511 if (image->alpha_trait != UndefinedPixelTrait)
1512 entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha*
1513 log2(density_x[x].direction[i].alpha));
1514 entropy_y.direction[i].red-=(density_y[x].direction[i].red*
1515 log2(density_y[x].direction[i].red));
1516 entropy_y.direction[i].green-=(density_y[x].direction[i].green*
1517 log2(density_y[x].direction[i].green));
1518 entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
1519 log2(density_y[x].direction[i].blue));
1520 if (image->colorspace == CMYKColorspace)
1521 entropy_y.direction[i].black-=(density_y[x].direction[i].black*
1522 log2(density_y[x].direction[i].black));
1523 if (image->alpha_trait != UndefinedPixelTrait)
1524 entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
1525 log2(density_y[x].direction[i].alpha));
1530 channel_features[RedPixelChannel].difference_variance[i]=
1531 (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1532 (variance.direction[i].red*variance.direction[i].red))/
1533 ((double) number_grays*number_grays*number_grays*number_grays);
1534 channel_features[GreenPixelChannel].difference_variance[i]=
1535 (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1536 (variance.direction[i].green*variance.direction[i].green))/
1537 ((double) number_grays*number_grays*number_grays*number_grays);
1538 channel_features[BluePixelChannel].difference_variance[i]=
1539 (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1540 (variance.direction[i].blue*variance.direction[i].blue))/
1541 ((double) number_grays*number_grays*number_grays*number_grays);
1542 if (image->colorspace == CMYKColorspace)
1543 channel_features[BlackPixelChannel].difference_variance[i]=
1544 (((double) number_grays*number_grays*sum_squares.direction[i].black)-
1545 (variance.direction[i].black*variance.direction[i].black))/
1546 ((double) number_grays*number_grays*number_grays*number_grays);
1547 if (image->alpha_trait != UndefinedPixelTrait)
1548 channel_features[AlphaPixelChannel].difference_variance[i]=
1549 (((double) number_grays*number_grays*sum_squares.direction[i].alpha)-
1550 (variance.direction[i].alpha*variance.direction[i].alpha))/
1551 ((double) number_grays*number_grays*number_grays*number_grays);
1555 channel_features[RedPixelChannel].measure_of_correlation_1[i]=
1556 (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1557 (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1558 entropy_x.direction[i].red : entropy_y.direction[i].red);
1559 channel_features[GreenPixelChannel].measure_of_correlation_1[i]=
1560 (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1561 (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1562 entropy_x.direction[i].green : entropy_y.direction[i].green);
1563 channel_features[BluePixelChannel].measure_of_correlation_1[i]=
1564 (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1565 (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1566 entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1567 if (image->colorspace == CMYKColorspace)
1568 channel_features[BlackPixelChannel].measure_of_correlation_1[i]=
1569 (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/
1570 (entropy_x.direction[i].black > entropy_y.direction[i].black ?
1571 entropy_x.direction[i].black : entropy_y.direction[i].black);
1572 if (image->alpha_trait != UndefinedPixelTrait)
1573 channel_features[AlphaPixelChannel].measure_of_correlation_1[i]=
1574 (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/
1575 (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ?
1576 entropy_x.direction[i].alpha : entropy_y.direction[i].alpha);
1577 channel_features[RedPixelChannel].measure_of_correlation_2[i]=
1578 (sqrt(fabs(1.0-exp(-2.0*(
double) (entropy_xy2.direction[i].red-
1579 entropy_xy.direction[i].red)))));
1580 channel_features[GreenPixelChannel].measure_of_correlation_2[i]=
1581 (sqrt(fabs(1.0-exp(-2.0*(
double) (entropy_xy2.direction[i].green-
1582 entropy_xy.direction[i].green)))));
1583 channel_features[BluePixelChannel].measure_of_correlation_2[i]=
1584 (sqrt(fabs(1.0-exp(-2.0*(
double) (entropy_xy2.direction[i].blue-
1585 entropy_xy.direction[i].blue)))));
1586 if (image->colorspace == CMYKColorspace)
1587 channel_features[BlackPixelChannel].measure_of_correlation_2[i]=
1588 (sqrt(fabs(1.0-exp(-2.0*(
double) (entropy_xy2.direction[i].black-
1589 entropy_xy.direction[i].black)))));
1590 if (image->alpha_trait != UndefinedPixelTrait)
1591 channel_features[AlphaPixelChannel].measure_of_correlation_2[i]=
1592 (sqrt(fabs(1.0-exp(-2.0*(
double) (entropy_xy2.direction[i].alpha-
1593 entropy_xy.direction[i].alpha)))));
1598#if defined(MAGICKCORE_OPENMP_SUPPORT)
1599 #pragma omp parallel for schedule(static) shared(status) \
1600 magick_number_threads(image,image,number_grays,1)
1602 for (i=0; i < 4; i++)
1607 for (z=0; z < (ssize_t) number_grays; z++)
1615 (void) memset(&pixel,0,
sizeof(pixel));
1616 for (y=0; y < (ssize_t) number_grays; y++)
1621 for (x=0; x < (ssize_t) number_grays; x++)
1626 if (((y-x) == z) || ((x-y) == z))
1628 pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1629 pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1630 pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1631 if (image->colorspace == CMYKColorspace)
1632 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
1633 if (image->alpha_trait != UndefinedPixelTrait)
1634 pixel.direction[i].alpha+=
1635 cooccurrence[x][y].direction[i].alpha;
1640 if ((fabs(density_x[z].direction[i].red) > MagickEpsilon) &&
1641 (fabs(density_y[x].direction[i].red) > MagickEpsilon))
1642 Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1643 cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1644 density_y[x].direction[i].red;
1645 if ((fabs(density_x[z].direction[i].green) > MagickEpsilon) &&
1646 (fabs(density_y[x].direction[i].red) > MagickEpsilon))
1647 Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1648 cooccurrence[y][x].direction[i].green/
1649 density_x[z].direction[i].green/density_y[x].direction[i].red;
1650 if ((fabs(density_x[z].direction[i].blue) > MagickEpsilon) &&
1651 (fabs(density_y[x].direction[i].blue) > MagickEpsilon))
1652 Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1653 cooccurrence[y][x].direction[i].blue/
1654 density_x[z].direction[i].blue/density_y[x].direction[i].blue;
1655 if (image->colorspace == CMYKColorspace)
1656 if ((fabs(density_x[z].direction[i].black) > MagickEpsilon) &&
1657 (fabs(density_y[x].direction[i].black) > MagickEpsilon))
1658 Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
1659 cooccurrence[y][x].direction[i].black/
1660 density_x[z].direction[i].black/density_y[x].direction[i].black;
1661 if (image->alpha_trait != UndefinedPixelTrait)
1662 if ((fabs(density_x[z].direction[i].alpha) > MagickEpsilon) &&
1663 (fabs(density_y[x].direction[i].alpha) > MagickEpsilon))
1664 Q[z][y].direction[i].alpha+=
1665 cooccurrence[z][x].direction[i].alpha*
1666 cooccurrence[y][x].direction[i].alpha/
1667 density_x[z].direction[i].alpha/
1668 density_y[x].direction[i].alpha;
1671 channel_features[RedPixelChannel].contrast[i]+=z*z*
1672 pixel.direction[i].red;
1673 channel_features[GreenPixelChannel].contrast[i]+=z*z*
1674 pixel.direction[i].green;
1675 channel_features[BluePixelChannel].contrast[i]+=z*z*
1676 pixel.direction[i].blue;
1677 if (image->colorspace == CMYKColorspace)
1678 channel_features[BlackPixelChannel].contrast[i]+=z*z*
1679 pixel.direction[i].black;
1680 if (image->alpha_trait != UndefinedPixelTrait)
1681 channel_features[AlphaPixelChannel].contrast[i]+=z*z*
1682 pixel.direction[i].alpha;
1688 channel_features[RedPixelChannel].maximum_correlation_coefficient[i]=
1690 channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]=
1692 channel_features[BluePixelChannel].maximum_correlation_coefficient[i]=
1694 if (image->colorspace == CMYKColorspace)
1695 channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]=
1697 if (image->alpha_trait != UndefinedPixelTrait)
1698 channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]=
1705 for (i=0; i < (ssize_t) number_grays; i++)
1711 for (i=0; i < (ssize_t) number_grays; i++)
1713 RelinquishMagickMemory(cooccurrence[i]);
1715 return(channel_features);
1758static inline double MagickRound(
double x)
1763 if ((x-floor(x)) < (ceil(x)-x))
1768static Image *RenderHoughLines(
const ImageInfo *image_info,
const size_t columns,
1771#define BoundingBox "viewbox"
1785 image=AcquireImage(image_info,exception);
1786 status=OpenBlob(image_info,image,ReadBinaryBlobMode,exception);
1787 if (status == MagickFalse)
1789 image=DestroyImageList(image);
1790 return((
Image *) NULL);
1792 image->columns=columns;
1794 draw_info=CloneDrawInfo(image_info,(
DrawInfo *) NULL);
1795 draw_info->affine.sx=image->resolution.x == 0.0 ? 1.0 : image->resolution.x/
1797 draw_info->affine.sy=image->resolution.y == 0.0 ? 1.0 : image->resolution.y/
1799 image->columns=CastDoubleToSizeT(draw_info->affine.sx*image->columns);
1800 image->rows=CastDoubleToSizeT(draw_info->affine.sy*image->rows);
1801 status=SetImageExtent(image,image->columns,image->rows,exception);
1802 if (status == MagickFalse)
1803 return(DestroyImageList(image));
1804 if (SetImageBackgroundColor(image,exception) == MagickFalse)
1806 image=DestroyImageList(image);
1807 return((
Image *) NULL);
1812 if (GetBlobStreamData(image) == (
unsigned char *) NULL)
1813 draw_info->primitive=FileToString(image->filename,~0UL,exception);
1816 draw_info->primitive=(
char *) AcquireQuantumMemory(1,(
size_t)
1817 GetBlobSize(image)+1);
1818 if (draw_info->primitive != (
char *) NULL)
1820 (void) memcpy(draw_info->primitive,GetBlobStreamData(image),
1821 (size_t) GetBlobSize(image));
1822 draw_info->primitive[GetBlobSize(image)]=
'\0';
1825 (void) DrawImage(image,draw_info,exception);
1826 draw_info=DestroyDrawInfo(draw_info);
1827 if (CloseBlob(image) == MagickFalse)
1828 image=DestroyImageList(image);
1829 return(GetFirstImageInList(image));
1832MagickExport
Image *HoughLineImage(
const Image *image,
const size_t width,
1833 const size_t height,
const size_t threshold,
ExceptionInfo *exception)
1835#define HoughLineImageTag "HoughLine/Image"
1841 message[MagickPathExtent],
1842 path[MagickPathExtent];
1851 *lines_image = NULL;
1882 assert(image != (
const Image *) NULL);
1883 assert(image->signature == MagickCoreSignature);
1885 assert(exception->signature == MagickCoreSignature);
1886 if (IsEventLogging() != MagickFalse)
1887 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
1888 accumulator_width=180;
1889 hough_height=((sqrt(2.0)*(double) (image->rows > image->columns ?
1890 image->rows : image->columns))/2.0);
1891 accumulator_height=(size_t) (2.0*hough_height);
1892 accumulator=AcquireMatrixInfo(accumulator_width,accumulator_height,
1893 sizeof(
double),exception);
1895 ThrowImageException(ResourceLimitError,
"MemoryAllocationFailed");
1896 if (NullMatrix(accumulator) == MagickFalse)
1898 accumulator=DestroyMatrixInfo(accumulator);
1899 ThrowImageException(ResourceLimitError,
"MemoryAllocationFailed");
1906 center.x=(double) image->columns/2.0;
1907 center.y=(double) image->rows/2.0;
1908 image_view=AcquireVirtualCacheView(image,exception);
1909 for (y=0; y < (ssize_t) image->rows; y++)
1917 if (status == MagickFalse)
1919 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1920 if (p == (Quantum *) NULL)
1925 for (x=0; x < (ssize_t) image->columns; x++)
1927 if (GetPixelIntensity(image,p) > ((double) QuantumRange/2.0))
1932 for (i=0; i < 180; i++)
1938 radius=(((double) x-center.x)*cos(DegreesToRadians((
double) i)))+
1939 (((
double) y-center.y)*sin(DegreesToRadians((
double) i)));
1940 (void) GetMatrixElement(accumulator,i,(ssize_t)
1941 MagickRound(radius+hough_height),&count);
1943 (void) SetMatrixElement(accumulator,i,(ssize_t)
1944 MagickRound(radius+hough_height),&count);
1947 p+=(ptrdiff_t) GetPixelChannels(image);
1949 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1954#if defined(MAGICKCORE_OPENMP_SUPPORT)
1958 proceed=SetImageProgress(image,CannyEdgeImageTag,progress,image->rows);
1959 if (proceed == MagickFalse)
1963 image_view=DestroyCacheView(image_view);
1964 if (status == MagickFalse)
1966 accumulator=DestroyMatrixInfo(accumulator);
1967 return((
Image *) NULL);
1972 file=AcquireUniqueFileResource(path);
1975 accumulator=DestroyMatrixInfo(accumulator);
1976 return((
Image *) NULL);
1978 (void) FormatLocaleString(message,MagickPathExtent,
1979 "# Hough line transform: %.20gx%.20g%+.20g\n",(
double) width,
1980 (double) height,(
double) threshold);
1981 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1983 (void) FormatLocaleString(message,MagickPathExtent,
1984 "viewbox 0 0 %.20g %.20g\n",(
double) image->columns,(double) image->rows);
1985 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1987 (void) FormatLocaleString(message,MagickPathExtent,
1988 "# x1,y1 x2,y2 # count angle distance\n");
1989 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1991 line_count=image->columns > image->rows ? image->columns/4 : image->rows/4;
1993 line_count=threshold;
1994 for (y=0; y < (ssize_t) accumulator_height; y++)
1999 for (x=0; x < (ssize_t) accumulator_width; x++)
2004 (void) GetMatrixElement(accumulator,x,y,&count);
2005 if (count >= (
double) line_count)
2020 for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
2025 for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
2027 if ((u != 0) || (v !=0))
2029 (void) GetMatrixElement(accumulator,x+u,y+v,&count);
2037 if (u < (ssize_t) (width/2))
2040 (void) GetMatrixElement(accumulator,x,y,&count);
2043 if ((x >= 45) && (x <= 135))
2049 line.y1=((double) (y-(accumulator_height/2.0))-((line.x1-
2050 (image->columns/2.0))*cos(DegreesToRadians((
double) x))))/
2051 sin(DegreesToRadians((
double) x))+(image->rows/2.0);
2052 line.x2=(double) image->columns;
2053 line.y2=((double) (y-(accumulator_height/2.0))-((line.x2-
2054 (image->columns/2.0))*cos(DegreesToRadians((
double) x))))/
2055 sin(DegreesToRadians((
double) x))+(image->rows/2.0);
2063 line.x1=((double) (y-(accumulator_height/2.0))-((line.y1-
2064 (image->rows/2.0))*sin(DegreesToRadians((
double) x))))/
2065 cos(DegreesToRadians((
double) x))+(image->columns/2.0);
2066 line.y2=(double) image->rows;
2067 line.x2=((double) (y-(accumulator_height/2.0))-((line.y2-
2068 (image->rows/2.0))*sin(DegreesToRadians((
double) x))))/
2069 cos(DegreesToRadians((
double) x))+(image->columns/2.0);
2071 (void) FormatLocaleString(message,MagickPathExtent,
2072 "line %g,%g %g,%g # %g %g %g\n",line.x1,line.y1,line.x2,line.y2,
2073 maxima,(
double) x,(double) y);
2074 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
2079 (void) close_utf8(file);
2083 image_info=AcquireImageInfo();
2084 image_info->background_color=image->background_color;
2085 (void) FormatLocaleString(image_info->filename,MagickPathExtent,
"%s",path);
2086 artifact=GetImageArtifact(image,
"background");
2087 if (artifact != (
const char *) NULL)
2088 (void) SetImageOption(image_info,
"background",artifact);
2089 artifact=GetImageArtifact(image,
"fill");
2090 if (artifact != (
const char *) NULL)
2091 (void) SetImageOption(image_info,
"fill",artifact);
2092 artifact=GetImageArtifact(image,
"stroke");
2093 if (artifact != (
const char *) NULL)
2094 (void) SetImageOption(image_info,
"stroke",artifact);
2095 artifact=GetImageArtifact(image,
"strokewidth");
2096 if (artifact != (
const char *) NULL)
2097 (void) SetImageOption(image_info,
"strokewidth",artifact);
2098 lines_image=RenderHoughLines(image_info,image->columns,image->rows,exception);
2099 artifact=GetImageArtifact(image,
"hough-lines:accumulator");
2100 if ((lines_image != (
Image *) NULL) &&
2101 (IsStringTrue(artifact) != MagickFalse))
2106 accumulator_image=MatrixToImage(accumulator,exception);
2107 if (accumulator_image != (
Image *) NULL)
2108 AppendImageToList(&lines_image,accumulator_image);
2113 accumulator=DestroyMatrixInfo(accumulator);
2114 image_info=DestroyImageInfo(image_info);
2115 (void) RelinquishUniqueFileResource(path);
2116 return(GetFirstImageInList(lines_image));
2159MagickExport
Image *MeanShiftImage(
const Image *image,
const size_t width,
2160 const size_t height,
const double color_distance,
ExceptionInfo *exception)
2162#define MaxMeanShiftIterations 100
2163#define MeanShiftImageTag "MeanShift/Image"
2182 assert(image != (
const Image *) NULL);
2183 assert(image->signature == MagickCoreSignature);
2185 assert(exception->signature == MagickCoreSignature);
2186 if (IsEventLogging() != MagickFalse)
2187 (void) LogMagickEvent(TraceEvent,GetMagickModule(),
"%s",image->filename);
2188 mean_image=CloneImage(image,0,0,MagickTrue,exception);
2189 if (mean_image == (
Image *) NULL)
2190 return((
Image *) NULL);
2191 if (SetImageStorageClass(mean_image,DirectClass,exception) == MagickFalse)
2193 mean_image=DestroyImage(mean_image);
2194 return((
Image *) NULL);
2198 image_view=AcquireVirtualCacheView(image,exception);
2199 pixel_view=AcquireVirtualCacheView(image,exception);
2200 mean_view=AcquireAuthenticCacheView(mean_image,exception);
2201#if defined(MAGICKCORE_OPENMP_SUPPORT)
2202 #pragma omp parallel for schedule(static) shared(status,progress) \
2203 magick_number_threads(mean_image,mean_image,mean_image->rows,1)
2205 for (y=0; y < (ssize_t) mean_image->rows; y++)
2216 if (status == MagickFalse)
2218 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2219 q=GetCacheViewAuthenticPixels(mean_view,0,y,mean_image->columns,1,
2221 if ((p == (
const Quantum *) NULL) || (q == (Quantum *) NULL))
2226 for (x=0; x < (ssize_t) mean_image->columns; x++)
2239 GetPixelInfo(image,&mean_pixel);
2240 GetPixelInfoPixel(image,p,&mean_pixel);
2241 mean_location.x=(double) x;
2242 mean_location.y=(double) y;
2243 for (i=0; i < MaxMeanShiftIterations; i++)
2261 GetPixelInfo(image,&sum_pixel);
2262 previous_location=mean_location;
2263 previous_pixel=mean_pixel;
2265 for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
2270 for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
2272 if ((v*v+u*u) <= (ssize_t) ((width/2)*(height/2)))
2277 status=GetOneCacheViewVirtualPixelInfo(pixel_view,(ssize_t)
2278 MagickRound(mean_location.x+u),(ssize_t) MagickRound(
2279 mean_location.y+v),&pixel,exception);
2280 distance=(mean_pixel.red-pixel.red)*(mean_pixel.red-pixel.red)+
2281 (mean_pixel.green-pixel.green)*(mean_pixel.green-pixel.green)+
2282 (mean_pixel.blue-pixel.blue)*(mean_pixel.blue-pixel.blue);
2283 if (distance <= (color_distance*color_distance))
2285 sum_location.x+=mean_location.x+u;
2286 sum_location.y+=mean_location.y+v;
2287 sum_pixel.red+=pixel.red;
2288 sum_pixel.green+=pixel.green;
2289 sum_pixel.blue+=pixel.blue;
2290 sum_pixel.alpha+=pixel.alpha;
2297 gamma=MagickSafeReciprocal((
double) count);
2298 mean_location.x=gamma*sum_location.x;
2299 mean_location.y=gamma*sum_location.y;
2300 mean_pixel.red=gamma*sum_pixel.red;
2301 mean_pixel.green=gamma*sum_pixel.green;
2302 mean_pixel.blue=gamma*sum_pixel.blue;
2303 mean_pixel.alpha=gamma*sum_pixel.alpha;
2304 distance=(mean_location.x-previous_location.x)*
2305 (mean_location.x-previous_location.x)+
2306 (mean_location.y-previous_location.y)*
2307 (mean_location.y-previous_location.y)+
2308 255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)*
2309 255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)+
2310 255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)*
2311 255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)+
2312 255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue)*
2313 255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue);
2314 if (distance <= 3.0)
2317 SetPixelRed(mean_image,ClampToQuantum(mean_pixel.red),q);
2318 SetPixelGreen(mean_image,ClampToQuantum(mean_pixel.green),q);
2319 SetPixelBlue(mean_image,ClampToQuantum(mean_pixel.blue),q);
2320 SetPixelAlpha(mean_image,ClampToQuantum(mean_pixel.alpha),q);
2321 p+=(ptrdiff_t) GetPixelChannels(image);
2322 q+=(ptrdiff_t) GetPixelChannels(mean_image);
2324 if (SyncCacheViewAuthenticPixels(mean_view,exception) == MagickFalse)
2326 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2331#if defined(MAGICKCORE_OPENMP_SUPPORT)
2335 proceed=SetImageProgress(image,MeanShiftImageTag,progress,image->rows);
2336 if (proceed == MagickFalse)
2340 mean_view=DestroyCacheView(mean_view);
2341 pixel_view=DestroyCacheView(pixel_view);
2342 image_view=DestroyCacheView(image_view);