1
///////////////////////////////////////////////////////////////////////////
3
// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
6
// All rights reserved.
8
// Redistribution and use in source and binary forms, with or without
9
// modification, are permitted provided that the following conditions are
11
// * Redistributions of source code must retain the above copyright
12
// notice, this list of conditions and the following disclaimer.
13
// * Redistributions in binary form must reproduce the above
14
// copyright notice, this list of conditions and the following disclaimer
15
// in the documentation and/or other materials provided with the
17
// * Neither the name of Industrial Light & Magic nor the names of
18
// its contributors may be used to endorse or promote products derived
19
// from this software without specific prior written permission.
21
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
///////////////////////////////////////////////////////////////////////////
37
//-----------------------------------------------------------------------------
39
// 16-bit Haar Wavelet encoding and decoding
41
// The source code in this file is derived from the encoding
42
// and decoding routines written by Christian Rouet for his
43
// PIZ image file format.
45
//-----------------------------------------------------------------------------
55
// Wavelet basis functions without modulo arithmetic; they produce
56
// the best compression ratios when the wavelet-transformed data are
57
// Huffman-encoded, but the wavelet transform works only for 14-bit
58
// data (untransformed data values must be less than (1 << 14)).
62
wenc14 (unsigned short a, unsigned short b,
63
unsigned short &l, unsigned short &h)
68
short ms = (as + bs) >> 1;
77
wdec14 (unsigned short l, unsigned short h,
78
unsigned short &a, unsigned short &b)
84
int ai = ls + (hi & 1) + (hi >> 1);
95
// Wavelet basis functions with modulo arithmetic; they work with full
96
// 16-bit data, but Huffman-encoding the wavelet-transformed data doesn't
97
// compress the data quite as well.
100
const int NBITS = 16;
101
const int A_OFFSET = 1 << (NBITS - 1);
102
const int M_OFFSET = 1 << (NBITS - 1);
103
const int MOD_MASK = (1 << NBITS) - 1;
107
wenc16 (unsigned short a, unsigned short b,
108
unsigned short &l, unsigned short &h)
110
int ao = (a + A_OFFSET) & MOD_MASK;
111
int m = ((ao + b) >> 1);
115
m = (m + M_OFFSET) & MOD_MASK;
125
wdec16 (unsigned short l, unsigned short h,
126
unsigned short &a, unsigned short &b)
130
int bb = (m - (d >> 1)) & MOD_MASK;
131
int aa = (d + bb - A_OFFSET) & MOD_MASK;
140
// 2D Wavelet encoding:
145
(unsigned short* in, // io: values are transformed in place
146
int nx, // i : x size
147
int ox, // i : x offset
148
int ny, // i : y size
149
int oy, // i : y offset
150
unsigned short mx) // i : maximum in[x][y] value
152
bool w14 = (mx < (1 << 14));
153
int n = (nx > ny)? ny: nx;
154
int p = 1; // == 1 << level
155
int p2 = 2; // == 1 << (level+1)
158
// Hierachical loop on smaller dimension n
163
unsigned short *py = in;
164
unsigned short *ey = in + oy * (ny - p2);
169
unsigned short i00,i01,i10,i11;
175
for (; py <= ey; py += oy2)
177
unsigned short *px = py;
178
unsigned short *ex = py + ox * (nx - p2);
184
for (; px <= ex; px += ox2)
186
unsigned short *p01 = px + ox1;
187
unsigned short *p10 = px + oy1;
188
unsigned short *p11 = p10 + ox1;
191
// 2D wavelet encoding
196
wenc14 (*px, *p01, i00, i01);
197
wenc14 (*p10, *p11, i10, i11);
198
wenc14 (i00, i10, *px, *p10);
199
wenc14 (i01, i11, *p01, *p11);
203
wenc16 (*px, *p01, i00, i01);
204
wenc16 (*p10, *p11, i10, i11);
205
wenc16 (i00, i10, *px, *p10);
206
wenc16 (i01, i11, *p01, *p11);
211
// Encode (1D) odd column (still in Y loop)
216
unsigned short *p10 = px + oy1;
219
wenc14 (*px, *p10, i00, *p10);
221
wenc16 (*px, *p10, i00, *p10);
228
// Encode (1D) odd line (must loop in X)
233
unsigned short *px = py;
234
unsigned short *ex = py + ox * (nx - p2);
236
for (; px <= ex; px += ox2)
238
unsigned short *p01 = px + ox1;
241
wenc14 (*px, *p01, i00, *p01);
243
wenc16 (*px, *p01, i00, *p01);
260
// 2D Wavelet decoding:
265
(unsigned short* in, // io: values are transformed in place
266
int nx, // i : x size
267
int ox, // i : x offset
268
int ny, // i : y size
269
int oy, // i : y offset
270
unsigned short mx) // i : maximum in[x][y] value
272
bool w14 = (mx < (1 << 14));
273
int n = (nx > ny)? ny: nx;
289
// Hierarchical loop on smaller dimension n
294
unsigned short *py = in;
295
unsigned short *ey = in + oy * (ny - p2);
300
unsigned short i00,i01,i10,i11;
306
for (; py <= ey; py += oy2)
308
unsigned short *px = py;
309
unsigned short *ex = py + ox * (nx - p2);
315
for (; px <= ex; px += ox2)
317
unsigned short *p01 = px + ox1;
318
unsigned short *p10 = px + oy1;
319
unsigned short *p11 = p10 + ox1;
322
// 2D wavelet decoding
327
wdec14 (*px, *p10, i00, i10);
328
wdec14 (*p01, *p11, i01, i11);
329
wdec14 (i00, i01, *px, *p01);
330
wdec14 (i10, i11, *p10, *p11);
334
wdec16 (*px, *p10, i00, i10);
335
wdec16 (*p01, *p11, i01, i11);
336
wdec16 (i00, i01, *px, *p01);
337
wdec16 (i10, i11, *p10, *p11);
342
// Decode (1D) odd column (still in Y loop)
347
unsigned short *p10 = px + oy1;
350
wdec14 (*px, *p10, i00, *p10);
352
wdec16 (*px, *p10, i00, *p10);
359
// Decode (1D) odd line (must loop in X)
364
unsigned short *px = py;
365
unsigned short *ex = py + ox * (nx - p2);
367
for (; px <= ex; px += ox2)
369
unsigned short *p01 = px + ox1;
372
wdec14 (*px, *p01, i00, *p01);
374
wdec16 (*px, *p01, i00, *p01);
1
///////////////////////////////////////////////////////////////////////////
3
// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
6
// All rights reserved.
8
// Redistribution and use in source and binary forms, with or without
9
// modification, are permitted provided that the following conditions are
11
// * Redistributions of source code must retain the above copyright
12
// notice, this list of conditions and the following disclaimer.
13
// * Redistributions in binary form must reproduce the above
14
// copyright notice, this list of conditions and the following disclaimer
15
// in the documentation and/or other materials provided with the
17
// * Neither the name of Industrial Light & Magic nor the names of
18
// its contributors may be used to endorse or promote products derived
19
// from this software without specific prior written permission.
21
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22
// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23
// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24
// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25
// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26
// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27
// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28
// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29
// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30
// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
///////////////////////////////////////////////////////////////////////////
37
//-----------------------------------------------------------------------------
39
// 16-bit Haar Wavelet encoding and decoding
41
// The source code in this file is derived from the encoding
42
// and decoding routines written by Christian Rouet for his
43
// PIZ image file format.
45
//-----------------------------------------------------------------------------
55
// Wavelet basis functions without modulo arithmetic; they produce
56
// the best compression ratios when the wavelet-transformed data are
57
// Huffman-encoded, but the wavelet transform works only for 14-bit
58
// data (untransformed data values must be less than (1 << 14)).
62
wenc14 (unsigned short a, unsigned short b,
63
unsigned short &l, unsigned short &h)
68
short ms = (as + bs) >> 1;
77
wdec14 (unsigned short l, unsigned short h,
78
unsigned short &a, unsigned short &b)
84
int ai = ls + (hi & 1) + (hi >> 1);
95
// Wavelet basis functions with modulo arithmetic; they work with full
96
// 16-bit data, but Huffman-encoding the wavelet-transformed data doesn't
97
// compress the data quite as well.
100
const int NBITS = 16;
101
const int A_OFFSET = 1 << (NBITS - 1);
102
const int M_OFFSET = 1 << (NBITS - 1);
103
const int MOD_MASK = (1 << NBITS) - 1;
107
wenc16 (unsigned short a, unsigned short b,
108
unsigned short &l, unsigned short &h)
110
int ao = (a + A_OFFSET) & MOD_MASK;
111
int m = ((ao + b) >> 1);
115
m = (m + M_OFFSET) & MOD_MASK;
125
wdec16 (unsigned short l, unsigned short h,
126
unsigned short &a, unsigned short &b)
130
int bb = (m - (d >> 1)) & MOD_MASK;
131
int aa = (d + bb - A_OFFSET) & MOD_MASK;
140
// 2D Wavelet encoding:
145
(unsigned short* in, // io: values are transformed in place
146
int nx, // i : x size
147
int ox, // i : x offset
148
int ny, // i : y size
149
int oy, // i : y offset
150
unsigned short mx) // i : maximum in[x][y] value
152
bool w14 = (mx < (1 << 14));
153
int n = (nx > ny)? ny: nx;
154
int p = 1; // == 1 << level
155
int p2 = 2; // == 1 << (level+1)
158
// Hierachical loop on smaller dimension n
163
unsigned short *py = in;
164
unsigned short *ey = in + oy * (ny - p2);
169
unsigned short i00,i01,i10,i11;
175
for (; py <= ey; py += oy2)
177
unsigned short *px = py;
178
unsigned short *ex = py + ox * (nx - p2);
184
for (; px <= ex; px += ox2)
186
unsigned short *p01 = px + ox1;
187
unsigned short *p10 = px + oy1;
188
unsigned short *p11 = p10 + ox1;
191
// 2D wavelet encoding
196
wenc14 (*px, *p01, i00, i01);
197
wenc14 (*p10, *p11, i10, i11);
198
wenc14 (i00, i10, *px, *p10);
199
wenc14 (i01, i11, *p01, *p11);
203
wenc16 (*px, *p01, i00, i01);
204
wenc16 (*p10, *p11, i10, i11);
205
wenc16 (i00, i10, *px, *p10);
206
wenc16 (i01, i11, *p01, *p11);
211
// Encode (1D) odd column (still in Y loop)
216
unsigned short *p10 = px + oy1;
219
wenc14 (*px, *p10, i00, *p10);
221
wenc16 (*px, *p10, i00, *p10);
228
// Encode (1D) odd line (must loop in X)
233
unsigned short *px = py;
234
unsigned short *ex = py + ox * (nx - p2);
236
for (; px <= ex; px += ox2)
238
unsigned short *p01 = px + ox1;
241
wenc14 (*px, *p01, i00, *p01);
243
wenc16 (*px, *p01, i00, *p01);
260
// 2D Wavelet decoding:
265
(unsigned short* in, // io: values are transformed in place
266
int nx, // i : x size
267
int ox, // i : x offset
268
int ny, // i : y size
269
int oy, // i : y offset
270
unsigned short mx) // i : maximum in[x][y] value
272
bool w14 = (mx < (1 << 14));
273
int n = (nx > ny)? ny: nx;
289
// Hierarchical loop on smaller dimension n
294
unsigned short *py = in;
295
unsigned short *ey = in + oy * (ny - p2);
300
unsigned short i00,i01,i10,i11;
306
for (; py <= ey; py += oy2)
308
unsigned short *px = py;
309
unsigned short *ex = py + ox * (nx - p2);
315
for (; px <= ex; px += ox2)
317
unsigned short *p01 = px + ox1;
318
unsigned short *p10 = px + oy1;
319
unsigned short *p11 = p10 + ox1;
322
// 2D wavelet decoding
327
wdec14 (*px, *p10, i00, i10);
328
wdec14 (*p01, *p11, i01, i11);
329
wdec14 (i00, i01, *px, *p01);
330
wdec14 (i10, i11, *p10, *p11);
334
wdec16 (*px, *p10, i00, i10);
335
wdec16 (*p01, *p11, i01, i11);
336
wdec16 (i00, i01, *px, *p01);
337
wdec16 (i10, i11, *p10, *p11);
342
// Decode (1D) odd column (still in Y loop)
347
unsigned short *p10 = px + oy1;
350
wdec14 (*px, *p10, i00, *p10);
352
wdec16 (*px, *p10, i00, *p10);
359
// Decode (1D) odd line (must loop in X)
364
unsigned short *px = py;
365
unsigned short *ex = py + ox * (nx - p2);
367
for (; px <= ex; px += ox2)
369
unsigned short *p01 = px + ox1;
372
wdec14 (*px, *p01, i00, *p01);
374
wdec16 (*px, *p01, i00, *p01);