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

// - Author: Leon Bottou, 08/1998

// From: Leon Bottou, 1/31/2002

// Lizardtech has split this file into a decoder and an encoder.

// Only superficial changes.  The meat is mine.

#define IW44IMAGE_IMPLIMENTATION /* */

#include "IW44Image.h"

#include "ZPCodec.h"

#include "GBitmap.h"

#include "GPixmap.h"

#include "IFFByteStream.h"

#include "GRect.h"

#include <stddef.h>

#include <stdlib.h>

#include <string.h>

#include <math.h>

#include "MMX.h"

#undef IWTRANSFORM_TIMER

#ifdef IWTRANSFORM_TIMER

#include "GOS.h"

#endif

#include <assert.h>

#include <string.h>

#include <math.h>

#ifndef NEED_DECODER_ONLY

#ifdef HAVE_NAMESPACES

namespace DJVU {

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

}

#endif

#endif

#define IWALLOCSIZE    4080

#define IWCODEC_MAJOR     1

#define IWCODEC_MINOR     2

#define DECIBEL_PRUNE   5.0

//////////////////////////////////////////////////////

// WAVELET DECOMPOSITION CONSTANTS

//////////////////////////////////////////////////////

// Parameters for IW44 wavelet.

// - iw_norm: norm of all wavelets (for db estimation)

// - iw_shift: scale applied before decomposition

static const float iw_norm[16] = {

  2.627989e+03F,

  1.832893e+02F, 1.832959e+02F, 5.114690e+01F,

  4.583344e+01F, 4.583462e+01F, 1.279225e+01F,

  1.149671e+01F, 1.149712e+01F, 3.218888e+00F,

  2.999281e+00F, 2.999476e+00F, 8.733161e-01F,

  1.074451e+00F, 1.074511e+00F, 4.289318e-01F

};

static const int iw_shift  = 6;

static const int iw_round  = (1<<(iw_shift-1));

static const struct { int start; int size; }  

bandbuckets[] = 

{

  // Code first bucket and number of buckets in each band

  { 0, 1 }, // -- band zero contains all lores info

  { 1, 1 }, { 2, 1 }, { 3, 1 }, 

  { 4, 4 }, { 8, 4 }, { 12,4 }, 

  { 16,16 }, { 32,16 }, { 48,16 }, 

};

/** IW44 encoded gray-level image.  This class provided functions for managing

    a gray level image represented as a collection of IW44 wavelet

    coefficients.  The coefficients are stored in a memory efficient data

    structure.  Member function \Ref{get_bitmap} renders an arbitrary segment

    of the image into a \Ref{GBitmap}.  Member functions \Ref{decode_iff} and

    \Ref{encode_iff} read and write DjVu IW44 files (see \Ref{IW44Image.h}).

    Both the copy constructor and the copy operator are declared as private

    members. It is therefore not possible to make multiple copies of instances

    of this class. */

class IWBitmap::Encode : public IWBitmap

{

public:

  /// Destructor

  virtual ~Encode(void);

  /** Null constructor.  Constructs an empty IWBitmap object. This object does

      not contain anything meaningful. You must call function \Ref{init},

      \Ref{decode_iff} or \Ref{decode_chunk} to populate the wavelet

      coefficient data structure. */

  Encode(void);

  /** Initializes an IWBitmap with image #bm#.  This constructor

      performs the wavelet decomposition of image #bm# and records the

      corresponding wavelet coefficient.  Argument #mask# is an optional

      bilevel image specifying the masked pixels (see \Ref{IW44Image.h}). */

  void init(const GBitmap &bm, const GP<GBitmap> mask=0);

  // CODER

  /** Encodes one data chunk into ByteStream #bs#.  Parameter #parms# controls

      how much data is generated.  The chunk data is written to ByteStream

      #bs# with no IFF header.  Successive calls to #encode_chunk# encode

      successive chunks.  You must call #close_codec# after encoding the last

      chunk of a file. */

  virtual int  encode_chunk(GP<ByteStream> gbs, const IWEncoderParms &parms);

