//C-  -*- C++ -*-

//C- -------------------------------------------------------------------

//C- DjVuLibre-3.5

//C- Copyright (c) 2002  Leon Bottou and Yann Le Cun.

//C- Copyright (c) 2001  AT&T

//C-

//C- This software is subject to, and may be distributed under, the

//C- GNU General Public License, either Version 2 of the license,

//C- or (at your option) any later version. The license should have

//C- accompanied the software or you may obtain a copy of the license

//C- from the Free Software Foundation at http://www.fsf.org .

//C-

//C- This program is distributed in the hope that it will be useful,

//C- but WITHOUT ANY WARRANTY; without even the implied warranty of

//C- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

//C- GNU General Public License for more details.

//C- 

//C- DjVuLibre-3.5 is derived from the DjVu(r) Reference Library from

//C- Lizardtech Software.  Lizardtech Software has authorized us to

//C- replace the original DjVu(r) Reference Library notice by the following

//C- text (see doc/lizard2002.djvu and doc/lizardtech2007.djvu):

//C-

//C-  ------------------------------------------------------------------

//C- | DjVu (r) Reference Library (v. 3.5)

//C- | Copyright (c) 1999-2001 LizardTech, Inc. All Rights Reserved.

//C- | The DjVu Reference Library is protected by U.S. Pat. No.

//C- | 6,058,214 and patents pending.

//C- |

//C- | This software is subject to, and may be distributed under, the

//C- | GNU General Public License, either Version 2 of the license,

//C- | or (at your option) any later version. The license should have

//C- | accompanied the software or you may obtain a copy of the license

//C- | from the Free Software Foundation at http://www.fsf.org .

//C- |

//C- | The computer code originally released by LizardTech under this

//C- | license and unmodified by other parties is deemed "the LIZARDTECH

//C- | ORIGINAL CODE."  Subject to any third party intellectual property

//C- | claims, LizardTech grants recipient a worldwide, royalty-free, 

//C- | non-exclusive license to make, use, sell, or otherwise dispose of 

//C- | the LIZARDTECH ORIGINAL CODE or of programs derived from the 

//C- | LIZARDTECH ORIGINAL CODE in compliance with the terms of the GNU 

//C- | General Public License.   This grant only confers the right to 

//C- | infringe patent claims underlying the LIZARDTECH ORIGINAL CODE to 

//C- | the extent such infringement is reasonably necessary to enable 

//C- | recipient to make, have made, practice, sell, or otherwise dispose 

//C- | of the LIZARDTECH ORIGINAL CODE (or portions thereof) and not to 

//C- | any greater extent that may be necessary to utilize further 

//C- | modifications or combinations.

//C- |

//C- | The LIZARDTECH ORIGINAL CODE is provided "AS IS" WITHOUT WARRANTY

//C- | OF ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED

//C- | TO ANY WARRANTY OF NON-INFRINGEMENT, OR ANY IMPLIED WARRANTY OF

//C- | MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.

//C- +------------------------------------------------------------------

#ifdef HAVE_CONFIG_H

# include "config.h"

#endif

#if NEED_GNUG_PRAGMAS

# pragma implementation

#endif

#include "GException.h"

#include "ByteStream.h"

#include "BSByteStream.h"

#include "DjVuPalette.h"

#include <stddef.h>

#include <stdlib.h>

#include <math.h>

#ifdef HAVE_NAMESPACES

namespace DJVU {

# ifdef NOT_DEFINED // Just to fool emacs c++ mode

}

#endif

#endif

#define CUBEBITS  4

#define CUBESIDE  (1<

#define CUBESIZE  (CUBESIDE*CUBESIDE*CUBESIDE)

#define RMUL 5

#define GMUL 9

#define BMUL 2

#define SMUL (RMUL+GMUL+BMUL)

#define MAXPALETTESIZE 65535 // Limit for a 16 bit unsigned read.

inline unsigned char 

umax(unsigned char a, unsigned char b) 

{ return (a>b) ? a : b; }

inline unsigned char 

umin(unsigned char a, unsigned char b) 

{ return (a>b) ? b : a; }

inline float 

fmin(float a, float b) 

{ return (a>b) ? b : a; }

// ------- DJVUPALETTE

DjVuPalette::DjVuPalette()

  : mask(0), hist(0), pmap(0)

{

}

DjVuPalette::~DjVuPalette()

{

  delete hist;

  delete pmap;

}

DjVuPalette& 

DjVuPalette::operator=(const DjVuPalette &ref)

{

  if (this != &ref)

    {

      delete hist;

      delete pmap;

      mask = 0;

      palette = ref.palette;

      colordata = ref.colordata;

    }

  return *this;

}

