root/magick/quantize.c

DEFINITIONS

1. AcquireQuantizeInfo
2. AssociateAlphaPixel
3. ColorToNodeId
4. IsSameColor
5. AssignImageColors
6. SetAssociatedAlpha
7. ClassifyImageColors
8. CloneQuantizeInfo
9. ClosestColor
10. CompressImageColormap
11. DefineImageColormap
12. DestroyCubeInfo
13. DestroyQuantizeInfo
14. DestroyPixelThreadSet
15. AcquirePixelThreadSet
16. CacheOffset
17. FloydSteinbergDither
18. Riemersma
19. RiemersmaDither
20. DitherImage
21. GetCubeInfo
22. GetNodeInfo
23. GetImageQuantizeError
24. GetQuantizeInfo
25. MagickRound
26. PosterizeImage
27. PosterizeImageChannel
28. PruneChild
29. PruneLevel
30. PruneToCubeDepth
31. DirectToColormapImage
32. QuantizeImage
33. QuantizeImages
34. QuantizeErrorFlatten
35. Reduce
36. MagickRealTypeCompare
37. ReduceImageColors
38. RemapImage
39. RemapImages
40. IntensityCompare
41. SetGrayscaleImage

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%           QQQ   U   U   AAA   N   N  TTTTT  IIIII   ZZZZZ  EEEEE            %
%          Q   Q  U   U  A   A  NN  N    T      I        ZZ  E                %
%          Q   Q  U   U  AAAAA  N N N    T      I      ZZZ   EEEEE            %
%          Q  QQ  U   U  A   A  N  NN    T      I     ZZ     E                %
%           QQQQ   UUU   A   A  N   N    T    IIIII   ZZZZZ  EEEEE            %
%                                                                             %
%                                                                             %
%    MagickCore Methods to Reduce the Number of Unique Colors in an Image     %
%                                                                             %
%                           Software Design                                   %
%                                Cristy                                       %
%                              July 1992                                      %
%                                                                             %
%                                                                             %
%  Copyright 1999-2016 ImageMagick Studio LLC, a non-profit organization      %
%  dedicated to making software imaging solutions freely available.           %
%                                                                             %
%  You may not use this file except in compliance with the License.  You may  %
%  obtain a copy of the License at                                            %
%                                                                             %
%    http://www.imagemagick.org/script/license.php                            %
%                                                                             %
%  Unless required by applicable law or agreed to in writing, software        %
%  distributed under the License is distributed on an "AS IS" BASIS,          %
%  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
%  See the License for the specific language governing permissions and        %
%  limitations under the License.                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Realism in computer graphics typically requires using 24 bits/pixel to
%  generate an image.  Yet many graphic display devices do not contain the
%  amount of memory necessary to match the spatial and color resolution of
%  the human eye.  The Quantize methods takes a 24 bit image and reduces
%  the number of colors so it can be displayed on raster device with less
%  bits per pixel.  In most instances, the quantized image closely
%  resembles the original reference image.
%
%  A reduction of colors in an image is also desirable for image
%  transmission and real-time animation.
%
%  QuantizeImage() takes a standard RGB or monochrome images and quantizes
%  them down to some fixed number of colors.
%
%  For purposes of color allocation, an image is a set of n pixels, where
%  each pixel is a point in RGB space.  RGB space is a 3-dimensional
%  vector space, and each pixel, Pi,  is defined by an ordered triple of
%  red, green, and blue coordinates, (Ri, Gi, Bi).
%
%  Each primary color component (red, green, or blue) represents an
%  intensity which varies linearly from 0 to a maximum value, Cmax, which
%  corresponds to full saturation of that color.  Color allocation is
%  defined over a domain consisting of the cube in RGB space with opposite
%  vertices at (0,0,0) and (Cmax, Cmax, Cmax).  QUANTIZE requires Cmax =
%  255.
%
%  The algorithm maps this domain onto a tree in which each node
%  represents a cube within that domain.  In the following discussion
%  these cubes are defined by the coordinate of two opposite vertices (vertex
%  nearest the origin in RGB space and the vertex farthest from the origin).
%
%  The tree's root node represents the entire domain, (0,0,0) through
%  (Cmax,Cmax,Cmax).  Each lower level in the tree is generated by
%  subdividing one node's cube into eight smaller cubes of equal size.
%  This corresponds to bisecting the parent cube with planes passing
%  through the midpoints of each edge.
%
%  The basic algorithm operates in three phases: Classification,
%  Reduction, and Assignment.  Classification builds a color description
%  tree for the image.  Reduction collapses the tree until the number it
%  represents, at most, the number of colors desired in the output image.
%  Assignment defines the output image's color map and sets each pixel's
%  color by restorage_class in the reduced tree.  Our goal is to minimize
%  the numerical discrepancies between the original colors and quantized
%  colors (quantization error).
%
%  Classification begins by initializing a color description tree of
%  sufficient depth to represent each possible input color in a leaf.
%  However, it is impractical to generate a fully-formed color description
%  tree in the storage_class phase for realistic values of Cmax.  If
%  colors components in the input image are quantized to k-bit precision,
%  so that Cmax= 2k-1, the tree would need k levels below the root node to
%  allow representing each possible input color in a leaf.  This becomes
%  prohibitive because the tree's total number of nodes is 1 +
%  sum(i=1, k, 8k).
%
%  A complete tree would require 19,173,961 nodes for k = 8, Cmax = 255.
%  Therefore, to avoid building a fully populated tree, QUANTIZE: (1)
%  Initializes data structures for nodes only as they are needed;  (2)
%  Chooses a maximum depth for the tree as a function of the desired
%  number of colors in the output image (currently log2(colormap size)).
%
%  For each pixel in the input image, storage_class scans downward from
%  the root of the color description tree.  At each level of the tree it
%  identifies the single node which represents a cube in RGB space
%  containing the pixel's color.  It updates the following data for each
%  such node:
%
%    n1: Number of pixels whose color is contained in the RGB cube which
%    this node represents;
%
%    n2: Number of pixels whose color is not represented in a node at
%    lower depth in the tree;  initially,  n2 = 0 for all nodes except
%    leaves of the tree.
%
%    Sr, Sg, Sb: Sums of the red, green, and blue component values for all
%    pixels not classified at a lower depth. The combination of these sums
%    and n2  will ultimately characterize the mean color of a set of
%    pixels represented by this node.
%
%    E: the distance squared in RGB space between each pixel contained
%    within a node and the nodes' center.  This represents the
%    quantization error for a node.
%
%  Reduction repeatedly prunes the tree until the number of nodes with n2
%  > 0 is less than or equal to the maximum number of colors allowed in
%  the output image.  On any given iteration over the tree, it selects
%  those nodes whose E count is minimal for pruning and merges their color
%  statistics upward. It uses a pruning threshold, Ep, to govern node
%  selection as follows:
%
%    Ep = 0
%    while number of nodes with (n2 > 0) > required maximum number of colors
%      prune all nodes such that E <= Ep
%      Set Ep to minimum E in remaining nodes
%
%  This has the effect of minimizing any quantization error when merging
%  two nodes together.
%
%  When a node to be pruned has offspring, the pruning procedure invokes
%  itself recursively in order to prune the tree from the leaves upward.
%  n2,  Sr, Sg,  and  Sb in a node being pruned are always added to the
%  corresponding data in that node's parent.  This retains the pruned
%  node's color characteristics for later averaging.
%
%  For each node, n2 pixels exist for which that node represents the
%  smallest volume in RGB space containing those pixel's colors.  When n2
%  > 0 the node will uniquely define a color in the output image. At the
%  beginning of reduction,  n2 = 0  for all nodes except a the leaves of
%  the tree which represent colors present in the input image.
%
%  The other pixel count, n1, indicates the total number of colors within
%  the cubic volume which the node represents.  This includes n1 - n2
%  pixels whose colors should be defined by nodes at a lower level in the
%  tree.
%
%  Assignment generates the output image from the pruned tree.  The output
%  image consists of two parts: (1)  A color map, which is an array of
%  color descriptions (RGB triples) for each color present in the output
%  image;  (2)  A pixel array, which represents each pixel as an index
%  into the color map array.
%
%  First, the assignment phase makes one pass over the pruned color
%  description tree to establish the image's color map.  For each node
%  with n2  > 0, it divides Sr, Sg, and Sb by n2 .  This produces the mean
%  color of all pixels that classify no lower than this node.  Each of
%  these colors becomes an entry in the color map.
%
%  Finally,  the assignment phase reclassifies each pixel in the pruned
%  tree to identify the deepest node containing the pixel's color.  The
%  pixel's value in the pixel array becomes the index of this node's mean
%  color in the color map.
%
%  This method is based on a similar algorithm written by Paul Raveling.
%
*/

/*
  Include declarations.
*/
#include "magick/studio.h"
#include "magick/attribute.h"
#include "magick/cache-view.h"
#include "magick/color.h"
#include "magick/color-private.h"
#include "magick/colormap.h"
#include "magick/colorspace.h"
#include "magick/colorspace-private.h"
#include "magick/enhance.h"
#include "magick/exception.h"
#include "magick/exception-private.h"
#include "magick/histogram.h"
#include "magick/image.h"
#include "magick/image-private.h"
#include "magick/list.h"
#include "magick/memory_.h"
#include "magick/monitor.h"
#include "magick/monitor-private.h"
#include "magick/option.h"
#include "magick/pixel-private.h"
#include "magick/quantize.h"
#include "magick/quantum.h"
#include "magick/resource_.h"
#include "magick/string_.h"
#include "magick/thread-private.h"

/*
  Define declarations.
*/
#if !defined(__APPLE__) && !defined(TARGET_OS_IPHONE)
#define CacheShift  2
#else
#define CacheShift  3
#endif
#define ErrorQueueLength  16
#define MaxNodes  266817
#define MaxTreeDepth  8
#define NodesInAList  1920

/*
  Typdef declarations.
*/
typedef struct _NodeInfo
{
  struct _NodeInfo
    *parent,
    *child[16];

  MagickSizeType
    number_unique;

  DoublePixelPacket
    total_color;

  MagickRealType
    quantize_error;

  size_t
    color_number,
    id,
    level;
} NodeInfo;

typedef struct _Nodes
{
  NodeInfo
    *nodes;

  struct _Nodes
    *next;
} Nodes;

typedef struct _CubeInfo
{
  NodeInfo
    *root;

  size_t
    colors,
    maximum_colors;

  ssize_t
    transparent_index;

  MagickSizeType
    transparent_pixels;

  DoublePixelPacket
    target;

  MagickRealType
    distance,
    pruning_threshold,
    next_threshold;

  size_t
    nodes,
    free_nodes,
    color_number;

  NodeInfo
    *next_node;

  Nodes
    *node_queue;

  MemoryInfo
    *memory_info;

  ssize_t
    *cache;

  DoublePixelPacket
    error[ErrorQueueLength];

  MagickRealType
    weights[ErrorQueueLength];

  QuantizeInfo
    *quantize_info;

  MagickBooleanType
    associate_alpha;

  ssize_t
    x,
    y;

  size_t
    depth;

  MagickOffsetType
    offset;

  MagickSizeType
    span;
} CubeInfo;

/*
  Method prototypes.
*/
static CubeInfo
  *GetCubeInfo(const QuantizeInfo *,const size_t,const size_t);

static NodeInfo
  *GetNodeInfo(CubeInfo *,const size_t,const size_t,NodeInfo *);

static MagickBooleanType
  AssignImageColors(Image *,CubeInfo *),
  ClassifyImageColors(CubeInfo *,const Image *,ExceptionInfo *),
  DitherImage(Image *,CubeInfo *),
  SetGrayscaleImage(Image *);

static size_t
  DefineImageColormap(Image *,CubeInfo *,NodeInfo *);