  /** Writes a gray level image into DjVu IW44 file.  This function creates a

      composite chunk (identifier #FORM:BM44#) composed of #nchunks# chunks

      (identifier #BM44#).  Data for each chunk is generated with

      #encode_chunk# using the corresponding parameters in array #parms#. */

  virtual void encode_iff(IFFByteStream &iff, int nchunks, const IWEncoderParms *parms);

  /** Resets the encoder/decoder state.  The first call to #decode_chunk# or

      #encode_chunk# initializes the coder for encoding or decoding.  Function

      #close_codec# must be called after processing the last chunk in order to

      reset the coder and release the associated memory. */

  virtual void close_codec(void);

protected:

  Codec::Encode *ycodec_enc;

};

/** IW44 encoded color image. This class provided functions for managing a

    color image represented as a collection of IW44 wavelet coefficients.  The

    coefficients are stored in a memory efficient data structure.  Member

    function \Ref{get_pixmap} renders an arbitrary segment of the image into a

    \Ref{GPixmap}.  Member functions \Ref{decode_iff} and \Ref{encode_iff}

    read and write DjVu IW44 files (see \Ref{IW44Image.h}).  Both the copy

    constructor and the copy operator are declared as private members. It is

    therefore not possible to make multiple copies of instances of this

    class. */

class IWPixmap::Encode : public IWPixmap

{

public:

  enum CRCBMode { 

    CRCBnone=IW44Image::CRCBnone, 

    CRCBhalf=IW44Image::CRCBhalf, 

    CRCBnormal=IW44Image::CRCBnormal, 

    CRCBfull=IW44Image::CRCBfull };

  /// Destructor

  virtual ~Encode(void);

  /** Null constructor.  Constructs an empty IWPixmap object. This object does

      not contain anything meaningful. You must call function \Ref{init},

      \Ref{decode_iff} or \Ref{decode_chunk} to populate the wavelet

      coefficient data structure. */

  Encode(void);

  /** Initializes an IWPixmap with color image #bm#.  This constructor

      performs the wavelet decomposition of image #bm# and records the

      corresponding wavelet coefficient.  Argument #mask# is an optional

      bilevel image specifying the masked pixels (see \Ref{IW44Image.h}).

      Argument #crcbmode# specifies how the chrominance information should be

      encoded (see \Ref{CRCBMode}). */

  void init(const GPixmap &bm, const GP<GBitmap> mask=0, CRCBMode crcbmode=CRCBnormal);

  // CODER

  /** Encodes one data chunk into ByteStream #bs#.  Parameter #parms# controls

      how much data is generated.  The chunk data is written to ByteStream

      #bs# with no IFF header.  Successive calls to #encode_chunk# encode

      successive chunks.  You must call #close_codec# after encoding the last

      chunk of a file. */

  virtual int  encode_chunk(GP<ByteStream> gbs, const IWEncoderParms &parms);

  /** Writes a color image into a DjVu IW44 file.  This function creates a

      composite chunk (identifier #FORM:PM44#) composed of #nchunks# chunks

      (identifier #PM44#).  Data for each chunk is generated with

      #encode_chunk# using the corresponding parameters in array #parms#. */

  virtual void encode_iff(IFFByteStream &iff, int nchunks, const IWEncoderParms *parms);

  /** Resets the encoder/decoder state.  The first call to #decode_chunk# or

      #encode_chunk# initializes the coder for encoding or decoding.  Function

      #close_codec# must be called after processing the last chunk in order to

      reset the coder and release the associated memory. */

  virtual void close_codec(void);

protected:

  Codec::Encode *ycodec_enc, *cbcodec_enc, *crcodec_enc;

};

class IW44Image::Map::Encode : public IW44Image::Map // DJVU_CLASS

{

public:

  Encode(const int w, const int h) : Map(w,h) {}

  // creation (from image)

  void create(const signed char *img8, int imgrowsize, 

              const signed char *msk8=0, int mskrowsize=0);

  // slash resolution

  void slashres(int res);

};

class IW44Image::Codec::Encode : public IW44Image::Codec

{

public:

  Encode(IW44Image::Map &map);

