source: trunk/libdjvu/IW44Image.cpp @ 206

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

DJVU plugin: djvulibre updated to version 3.5.19

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