Support -kmeans command-line option
diff --git a/MagickCore/quantize.c b/MagickCore/quantize.c
index 58f948c..410b573 100644
--- a/MagickCore/quantize.c
+++ b/MagickCore/quantize.c
@@ -182,6 +182,7 @@
#include "MagickCore/colormap.h"
#include "MagickCore/colorspace.h"
#include "MagickCore/colorspace-private.h"
+#include "MagickCore/compare.h"
#include "MagickCore/enhance.h"
#include "MagickCore/exception.h"
#include "MagickCore/exception-private.h"
@@ -2305,6 +2306,267 @@
% %
% %
% %
+% K m e a n s I m a g e %
+% %
+% %
+% %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% KmeansImage() applies k-means color reduction to an image. This is a
+% colorspace clustering or segmentation technique.
+%
+% The format of the KmeansImage method is:
+%
+% MagickBooleanType KmeansImage(Image *image,const size_t number_colors,
+% const size_t max_iterations,const double max_distortion,
+% ExceptionInfo *exception)
+%
+% A description of each parameter follows:
+%
+% o image: the image.
+%
+% o number_colors: number of colors to use as seeds.
+%
+% o max_iterations: maximum number of iterations while converging.
+%
+% o max_distortion: the maximum quantization distortion.
+%
+% o exception: return any errors or warnings in this structure.
+%
+*/
+
+static inline double KmeansDistance(const Image *magick_restrict image,
+ const Quantum *magick_restrict p,const PixelInfo *magick_restrict q)
+{
+ double
+ distance,
+ gamma,
+ pixel;
+
+ gamma=1.0;
+ distance=0.0;
+ if (image->alpha_trait != UndefinedPixelTrait)
+ {
+ pixel=GetPixelAlpha(image,p)-q->alpha;
+ distance+=gamma*QuantumScale*pixel*pixel;
+ gamma=(QuantumScale*GetPixelAlpha(image,p))*(QuantumScale*q->alpha);
+ }
+ distance=0.0;
+ if (image->colorspace == CMYKColorspace)
+ {
+ pixel=GetPixelBlack(image,p)-q->black;
+ distance+=gamma*QuantumScale*pixel*pixel;
+ gamma*=(double) (QuantumScale*(QuantumRange-GetPixelBlack(image,p)));
+ gamma*=(double) (QuantumScale*(QuantumRange-q->black));
+ }
+ distance*=3.0;
+ pixel=GetPixelRed(image,p)-q->red;
+ if ((image->colorspace == HSLColorspace) ||
+ (image->colorspace == HSBColorspace) ||
+ (image->colorspace == HWBColorspace))
+ {
+ if (fabs((double) pixel) > (QuantumRange/2))
+ pixel-=QuantumRange;
+ pixel*=2.0;
+ }
+ distance+=gamma*QuantumScale*pixel*pixel;
+ pixel=GetPixelGreen(image,p)-q->green;
+ distance+=gamma*QuantumScale*pixel*pixel;
+ pixel=GetPixelBlue(image,p)-q->blue;
+ distance+=gamma*QuantumScale*pixel*pixel;
+ return(distance);
+}
+
+MagickExport MagickBooleanType KmeansImage(Image *image,
+ const size_t number_colors,const size_t max_iterations,
+ const double max_distortion,ExceptionInfo *exception)
+{
+#define KmeansImageTag "Kmeans/Image"
+
+ double
+ last_distortion;
+
+ Image
+ *kmeans_image;
+
+ MagickBooleanType
+ status;
+
+ QuantizeInfo
+ *quantize_info;
+
+ register ssize_t
+ n;
+
+ /*
+ Initialize the initial set of k means.
+ */
+ assert(image != (Image *) NULL);
+ assert(image->signature == MagickCoreSignature);
+ if (image->debug != MagickFalse)
+ (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
+ assert(exception != (ExceptionInfo *) NULL);
+ assert(exception->signature == MagickCoreSignature);
+ kmeans_image=CloneImage(image,0,0,MagickTrue,exception);
+ if (kmeans_image == (Image *) NULL)
+ return(MagickFalse);
+ quantize_info=AcquireQuantizeInfo((ImageInfo *) NULL);
+ quantize_info->number_colors=number_colors;
+ quantize_info->dither_method=NoDitherMethod;
+ status=QuantizeImage(quantize_info,kmeans_image,exception);
+ quantize_info=DestroyQuantizeInfo(quantize_info);
+ if (status == MagickFalse)
+ return(status);
+ if (AcquireImageColormap(image,kmeans_image->colors,exception) == MagickFalse)
+ {
+ kmeans_image=DestroyImage(kmeans_image);
+ return(MagickFalse);
+ }
+ (void) memcpy(image->colormap,kmeans_image->colormap,kmeans_image->colors*
+ sizeof(*kmeans_image->colormap));
+ /*
+ Iterative refinement.
+ */
+ last_distortion=0.0;
+ for (n=0; n < max_iterations; n++)
+ {
+ CacheView
+ *image_view;
+
+ double
+ distortion;
+
+ PixelInfo
+ *kmeans_colormap;
+
+ register ssize_t
+ i;
+
+ ssize_t
+ y;
+
+ kmeans_colormap=kmeans_image->colormap;
+ for (i=0; i < (ssize_t) kmeans_image->colors; i++)
+ {
+ kmeans_colormap[i].red=0.0;
+ kmeans_colormap[i].green=0.0;
+ kmeans_colormap[i].blue=0.0;
+ kmeans_colormap[i].alpha=0.0;
+ kmeans_colormap[i].black=0.0;
+ kmeans_colormap[i].count=0.0;
+ kmeans_colormap[i].distortion=0.0;
+ }
+ image_view=AcquireAuthenticCacheView(image,exception);
+#if defined(MAGICKCORE_OPENMP_SUPPORT)
+ #pragma omp parallel for schedule(static) shared(status) \
+ magick_number_threads(image,image,image->rows,1)
+#endif
+ for (y=0; y < (ssize_t) image->rows; y++)
+ {
+ register Quantum
+ *magick_restrict q;
+
+ register ssize_t
+ x;
+
+ if (status == MagickFalse)
+ continue;
+ q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
+ if (q == (Quantum *) NULL)
+ {
+ status=MagickFalse;
+ continue;
+ }
+ for (x=0; x < (ssize_t) image->columns; x++)
+ {
+ double
+ min_distance;
+
+ register ssize_t
+ i;
+
+ ssize_t
+ j;
+
+ /*
+ Assign each pixel whose mean has the least squared Euclidean distance.
+ */
+ j=0;
+ min_distance=KmeansDistance(image,q,image->colormap+0);
+ for (i=1; i < (ssize_t) image->colors; i++)
+ {
+ double
+ distance;
+
+ if (min_distance <= MagickEpsilon)
+ break;
+ distance=KmeansDistance(image,q,image->colormap+i);
+ if (distance < min_distance)
+ {
+ min_distance=distance;
+ j=i;
+ }
+ }
+ kmeans_colormap[j].red+=QuantumScale*GetPixelRed(image,q);
+ kmeans_colormap[j].green+=QuantumScale*GetPixelGreen(image,q);
+ kmeans_colormap[j].blue+=QuantumScale*GetPixelBlue(image,q);
+ if (kmeans_image->alpha_trait != BlendPixelTrait)
+ kmeans_colormap[j].alpha+=QuantumScale*GetPixelAlpha(image,q);
+ if (kmeans_image->colorspace == CMYKColorspace)
+ kmeans_colormap[j].black+=QuantumScale*GetPixelBlack(image,q);
+ kmeans_colormap[j].count++;
+ kmeans_colormap[j].distortion+=min_distance;
+ SetPixelIndex(image,j,q);
+ q+=GetPixelChannels(image);
+ }
+ }
+ /*
+ Calculate the new means (centroids) of the pixels in the new clusters.
+ */
+ distortion=0.0;
+ for (i=0; i < (ssize_t) kmeans_image->colors; i++)
+ {
+ image->colormap[i].red=QuantumRange*kmeans_colormap[i].red/(double)
+ kmeans_colormap[i].count;
+ image->colormap[i].green=QuantumRange*kmeans_colormap[i].green/(double)
+ kmeans_colormap[i].count;
+ image->colormap[i].blue=QuantumRange*kmeans_colormap[i].blue/(double)
+ kmeans_colormap[i].count;
+ if (kmeans_image->alpha_trait != BlendPixelTrait)
+ image->colormap[i].alpha=QuantumRange*kmeans_colormap[i].alpha/(double)
+ kmeans_colormap[i].count;
+ if (kmeans_image->colorspace == CMYKColorspace)
+ image->colormap[i].black=QuantumRange*kmeans_colormap[i].black/(double)
+ kmeans_colormap[i].count;
+ distortion+=kmeans_colormap[i].distortion/kmeans_image->colors;
+ }
+ image_view=DestroyCacheView(image_view);
+ if (fabs(distortion-last_distortion) <= max_distortion)
+ break;
+ last_distortion=distortion;
+ if (image->progress_monitor != (MagickProgressMonitor) NULL)
+ {
+ MagickBooleanType
+ proceed;
+
+ proceed=SetImageProgress(image,KmeansImageTag,n,max_iterations);
+ if (proceed == MagickFalse)
+ status=MagickFalse;
+ }
+ }
+ if (image->progress_monitor != (MagickProgressMonitor) NULL)
+ (void) SetImageProgress(image,KmeansImageTag,max_iterations+1,
+ max_iterations);
+ if (status != MagickFalse)
+ SyncImage(image,exception);
+ return(status);
+}
+
+/*
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% %
+% %
+% %
% P o s t e r i z e I m a g e %
% %
% %