DjVuPalette::DjVuPalette(const DjVuPalette &ref)

  : mask(0), hist(0), pmap(0)

{

  this->operator=(ref);

}

// -------- HISTOGRAM ALLOCATION

void

DjVuPalette::allocate_hist()

{

  if (! hist)

    {

      hist = new GMap<int,int>;

      mask = 0;

    }

  else

    {

      GMap<int,int> *old = hist;

      hist = new GMap<int,int>;

      mask = (mask<<1)|(0x010101);

      for (GPosition p = *old; p; ++p)

        {

          int k = old->key(p);

          int w = (*old)[p];

          (*hist)[k | mask] += w;

        }

      delete old;

    }

}

// -------- PALETTE COMPUTATION

#ifndef NEED_DECODER_ONLY

struct PData

{

  unsigned char p[3];

  int w;

};

struct PBox 

{

  PData *data;

  int colors;

  int boxsize;

  int sum;

};

int

DjVuPalette::bcomp (const void *a, const void *b)

{

  return ((PData*)a)->p[0] - ((PData*)b)->p[0];

}

int

DjVuPalette::gcomp (const void *a, const void *b)

{

  return ((PData*)a)->p[1] - ((PData*)b)->p[1];

}

int

DjVuPalette::rcomp (const void *a, const void *b)

{

  return ((PData*)a)->p[2] - ((PData*)b)->p[2];

}

int

DjVuPalette::lcomp (const void *a, const void *b)

{

  unsigned char *aa = ((PColor*)a)->p;

  unsigned char *bb = ((PColor*)b)->p;

  if (aa[3] != bb[3])

    return aa[3]-bb[3];

  else if (aa[2] != bb[2])

    return aa[2]-bb[2];

  else if (aa[1] != bb[1])

    return aa[1]=bb[1];

  else

    return aa[0]-bb[0];

}

int

DjVuPalette::compute_palette(int maxcolors, int minboxsize)