  // Coding

  virtual int code_slice(ZPCodec &zp;;

  float estimate_decibel(float frac);

  // Data

  void encode_buckets(ZPCodec &zp, int bit, int band,

    IW44Image::Block &blk, IW44Image::Block &eblk, int fbucket, int nbucket);

  int encode_prepare(int band, int fbucket, int nbucket, IW44Image::Block &blk, IW44Image::Block &eblk);

  IW44Image::Map emap;

};

IW44Image::Codec::Encode::Encode(IW44Image::Map &map)

: Codec(map), emap(map.iw,map.ih) {}

//////////////////////////////////////////////////////

/** IW44Image::Transform::Encode

*/

class IW44Image::Transform::Encode : IW44Image::Transform

{

 public:

 // WAVELET TRANSFORM

  /** Forward transform. */

  static void forward(short *p, int w, int h, int rowsize, int begin, int end);

  // COLOR TRANSFORM

  /** Extracts Y */

  static void RGB_to_Y(const GPixel *p, int w, int h, int rowsize, 

                       signed char *out, int outrowsize);

  /** Extracts Cb */

  static void RGB_to_Cb(const GPixel *p, int w, int h, int rowsize, 

                        signed char *out, int outrowsize);

  /** Extracts Cr */

  static void RGB_to_Cr(const GPixel *p, int w, int h, int rowsize, 

                        signed char *out, int outrowsize);

};

//////////////////////////////////////////////////////

// MMX IMPLEMENTATION HELPERS

//////////////////////////////////////////////////////

// Note:

// MMX implementation for vertical transforms only.

// Speedup is basically related to faster memory transfer

// The IW44 transform is not CPU bound, it is memory bound.

#ifdef MMX

static const short w9[]  = {9,9,9,9};

static const short w1[]  = {1,1,1,1};

static const int   d8[]  = {8,8};

static const int   d16[] = {16,16};

static inline void

mmx_fv_1 ( short* &q, short* e, int s, int s3 )

{

  while (q

    {

      int a = (int)q[-s] + (int)q[s];

      int b = (int)q[-s3] + (int)q[s3];

      *q -= (((a<<3)+a-b+8)>>4);

      q++;

    }

  while (q+3

    {

      MMXar( movq,       q-s,mm0);  // MM0=[ b3, b2, b1, b0 ]

      MMXar( movq,       q+s,mm2);  // MM2=[ c3, c2, c1, c0 ]

      MMXrr( movq,       mm0,mm1);  

      MMXrr( punpcklwd,  mm2,mm0);  // MM0=[ c1, b1, c0, b0 ]

      MMXrr( punpckhwd,  mm2,mm1);  // MM1=[ c3, b3, c2, b2 ]

      MMXar( pmaddwd,    w9,mm0);   // MM0=[ (c1+b1)*9, (c0+b0)*9 ]

      MMXar( pmaddwd,    w9,mm1);   // MM1=[ (c3+b3)*9, (c2+b2)*9 ]

      MMXar( movq,       q-s3,mm2);

      MMXar( movq,       q+s3,mm4);

      MMXrr( movq,       mm2,mm3);

      MMXrr( punpcklwd,  mm4,mm2);  // MM2=[ d1, a1, d0, a0 ]

      MMXrr( punpckhwd,  mm4,mm3);  // MM3=[ d3, a3, d2, a2 ]

      MMXar( pmaddwd,    w1,mm2);   // MM2=[ (a1+d1)*1, (a0+d0)*1 ]

      MMXar( pmaddwd,    w1,mm3);   // MM3=[ (a3+d3)*1, (a2+d2)*1 ]

      MMXar( paddd,      d8,mm0);

      MMXar( paddd,      d8,mm1);

      MMXrr( psubd,      mm2,mm0);  // MM0=[ (c1+b1)*9-a1-d1+8, ...

      MMXrr( psubd,      mm3,mm1);  // MM1=[ (c3+b3)*9-a3-d3+8, ...

      MMXir( psrad,      4,mm0);

      MMXar( movq,       q,mm7);    // MM7=[ p3,p2,p1,p0 ]

      MMXir( psrad,      4,mm1);

      MMXrr( packssdw,   mm1,mm0);  // MM0=[ x3,x2,x1,x0 ]

      MMXrr( psubw,      mm0,mm7);  // MM7=[ p3-x3, p2-x2, ... ]

      MMXra( movq,       mm7,q);

#if defined(_MSC_VER) && defined(_DEBUG)

      MMXemms;

#endif

      q += 4;

    }

}

static inline void

mmx_fv_2 ( short* &q, short* e, int s, int s3 )

{

  while (q

    {

      int a = (int)q[-s] + (int)q[s];

      int b = (int)q[-s3] + (int)q[s3];

      *q += (((a<<3)+a-b+16)>>5);

      q ++;

    }

  while (q+3

    {

      MMXar( movq,       q-s,mm0);  // MM0=[ b3, b2, b1, b0 ]

      MMXar( movq,       q+s,mm2);  // MM2=[ c3, c2, c1, c0 ]

      MMXrr( movq,       mm0,mm1);  

      MMXrr( punpcklwd,  mm2,mm0);  // MM0=[ c1, b1, c0, b0 ]

      MMXrr( punpckhwd,  mm2,mm1);  // MM1=[ c3, b3, c2, b2 ]

      MMXar( pmaddwd,    w9,mm0);   // MM0=[ (c1+b1)*9, (c0+b0)*9 ]

      MMXar( pmaddwd,    w9,mm1);   // MM1=[ (c3+b3)*9, (c2+b2)*9 ]

      MMXar( movq,       q-s3,mm2);

      MMXar( movq,       q+s3,mm4);

      MMXrr( movq,       mm2,mm3);

      MMXrr( punpcklwd,  mm4,mm2);  // MM2=[ d1, a1, d0, a0 ]

      MMXrr( punpckhwd,  mm4,mm3);  // MM3=[ d3, a3, d2, a2 ]

      MMXar( pmaddwd,    w1,mm2);   // MM2=[ (a1+d1)*1, (a0+d0)*1 ]

      MMXar( pmaddwd,    w1,mm3);   // MM3=[ (a3+d3)*1, (a2+d2)*1 ]

      MMXar( paddd,      d16,mm0);

      MMXar( paddd,      d16,mm1);

      MMXrr( psubd,      mm2,mm0);  // MM0=[ (c1+b1)*9-a1-d1+8, ...

      MMXrr( psubd,      mm3,mm1);  // MM1=[ (c3+b3)*9-a3-d3+8, ...

      MMXir( psrad,      5,mm0);

      MMXar( movq,       q,mm7);    // MM7=[ p3,p2,p1,p0 ]

      MMXir( psrad,      5,mm1);

      MMXrr( packssdw,   mm1,mm0);  // MM0=[ x3,x2,x1,x0 ]

      MMXrr( paddw,      mm0,mm7);  // MM7=[ p3+x3, p2+x2, ... ]

      MMXra( movq,       mm7,q);

#if defined(_MSC_VER) && defined(_DEBUG)

      MMXemms;

#endif

      q += 4;

    }

}

#endif /* MMX */

//////////////////////////////////////////////////////

// NEW FILTERS

//////////////////////////////////////////////////////

static void 

filter_fv(short *p, int w, int h, int rowsize, int scale)

{

  int y = 0;

  int s = scale*rowsize;

  int s3 = s+s+s;

  h = ((h-1)/scale)+1;

  y += 1;

  p += s;

  while (y-3 < h)

    {

      // 1-Delta

      {

        short *q = p;

        short *e = q+w;

        if (y>=3 && y+3

          {

            // Generic case

#ifdef MMX

            if (scale==1 && MMXControl::mmxflag>0)

              mmx_fv_1(q, e, s, s3);

#endif

            while (q

              {

                int a = (int)q[-s] + (int)q[s];

                int b = (int)q[-s3] + (int)q[s3];

                *q -= (((a<<3)+a-b+8)>>4);

                q += scale;

              }

          }

        else if (y

          {

            // Special cases

            short *q1 = (y+1

            while (q

              {

                int a = (int)q[-s] + (int)(*q1);

                *q -= ((a+1)>>1);

                q += scale;

                q1 += scale;

              }

          }

      }

      // 2-Update

      {

        short *q = p-s3;

        short *e = q+w;

        if (y>=6 && y

          {

            // Generic case

#ifdef MMX

            if (scale==1 && MMXControl::mmxflag>0)

              mmx_fv_2(q, e, s, s3);

#endif

            while (q

              {

                int a = (int)q[-s] + (int)q[s];

                int b = (int)q[-s3] + (int)q[s3];

                *q += (((a<<3)+a-b+16)>>5);

                q += scale;

              }

          }

        else if (y>=3)

          {

            // Special cases

            short *q1 = (y-2

            short *q3 = (y

            if (y>=6)

              {

                while (q

                  {

                    int a = (int)q[-s] + (q1 ? (int)(*q1) : 0);

                    int b = (int)q[-s3] + (q3 ? (int)(*q3) : 0);

                    *q += (((a<<3)+a-b+16)>>5);

                    q += scale;

                    if (q1) q1 += scale;

                    if (q3) q3 += scale;

                  }

              }

            else if (y>=4)

              {

                while (q

                  {

                    int a = (int)q[-s] + (q1 ? (int)(*q1) : 0);

                    int b = (q3 ? (int)(*q3) : 0);

                    *q += (((a<<3)+a-b+16)>>5);

                    q += scale;

                    if (q1) q1 += scale;

                    if (q3) q3 += scale;

                  }

              }

            else

              {

                while (q

                  {

                    int a = (q1 ? (int)(*q1) : 0);

                    int b = (q3 ? (int)(*q3) : 0);

                    *q += (((a<<3)+a-b+16)>>5);

                    q += scale;

                    if (q1) q1 += scale;

                    if (q3) q3 += scale;

                  }

              }

          }

      }

      y += 2;

      p += s+s;

    }

}

static void 

filter_fh(short *p, int w, int h, int rowsize, int scale)

{

  int y = 0;

  int s = scale;

  int s3 = s+s+s;

  rowsize *= scale;

  while (y

    {

      short *q = p+s;

      short *e = p+w;

      int a0=0, a1=0, a2=0, a3=0;

      int b0=0, b1=0, b2=0, b3=0;

      if (q < e)

        {

          // Special case: x=1

          a1 = a2 = a3 = q[-s];

          if (q+s

            a2 = q[s];

          if (q+s3

            a3 = q[s3];

          b3 = q[0] - ((a1+a2+1)>>1);

          q[0] = b3;

          q += s+s;

        }

      while (q+s3 < e)

        {

          // Generic case

          a0=a1; 

          a1=a2; 

          a2=a3;

          a3=q[s3];

          b0=b1; 

          b1=b2; 

          b2=b3;

          b3 = q[0] - ((((a1+a2)<<3)+(a1+a2)-a0-a3+8) >> 4);

          q[0] = b3;

          q[-s3] = q[-s3] + ((((b1+b2)<<3)+(b1+b2)-b0-b3+16) >> 5);

          q += s+s;

        }

      while (q < e)

        {

          // Special case: w-3 <= x < w

          a1=a2; 

          a2=a3;

          b0=b1; 

          b1=b2; 

          b2=b3;

          b3 = q[0] - ((a1+a2+1)>>1);

          q[0] = b3;

          q[-s3] = q[-s3] + ((((b1+b2)<<3)+(b1+b2)-b0-b3+16) >> 5);

          q += s+s;

        }

      while (q-s3 < e)

        {

          // Special case  w <= x < w+3

          b0=b1; 

          b1=b2; 

          b2=b3;

          b3=0;

          if (q-s3 >= p)

            q[-s3] = q[-s3] + ((((b1+b2)<<3)+(b1+b2)-b0-b3+16) >> 5);

          q += s+s;

        }

      y += scale;

      p += rowsize;

    }

}

//////////////////////////////////////////////////////

// WAVELET TRANSFORM 

//////////////////////////////////////////////////////

//----------------------------------------------------

// Function for applying bidimensional IW44 between 

// scale intervals begin(inclusive) and end(exclusive)

void

IW44Image::Transform::Encode::forward(short *p, int w, int h, int rowsize, int begin, int end)

{ 

  // PREPARATION

  filter_begin(w,h);

  // LOOP ON SCALES

  for (int scale=begin; scale

    {

#ifdef IWTRANSFORM_TIMER

      int tv,th;

      th = tv = GOS::ticks();

#endif

      filter_fh(p, w, h, rowsize, scale);

#ifdef IWTRANSFORM_TIMER

      th = GOS::ticks();

      tv = th - tv;

#endif

      filter_fv(p, w, h, rowsize, scale);

#ifdef IWTRANSFORM_TIMER

      th = GOS::ticks()-th;

      DjVuPrintErrorUTF8("forw%d\tv=%dms h=%dms\n", scale,th,tv);

#endif

    }

  // TERMINATE

  filter_end();

}

//////////////////////////////////////////////////////

// COLOR TRANSFORM 

//////////////////////////////////////////////////////

static const float 

rgb_to_ycc[3][3] = 

{ { 0.304348F,  0.608696F,  0.086956F },      

  { 0.463768F, -0.405797F, -0.057971F },

  {-0.173913F, -0.347826F,  0.521739F } };

/* Extracts Y */

void 

IW44Image::Transform::Encode::RGB_to_Y(const GPixel *p, int w, int h, int rowsize, 

                      signed char *out, int outrowsize)

{

  int rmul[256], gmul[256], bmul[256];

  for (int k=0; k<256; k++)

    {

      rmul[k] = (int)(k*0x10000*rgb_to_ycc[0][0]);

      gmul[k] = (int)(k*0x10000*rgb_to_ycc[0][1]);

      bmul[k] = (int)(k*0x10000*rgb_to_ycc[0][2]);

    }

  for (int i=0; i

    {

      const GPixel *p2 = p;

      signed char *out2 = out;

      for (int j=0; j

        {

          int y = rmul[p2->r] + gmul[p2->g] + bmul[p2->b] + 32768;

          *out2 = (y>>16) - 128;

        }

    }

}

#ifdef min

#undef min

#endif

static inline int min(const int x,const int y) {return (x

#ifdef max

#undef max

#endif

static inline int max(const int x,const int y) {return (x>y)?x:y;}

/* Extracts Cb */

void 

IW44Image::Transform::Encode::RGB_to_Cb(const GPixel *p, int w, int h, int rowsize, 

                       signed char *out, int outrowsize)

{

  int rmul[256], gmul[256], bmul[256];

  for (int k=0; k<256; k++)

    {

      rmul[k] = (int)(k*0x10000*rgb_to_ycc[2][0]);

      gmul[k] = (int)(k*0x10000*rgb_to_ycc[2][1]);

      bmul[k] = (int)(k*0x10000*rgb_to_ycc[2][2]);

    }

  for (int i=0; i

    {

      const GPixel *p2 = p;

      signed char *out2 = out;

      for (int j=0; j

        {

          int c = rmul[p2->r] + gmul[p2->g] + bmul[p2->b] + 32768;

          *out2 = max(-128, min(127, c>>16));

        }

    }

}

/* Extracts Cr */

void 

IW44Image::Transform::Encode::RGB_to_Cr(const GPixel *p, int w, int h, int rowsize, 

                       signed char *out, int outrowsize)

{

  int rmul[256], gmul[256], bmul[256];

  for (int k=0; k<256; k++)

    {

      rmul[k] = (int)((k*0x10000)*rgb_to_ycc[1][0]);

      gmul[k] = (int)((k*0x10000)*rgb_to_ycc[1][1]);

      bmul[k] = (int)((k*0x10000)*rgb_to_ycc[1][2]);

    }

  for (int i=0; i

    {

      const GPixel *p2 = p;

      signed char *out2 = out;

      for (int j=0; j

        {

          int c = rmul[p2->r] + gmul[p2->g] + bmul[p2->b] + 32768;

          *out2 = max(-128, min(127, c>>16));

        }

    }

}

//////////////////////////////////////////////////////

// MASKING DECOMPOSITION

//////////////////////////////////////////////////////

//----------------------------------------------------

// Function for applying bidimensional IW44 between 

// scale intervals begin(inclusive) and end(exclusive)

// with a MASK bitmap

static void

interpolate_mask(short *data16, int w, int h, int rowsize,

                 const signed char *mask8, int mskrowsize)

{

  int i,j;

  // count masked bits

  short *count;

  GPBuffer<short> gcount(count,w*h);

  short *cp = count;

  for (i=0; i

    for (j=0; j

      cp[j] = (mask8[j] ? 0 : 0x1000);

  // copy image

  short *sdata;

  GPBuffer<short> gsdata(sdata,w*h);

  short *p = sdata;

  short *q = data16;

  for (i=0; i

    for (j=0; j

      p[j] = q[j];

  // iterate over resolutions

  int split = 1;

  int scale = 2;

  int again = 1;

  while (again && scale

    {

      again = 0;

      p = data16;

      q = sdata;

      cp = count;

      // iterate over block

      for (i=0; i

        for (j=0; j

          {

            int ii, jj;

            int gotz = 0;

            int gray = 0;

            int npix = 0;

            short *cpp = cp;

            short *qq = q;

            // look around when square goes beyond border

            int istart = i;

            if (istart+split>h)

              {

                istart -= scale;

                cpp -= w*scale;

                qq -= w*scale;

              }

            int jstart = j;

            if (jstart+split>w)

              jstart -= scale;

            // compute gray level

            for (ii=istart; ii

              for (jj=jstart; jj

                {

                  if (cpp[jj]>0) 

                    {

                      npix += cpp[jj];

                      gray += cpp[jj] * qq[jj];

                    } 

                  else if (ii>=i && jj>=j)

                    {

                      gotz = 1;

                    }

                }

            // process result

            if (npix == 0)

              {

                // continue to next resolution

                again = 1;

                cp[j] = 0;

              }

            else

              {

                gray = gray / npix;

                // check whether initial image require fix

                if (gotz)

                  {

                    cpp = cp;

                    qq = p;

                    for (ii=i; ii

                      for (jj=j; jj

                        if (cpp[jj] == 0)

                          {

                            qq[jj] = gray;

                            cpp[jj] = 1;

                          }

                  }

                // store average for next iteration

                cp[j] = npix>>2;

                q[j] = gray;

              }

          }

      // double resolution

      split = scale;

      scale = scale+scale;

    }

}

static void

forward_mask(short *data16, int w, int h, int rowsize, int begin, int end,

             const signed char *mask8, int mskrowsize )

{

  int i,j;

  signed char *m;

  short *p;

  short *d;

  // Allocate buffers

  short *sdata;

  GPBuffer<short> gsdata(sdata,w*h);

  signed char *smask;

  GPBuffer<signed char> gsmask(smask,w*h);

  // Copy mask

  m = smask;

  for (i=0; i

    memcpy((void*)m, (void*)mask8, w);

  // Loop over scale

  for (int scale=begin; scale

    {

      // Copy data into sdata buffer

      p = data16;

      d = sdata;

      for (i=0; i

        {

          for (j=0; j

            d[j] = p[j];

          p += rowsize * scale;

          d += w * scale;

        }

      // Decompose

      IW44Image::Transform::Encode::forward(sdata, w, h, w, scale, scale+scale);

      // Cancel masked coefficients

      d = sdata;

      m = smask;

      for (i=0; i

        {

          for (j=scale; j

            if (m[j])

              d[j] = 0;

          d += w * scale;

          m += w * scale;

          if (i+scale < h)

            {

              for (j=0; j

                if (m[j])

                  d[j] = 0;

              d += w * scale;

              m += w * scale;

            }

        }

      // Reconstruct

      IW44Image::Transform::Decode::backward(sdata, w, h, w, scale+scale, scale);

      // Correct visible pixels

      p = data16;

      d = sdata;

      m = smask;

      for (i=0; i

        {