/*

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                                                                             %

%                                                                             %

%                                                                             %

%           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 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                                            %

%                                                                             %

%    https://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 "MagickCore/studio.h"

#include "MagickCore/artifact.h"

#include "MagickCore/attribute.h"

#include "MagickCore/cache-view.h"

#include "MagickCore/color.h"

#include "MagickCore/color-private.h"

#include "MagickCore/colormap.h"

#include "MagickCore/colorspace.h"

#include "MagickCore/colorspace-private.h"

#include "MagickCore/compare.h"

#include "MagickCore/enhance.h"

#include "MagickCore/exception.h"

#include "MagickCore/exception-private.h"

#include "MagickCore/histogram.h"

#include "MagickCore/image.h"

#include "MagickCore/image-private.h"

#include "MagickCore/list.h"

#include "MagickCore/memory_.h"

#include "MagickCore/memory-private.h"

#include "MagickCore/monitor.h"

#include "MagickCore/monitor-private.h"

#include "MagickCore/option.h"

#include "MagickCore/pixel-accessor.h"

#include "MagickCore/property.h"

#include "MagickCore/quantize.h"

#include "MagickCore/quantum.h"

#include "MagickCore/quantum-private.h"

#include "MagickCore/random_.h"

#include "MagickCore/resource_.h"

#include "MagickCore/string_.h"

#include "MagickCore/string-private.h"

#include "MagickCore/thread-private.h"

/*

  Define declarations.

*/

#if !defined(__APPLE__) && !defined(TARGET_OS_IPHONE)

#define CacheShift  2

#else

#define CacheShift  3

#endif

#define ErrorQueueLength  16

#define ErrorRelativeWeight  MagickSafeReciprocal(16)

#define MaxQNodes  266817

#define MaxTreeDepth  8

#define QNodesInAList  1920

/*

  Typedef declarations.

*/

typedef struct _DoublePixelPacket

{

  double

    red,

    green,

    blue,

    alpha;

} DoublePixelPacket;

typedef struct _QNodeInfo

{

  struct _QNodeInfo

    *parent,

    *child[16];

  MagickSizeType

    number_unique;

  DoublePixelPacket

    total_color;

  double

    quantize_error;

  size_t

    color_number,

    id,

    level;

} QNodeInfo;

typedef struct _QNodes

{

  QNodeInfo

    *nodes;

  struct _QNodes

    *next;

} QNodes;

typedef struct _QCubeInfo

{

  QNodeInfo

    *root;

  size_t

    colors,

    maximum_colors;

  ssize_t

    transparent_index;

  MagickSizeType

    transparent_pixels;

  DoublePixelPacket

    target;

  double

    distance,

    pruning_threshold,

    next_threshold;

  size_t

    nodes,

    free_nodes,

    color_number;

  QNodeInfo

    *next_node;

  QNodes

    *node_queue;

  MemoryInfo

    *memory_info;

  ssize_t

    *cache;

  DoublePixelPacket

    error[ErrorQueueLength];

  double

    diffusion,

    weights[ErrorQueueLength];

  QuantizeInfo

    *quantize_info;

  MagickBooleanType

    associate_alpha;

  ssize_t

    x,

    y;

  size_t

    depth;

  MagickOffsetType

    offset;

  MagickSizeType

    span;

} QCubeInfo;

/*

  Method prototypes.

*/

static QCubeInfo

  *GetQCubeInfo(const QuantizeInfo *,const size_t,const size_t);

static QNodeInfo

  *GetQNodeInfo(QCubeInfo *,const size_t,const size_t,QNodeInfo *);

static MagickBooleanType

  AssignImageColors(Image *,QCubeInfo *,ExceptionInfo *),

  ClassifyImageColors(QCubeInfo *,const Image *,ExceptionInfo *),

  DitherImage(Image *,QCubeInfo *,ExceptionInfo *),

  SetGrayscaleImage(Image *,ExceptionInfo *),

  SetImageColormap(Image *,QCubeInfo *,ExceptionInfo *);

static void

  ClosestColor(const Image *,QCubeInfo *,const QNodeInfo *),

  DefineImageColormap(Image *,QCubeInfo *,QNodeInfo *),

  DestroyQCubeInfo(QCubeInfo *),

  PruneLevel(QCubeInfo *,const QNodeInfo *),

  PruneToCubeDepth(QCubeInfo *,const QNodeInfo *),

  ReduceImageColors(const Image *,QCubeInfo *);