{

  if (!hist)

    G_THROW( ERR_MSG("DjVuPalette.no_color") );

  if (maxcolors<1 || maxcolors>MAXPALETTESIZE)

    G_THROW( ERR_MSG("DjVuPalette.many_colors") );

  // Paul Heckbert: "Color Image Quantization for Frame Buffer Display", 

  // SIGGRAPH '82 Proceedings, page 297.  (also in ppmquant)

  // Collect histogram colors

  int sum = 0;

  int ncolors = 0;

  GTArray<PData> pdata;

  { // extra nesting for windows

    for (GPosition p = *hist; p; ++p)

    {

      pdata.touch(ncolors);

      PData &data = pdata[ncolors++];

      int k = hist->key(p);

      data.p[0] = (k>>16) & 0xff;

      data.p[1] = (k>>8) & 0xff;

      data.p[2] = (k) & 0xff;

      data.w = (*hist)[p];

      sum += data.w;

    }

  }

  // Create first box

  GList<PBox> boxes;

  PBox newbox;

  newbox.data = pdata;

  newbox.colors = ncolors;

  newbox.boxsize = 256;

  newbox.sum = sum;

  boxes.append(newbox);

  // Repeat spliting boxes

  while (boxes.size() < maxcolors)

    {

      // Find suitable box

      GPosition p;

      for (p=boxes; p; ++p)

        if (boxes[p].colors>=2 && boxes[p].boxsize>minboxsize) 

          break;

      if (! p)

        break;

      // Find box boundaries

      PBox &splitbox = boxes[p];

      unsigned char pmax[3];

      unsigned char pmin[3];

      pmax[0] = pmin[0] = splitbox.data->p[0];

      pmax[1] = pmin[1] = splitbox.data->p[1];

      pmax[2] = pmin[2] = splitbox.data->p[2];

      { // extra nesting for windows

        for (int j=1; j

        {

          pmax[0] = umax(pmax[0], splitbox.data[j].p[0]);

          pmax[1] = umax(pmax[1], splitbox.data[j].p[1]);

          pmax[2] = umax(pmax[2], splitbox.data[j].p[2]);

          pmin[0] = umin(pmin[0], splitbox.data[j].p[0]);

          pmin[1] = umin(pmin[1], splitbox.data[j].p[1]);

          pmin[2] = umin(pmin[2], splitbox.data[j].p[2]);

        }

      }

      // Determine split direction and sort

      int bl = pmax[0]-pmin[0]; 

      int gl = pmax[1]-pmin[1];

      int rl = pmax[2]-pmin[2];

      splitbox.boxsize = (bl>gl ? (rl>bl ? rl : bl) : (rl>gl ? rl : gl));

      if (splitbox.boxsize <= minboxsize)

        continue;

      if (gl == splitbox.boxsize)

        qsort(splitbox.data, splitbox.colors, sizeof(PData), gcomp);

      else if (rl == splitbox.boxsize)

        qsort(splitbox.data, splitbox.colors, sizeof(PData), rcomp);

      else

        qsort(splitbox.data, splitbox.colors, sizeof(PData), bcomp);

      // Find median

      int lowercolors = 0;

      int lowersum = 0;

      while (lowercolors

        lowersum += splitbox.data[lowercolors++].w;

      // Compute new boxes

      newbox.data = splitbox.data + lowercolors;

      newbox.colors = splitbox.colors - lowercolors;

      newbox.sum = splitbox.sum - lowersum;

      splitbox.colors = lowercolors;

      splitbox.sum = lowersum;

      // Insert boxes at proper location

      GPosition q;

      for (q=p; q; ++q)

        if (boxes[q].sum < newbox.sum)

          break;

      boxes.insert_before(q, newbox);

      for (q=p; q; ++q)

        if (boxes[q].sum < splitbox.sum)

          break;

      boxes.insert_before(q, boxes, p);

    }

  // Fill palette array

  ncolors = 0;

  palette.empty();

  palette.resize(0,boxes.size()-1);

  { // extra nesting for windows

    for (GPosition p=boxes; p; ++p)

    {

      PBox &box = boxes[p];

      // Compute box representative color

      float bsum = 0;

      float gsum = 0;

      float rsum = 0;

      for (int j=0; j

        {

          float w = (float)box.data[j].w;

          bsum += box.data[j].p[0] * w;

          gsum += box.data[j].p[1] * w;

          rsum += box.data[j].p[2] * w;

        }

      PColor &color = palette[ncolors++];

      color.p[0] = (unsigned char) fmin(255, bsum/box.sum);

      color.p[1] = (unsigned char) fmin(255, gsum/box.sum);

      color.p[2] = (unsigned char) fmin(255, rsum/box.sum);

      color.p[3] = ( color.p[0]*BMUL + color.p[1]*GMUL + color.p[2]*RMUL) / SMUL;

    }

  }

  // Save dominant color

  PColor dcolor = palette[0];

  // Sort palette colors in luminance order

  qsort((PColor*)palette, ncolors, sizeof(PColor), lcomp);

  // Clear invalid data

  colordata.empty();

  delete pmap;

  pmap = 0;

  // Return dominant color

  return color_to_index_slow(dcolor.p);

}

int 

DjVuPalette::compute_pixmap_palette(const GPixmap &pm, int ncolors, int minboxsize)

{

  // Prepare histogram

  histogram_clear();

  { // extra nesting for windows

    for (int j=0; j<(int)pm.rows(); j++)

    {

      const GPixel *p = pm[j];

      for (int i=0; i<(int)pm.columns(); i++)

        histogram_add(p[i], 1);

    }

  }

  // Compute palette

  return compute_palette(ncolors, minboxsize);

}

#endif

// -------- QUANTIZATION

void

DjVuPalette::allocate_pmap()

{

  if (! pmap)

    pmap = new GMap<int,int>;

}

int 

DjVuPalette::color_to_index_slow(const unsigned char *bgr)

{

  PColor *pal = palette;

  const int ncolors = palette.size();

  if (! ncolors)

    G_THROW( ERR_MSG("DjVuPalette.not_init") );

  // Should be able to do better

  int found = 0;

  int founddist = 3*256*256;

  { // extra nesting for windows

    for (int i=0; i

    {

      int bd = bgr[0] - pal[i].p[0];

      int gd = bgr[1] - pal[i].p[1];

      int rd = bgr[2] - pal[i].p[2];

      int dist = (bd*bd)+(gd*gd)+(rd*rd);

      if (dist < founddist)

        {

          found = i;

          founddist = dist;

        }

    }

  }

  // Store in pmap

  if (pmap && pmap->size()<0x8000)

    {

      int key = (bgr[0]<<16)|(bgr[1]<<8)|(bgr[2]);

      (*pmap)[key] = found;

    }

  // Return

  return found;

}

#ifndef NEED_DECODER_ONLY

void 

DjVuPalette::quantize(GPixmap &pm)

{

  { // extra nesting for windows

    for (int j=0; j<(int)pm.rows(); j++)

    {

      GPixel *p = pm[j];

      for (int i=0; i<(int)pm.columns(); i++)

        index_to_color(color_to_index(p[i]), p[i]);

    }

  }

}

int 

DjVuPalette::compute_palette_and_quantize(GPixmap &pm, int maxcolors, int minboxsize)

{

  int result = compute_pixmap_palette(pm, maxcolors, minboxsize);

  quantize(pm);

  return result;

}

void 

DjVuPalette::color_correct(double corr)

{

  const int palettesize = palette.size();

  if (palettesize > 0)

    {

      // Copy colors

      int i;

      GTArray<GPixel> pix(0,palettesize-1);

      GPixel *r = pix;

      PColor *q = palette;

      for (i=0; i

        {

          r[i].b = q[i].p[0];

          r[i].g = q[i].p[1];

          r[i].r = q[i].p[2];

        }

      // Apply color correction

      GPixmap::color_correct(corr, r, palettesize);

      // Restore colors

      for (i=0; i

        {

          q[i].p[0] = r[i].b;

          q[i].p[1] = r[i].g;

          q[i].p[2] = r[i].r;

        }

    }

}

#endif

// -------- ENCODE AND DECODE

#define DJVUPALETTEVERSION 0

void

DjVuPalette::encode_rgb_entries(ByteStream &bs) const

{

  const int palettesize = palette.size();

  { // extra nesting for windows

    for (int c=0; c

    {

      unsigned char p[3];

      p[2] = palette[c].p[0];

      p[1] = palette[c].p[1];

      p[0] = palette[c].p[2];

      bs.writall((const void*)p, 3);

    }

  }

}

void 

DjVuPalette::encode(GP<ByteStream> gbs) const

{

  ByteStream &bs=*gbs;

  const int palettesize = palette.size();

  const int datasize = colordata.size();

  // Code version number

  int version = DJVUPALETTEVERSION;

  if (datasize>0) version |= 0x80;

  bs.write8(version);

  // Code palette

  bs.write16(palettesize);

  { // extra nesting for windows

    for (int c=0; c

    {

      unsigned char p[3];

      p[0] = palette[c].p[0];

      p[1] = palette[c].p[1];

      p[2] = palette[c].p[2];

      bs.writall((const void*)p, 3);

    }

  }

  // Code colordata

  if (datasize > 0)

    {

      bs.write24(datasize);

      GP<ByteStream> gbsb=BSByteStream::create(gbs, 50);

      ByteStream &bsb=*gbsb;

      for (int d=0; d

        bsb.write16(colordata[d]);

    }

}

void 

DjVuPalette::decode_rgb_entries(ByteStream &bs, const int palettesize)

{

  palette.resize(0,palettesize-1);

  { // extra nesting for windows

    for (int c=0; c

    {

      unsigned char p[3];

      bs.readall((void*)p, 3);

      palette[c].p[0] = p[2];

      palette[c].p[1] = p[1];

      palette[c].p[2] = p[0];

      palette[c].p[3] = (p[0]*BMUL+p[1]*GMUL+p[2]*RMUL)/SMUL;

    }

  }

}

void 

DjVuPalette::decode(GP<ByteStream> gbs)

{

  ByteStream &bs=*gbs;

  // Make sure that everything is clear

  delete hist;

  delete pmap;

  hist = 0;

  pmap = 0;

  mask = 0;

  // Code version

  int version = bs.read8();

  if ( (version & 0x7f) != DJVUPALETTEVERSION)

    G_THROW( ERR_MSG("DjVuPalette.bad_version") );

  // Code palette

  const int palettesize = bs.read16();

  if (palettesize<0 || palettesize>MAXPALETTESIZE)

    G_THROW( ERR_MSG("DjVuPalette.bad_palette") );

  palette.resize(0,palettesize-1);

  { // extra nesting for windows

    for (int c=0; c

    {

      unsigned char p[3];

      bs.readall((void*)p, 3);

      palette[c].p[0] = p[0];

      palette[c].p[1] = p[1];

      palette[c].p[2] = p[2];

      palette[c].p[3] = (p[0]*BMUL+p[1]*GMUL+p[2]*RMUL)/SMUL;

    }

  }

  // Code data

  if (version & 0x80)

    {

      int datasize = bs.read24();

      if (datasize<0)

        G_THROW( ERR_MSG("DjVuPalette.bad_palette") );

      colordata.resize(0,datasize-1);

      GP<ByteStream> gbsb=BSByteStream::create(gbs);

      ByteStream &bsb=*gbsb;

      { // extra nesting for windows

        for (int d=0; d

        {

          short s = bsb.read16();

          if (s<0 || s>=palettesize)

            G_THROW( ERR_MSG("DjVuPalette.bad_palette") );        

          colordata[d] = s;

        }

      }

    }

}

#ifdef HAVE_NAMESPACES

}

# ifndef NOT_USING_DJVU_NAMESPACE

using namespace DJVU;

# endif

#endif