source: trunk/libdjvu/IW44EncodeCodec.cpp @ 426

Last change on this file since 426 was 280, checked in by rbri, 12 years ago

DJVU plugin: djvulibre updated to version 3.5.22

File size: 52.9 KB
Line 
1//C-  -*- C++ -*-
2//C- -------------------------------------------------------------------
3//C- DjVuLibre-3.5
4//C- Copyright (c) 2002  Leon Bottou and Yann Le Cun.
5//C- Copyright (c) 2001  AT&T
6//C-
7//C- This software is subject to, and may be distributed under, the
8//C- GNU General Public License, either Version 2 of the license,
9//C- or (at your option) any later version. The license should have
10//C- accompanied the software or you may obtain a copy of the license
11//C- from the Free Software Foundation at http://www.fsf.org .
12//C-
13//C- This program is distributed in the hope that it will be useful,
14//C- but WITHOUT ANY WARRANTY; without even the implied warranty of
15//C- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16//C- GNU General Public License for more details.
17//C-
18//C- DjVuLibre-3.5 is derived from the DjVu(r) Reference Library from
19//C- Lizardtech Software.  Lizardtech Software has authorized us to
20//C- replace the original DjVu(r) Reference Library notice by the following
21//C- text (see doc/lizard2002.djvu and doc/lizardtech2007.djvu):
22//C-
23//C-  ------------------------------------------------------------------
24//C- | DjVu (r) Reference Library (v. 3.5)
25//C- | Copyright (c) 1999-2001 LizardTech, Inc. All Rights Reserved.
26//C- | The DjVu Reference Library is protected by U.S. Pat. No.
27//C- | 6,058,214 and patents pending.
28//C- |
29//C- | This software is subject to, and may be distributed under, the
30//C- | GNU General Public License, either Version 2 of the license,
31//C- | or (at your option) any later version. The license should have
32//C- | accompanied the software or you may obtain a copy of the license
33//C- | from the Free Software Foundation at http://www.fsf.org .
34//C- |
35//C- | The computer code originally released by LizardTech under this
36//C- | license and unmodified by other parties is deemed "the LIZARDTECH
37//C- | ORIGINAL CODE."  Subject to any third party intellectual property
38//C- | claims, LizardTech grants recipient a worldwide, royalty-free,
39//C- | non-exclusive license to make, use, sell, or otherwise dispose of
40//C- | the LIZARDTECH ORIGINAL CODE or of programs derived from the
41//C- | LIZARDTECH ORIGINAL CODE in compliance with the terms of the GNU
42//C- | General Public License.   This grant only confers the right to
43//C- | infringe patent claims underlying the LIZARDTECH ORIGINAL CODE to
44//C- | the extent such infringement is reasonably necessary to enable
45//C- | recipient to make, have made, practice, sell, or otherwise dispose
46//C- | of the LIZARDTECH ORIGINAL CODE (or portions thereof) and not to
47//C- | any greater extent that may be necessary to utilize further
48//C- | modifications or combinations.
49//C- |
50//C- | The LIZARDTECH ORIGINAL CODE is provided "AS IS" WITHOUT WARRANTY
51//C- | OF ANY KIND, EITHER EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED
52//C- | TO ANY WARRANTY OF NON-INFRINGEMENT, OR ANY IMPLIED WARRANTY OF
53//C- | MERCHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE.
54//C- +------------------------------------------------------------------
55//
56// $Id: IW44EncodeCodec.cpp,v 1.12 2007/03/25 20:48:32 leonb Exp $
57// $Name: release_3_5_22 $
58
59#ifdef HAVE_CONFIG_H
60# include "config.h"
61#endif
62#if NEED_GNUG_PRAGMAS
63# pragma implementation
64#endif
65
66// - Author: Leon Bottou, 08/1998
67// From: Leon Bottou, 1/31/2002
68// Lizardtech has split this file into a decoder and an encoder.
69// Only superficial changes.  The meat is mine.
70
71#define IW44IMAGE_IMPLIMENTATION /* */
72
73
74
75#include "IW44Image.h"
76#include "ZPCodec.h"
77#include "GBitmap.h"
78#include "GPixmap.h"
79#include "IFFByteStream.h"
80#include "GRect.h"
81
82#include <stdlib.h>
83#include <string.h>
84#include <math.h>
85#include "MMX.h"
86#undef IWTRANSFORM_TIMER
87#ifdef IWTRANSFORM_TIMER
88#include "GOS.h"
89#endif
90
91#include <assert.h>
92#include <string.h>
93#include <math.h>
94
95#ifndef NEED_DECODER_ONLY
96
97
98#ifdef HAVE_NAMESPACES
99namespace DJVU {
100# ifdef NOT_DEFINED // Just to fool emacs c++ mode
101}
102#endif
103#endif
104
105#define IWALLOCSIZE    4080
106#define IWCODEC_MAJOR     1
107#define IWCODEC_MINOR     2
108#define DECIBEL_PRUNE   5.0
109
110
111//////////////////////////////////////////////////////
112// WAVELET DECOMPOSITION CONSTANTS
113//////////////////////////////////////////////////////
114
115// Parameters for IW44 wavelet.
116// - iw_norm: norm of all wavelets (for db estimation)
117// - iw_shift: scale applied before decomposition
118
119
120static const float iw_norm[16] = {
121  2.627989e+03F,
122  1.832893e+02F, 1.832959e+02F, 5.114690e+01F,
123  4.583344e+01F, 4.583462e+01F, 1.279225e+01F,
124  1.149671e+01F, 1.149712e+01F, 3.218888e+00F,
125  2.999281e+00F, 2.999476e+00F, 8.733161e-01F,
126  1.074451e+00F, 1.074511e+00F, 4.289318e-01F
127};
128
129static const int iw_shift  = 6;
130static const int iw_round  = (1<<(iw_shift-1));
131
132static const struct { int start; int size; } 
133bandbuckets[] = 
134{
135  // Code first bucket and number of buckets in each band
136  { 0, 1 }, // -- band zero contains all lores info
137  { 1, 1 }, { 2, 1 }, { 3, 1 }, 
138  { 4, 4 }, { 8, 4 }, { 12,4 }, 
139  { 16,16 }, { 32,16 }, { 48,16 }, 
140};
141
142
143/** IW44 encoded gray-level image.  This class provided functions for managing
144    a gray level image represented as a collection of IW44 wavelet
145    coefficients.  The coefficients are stored in a memory efficient data
146    structure.  Member function \Ref{get_bitmap} renders an arbitrary segment
147    of the image into a \Ref{GBitmap}.  Member functions \Ref{decode_iff} and
148    \Ref{encode_iff} read and write DjVu IW44 files (see \Ref{IW44Image.h}).
149    Both the copy constructor and the copy operator are declared as private
150    members. It is therefore not possible to make multiple copies of instances
151    of this class. */
152
153class IWBitmap::Encode : public IWBitmap
154{
155public:
156  /// Destructor
157  virtual ~Encode(void);
158  /** Null constructor.  Constructs an empty IWBitmap object. This object does
159      not contain anything meaningful. You must call function \Ref{init},
160      \Ref{decode_iff} or \Ref{decode_chunk} to populate the wavelet
161      coefficient data structure. */
162  Encode(void);
163  /** Initializes an IWBitmap with image #bm#.  This constructor
164      performs the wavelet decomposition of image #bm# and records the
165      corresponding wavelet coefficient.  Argument #mask# is an optional
166      bilevel image specifying the masked pixels (see \Ref{IW44Image.h}). */
167  void init(const GBitmap &bm, const GP<GBitmap> mask=0);
168  // CODER
169  /** Encodes one data chunk into ByteStream #bs#.  Parameter #parms# controls
170      how much data is generated.  The chunk data is written to ByteStream
171      #bs# with no IFF header.  Successive calls to #encode_chunk# encode
172      successive chunks.  You must call #close_codec# after encoding the last
173      chunk of a file. */
174  virtual int  encode_chunk(GP<ByteStream> gbs, const IWEncoderParms &parms);
175  /** Writes a gray level image into DjVu IW44 file.  This function creates a
176      composite chunk (identifier #FORM:BM44#) composed of #nchunks# chunks
177      (identifier #BM44#).  Data for each chunk is generated with
178      #encode_chunk# using the corresponding parameters in array #parms#. */
179  virtual void encode_iff(IFFByteStream &iff, int nchunks, const IWEncoderParms *parms);
180  /** Resets the encoder/decoder state.  The first call to #decode_chunk# or
181      #encode_chunk# initializes the coder for encoding or decoding.  Function
182      #close_codec# must be called after processing the last chunk in order to
183      reset the coder and release the associated memory. */
184  virtual void close_codec(void);
185protected:
186  Codec::Encode *ycodec_enc;
187};
188
189/** IW44 encoded color image. This class provided functions for managing a
190    color image represented as a collection of IW44 wavelet coefficients.  The
191    coefficients are stored in a memory efficient data structure.  Member
192    function \Ref{get_pixmap} renders an arbitrary segment of the image into a
193    \Ref{GPixmap}.  Member functions \Ref{decode_iff} and \Ref{encode_iff}
194    read and write DjVu IW44 files (see \Ref{IW44Image.h}).  Both the copy
195    constructor and the copy operator are declared as private members. It is
196    therefore not possible to make multiple copies of instances of this
197    class. */
198
199class IWPixmap::Encode : public IWPixmap
200{
201public:
202  enum CRCBMode { 
203    CRCBnone=IW44Image::CRCBnone, 
204    CRCBhalf=IW44Image::CRCBhalf, 
205    CRCBnormal=IW44Image::CRCBnormal, 
206    CRCBfull=IW44Image::CRCBfull };
207  /// Destructor
208  virtual ~Encode(void);
209  /** Null constructor.  Constructs an empty IWPixmap object. This object does
210      not contain anything meaningful. You must call function \Ref{init},
211      \Ref{decode_iff} or \Ref{decode_chunk} to populate the wavelet
212      coefficient data structure. */
213  Encode(void);
214  /** Initializes an IWPixmap with color image #bm#.  This constructor
215      performs the wavelet decomposition of image #bm# and records the
216      corresponding wavelet coefficient.  Argument #mask# is an optional
217      bilevel image specifying the masked pixels (see \Ref{IW44Image.h}).
218      Argument #crcbmode# specifies how the chrominance information should be
219      encoded (see \Ref{CRCBMode}). */
220  void init(const GPixmap &bm, const GP<GBitmap> mask=0, CRCBMode crcbmode=CRCBnormal);
221  // CODER
222  /** Encodes one data chunk into ByteStream #bs#.  Parameter #parms# controls
223      how much data is generated.  The chunk data is written to ByteStream
224      #bs# with no IFF header.  Successive calls to #encode_chunk# encode
225      successive chunks.  You must call #close_codec# after encoding the last
226      chunk of a file. */
227  virtual int  encode_chunk(GP<ByteStream> gbs, const IWEncoderParms &parms);
228  /** Writes a color image into a DjVu IW44 file.  This function creates a
229      composite chunk (identifier #FORM:PM44#) composed of #nchunks# chunks
230      (identifier #PM44#).  Data for each chunk is generated with
231      #encode_chunk# using the corresponding parameters in array #parms#. */
232  virtual void encode_iff(IFFByteStream &iff, int nchunks, const IWEncoderParms *parms);
233  /** Resets the encoder/decoder state.  The first call to #decode_chunk# or
234      #encode_chunk# initializes the coder for encoding or decoding.  Function
235      #close_codec# must be called after processing the last chunk in order to
236      reset the coder and release the associated memory. */
237  virtual void close_codec(void);
238protected:
239  Codec::Encode *ycodec_enc, *cbcodec_enc, *crcodec_enc;
240};
241
242class IW44Image::Map::Encode : public IW44Image::Map // DJVU_CLASS
243{
244public:
245  Encode(const int w, const int h) : Map(w,h) {}
246  // creation (from image)
247  void create(const signed char *img8, int imgrowsize, 
248              const signed char *msk8=0, int mskrowsize=0);
249  // slash resolution
250  void slashres(int res);
251};
252
253class IW44Image::Codec::Encode : public IW44Image::Codec
254{
255public:
256  Encode(IW44Image::Map &map);
257  // Coding
258  virtual int code_slice(ZPCodec &zp);
259  float estimate_decibel(float frac);
260  // Data
261  void encode_buckets(ZPCodec &zp, int bit, int band,
262    IW44Image::Block &blk, IW44Image::Block &eblk, int fbucket, int nbucket);
263  int encode_prepare(int band, int fbucket, int nbucket, IW44Image::Block &blk, IW44Image::Block &eblk);
264  IW44Image::Map emap;
265};
266
267IW44Image::Codec::Encode::Encode(IW44Image::Map &map)
268: Codec(map), emap(map.iw,map.ih) {}
269
270//////////////////////////////////////////////////////
271/** IW44Image::Transform::Encode
272*/
273
274class IW44Image::Transform::Encode : IW44Image::Transform
275{
276 public:
277 // WAVELET TRANSFORM
278  /** Forward transform. */
279  static void forward(short *p, int w, int h, int rowsize, int begin, int end);
280 
281  // COLOR TRANSFORM
282  /** Extracts Y */
283  static void RGB_to_Y(const GPixel *p, int w, int h, int rowsize, 
284                       signed char *out, int outrowsize);
285  /** Extracts Cb */
286  static void RGB_to_Cb(const GPixel *p, int w, int h, int rowsize, 
287                        signed char *out, int outrowsize);
288  /** Extracts Cr */
289  static void RGB_to_Cr(const GPixel *p, int w, int h, int rowsize, 
290                        signed char *out, int outrowsize);
291};
292
293
294//////////////////////////////////////////////////////
295// MMX IMPLEMENTATION HELPERS
296//////////////////////////////////////////////////////
297
298
299// Note:
300// MMX implementation for vertical transforms only.
301// Speedup is basically related to faster memory transfer
302// The IW44 transform is not CPU bound, it is memory bound.
303
304#ifdef MMX
305
306static const short w9[]  = {9,9,9,9};
307static const short w1[]  = {1,1,1,1};
308static const int   d8[]  = {8,8};
309static const int   d16[] = {16,16};
310
311
312static inline void
313mmx_fv_1 ( short* &q, short* e, int s, int s3 )
314{
315  while (q<e && (((long)q)&0x7))
316    {
317      int a = (int)q[-s] + (int)q[s];
318      int b = (int)q[-s3] + (int)q[s3];
319      *q -= (((a<<3)+a-b+8)>>4);
320      q++;
321    }
322  while (q+3<e)
323    {
324      MMXar( movq,       q-s,mm0);  // MM0=[ b3, b2, b1, b0 ]
325      MMXar( movq,       q+s,mm2);  // MM2=[ c3, c2, c1, c0 ]
326      MMXrr( movq,       mm0,mm1); 
327      MMXrr( punpcklwd,  mm2,mm0);  // MM0=[ c1, b1, c0, b0 ]
328      MMXrr( punpckhwd,  mm2,mm1);  // MM1=[ c3, b3, c2, b2 ]
329      MMXar( pmaddwd,    w9,mm0);   // MM0=[ (c1+b1)*9, (c0+b0)*9 ]
330      MMXar( pmaddwd,    w9,mm1);   // MM1=[ (c3+b3)*9, (c2+b2)*9 ]
331      MMXar( movq,       q-s3,mm2);
332      MMXar( movq,       q+s3,mm4);
333      MMXrr( movq,       mm2,mm3);
334      MMXrr( punpcklwd,  mm4,mm2);  // MM2=[ d1, a1, d0, a0 ]
335      MMXrr( punpckhwd,  mm4,mm3);  // MM3=[ d3, a3, d2, a2 ]
336      MMXar( pmaddwd,    w1,mm2);   // MM2=[ (a1+d1)*1, (a0+d0)*1 ]
337      MMXar( pmaddwd,    w1,mm3);   // MM3=[ (a3+d3)*1, (a2+d2)*1 ]
338      MMXar( paddd,      d8,mm0);
339      MMXar( paddd,      d8,mm1);
340      MMXrr( psubd,      mm2,mm0);  // MM0=[ (c1+b1)*9-a1-d1+8, ...
341      MMXrr( psubd,      mm3,mm1);  // MM1=[ (c3+b3)*9-a3-d3+8, ...
342      MMXir( psrad,      4,mm0);
343      MMXar( movq,       q,mm7);    // MM7=[ p3,p2,p1,p0 ]
344      MMXir( psrad,      4,mm1);
345      MMXrr( packssdw,   mm1,mm0);  // MM0=[ x3,x2,x1,x0 ]
346      MMXrr( psubw,      mm0,mm7);  // MM7=[ p3-x3, p2-x2, ... ]
347      MMXra( movq,       mm7,q);
348#if defined(_MSC_VER) && defined(_DEBUG)
349      MMXemms;
350#endif
351      q += 4;
352    }
353}
354
355static inline void
356mmx_fv_2 ( short* &q, short* e, int s, int s3 )
357{
358  while (q<e && (((long)q)&0x7))
359    {
360      int a = (int)q[-s] + (int)q[s];
361      int b = (int)q[-s3] + (int)q[s3];
362      *q += (((a<<3)+a-b+16)>>5);
363      q ++;
364    }
365  while (q+3<e)
366    {
367      MMXar( movq,       q-s,mm0);  // MM0=[ b3, b2, b1, b0 ]
368      MMXar( movq,       q+s,mm2);  // MM2=[ c3, c2, c1, c0 ]
369      MMXrr( movq,       mm0,mm1); 
370      MMXrr( punpcklwd,  mm2,mm0);  // MM0=[ c1, b1, c0, b0 ]
371      MMXrr( punpckhwd,  mm2,mm1);  // MM1=[ c3, b3, c2, b2 ]
372      MMXar( pmaddwd,    w9,mm0);   // MM0=[ (c1+b1)*9, (c0+b0)*9 ]
373      MMXar( pmaddwd,    w9,mm1);   // MM1=[ (c3+b3)*9, (c2+b2)*9 ]
374      MMXar( movq,       q-s3,mm2);
375      MMXar( movq,       q+s3,mm4);
376      MMXrr( movq,       mm2,mm3);
377      MMXrr( punpcklwd,  mm4,mm2);  // MM2=[ d1, a1, d0, a0 ]
378      MMXrr( punpckhwd,  mm4,mm3);  // MM3=[ d3, a3, d2, a2 ]
379      MMXar( pmaddwd,    w1,mm2);   // MM2=[ (a1+d1)*1, (a0+d0)*1 ]
380      MMXar( pmaddwd,    w1,mm3);   // MM3=[ (a3+d3)*1, (a2+d2)*1 ]
381      MMXar( paddd,      d16,mm0);
382      MMXar( paddd,      d16,mm1);
383      MMXrr( psubd,      mm2,mm0);  // MM0=[ (c1+b1)*9-a1-d1+8, ...
384      MMXrr( psubd,      mm3,mm1);  // MM1=[ (c3+b3)*9-a3-d3+8, ...
385      MMXir( psrad,      5,mm0);
386      MMXar( movq,       q,mm7);    // MM7=[ p3,p2,p1,p0 ]
387      MMXir( psrad,      5,mm1);
388      MMXrr( packssdw,   mm1,mm0);  // MM0=[ x3,x2,x1,x0 ]
389      MMXrr( paddw,      mm0,mm7);  // MM7=[ p3+x3, p2+x2, ... ]
390      MMXra( movq,       mm7,q);
391#if defined(_MSC_VER) && defined(_DEBUG)
392      MMXemms;
393#endif
394      q += 4;
395    }
396}
397
398#endif /* MMX */
399
400//////////////////////////////////////////////////////
401// NEW FILTERS
402//////////////////////////////////////////////////////
403
404static void 
405filter_fv(short *p, int w, int h, int rowsize, int scale)
406{
407  int y = 0;
408  int s = scale*rowsize;
409  int s3 = s+s+s;
410  h = ((h-1)/scale)+1;
411  y += 1;
412  p += s;
413  while (y-3 < h)
414    {
415      // 1-Delta
416      {
417        short *q = p;
418        short *e = q+w;
419        if (y>=3 && y+3<h)
420          {
421            // Generic case
422#ifdef MMX
423            if (scale==1 && MMXControl::mmxflag>0)
424              mmx_fv_1(q, e, s, s3);
425#endif
426            while (q<e)
427              {
428                int a = (int)q[-s] + (int)q[s];
429                int b = (int)q[-s3] + (int)q[s3];
430                *q -= (((a<<3)+a-b+8)>>4);
431                q += scale;
432              }
433          }
434        else if (y<h)
435          {
436            // Special cases
437            short *q1 = (y+1<h ? q+s : q-s);
438            while (q<e)
439              {
440                int a = (int)q[-s] + (int)(*q1);
441                *q -= ((a+1)>>1);
442                q += scale;
443                q1 += scale;
444              }
445          }
446      }
447      // 2-Update
448      {
449        short *q = p-s3;
450        short *e = q+w;
451        if (y>=6 && y<h)
452          {
453            // Generic case
454#ifdef MMX
455            if (scale==1 && MMXControl::mmxflag>0)
456              mmx_fv_2(q, e, s, s3);
457#endif
458            while (q<e)
459              {
460                int a = (int)q[-s] + (int)q[s];
461                int b = (int)q[-s3] + (int)q[s3];
462                *q += (((a<<3)+a-b+16)>>5);
463                q += scale;
464              }
465          }
466        else if (y>=3)
467          {
468            // Special cases
469            short *q1 = (y-2<h ? q+s : 0);
470            short *q3 = (y<h ? q+s3 : 0);
471            if (y>=6)
472              {
473                while (q<e)
474                  {
475                    int a = (int)q[-s] + (q1 ? (int)(*q1) : 0);
476                    int b = (int)q[-s3] + (q3 ? (int)(*q3) : 0);
477                    *q += (((a<<3)+a-b+16)>>5);
478                    q += scale;
479                    if (q1) q1 += scale;
480                    if (q3) q3 += scale;
481                  }
482              }
483            else if (y>=4)
484              {
485                while (q<e)
486                  {
487                    int a = (int)q[-s] + (q1 ? (int)(*q1) : 0);
488                    int b = (q3 ? (int)(*q3) : 0);
489                    *q += (((a<<3)+a-b+16)>>5);
490                    q += scale;
491                    if (q1) q1 += scale;
492                    if (q3) q3 += scale;
493                  }
494              }
495            else
496              {
497                while (q<e)
498                  {
499                    int a = (q1 ? (int)(*q1) : 0);
500                    int b = (q3 ? (int)(*q3) : 0);
501                    *q += (((a<<3)+a-b+16)>>5);
502                    q += scale;
503                    if (q1) q1 += scale;
504                    if (q3) q3 += scale;
505                  }
506              }
507          }
508      }
509      y += 2;
510      p += s+s;
511    }
512}
513
514static void 
515filter_fh(short *p, int w, int h, int rowsize, int scale)
516{
517  int y = 0;
518  int s = scale;
519  int s3 = s+s+s;
520  rowsize *= scale;
521  while (y<h)
522    {
523      short *q = p+s;
524      short *e = p+w;
525      int a0=0, a1=0, a2=0, a3=0;
526      int b0=0, b1=0, b2=0, b3=0;
527      if (q < e)
528        {
529          // Special case: x=1
530          a1 = a2 = a3 = q[-s];
531          if (q+s<e)
532            a2 = q[s];
533          if (q+s3<e)
534            a3 = q[s3];
535          b3 = q[0] - ((a1+a2+1)>>1);
536          q[0] = b3;
537          q += s+s;
538        }
539      while (q+s3 < e)
540        {
541          // Generic case
542          a0=a1; 
543          a1=a2; 
544          a2=a3;
545          a3=q[s3];
546          b0=b1; 
547          b1=b2; 
548          b2=b3;
549          b3 = q[0] - ((((a1+a2)<<3)+(a1+a2)-a0-a3+8) >> 4);
550          q[0] = b3;
551          q[-s3] = q[-s3] + ((((b1+b2)<<3)+(b1+b2)-b0-b3+16) >> 5);
552          q += s+s;
553        }
554      while (q < e)
555        {
556          // Special case: w-3 <= x < w
557          a1=a2; 
558          a2=a3;
559          b0=b1; 
560          b1=b2; 
561          b2=b3;
562          b3 = q[0] - ((a1+a2+1)>>1);
563          q[0] = b3;
564          q[-s3] = q[-s3] + ((((b1+b2)<<3)+(b1+b2)-b0-b3+16) >> 5);
565          q += s+s;
566        }
567      while (q-s3 < e)
568        {
569          // Special case  w <= x < w+3
570          b0=b1; 
571          b1=b2; 
572          b2=b3;
573          b3=0;
574          if (q-s3 >= p)
575            q[-s3] = q[-s3] + ((((b1+b2)<<3)+(b1+b2)-b0-b3+16) >> 5);
576          q += s+s;
577        }
578      y += scale;
579      p += rowsize;
580    }
581}
582
583
584//////////////////////////////////////////////////////
585// WAVELET TRANSFORM
586//////////////////////////////////////////////////////
587
588
589//----------------------------------------------------
590// Function for applying bidimensional IW44 between
591// scale intervals begin(inclusive) and end(exclusive)
592
593void
594IW44Image::Transform::Encode::forward(short *p, int w, int h, int rowsize, int begin, int end)
595{ 
596
597  // PREPARATION
598  filter_begin(w,h);
599  // LOOP ON SCALES
600  for (int scale=begin; scale<end; scale<<=1)
601    {
602#ifdef IWTRANSFORM_TIMER
603      int tv,th;
604      th = tv = GOS::ticks();
605#endif
606      filter_fh(p, w, h, rowsize, scale);
607#ifdef IWTRANSFORM_TIMER
608      th = GOS::ticks();
609      tv = th - tv;
610#endif
611      filter_fv(p, w, h, rowsize, scale);
612#ifdef IWTRANSFORM_TIMER
613      th = GOS::ticks()-th;
614      DjVuPrintErrorUTF8("forw%d\tv=%dms h=%dms\n", scale,th,tv);
615#endif
616    }
617  // TERMINATE
618  filter_end();
619}
620
621//////////////////////////////////////////////////////
622// COLOR TRANSFORM
623//////////////////////////////////////////////////////
624
625static const float 
626rgb_to_ycc[3][3] = 
627{ { 0.304348F,  0.608696F,  0.086956F },     
628  { 0.463768F, -0.405797F, -0.057971F },
629  {-0.173913F, -0.347826F,  0.521739F } };
630
631
632/* Extracts Y */
633void 
634IW44Image::Transform::Encode::RGB_to_Y(const GPixel *p, int w, int h, int rowsize, 
635                      signed char *out, int outrowsize)
636{
637  int rmul[256], gmul[256], bmul[256];
638  for (int k=0; k<256; k++)
639    {
640      rmul[k] = (int)(k*0x10000*rgb_to_ycc[0][0]);
641      gmul[k] = (int)(k*0x10000*rgb_to_ycc[0][1]);
642      bmul[k] = (int)(k*0x10000*rgb_to_ycc[0][2]);
643    }
644  for (int i=0; i<h; i++, p+=rowsize, out+=outrowsize)
645    {
646      const GPixel *p2 = p;
647      signed char *out2 = out;
648      for (int j=0; j<w; j++,p2++,out2++)
649        {
650          int y = rmul[p2->r] + gmul[p2->g] + bmul[p2->b] + 32768;
651          *out2 = (y>>16) - 128;
652        }
653    }
654}
655
656#ifdef min
657#undef min
658#endif
659static inline int min(const int x,const int y) {return (x<y)?x:y;}
660#ifdef max
661#undef max
662#endif
663static inline int max(const int x,const int y) {return (x>y)?x:y;}
664
665/* Extracts Cb */
666void 
667IW44Image::Transform::Encode::RGB_to_Cb(const GPixel *p, int w, int h, int rowsize, 
668                       signed char *out, int outrowsize)
669{
670  int rmul[256], gmul[256], bmul[256];
671  for (int k=0; k<256; k++)
672    {
673      rmul[k] = (int)(k*0x10000*rgb_to_ycc[2][0]);
674      gmul[k] = (int)(k*0x10000*rgb_to_ycc[2][1]);
675      bmul[k] = (int)(k*0x10000*rgb_to_ycc[2][2]);
676    }
677  for (int i=0; i<h; i++, p+=rowsize, out+=outrowsize)
678    {
679      const GPixel *p2 = p;
680      signed char *out2 = out;
681      for (int j=0; j<w; j++,p2++,out2++)
682        {
683          int c = rmul[p2->r] + gmul[p2->g] + bmul[p2->b] + 32768;
684          *out2 = max(-128, min(127, c>>16));
685        }
686    }
687}
688
689/* Extracts Cr */
690void 
691IW44Image::Transform::Encode::RGB_to_Cr(const GPixel *p, int w, int h, int rowsize, 
692                       signed char *out, int outrowsize)
693{
694  int rmul[256], gmul[256], bmul[256];
695  for (int k=0; k<256; k++)
696    {
697      rmul[k] = (int)((k*0x10000)*rgb_to_ycc[1][0]);
698      gmul[k] = (int)((k*0x10000)*rgb_to_ycc[1][1]);
699      bmul[k] = (int)((k*0x10000)*rgb_to_ycc[1][2]);
700    }
701  for (int i=0; i<h; i++, p+=rowsize, out+=outrowsize)
702    {
703      const GPixel *p2 = p;
704      signed char *out2 = out;
705      for (int j=0; j<w; j++,p2++,out2++)
706        {
707          int c = rmul[p2->r] + gmul[p2->g] + bmul[p2->b] + 32768;
708          *out2 = max(-128, min(127, c>>16));
709        }
710    }
711}
712
713
714//////////////////////////////////////////////////////
715// MASKING DECOMPOSITION
716//////////////////////////////////////////////////////
717
718//----------------------------------------------------
719// Function for applying bidimensional IW44 between
720// scale intervals begin(inclusive) and end(exclusive)
721// with a MASK bitmap
722
723
724static void
725interpolate_mask(short *data16, int w, int h, int rowsize,
726                 const signed char *mask8, int mskrowsize)
727{
728  int i,j;
729  // count masked bits
730  short *count;
731  GPBuffer<short> gcount(count,w*h);
732  short *cp = count;
733  for (i=0; i<h; i++, cp+=w, mask8+=mskrowsize)
734    for (j=0; j<w; j++)
735      cp[j] = (mask8[j] ? 0 : 0x1000);
736  // copy image
737  short *sdata;
738  GPBuffer<short> gsdata(sdata,w*h);
739  short *p = sdata;
740  short *q = data16;
741  for (i=0; i<h; i++, p+=w, q+=rowsize)
742    for (j=0; j<w; j++)
743      p[j] = q[j];
744  // iterate over resolutions
745  int split = 1;
746  int scale = 2;
747  int again = 1;
748  while (again && scale<w && scale<h)
749    {
750      again = 0;
751      p = data16;
752      q = sdata;
753      cp = count;
754      // iterate over block
755      for (i=0; i<h; i+=scale, cp+=w*scale, q+=w*scale, p+=rowsize*scale)
756        for (j=0; j<w; j+=scale)
757          {
758            int ii, jj;
759            int gotz = 0;
760            int gray = 0;
761            int npix = 0;
762            short *cpp = cp;
763            short *qq = q;
764            // look around when square goes beyond border
765            int istart = i;
766            if (istart+split>h)
767              {
768                istart -= scale;
769                cpp -= w*scale;
770                qq -= w*scale;
771              }
772            int jstart = j;
773            if (jstart+split>w)
774              jstart -= scale;
775            // compute gray level
776            for (ii=istart; ii<i+scale && ii<h; ii+=split, cpp+=w*split, qq+=w*split)
777              for (jj=jstart; jj<j+scale && jj<w; jj+=split)
778                {
779                  if (cpp[jj]>0) 
780                    {
781                      npix += cpp[jj];
782                      gray += cpp[jj] * qq[jj];
783                    } 
784                  else if (ii>=i && jj>=j)
785                    {
786                      gotz = 1;
787                    }
788                }
789            // process result
790            if (npix == 0)
791              {
792                // continue to next resolution
793                again = 1;
794                cp[j] = 0;
795              }
796            else
797              {
798                gray = gray / npix;
799                // check whether initial image require fix
800                if (gotz)
801                  {
802                    cpp = cp;
803                    qq = p;
804                    for (ii=i; ii<i+scale && ii<h; ii+=1, cpp+=w, qq+=rowsize)
805                      for (jj=j; jj<j+scale && jj<w; jj+=1)
806                        if (cpp[jj] == 0)
807                          {
808                            qq[jj] = gray;
809                            cpp[jj] = 1;
810                          }
811                  }
812                // store average for next iteration
813                cp[j] = npix>>2;
814                q[j] = gray;
815              }
816          }
817      // double resolution
818      split = scale;
819      scale = scale+scale;
820    }
821}
822
823
824static void
825forward_mask(short *data16, int w, int h, int rowsize, int begin, int end,
826             const signed char *mask8, int mskrowsize )
827{
828  int i,j;
829  signed char *m;
830  short *p;
831  short *d;
832  // Allocate buffers
833  short *sdata;
834  GPBuffer<short> gsdata(sdata,w*h);
835  signed char *smask;
836  GPBuffer<signed char> gsmask(smask,w*h);
837  // Copy mask
838  m = smask;
839  for (i=0; i<h; i+=1, m+=w, mask8+=mskrowsize)
840    memcpy((void*)m, (void*)mask8, w);
841  // Loop over scale
842  for (int scale=begin; scale<end; scale<<=1)
843    {
844      // Copy data into sdata buffer
845      p = data16;
846      d = sdata;
847      for (i=0; i<h; i+=scale)
848        {
849          for (j=0; j<w; j+=scale)
850            d[j] = p[j];
851          p += rowsize * scale;
852          d += w * scale;
853        }
854      // Decompose
855      IW44Image::Transform::Encode::forward(sdata, w, h, w, scale, scale+scale);
856      // Cancel masked coefficients
857      d = sdata;
858      m = smask;
859      for (i=0; i<h; i+=scale+scale)
860        {
861          for (j=scale; j<w; j+=scale+scale)
862            if (m[j])
863              d[j] = 0;
864          d += w * scale;
865          m += w * scale;
866          if (i+scale < h)
867            {
868              for (j=0; j<w; j+=scale)
869                if (m[j])
870                  d[j] = 0;
871              d += w * scale;
872              m += w * scale;
873            }
874        }
875      // Reconstruct
876      IW44Image::Transform::Decode::backward(sdata, w, h, w, scale+scale, scale);
877      // Correct visible pixels
878      p = data16;
879      d = sdata;
880      m = smask;
881      for (i=0; i<h; i+=scale)
882        {
883          for (j=0; j<w; j+=scale)
884            if (! m[j])
885              d[j] = p[j];
886          p += rowsize*scale;
887          m += w*scale;
888          d += w*scale;
889        }
890      // Decompose again (no need to iterate actually!)
891      IW44Image::Transform::Encode::forward(sdata, w, h, w, scale, scale+scale);
892      // Copy coefficients from sdata buffer
893      p = data16;
894      d = sdata;
895      for (i=0; i<h; i+=scale)
896        {
897          for (j=0; j<w; j+=scale)
898            p[j] = d[j];
899          p += rowsize * scale;
900          d += w * scale;
901        }
902      // Compute new mask for next scale
903      m = smask;
904      signed char *m0 = m;
905      signed char *m1 = m;
906      for (i=0; i<h; i+=scale+scale)
907        {
908          m0 = m1;
909          if (i+scale < h)
910            m1 = m + w*scale;
911          for (j=0; j<w; j+=scale+scale)
912            if (m[j] && m0[j] && m1[j] && (j<=0 || m[j-scale]) && (j+scale>=w || m[j+scale]))
913              m[j] = 1;
914            else
915              m[j] = 0;
916          m = m1 + w*scale;
917        }
918    }
919  // Free buffers
920}
921
922void 
923IW44Image::Map::Encode::create(const signed char *img8, int imgrowsize, 
924               const signed char *msk8, int mskrowsize )
925{
926  int i, j;
927  // Progress
928  DJVU_PROGRESS_TASK(transf,"create iw44 map",3);
929  // Allocate decomposition buffer
930  short *data16;
931  GPBuffer<short> gdata16(data16,bw*bh);
932  // Copy pixels
933  short *p = data16;
934  const signed char *row = img8;
935  for (i=0; i<ih; i++)
936    {
937      for (j=0; j<iw; j++)
938        *p++ = (int)(row[j]) << iw_shift;
939      row += imgrowsize;
940      for (j=iw; j<bw; j++)
941        *p++ = 0;
942    }
943  for (i=ih; i<bh; i++)
944    for (j=0; j<bw; j++)
945      *p++ = 0;
946  // Handle bitmask
947  if (msk8)
948    {
949      // Interpolate pixels below mask
950      DJVU_PROGRESS_RUN(transf, 1);
951      interpolate_mask(data16, iw, ih, bw, msk8, mskrowsize);
952      // Multiscale iterative masked decomposition
953      DJVU_PROGRESS_RUN(transf, 3);
954      forward_mask(data16, iw, ih, bw, 1, 32, msk8, mskrowsize);
955    }
956  else
957    {
958      // Perform traditional decomposition
959      DJVU_PROGRESS_RUN(transf, 3);
960      IW44Image::Transform::Encode::forward(data16, iw, ih, bw, 1, 32);
961    }
962  // Copy coefficient into blocks
963  p = data16;
964  IW44Image::Block *block = blocks;
965  for (i=0; i<bh; i+=32)
966    {
967      for (j=0; j<bw; j+=32)
968        {
969          short liftblock[1024];
970          // transfer coefficients at (p+j) into aligned block
971          short *pp = p + j;
972          short *pl = liftblock;
973          for (int ii=0; ii<32; ii++, pp+=bw)
974            for (int jj=0; jj<32; jj++) 
975              *pl++ = pp[jj];
976          // transfer into IW44Image::Block (apply zigzag and scaling)
977          block->read_liftblock(liftblock, this);
978          block++;
979        }
980      // next row of blocks
981      p += 32*bw;
982    }
983}
984
985void 
986IW44Image::Map::Encode::slashres(int res)
987{
988  int minbucket = 1;
989  if (res < 2)
990    return;
991  else if (res < 4)
992    minbucket=16;
993  else if (res < 8)
994    minbucket=4;
995  for (int blockno=0; blockno<nb; blockno++)
996    for (int buckno=minbucket; buckno<64; buckno++)
997      blocks[blockno].zero(buckno);
998}
999
1000// encode_prepare
1001// -- compute the states prior to encoding the buckets
1002int
1003IW44Image::Codec::Encode::encode_prepare(int band, int fbucket, int nbucket, IW44Image::Block &blk, IW44Image::Block &eblk)
1004{
1005  int bbstate = 0;
1006  // compute state of all coefficients in all buckets
1007  if (band) 
1008    {
1009      // Band other than zero
1010      int thres = quant_hi[band];
1011      char *cstate = coeffstate;
1012      for (int buckno=0; buckno<nbucket; buckno++, cstate+=16)
1013        {
1014          const short *pcoeff = blk.data(fbucket+buckno);
1015          const short *epcoeff = eblk.data(fbucket+buckno);
1016          int bstatetmp = 0;
1017          if (! pcoeff)
1018            {
1019              bstatetmp = UNK;
1020              // cstate[i] is not used and does not need initialization
1021            }
1022          else if (! epcoeff)
1023            {
1024              for (int i=0; i<16; i++)
1025                {
1026                  int cstatetmp = UNK;
1027                  if  ((int)(pcoeff[i])>=thres || (int)(pcoeff[i])<=-thres)
1028                    cstatetmp = NEW|UNK;
1029                  cstate[i] = cstatetmp;
1030                  bstatetmp |= cstatetmp;
1031                }
1032            }
1033          else
1034            {
1035              for (int i=0; i<16; i++)
1036                {
1037                  int cstatetmp = UNK;
1038                  if (epcoeff[i])
1039                    cstatetmp = ACTIVE;
1040                  else if  ((int)(pcoeff[i])>=thres || (int)(pcoeff[i])<=-thres)
1041                    cstatetmp = NEW|UNK;
1042                  cstate[i] = cstatetmp;
1043                  bstatetmp |= cstatetmp;
1044                }
1045            }
1046          bucketstate[buckno] = bstatetmp;
1047          bbstate |= bstatetmp;
1048        }
1049    }
1050  else
1051    {
1052      // Band zero ( fbucket==0 implies band==zero and nbucket==1 )
1053      const short *pcoeff = blk.data(0, &map);
1054      const short *epcoeff = eblk.data(0, &emap);
1055      char *cstate = coeffstate;
1056      for (int i=0; i<16; i++)
1057        {
1058          int thres = quant_lo[i];
1059          int cstatetmp = cstate[i];
1060          if (cstatetmp != ZERO)
1061            {
1062              cstatetmp = UNK;
1063              if (epcoeff[i])
1064                cstatetmp = ACTIVE;
1065              else if ((int)(pcoeff[i])>=thres || (int)(pcoeff[i])<=-thres)
1066                cstatetmp = NEW|UNK;
1067            }
1068          cstate[i] = cstatetmp;
1069          bbstate |= cstatetmp;
1070        }
1071      bucketstate[0] = bbstate;
1072    }
1073  return bbstate;
1074}
1075
1076// encode_buckets
1077// -- code a sequence of buckets in a given block
1078void
1079IW44Image::Codec::Encode::encode_buckets(ZPCodec &zp, int bit, int band, 
1080                         IW44Image::Block &blk, IW44Image::Block &eblk,
1081                         int fbucket, int nbucket)
1082{
1083  // compute state of all coefficients in all buckets
1084  int bbstate = encode_prepare(band, fbucket, nbucket, blk, eblk);
1085
1086  // code root bit
1087  if ((nbucket<16) || (bbstate&ACTIVE))
1088    {
1089      bbstate |= NEW;
1090    }
1091  else if (bbstate & UNK)
1092    {
1093      zp.encoder( (bbstate&NEW) ? 1 : 0 , ctxRoot);
1094#ifdef TRACE
1095      DjVuPrintMessage("bbstate[bit=%d,band=%d] = %d\n", bit, band, bbstate);
1096#endif
1097    }
1098 
1099  // code bucket bits
1100  if (bbstate & NEW)
1101    for (int buckno=0; buckno<nbucket; buckno++)
1102      {
1103        // Code bucket bit
1104        if (bucketstate[buckno] & UNK)
1105          {
1106            // Context
1107            int ctx = 0;
1108#ifndef NOCTX_BUCKET_UPPER
1109            if (band>0)
1110              {
1111                int k = (fbucket+buckno)<<2;
1112                const short *b = eblk.data(k>>4);
1113                if (b)
1114                  {
1115                    k = k & 0xf;
1116                    if (b[k])
1117                      ctx += 1;
1118                    if (b[k+1])
1119                      ctx += 1;
1120                    if (b[k+2])
1121                      ctx += 1;
1122                    if (ctx<3 && b[k+3])
1123                      ctx += 1;
1124                  }
1125              }
1126#endif
1127#ifndef NOCTX_BUCKET_ACTIVE
1128            if (bbstate & ACTIVE)
1129              ctx |= 4; 
1130#endif
1131            // Code
1132            zp.encoder( (bucketstate[buckno]&NEW) ? 1 : 0, ctxBucket[band][ctx] );
1133#ifdef TRACE
1134            DjVuPrintMessage("  bucketstate[bit=%d,band=%d,buck=%d] = %d\n", 
1135                   bit, band, buckno, bucketstate[buckno] & ~ZERO);
1136#endif
1137          }
1138      }
1139 
1140  // code new active coefficient (with their sign)
1141  if (bbstate & NEW)
1142    {
1143      int thres = quant_hi[band];
1144      char *cstate = coeffstate;
1145      for (int buckno=0; buckno<nbucket; buckno++, cstate+=16)
1146        if (bucketstate[buckno] & NEW)
1147          {
1148            int i;
1149#ifndef NOCTX_EXPECT
1150            int gotcha = 0;
1151            const int maxgotcha = 7;
1152            for (i=0; i<16; i++)
1153              if (cstate[i] & UNK)
1154                gotcha += 1;
1155#endif
1156            const short *pcoeff = blk.data(fbucket+buckno);
1157            short *epcoeff = eblk.data(fbucket+buckno, &emap);
1158            // iterate within bucket
1159            for (i=0; i<16; i++)
1160              {
1161                if (cstate[i] & UNK)
1162                  {
1163                    // Prepare context
1164                    int ctx = 0;
1165#ifndef NOCTX_EXPECT
1166                    if (gotcha>=maxgotcha)
1167                      ctx = maxgotcha;
1168                    else
1169                      ctx = gotcha;
1170#endif
1171#ifndef NOCTX_ACTIVE
1172                    if (bucketstate[buckno] & ACTIVE)
1173                      ctx |= 8;
1174#endif
1175                    // Code
1176                    zp.encoder( (cstate[i]&NEW) ? 1 : 0, ctxStart[ctx] );
1177                    if (cstate[i] & NEW)
1178                      {
1179                        // Code sign
1180                        zp.IWencoder( (pcoeff[i]<0) ? 1 : 0 );
1181                        // Set encoder state
1182                        if (band==0)
1183                          thres = quant_lo[i];
1184                        epcoeff[i] = thres + (thres>>1);
1185                      }
1186#ifndef NOCTX_EXPECT
1187                    if (cstate[i] & NEW)
1188                      gotcha = 0;
1189                    else if (gotcha > 0)
1190                      gotcha -= 1;
1191#endif
1192#ifdef TRACE
1193                    DjVuPrintMessage("    coeffstate[bit=%d,band=%d,buck=%d,c=%d] = %d\n", 
1194                           bit, band, buckno, i, cstate[i]);
1195#endif
1196                  }
1197              }
1198          }
1199    }
1200
1201  // code mantissa bits
1202  if (bbstate & ACTIVE)
1203    {
1204      int thres = quant_hi[band];
1205      char *cstate = coeffstate;
1206      for (int buckno=0; buckno<nbucket; buckno++, cstate+=16)
1207        if (bucketstate[buckno] & ACTIVE)
1208          {
1209            const short *pcoeff = blk.data(fbucket+buckno);
1210            short *epcoeff = eblk.data(fbucket+buckno, &emap);
1211            for (int i=0; i<16; i++)
1212              if (cstate[i] & ACTIVE)
1213                {
1214                  // get coefficient
1215                  int coeff = pcoeff[i];
1216                  int ecoeff = epcoeff[i];
1217                  if (coeff < 0)
1218                    coeff = -coeff;
1219                  // get band zero thresholds
1220                  if (band == 0)
1221                    thres = quant_lo[i];
1222                  // compute mantissa bit
1223                  int pix = 0;
1224                  if (coeff >= ecoeff)
1225                    pix = 1;
1226                  // encode second or lesser mantissa bit
1227                  if (ecoeff <= 3*thres)
1228                    zp.encoder(pix, ctxMant);                     
1229                  else
1230                                          zp.IWencoder(!!pix);
1231                  // adjust epcoeff
1232                  epcoeff[i] = ecoeff - (pix ? 0 : thres) + (thres>>1);
1233                }
1234          }
1235    }
1236}
1237
1238// IW44Image::Codec::estimate_decibel
1239// -- estimate encoding error (after code_slice) in decibels.
1240float
1241IW44Image::Codec::Encode::estimate_decibel(float frac)
1242{
1243  int i,j;
1244  const float *q;
1245  // Fill norm arrays
1246  float norm_lo[16];
1247  float norm_hi[10];
1248  // -- lo coefficients
1249  q = iw_norm;
1250  for (i=j=0; i<4; j++)
1251    norm_lo[i++] = *q++;
1252  for (j=0; j<4; j++)
1253    norm_lo[i++] = *q;
1254  q += 1;
1255  for (j=0; j<4; j++)
1256    norm_lo[i++] = *q;
1257  q += 1;
1258  for (j=0; j<4; j++)
1259    norm_lo[i++] = *q;
1260  q += 1;
1261  // -- hi coefficients
1262  norm_hi[0] = 0;
1263  for (j=1; j<10; j++)
1264    norm_hi[j] = *q++;
1265  // Initialize mse array
1266  float *xmse;
1267  GPBuffer<float> gxmse(xmse,map.nb);
1268  // Compute mse in each block
1269  for (int blockno=0; blockno<map.nb; blockno++)
1270    {
1271      float mse = 0;
1272      // Iterate over bands
1273      for (int bandno=0; bandno<10; bandno++)
1274        {
1275          int fbucket = bandbuckets[bandno].start;
1276          int nbucket = bandbuckets[bandno].size;
1277          IW44Image::Block &blk = map.blocks[blockno];
1278          IW44Image::Block &eblk = emap.blocks[blockno];
1279          float norm = norm_hi[bandno];
1280          for (int buckno=0; buckno<nbucket; buckno++)
1281            {
1282              const short *pcoeff = blk.data(fbucket+buckno);
1283              const short *epcoeff = eblk.data(fbucket+buckno);
1284              if (pcoeff)
1285                {
1286                  if (epcoeff)
1287                    {
1288                      for (i=0; i<16; i++)
1289                        {
1290                          if (bandno == 0)
1291                            norm = norm_lo[i];
1292                          float delta = (float)(pcoeff[i]<0 ? -pcoeff[i] : pcoeff[i]);
1293                          delta = delta - epcoeff[i];
1294                          mse = mse + norm * delta * delta;
1295                        }
1296                    }
1297                  else
1298                    {
1299                      for (i=0; i<16; i++)
1300                        {
1301                          if (bandno == 0)
1302                            norm = norm_lo[i];
1303                          float delta = (float)(pcoeff[i]);
1304                          mse = mse + norm * delta * delta;
1305                        }
1306                    }
1307                }
1308            }
1309        }
1310      xmse[blockno] = mse / 1024;
1311    }
1312  // Compute partition point
1313  int n = 0;
1314  int m = map.nb - 1;
1315  int p = (int)floor(m*(1.0-frac)+0.5);
1316  p = (p>m ? m : (p<0 ? 0 : p));
1317  float pivot = 0;
1318  // Partition array
1319  while (n < p)
1320    {
1321      int l = n;
1322      int h = m;
1323      if (xmse[l] > xmse[h]) { float tmp=xmse[l]; xmse[l]=xmse[h]; xmse[h]=tmp; }
1324      pivot = xmse[(l+h)/2];
1325      if (pivot < xmse[l]) { float tmp=pivot; pivot=xmse[l]; xmse[l]=tmp; }
1326      if (pivot > xmse[h]) { float tmp=pivot; pivot=xmse[h]; xmse[h]=tmp; }
1327      while (l < h)
1328        {
1329          if (xmse[l] > xmse[h]) { float tmp=xmse[l]; xmse[l]=xmse[h]; xmse[h]=tmp; }
1330          while (xmse[l]<pivot || (xmse[l]==pivot && l<h)) l++;
1331          while (xmse[h]>pivot) h--;
1332        }
1333      if (p>=l) 
1334        n = l;
1335      else 
1336        m = l-1;
1337    }
1338  // Compute average mse
1339  float mse = 0;
1340  for (i=p; i<map.nb; i++)
1341    mse = mse + xmse[i];
1342  mse = mse / (map.nb - p);
1343  // Return
1344  float factor = 255 << iw_shift;
1345  float decibel = (float)(10.0 * log ( factor * factor / mse ) / 2.302585125);
1346  return decibel;
1347}
1348
1349
1350
1351
1352//////////////////////////////////////////////////////
1353// IW44IMAGE ENCODING ROUTINES
1354//////////////////////////////////////////////////////
1355
1356
1357void 
1358IW44Image::PrimaryHeader::encode(GP<ByteStream> gbs)
1359{
1360  gbs->write8(serial);
1361  gbs->write8(slices);
1362}
1363
1364void 
1365IW44Image::SecondaryHeader::encode(GP<ByteStream> gbs)
1366{
1367  gbs->write8(major);
1368  gbs->write8(minor);
1369}
1370
1371void 
1372IW44Image::TertiaryHeader::encode(GP<ByteStream> gbs)
1373{
1374  gbs->write8(xhi);
1375  gbs->write8(xlo);
1376  gbs->write8(yhi);
1377  gbs->write8(ylo);
1378  gbs->write8(crcbdelay);
1379}
1380
1381
1382
1383GP<IW44Image>
1384IW44Image::create_encode(const ImageType itype)
1385{
1386  switch(itype)
1387  {
1388  case COLOR:
1389    return new IWPixmap::Encode();
1390  case GRAY:
1391    return new IWBitmap::Encode();
1392  default:
1393    return 0;
1394  }
1395}
1396
1397GP<IW44Image>
1398IW44Image::create_encode(const GBitmap &bm, const GP<GBitmap> mask)
1399{
1400  IWBitmap::Encode *bit=new IWBitmap::Encode();
1401  GP<IW44Image> retval=bit;
1402  bit->init(bm, mask);
1403  return retval;
1404}
1405
1406
1407IWBitmap::Encode::Encode(void)
1408: IWBitmap(), ycodec_enc(0)
1409{}
1410
1411IWBitmap::Encode::~Encode()
1412{
1413  close_codec();
1414}
1415
1416void
1417IWBitmap::Encode::init(const GBitmap &bm, const GP<GBitmap> gmask)
1418{
1419  // Free
1420  close_codec();
1421  delete ymap;
1422  ymap = 0;
1423  // Init
1424  int i, j;
1425  int w = bm.columns();
1426  int h = bm.rows();
1427  int g = bm.get_grays()-1;
1428  signed char *buffer;
1429  GPBuffer<signed char> gbuffer(buffer,w*h);
1430  // Prepare gray level conversion table
1431  signed char  bconv[256];
1432  for (i=0; i<256; i++)
1433    bconv[i] = max(0,min(255,i*255/g)) - 128;
1434  // Perform decomposition
1435  // Prepare mask information
1436  const signed char *msk8 = 0;
1437  int mskrowsize = 0;
1438  GBitmap *mask=gmask;
1439  if (gmask)
1440  {
1441    msk8 = (const signed char*)((*mask)[0]);
1442    mskrowsize = mask->rowsize();
1443  }
1444  // Prepare a buffer of signed bytes
1445  for (i=0; i<h; i++)
1446    {
1447      signed char *bufrow = buffer + i*w;
1448      const unsigned char *bmrow = bm[i];
1449      for (j=0; j<w; j++)
1450        bufrow[j] = bconv[bmrow[j]];
1451    }
1452  // Create map
1453  Map::Encode *eymap=new Map::Encode(w,h);
1454  ymap = eymap;
1455  eymap->create(buffer, w, msk8, mskrowsize);
1456}
1457
1458void 
1459IWBitmap::Encode::close_codec(void)
1460{
1461  delete ycodec_enc;
1462  ycodec_enc = 0;
1463  IWBitmap::close_codec();
1464}
1465
1466int 
1467IWBitmap::Encode::encode_chunk(GP<ByteStream> gbs, const IWEncoderParms &parm)
1468{
1469  // Check
1470  if (parm.slices==0 && parm.bytes==0 && parm.decibels==0)
1471    G_THROW( ERR_MSG("IW44Image.need_stop") );
1472  if (! ymap)
1473    G_THROW( ERR_MSG("IW44Image.empty_object") );
1474  // Open codec
1475  if (!ycodec_enc)
1476  {
1477    cslice = cserial = cbytes = 0;
1478    ycodec_enc = new Codec::Encode(*ymap);
1479  }
1480  // Adjust cbytes
1481  cbytes += sizeof(struct IW44Image::PrimaryHeader);
1482  if (cserial == 0)
1483    cbytes += sizeof(struct IW44Image::SecondaryHeader) + sizeof(struct IW44Image::TertiaryHeader);
1484  // Prepare zcoded slices
1485  int flag = 1;
1486  int nslices = 0;
1487  GP<ByteStream> gmbs=ByteStream::create();
1488  ByteStream &mbs=*gmbs;
1489  DJVU_PROGRESS_TASK(chunk,"encode chunk",parm.slices-cslice);
1490  {
1491    float estdb = -1.0;
1492    GP<ZPCodec> gzp=ZPCodec::create(gmbs, true, true);
1493    ZPCodec &zp=*gzp;
1494    while (flag)
1495      {
1496        if (parm.decibels>0  && estdb>=parm.decibels)
1497          break;
1498        if (parm.bytes>0  && mbs.tell()+cbytes>=parm.bytes)
1499          break;
1500        if (parm.slices>0 && nslices+cslice>=parm.slices)
1501          break;
1502        DJVU_PROGRESS_RUN(chunk, (1+nslices-cslice)|0xf);
1503        flag = ycodec_enc->code_slice(zp);
1504        if (flag && parm.decibels>0.0)
1505          if (ycodec_enc->curband==0 || estdb>=parm.decibels-DECIBEL_PRUNE)
1506            estdb = ycodec_enc->estimate_decibel(db_frac);
1507        nslices++;
1508      }
1509  }
1510  // Write primary header
1511  struct IW44Image::PrimaryHeader primary;
1512  primary.serial = cserial;
1513  primary.slices = nslices;
1514  primary.encode(gbs);
1515  // Write auxilliary headers
1516  if (cserial == 0)
1517    {
1518      struct IW44Image::SecondaryHeader secondary;
1519      secondary.major = IWCODEC_MAJOR + 0x80;
1520      secondary.minor = IWCODEC_MINOR;
1521      secondary.encode(gbs);
1522      struct IW44Image::TertiaryHeader tertiary;
1523      tertiary.xhi = (ymap->iw >> 8) & 0xff;
1524      tertiary.xlo = (ymap->iw >> 0) & 0xff;
1525      tertiary.yhi = (ymap->ih >> 8) & 0xff;
1526      tertiary.ylo = (ymap->ih >> 0) & 0xff;
1527      tertiary.crcbdelay = 0;
1528      tertiary.encode(gbs);
1529    }
1530  // Write slices
1531  mbs.seek(0);
1532  gbs->copy(mbs);
1533  // Return
1534  cbytes  += mbs.tell();
1535  cslice  += nslices;
1536  cserial += 1;
1537  return flag;
1538}
1539
1540void 
1541IWBitmap::Encode::encode_iff(IFFByteStream &iff, int nchunks, const IWEncoderParms *parms)
1542{
1543  if (ycodec_enc)
1544    G_THROW( ERR_MSG("IW44Image.left_open1") );
1545  int flag = 1;
1546  iff.put_chunk("FORM:BM44", 1);
1547  DJVU_PROGRESS_TASK(iff,"encode iff chunk",nchunks);
1548  for (int i=0; flag && i<nchunks; i++)
1549    {
1550      DJVU_PROGRESS_RUN(iff,i+1);
1551      iff.put_chunk("BM44");
1552      flag = encode_chunk(iff.get_bytestream(),parms[i]);
1553      iff.close_chunk();
1554    }
1555  iff.close_chunk();
1556  close_codec();
1557}
1558
1559GP<IW44Image>
1560IW44Image::create_encode(
1561  const GPixmap &pm, const GP<GBitmap> gmask, CRCBMode crcbmode)
1562{
1563  IWPixmap::Encode *pix=new IWPixmap::Encode();
1564  GP<IW44Image> retval=pix;
1565  pix->init(pm, gmask,(IWPixmap::Encode::CRCBMode)crcbmode);
1566  return retval;
1567}
1568
1569IWPixmap::Encode::Encode(void)
1570: IWPixmap(), ycodec_enc(0), cbcodec_enc(0), crcodec_enc(0)
1571{}
1572
1573IWPixmap::Encode::~Encode()
1574{
1575  close_codec();
1576}
1577
1578void
1579IWPixmap::Encode::init(const GPixmap &pm, const GP<GBitmap> gmask, CRCBMode crcbmode)
1580{
1581  /* Free */
1582  close_codec();
1583  delete ymap;
1584  delete cbmap;
1585  delete crmap;
1586  ymap = cbmap = crmap = 0;
1587  /* Create */
1588  int w = pm.columns();
1589  int h = pm.rows();
1590  signed char *buffer;
1591  GPBuffer<signed char> gbuffer(buffer,w*h);
1592  // Create maps
1593  Map::Encode *eymap = new Map::Encode(w,h);
1594  ymap = eymap;
1595  // Handle CRCB mode
1596  switch (crcbmode) 
1597    {
1598    case CRCBnone:   crcb_half=1; crcb_delay=-1; break;
1599    case CRCBhalf:   crcb_half=1; crcb_delay=10; break;       
1600    case CRCBnormal: crcb_half=0; crcb_delay=10; break;
1601    case CRCBfull:   crcb_half=0; crcb_delay= 0; break;
1602    }
1603  // Prepare mask information
1604  const signed char *msk8 = 0;
1605  int mskrowsize = 0;
1606  GBitmap *mask=gmask;
1607  if (mask)
1608  {
1609    msk8 = (signed char const *)((*mask)[0]);
1610    mskrowsize = mask->rowsize();
1611  }
1612  // Fill buffer with luminance information
1613  DJVU_PROGRESS_TASK(create,"initialize pixmap",3);
1614  DJVU_PROGRESS_RUN(create,(crcb_delay>=0 ? 1 : 3));
1615  Transform::Encode::RGB_to_Y(pm[0], w, h, pm.rowsize(), buffer, w);
1616  if (crcb_delay < 0)
1617    {
1618      // Stupid inversion for gray images
1619      signed char *e = buffer + w*h;
1620      for (signed char *b=buffer; b<e; b++)
1621        *b = 255 - *b;
1622    }
1623  // Create YMAP
1624  eymap->create(buffer, w, msk8, mskrowsize);
1625  // Create chrominance maps
1626  if (crcb_delay >= 0)
1627    {
1628      Map::Encode *ecbmap = new Map::Encode(w,h);
1629      cbmap = ecbmap;
1630      Map::Encode *ecrmap = new Map::Encode(w,h);
1631      crmap = ecrmap;
1632      // Process CB information
1633      DJVU_PROGRESS_RUN(create,2);
1634      Transform::Encode::RGB_to_Cb(pm[0], w, h, pm.rowsize(), buffer, w);
1635      ecbmap->create(buffer, w, msk8, mskrowsize);
1636      // Process CR information
1637      DJVU_PROGRESS_RUN(create,3);
1638      Transform::Encode::RGB_to_Cr(pm[0], w, h, pm.rowsize(), buffer, w); 
1639      ecrmap->create(buffer, w, msk8, mskrowsize);
1640      // Perform chrominance reduction (CRCBhalf)
1641      if (crcb_half)
1642        {
1643          ecbmap->slashres(2);
1644          ecrmap->slashres(2);
1645        }
1646    }
1647}
1648
1649void 
1650IWPixmap::Encode::encode_iff(IFFByteStream &iff, int nchunks, const IWEncoderParms *parms)
1651{
1652  if (ycodec_enc)
1653    G_THROW( ERR_MSG("IW44Image.left_open3") );
1654  int flag = 1;
1655  iff.put_chunk("FORM:PM44", 1);
1656  DJVU_PROGRESS_TASK(iff,"encode pixmap chunk", nchunks);
1657  for (int i=0; flag && i<nchunks; i++)
1658    {
1659      DJVU_PROGRESS_RUN(iff,i+1);
1660      iff.put_chunk("PM44");
1661      flag = encode_chunk(iff.get_bytestream(), parms[i]);
1662      iff.close_chunk();
1663    }
1664  iff.close_chunk();
1665  close_codec();
1666}
1667
1668void 
1669IWPixmap::Encode::close_codec(void)
1670{
1671  delete ycodec_enc;
1672  delete cbcodec_enc;
1673  delete crcodec_enc;
1674  ycodec_enc = crcodec_enc = cbcodec_enc = 0;
1675  IWPixmap::close_codec();
1676}
1677
1678int 
1679IWPixmap::Encode::encode_chunk(GP<ByteStream> gbs, const IWEncoderParms &parm)
1680{
1681  // Check
1682  if (parm.slices==0 && parm.bytes==0 && parm.decibels==0)
1683    G_THROW( ERR_MSG("IW44Image.need_stop2") );
1684  if (!ymap)
1685    G_THROW( ERR_MSG("IW44Image.empty_object2") );
1686  // Open
1687  if (!ycodec_enc)
1688  {
1689    cslice = cserial = cbytes = 0;
1690    ycodec_enc = new Codec::Encode(*ymap);
1691    if (crmap && cbmap)
1692    {
1693      cbcodec_enc = new Codec::Encode(*cbmap);
1694      crcodec_enc = new Codec::Encode(*crmap);
1695    }
1696  }
1697
1698  // Adjust cbytes
1699  cbytes += sizeof(struct IW44Image::PrimaryHeader);
1700  if (cserial == 0)
1701    cbytes += sizeof(struct IW44Image::SecondaryHeader) + sizeof(struct IW44Image::TertiaryHeader);
1702  // Prepare zcodec slices
1703  int flag = 1;
1704  int nslices = 0;
1705  GP<ByteStream> gmbs=ByteStream::create();
1706  ByteStream &mbs=*gmbs;
1707  DJVU_PROGRESS_TASK(chunk, "encode pixmap chunk", parm.slices-cslice);
1708  {
1709    float estdb = -1.0;
1710    GP<ZPCodec> gzp=ZPCodec::create(gmbs, true, true);
1711    ZPCodec &zp=*gzp;
1712    while (flag)
1713      {
1714        if (parm.decibels>0  && estdb>=parm.decibels)
1715          break;
1716        if (parm.bytes>0  && mbs.tell()+cbytes>=parm.bytes)
1717          break;
1718        if (parm.slices>0 && nslices+cslice>=parm.slices)
1719          break;
1720        DJVU_PROGRESS_RUN(chunk,(1+nslices-cslice)|0xf);
1721        flag = ycodec_enc->code_slice(zp);
1722        if (flag && parm.decibels>0)
1723          if (ycodec_enc->curband==0 || estdb>=parm.decibels-DECIBEL_PRUNE)
1724            estdb = ycodec_enc->estimate_decibel(db_frac);
1725        if (crcodec_enc && cbcodec_enc && cslice+nslices>=crcb_delay)
1726          {
1727            flag |= cbcodec_enc->code_slice(zp);
1728            flag |= crcodec_enc->code_slice(zp);
1729          }
1730        nslices++;
1731      }
1732  }
1733  // Write primary header
1734  struct IW44Image::PrimaryHeader primary;
1735  primary.serial = cserial;
1736  primary.slices = nslices;
1737  primary.encode(gbs);
1738  // Write secondary header
1739  if (cserial == 0)
1740    {
1741      struct IW44Image::SecondaryHeader secondary;
1742      secondary.major = IWCODEC_MAJOR;
1743      secondary.minor = IWCODEC_MINOR;
1744      if (! (crmap && cbmap))
1745        secondary.major |= 0x80;
1746      secondary.encode(gbs);
1747      struct IW44Image::TertiaryHeader tertiary;
1748      tertiary.xhi = (ymap->iw >> 8) & 0xff;
1749      tertiary.xlo = (ymap->iw >> 0) & 0xff;
1750      tertiary.yhi = (ymap->ih >> 8) & 0xff;
1751      tertiary.ylo = (ymap->ih >> 0) & 0xff;
1752      tertiary.crcbdelay = (crcb_half ? 0x00 : 0x80);
1753      tertiary.crcbdelay |= (crcb_delay>=0 ? crcb_delay : 0x00);
1754      tertiary.encode(gbs);
1755    }
1756  // Write slices
1757  mbs.seek(0);
1758  gbs->copy(mbs);
1759  // Return
1760  cbytes  += mbs.tell();
1761  cslice  += nslices;
1762  cserial += 1;
1763  return flag;
1764}
1765
1766// code_slice
1767// -- read/write a slice of datafile
1768
1769int
1770IW44Image::Codec::Encode::code_slice(ZPCodec &zp)
1771{
1772  // Check that code_slice can still run
1773  if (curbit < 0)
1774    return 0;
1775  // Perform coding
1776  if (! is_null_slice(curbit, curband))
1777    {
1778      for (int blockno=0; blockno<map.nb; blockno++)
1779        {
1780          const int fbucket = bandbuckets[curband].start;
1781          const int nbucket = bandbuckets[curband].size;
1782          encode_buckets(zp, curbit, curband, 
1783                         map.blocks[blockno], emap.blocks[blockno], 
1784                         fbucket, nbucket);
1785        }
1786    }
1787  return finish_code_slice(zp);
1788}
1789
1790
1791
1792#ifdef HAVE_NAMESPACES
1793}
1794# ifndef NOT_USING_DJVU_NAMESPACE
1795using namespace DJVU;
1796# endif
1797#endif
1798#endif // NEED_DECODER_ONLY
1799
Note: See TracBrowser for help on using the repository browser.