1 ///////////////////////////////////////////////////////////////////////////
2 //
3 // Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
4 // Digital Ltd. LLC
5 //
6 // All rights reserved.
7 //
8 // Redistribution and use in source and binary forms, with or without
9 // modification, are permitted provided that the following conditions are
10 // met:
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
16 // distribution.
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.
20 //
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.
32 //
33 ///////////////////////////////////////////////////////////////////////////
34
35
36
37 //-----------------------------------------------------------------------------
38 //
39 // 16-bit Haar Wavelet encoding and decoding
40 //
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.
44 //
45 //-----------------------------------------------------------------------------
46
47
48 #include <ImfWav.h>
49 #include "ImfNamespace.h"
50
51 OPENEXR_IMF_INTERNAL_NAMESPACE_SOURCE_ENTER
52 namespace {
53
54
55 //
56 // Wavelet basis functions without modulo arithmetic; they produce
57 // the best compression ratios when the wavelet-transformed data are
58 // Huffman-encoded, but the wavelet transform works only for 14-bit
59 // data (untransformed data values must be less than (1 << 14)).
60 //
61
62 inline void
wenc14(unsigned short a,unsigned short b,unsigned short & l,unsigned short & h)63 wenc14 (unsigned short a, unsigned short b,
64 unsigned short &l, unsigned short &h)
65 {
66 short as = a;
67 short bs = b;
68
69 short ms = (as + bs) >> 1;
70 short ds = as - bs;
71
72 l = ms;
73 h = ds;
74 }
75
76
77 inline void
wdec14(unsigned short l,unsigned short h,unsigned short & a,unsigned short & b)78 wdec14 (unsigned short l, unsigned short h,
79 unsigned short &a, unsigned short &b)
80 {
81 short ls = l;
82 short hs = h;
83
84 int hi = hs;
85 int ai = ls + (hi & 1) + (hi >> 1);
86
87 short as = ai;
88 short bs = ai - hi;
89
90 a = as;
91 b = bs;
92 }
93
94
95 //
96 // Wavelet basis functions with modulo arithmetic; they work with full
97 // 16-bit data, but Huffman-encoding the wavelet-transformed data doesn't
98 // compress the data quite as well.
99 //
100
101 const int NBITS = 16;
102 const int A_OFFSET = 1 << (NBITS - 1);
103 const int M_OFFSET = 1 << (NBITS - 1);
104 const int MOD_MASK = (1 << NBITS) - 1;
105
106
107 inline void
wenc16(unsigned short a,unsigned short b,unsigned short & l,unsigned short & h)108 wenc16 (unsigned short a, unsigned short b,
109 unsigned short &l, unsigned short &h)
110 {
111 int ao = (a + A_OFFSET) & MOD_MASK;
112 int m = ((ao + b) >> 1);
113 int d = ao - b;
114
115 if (d < 0)
116 m = (m + M_OFFSET) & MOD_MASK;
117
118 d &= MOD_MASK;
119
120 l = m;
121 h = d;
122 }
123
124
125 inline void
wdec16(unsigned short l,unsigned short h,unsigned short & a,unsigned short & b)126 wdec16 (unsigned short l, unsigned short h,
127 unsigned short &a, unsigned short &b)
128 {
129 int m = l;
130 int d = h;
131 int bb = (m - (d >> 1)) & MOD_MASK;
132 int aa = (d + bb - A_OFFSET) & MOD_MASK;
133 b = bb;
134 a = aa;
135 }
136
137 } // namespace
138
139
140 //
141 // 2D Wavelet encoding:
142 //
143
144 void
wav2Encode(unsigned short * in,int nx,int ox,int ny,int oy,unsigned short mx)145 wav2Encode
146 (unsigned short* in, // io: values are transformed in place
147 int nx, // i : x size
148 int ox, // i : x offset
149 int ny, // i : y size
150 int oy, // i : y offset
151 unsigned short mx) // i : maximum in[x][y] value
152 {
153 bool w14 = (mx < (1 << 14));
154 int n = (nx > ny)? ny: nx;
155 int p = 1; // == 1 << level
156 int p2 = 2; // == 1 << (level+1)
157
158 //
159 // Hierachical loop on smaller dimension n
160 //
161
162 while (p2 <= n)
163 {
164 unsigned short *py = in;
165 unsigned short *ey = in + oy * (ny - p2);
166 int oy1 = oy * p;
167 int oy2 = oy * p2;
168 int ox1 = ox * p;
169 int ox2 = ox * p2;
170 unsigned short i00,i01,i10,i11;
171
172 //
173 // Y loop
174 //
175
176 for (; py <= ey; py += oy2)
177 {
178 unsigned short *px = py;
179 unsigned short *ex = py + ox * (nx - p2);
180
181 //
182 // X loop
183 //
184
185 for (; px <= ex; px += ox2)
186 {
187 unsigned short *p01 = px + ox1;
188 unsigned short *p10 = px + oy1;
189 unsigned short *p11 = p10 + ox1;
190
191 //
192 // 2D wavelet encoding
193 //
194
195 if (w14)
196 {
197 wenc14 (*px, *p01, i00, i01);
198 wenc14 (*p10, *p11, i10, i11);
199 wenc14 (i00, i10, *px, *p10);
200 wenc14 (i01, i11, *p01, *p11);
201 }
202 else
203 {
204 wenc16 (*px, *p01, i00, i01);
205 wenc16 (*p10, *p11, i10, i11);
206 wenc16 (i00, i10, *px, *p10);
207 wenc16 (i01, i11, *p01, *p11);
208 }
209 }
210
211 //
212 // Encode (1D) odd column (still in Y loop)
213 //
214
215 if (nx & p)
216 {
217 unsigned short *p10 = px + oy1;
218
219 if (w14)
220 wenc14 (*px, *p10, i00, *p10);
221 else
222 wenc16 (*px, *p10, i00, *p10);
223
224 *px= i00;
225 }
226 }
227
228 //
229 // Encode (1D) odd line (must loop in X)
230 //
231
232 if (ny & p)
233 {
234 unsigned short *px = py;
235 unsigned short *ex = py + ox * (nx - p2);
236
237 for (; px <= ex; px += ox2)
238 {
239 unsigned short *p01 = px + ox1;
240
241 if (w14)
242 wenc14 (*px, *p01, i00, *p01);
243 else
244 wenc16 (*px, *p01, i00, *p01);
245
246 *px= i00;
247 }
248 }
249
250 //
251 // Next level
252 //
253
254 p = p2;
255 p2 <<= 1;
256 }
257 }
258
259
260 //
261 // 2D Wavelet decoding:
262 //
263
264 void
wav2Decode(unsigned short * in,int nx,int ox,int ny,int oy,unsigned short mx)265 wav2Decode
266 (unsigned short* in, // io: values are transformed in place
267 int nx, // i : x size
268 int ox, // i : x offset
269 int ny, // i : y size
270 int oy, // i : y offset
271 unsigned short mx) // i : maximum in[x][y] value
272 {
273 bool w14 = (mx < (1 << 14));
274 int n = (nx > ny)? ny: nx;
275 int p = 1;
276 int p2;
277
278 //
279 // Search max level
280 //
281
282 while (p <= n)
283 p <<= 1;
284
285 p >>= 1;
286 p2 = p;
287 p >>= 1;
288
289 //
290 // Hierarchical loop on smaller dimension n
291 //
292
293 while (p >= 1)
294 {
295 unsigned short *py = in;
296 unsigned short *ey = in + oy * (ny - p2);
297 int oy1 = oy * p;
298 int oy2 = oy * p2;
299 int ox1 = ox * p;
300 int ox2 = ox * p2;
301 unsigned short i00,i01,i10,i11;
302
303 //
304 // Y loop
305 //
306
307 for (; py <= ey; py += oy2)
308 {
309 unsigned short *px = py;
310 unsigned short *ex = py + ox * (nx - p2);
311
312 //
313 // X loop
314 //
315
316 for (; px <= ex; px += ox2)
317 {
318 unsigned short *p01 = px + ox1;
319 unsigned short *p10 = px + oy1;
320 unsigned short *p11 = p10 + ox1;
321
322 //
323 // 2D wavelet decoding
324 //
325
326 if (w14)
327 {
328 wdec14 (*px, *p10, i00, i10);
329 wdec14 (*p01, *p11, i01, i11);
330 wdec14 (i00, i01, *px, *p01);
331 wdec14 (i10, i11, *p10, *p11);
332 }
333 else
334 {
335 wdec16 (*px, *p10, i00, i10);
336 wdec16 (*p01, *p11, i01, i11);
337 wdec16 (i00, i01, *px, *p01);
338 wdec16 (i10, i11, *p10, *p11);
339 }
340 }
341
342 //
343 // Decode (1D) odd column (still in Y loop)
344 //
345
346 if (nx & p)
347 {
348 unsigned short *p10 = px + oy1;
349
350 if (w14)
351 wdec14 (*px, *p10, i00, *p10);
352 else
353 wdec16 (*px, *p10, i00, *p10);
354
355 *px= i00;
356 }
357 }
358
359 //
360 // Decode (1D) odd line (must loop in X)
361 //
362
363 if (ny & p)
364 {
365 unsigned short *px = py;
366 unsigned short *ex = py + ox * (nx - p2);
367
368 for (; px <= ex; px += ox2)
369 {
370 unsigned short *p01 = px + ox1;
371
372 if (w14)
373 wdec14 (*px, *p01, i00, *p01);
374 else
375 wdec16 (*px, *p01, i00, *p01);
376
377 *px= i00;
378 }
379 }
380
381 //
382 // Next level
383 //
384
385 p2 = p;
386 p >>= 1;
387 }
388 }
389
390
391 OPENEXR_IMF_INTERNAL_NAMESPACE_SOURCE_EXIT
392