source: trunk/libdjvu/IW44Image.cpp @ 15

Last change on this file since 15 was 15, checked in by Eugene Romanenko, 15 years ago

needed libs update

File size: 52.3 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: IW44Image.cpp,v 1.11 2004/08/06 15:11:29 leonb Exp $
55// $Name: release_3_5_16 $
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
66// From: Leon Bottou, 1/31/2002
67// Lizardtech has split this file into a decoder and an encoder.
68// Only superficial changes.  The meat is mine.
69
70#define IW44IMAGE_IMPLIMENTATION /* */
71// -------------------^  not my spelling mistake (Leon Bottou)
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
94#ifdef HAVE_NAMESPACES
95namespace DJVU {
96# ifdef NOT_DEFINED // Just to fool emacs c++ mode
97}
98#endif
99#endif
100
101
102#define IWALLOCSIZE    4080
103#define IWCODEC_MAJOR     1
104#define IWCODEC_MINOR     2
105#define DECIBEL_PRUNE   5.0
106
107
108//////////////////////////////////////////////////////
109// WAVELET DECOMPOSITION CONSTANTS
110//////////////////////////////////////////////////////
111
112// Parameters for IW44 wavelet.
113// - iw_quant: quantization for all 16 sub-bands
114// - iw_norm: norm of all wavelets (for db estimation)
115// - iw_border: pixel border required to run filters
116// - iw_shift: scale applied before decomposition
117
118
119static const int iw_quant[16] = {
120  0x004000, 
121  0x008000, 0x008000, 0x010000,
122  0x010000, 0x010000, 0x020000,
123  0x020000, 0x020000, 0x040000,
124  0x040000, 0x040000, 0x080000, 
125  0x040000, 0x040000, 0x080000
126};
127
128static const float iw_norm[16] = {
129  2.627989e+03F,
130  1.832893e+02F, 1.832959e+02F, 5.114690e+01F,
131  4.583344e+01F, 4.583462e+01F, 1.279225e+01F,
132  1.149671e+01F, 1.149712e+01F, 3.218888e+00F,
133  2.999281e+00F, 2.999476e+00F, 8.733161e-01F,
134  1.074451e+00F, 1.074511e+00F, 4.289318e-01F
135};
136
137static const int iw_border = 3;
138static const int iw_shift  = 6;
139static const int iw_round  = (1<<(iw_shift-1));
140
141class IW44Image::Codec::Decode : public IW44Image::Codec
142{
143public:
144  // Construction
145  Decode(IW44Image::Map &map) : Codec(map) {}
146  // Coding
147  virtual int code_slice(ZPCodec &zp);
148};
149
150//////////////////////////////////////////////////////
151// MMX IMPLEMENTATION HELPERS
152//////////////////////////////////////////////////////
153
154
155// Note:
156// MMX implementation for vertical transforms only.
157// Speedup is basically related to faster memory transfer
158// The IW44 transform is not CPU bound, it is memory bound.
159
160#ifdef MMX
161
162static const short w9[]  = {9,9,9,9};
163static const short w1[]  = {1,1,1,1};
164static const int   d8[]  = {8,8};
165static const int   d16[] = {16,16};
166
167static void
168mmx_bv_1 ( short* &q, short* e, int s, int s3 )
169{
170  while (q<e && (((long)q)&0x7))
171    {
172      int a = (int)q[-s] + (int)q[s];
173      int b = (int)q[-s3] + (int)q[s3];
174      *q -= (((a<<3)+a-b+16)>>5);
175      q ++;
176  }
177  while (q+3 < e)
178    {
179      MMXar( movq,       q-s,mm0);  // MM0=[ b3, b2, b1, b0 ]
180      MMXar( movq,       q+s,mm2);  // MM2=[ c3, c2, c1, c0 ]
181      MMXrr( movq,       mm0,mm1); 
182      MMXrr( punpcklwd,  mm2,mm0);  // MM0=[ c1, b1, c0, b0 ]
183      MMXrr( punpckhwd,  mm2,mm1);  // MM1=[ c3, b3, c2, b2 ]
184      MMXar( pmaddwd,    w9,mm0);   // MM0=[ (c1+b1)*9, (c0+b0)*9 ]
185      MMXar( pmaddwd,    w9,mm1);   // MM1=[ (c3+b3)*9, (c2+b2)*9 ]
186      MMXar( movq,       q-s3,mm2);
187      MMXar( movq,       q+s3,mm4);
188      MMXrr( movq,       mm2,mm3);
189      MMXrr( punpcklwd,  mm4,mm2);  // MM2=[ d1, a1, d0, a0 ]
190      MMXrr( punpckhwd,  mm4,mm3);  // MM3=[ d3, a3, d2, a2 ]
191      MMXar( pmaddwd,    w1,mm2);   // MM2=[ (a1+d1)*1, (a0+d0)*1 ]
192      MMXar( pmaddwd,    w1,mm3);   // MM3=[ (a3+d3)*1, (a2+d2)*1 ]
193      MMXar( paddd,      d16,mm0);
194      MMXar( paddd,      d16,mm1);
195      MMXrr( psubd,      mm2,mm0);  // MM0=[ (c1+b1)*9-a1-d1+8, ...
196      MMXrr( psubd,      mm3,mm1);  // MM1=[ (c3+b3)*9-a3-d3+8, ...
197      MMXir( psrad,      5,mm0);
198      MMXar( movq,       q,mm7);    // MM7=[ p3,p2,p1,p0 ]
199      MMXir( psrad,      5,mm1);
200      MMXrr( packssdw,   mm1,mm0);  // MM0=[ x3,x2,x1,x0 ]
201      MMXrr( psubw,      mm0,mm7);  // MM7=[ p3-x3, p2-x2, ... ]
202      MMXra( movq,       mm7,q);
203#if defined(_MSC_VER) && defined(_DEBUG)
204      MMXemms;
205#endif
206      q += 4;
207    }
208}
209
210
211static void
212mmx_bv_2 ( short* &q, short* e, int s, int s3 )
213{
214  while (q<e && (((long)q)&0x7))
215    {
216      int a = (int)q[-s] + (int)q[s];
217      int b = (int)q[-s3] + (int)q[s3];
218      *q += (((a<<3)+a-b+8)>>4);
219      q ++;
220    }
221  while (q+3 < e)
222    {
223      MMXar( movq,       q-s,mm0);  // MM0=[ b3, b2, b1, b0 ]
224      MMXar( movq,       q+s,mm2);  // MM2=[ c3, c2, c1, c0 ]
225      MMXrr( movq,       mm0,mm1); 
226      MMXrr( punpcklwd,  mm2,mm0);  // MM0=[ c1, b1, c0, b0 ]
227      MMXrr( punpckhwd,  mm2,mm1);  // MM1=[ c3, b3, c2, b2 ]
228      MMXar( pmaddwd,    w9,mm0);   // MM0=[ (c1+b1)*9, (c0+b0)*9 ]
229      MMXar( pmaddwd,    w9,mm1);   // MM1=[ (c3+b3)*9, (c2+b2)*9 ]
230      MMXar( movq,       q-s3,mm2);
231      MMXar( movq,       q+s3,mm4);
232      MMXrr( movq,       mm2,mm3);
233      MMXrr( punpcklwd,  mm4,mm2);  // MM2=[ d1, a1, d0, a0 ]
234      MMXrr( punpckhwd,  mm4,mm3);  // MM3=[ d3, a3, d2, a2 ]
235      MMXar( pmaddwd,    w1,mm2);   // MM2=[ (a1+d1)*1, (a0+d0)*1 ]
236      MMXar( pmaddwd,    w1,mm3);   // MM3=[ (a3+d3)*1, (a2+d2)*1 ]
237      MMXar( paddd,      d8,mm0);
238      MMXar( paddd,      d8,mm1);
239      MMXrr( psubd,      mm2,mm0);  // MM0=[ (c1+b1)*9-a1-d1+8, ...
240      MMXrr( psubd,      mm3,mm1);  // MM1=[ (c3+b3)*9-a3-d3+8, ...
241      MMXir( psrad,      4,mm0);
242      MMXar( movq,       q,mm7);    // MM7=[ p3,p2,p1,p0 ]
243      MMXir( psrad,      4,mm1);
244      MMXrr( packssdw,   mm1,mm0);  // MM0=[ x3,x2,x1,x0 ]
245      MMXrr( paddw,      mm0,mm7);  // MM7=[ p3+x3, p2+x2, ... ]
246      MMXra( movq,       mm7,q);
247#if defined(_MSC_VER) && defined(_DEBUG)
248      MMXemms;
249#endif
250      q += 4;
251    }
252}
253#endif /* MMX */
254
255static void 
256filter_bv(short *p, int w, int h, int rowsize, int scale)
257{
258  int y = 0;
259  int s = scale*rowsize;
260  int s3 = s+s+s;
261  h = ((h-1)/scale)+1;
262  while (y-3 < h)
263    {
264      // 1-Lifting
265      {
266        short *q = p;
267        short *e = q+w;
268        if (y>=3 && y+3<h)
269          {
270            // Generic case
271#ifdef MMX
272            if (scale==1 && MMXControl::mmxflag>0)
273              mmx_bv_1(q, e, s, s3);
274#endif
275            while (q<e)
276              {
277                int a = (int)q[-s] + (int)q[s];
278                int b = (int)q[-s3] + (int)q[s3];
279                *q -= (((a<<3)+a-b+16)>>5);
280                q += scale;
281              }
282          }
283        else if (y<h)
284          {
285            // Special cases
286            short *q1 = (y+1<h ? q+s : 0);
287            short *q3 = (y+3<h ? q+s3 : 0);
288            if (y>=3)
289              {
290                while (q<e)
291                  {
292                    int a = (int)q[-s] + (q1 ? (int)(*q1) : 0);
293                    int b = (int)q[-s3] + (q3 ? (int)(*q3) : 0);
294                    *q -= (((a<<3)+a-b+16)>>5);
295                    q += scale;
296                    if (q1) q1 += scale;
297                    if (q3) q3 += scale;
298                  }
299              }
300            else if (y>=1)
301              {
302                while (q<e)
303                  {
304                    int a = (int)q[-s] + (q1 ? (int)(*q1) : 0);
305                    int b = (q3 ? (int)(*q3) : 0);
306                    *q -= (((a<<3)+a-b+16)>>5);
307                    q += scale;
308                    if (q1) q1 += scale;
309                    if (q3) q3 += scale;
310                  }
311              }
312            else
313              {
314                while (q<e)
315                  {
316                    int a = (q1 ? (int)(*q1) : 0);
317                    int b = (q3 ? (int)(*q3) : 0);
318                    *q -= (((a<<3)+a-b+16)>>5);
319                    q += scale;
320                    if (q1) q1 += scale;
321                    if (q3) q3 += scale;
322                  }
323              }
324          }
325      }
326      // 2-Interpolation
327      {
328        short *q = p-s3;
329        short *e = q+w;
330        if (y>=6 && y<h)
331          {
332            // Generic case
333#ifdef MMX
334            if (scale==1 && MMXControl::mmxflag>0)
335              mmx_bv_2(q, e, s, s3);
336#endif
337            while (q<e)
338              {
339                int a = (int)q[-s] + (int)q[s];
340                int b = (int)q[-s3] + (int)q[s3];
341                *q += (((a<<3)+a-b+8)>>4);
342                q += scale;
343              }
344          }
345        else if (y>=3)
346          {
347            // Special cases
348            short *q1 = (y-2<h ? q+s : q-s);
349            while (q<e)
350              {
351                int a = (int)q[-s] + (int)(*q1);
352                *q += ((a+1)>>1);
353                q += scale;
354                q1 += scale;
355              }
356          }
357      }
358      y += 2;
359      p += s+s;
360    }
361}
362
363static void 
364filter_bh(short *p, int w, int h, int rowsize, int scale)
365{
366  int y = 0;
367  int s = scale;
368  int s3 = s+s+s;
369  rowsize *= scale;
370  while (y<h)
371    {
372      short *q = p;
373      short *e = p+w;
374      int a0=0, a1=0, a2=0, a3=0;
375      int b0=0, b1=0, b2=0, b3=0;
376      if (q<e)
377        {
378          // Special case:  x=0
379          if (q+s < e)
380            a2 = q[s];
381          if (q+s3 < e)
382            a3 = q[s3];
383          b2 = b3 = q[0] - ((((a1+a2)<<3)+(a1+a2)-a0-a3+16) >> 5);
384          q[0] = b3;
385          q += s+s;
386        }
387      if (q<e)
388        {
389          // Special case:  x=2
390          a0 = a1;
391          a1 = a2;
392          a2 = a3;
393          if (q+s3 < e)
394            a3 = q[s3];
395          b3 = q[0] - ((((a1+a2)<<3)+(a1+a2)-a0-a3+16) >> 5);
396          q[0] = b3;
397          q += s+s;
398        }
399      if (q<e)
400        {
401          // Special case:  x=4
402          b1 = b2;
403          b2 = b3;
404          a0 = a1;
405          a1 = a2;
406          a2 = a3;
407          if (q+s3 < e)
408            a3 = q[s3];
409          b3 = q[0] - ((((a1+a2)<<3)+(a1+a2)-a0-a3+16) >> 5);
410          q[0] = b3;
411          q[-s3] = q[-s3] + ((b1+b2+1)>>1);
412          q += s+s;
413        }
414      while (q+s3 < e)
415        {
416          // Generic case
417          a0=a1; 
418          a1=a2; 
419          a2=a3;
420          a3=q[s3];
421          b0=b1; 
422          b1=b2; 
423          b2=b3;
424          b3 = q[0] - ((((a1+a2)<<3)+(a1+a2)-a0-a3+16) >> 5);
425          q[0] = b3;
426          q[-s3] = q[-s3] + ((((b1+b2)<<3)+(b1+b2)-b0-b3+8) >> 4);
427          q += s+s;
428        }
429      while (q < e)
430        {
431          // Special case:  w-3 <= x < w
432          a0=a1;
433          a1=a2; 
434          a2=a3;
435          a3=0;
436          b0=b1; 
437          b1=b2; 
438          b2=b3;
439          b3 = q[0] - ((((a1+a2)<<3)+(a1+a2)-a0-a3+16) >> 5);
440          q[0] = b3;
441          q[-s3] = q[-s3] + ((((b1+b2)<<3)+(b1+b2)-b0-b3+8) >> 4);
442          q += s+s;
443        }
444      while (q-s3 < e)
445        {
446          // Special case  w <= x < w+3
447          b0=b1; 
448          b1=b2; 
449          b2=b3;
450          if (q-s3 >= p)
451            q[-s3] = q[-s3] + ((b1+b2+1)>>1);
452          q += s+s;
453        }
454      y += scale;
455      p += rowsize;
456    }
457}
458
459
460//////////////////////////////////////////////////////
461// REPRESENTATION OF WAVELET DECOMPOSED IMAGES
462//////////////////////////////////////////////////////
463
464
465
466//---------------------------------------------------------------
467// Zig zag location in a 1024 liftblock.
468// These numbers have been generated with the following program:
469//
470// int x=0, y=0;
471// for (int i=0; i<5; i++) {
472//   x = (x<<1) | (n&1);  n >>= 1;
473//   y = (y<<1) | (n&1);  n >>= 1;
474// }
475
476
477static int zigzagloc[1024] = {
478   0,  16, 512, 528,   8,  24, 520, 536, 256, 272, 768, 784, 264, 280, 776, 792,
479   4,  20, 516, 532,  12,  28, 524, 540, 260, 276, 772, 788, 268, 284, 780, 796,
480 128, 144, 640, 656, 136, 152, 648, 664, 384, 400, 896, 912, 392, 408, 904, 920,
481 132, 148, 644, 660, 140, 156, 652, 668, 388, 404, 900, 916, 396, 412, 908, 924,
482   2,  18, 514, 530,  10,  26, 522, 538, 258, 274, 770, 786, 266, 282, 778, 794,
483   6,  22, 518, 534,  14,  30, 526, 542, 262, 278, 774, 790, 270, 286, 782, 798,
484 130, 146, 642, 658, 138, 154, 650, 666, 386, 402, 898, 914, 394, 410, 906, 922,
485 134, 150, 646, 662, 142, 158, 654, 670, 390, 406, 902, 918, 398, 414, 910, 926,
486  64,  80, 576, 592,  72,  88, 584, 600, 320, 336, 832, 848, 328, 344, 840, 856,
487  68,  84, 580, 596,  76,  92, 588, 604, 324, 340, 836, 852, 332, 348, 844, 860,
488 192, 208, 704, 720, 200, 216, 712, 728, 448, 464, 960, 976, 456, 472, 968, 984,
489 196, 212, 708, 724, 204, 220, 716, 732, 452, 468, 964, 980, 460, 476, 972, 988,
490  66,  82, 578, 594,  74,  90, 586, 602, 322, 338, 834, 850, 330, 346, 842, 858,
491  70,  86, 582, 598,  78,  94, 590, 606, 326, 342, 838, 854, 334, 350, 846, 862,
492 194, 210, 706, 722, 202, 218, 714, 730, 450, 466, 962, 978, 458, 474, 970, 986,
493 198, 214, 710, 726, 206, 222, 718, 734, 454, 470, 966, 982, 462, 478, 974, 990, // 255
494   1,  17, 513, 529,   9,  25, 521, 537, 257, 273, 769, 785, 265, 281, 777, 793,
495   5,  21, 517, 533,  13,  29, 525, 541, 261, 277, 773, 789, 269, 285, 781, 797,
496 129, 145, 641, 657, 137, 153, 649, 665, 385, 401, 897, 913, 393, 409, 905, 921,
497 133, 149, 645, 661, 141, 157, 653, 669, 389, 405, 901, 917, 397, 413, 909, 925,
498   3,  19, 515, 531,  11,  27, 523, 539, 259, 275, 771, 787, 267, 283, 779, 795,
499   7,  23, 519, 535,  15,  31, 527, 543, 263, 279, 775, 791, 271, 287, 783, 799,
500 131, 147, 643, 659, 139, 155, 651, 667, 387, 403, 899, 915, 395, 411, 907, 923,
501 135, 151, 647, 663, 143, 159, 655, 671, 391, 407, 903, 919, 399, 415, 911, 927,
502  65,  81, 577, 593,  73,  89, 585, 601, 321, 337, 833, 849, 329, 345, 841, 857,
503  69,  85, 581, 597,  77,  93, 589, 605, 325, 341, 837, 853, 333, 349, 845, 861,
504 193, 209, 705, 721, 201, 217, 713, 729, 449, 465, 961, 977, 457, 473, 969, 985,
505 197, 213, 709, 725, 205, 221, 717, 733, 453, 469, 965, 981, 461, 477, 973, 989,
506  67,  83, 579, 595,  75,  91, 587, 603, 323, 339, 835, 851, 331, 347, 843, 859,
507  71,  87, 583, 599,  79,  95, 591, 607, 327, 343, 839, 855, 335, 351, 847, 863,
508 195, 211, 707, 723, 203, 219, 715, 731, 451, 467, 963, 979, 459, 475, 971, 987,
509 199, 215, 711, 727, 207, 223, 719, 735, 455, 471, 967, 983, 463, 479, 975, 991, // 511
510  32,  48, 544, 560,  40,  56, 552, 568, 288, 304, 800, 816, 296, 312, 808, 824,
511  36,  52, 548, 564,  44,  60, 556, 572, 292, 308, 804, 820, 300, 316, 812, 828,
512 160, 176, 672, 688, 168, 184, 680, 696, 416, 432, 928, 944, 424, 440, 936, 952,
513 164, 180, 676, 692, 172, 188, 684, 700, 420, 436, 932, 948, 428, 444, 940, 956,
514  34,  50, 546, 562,  42,  58, 554, 570, 290, 306, 802, 818, 298, 314, 810, 826,
515  38,  54, 550, 566,  46,  62, 558, 574, 294, 310, 806, 822, 302, 318, 814, 830,
516 162, 178, 674, 690, 170, 186, 682, 698, 418, 434, 930, 946, 426, 442, 938, 954,
517 166, 182, 678, 694, 174, 190, 686, 702, 422, 438, 934, 950, 430, 446, 942, 958,
518  96, 112, 608, 624, 104, 120, 616, 632, 352, 368, 864, 880, 360, 376, 872, 888,
519 100, 116, 612, 628, 108, 124, 620, 636, 356, 372, 868, 884, 364, 380, 876, 892,
520 224, 240, 736, 752, 232, 248, 744, 760, 480, 496, 992,1008, 488, 504,1000,1016,
521 228, 244, 740, 756, 236, 252, 748, 764, 484, 500, 996,1012, 492, 508,1004,1020,
522  98, 114, 610, 626, 106, 122, 618, 634, 354, 370, 866, 882, 362, 378, 874, 890,
523 102, 118, 614, 630, 110, 126, 622, 638, 358, 374, 870, 886, 366, 382, 878, 894,
524 226, 242, 738, 754, 234, 250, 746, 762, 482, 498, 994,1010, 490, 506,1002,1018,
525 230, 246, 742, 758, 238, 254, 750, 766, 486, 502, 998,1014, 494, 510,1006,1022, // 767
526  33,  49, 545, 561,  41,  57, 553, 569, 289, 305, 801, 817, 297, 313, 809, 825,
527  37,  53, 549, 565,  45,  61, 557, 573, 293, 309, 805, 821, 301, 317, 813, 829,
528 161, 177, 673, 689, 169, 185, 681, 697, 417, 433, 929, 945, 425, 441, 937, 953,
529 165, 181, 677, 693, 173, 189, 685, 701, 421, 437, 933, 949, 429, 445, 941, 957,
530  35,  51, 547, 563,  43,  59, 555, 571, 291, 307, 803, 819, 299, 315, 811, 827,
531  39,  55, 551, 567,  47,  63, 559, 575, 295, 311, 807, 823, 303, 319, 815, 831,
532 163, 179, 675, 691, 171, 187, 683, 699, 419, 435, 931, 947, 427, 443, 939, 955,
533 167, 183, 679, 695, 175, 191, 687, 703, 423, 439, 935, 951, 431, 447, 943, 959,
534  97, 113, 609, 625, 105, 121, 617, 633, 353, 369, 865, 881, 361, 377, 873, 889,
535 101, 117, 613, 629, 109, 125, 621, 637, 357, 373, 869, 885, 365, 381, 877, 893,
536 225, 241, 737, 753, 233, 249, 745, 761, 481, 497, 993,1009, 489, 505,1001,1017,
537 229, 245, 741, 757, 237, 253, 749, 765, 485, 501, 997,1013, 493, 509,1005,1021,
538  99, 115, 611, 627, 107, 123, 619, 635, 355, 371, 867, 883, 363, 379, 875, 891,
539 103, 119, 615, 631, 111, 127, 623, 639, 359, 375, 871, 887, 367, 383, 879, 895,
540 227, 243, 739, 755, 235, 251, 747, 763, 483, 499, 995,1011, 491, 507,1003,1019,
541 231, 247, 743, 759, 239, 255, 751, 767, 487, 503, 999,1015, 495, 511,1007,1023, // 1023
542};
543
544//---------------------------------------------------------------
545// *** Class IW44Image::Alloc [declaration]
546
547struct IW44Image::Alloc // DJVU_CLASS
548{
549  Alloc *next;
550  short data[IWALLOCSIZE];
551};
552
553//---------------------------------------------------------------
554// *** Class IW44Image::Block [implementation]
555
556
557IW44Image::Block::Block(void)
558{
559  pdata[0] = pdata[1] = pdata[2] = pdata[3] = 0;
560}
561
562void 
563IW44Image::Block::zero(int n)
564{
565  if (pdata[n>>4])
566    pdata[n>>4][n&15] = 0;
567}
568
569void 
570IW44Image::Block::read_liftblock(const short *coeff, IW44Image::Map *map)
571{
572  int n=0;
573  for (int n1=0; n1<64; n1++)
574    {
575      short *d = data(n1,map);
576      for (int n2=0; n2<16; n2++,n++)
577        d[n2] = coeff[zigzagloc[n]];
578    }
579}
580
581void 
582IW44Image::Block::write_liftblock(short *coeff, int bmin, int bmax) const
583{
584  int n = bmin<<4;
585  memset(coeff, 0, 1024*sizeof(short));
586  for (int n1=bmin; n1<bmax; n1++)
587    {
588      const short *d = data(n1);
589      if (d == 0)
590        n += 16;
591      else
592        for (int n2=0; n2<16; n2++,n++)
593          coeff[zigzagloc[n]] = d[n2];
594    }
595}
596
597//---------------------------------------------------------------
598// *** Class IW44Image::Map [implementation]
599
600
601IW44Image::Map::Map(int w, int h)
602  :  blocks(0), iw(w), ih(h), chain(0)
603{
604  bw = (w+0x20-1) & ~0x1f;
605  bh = (h+0x20-1) & ~0x1f;
606  nb = (bw * bh) / (32 * 32);
607  blocks = new IW44Image::Block[nb];
608  top = IWALLOCSIZE;
609}
610
611IW44Image::Map::~Map()
612{
613  while (chain)
614    {
615      IW44Image::Alloc *next = chain->next;
616      delete chain;
617      chain = next;
618    }
619  delete [] blocks;
620}
621
622short *
623IW44Image::Map::alloc(int n)
624{
625  if (top+n > IWALLOCSIZE)
626    {
627      IW44Image::Alloc *n = new IW44Image::Alloc;
628      n->next = chain;
629      chain = n;
630      top = 0;
631    }
632  short *ans = chain->data + top;
633  top += n;
634  memset((void*)ans, 0, sizeof(short)*n);
635  return ans;
636}
637
638short **
639IW44Image::Map::allocp(int n)
640{
641  // Allocate enough room for pointers plus alignment
642  short *p = alloc( (n+1) * sizeof(short*) / sizeof(short) );
643  // Align on pointer size
644  while ( ((long)p) & (sizeof(short*)-1) )
645    p += 1;
646  // Cast and return
647  return (short**)p;
648}
649
650int 
651IW44Image::Map::get_bucket_count(void) const
652{
653  int buckets = 0;
654  for (int blockno=0; blockno<nb; blockno++)
655    for (int buckno=0; buckno<64; buckno++)
656      if (blocks[blockno].data(buckno))
657        buckets += 1;
658  return buckets;
659}
660
661unsigned int 
662IW44Image::Map::get_memory_usage(void) const
663{
664  unsigned int usage = sizeof(Map);
665  usage += sizeof(IW44Image::Block) * nb;
666  for (IW44Image::Alloc *n = chain; n; n=n->next)
667    usage += sizeof(IW44Image::Alloc);
668  return usage;
669}
670
671
672
673
674void 
675IW44Image::Map::image(signed char *img8, int rowsize, int pixsep, int fast)
676{
677  // Allocate reconstruction buffer
678  short *data16;
679  GPBuffer<short> gdata16(data16,bw*bh);
680  // Copy coefficients
681  int i;
682  short *p = data16;
683  const IW44Image::Block *block = blocks;
684  for (i=0; i<bh; i+=32)
685    {
686      for (int j=0; j<bw; j+=32)
687        {
688          short liftblock[1024];
689          // transfer into IW44Image::Block (apply zigzag and scaling)
690          block->write_liftblock(liftblock);
691          block++;
692          // transfer into coefficient matrix at (p+j)
693          short *pp = p + j;
694          short *pl = liftblock;
695          for (int ii=0; ii<32; ii++, pp+=bw,pl+=32)
696            memcpy((void*)pp, (void*)pl, 32*sizeof(short));
697        }
698      // next row of blocks
699      p += 32*bw;
700    }
701  // Reconstruction
702  if (fast)
703    {
704      IW44Image::Transform::Decode::backward(data16, iw, ih, bw, 32, 2); 
705      p = data16;
706      for (i=0; i<bh; i+=2,p+=bw)
707        for (int jj=0; jj<bw; jj+=2,p+=2)
708          p[bw] = p[bw+1] = p[1] = p[0];
709    }
710  else
711    {
712      IW44Image::Transform::Decode::backward(data16, iw, ih, bw, 32, 1); 
713    }
714  // Copy result into image
715  p = data16;
716  signed char *row = img8; 
717  for (i=0; i<ih; i++)
718    {
719      signed char *pix = row;
720      for (int j=0; j<iw; j+=1,pix+=pixsep)
721        {
722          int x = (p[j] + iw_round) >> iw_shift;
723          if (x < -128)
724            x = -128;
725          else if (x > 127)
726            x = 127;
727          *pix = x;
728        }
729      row += rowsize;
730      p += bw;
731    }
732}
733
734void 
735IW44Image::Map::image(int subsample, const GRect &rect, 
736              signed char *img8, int rowsize, int pixsep, int fast)
737{
738  int i;
739  // Compute number of decomposition levels
740  int nlevel = 0;
741  while (nlevel<5 && (32>>nlevel)>subsample)
742    nlevel += 1;
743  int boxsize = 1<<nlevel;
744  // Parameter check
745  if (subsample!=(32>>nlevel))
746    G_THROW( ERR_MSG("IW44Image.sample_factor") );
747  if (rect.isempty())
748    G_THROW( ERR_MSG("IW44Image.empty_rect") );   
749  GRect irect(0,0,(iw+subsample-1)/subsample,(ih+subsample-1)/subsample);
750  if (rect.xmin<0 || rect.ymin<0 || rect.xmax>irect.xmax || rect.ymax>irect.ymax)
751    G_THROW( ERR_MSG("IW44Image.bad_rect") );
752  // Multiresolution rectangles
753  // -- needed[i] tells which coeffs are required for the next step
754  // -- recomp[i] tells which coeffs need to be computed at this level
755  GRect needed[8];
756  GRect recomp[8];
757  int r = 1;
758  needed[nlevel] = rect;
759  recomp[nlevel] = rect;
760  for (i=nlevel-1; i>=0; i--)
761    {
762      needed[i] = recomp[i+1];
763      needed[i].inflate(iw_border*r, iw_border*r);
764      needed[i].intersect(needed[i], irect);
765      r += r;
766      recomp[i].xmin = (needed[i].xmin + r-1) & ~(r-1);
767      recomp[i].xmax = (needed[i].xmax) & ~(r-1);
768      recomp[i].ymin = (needed[i].ymin + r-1) & ~(r-1);
769      recomp[i].ymax = (needed[i].ymax) & ~(r-1);
770    }
771  // Working rectangle
772  // -- a rectangle large enough to hold all the data
773  GRect work;
774  work.xmin = (needed[0].xmin) & ~(boxsize-1);
775  work.ymin = (needed[0].ymin) & ~(boxsize-1);
776  work.xmax = ((needed[0].xmax-1) & ~(boxsize-1) ) + boxsize;
777  work.ymax = ((needed[0].ymax-1) & ~(boxsize-1) ) + boxsize;
778  // -- allocate work buffer
779  int dataw = work.xmax - work.xmin;     // Note: cannot use inline width() or height()
780  int datah = work.ymax - work.ymin;     // because Intel C++ compiler optimizes it wrong !
781  short *data;
782  GPBuffer<short> gdata(data,dataw*datah);
783  // Fill working rectangle
784  // -- loop over liftblocks rows
785  short *ldata = data;
786  int blkw = (bw>>5);
787  const IW44Image::Block *lblock = blocks + (work.ymin>>nlevel)*blkw + (work.xmin>>nlevel);
788  for (int by=work.ymin; by<work.ymax; by+=boxsize)
789    {
790      // -- loop over liftblocks in row
791      const IW44Image::Block *block = lblock;
792      short *rdata = ldata;
793      for (int bx=work.xmin; bx<work.xmax; bx+=boxsize)       
794        {
795          // -- decide how much to load
796          int mlevel = nlevel;
797          if (nlevel>2)
798            if (bx+31<needed[2].xmin || bx>needed[2].xmax ||
799                by+31<needed[2].ymin || by>needed[2].ymax )
800              mlevel = 2;
801          int bmax   = ((1<<(mlevel+mlevel))+15)>>4;
802          int ppinc  = (1<<(nlevel-mlevel));
803          int ppmod1 = (dataw<<(nlevel-mlevel));
804          int ttmod0 = (32 >> mlevel);
805          int ttmod1 = (ttmod0 << 5);
806          // -- get current block
807          short liftblock[1024];
808          block->write_liftblock(liftblock, 0, bmax );
809          // -- copy liftblock into image
810          short *tt = liftblock;
811          short *pp = rdata;
812          for (int ii=0; ii<boxsize; ii+=ppinc,pp+=ppmod1,tt+=ttmod1-32)
813            for (int jj=0; jj<boxsize; jj+=ppinc,tt+=ttmod0)
814              pp[jj] = *tt;
815          // -- next block in row
816          rdata += boxsize;
817          block += 1;
818        }
819      // -- next row of blocks
820      ldata += dataw << nlevel;
821      lblock += blkw;
822    }
823  // Perform reconstruction
824  r = boxsize;
825  for (i=0; i<nlevel; i++)
826    {
827      GRect comp = needed[i];
828      comp.xmin = comp.xmin & ~(r-1);
829      comp.ymin = comp.ymin & ~(r-1);
830      comp.translate(-work.xmin, -work.ymin);
831      // Fast mode shortcuts finer resolution
832      if (fast && i>=4) 
833        {
834          short *pp = data + comp.ymin*dataw;
835          for (int ii=comp.ymin; ii<comp.ymax; ii+=2, pp+=dataw+dataw)
836            for (int jj=comp.xmin; jj<comp.xmax; jj+=2)
837              pp[jj+dataw] = pp[jj+dataw+1] = pp[jj+1] = pp[jj];
838          break;
839        }
840      else
841        {
842          short *pp = data + comp.ymin*dataw + comp.xmin;
843          IW44Image::Transform::Decode::backward(pp, comp.width(), comp.height(), dataw, r, r>>1);
844        }
845      r = r>>1;
846    }
847  // Copy result into image
848  GRect nrect = rect;
849  nrect.translate(-work.xmin, -work.ymin);
850  short *p = data + nrect.ymin*dataw;
851  signed char *row = img8; 
852  for (i=nrect.ymin; i<nrect.ymax; i++)
853    {
854      int j;
855      signed char *pix = row;
856      for (j=nrect.xmin; j<nrect.xmax; j+=1,pix+=pixsep)
857        {
858          int x = (p[j] + iw_round) >> iw_shift;
859          if (x < -128)
860            x = -128;
861          else if (x > 127)
862            x = 127;
863          *pix = x;
864        }
865      row += rowsize;
866      p += dataw;
867    }
868}
869
870
871
872
873//////////////////////////////////////////////////////
874// ENCODING/DECODING WAVELET COEFFICIENTS
875//    USING HIERARCHICAL SET DIFFERENCE
876//////////////////////////////////////////////////////
877
878
879//-----------------------------------------------
880// Class IW44Image::Codec [implementation]
881// Maintains information shared while encoding or decoding
882
883
884// Constant
885
886static const struct { int start; int size; } 
887bandbuckets[] = 
888{
889  // Code first bucket and number of buckets in each band
890  { 0, 1 }, // -- band zero contains all lores info
891  { 1, 1 }, { 2, 1 }, { 3, 1 }, 
892  { 4, 4 }, { 8, 4 }, { 12,4 }, 
893  { 16,16 }, { 32,16 }, { 48,16 }, 
894};
895
896
897// IW44Image::Codec constructor
898
899IW44Image::Codec::Codec(IW44Image::Map &xmap)
900  : map(xmap), 
901    curband(0),
902    curbit(1)
903{
904  // Initialize quantification
905  int  j;
906  int  i = 0;
907  const int *q = iw_quant;
908  // -- lo coefficients
909  for (j=0; i<4; j++)
910    quant_lo[i++] = *q++;
911  for (j=0; j<4; j++)
912    quant_lo[i++] = *q;
913  q += 1;
914  for (j=0; j<4; j++)
915    quant_lo[i++] = *q;
916  q += 1;
917  for (j=0; j<4; j++)
918    quant_lo[i++] = *q;
919  q += 1;
920  // -- hi coefficients
921  quant_hi[0] = 0;
922  for (j=1; j<10; j++)
923    quant_hi[j] = *q++;
924  // Initialize coding contexts
925  memset((void*)ctxStart, 0, sizeof(ctxStart));
926  memset((void*)ctxBucket, 0, sizeof(ctxBucket));
927  ctxMant = 0;
928  ctxRoot = 0;
929}
930
931
932// IW44Image::Codec destructor
933
934IW44Image::Codec::~Codec() {}
935
936// is_null_slice
937// -- check if data can be produced for this band/mask
938// -- also fills the sure_zero array
939
940int 
941IW44Image::Codec::is_null_slice(int bit, int band)
942{
943  if (band == 0)
944    {
945      int is_null = 1;
946      for (int i=0; i<16; i++) 
947        {
948          int threshold = quant_lo[i];
949          coeffstate[i] = ZERO;
950          if (threshold>0 && threshold<0x8000)
951            {
952              coeffstate[i] = UNK;
953              is_null = 0;
954            }
955        }
956      return is_null;
957    }
958  else
959    {
960      int threshold = quant_hi[band];
961      return (! (threshold>0 && threshold<0x8000));
962    }
963}
964
965
966// code_slice
967// -- read/write a slice of datafile
968
969int
970IW44Image::Codec::Decode::code_slice(ZPCodec &zp)
971{
972  // Check that code_slice can still run
973  if (curbit < 0)
974    return 0;
975  // Perform coding
976  if (! is_null_slice(curbit, curband))
977    {
978      for (int blockno=0; blockno<map.nb; blockno++)
979        {
980          int fbucket = bandbuckets[curband].start;
981          int nbucket = bandbuckets[curband].size;
982          decode_buckets(zp, curbit, curband, 
983                           map.blocks[blockno], 
984                           fbucket, nbucket);
985        }
986    }
987  return finish_code_slice(zp);
988}
989
990// code_slice
991// -- read/write a slice of datafile
992
993int
994IW44Image::Codec::finish_code_slice(ZPCodec &zp)
995{
996  // Reduce quantization threshold
997  quant_hi[curband] = quant_hi[curband] >> 1;
998  if (curband == 0)
999    for (int i=0; i<16; i++) 
1000      quant_lo[i] = quant_lo[i] >> 1;
1001  // Proceed to the next slice
1002  if (++curband >= (int)(sizeof(bandbuckets)/sizeof(bandbuckets[0])))
1003    {
1004      curband = 0;
1005      curbit += 1;
1006      if (quant_hi[(sizeof(bandbuckets)/sizeof(bandbuckets[0]))-1] == 0)
1007        {
1008          // All quantization thresholds are null
1009          curbit = -1;
1010          return 0;
1011        }
1012    }
1013  return 1;
1014}
1015
1016// decode_prepare
1017// -- prepare the states before decoding buckets
1018
1019int
1020IW44Image::Codec::decode_prepare(int fbucket, int nbucket, IW44Image::Block &blk)
1021{ 
1022  int bbstate = 0;
1023  char *cstate = coeffstate;
1024  if (fbucket)
1025    {
1026      // Band other than zero
1027      for (int buckno=0; buckno<nbucket; buckno++, cstate+=16)
1028        {
1029          int bstatetmp = 0;
1030          const short *pcoeff = blk.data(fbucket+buckno);
1031          if (! pcoeff)
1032            {
1033              // cstate[0..15] will be filled later
1034              bstatetmp = UNK;
1035            }
1036          else
1037            {
1038              for (int i=0; i<16; i++)
1039                {
1040                  int cstatetmp = UNK;
1041                  if (pcoeff[i])
1042                    cstatetmp = ACTIVE;
1043                  cstate[i] = cstatetmp;
1044                  bstatetmp |= cstatetmp;
1045                }
1046            }
1047          bucketstate[buckno] = bstatetmp;
1048          bbstate |= bstatetmp;
1049        }
1050    }
1051  else
1052    {
1053      // Band zero ( fbucket==0 implies band==zero and nbucket==1 )
1054      const short *pcoeff = blk.data(0);
1055      if (! pcoeff)
1056        {
1057          // cstate[0..15] will be filled later
1058          bbstate = UNK;     
1059        }
1060      else
1061        {
1062          for (int i=0; i<16; i++)
1063            {
1064              int cstatetmp = cstate[i];
1065              if (cstatetmp != ZERO)
1066                {
1067                  cstatetmp = UNK;
1068                  if (pcoeff[i])
1069                    cstatetmp = ACTIVE;
1070                }
1071              cstate[i] = cstatetmp;
1072              bbstate |= cstatetmp;
1073            }
1074        }
1075      bucketstate[0] = bbstate;
1076    }
1077  return bbstate;
1078}
1079
1080
1081// decode_buckets
1082// -- code a sequence of buckets in a given block
1083
1084void
1085IW44Image::Codec::decode_buckets(ZPCodec &zp, int bit, int band, 
1086                         IW44Image::Block &blk,
1087                         int fbucket, int nbucket)
1088{
1089  // compute state of all coefficients in all buckets
1090  int bbstate = decode_prepare(fbucket, nbucket, blk);
1091  // code root bit
1092  if ((nbucket<16) || (bbstate&ACTIVE))
1093    {
1094      bbstate |= NEW;
1095    }
1096  else if (bbstate & UNK)
1097    {
1098      if (zp.decoder(ctxRoot))
1099        bbstate |= NEW;
1100#ifdef TRACE
1101      DjVuPrintMessage("bbstate[bit=%d,band=%d] = %d\n", bit, band, bbstate);
1102#endif
1103    }
1104 
1105  // code bucket bits
1106  if (bbstate & NEW)
1107    for (int buckno=0; buckno<nbucket; buckno++)
1108      {
1109        // Code bucket bit
1110        if (bucketstate[buckno] & UNK)
1111          {
1112            // Context
1113            int ctx = 0;
1114#ifndef NOCTX_BUCKET_UPPER
1115            if (band>0)
1116              {
1117                int k = (fbucket+buckno)<<2;
1118                const short *b = blk.data(k>>4);
1119                if (b)
1120                  {
1121                    k = k & 0xf;
1122                    if (b[k])
1123                      ctx += 1;
1124                    if (b[k+1])
1125                      ctx += 1;
1126                    if (b[k+2])
1127                      ctx += 1;
1128                    if (ctx<3 && b[k+3])
1129                      ctx += 1;
1130                  }
1131              }
1132#endif // NOCTX_BUCKET_UPPER
1133#ifndef NOCTX_BUCKET_ACTIVE
1134            if (bbstate & ACTIVE)
1135              ctx |= 4; 
1136#endif
1137            // Code
1138            if (zp.decoder( ctxBucket[band][ctx] ))
1139              bucketstate[buckno] |= NEW;
1140#ifdef TRACE
1141            DjVuPrintMessage("  bucketstate[bit=%d,band=%d,buck=%d] = %d\n", 
1142                   bit, band, buckno, bucketstate[buckno]);
1143#endif
1144          }
1145      }
1146
1147  // code new active coefficient (with their sign)
1148  if (bbstate & NEW)
1149    {
1150      int thres = quant_hi[band];
1151      char *cstate = coeffstate;
1152      for (int buckno=0; buckno<nbucket; buckno++, cstate+=16)
1153        if (bucketstate[buckno] & NEW)
1154          {
1155            int i;
1156            short *pcoeff = (short*)blk.data(fbucket+buckno);
1157            if (!pcoeff)
1158              {
1159                pcoeff = blk.data(fbucket+buckno, &map);
1160                // time to fill cstate[0..15]
1161                if (fbucket == 0) // band zero
1162                  {
1163                    for (i=0; i<16; i++)
1164                      if (cstate[i] != ZERO)
1165                        cstate[i] = UNK;
1166                  }
1167                else
1168                  {
1169                    for (i=0; i<16; i++)
1170                      cstate[i] = UNK;
1171                  }
1172              }
1173#ifndef NOCTX_EXPECT
1174            int gotcha = 0;
1175            const int maxgotcha = 7;
1176            for (i=0; i<16; i++)
1177              if (cstate[i] & UNK)
1178                gotcha += 1;
1179#endif
1180            for (i=0; i<16; i++)
1181              {
1182                if (cstate[i] & UNK)
1183                  {
1184                    // find lores threshold
1185                    if (band == 0)
1186                      thres = quant_lo[i];
1187                    // prepare context
1188                    int ctx = 0;
1189#ifndef NOCTX_EXPECT
1190                    if (gotcha>=maxgotcha)
1191                      ctx = maxgotcha;
1192                    else
1193                      ctx = gotcha;
1194#endif
1195#ifndef NOCTX_ACTIVE
1196                    if (bucketstate[buckno] & ACTIVE)
1197                      ctx |= 8;
1198#endif
1199                    // code difference bit
1200                    if (zp.decoder( ctxStart[ctx] ))
1201                      {
1202                        cstate[i] |= NEW;
1203                        int halfthres = thres>>1;
1204                        int coeff = thres+halfthres-(halfthres>>2);
1205                        if (zp.IWdecoder())
1206                          pcoeff[i] = -coeff;
1207                        else
1208                          pcoeff[i] = coeff;
1209                      }
1210#ifndef NOCTX_EXPECT
1211                    if (cstate[i] & NEW)
1212                      gotcha = 0;
1213                    else if (gotcha > 0)
1214                      gotcha -= 1;
1215#endif
1216#ifdef TRACE
1217                    DjVuPrintMessage("    coeffstate[bit=%d,band=%d,buck=%d,c=%d] = %d\n", 
1218                           bit, band, buckno, i, cstate[i]);
1219#endif
1220                  }
1221              }
1222          }
1223    }
1224 
1225  // code mantissa bits
1226  if (bbstate & ACTIVE)
1227    {
1228      int thres = quant_hi[band];
1229      char *cstate = coeffstate;
1230      for (int buckno=0; buckno<nbucket; buckno++, cstate+=16)
1231        if (bucketstate[buckno] & ACTIVE)
1232          {
1233            short *pcoeff = (short*)blk.data(fbucket+buckno);
1234            for (int i=0; i<16; i++)
1235              if (cstate[i] & ACTIVE)
1236                {
1237                  int coeff = pcoeff[i];
1238                  if (coeff < 0)
1239                    coeff = -coeff;
1240                  // find lores threshold
1241                  if (band == 0)
1242                    thres = quant_lo[i];
1243                  // adjust coefficient
1244                  if (coeff <= 3*thres)
1245                    {
1246                      // second mantissa bit
1247                      coeff = coeff + (thres>>2);
1248                      if (zp.decoder(ctxMant))
1249                        coeff = coeff + (thres>>1);
1250                      else
1251                        coeff = coeff - thres + (thres>>1);
1252                    }
1253                  else
1254                    {
1255                      if (zp.IWdecoder())
1256                        coeff = coeff + (thres>>1);
1257                      else
1258                        coeff = coeff - thres + (thres>>1);
1259                    }
1260                  // store coefficient
1261                  if (pcoeff[i] > 0)
1262                    pcoeff[i] = coeff;
1263                  else
1264                    pcoeff[i] = -coeff;
1265                }
1266          }
1267    }
1268}
1269
1270
1271//////////////////////////////////////////////////////
1272// UTILITIES
1273//////////////////////////////////////////////////////
1274
1275
1276#ifdef min
1277#undef min
1278#endif
1279static inline int
1280min(const int x, const int y)
1281{
1282  return (x <= y) ? x : y;
1283}
1284
1285#ifdef max
1286#undef max
1287#endif
1288static inline int
1289max(const int x, const int y)
1290{
1291  return (y <= x) ? x : y;
1292}
1293
1294
1295void 
1296IW44Image::PrimaryHeader::decode(GP<ByteStream> gbs)
1297{
1298  serial = gbs->read8();
1299  slices = gbs->read8();
1300}
1301
1302void 
1303IW44Image::SecondaryHeader::decode(GP<ByteStream> gbs)
1304{
1305  major = gbs->read8();
1306  minor = gbs->read8();
1307}
1308
1309void 
1310IW44Image::TertiaryHeader::decode(GP<ByteStream> gbs, int major, int minor)
1311{
1312  xhi = gbs->read8();
1313  xlo = gbs->read8();
1314  yhi = gbs->read8();
1315  ylo = gbs->read8();
1316  crcbdelay = 0;
1317  if (major== 1 && minor>=2)
1318    crcbdelay = gbs->read8();
1319}
1320
1321
1322
1323//////////////////////////////////////////////////////
1324// CLASS IW44Image
1325//////////////////////////////////////////////////////
1326
1327IW44Image::IW44Image(void)
1328  : db_frac(1.0),
1329    ymap(0), cbmap(0), crmap(0),
1330    cslice(0), cserial(0), cbytes(0)
1331{}
1332
1333IW44Image::~IW44Image()
1334{
1335  delete ymap;
1336  delete cbmap;
1337  delete crmap;
1338}
1339
1340GP<IW44Image>
1341IW44Image::create_decode(const ImageType itype)
1342{
1343  switch(itype)
1344  {
1345  case COLOR:
1346    return new IWPixmap();
1347  case GRAY:
1348    return new IWBitmap();
1349  default:
1350    return 0;
1351  }
1352}
1353
1354int
1355IW44Image::encode_chunk(GP<ByteStream>, const IWEncoderParms &)
1356{
1357  G_THROW( ERR_MSG("IW44Image.codec_open2") );
1358  return 0;
1359}
1360
1361void 
1362IW44Image::encode_iff(IFFByteStream &, int nchunks, const IWEncoderParms *)
1363{
1364  G_THROW( ERR_MSG("IW44Image.codec_open2") );
1365}
1366
1367 
1368void 
1369IWBitmap::close_codec(void)
1370{
1371  delete ycodec;
1372  ycodec = 0;
1373  cslice = cbytes = cserial = 0;
1374}
1375
1376void 
1377IWPixmap::close_codec(void)
1378{
1379  delete ycodec;
1380  delete cbcodec;
1381  delete crcodec;
1382  ycodec = crcodec = cbcodec = 0;
1383  cslice = cbytes = cserial = 0;
1384}
1385
1386int 
1387IW44Image::get_width(void) const
1388{
1389  return (ymap)?(ymap->iw):0;
1390}
1391
1392int 
1393IW44Image::get_height(void) const
1394{
1395  return (ymap)?(ymap->ih):0;
1396}
1397
1398
1399//////////////////////////////////////////////////////
1400// CLASS IWBITMAP
1401//////////////////////////////////////////////////////
1402
1403IWBitmap::IWBitmap(void )
1404: IW44Image(), ycodec(0)
1405{}
1406
1407IWBitmap::~IWBitmap()
1408{
1409  close_codec();
1410}
1411
1412int
1413IWBitmap::get_percent_memory(void) const
1414{
1415  int buckets = 0;
1416  int maximum = 0;
1417  if (ymap) 
1418    {
1419      buckets += ymap->get_bucket_count();
1420      maximum += 64 * ymap->nb;
1421    }
1422  return 100*buckets/ (maximum ? maximum : 1);
1423}
1424
1425unsigned int
1426IWBitmap::get_memory_usage(void) const
1427{
1428  unsigned int usage = sizeof(GBitmap);
1429  if (ymap)
1430    usage += ymap->get_memory_usage();
1431  return usage;
1432}
1433
1434
1435GP<GBitmap> 
1436IWBitmap::get_bitmap(void)
1437{
1438  // Check presence of data
1439  if (ymap == 0)
1440    return 0;
1441  // Perform wavelet reconstruction
1442  int w = ymap->iw;
1443  int h = ymap->ih;
1444  GP<GBitmap> pbm = GBitmap::create(h, w);
1445  ymap->image((signed char*)(*pbm)[0],pbm->rowsize());
1446  // Shift image data
1447  for (int i=0; i<h; i++)
1448    {
1449      unsigned char *urow = (*pbm)[i];
1450      signed char *srow = (signed char*)urow;
1451      for (int j=0; j<w; j++)
1452        urow[j] = (int)(srow[j]) + 128;
1453    }
1454  pbm->set_grays(256);
1455  return pbm;
1456}
1457
1458
1459GP<GBitmap>
1460IWBitmap::get_bitmap(int subsample, const GRect &rect)
1461{
1462  if (ymap == 0)
1463    return 0;
1464  // Allocate bitmap
1465  int w = rect.width();
1466  int h = rect.height();
1467  GP<GBitmap> pbm = GBitmap::create(h,w);
1468  ymap->image(subsample, rect, (signed char*)(*pbm)[0],pbm->rowsize());
1469  // Shift image data
1470  for (int i=0; i<h; i++)
1471    {
1472      unsigned char *urow = (*pbm)[i];
1473      signed char *srow = (signed char*)urow;
1474      for (int j=0; j<w; j++)
1475        urow[j] = (int)(srow[j]) + 128;
1476    }
1477  pbm->set_grays(256);
1478  return pbm;
1479}
1480
1481
1482int
1483IWBitmap::decode_chunk(GP<ByteStream> gbs)
1484{
1485  // Open
1486  if (! ycodec)
1487  {
1488    cslice = cserial = 0;
1489    delete ymap;
1490    ymap = 0;
1491  }
1492  // Read primary header
1493  struct IW44Image::PrimaryHeader primary;
1494  primary.decode(gbs);
1495  if (primary.serial != cserial)
1496    G_THROW( ERR_MSG("IW44Image.wrong_serial") );
1497  int nslices = cslice + primary.slices;
1498  // Read auxilliary headers
1499  if (cserial == 0)
1500    {
1501      struct IW44Image::SecondaryHeader secondary;
1502      secondary.decode(gbs);
1503      if ((secondary.major & 0x7f) != IWCODEC_MAJOR)
1504        G_THROW( ERR_MSG("IW44Image.incompat_codec") );
1505      if (secondary.minor > IWCODEC_MINOR)
1506        G_THROW( ERR_MSG("IW44Image.recent_codec") );
1507      // Read tertiary header
1508      struct IW44Image::TertiaryHeader tertiary;
1509      tertiary.decode(gbs, secondary.major & 0x7f, secondary.minor);
1510      if (! (secondary.major & 0x80))
1511        G_THROW( ERR_MSG("IW44Image.has_color") );
1512      // Create ymap and ycodec
1513      int w = (tertiary.xhi << 8) | tertiary.xlo;
1514      int h = (tertiary.yhi << 8) | tertiary.ylo;
1515      assert(! ymap);
1516      ymap = new Map(w, h);
1517      assert(! ycodec);
1518      ycodec = new Codec::Decode(*ymap);
1519    }
1520  // Read data
1521  assert(ymap);
1522  assert(ycodec);
1523  GP<ZPCodec> gzp=ZPCodec::create(gbs, false, true);
1524  ZPCodec &zp=*gzp;
1525  int flag = 1;
1526  while (flag && cslice<nslices)
1527    {
1528      flag = ycodec->code_slice(zp);
1529      cslice++;
1530    }
1531  // Return
1532  cserial += 1;
1533  return nslices;
1534}
1535
1536void 
1537IWBitmap::parm_dbfrac(float frac)
1538{
1539  if (frac>0 && frac<=1)
1540    db_frac = frac;
1541  else
1542    G_THROW( ERR_MSG("IW44Image.param_range") );
1543}
1544
1545
1546int 
1547IWBitmap::get_serial(void)
1548{
1549  return cserial;
1550}
1551
1552void 
1553IWBitmap::decode_iff(IFFByteStream &iff, int maxchunks)
1554{
1555  if (ycodec)
1556    G_THROW( ERR_MSG("IW44Image.left_open2") );
1557  GUTF8String chkid;
1558  iff.get_chunk(chkid);
1559  if (chkid != "FORM:BM44")
1560    G_THROW( ERR_MSG("IW44Image.corrupt_BM44") );
1561  while (--maxchunks>=0 && iff.get_chunk(chkid))
1562    {
1563      if (chkid == "BM44")
1564        decode_chunk(iff.get_bytestream());
1565      iff.close_chunk();
1566    }
1567  iff.close_chunk();
1568  close_codec();
1569}
1570
1571
1572
1573
1574//////////////////////////////////////////////////////
1575// CLASS IWENCODERPARMS
1576//////////////////////////////////////////////////////
1577
1578
1579IWEncoderParms::IWEncoderParms(void)
1580{
1581  // Zero represent default values
1582  memset((void*)this, 0, sizeof(IWEncoderParms));
1583}
1584
1585
1586
1587
1588
1589//////////////////////////////////////////////////////
1590// CLASS IWPIXMAP
1591//////////////////////////////////////////////////////
1592
1593
1594IWPixmap::IWPixmap(void)
1595: IW44Image(), crcb_delay(10), crcb_half(0), ycodec(0), cbcodec(0), crcodec(0)
1596{}
1597
1598IWPixmap::~IWPixmap()
1599{
1600  close_codec();
1601}
1602
1603int
1604IWPixmap::get_percent_memory(void) const
1605{
1606  int buckets = 0;
1607  int maximum = 0;
1608  if (ymap)
1609    {
1610      buckets += ymap->get_bucket_count();
1611      maximum += 64*ymap->nb;
1612    }
1613  if (cbmap)
1614    {
1615      buckets += cbmap->get_bucket_count();
1616      maximum += 64*cbmap->nb;
1617    }
1618  if (crmap)
1619    {
1620      buckets += crmap->get_bucket_count();
1621      maximum += 64*crmap->nb;
1622    }
1623  return 100*buckets/ (maximum ? maximum : 1);
1624}
1625
1626unsigned int
1627IWPixmap::get_memory_usage(void) const
1628{
1629  unsigned int usage = sizeof(GPixmap);
1630  if (ymap)
1631    usage += ymap->get_memory_usage();
1632  if (cbmap)
1633    usage += cbmap->get_memory_usage();
1634  if (crmap)
1635    usage += crmap->get_memory_usage();
1636  return usage;
1637}
1638
1639
1640GP<GPixmap> 
1641IWPixmap::get_pixmap(void)
1642{
1643  // Check presence of data
1644  if (ymap == 0)
1645    return 0;
1646  // Allocate pixmap
1647  int w = ymap->iw;
1648  int h = ymap->ih;
1649  GP<GPixmap> ppm = GPixmap::create(h, w);
1650  // Perform wavelet reconstruction
1651  signed char *ptr = (signed char*) (*ppm)[0];
1652  int rowsep = ppm->rowsize() * sizeof(GPixel);
1653  int pixsep = sizeof(GPixel);
1654  ymap->image(ptr, rowsep, pixsep);
1655  if (crmap && cbmap && crcb_delay >= 0)
1656  {
1657    cbmap->image(ptr+1, rowsep, pixsep, crcb_half);
1658    crmap->image(ptr+2, rowsep, pixsep, crcb_half);
1659  }
1660  // Convert image data to RGB
1661  if (crmap && cbmap && crcb_delay >= 0)
1662    {
1663      Transform::Decode::YCbCr_to_RGB((*ppm)[0], w, h, ppm->rowsize());
1664    }
1665  else
1666    {
1667      for (int i=0; i<h; i++)
1668        {
1669          GPixel *pixrow = (*ppm)[i];
1670          for (int j=0; j<w; j++, pixrow++)
1671            pixrow->b = pixrow->g = pixrow->r
1672              = 127 - (int)(((signed char*)pixrow)[0]);
1673        }
1674    }
1675  // Return
1676  return ppm;
1677}
1678
1679
1680
1681GP<GPixmap>
1682IWPixmap::get_pixmap(int subsample, const GRect &rect)
1683{
1684  if (ymap == 0)
1685    return 0;
1686  // Allocate
1687  int w = rect.width();
1688  int h = rect.height();
1689  GP<GPixmap> ppm = GPixmap::create(h,w);
1690  // Perform wavelet reconstruction
1691  signed char *ptr = (signed char*) (*ppm)[0];
1692  int rowsep = ppm->rowsize() * sizeof(GPixel);
1693  int pixsep = sizeof(GPixel);
1694  ymap->image(subsample, rect, ptr, rowsep, pixsep);
1695  if (crmap && cbmap && crcb_delay >= 0)
1696  {
1697    cbmap->image(subsample, rect, ptr+1, rowsep, pixsep, crcb_half);
1698    crmap->image(subsample, rect, ptr+2, rowsep, pixsep, crcb_half);
1699  }
1700  // Convert image data to RGB
1701  if (crmap && cbmap && crcb_delay >= 0)
1702    {
1703      Transform::Decode::YCbCr_to_RGB((*ppm)[0], w, h, ppm->rowsize());
1704    }
1705  else
1706    {
1707      for (int i=0; i<h; i++)
1708        {
1709          GPixel *pixrow = (*ppm)[i];
1710          for (int j=0; j<w; j++, pixrow++)
1711            pixrow->b = pixrow->g = pixrow->r
1712              = 127 - (int)(((signed char*)pixrow)[0]);
1713        }
1714    }
1715  // Return
1716  return ppm;
1717}
1718
1719
1720int
1721IWPixmap::decode_chunk(GP<ByteStream> gbs)
1722{
1723  // Open
1724  if (! ycodec)
1725  {
1726      cslice = cserial = 0;
1727      delete ymap;
1728      ymap = 0;
1729  }
1730
1731  // Read primary header
1732  struct IW44Image::PrimaryHeader primary;
1733  primary.decode(gbs);
1734  if (primary.serial != cserial)
1735    G_THROW( ERR_MSG("IW44Image.wrong_serial2") );
1736  int nslices = cslice + primary.slices;
1737  // Read secondary header
1738  if (cserial == 0)
1739    {
1740      struct IW44Image::SecondaryHeader secondary;
1741      secondary.decode(gbs);
1742      if ((secondary.major & 0x7f) != IWCODEC_MAJOR)
1743        G_THROW( ERR_MSG("IW44Image.incompat_codec2") );
1744      if (secondary.minor > IWCODEC_MINOR)
1745        G_THROW( ERR_MSG("IW44Image.recent_codec2") );
1746      // Read tertiary header
1747      struct IW44Image::TertiaryHeader tertiary;
1748      tertiary.decode(gbs, secondary.major & 0x7f, secondary.minor);
1749      // Handle header information
1750      int w = (tertiary.xhi << 8) | tertiary.xlo;
1751      int h = (tertiary.yhi << 8) | tertiary.ylo;
1752      crcb_delay = 0;
1753      crcb_half = 0;
1754      if (secondary.minor>=2)
1755        crcb_delay = tertiary.crcbdelay & 0x7f;
1756      if (secondary.minor>=2)
1757        crcb_half = (tertiary.crcbdelay & 0x80 ? 0 : 1);
1758      if (secondary.major & 0x80)
1759        crcb_delay = -1;
1760      // Create ymap and ycodec   
1761      assert(! ymap);
1762      assert(! ycodec);
1763      ymap = new Map(w, h);
1764      ycodec = new Codec::Decode(*ymap);
1765      if (crcb_delay >= 0)
1766        {
1767          cbmap = new Map(w, h);
1768          crmap = new Map(w, h);
1769          cbcodec = new Codec::Decode(*cbmap);
1770          crcodec = new Codec::Decode(*crmap);
1771        }
1772    }
1773  // Read data
1774  assert(ymap);
1775  assert(ycodec);
1776  GP<ZPCodec> gzp=ZPCodec::create(gbs, false, true);
1777  ZPCodec &zp=*gzp;
1778  int flag = 1;
1779  while (flag && cslice<nslices)
1780    {
1781      flag = ycodec->code_slice(zp);
1782      if (crcodec && cbcodec && crcb_delay<=cslice)
1783        {
1784          flag |= cbcodec->code_slice(zp);
1785          flag |= crcodec->code_slice(zp);
1786        }
1787      cslice++;
1788    }
1789  // Return
1790  cserial += 1;
1791  return nslices;
1792}
1793
1794
1795int 
1796IWPixmap::parm_crcbdelay(const int parm)
1797{
1798  if (parm >= 0)
1799    crcb_delay = parm;
1800  return crcb_delay;
1801}
1802
1803void 
1804IWPixmap::parm_dbfrac(float frac)
1805{
1806  if (frac>0 && frac<=1)
1807    db_frac = frac;
1808  else
1809    G_THROW( ERR_MSG("IW44Image.param_range2") );
1810}
1811
1812int 
1813IWPixmap::get_serial(void)
1814{
1815  return cserial;
1816}
1817
1818
1819void 
1820IWPixmap::decode_iff(IFFByteStream &iff, int maxchunks)
1821{
1822  if (ycodec)
1823    G_THROW( ERR_MSG("IW44Image.left_open4") );
1824  GUTF8String chkid;
1825  iff.get_chunk(chkid);
1826  if (chkid!="FORM:PM44" && chkid!="FORM:BM44")
1827    G_THROW( ERR_MSG("IW44Image.corrupt_BM44_2") );
1828  while (--maxchunks>=0 && iff.get_chunk(chkid))
1829    {
1830      if (chkid=="PM44" || chkid=="BM44")
1831        decode_chunk(iff.get_bytestream());
1832      iff.close_chunk();
1833    }
1834  iff.close_chunk();
1835  close_codec();
1836}
1837
1838//////////////////////////////////////////////////////
1839// NEW FILTERS
1840//////////////////////////////////////////////////////
1841
1842void
1843IW44Image::Transform::filter_begin(int w, int h)
1844{
1845  if (MMXControl::mmxflag < 0) 
1846    MMXControl::enable_mmx();
1847}
1848
1849
1850void
1851IW44Image::Transform::filter_end(void)
1852{
1853#ifdef MMX
1854  if (MMXControl::mmxflag > 0)
1855    MMXemms;
1856#endif
1857}
1858
1859
1860//////////////////////////////////////////////////////
1861// WAVELET TRANSFORM
1862//////////////////////////////////////////////////////
1863
1864
1865//----------------------------------------------------
1866// Function for applying bidimensional IW44 between
1867// scale intervals begin(inclusive) and end(exclusive)
1868
1869void
1870IW44Image::Transform::Decode::backward(short *p, int w, int h, int rowsize, int begin, int end)
1871{ 
1872  // PREPARATION
1873  filter_begin(w,h);
1874  // LOOP ON SCALES
1875  for (int scale=begin>>1; scale>=end; scale>>=1)
1876    {
1877#ifdef IWTRANSFORM_TIMER
1878      int tv,th;
1879      th = tv = GOS::ticks();
1880#endif
1881      filter_bv(p, w, h, rowsize, scale);
1882#ifdef IWTRANSFORM_TIMER
1883      th = GOS::ticks();
1884      tv = th - tv;
1885#endif
1886      filter_bh(p, w, h, rowsize, scale);
1887#ifdef IWTRANSFORM_TIMER
1888      th = GOS::ticks()-th;
1889      DjVuPrintErrorUTF8("back%d\tv=%dms h=%dms\n", scale,tv,th);
1890#endif
1891    }
1892  // TERMINATE
1893  filter_end();
1894}
1895 
1896
1897
1898
1899//////////////////////////////////////////////////////
1900// COLOR TRANSFORM
1901//////////////////////////////////////////////////////
1902
1903/* Converts YCbCr to RGB. */
1904void 
1905IW44Image::Transform::Decode::YCbCr_to_RGB(GPixel *p, int w, int h, int rowsize)
1906{
1907  for (int i=0; i<h; i++,p+=rowsize)
1908    {
1909      GPixel *q = p;
1910      for (int j=0; j<w; j++,q++)
1911        {
1912          signed char y = ((signed char*)q)[0];
1913          signed char b = ((signed char*)q)[1];
1914          signed char r = ((signed char*)q)[2];
1915          // This is the Pigeon transform
1916          int t1 = b >> 2 ; 
1917          int t2 = r + (r >> 1);
1918          int t3 = y + 128 - t1;
1919          int tr = y + 128 + t2;
1920          int tg = t3 - (t2 >> 1);
1921          int tb = t3 + (b << 1);
1922          q->r = max(0,min(255,tr));
1923          q->g = max(0,min(255,tg));
1924          q->b = max(0,min(255,tb));
1925        }
1926    }
1927}
1928
1929
1930#ifdef HAVE_NAMESPACES
1931}
1932# ifndef NOT_USING_DJVU_NAMESPACE
1933using namespace DJVU;
1934# endif
1935#endif
Note: See TracBrowser for help on using the repository browser.