static void
  ClosestColor(const Image *,CubeInfo *,const NodeInfo *),
  DestroyCubeInfo(CubeInfo *),
  PruneLevel(CubeInfo *,const NodeInfo *),
  PruneToCubeDepth(CubeInfo *,const NodeInfo *),
  ReduceImageColors(const Image *,CubeInfo *);

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   A c q u i r e Q u a n t i z e I n f o                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AcquireQuantizeInfo() allocates the QuantizeInfo structure.
%
%  The format of the AcquireQuantizeInfo method is:
%
%      QuantizeInfo *AcquireQuantizeInfo(const ImageInfo *image_info)
%
%  A description of each parameter follows:
%
%    o image_info: the image info.
%
*/
MagickExport QuantizeInfo *AcquireQuantizeInfo(const ImageInfo *image_info)
{
  QuantizeInfo
    *quantize_info;

  quantize_info=(QuantizeInfo *) AcquireMagickMemory(sizeof(*quantize_info));
  if (quantize_info == (QuantizeInfo *) NULL)
    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
  GetQuantizeInfo(quantize_info);
  if (image_info != (ImageInfo *) NULL)
    {
      const char
        *option;

      quantize_info->dither=image_info->dither;
      option=GetImageOption(image_info,"dither");
      if (option != (const char *) NULL)
        quantize_info->dither_method=(DitherMethod) ParseCommandOption(
          MagickDitherOptions,MagickFalse,option);
      quantize_info->measure_error=image_info->verbose;
    }
  return(quantize_info);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   A s s i g n I m a g e C o l o r s                                         %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  AssignImageColors() generates the output image from the pruned tree.  The
%  output image consists of two parts: (1)  A color map, which is an array
%  of color descriptions (RGB triples) for each color present in the
%  output image;  (2)  A pixel array, which represents each pixel as an
%  index into the color map array.
%
%  First, the assignment phase makes one pass over the pruned color
%  description tree to establish the image's color map.  For each node
%  with n2  > 0, it divides Sr, Sg, and Sb by n2 .  This produces the mean
%  color of all pixels that classify no lower than this node.  Each of
%  these colors becomes an entry in the color map.
%
%  Finally,  the assignment phase reclassifies each pixel in the pruned
%  tree to identify the deepest node containing the pixel's color.  The
%  pixel's value in the pixel array becomes the index of this node's mean
%  color in the color map.
%
%  The format of the AssignImageColors() method is:
%
%      MagickBooleanType AssignImageColors(Image *image,CubeInfo *cube_info)
%
%  A description of each parameter follows.
%
%    o image: the image.
%
%    o cube_info: A pointer to the Cube structure.
%
*/

static inline void AssociateAlphaPixel(const CubeInfo *cube_info,
  const PixelPacket *pixel,DoublePixelPacket *alpha_pixel)
{
  MagickRealType
    alpha;

  if ((cube_info->associate_alpha == MagickFalse) ||
      (pixel->opacity == OpaqueOpacity))
    {
      alpha_pixel->red=(MagickRealType) GetPixelRed(pixel);
      alpha_pixel->green=(MagickRealType) GetPixelGreen(pixel);
      alpha_pixel->blue=(MagickRealType) GetPixelBlue(pixel);
      alpha_pixel->opacity=(MagickRealType) GetPixelOpacity(pixel);
      return;
    }
  alpha=(MagickRealType) (QuantumScale*(QuantumRange-GetPixelOpacity(pixel)));
  alpha_pixel->red=alpha*GetPixelRed(pixel);
  alpha_pixel->green=alpha*GetPixelGreen(pixel);
  alpha_pixel->blue=alpha*GetPixelBlue(pixel);
  alpha_pixel->opacity=(MagickRealType) GetPixelOpacity(pixel);
}

static inline size_t ColorToNodeId(const CubeInfo *cube_info,
  const DoublePixelPacket *pixel,size_t index)
{
  size_t
    id;

  id=(size_t) (((ScaleQuantumToChar(ClampPixel(GetPixelRed(pixel))) >> index) &
    0x01) | ((ScaleQuantumToChar(ClampPixel(GetPixelGreen(pixel))) >> index) &
    0x01) << 1 | ((ScaleQuantumToChar(ClampPixel(GetPixelBlue(pixel))) >>
    index) & 0x01) << 2);
  if (cube_info->associate_alpha != MagickFalse)
    id|=((ScaleQuantumToChar(ClampPixel(GetPixelOpacity(pixel))) >> index) &
      0x1) << 3;
  return(id);
}

static inline MagickBooleanType IsSameColor(const Image *image,
  const PixelPacket *p,const PixelPacket *q)
{
  if ((GetPixelRed(p) != GetPixelRed(q)) ||
      (GetPixelGreen(p) != GetPixelGreen(q)) ||
      (GetPixelBlue(p) != GetPixelBlue(q)))
    return(MagickFalse);
  if ((image->matte != MagickFalse) &&
      (GetPixelOpacity(p) != GetPixelOpacity(q)))
    return(MagickFalse);
  return(MagickTrue);
}

static MagickBooleanType AssignImageColors(Image *image,CubeInfo *cube_info)
{
#define AssignImageTag  "Assign/Image"

  ssize_t
    y;

  /*
    Allocate image colormap.
  */
  if ((cube_info->quantize_info->colorspace != UndefinedColorspace) &&
      (cube_info->quantize_info->colorspace != CMYKColorspace))
    (void) TransformImageColorspace(image,cube_info->quantize_info->colorspace);
  else
    if (IssRGBCompatibleColorspace(image->colorspace) == MagickFalse)
      (void) TransformImageColorspace(image,sRGBColorspace);
  if (AcquireImageColormap(image,cube_info->colors) == MagickFalse)
    ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
      image->filename);
  image->colors=0;
  cube_info->transparent_pixels=0;
  cube_info->transparent_index=(-1);
  (void) DefineImageColormap(image,cube_info,cube_info->root);
  /*
    Create a reduced color image.
  */
  if ((cube_info->quantize_info->dither != MagickFalse) &&
      (cube_info->quantize_info->dither_method != NoDitherMethod))
    (void) DitherImage(image,cube_info);
  else
    {
      CacheView
        *image_view;

      ExceptionInfo
        *exception;

      MagickBooleanType
        status;

      status=MagickTrue;
      exception=(&image->exception);
      image_view=AcquireAuthenticCacheView(image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
      #pragma omp parallel for schedule(static,4) shared(status) \
        magick_threads(image,image,image->rows,1)
#endif
      for (y=0; y < (ssize_t) image->rows; y++)
      {
        CubeInfo
          cube;

        register IndexPacket
          *magick_restrict indexes;

        register PixelPacket
          *magick_restrict q;

        register ssize_t
          x;

        ssize_t
          count;

        if (status == MagickFalse)
          continue;
        q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,
          exception);
        if (q == (PixelPacket *) NULL)
          {
            status=MagickFalse;
            continue;
          }
        indexes=GetCacheViewAuthenticIndexQueue(image_view);
        cube=(*cube_info);
        for (x=0; x < (ssize_t) image->columns; x+=count)
        {
          DoublePixelPacket
            pixel;

          register const NodeInfo
            *node_info;

          register ssize_t
            i;

          size_t
            id,
            index;

          /*
            Identify the deepest node containing the pixel's color.
          */
          for (count=1; (x+count) < (ssize_t) image->columns; count++)
            if (IsSameColor(image,q,q+count) == MagickFalse)
              break;
          AssociateAlphaPixel(&cube,q,&pixel);
          node_info=cube.root;
          for (index=MaxTreeDepth-1; (ssize_t) index > 0; index--)
          {
            id=ColorToNodeId(&cube,&pixel,index);
            if (node_info->child[id] == (NodeInfo *) NULL)
              break;
            node_info=node_info->child[id];
          }
          /*
            Find closest color among siblings and their children.
          */
          cube.target=pixel;
          cube.distance=(MagickRealType) (4.0*(QuantumRange+1.0)*
            (QuantumRange+1.0)+1.0);
          ClosestColor(image,&cube,node_info->parent);
          index=cube.color_number;
          for (i=0; i < (ssize_t) count; i++)
          {
            if (image->storage_class == PseudoClass)
              SetPixelIndex(indexes+x+i,index);
            if (cube.quantize_info->measure_error == MagickFalse)
              {
                SetPixelRgb(q,image->colormap+index);
                if (cube.associate_alpha != MagickFalse)
                  SetPixelOpacity(q,image->colormap[index].opacity);
              }
            q++;
          }
        }
        if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
          status=MagickFalse;
        if (image->progress_monitor != (MagickProgressMonitor) NULL)
          {
            MagickBooleanType
              proceed;

#if defined(MAGICKCORE_OPENMP_SUPPORT)
            #pragma omp critical (MagickCore_AssignImageColors)
#endif
            proceed=SetImageProgress(image,AssignImageTag,(MagickOffsetType) y,
              image->rows);
            if (proceed == MagickFalse)
              status=MagickFalse;
          }
      }
      image_view=DestroyCacheView(image_view);
    }
  if (cube_info->quantize_info->measure_error != MagickFalse)
    (void) GetImageQuantizeError(image);
  if ((cube_info->quantize_info->number_colors == 2) &&
      (cube_info->quantize_info->colorspace == GRAYColorspace))
    {
      double
        intensity;

      /*
        Monochrome image.
      */
      intensity=0.0;
      if (GetPixelLuma(image,image->colormap+0) > 
          GetPixelLuma(image,image->colormap+1))
        intensity=(double) QuantumRange;
      image->colormap[0].red=intensity;
      image->colormap[0].green=intensity;
      image->colormap[0].blue=intensity;
      image->colormap[1].red=(double) QuantumRange-intensity;
      image->colormap[1].green=(double) QuantumRange-intensity;
      image->colormap[1].blue=(double) QuantumRange-intensity;
    }
  (void) SyncImage(image);
  if ((cube_info->quantize_info->colorspace != UndefinedColorspace) &&
      (cube_info->quantize_info->colorspace != CMYKColorspace))
    (void) TransformImageColorspace((Image *) image,sRGBColorspace);
  return(MagickTrue);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   C l a s s i f y I m a g e C o l o r s                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ClassifyImageColors() begins by initializing a color description tree
%  of sufficient depth to represent each possible input color in a leaf.
%  However, it is impractical to generate a fully-formed color
%  description tree in the storage_class phase for realistic values of
%  Cmax.  If colors components in the input image are quantized to k-bit
%  precision, so that Cmax= 2k-1, the tree would need k levels below the
%  root node to allow representing each possible input color in a leaf.
%  This becomes prohibitive because the tree's total number of nodes is
%  1 + sum(i=1,k,8k).
%
%  A complete tree would require 19,173,961 nodes for k = 8, Cmax = 255.
%  Therefore, to avoid building a fully populated tree, QUANTIZE: (1)
%  Initializes data structures for nodes only as they are needed;  (2)
%  Chooses a maximum depth for the tree as a function of the desired
%  number of colors in the output image (currently log2(colormap size)).
%
%  For each pixel in the input image, storage_class scans downward from
%  the root of the color description tree.  At each level of the tree it
%  identifies the single node which represents a cube in RGB space
%  containing It updates the following data for each such node:
%
%    n1 : Number of pixels whose color is contained in the RGB cube
%    which this node represents;
%
%    n2 : Number of pixels whose color is not represented in a node at
%    lower depth in the tree;  initially,  n2 = 0 for all nodes except
%    leaves of the tree.
%
%    Sr, Sg, Sb : Sums of the red, green, and blue component values for
%    all pixels not classified at a lower depth. The combination of
%    these sums and n2 will ultimately characterize the mean color of a
%    set of pixels represented by this node.
%
%    E: the distance squared in RGB space between each pixel contained
%    within a node and the nodes' center.  This represents the quantization
%    error for a node.
%
%  The format of the ClassifyImageColors() method is:
%
%      MagickBooleanType ClassifyImageColors(CubeInfo *cube_info,
%        const Image *image,ExceptionInfo *exception)
%
%  A description of each parameter follows.
%
%    o cube_info: A pointer to the Cube structure.
%
%    o image: the image.
%
*/

static inline void SetAssociatedAlpha(const Image *image,CubeInfo *cube_info)
{
  MagickBooleanType
    associate_alpha;

  associate_alpha=image->matte;
  if ((cube_info->quantize_info->number_colors == 2) &&
      (cube_info->quantize_info->colorspace == GRAYColorspace))
    associate_alpha=MagickFalse;
  cube_info->associate_alpha=associate_alpha;
}

static MagickBooleanType ClassifyImageColors(CubeInfo *cube_info,
  const Image *image,ExceptionInfo *exception)
{
#define ClassifyImageTag  "Classify/Image"

  CacheView
    *image_view;

  DoublePixelPacket
    error,
    mid,
    midpoint,
    pixel;

  MagickBooleanType
    proceed;

  MagickRealType
    bisect;

  NodeInfo
    *node_info;

  size_t
    count,
    id,
    index,
    level;

  ssize_t
    y;

  /*
    Classify the first cube_info->maximum_colors colors to a tree depth of 8.
  */
  SetAssociatedAlpha(image,cube_info);
  if ((cube_info->quantize_info->colorspace != UndefinedColorspace) &&
      (cube_info->quantize_info->colorspace != CMYKColorspace))
    (void) TransformImageColorspace((Image *) image,
      cube_info->quantize_info->colorspace);
  else
    if (IssRGBCompatibleColorspace(image->colorspace) == MagickFalse)
      (void) TransformImageColorspace((Image *) image,sRGBColorspace);
  midpoint.red=(MagickRealType) QuantumRange/2.0;
  midpoint.green=(MagickRealType) QuantumRange/2.0;
  midpoint.blue=(MagickRealType) QuantumRange/2.0;
  midpoint.opacity=(MagickRealType) QuantumRange/2.0;
  error.opacity=0.0;
  image_view=AcquireVirtualCacheView(image,exception);
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    register const PixelPacket
      *magick_restrict p;

    register ssize_t
      x;

    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
    if (p == (const PixelPacket *) NULL)
      break;
    if (cube_info->nodes > MaxNodes)
      {
        /*
          Prune one level if the color tree is too large.
        */
        PruneLevel(cube_info,cube_info->root);
        cube_info->depth--;
      }
    for (x=0; x < (ssize_t) image->columns; x+=(ssize_t) count)
    {
      /*
        Start at the root and descend the color cube tree.
      */
      for (count=1; (x+(ssize_t) count) < (ssize_t) image->columns; count++)
        if (IsSameColor(image,p,p+count) == MagickFalse)
          break;
      AssociateAlphaPixel(cube_info,p,&pixel);
      index=MaxTreeDepth-1;
      bisect=((MagickRealType) QuantumRange+1.0)/2.0;
      mid=midpoint;
      node_info=cube_info->root;
      for (level=1; level <= MaxTreeDepth; level++)
      {
        double
          distance;

        bisect*=0.5;
        id=ColorToNodeId(cube_info,&pixel,index);
        mid.red+=(id & 1) != 0 ? bisect : -bisect;
        mid.green+=(id & 2) != 0 ? bisect : -bisect;
        mid.blue+=(id & 4) != 0 ? bisect : -bisect;
        mid.opacity+=(id & 8) != 0 ? bisect : -bisect;
        if (node_info->child[id] == (NodeInfo *) NULL)
          {
            /*
              Set colors of new node to contain pixel.
            */
            node_info->child[id]=GetNodeInfo(cube_info,id,level,node_info);
            if (node_info->child[id] == (NodeInfo *) NULL)
              {
                (void) ThrowMagickException(exception,GetMagickModule(),
                  ResourceLimitError,"MemoryAllocationFailed","`%s'",
                  image->filename);
                continue;
              }
            if (level == MaxTreeDepth)
              cube_info->colors++;
          }
        /*
          Approximate the quantization error represented by this node.
        */
        node_info=node_info->child[id];
        error.red=QuantumScale*(pixel.red-mid.red);
        error.green=QuantumScale*(pixel.green-mid.green);
        error.blue=QuantumScale*(pixel.blue-mid.blue);
        if (cube_info->associate_alpha != MagickFalse)
          error.opacity=QuantumScale*(pixel.opacity-mid.opacity);
        distance=(double) (error.red*error.red+error.green*error.green+
          error.blue*error.blue+error.opacity*error.opacity);
        if (IsNaN(distance))
          distance=0.0;
        node_info->quantize_error+=count*sqrt(distance);
        cube_info->root->quantize_error+=node_info->quantize_error;
        index--;
      }
      /*
        Sum RGB for this leaf for later derivation of the mean cube color.
      */
      node_info->number_unique+=count;
      node_info->total_color.red+=count*QuantumScale*ClampPixel(pixel.red);
      node_info->total_color.green+=count*QuantumScale*ClampPixel(pixel.green);
      node_info->total_color.blue+=count*QuantumScale*ClampPixel(pixel.blue);
      if (cube_info->associate_alpha != MagickFalse)
        node_info->total_color.opacity+=count*QuantumScale*
          ClampPixel(pixel.opacity);
      else
        node_info->total_color.opacity+=count*QuantumScale*
          ClampPixel(OpaqueOpacity);
      p+=count;
    }
    if (cube_info->colors > cube_info->maximum_colors)
      {
        PruneToCubeDepth(cube_info,cube_info->root);
        break;
      }
    proceed=SetImageProgress(image,ClassifyImageTag,(MagickOffsetType) y,
      image->rows);
    if (proceed == MagickFalse)
      break;
  }
  for (y++; y < (ssize_t) image->rows; y++)
  {
    register const PixelPacket
      *magick_restrict p;

    register ssize_t
      x;

    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
    if (p == (const PixelPacket *) NULL)
      break;
    if (cube_info->nodes > MaxNodes)
      {
        /*
          Prune one level if the color tree is too large.
        */
        PruneLevel(cube_info,cube_info->root);
        cube_info->depth--;
      }
    for (x=0; x < (ssize_t) image->columns; x+=(ssize_t) count)
    {
      /*
        Start at the root and descend the color cube tree.
      */
      for (count=1; (x+(ssize_t) count) < (ssize_t) image->columns; count++)
        if (IsSameColor(image,p,p+count) == MagickFalse)
          break;
      AssociateAlphaPixel(cube_info,p,&pixel);
      index=MaxTreeDepth-1;
      bisect=((MagickRealType) QuantumRange+1.0)/2.0;
      mid=midpoint;
      node_info=cube_info->root;
      for (level=1; level <= cube_info->depth; level++)
      {
        double
          distance;

        bisect*=0.5;
        id=ColorToNodeId(cube_info,&pixel,index);
        mid.red+=(id & 1) != 0 ? bisect : -bisect;
        mid.green+=(id & 2) != 0 ? bisect : -bisect;
        mid.blue+=(id & 4) != 0 ? bisect : -bisect;
        mid.opacity+=(id & 8) != 0 ? bisect : -bisect;
        if (node_info->child[id] == (NodeInfo *) NULL)
          {
            /*
              Set colors of new node to contain pixel.
            */
            node_info->child[id]=GetNodeInfo(cube_info,id,level,node_info);
            if (node_info->child[id] == (NodeInfo *) NULL)
              {
                (void) ThrowMagickException(exception,GetMagickModule(),
                  ResourceLimitError,"MemoryAllocationFailed","%s",
                  image->filename);
                  continue;
               }
            if (level == cube_info->depth)
              cube_info->colors++;
          }
        /*
          Approximate the quantization error represented by this node.
        */
        node_info=node_info->child[id];
        error.red=QuantumScale*(pixel.red-mid.red);
        error.green=QuantumScale*(pixel.green-mid.green);
        error.blue=QuantumScale*(pixel.blue-mid.blue);
        if (cube_info->associate_alpha != MagickFalse)
          error.opacity=QuantumScale*(pixel.opacity-mid.opacity);
        distance=(double) (error.red*error.red+error.green*error.green+
          error.blue*error.blue+error.opacity*error.opacity);
        if (IsNaN(distance))
          distance=0.0;
        node_info->quantize_error+=count*sqrt(distance);
        cube_info->root->quantize_error+=node_info->quantize_error;
        index--;
      }
      /*
        Sum RGB for this leaf for later derivation of the mean cube color.
      */
      node_info->number_unique+=count;
      node_info->total_color.red+=count*QuantumScale*ClampPixel(pixel.red);
      node_info->total_color.green+=count*QuantumScale*ClampPixel(pixel.green);
      node_info->total_color.blue+=count*QuantumScale*ClampPixel(pixel.blue);
      if (cube_info->associate_alpha != MagickFalse)
        node_info->total_color.opacity+=count*QuantumScale*ClampPixel(
          pixel.opacity);
      else
        node_info->total_color.opacity+=count*QuantumScale*
          ClampPixel(OpaqueOpacity);
      p+=count;
    }
    proceed=SetImageProgress(image,ClassifyImageTag,(MagickOffsetType) y,
      image->rows);
    if (proceed == MagickFalse)
      break;
  }
  image_view=DestroyCacheView(image_view);
  if ((cube_info->quantize_info->colorspace != UndefinedColorspace) &&
      (cube_info->quantize_info->colorspace != CMYKColorspace))
    (void) TransformImageColorspace((Image *) image,sRGBColorspace);
  return(y < (ssize_t) image->rows ? MagickFalse : MagickTrue);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   C l o n e Q u a n t i z e I n f o                                         %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  CloneQuantizeInfo() makes a duplicate of the given quantize info structure,
%  or if quantize info is NULL, a new one.
%
%  The format of the CloneQuantizeInfo method is:
%
%      QuantizeInfo *CloneQuantizeInfo(const QuantizeInfo *quantize_info)
%
%  A description of each parameter follows:
%
%    o clone_info: Method CloneQuantizeInfo returns a duplicate of the given
%      quantize info, or if image info is NULL a new one.
%
%    o quantize_info: a structure of type info.
%
*/
MagickExport QuantizeInfo *CloneQuantizeInfo(const QuantizeInfo *quantize_info)
{
  QuantizeInfo
    *clone_info;

  clone_info=(QuantizeInfo *) AcquireMagickMemory(sizeof(*clone_info));
  if (clone_info == (QuantizeInfo *) NULL)
    ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
  GetQuantizeInfo(clone_info);
  if (quantize_info == (QuantizeInfo *) NULL)
    return(clone_info);
  clone_info->number_colors=quantize_info->number_colors;
  clone_info->tree_depth=quantize_info->tree_depth;
  clone_info->dither=quantize_info->dither;
  clone_info->dither_method=quantize_info->dither_method;
  clone_info->colorspace=quantize_info->colorspace;
  clone_info->measure_error=quantize_info->measure_error;
  return(clone_info);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   C l o s e s t C o l o r                                                   %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ClosestColor() traverses the color cube tree at a particular node and
%  determines which colormap entry best represents the input color.
%
%  The format of the ClosestColor method is:
%
%      void ClosestColor(const Image *image,CubeInfo *cube_info,
%        const NodeInfo *node_info)
%
%  A description of each parameter follows.
%
%    o image: the image.
%
%    o cube_info: A pointer to the Cube structure.
%
%    o node_info: the address of a structure of type NodeInfo which points to a
%      node in the color cube tree that is to be pruned.
%
*/
static void ClosestColor(const Image *image,CubeInfo *cube_info,
  const NodeInfo *node_info)
{
  register ssize_t
    i;

  size_t
    number_children;

  /*
    Traverse any children.
  */
  number_children=cube_info->associate_alpha == MagickFalse ? 8UL : 16UL;
  for (i=0; i < (ssize_t) number_children; i++)
    if (node_info->child[i] != (NodeInfo *) NULL)
      ClosestColor(image,cube_info,node_info->child[i]);
  if (node_info->number_unique != 0)
    {
      MagickRealType
        pixel;

      register DoublePixelPacket
        *magick_restrict q;

      register MagickRealType
        alpha,
        beta,
        distance;

      register PixelPacket
        *magick_restrict p;

      /*
        Determine if this color is "closest".
      */
      p=image->colormap+node_info->color_number;
      q=(&cube_info->target);
      alpha=1.0;
      beta=1.0;
      if (cube_info->associate_alpha != MagickFalse)
        {
          alpha=(MagickRealType) (QuantumScale*GetPixelAlpha(p));
          beta=(MagickRealType) (QuantumScale*GetPixelAlpha(q));
        }
      pixel=alpha*GetPixelRed(p)-beta*GetPixelRed(q);
      distance=pixel*pixel;
      if (distance <= cube_info->distance)
        {
          pixel=alpha*GetPixelGreen(p)-beta*GetPixelGreen(q);
          distance+=pixel*pixel;
          if (distance <= cube_info->distance)
            {
              pixel=alpha*GetPixelBlue(p)-beta*GetPixelBlue(q);
              distance+=pixel*pixel;
              if (distance <= cube_info->distance)
                {
                  if (cube_info->associate_alpha != MagickFalse)
                    {
                      pixel=GetPixelAlpha(p)-GetPixelAlpha(q);
                      distance+=pixel*pixel;
                    }
                  if (distance <= cube_info->distance)
                    {
                      cube_info->distance=distance;
                      cube_info->color_number=node_info->color_number;
                    }
                }
            }
        }
    }
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   C o m p r e s s I m a g e C o l o r m a p                                 %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  CompressImageColormap() compresses an image colormap by removing any
%  duplicate or unused color entries.
%
%  The format of the CompressImageColormap method is:
%
%      MagickBooleanType CompressImageColormap(Image *image)
%
%  A description of each parameter follows:
%
%    o image: the image.
%
*/
MagickExport MagickBooleanType CompressImageColormap(Image *image)
{
  QuantizeInfo
    quantize_info;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  if (IsPaletteImage(image,&image->exception) == MagickFalse)
    return(MagickFalse);
  GetQuantizeInfo(&quantize_info);
  quantize_info.number_colors=image->colors;
  quantize_info.tree_depth=MaxTreeDepth;
  return(QuantizeImage(&quantize_info,image));
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   D e f i n e I m a g e C o l o r m a p                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  DefineImageColormap() traverses the color cube tree and notes each colormap
%  entry.  A colormap entry is any node in the color cube tree where the
%  of unique colors is not zero.  DefineImageColormap() returns the number of
%  colors in the image colormap.
%
%  The format of the DefineImageColormap method is:
%
%      size_t DefineImageColormap(Image *image,CubeInfo *cube_info,
%        NodeInfo *node_info)
%
%  A description of each parameter follows.
%
%    o image: the image.
%
%    o cube_info: A pointer to the Cube structure.
%
%    o node_info: the address of a structure of type NodeInfo which points to a
%      node in the color cube tree that is to be pruned.
%
*/
static size_t DefineImageColormap(Image *image,CubeInfo *cube_info,
  NodeInfo *node_info)
{
  register ssize_t
    i;

  size_t
    number_children;

  /*
    Traverse any children.
  */
  number_children=cube_info->associate_alpha == MagickFalse ? 8UL : 16UL;
  for (i=0; i < (ssize_t) number_children; i++)
    if (node_info->child[i] != (NodeInfo *) NULL)
      (void) DefineImageColormap(image,cube_info,node_info->child[i]);
  if (node_info->number_unique != 0)
    {
      register MagickRealType
        alpha;

      register PixelPacket
        *magick_restrict q;

      /*
        Colormap entry is defined by the mean color in this cube.
      */
      q=image->colormap+image->colors;
      alpha=(MagickRealType) ((MagickOffsetType) node_info->number_unique);
      alpha=PerceptibleReciprocal(alpha);
      if (cube_info->associate_alpha == MagickFalse)
        {
          SetPixelRed(q,ClampToQuantum((MagickRealType) (alpha*
            QuantumRange*node_info->total_color.red)));
          SetPixelGreen(q,ClampToQuantum((MagickRealType) (alpha*
            QuantumRange*node_info->total_color.green)));
          SetPixelBlue(q,ClampToQuantum((MagickRealType) (alpha*
            QuantumRange*node_info->total_color.blue)));
          SetPixelOpacity(q,OpaqueOpacity);
        }
      else
        {
          MagickRealType
            opacity;

          opacity=(MagickRealType) (alpha*QuantumRange*
            node_info->total_color.opacity);
          SetPixelOpacity(q,ClampToQuantum(opacity));
          if (q->opacity == OpaqueOpacity)
            {
              SetPixelRed(q,ClampToQuantum((MagickRealType) (alpha*
                QuantumRange*node_info->total_color.red)));
              SetPixelGreen(q,ClampToQuantum((MagickRealType) (alpha*
                QuantumRange*node_info->total_color.green)));
              SetPixelBlue(q,ClampToQuantum((MagickRealType) (alpha*
                QuantumRange*node_info->total_color.blue)));
            }
          else
            {
              double
                gamma;

              gamma=(double) (QuantumScale*(QuantumRange-(double) q->opacity));
              gamma=PerceptibleReciprocal(gamma);
              SetPixelRed(q,ClampToQuantum((MagickRealType) (alpha*
                gamma*QuantumRange*node_info->total_color.red)));
              SetPixelGreen(q,ClampToQuantum((MagickRealType) (alpha*
                gamma*QuantumRange*node_info->total_color.green)));
              SetPixelBlue(q,ClampToQuantum((MagickRealType) (alpha*
                gamma*QuantumRange*node_info->total_color.blue)));
              if (node_info->number_unique > cube_info->transparent_pixels)
                {
                  cube_info->transparent_pixels=node_info->number_unique;
                  cube_info->transparent_index=(ssize_t) image->colors;
                }
            }
        }
      node_info->color_number=image->colors++;
    }
  return(image->colors);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   D e s t r o y C u b e I n f o                                             %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  DestroyCubeInfo() deallocates memory associated with an image.
%
%  The format of the DestroyCubeInfo method is:
%
%      DestroyCubeInfo(CubeInfo *cube_info)
%
%  A description of each parameter follows:
%
%    o cube_info: the address of a structure of type CubeInfo.
%
*/
static void DestroyCubeInfo(CubeInfo *cube_info)
{
  register Nodes
    *nodes;

  /*
    Release color cube tree storage.
  */
  do
  {
    nodes=cube_info->node_queue->next;
    cube_info->node_queue->nodes=(NodeInfo *) RelinquishMagickMemory(
      cube_info->node_queue->nodes);
    cube_info->node_queue=(Nodes *) RelinquishMagickMemory(
      cube_info->node_queue);
    cube_info->node_queue=nodes;
  } while (cube_info->node_queue != (Nodes *) NULL);
  if (cube_info->memory_info != (MemoryInfo *) NULL)
    cube_info->memory_info=RelinquishVirtualMemory(cube_info->memory_info);
  cube_info->quantize_info=DestroyQuantizeInfo(cube_info->quantize_info);
  cube_info=(CubeInfo *) RelinquishMagickMemory(cube_info);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   D e s t r o y Q u a n t i z e I n f o                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  DestroyQuantizeInfo() deallocates memory associated with an QuantizeInfo
%  structure.
%
%  The format of the DestroyQuantizeInfo method is:
%
%      QuantizeInfo *DestroyQuantizeInfo(QuantizeInfo *quantize_info)
%
%  A description of each parameter follows:
%
%    o quantize_info: Specifies a pointer to an QuantizeInfo structure.
%
*/
MagickExport QuantizeInfo *DestroyQuantizeInfo(QuantizeInfo *quantize_info)
{
  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
  assert(quantize_info != (QuantizeInfo *) NULL);
  assert(quantize_info->signature == MagickSignature);
  quantize_info->signature=(~MagickSignature);
  quantize_info=(QuantizeInfo *) RelinquishMagickMemory(quantize_info);
  return(quantize_info);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   D i t h e r I m a g e                                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  DitherImage() distributes the difference between an original image and
%  the corresponding color reduced algorithm to neighboring pixels using
%  serpentine-scan Floyd-Steinberg error diffusion. DitherImage returns
%  MagickTrue if the image is dithered otherwise MagickFalse.
%
%  The format of the DitherImage method is:
%
%      MagickBooleanType DitherImage(Image *image,CubeInfo *cube_info)
%
%  A description of each parameter follows.
%
%    o image: the image.
%
%    o cube_info: A pointer to the Cube structure.
%
*/

static DoublePixelPacket **DestroyPixelThreadSet(DoublePixelPacket **pixels)
{
  register ssize_t
    i;

  assert(pixels != (DoublePixelPacket **) NULL);
  for (i=0; i < (ssize_t) GetMagickResourceLimit(ThreadResource); i++)
    if (pixels[i] != (DoublePixelPacket *) NULL)
      pixels[i]=(DoublePixelPacket *) RelinquishMagickMemory(pixels[i]);
  pixels=(DoublePixelPacket **) RelinquishMagickMemory(pixels);
  return(pixels);
}

static DoublePixelPacket **AcquirePixelThreadSet(const size_t count)
{
  DoublePixelPacket
    **pixels;

  register ssize_t
    i;

  size_t
    number_threads;

  number_threads=(size_t) GetMagickResourceLimit(ThreadResource);
  pixels=(DoublePixelPacket **) AcquireQuantumMemory(number_threads,
    sizeof(*pixels));
  if (pixels == (DoublePixelPacket **) NULL)
    return((DoublePixelPacket **) NULL);
  (void) ResetMagickMemory(pixels,0,number_threads*sizeof(*pixels));
  for (i=0; i < (ssize_t) number_threads; i++)
  {
    pixels[i]=(DoublePixelPacket *) AcquireQuantumMemory(count,
      2*sizeof(**pixels));
    if (pixels[i] == (DoublePixelPacket *) NULL)
      return(DestroyPixelThreadSet(pixels));
  }
  return(pixels);
}

static inline ssize_t CacheOffset(CubeInfo *cube_info,
  const DoublePixelPacket *pixel)
{
#define RedShift(pixel) (((pixel) >> CacheShift) << (0*(8-CacheShift)))
#define GreenShift(pixel) (((pixel) >> CacheShift) << (1*(8-CacheShift)))
#define BlueShift(pixel) (((pixel) >> CacheShift) << (2*(8-CacheShift)))
#define AlphaShift(pixel) (((pixel) >> CacheShift) << (3*(8-CacheShift)))

  ssize_t
    offset;

  offset=(ssize_t) (RedShift(ScaleQuantumToChar(ClampPixel(pixel->red))) |
    GreenShift(ScaleQuantumToChar(ClampPixel(pixel->green))) |
    BlueShift(ScaleQuantumToChar(ClampPixel(pixel->blue))));
  if (cube_info->associate_alpha != MagickFalse)
    offset|=AlphaShift(ScaleQuantumToChar(ClampPixel(pixel->opacity)));
  return(offset);
}

static MagickBooleanType FloydSteinbergDither(Image *image,CubeInfo *cube_info)
{
#define DitherImageTag  "Dither/Image"

  CacheView
    *image_view;

  DoublePixelPacket
    **pixels;

  ExceptionInfo
    *exception;

  MagickBooleanType
    status;

  ssize_t
    y;

  /*
    Distribute quantization error using Floyd-Steinberg.
  */
  pixels=AcquirePixelThreadSet(image->columns);
  if (pixels == (DoublePixelPacket **) NULL)
    return(MagickFalse);
  exception=(&image->exception);
  status=MagickTrue;
  image_view=AcquireAuthenticCacheView(image,exception);
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    const int
      id = GetOpenMPThreadId();

    CubeInfo
      cube;

    DoublePixelPacket
      *current,
      *previous;

    register IndexPacket
      *magick_restrict indexes;

    register PixelPacket
      *magick_restrict q;

    register ssize_t
      x;

    size_t
      index;

    ssize_t
      v;

    if (status == MagickFalse)
      continue;
    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
    if (q == (PixelPacket *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    indexes=GetCacheViewAuthenticIndexQueue(image_view);
    cube=(*cube_info);
    current=pixels[id]+(y & 0x01)*image->columns;
    previous=pixels[id]+((y+1) & 0x01)*image->columns;
    v=(ssize_t) ((y & 0x01) ? -1 : 1);
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      DoublePixelPacket
        color,
        pixel;

      register ssize_t
        i;

      ssize_t
        u;

      u=(y & 0x01) ? (ssize_t) image->columns-1-x : x;
      AssociateAlphaPixel(&cube,q+u,&pixel);
      if (x > 0)
        {
          pixel.red+=7*current[u-v].red/16;
          pixel.green+=7*current[u-v].green/16;
          pixel.blue+=7*current[u-v].blue/16;
          if (cube.associate_alpha != MagickFalse)
            pixel.opacity+=7*current[u-v].opacity/16;
        }
      if (y > 0)
        {
          if (x < (ssize_t) (image->columns-1))
            {
              pixel.red+=previous[u+v].red/16;
              pixel.green+=previous[u+v].green/16;
              pixel.blue+=previous[u+v].blue/16;
              if (cube.associate_alpha != MagickFalse)
                pixel.opacity+=previous[u+v].opacity/16;
            }
          pixel.red+=5*previous[u].red/16;
          pixel.green+=5*previous[u].green/16;
          pixel.blue+=5*previous[u].blue/16;
          if (cube.associate_alpha != MagickFalse)
            pixel.opacity+=5*previous[u].opacity/16;
          if (x > 0)
            {
              pixel.red+=3*previous[u-v].red/16;
              pixel.green+=3*previous[u-v].green/16;
              pixel.blue+=3*previous[u-v].blue/16;
              if (cube.associate_alpha != MagickFalse)
                pixel.opacity+=3*previous[u-v].opacity/16;
            }
        }
      pixel.red=(MagickRealType) ClampPixel(pixel.red);
      pixel.green=(MagickRealType) ClampPixel(pixel.green);
      pixel.blue=(MagickRealType) ClampPixel(pixel.blue);
      if (cube.associate_alpha != MagickFalse)
        pixel.opacity=(MagickRealType) ClampPixel(pixel.opacity);
      i=CacheOffset(&cube,&pixel);
      if (cube.cache[i] < 0)
        {
          register NodeInfo
            *node_info;

          register size_t
            id;

          /*
            Identify the deepest node containing the pixel's color.
          */
          node_info=cube.root;
          for (index=MaxTreeDepth-1; (ssize_t) index > 0; index--)
          {
            id=ColorToNodeId(&cube,&pixel,index);
            if (node_info->child[id] == (NodeInfo *) NULL)
              break;
            node_info=node_info->child[id];
          }
          /*
            Find closest color among siblings and their children.
          */
          cube.target=pixel;
          cube.distance=(MagickRealType) (4.0*(QuantumRange+1.0)*(QuantumRange+
            1.0)+1.0);
          ClosestColor(image,&cube,node_info->parent);
          cube.cache[i]=(ssize_t) cube.color_number;
        }
      /*
        Assign pixel to closest colormap entry.
      */
      index=(size_t) cube.cache[i];
      if (image->storage_class == PseudoClass)
        SetPixelIndex(indexes+u,index);
      if (cube.quantize_info->measure_error == MagickFalse)
        {
          SetPixelRgb(q+u,image->colormap+index);
          if (cube.associate_alpha != MagickFalse)
            SetPixelOpacity(q+u,image->colormap[index].opacity);
        }
      if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
        status=MagickFalse;
      /*
        Store the error.
      */
      AssociateAlphaPixel(&cube,image->colormap+index,&color);
      current[u].red=pixel.red-color.red;
      current[u].green=pixel.green-color.green;
      current[u].blue=pixel.blue-color.blue;
      if (cube.associate_alpha != MagickFalse)
        current[u].opacity=pixel.opacity-color.opacity;
      if (image->progress_monitor != (MagickProgressMonitor) NULL)
        {
          MagickBooleanType
            proceed;

          proceed=SetImageProgress(image,DitherImageTag,(MagickOffsetType) y,
            image->rows);
          if (proceed == MagickFalse)
            status=MagickFalse;
        }
    }
  }
  image_view=DestroyCacheView(image_view);
  pixels=DestroyPixelThreadSet(pixels);
  return(MagickTrue);
}

static MagickBooleanType
  RiemersmaDither(Image *,CacheView *,CubeInfo *,const unsigned int);

static void Riemersma(Image *image,CacheView *image_view,CubeInfo *cube_info,
  const size_t level,const unsigned int direction)
{
  if (level == 1)
    switch (direction)
    {
      case WestGravity:
      {
        (void) RiemersmaDither(image,image_view,cube_info,EastGravity);
        (void) RiemersmaDither(image,image_view,cube_info,SouthGravity);
        (void) RiemersmaDither(image,image_view,cube_info,WestGravity);
        break;
      }
      case EastGravity:
      {
        (void) RiemersmaDither(image,image_view,cube_info,WestGravity);
        (void) RiemersmaDither(image,image_view,cube_info,NorthGravity);
        (void) RiemersmaDither(image,image_view,cube_info,EastGravity);
        break;
      }
      case NorthGravity:
      {
        (void) RiemersmaDither(image,image_view,cube_info,SouthGravity);
        (void) RiemersmaDither(image,image_view,cube_info,EastGravity);
        (void) RiemersmaDither(image,image_view,cube_info,NorthGravity);
        break;
      }
      case SouthGravity:
      {
        (void) RiemersmaDither(image,image_view,cube_info,NorthGravity);
        (void) RiemersmaDither(image,image_view,cube_info,WestGravity);
        (void) RiemersmaDither(image,image_view,cube_info,SouthGravity);
        break;
      }
      default:
        break;
    }
  else
    switch (direction)
    {
      case WestGravity:
      {
        Riemersma(image,image_view,cube_info,level-1,NorthGravity);
        (void) RiemersmaDither(image,image_view,cube_info,EastGravity);
        Riemersma(image,image_view,cube_info,level-1,WestGravity);
        (void) RiemersmaDither(image,image_view,cube_info,SouthGravity);
        Riemersma(image,image_view,cube_info,level-1,WestGravity);
        (void) RiemersmaDither(image,image_view,cube_info,WestGravity);
        Riemersma(image,image_view,cube_info,level-1,SouthGravity);
        break;
      }
      case EastGravity:
      {
        Riemersma(image,image_view,cube_info,level-1,SouthGravity);
        (void) RiemersmaDither(image,image_view,cube_info,WestGravity);
        Riemersma(image,image_view,cube_info,level-1,EastGravity);
        (void) RiemersmaDither(image,image_view,cube_info,NorthGravity);
        Riemersma(image,image_view,cube_info,level-1,EastGravity);
        (void) RiemersmaDither(image,image_view,cube_info,EastGravity);
        Riemersma(image,image_view,cube_info,level-1,NorthGravity);
        break;
      }
      case NorthGravity:
      {
        Riemersma(image,image_view,cube_info,level-1,WestGravity);
        (void) RiemersmaDither(image,image_view,cube_info,SouthGravity);
        Riemersma(image,image_view,cube_info,level-1,NorthGravity);
        (void) RiemersmaDither(image,image_view,cube_info,EastGravity);
        Riemersma(image,image_view,cube_info,level-1,NorthGravity);
        (void) RiemersmaDither(image,image_view,cube_info,NorthGravity);
        Riemersma(image,image_view,cube_info,level-1,EastGravity);
        break;
      }
      case SouthGravity:
      {
        Riemersma(image,image_view,cube_info,level-1,EastGravity);
        (void) RiemersmaDither(image,image_view,cube_info,NorthGravity);
        Riemersma(image,image_view,cube_info,level-1,SouthGravity);
        (void) RiemersmaDither(image,image_view,cube_info,WestGravity);
        Riemersma(image,image_view,cube_info,level-1,SouthGravity);
        (void) RiemersmaDither(image,image_view,cube_info,SouthGravity);
        Riemersma(image,image_view,cube_info,level-1,WestGravity);
        break;
      }
      default:
        break;
    }
}

static MagickBooleanType RiemersmaDither(Image *image,CacheView *image_view,
  CubeInfo *cube_info,const unsigned int direction)
{
#define DitherImageTag  "Dither/Image"

  DoublePixelPacket
    color,
    pixel;

  MagickBooleanType
    proceed;

  register CubeInfo
    *p;

  size_t
    index;

  p=cube_info;
  if ((p->x >= 0) && (p->x < (ssize_t) image->columns) &&
      (p->y >= 0) && (p->y < (ssize_t) image->rows))
    {
      ExceptionInfo
        *exception;

      register IndexPacket
        *magick_restrict indexes;

      register PixelPacket
        *magick_restrict q;

      register ssize_t
        i;

      /*
        Distribute error.
      */
      exception=(&image->exception);
      q=GetCacheViewAuthenticPixels(image_view,p->x,p->y,1,1,exception);
      if (q == (PixelPacket *) NULL)
        return(MagickFalse);
      indexes=GetCacheViewAuthenticIndexQueue(image_view);
      AssociateAlphaPixel(cube_info,q,&pixel);
      for (i=0; i < ErrorQueueLength; i++)
      {
        pixel.red+=p->weights[i]*p->error[i].red;
        pixel.green+=p->weights[i]*p->error[i].green;
        pixel.blue+=p->weights[i]*p->error[i].blue;
        if (cube_info->associate_alpha != MagickFalse)
          pixel.opacity+=p->weights[i]*p->error[i].opacity;
      }
      pixel.red=(MagickRealType) ClampPixel(pixel.red);
      pixel.green=(MagickRealType) ClampPixel(pixel.green);
      pixel.blue=(MagickRealType) ClampPixel(pixel.blue);
      if (cube_info->associate_alpha != MagickFalse)
        pixel.opacity=(MagickRealType) ClampPixel(pixel.opacity);
      i=CacheOffset(cube_info,&pixel);
      if (p->cache[i] < 0)
        {
          register NodeInfo
            *node_info;

          register size_t
            id;

          /*
            Identify the deepest node containing the pixel's color.
          */
          node_info=p->root;
          for (index=MaxTreeDepth-1; (ssize_t) index > 0; index--)
          {
            id=ColorToNodeId(cube_info,&pixel,index);
            if (node_info->child[id] == (NodeInfo *) NULL)
              break;
            node_info=node_info->child[id];
          }
          /*
            Find closest color among siblings and their children.
          */
          p->target=pixel;
          p->distance=(MagickRealType) (4.0*(QuantumRange+1.0)*((MagickRealType)
            QuantumRange+1.0)+1.0);
          ClosestColor(image,p,node_info->parent);
          p->cache[i]=(ssize_t) p->color_number;
        }
      /*
        Assign pixel to closest colormap entry.
      */
      index=(size_t) (1*p->cache[i]);
      if (image->storage_class == PseudoClass)
        *indexes=(IndexPacket) index;
      if (cube_info->quantize_info->measure_error == MagickFalse)
        {
          SetPixelRgb(q,image->colormap+index);
          if (cube_info->associate_alpha != MagickFalse)
            SetPixelOpacity(q,image->colormap[index].opacity);
        }
      if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
        return(MagickFalse);
      /*
        Propagate the error as the last entry of the error queue.
      */
      (void) CopyMagickMemory(p->error,p->error+1,(ErrorQueueLength-1)*
        sizeof(p->error[0]));
      AssociateAlphaPixel(cube_info,image->colormap+index,&color);
      p->error[ErrorQueueLength-1].red=pixel.red-color.red;
      p->error[ErrorQueueLength-1].green=pixel.green-color.green;
      p->error[ErrorQueueLength-1].blue=pixel.blue-color.blue;
      if (cube_info->associate_alpha != MagickFalse)
        p->error[ErrorQueueLength-1].opacity=pixel.opacity-color.opacity;
      proceed=SetImageProgress(image,DitherImageTag,p->offset,p->span);
      if (proceed == MagickFalse)
        return(MagickFalse);
      p->offset++;
    }
  switch (direction)
  {
    case WestGravity: p->x--; break;
    case EastGravity: p->x++; break;
    case NorthGravity: p->y--; break;
    case SouthGravity: p->y++; break;
  }
  return(MagickTrue);
}

static MagickBooleanType DitherImage(Image *image,CubeInfo *cube_info)
{
  CacheView
    *image_view;

  MagickBooleanType
    status;

  register ssize_t
    i;

  size_t
    depth;

  if (cube_info->quantize_info->dither_method != RiemersmaDitherMethod)
    return(FloydSteinbergDither(image,cube_info));
  /*
    Distribute quantization error along a Hilbert curve.
  */
  (void) ResetMagickMemory(cube_info->error,0,ErrorQueueLength*
    sizeof(*cube_info->error));
  cube_info->x=0;
  cube_info->y=0;
  i=MagickMax((ssize_t) image->columns,(ssize_t) image->rows);
  for (depth=1; i != 0; depth++)
    i>>=1;
  if ((ssize_t) (1L << depth) < MagickMax((ssize_t) image->columns,(ssize_t) image->rows))
    depth++;
  cube_info->offset=0;
  cube_info->span=(MagickSizeType) image->columns*image->rows;
  image_view=AcquireAuthenticCacheView(image,&image->exception);
  if (depth > 1)
    Riemersma(image,image_view,cube_info,depth-1,NorthGravity);
  status=RiemersmaDither(image,image_view,cube_info,ForgetGravity);
  image_view=DestroyCacheView(image_view);
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   G e t C u b e I n f o                                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GetCubeInfo() initialize the Cube data structure.
%
%  The format of the GetCubeInfo method is:
%
%      CubeInfo GetCubeInfo(const QuantizeInfo *quantize_info,
%        const size_t depth,const size_t maximum_colors)
%
%  A description of each parameter follows.
%
%    o quantize_info: Specifies a pointer to an QuantizeInfo structure.
%
%    o depth: Normally, this integer value is zero or one.  A zero or
%      one tells Quantize to choose a optimal tree depth of Log4(number_colors).
%      A tree of this depth generally allows the best representation of the
%      reference image with the least amount of memory and the fastest
%      computational speed.  In some cases, such as an image with low color
%      dispersion (a few number of colors), a value other than
%      Log4(number_colors) is required.  To expand the color tree completely,
%      use a value of 8.
%
%    o maximum_colors: maximum colors.
%
*/
static CubeInfo *GetCubeInfo(const QuantizeInfo *quantize_info,
  const size_t depth,const size_t maximum_colors)
{
  CubeInfo
    *cube_info;

  MagickRealType
    sum,
    weight;

  register ssize_t
    i;

  size_t
    length;

  /*
    Initialize tree to describe color cube_info.
  */
  cube_info=(CubeInfo *) AcquireMagickMemory(sizeof(*cube_info));
  if (cube_info == (CubeInfo *) NULL)
    return((CubeInfo *) NULL);
  (void) ResetMagickMemory(cube_info,0,sizeof(*cube_info));
  cube_info->depth=depth;
  if (cube_info->depth > MaxTreeDepth)
    cube_info->depth=MaxTreeDepth;
  if (cube_info->depth < 2)
    cube_info->depth=2;
  cube_info->maximum_colors=maximum_colors;
  /*
    Initialize root node.
  */
  cube_info->root=GetNodeInfo(cube_info,0,0,(NodeInfo *) NULL);
  if (cube_info->root == (NodeInfo *) NULL)
    return((CubeInfo *) NULL);
  cube_info->root->parent=cube_info->root;
  cube_info->quantize_info=CloneQuantizeInfo(quantize_info);
  if (cube_info->quantize_info->dither == MagickFalse)
    return(cube_info);
  /*
    Initialize dither resources.
  */
  length=(size_t) (1UL << (4*(8-CacheShift)));
  cube_info->memory_info=AcquireVirtualMemory(length,sizeof(*cube_info->cache));
  if (cube_info->memory_info == (MemoryInfo *) NULL)
    return((CubeInfo *) NULL);
  cube_info->cache=(ssize_t *) GetVirtualMemoryBlob(cube_info->memory_info);
  /*
    Initialize color cache.
  */
  (void) ResetMagickMemory(cube_info->cache,(-1),sizeof(*cube_info->cache)*
    length);
  /*
    Distribute weights along a curve of exponential decay.
  */
  weight=1.0;
  for (i=0; i < ErrorQueueLength; i++)
  {
    cube_info->weights[ErrorQueueLength-i-1]=PerceptibleReciprocal(weight);
    weight*=exp(log(((double) QuantumRange+1.0))/(ErrorQueueLength-1.0));
  }
  /*
    Normalize the weighting factors.
  */
  weight=0.0;
  for (i=0; i < ErrorQueueLength; i++)
    weight+=cube_info->weights[i];
  sum=0.0;
  for (i=0; i < ErrorQueueLength; i++)
  {
    cube_info->weights[i]/=weight;
    sum+=cube_info->weights[i];
  }
  cube_info->weights[0]+=1.0-sum;
  return(cube_info);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   G e t N o d e I n f o                                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GetNodeInfo() allocates memory for a new node in the color cube tree and
%  presets all fields to zero.
%
%  The format of the GetNodeInfo method is:
%
%      NodeInfo *GetNodeInfo(CubeInfo *cube_info,const size_t id,
%        const size_t level,NodeInfo *parent)
%
%  A description of each parameter follows.
%
%    o node: The GetNodeInfo method returns a pointer to a queue of nodes.
%
%    o id: Specifies the child number of the node.
%
%    o level: Specifies the level in the storage_class the node resides.
%
*/
static NodeInfo *GetNodeInfo(CubeInfo *cube_info,const size_t id,
  const size_t level,NodeInfo *parent)
{
  NodeInfo
    *node_info;

  if (cube_info->free_nodes == 0)
    {
      Nodes
        *nodes;

      /*
        Allocate a new queue of nodes.
      */
      nodes=(Nodes *) AcquireMagickMemory(sizeof(*nodes));
      if (nodes == (Nodes *) NULL)
        return((NodeInfo *) NULL);
      nodes->nodes=(NodeInfo *) AcquireQuantumMemory(NodesInAList,
        sizeof(*nodes->nodes));
      if (nodes->nodes == (NodeInfo *) NULL)
        return((NodeInfo *) NULL);
      nodes->next=cube_info->node_queue;
      cube_info->node_queue=nodes;
      cube_info->next_node=nodes->nodes;
      cube_info->free_nodes=NodesInAList;
    }
  cube_info->nodes++;
  cube_info->free_nodes--;
  node_info=cube_info->next_node++;
  (void) ResetMagickMemory(node_info,0,sizeof(*node_info));
  node_info->parent=parent;
  node_info->id=id;
  node_info->level=level;
  return(node_info);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%  G e t I m a g e Q u a n t i z e E r r o r                                  %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GetImageQuantizeError() measures the difference between the original
%  and quantized images.  This difference is the total quantization error.
%  The error is computed by summing over all pixels in an image the distance
%  squared in RGB space between each reference pixel value and its quantized
%  value.  These values are computed:
%
%    o mean_error_per_pixel:  This value is the mean error for any single
%      pixel in the image.
%
%    o normalized_mean_square_error:  This value is the normalized mean
%      quantization error for any single pixel in the image.  This distance
%      measure is normalized to a range between 0 and 1.  It is independent
%      of the range of red, green, and blue values in the image.
%
%    o normalized_maximum_square_error:  Thsi value is the normalized
%      maximum quantization error for any single pixel in the image.  This
%      distance measure is normalized to a range between 0 and 1.  It is
%      independent of the range of red, green, and blue values in your image.
%
%  The format of the GetImageQuantizeError method is:
%
%      MagickBooleanType GetImageQuantizeError(Image *image)
%
%  A description of each parameter follows.
%
%    o image: the image.
%
*/
MagickExport MagickBooleanType GetImageQuantizeError(Image *image)
{
  CacheView
    *image_view;

  ExceptionInfo
    *exception;

  IndexPacket
    *indexes;

  MagickRealType
    alpha,
    area,
    beta,
    distance,
    maximum_error,
    mean_error,
    mean_error_per_pixel;

  size_t
    index;

  ssize_t
    y;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  image->total_colors=GetNumberColors(image,(FILE *) NULL,&image->exception);
  (void) ResetMagickMemory(&image->error,0,sizeof(image->error));
  if (image->storage_class == DirectClass)
    return(MagickTrue);
  alpha=1.0;
  beta=1.0;
  area=3.0*image->columns*image->rows;
  maximum_error=0.0;
  mean_error_per_pixel=0.0;
  mean_error=0.0;
  exception=(&image->exception);
  image_view=AcquireVirtualCacheView(image,exception);
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    register const PixelPacket
      *magick_restrict p;

    register ssize_t
      x;

    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
    if (p == (const PixelPacket *) NULL)
      break;
    indexes=GetCacheViewAuthenticIndexQueue(image_view);
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      index=1UL*GetPixelIndex(indexes+x);
      if (image->matte != MagickFalse)
        {
          alpha=(MagickRealType) (QuantumScale*(GetPixelAlpha(p)));
          beta=(MagickRealType) (QuantumScale*(QuantumRange-
            image->colormap[index].opacity));
        }
      distance=fabs((double) (alpha*GetPixelRed(p)-beta*
        image->colormap[index].red));
      mean_error_per_pixel+=distance;
      mean_error+=distance*distance;
      if (distance > maximum_error)
        maximum_error=distance;
      distance=fabs((double) (alpha*GetPixelGreen(p)-beta*
        image->colormap[index].green));
      mean_error_per_pixel+=distance;
      mean_error+=distance*distance;
      if (distance > maximum_error)
        maximum_error=distance;
      distance=fabs((double) (alpha*GetPixelBlue(p)-beta*
        image->colormap[index].blue));
      mean_error_per_pixel+=distance;
      mean_error+=distance*distance;
      if (distance > maximum_error)
        maximum_error=distance;
      p++;
    }
  }
  image_view=DestroyCacheView(image_view);
  image->error.mean_error_per_pixel=(double) mean_error_per_pixel/area;
  image->error.normalized_mean_error=(double) QuantumScale*QuantumScale*
    mean_error/area;
  image->error.normalized_maximum_error=(double) QuantumScale*maximum_error;
  return(MagickTrue);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   G e t Q u a n t i z e I n f o                                             %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  GetQuantizeInfo() initializes the QuantizeInfo structure.
%
%  The format of the GetQuantizeInfo method is:
%
%      GetQuantizeInfo(QuantizeInfo *quantize_info)
%
%  A description of each parameter follows:
%
%    o quantize_info: Specifies a pointer to a QuantizeInfo structure.
%
*/
MagickExport void GetQuantizeInfo(QuantizeInfo *quantize_info)
{
  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"...");
  assert(quantize_info != (QuantizeInfo *) NULL);
  (void) ResetMagickMemory(quantize_info,0,sizeof(*quantize_info));
  quantize_info->number_colors=256;
  quantize_info->dither=MagickTrue;
  quantize_info->dither_method=RiemersmaDitherMethod;
  quantize_info->colorspace=UndefinedColorspace;
  quantize_info->measure_error=MagickFalse;
  quantize_info->signature=MagickSignature;
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%     P o s t e r i z e I m a g e C h a n n e l                               %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  PosterizeImage() reduces the image to a limited number of colors for a
%  "poster" effect.
%
%  The format of the PosterizeImage method is:
%
%      MagickBooleanType PosterizeImage(Image *image,const size_t levels,
%        const MagickBooleanType dither)
%      MagickBooleanType PosterizeImageChannel(Image *image,
%        const ChannelType channel,const size_t levels,
%        const MagickBooleanType dither)
%
%  A description of each parameter follows:
%
%    o image: Specifies a pointer to an Image structure.
%
%    o levels: Number of color levels allowed in each channel.  Very low values
%      (2, 3, or 4) have the most visible effect.
%
%    o dither: Set this integer value to something other than zero to dither
%      the mapped image.
%
*/

static inline double MagickRound(double x)
{
  /*
    Round the fraction to nearest integer.
  */
  if ((x-floor(x)) < (ceil(x)-x))
    return(floor(x));
  return(ceil(x));
}

MagickExport MagickBooleanType PosterizeImage(Image *image,const size_t levels,
  const MagickBooleanType dither)
{
  MagickBooleanType
    status;

  status=PosterizeImageChannel(image,DefaultChannels,levels,dither);
  return(status);
}

MagickExport MagickBooleanType PosterizeImageChannel(Image *image,
  const ChannelType channel,const size_t levels,const MagickBooleanType dither)
{
#define PosterizeImageTag  "Posterize/Image"
#define PosterizePixel(pixel) (Quantum) (QuantumRange*(MagickRound( \
  QuantumScale*pixel*(levels-1)))/MagickMax((ssize_t) levels-1,1))

  CacheView
    *image_view;

  ExceptionInfo
    *exception;

  MagickBooleanType
    status;

  MagickOffsetType
    progress;

  QuantizeInfo
    *quantize_info;

  register ssize_t
    i;

  ssize_t
    y;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  if (image->storage_class == PseudoClass)
#if defined(MAGICKCORE_OPENMP_SUPPORT)
    #pragma omp parallel for schedule(static,4) shared(progress,status) \
      magick_threads(image,image,1,1)
#endif
    for (i=0; i < (ssize_t) image->colors; i++)
    {
      /*
        Posterize colormap.
      */
      if ((channel & RedChannel) != 0)
        image->colormap[i].red=PosterizePixel(image->colormap[i].red);
      if ((channel & GreenChannel) != 0)
        image->colormap[i].green=PosterizePixel(image->colormap[i].green);
      if ((channel & BlueChannel) != 0)
        image->colormap[i].blue=PosterizePixel(image->colormap[i].blue);
      if ((channel & OpacityChannel) != 0)
        image->colormap[i].opacity=PosterizePixel(image->colormap[i].opacity);
    }
  /*
    Posterize image.
  */
  status=MagickTrue;
  progress=0;
  exception=(&image->exception);
  image_view=AcquireAuthenticCacheView(image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static,4) shared(progress,status) \
    magick_threads(image,image,image->rows,1)
#endif
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    register IndexPacket
      *magick_restrict indexes;

    register PixelPacket
      *magick_restrict q;

    register ssize_t
      x;

    if (status == MagickFalse)
      continue;
    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
    if (q == (PixelPacket *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    indexes=GetCacheViewAuthenticIndexQueue(image_view);
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      if ((channel & RedChannel) != 0)
        SetPixelRed(q,PosterizePixel(GetPixelRed(q)));
      if ((channel & GreenChannel) != 0)
        SetPixelGreen(q,PosterizePixel(GetPixelGreen(q)));
      if ((channel & BlueChannel) != 0)
        SetPixelBlue(q,PosterizePixel(GetPixelBlue(q)));
      if (((channel & OpacityChannel) != 0) &&
          (image->matte != MagickFalse))
        SetPixelOpacity(q,PosterizePixel(GetPixelOpacity(q)));
      if (((channel & IndexChannel) != 0) &&
          (image->colorspace == CMYKColorspace))
        SetPixelIndex(indexes+x,PosterizePixel(GetPixelIndex(indexes+x)));
      q++;
    }
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      status=MagickFalse;
    if (image->progress_monitor != (MagickProgressMonitor) NULL)
      {
        MagickBooleanType
          proceed;

#if defined(MAGICKCORE_OPENMP_SUPPORT)
        #pragma omp critical (MagickCore_PosterizeImageChannel)
#endif
        proceed=SetImageProgress(image,PosterizeImageTag,progress++,
          image->rows);
        if (proceed == MagickFalse)
          status=MagickFalse;
      }
  }
  image_view=DestroyCacheView(image_view);
  quantize_info=AcquireQuantizeInfo((ImageInfo *) NULL);
  quantize_info->number_colors=(size_t) MagickMin((ssize_t) levels*levels*
    levels,MaxColormapSize+1);
  quantize_info->dither=dither;
  quantize_info->tree_depth=MaxTreeDepth;
  status=QuantizeImage(quantize_info,image);
  quantize_info=DestroyQuantizeInfo(quantize_info);
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   P r u n e C h i l d                                                       %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  PruneChild() deletes the given node and merges its statistics into its
%  parent.
%
%  The format of the PruneSubtree method is:
%
%      PruneChild(CubeInfo *cube_info,const NodeInfo *node_info)
%
%  A description of each parameter follows.
%
%    o cube_info: A pointer to the Cube structure.
%
%    o node_info: pointer to node in color cube tree that is to be pruned.
%
*/
static void PruneChild(CubeInfo *cube_info,const NodeInfo *node_info)
{
  NodeInfo
    *parent;

  register ssize_t
    i;

  size_t
    number_children;

  /*
    Traverse any children.
  */
  number_children=cube_info->associate_alpha == MagickFalse ? 8UL : 16UL;
  for (i=0; i < (ssize_t) number_children; i++)
    if (node_info->child[i] != (NodeInfo *) NULL)
      PruneChild(cube_info,node_info->child[i]);
  /*
    Merge color statistics into parent.
  */
  parent=node_info->parent;
  parent->number_unique+=node_info->number_unique;
  parent->total_color.red+=node_info->total_color.red;
  parent->total_color.green+=node_info->total_color.green;
  parent->total_color.blue+=node_info->total_color.blue;
  parent->total_color.opacity+=node_info->total_color.opacity;
  parent->child[node_info->id]=(NodeInfo *) NULL;
  cube_info->nodes--;
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+  P r u n e L e v e l                                                        %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  PruneLevel() deletes all nodes at the bottom level of the color tree merging
%  their color statistics into their parent node.
%
%  The format of the PruneLevel method is:
%
%      PruneLevel(CubeInfo *cube_info,const NodeInfo *node_info)
%
%  A description of each parameter follows.
%
%    o cube_info: A pointer to the Cube structure.
%
%    o node_info: pointer to node in color cube tree that is to be pruned.
%
*/
static void PruneLevel(CubeInfo *cube_info,const NodeInfo *node_info)
{
  register ssize_t
    i;

  size_t
    number_children;

  /*
    Traverse any children.
  */
  number_children=cube_info->associate_alpha == MagickFalse ? 8UL : 16UL;
  for (i=0; i < (ssize_t) number_children; i++)
    if (node_info->child[i] != (NodeInfo *) NULL)
      PruneLevel(cube_info,node_info->child[i]);
  if (node_info->level == cube_info->depth)
    PruneChild(cube_info,node_info);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+  P r u n e T o C u b e D e p t h                                            %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  PruneToCubeDepth() deletes any nodes at a depth greater than
%  cube_info->depth while merging their color statistics into their parent
%  node.
%
%  The format of the PruneToCubeDepth method is:
%
%      PruneToCubeDepth(CubeInfo *cube_info,const NodeInfo *node_info)
%
%  A description of each parameter follows.
%
%    o cube_info: A pointer to the Cube structure.
%
%    o node_info: pointer to node in color cube tree that is to be pruned.
%
*/
static void PruneToCubeDepth(CubeInfo *cube_info,const NodeInfo *node_info)
{
  register ssize_t
    i;

  size_t
    number_children;

  /*
    Traverse any children.
  */
  number_children=cube_info->associate_alpha == MagickFalse ? 8UL : 16UL;
  for (i=0; i < (ssize_t) number_children; i++)
    if (node_info->child[i] != (NodeInfo *) NULL)
      PruneToCubeDepth(cube_info,node_info->child[i]);
  if (node_info->level > cube_info->depth)
    PruneChild(cube_info,node_info);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%  Q u a n t i z e I m a g e                                                  %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  QuantizeImage() analyzes the colors within a reference image and chooses a
%  fixed number of colors to represent the image.  The goal of the algorithm
%  is to minimize the color difference between the input and output image while
%  minimizing the processing time.
%
%  The format of the QuantizeImage method is:
%
%      MagickBooleanType QuantizeImage(const QuantizeInfo *quantize_info,
%        Image *image)
%
%  A description of each parameter follows:
%
%    o quantize_info: Specifies a pointer to an QuantizeInfo structure.
%
%    o image: the image.
%
*/

static MagickBooleanType DirectToColormapImage(Image *image,
  ExceptionInfo *exception)
{
  CacheView
    *image_view;

  MagickBooleanType
    status;

  register ssize_t
    i;

  size_t
    number_colors;

  ssize_t
    y;

  status=MagickTrue;
  number_colors=(size_t) (image->columns*image->rows);
  if (AcquireImageColormap(image,number_colors) == MagickFalse)
    ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
      image->filename);
  if (image->colors != number_colors)
    return(MagickFalse);
  i=0;
  image_view=AcquireAuthenticCacheView(image,exception);
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    MagickBooleanType
      proceed;

    register IndexPacket
      *magick_restrict indexes;

    register PixelPacket
      *magick_restrict q;

    register ssize_t
      x;

    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
    if (q == (const PixelPacket *) NULL)
      break;
    indexes=GetCacheViewAuthenticIndexQueue(image_view);
    for (x=0; x < (ssize_t) image->columns; x++)
    {
      image->colormap[i].red=GetPixelRed(q);
      image->colormap[i].green=GetPixelGreen(q);
      image->colormap[i].blue=GetPixelBlue(q);
      image->colormap[i].opacity=GetPixelOpacity(q);
      SetPixelIndex(indexes+x,i);
      i++;
      q++;
    }
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      break;
    proceed=SetImageProgress(image,AssignImageTag,(MagickOffsetType) y,
      image->rows);
    if (proceed == MagickFalse)
      status=MagickFalse;
  }
  image_view=DestroyCacheView(image_view);
  return(status);
}

MagickExport MagickBooleanType QuantizeImage(const QuantizeInfo *quantize_info,
  Image *image)
{
  CubeInfo
    *cube_info;

  MagickBooleanType
    status;

  size_t
    depth,
    maximum_colors;

  assert(quantize_info != (const QuantizeInfo *) NULL);
  assert(quantize_info->signature == MagickSignature);
  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  maximum_colors=quantize_info->number_colors;
  if (maximum_colors == 0)
    maximum_colors=MaxColormapSize;
  if (maximum_colors > MaxColormapSize)
    maximum_colors=MaxColormapSize;
  if (image->matte == MagickFalse)
    {
      if ((image->columns*image->rows) <= maximum_colors)
        (void) DirectToColormapImage(image,&image->exception);
      if (SetImageGray(image,&image->exception) != MagickFalse)
        (void) SetGrayscaleImage(image);
    }
  if ((image->storage_class == PseudoClass) &&
      (image->colors <= maximum_colors))
    {
      if ((quantize_info->colorspace != UndefinedColorspace) &&
          (quantize_info->colorspace != CMYKColorspace))
        (void) TransformImageColorspace(image,quantize_info->colorspace);
      return(MagickTrue);
    }
  depth=quantize_info->tree_depth;
  if (depth == 0)
    {
      size_t
        colors;

      /*
        Depth of color tree is: Log4(colormap size)+2.
      */
      colors=maximum_colors;
      for (depth=1; colors != 0; depth++)
        colors>>=2;
      if ((quantize_info->dither != MagickFalse) && (depth > 2))
        depth--;
      if ((image->matte != MagickFalse) && (depth > 5))
        depth--;
      if (SetImageGray(image,&image->exception) != MagickFalse)
        depth=MaxTreeDepth;
    }
  /*
    Initialize color cube.
  */
  cube_info=GetCubeInfo(quantize_info,depth,maximum_colors);
  if (cube_info == (CubeInfo *) NULL)
    ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
      image->filename);
  status=ClassifyImageColors(cube_info,image,&image->exception);
  if (status != MagickFalse)
    {
      /*
        Reduce the number of colors in the image if it contains more than the
        maximum, otherwise we can disable dithering to improve the performance.
      */
      if (cube_info->colors > cube_info->maximum_colors)
        ReduceImageColors(image,cube_info);
      else
        cube_info->quantize_info->dither_method=NoDitherMethod;
      status=AssignImageColors(image,cube_info);
    }
  DestroyCubeInfo(cube_info);
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   Q u a n t i z e I m a g e s                                               %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  QuantizeImages() analyzes the colors within a set of reference images and
%  chooses a fixed number of colors to represent the set.  The goal of the
%  algorithm is to minimize the color difference between the input and output
%  images while minimizing the processing time.
%
%  The format of the QuantizeImages method is:
%
%      MagickBooleanType QuantizeImages(const QuantizeInfo *quantize_info,
%        Image *images)
%
%  A description of each parameter follows:
%
%    o quantize_info: Specifies a pointer to an QuantizeInfo structure.
%
%    o images: Specifies a pointer to a list of Image structures.
%
*/
MagickExport MagickBooleanType QuantizeImages(const QuantizeInfo *quantize_info,
  Image *images)
{
  CubeInfo
    *cube_info;

  Image
    *image;

  MagickBooleanType
    proceed,
    status;

  MagickProgressMonitor
    progress_monitor;

  register ssize_t
    i;

  size_t
    depth,
    maximum_colors,
    number_images;

  assert(quantize_info != (const QuantizeInfo *) NULL);
  assert(quantize_info->signature == MagickSignature);
  assert(images != (Image *) NULL);
  assert(images->signature == MagickSignature);
  if (images->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
  if (GetNextImageInList(images) == (Image *) NULL)
    {
      /*
        Handle a single image with QuantizeImage.
      */
      status=QuantizeImage(quantize_info,images);
      return(status);
    }
  status=MagickFalse;
  maximum_colors=quantize_info->number_colors;
  if (maximum_colors == 0)
    maximum_colors=MaxColormapSize;
  if (maximum_colors > MaxColormapSize)
    maximum_colors=MaxColormapSize;
  depth=quantize_info->tree_depth;
  if (depth == 0)
    {
      size_t
        colors;

      /*
        Depth of color tree is: Log4(colormap size)+2.
      */
      colors=maximum_colors;
      for (depth=1; colors != 0; depth++)
        colors>>=2;
      if (quantize_info->dither != MagickFalse)
        depth--;
    }
  /*
    Initialize color cube.
  */
  cube_info=GetCubeInfo(quantize_info,depth,maximum_colors);
  if (cube_info == (CubeInfo *) NULL)
    {
      (void) ThrowMagickException(&images->exception,GetMagickModule(),
        ResourceLimitError,"MemoryAllocationFailed","`%s'",images->filename);
      return(MagickFalse);
    }
  number_images=GetImageListLength(images);
  image=images;
  for (i=0; image != (Image *) NULL; i++)
  {
    progress_monitor=SetImageProgressMonitor(image,(MagickProgressMonitor) NULL,
      image->client_data);
    status=ClassifyImageColors(cube_info,image,&image->exception);
    if (status == MagickFalse)
      break;
    (void) SetImageProgressMonitor(image,progress_monitor,image->client_data);
    proceed=SetImageProgress(image,AssignImageTag,(MagickOffsetType) i,
      number_images);
    if (proceed == MagickFalse)
      break;
    image=GetNextImageInList(image);
  }
  if (status != MagickFalse)
    {
      /*
        Reduce the number of colors in an image sequence.
      */
      ReduceImageColors(images,cube_info);
      image=images;
      for (i=0; image != (Image *) NULL; i++)
      {
        progress_monitor=SetImageProgressMonitor(image,(MagickProgressMonitor)
          NULL,image->client_data);
        status=AssignImageColors(image,cube_info);
        if (status == MagickFalse)
          break;
        (void) SetImageProgressMonitor(image,progress_monitor,
          image->client_data);
        proceed=SetImageProgress(image,AssignImageTag,(MagickOffsetType) i,
          number_images);
        if (proceed == MagickFalse)
          break;
        image=GetNextImageInList(image);
      }
    }
  DestroyCubeInfo(cube_info);
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   Q u a n t i z e E r r o r F l a t t e n                                   %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  QuantizeErrorFlatten() traverses the color cube and flattens the quantization
%  error into a sorted 1D array.  This accelerates the color reduction process.
%
%  Contributed by Yoya.
%
%  The format of the QuantizeErrorFlatten method is:
%
%      size_t QuantizeErrorFlatten(const CubeInfo *cube_info,
%        const NodeInfo *node_info,const ssize_t offset,
%        MagickRealType *quantize_error)
%
%  A description of each parameter follows.
%
%    o cube_info: A pointer to the Cube structure.
%
%    o node_info: pointer to node in color cube tree that is current pointer.
%
%    o offset: quantize error offset.
%
%    o quantize_error: the quantization error vector.
%
*/
static size_t QuantizeErrorFlatten(const CubeInfo *cube_info,
  const NodeInfo *node_info,const ssize_t offset,
  MagickRealType *quantize_error)
{
  register ssize_t
    i;

  size_t
    n,
    number_children;

  if (offset >= (ssize_t) cube_info->nodes)
    return(0);
  quantize_error[offset]=node_info->quantize_error;
  n=1;
  number_children=cube_info->associate_alpha == MagickFalse ? 8UL : 16UL;
  for (i=0; i < (ssize_t) number_children ; i++)
    if (node_info->child[i] != (NodeInfo *) NULL)
      n+=QuantizeErrorFlatten(cube_info,node_info->child[i],offset+n,
        quantize_error);
  return(n);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   R e d u c e                                                               %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Reduce() traverses the color cube tree and prunes any node whose
%  quantization error falls below a particular threshold.
%
%  The format of the Reduce method is:
%
%      Reduce(CubeInfo *cube_info,const NodeInfo *node_info)
%
%  A description of each parameter follows.
%
%    o cube_info: A pointer to the Cube structure.
%
%    o node_info: pointer to node in color cube tree that is to be pruned.
%
*/
static void Reduce(CubeInfo *cube_info,const NodeInfo *node_info)
{
  register ssize_t
    i;

  size_t
    number_children;

  /*
    Traverse any children.
  */
  number_children=cube_info->associate_alpha == MagickFalse ? 8UL : 16UL;
  for (i=0; i < (ssize_t) number_children; i++)
    if (node_info->child[i] != (NodeInfo *) NULL)
      Reduce(cube_info,node_info->child[i]);
  if (node_info->quantize_error <= cube_info->pruning_threshold)
    PruneChild(cube_info,node_info);
  else
    {
      /*
        Find minimum pruning threshold.
      */
      if (node_info->number_unique > 0)
        cube_info->colors++;
      if (node_info->quantize_error < cube_info->next_threshold)
        cube_info->next_threshold=node_info->quantize_error;
    }
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
+   R e d u c e I m a g e C o l o r s                                         %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  ReduceImageColors() repeatedly prunes the tree until the number of nodes
%  with n2 > 0 is less than or equal to the maximum number of colors allowed
%  in the output image.  On any given iteration over the tree, it selects
%  those nodes whose E value is minimal for pruning and merges their
%  color statistics upward. It uses a pruning threshold, Ep, to govern
%  node selection as follows:
%
%    Ep = 0
%    while number of nodes with (n2 > 0) > required maximum number of colors
%      prune all nodes such that E <= Ep
%      Set Ep to minimum E in remaining nodes
%
%  This has the effect of minimizing any quantization error when merging
%  two nodes together.
%
%  When a node to be pruned has offspring, the pruning procedure invokes
%  itself recursively in order to prune the tree from the leaves upward.
%  n2,  Sr, Sg,  and  Sb in a node being pruned are always added to the
%  corresponding data in that node's parent.  This retains the pruned
%  node's color characteristics for later averaging.
%
%  For each node, n2 pixels exist for which that node represents the
%  smallest volume in RGB space containing those pixel's colors.  When n2
%  > 0 the node will uniquely define a color in the output image. At the
%  beginning of reduction,  n2 = 0  for all nodes except a the leaves of
%  the tree which represent colors present in the input image.
%
%  The other pixel count, n1, indicates the total number of colors
%  within the cubic volume which the node represents.  This includes n1 -
%  n2  pixels whose colors should be defined by nodes at a lower level in
%  the tree.
%
%  The format of the ReduceImageColors method is:
%
%      ReduceImageColors(const Image *image,CubeInfo *cube_info)
%
%  A description of each parameter follows.
%
%    o image: the image.
%
%    o cube_info: A pointer to the Cube structure.
%
*/

static int MagickRealTypeCompare(const void *error_p,const void *error_q)
{
  MagickRealType
    *p,
    *q;

  p=(MagickRealType *) error_p;
  q=(MagickRealType *) error_q;
  if (*p > *q)
    return(1);
  if (fabs((double) (*q-*p)) <= MagickEpsilon)
    return(0);
  return(-1);
}

static void ReduceImageColors(const Image *image,CubeInfo *cube_info)
{
#define ReduceImageTag  "Reduce/Image"

  MagickBooleanType
    proceed;

  MagickOffsetType
    offset;

  size_t
    span;

  cube_info->next_threshold=0.0;
  if (cube_info->colors > cube_info->maximum_colors)
    {
      MagickRealType
        *quantize_error;

      /*
        Enable rapid reduction of the number of unique colors.
      */
      quantize_error=(MagickRealType *) AcquireQuantumMemory(cube_info->nodes,
        sizeof(*quantize_error));
      if (quantize_error != (MagickRealType *) NULL)
        {
          (void) QuantizeErrorFlatten(cube_info,cube_info->root,0,
            quantize_error);
          qsort(quantize_error,cube_info->nodes,sizeof(MagickRealType),
            MagickRealTypeCompare);
          if (cube_info->nodes > (110*(cube_info->maximum_colors+1)/100))
            cube_info->next_threshold=quantize_error[cube_info->nodes-110*
              (cube_info->maximum_colors+1)/100];
          quantize_error=(MagickRealType *) RelinquishMagickMemory(
            quantize_error);
        }
  }
  for (span=cube_info->colors; cube_info->colors > cube_info->maximum_colors; )
  {
    cube_info->pruning_threshold=cube_info->next_threshold;
    cube_info->next_threshold=cube_info->root->quantize_error-1;
    cube_info->colors=0;
    Reduce(cube_info,cube_info->root);
    offset=(MagickOffsetType) span-cube_info->colors;
    proceed=SetImageProgress(image,ReduceImageTag,offset,span-
      cube_info->maximum_colors+1);
    if (proceed == MagickFalse)
      break;
  }
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   R e m a p I m a g e                                                       %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  RemapImage() replaces the colors of an image with the closest color from
%  a reference image.
%
%  The format of the RemapImage method is:
%
%      MagickBooleanType RemapImage(const QuantizeInfo *quantize_info,
%        Image *image,const Image *remap_image)
%
%  A description of each parameter follows:
%
%    o quantize_info: Specifies a pointer to an QuantizeInfo structure.
%
%    o image: the image.
%
%    o remap_image: the reference image.
%
*/
MagickExport MagickBooleanType RemapImage(const QuantizeInfo *quantize_info,
  Image *image,const Image *remap_image)
{
  CubeInfo
    *cube_info;

  MagickBooleanType
    status;

  /*
    Initialize color cube.
  */
  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
  assert(remap_image != (Image *) NULL);
  assert(remap_image->signature == MagickSignature);
  cube_info=GetCubeInfo(quantize_info,MaxTreeDepth,
    quantize_info->number_colors);
  if (cube_info == (CubeInfo *) NULL)
    ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
      image->filename);
  status=ClassifyImageColors(cube_info,remap_image,&image->exception);
  if (status != MagickFalse)
    {
      /*
        Classify image colors from the reference image.
      */
      cube_info->quantize_info->number_colors=cube_info->colors;
      status=AssignImageColors(image,cube_info);
    }
  DestroyCubeInfo(cube_info);
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   R e m a p I m a g e s                                                     %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  RemapImages() replaces the colors of a sequence of images with the
%  closest color from a reference image.
%
%  The format of the RemapImage method is:
%
%      MagickBooleanType RemapImages(const QuantizeInfo *quantize_info,
%        Image *images,Image *remap_image)
%
%  A description of each parameter follows:
%
%    o quantize_info: Specifies a pointer to an QuantizeInfo structure.
%
%    o images: the image sequence.
%
%    o remap_image: the reference image.
%
*/
MagickExport MagickBooleanType RemapImages(const QuantizeInfo *quantize_info,
  Image *images,const Image *remap_image)
{
  CubeInfo
    *cube_info;

  Image
    *image;

  MagickBooleanType
    status;

  assert(images != (Image *) NULL);
  assert(images->signature == MagickSignature);
  if (images->debug != MagickFalse)
    (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename);
  image=images;
  if (remap_image == (Image *) NULL)
    {
      /*
        Create a global colormap for an image sequence.
      */
      status=QuantizeImages(quantize_info,images);
      return(status);
    }
  /*
    Classify image colors from the reference image.
  */
  cube_info=GetCubeInfo(quantize_info,MaxTreeDepth,
    quantize_info->number_colors);
  if (cube_info == (CubeInfo *) NULL)
    ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
      image->filename);
  status=ClassifyImageColors(cube_info,remap_image,&image->exception);
  if (status != MagickFalse)
    {
      /*
        Classify image colors from the reference image.
      */
      cube_info->quantize_info->number_colors=cube_info->colors;
      image=images;
      for ( ; image != (Image *) NULL; image=GetNextImageInList(image))
      {
        status=AssignImageColors(image,cube_info);
        if (status == MagickFalse)
          break;
      }
    }
  DestroyCubeInfo(cube_info);
  return(status);
}

/*
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                             %
%                                                                             %
%                                                                             %
%   S e t G r a y s c a l e I m a g e                                         %
%                                                                             %
%                                                                             %
%                                                                             %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  SetGrayscaleImage() converts an image to a PseudoClass grayscale image.
%
%  The format of the SetGrayscaleImage method is:
%
%      MagickBooleanType SetGrayscaleImage(Image *image)
%
%  A description of each parameter follows:
%
%    o image: The image.
%
*/

#if defined(__cplusplus) || defined(c_plusplus)
extern "C" {
#endif

static int IntensityCompare(const void *x,const void *y)
{
  PixelPacket
    *color_1,
    *color_2;

  int
    intensity;

  color_1=(PixelPacket *) x;
  color_2=(PixelPacket *) y;
  intensity=PixelPacketIntensity(color_1)-(int) PixelPacketIntensity(color_2);
  return((int) intensity);
}

#if defined(__cplusplus) || defined(c_plusplus)
}
#endif

static MagickBooleanType SetGrayscaleImage(Image *image)
{
  CacheView
    *image_view;

  ExceptionInfo
    *exception;

  MagickBooleanType
    status;

  PixelPacket
    *colormap;

  register ssize_t
    i;

  ssize_t
    *colormap_index,
    j,
    y;

  assert(image != (Image *) NULL);
  assert(image->signature == MagickSignature);
  if (image->type != GrayscaleType)
    (void) TransformImageColorspace(image,GRAYColorspace);
  colormap_index=(ssize_t *) AcquireQuantumMemory(MaxColormapSize,
    sizeof(*colormap_index));
  if (colormap_index == (ssize_t *) NULL)
    ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
      image->filename);
  if (image->storage_class != PseudoClass)
    {
      ExceptionInfo
        *exception;

      (void) ResetMagickMemory(colormap_index,(-1),MaxColormapSize*
        sizeof(*colormap_index));
      if (AcquireImageColormap(image,MaxColormapSize) == MagickFalse)
        ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
          image->filename);
      image->colors=0;
      status=MagickTrue;
      exception=(&image->exception);
      image_view=AcquireAuthenticCacheView(image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
      #pragma omp parallel for schedule(static,4) shared(status) \
        magick_threads(image,image,image->rows,1)
#endif
      for (y=0; y < (ssize_t) image->rows; y++)
      {
        register IndexPacket
          *magick_restrict indexes;

        register PixelPacket
          *magick_restrict q;

        register ssize_t
          x;

        if (status == MagickFalse)
          continue;
        q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,
          exception);
        if (q == (PixelPacket *) NULL)
          {
            status=MagickFalse;
            continue;
          }
        indexes=GetCacheViewAuthenticIndexQueue(image_view);
        for (x=0; x < (ssize_t) image->columns; x++)
        {
          register size_t
            intensity;

          intensity=ScaleQuantumToMap(GetPixelRed(q));
          if (colormap_index[intensity] < 0)
            {
#if defined(MAGICKCORE_OPENMP_SUPPORT)
              #pragma omp critical (MagickCore_SetGrayscaleImage)
#endif
              if (colormap_index[intensity] < 0)
                {
                  colormap_index[intensity]=(ssize_t) image->colors;
                  image->colormap[image->colors].red=GetPixelRed(q);
                  image->colormap[image->colors].green=GetPixelGreen(q);
                  image->colormap[image->colors].blue=GetPixelBlue(q);
                  image->colors++;
               }
            }
          SetPixelIndex(indexes+x,colormap_index[intensity]);
          q++;
        }
        if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
          status=MagickFalse;
      }
      image_view=DestroyCacheView(image_view);
    }
  for (i=0; i < (ssize_t) image->colors; i++)
    image->colormap[i].opacity=(unsigned short) i;
  qsort((void *) image->colormap,image->colors,sizeof(PixelPacket),
    IntensityCompare);
  colormap=(PixelPacket *) AcquireQuantumMemory(image->colors,
    sizeof(*colormap));
  if (colormap == (PixelPacket *) NULL)
    ThrowBinaryException(ResourceLimitError,"MemoryAllocationFailed",
      image->filename);
  j=0;
  colormap[j]=image->colormap[0];
  for (i=0; i < (ssize_t) image->colors; i++)
  {
    if (IsSameColor(image,&colormap[j],&image->colormap[i]) == MagickFalse)
      {
        j++;
        colormap[j]=image->colormap[i];
      }
    colormap_index[(ssize_t) image->colormap[i].opacity]=j;
  }
  image->colors=(size_t) (j+1);
  image->colormap=(PixelPacket *) RelinquishMagickMemory(image->colormap);
  image->colormap=colormap;
  status=MagickTrue;
  exception=(&image->exception);
  image_view=AcquireAuthenticCacheView(image,exception);
#if defined(MAGICKCORE_OPENMP_SUPPORT)
  #pragma omp parallel for schedule(static,4) shared(status) \
    magick_threads(image,image,image->rows,1)
#endif
  for (y=0; y < (ssize_t) image->rows; y++)
  {
    register IndexPacket
      *magick_restrict indexes;

    register const PixelPacket
      *magick_restrict q;

    register ssize_t
      x;

    if (status == MagickFalse)
      continue;
    q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,exception);
    if (q == (PixelPacket *) NULL)
      {
        status=MagickFalse;
        continue;
      }
    indexes=GetCacheViewAuthenticIndexQueue(image_view);
    for (x=0; x < (ssize_t) image->columns; x++)
      SetPixelIndex(indexes+x,colormap_index[ScaleQuantumToMap(GetPixelIndex(
        indexes+x))]);
    if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)
      status=MagickFalse;
  }
  image_view=DestroyCacheView(image_view);
  colormap_index=(ssize_t *) RelinquishMagickMemory(colormap_index);
  image->type=GrayscaleType;
  if (SetImageMonochrome(image,&image->exception) != MagickFalse)
    image->type=BilevelType;
  return(status);
}