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