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