source: trunk/libdjvu/IW44EncodeCodec.cpp @ 136

Last change on this file since 136 was 17, checked in by Eugene Romanenko, 16 years ago

update makefiles, remove absolute paths, update djvulibre to version 3.5.17

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