/*

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                                                                             %

%                                                                             %

%                                                                             %

%   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 *) AcquireCriticalMemory(sizeof(*quantize_info));

  GetQuantizeInfo(quantize_info);

  if (image_info != (ImageInfo *) NULL)

    {

      const char

        *option;

      quantize_info->dither_method=image_info->dither == MagickFalse ?

        NoDitherMethod : RiemersmaDitherMethod;

      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,QCubeInfo *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 Image *image,

  const QCubeInfo *cube_info,const Quantum *pixel,

  DoublePixelPacket *alpha_pixel)

{

  double

    alpha;

  if ((cube_info->associate_alpha == MagickFalse) ||

      (GetPixelAlpha(image,pixel) == OpaqueAlpha))

    {

      alpha_pixel->red=(double) GetPixelRed(image,pixel);

      alpha_pixel->green=(double) GetPixelGreen(image,pixel);

      alpha_pixel->blue=(double) GetPixelBlue(image,pixel);

      alpha_pixel->alpha=(double) GetPixelAlpha(image,pixel);

      return;

    }

  alpha=QuantumScale*(double) GetPixelAlpha(image,pixel);

  alpha_pixel->red=alpha*(double) GetPixelRed(image,pixel);

  alpha_pixel->green=alpha*(double) GetPixelGreen(image,pixel);

  alpha_pixel->blue=alpha*(double) GetPixelBlue(image,pixel);

  alpha_pixel->alpha=(double) GetPixelAlpha(image,pixel);

}

static inline void AssociateAlphaPixelInfo(const QCubeInfo *cube_info,

  const PixelInfo *pixel,DoublePixelPacket *alpha_pixel)

{

  double

    alpha;

  if ((cube_info->associate_alpha == MagickFalse) ||

      (pixel->alpha == (double) OpaqueAlpha))

    {

      alpha_pixel->red=(double) pixel->red;

      alpha_pixel->green=(double) pixel->green;

      alpha_pixel->blue=(double) pixel->blue;

      alpha_pixel->alpha=(double) pixel->alpha;

      return;

    }

  alpha=(double) (QuantumScale*pixel->alpha);

  alpha_pixel->red=alpha*pixel->red;

  alpha_pixel->green=alpha*pixel->green;

  alpha_pixel->blue=alpha*pixel->blue;

  alpha_pixel->alpha=(double) pixel->alpha;

}

static inline size_t ColorToQNodeId(const QCubeInfo *cube_info,

  const DoublePixelPacket *pixel,size_t index)

{

  size_t

    id;

  id=(size_t) (((ScaleQuantumToChar(ClampPixel(pixel->red)) >> index) & 0x01) |

    ((ScaleQuantumToChar(ClampPixel(pixel->green)) >> index) & 0x01) << 1 |

    ((ScaleQuantumToChar(ClampPixel(pixel->blue)) >> index) & 0x01) << 2);

  if (cube_info->associate_alpha != MagickFalse)

    id|=((((size_t) ScaleQuantumToChar(ClampPixel(pixel->alpha)) >> index) &

      0x1) << 3);

  return(id);

}

static MagickBooleanType AssignImageColors(Image *image,QCubeInfo *cube_info,

  ExceptionInfo *exception)

{

#define AssignImageTag  "Assign/Image"

  ColorspaceType

    colorspace;

  ssize_t

    y;

  /*

    Allocate image colormap.

  */

  colorspace=image->colorspace;

  if (cube_info->quantize_info->colorspace != UndefinedColorspace)

    (void) TransformImageColorspace(image,cube_info->quantize_info->colorspace,

      exception);

  cube_info->transparent_pixels=0;

  cube_info->transparent_index=(-1);

  if (SetImageColormap(image,cube_info,exception) == MagickFalse)

    return(MagickFalse);

  /*

    Create a reduced color image.

  */

  if (cube_info->quantize_info->dither_method != NoDitherMethod)

    (void) DitherImage(image,cube_info,exception);

  else

    {

      CacheView

        *image_view;

      MagickBooleanType

        status;

      status=MagickTrue;

      image_view=AcquireAuthenticCacheView(image,exception);

#if defined(MAGICKCORE_OPENMP_SUPPORT)

      #pragma omp parallel for schedule(static) shared(status) \

        magick_number_threads(image,image,image->rows,1)

#endif

      for (y=0; y < (ssize_t) image->rows; y++)

      {

        QCubeInfo

          cube;

        Quantum

          *magick_restrict q;

        ssize_t

          count,

          x;

        if (status == MagickFalse)

          continue;

        q=GetCacheViewAuthenticPixels(image_view,0,y,image->columns,1,

          exception);

        if (q == (Quantum *) NULL)

          {

            status=MagickFalse;

            continue;

          }

        cube=(*cube_info);

        for (x=0; x < (ssize_t) image->columns; x+=count)

        {

          DoublePixelPacket

            pixel;

          const QNodeInfo

            *node_info;

          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++)

          {

            PixelInfo

              packet;

            GetPixelInfoPixel(image,q+count*(ssize_t) GetPixelChannels(image),

              &packet);

            if (IsPixelEquivalent(image,q,&packet) == MagickFalse)

              break;

          }

          AssociateAlphaPixel(image,&cube,q,&pixel);

          node_info=cube.root;

          for (index=MaxTreeDepth-1; (ssize_t) index > 0; index--)

          {

            id=ColorToQNodeId(&cube,&pixel,index);

            if (node_info->child[id] == (QNodeInfo *) NULL)

              break;

            node_info=node_info->child[id];

          }

          /*

            Find closest color among siblings and their children.

          */

          cube.target=pixel;

          cube.distance=(double) (4.0*((double) QuantumRange+1.0)*

            ((double) 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(image,(Quantum) index,q);

            if (cube.quantize_info->measure_error == MagickFalse)

              {

                SetPixelRed(image,ClampToQuantum(

                  image->colormap[index].red),q);

                SetPixelGreen(image,ClampToQuantum(

                  image->colormap[index].green),q);

                SetPixelBlue(image,ClampToQuantum(

                  image->colormap[index].blue),q);

                if (cube.associate_alpha != MagickFalse)

                  SetPixelAlpha(image,ClampToQuantum(

                    image->colormap[index].alpha),q);

              }

            q+=(ptrdiff_t) GetPixelChannels(image);

          }

        }

        if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse)

          status=MagickFalse;

        if (image->progress_monitor != (MagickProgressMonitor) NULL)

          {

            MagickBooleanType

              proceed;

            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,exception);

  if ((cube_info->quantize_info->number_colors == 2) &&

      (IsGrayColorspace(cube_info->quantize_info->colorspace)))

    {

      double

        intensity;

      /*

        Monochrome image.

      */

      intensity=GetPixelInfoLuma(image->colormap+0) < (double)

        QuantumRange/2.0 ? 0.0 : (double) QuantumRange;

      if (image->colors > 1)

        {

          intensity=0.0;

          if (GetPixelInfoLuma(image->colormap+0) >

              GetPixelInfoLuma(image->colormap+1))

            intensity=(double) QuantumRange;

        }

      image->colormap[0].red=intensity;

      image->colormap[0].green=intensity;

      image->colormap[0].blue=intensity;

      if (image->colors > 1)

        {

          image->colormap[1].red=(double) QuantumRange-intensity;

          image->colormap[1].green=(double) QuantumRange-intensity;

          image->colormap[1].blue=(double) QuantumRange-intensity;

        }

    }

  (void) SyncImage(image,exception);

  if ((cube_info->quantize_info->colorspace != UndefinedColorspace) &&

      (IssRGBCompatibleColorspace(colorspace) == MagickFalse))

    (void) TransformImageColorspace(image,colorspace,exception);

  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(QCubeInfo *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,QCubeInfo *cube_info)

{

  MagickBooleanType

    associate_alpha;

  associate_alpha=image->alpha_trait != UndefinedPixelTrait ? MagickTrue :

    MagickFalse;

  if ((cube_info->quantize_info->number_colors == 2) &&

      ((cube_info->quantize_info->colorspace == LinearGRAYColorspace) ||

       (cube_info->quantize_info->colorspace == GRAYColorspace)))

    associate_alpha=MagickFalse;

  cube_info->associate_alpha=associate_alpha;

}

static MagickBooleanType ClassifyImageColors(QCubeInfo *cube_info,

  const Image *image,ExceptionInfo *exception)

{

#define ClassifyImageTag  "Classify/Image"

  CacheView

    *image_view;

  double

    bisect;

  DoublePixelPacket

    error,

    mid,

    midpoint,

    pixel;

  MagickBooleanType

    proceed;

  QNodeInfo

    *node_info;

  size_t

    id,

    index,

    level;

  ssize_t

    count,

    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 != image->colorspace)

    {

      if ((cube_info->quantize_info->colorspace != UndefinedColorspace) &&

          (cube_info->quantize_info->colorspace != CMYKColorspace))

        (void) TransformImageColorspace((Image *) image,

          cube_info->quantize_info->colorspace,exception);

      else

        if (IssRGBCompatibleColorspace(image->colorspace) == MagickFalse)

          (void) TransformImageColorspace((Image *) image,sRGBColorspace,

            exception);

    }

  midpoint.red=(double) QuantumRange/2.0;

  midpoint.green=(double) QuantumRange/2.0;

  midpoint.blue=(double) QuantumRange/2.0;

  midpoint.alpha=(double) QuantumRange/2.0;

  error.alpha=0.0;

  image_view=AcquireVirtualCacheView(image,exception);

  for (y=0; y < (ssize_t) image->rows; y++)

  {

    const Quantum

      *magick_restrict p;

    ssize_t

      x;

    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);

    if (p == (const Quantum *) NULL)

      break;

    if (cube_info->nodes > MaxQNodes)

      {

        /*

          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++)

      {

        PixelInfo

          packet;

        GetPixelInfoPixel(image,p+count*(ssize_t) GetPixelChannels(image),

          &packet);

        if (IsPixelEquivalent(image,p,&packet) == MagickFalse)

          break;

      }

      AssociateAlphaPixel(image,cube_info,p,&pixel);

      index=MaxTreeDepth-1;

      bisect=((double) QuantumRange+1.0)/2.0;

      mid=midpoint;

      node_info=cube_info->root;

      for (level=1; level <= MaxTreeDepth; level++)

      {

        double

          distance;

        bisect*=0.5;

        id=ColorToQNodeId(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.alpha+=(id & 8) != 0 ? bisect : -bisect;

        if (node_info->child[id] == (QNodeInfo *) NULL)

          {

            /*

              Set colors of new node to contain pixel.

            */

            node_info->child[id]=GetQNodeInfo(cube_info,id,level,node_info);

            if (node_info->child[id] == (QNodeInfo *) 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.alpha=QuantumScale*(pixel.alpha-mid.alpha);

        distance=(double) (error.red*error.red+error.green*error.green+

          error.blue*error.blue+error.alpha*error.alpha);

        if (IsNaN(distance) != 0)

          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=(size_t) ((ssize_t) node_info->number_unique+

        count);

      node_info->total_color.red+=count*QuantumScale*(double)

        ClampPixel(pixel.red);

      node_info->total_color.green+=count*QuantumScale*(double)

        ClampPixel(pixel.green);

      node_info->total_color.blue+=count*QuantumScale*(double)

        ClampPixel(pixel.blue);

      if (cube_info->associate_alpha != MagickFalse)

        node_info->total_color.alpha+=count*QuantumScale*(double)

          ClampPixel(pixel.alpha);

      else

        node_info->total_color.alpha+=count*QuantumScale*(double)

          ClampPixel((double) OpaqueAlpha);

      p+=(ptrdiff_t) count*(ssize_t) GetPixelChannels(image);

    }

    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++)

  {

    const Quantum

      *magick_restrict p;

    ssize_t

      x;

    p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);

    if (p == (const Quantum *) NULL)

      break;

    if (cube_info->nodes > MaxQNodes)

      {

        /*

          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++)

      {

        PixelInfo

          packet;

        GetPixelInfoPixel(image,p+count*(ssize_t) GetPixelChannels(image),

          &packet);

        if (IsPixelEquivalent(image,p,&packet) == MagickFalse)

          break;

      }

      AssociateAlphaPixel(image,cube_info,p,&pixel);

      index=MaxTreeDepth-1;

      bisect=((double) 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=ColorToQNodeId(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.alpha+=(id & 8) != 0 ? bisect : -bisect;

        if (node_info->child[id] == (QNodeInfo *) NULL)

          {

            /*

              Set colors of new node to contain pixel.

            */

            node_info->child[id]=GetQNodeInfo(cube_info,id,level,node_info);

            if (node_info->child[id] == (QNodeInfo *) 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.alpha=QuantumScale*(pixel.alpha-mid.alpha);

        distance=(double) (error.red*error.red+error.green*error.green+

          error.blue*error.blue+error.alpha*error.alpha);

        if (IsNaN(distance) != 0)

          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=(size_t) ((ssize_t) node_info->number_unique+

        count);

      node_info->total_color.red+=count*QuantumScale*(double)

        ClampPixel(pixel.red);

      node_info->total_color.green+=count*QuantumScale*(double)

        ClampPixel(pixel.green);

      node_info->total_color.blue+=count*QuantumScale*(double)

        ClampPixel(pixel.blue);

      if (cube_info->associate_alpha != MagickFalse)

        node_info->total_color.alpha+=count*QuantumScale*(double)

          ClampPixel(pixel.alpha);

      else

        node_info->total_color.alpha+=count*QuantumScale*(double)

          ClampPixel((MagickRealType) OpaqueAlpha);

      p+=(ptrdiff_t) count*(ssize_t) GetPixelChannels(image);

    }

    proceed=SetImageProgress(image,ClassifyImageTag,(MagickOffsetType) y,

      image->rows);

    if (proceed == MagickFalse)

      break;

  }

  image_view=DestroyCacheView(image_view);

  if (cube_info->quantize_info->colorspace != image->colorspace)

    if ((cube_info->quantize_info->colorspace != UndefinedColorspace) &&

        (cube_info->quantize_info->colorspace != CMYKColorspace))

      (void) TransformImageColorspace((Image *) image,sRGBColorspace,exception);

  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;