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
50 namespace Imf {
51 namespace {
52
53
54 //
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)).
59 //
60
61 inline void
wenc14(unsigned short a,unsigned short b,unsigned short & l,unsigned short & h)62 wenc14 (unsigned short a, unsigned short b,
63 unsigned short &l, unsigned short &h)
64 {
65 short as = a;
66 short bs = b;
67
68 short ms = (as + bs) >> 1;
69 short ds = as - bs;
70
71 l = ms;
72 h = ds;
73 }
74
75
76 inline void
wdec14(unsigned short l,unsigned short h,unsigned short & a,unsigned short & b)77 wdec14 (unsigned short l, unsigned short h,
78 unsigned short &a, unsigned short &b)
79 {
80 short ls = l;
81 short hs = h;
82
83 int hi = hs;
84 int ai = ls + (hi & 1) + (hi >> 1);
85
86 short as = ai;
87 short bs = ai - hi;
88
89 a = as;
90 b = bs;
91 }
92
93
94 //
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.
98 //
99
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;
104
105
106 inline void
wenc16(unsigned short a,unsigned short b,unsigned short & l,unsigned short & h)107 wenc16 (unsigned short a, unsigned short b,
108 unsigned short &l, unsigned short &h)
109 {
110 int ao = (a + A_OFFSET) & MOD_MASK;
111 int m = ((ao + b) >> 1);
112 int d = ao - b;
113
114 if (d < 0)
115 m = (m + M_OFFSET) & MOD_MASK;
116
117 d &= MOD_MASK;
118
119 l = m;
120 h = d;
121 }
122
123
124 inline void
wdec16(unsigned short l,unsigned short h,unsigned short & a,unsigned short & b)125 wdec16 (unsigned short l, unsigned short h,
126 unsigned short &a, unsigned short &b)
127 {
128 int m = l;
129 int d = h;
130 int bb = (m - (d >> 1)) & MOD_MASK;
131 int aa = (d + bb - A_OFFSET) & MOD_MASK;
132 b = bb;
133 a = aa;
134 }
135
136 } // namespace
137
138
139 //
140 // 2D Wavelet encoding:
141 //
142
143 void
wav2Encode(unsigned short * in,int nx,int ox,int ny,int oy,unsigned short mx)144 wav2Encode
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
151 {
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)
156
157 //
158 // Hierachical loop on smaller dimension n
159 //
160
161 while (p2 <= n)
162 {
163 unsigned short *py = in;
164 unsigned short *ey = in + oy * (ny - p2);
165 int oy1 = oy * p;
166 int oy2 = oy * p2;
167 int ox1 = ox * p;
168 int ox2 = ox * p2;
169 unsigned short i00,i01,i10,i11;
170
171 //
172 // Y loop
173 //
174
175 for (; py <= ey; py += oy2)
176 {
177 unsigned short *px = py;
178 unsigned short *ex = py + ox * (nx - p2);
179
180 //
181 // X loop
182 //
183
184 for (; px <= ex; px += ox2)
185 {
186 unsigned short *p01 = px + ox1;
187 unsigned short *p10 = px + oy1;
188 unsigned short *p11 = p10 + ox1;
189
190 //
191 // 2D wavelet encoding
192 //
193
194 if (w14)
195 {
196 wenc14 (*px, *p01, i00, i01);
197 wenc14 (*p10, *p11, i10, i11);
198 wenc14 (i00, i10, *px, *p10);
199 wenc14 (i01, i11, *p01, *p11);
200 }
201 else
202 {
203 wenc16 (*px, *p01, i00, i01);
204 wenc16 (*p10, *p11, i10, i11);
205 wenc16 (i00, i10, *px, *p10);
206 wenc16 (i01, i11, *p01, *p11);
207 }
208 }
209
210 //
211 // Encode (1D) odd column (still in Y loop)
212 //
213
214 if (nx & p)
215 {
216 unsigned short *p10 = px + oy1;
217
218 if (w14)
219 wenc14 (*px, *p10, i00, *p10);
220 else
221 wenc16 (*px, *p10, i00, *p10);
222
223 *px= i00;
224 }
225 }
226
227 //
228 // Encode (1D) odd line (must loop in X)
229 //
230
231 if (ny & p)
232 {
233 unsigned short *px = py;
234 unsigned short *ex = py + ox * (nx - p2);
235
236 for (; px <= ex; px += ox2)
237 {
238 unsigned short *p01 = px + ox1;
239
240 if (w14)
241 wenc14 (*px, *p01, i00, *p01);
242 else
243 wenc16 (*px, *p01, i00, *p01);
244
245 *px= i00;
246 }
247 }
248
249 //
250 // Next level
251 //
252
253 p = p2;
254 p2 <<= 1;
255 }
256 }
257
258
259 //
260 // 2D Wavelet decoding:
261 //
262
263 void
wav2Decode(unsigned short * in,int nx,int ox,int ny,int oy,unsigned short mx)264 wav2Decode
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
271 {
272 bool w14 = (mx < (1 << 14));
273 int n = (nx > ny)? ny: nx;
274 int p = 1;
275 int p2;
276
277 //
278 // Search max level
279 //
280
281 while (p <= n)
282 p <<= 1;
283
284 p >>= 1;
285 p2 = p;
286 p >>= 1;
287
288 //
289 // Hierarchical loop on smaller dimension n
290 //
291
292 while (p >= 1)
293 {
294 unsigned short *py = in;
295 unsigned short *ey = in + oy * (ny - p2);
296 int oy1 = oy * p;
297 int oy2 = oy * p2;
298 int ox1 = ox * p;
299 int ox2 = ox * p2;
300 unsigned short i00,i01,i10,i11;
301
302 //
303 // Y loop
304 //
305
306 for (; py <= ey; py += oy2)
307 {
308 unsigned short *px = py;
309 unsigned short *ex = py + ox * (nx - p2);
310
311 //
312 // X loop
313 //
314
315 for (; px <= ex; px += ox2)
316 {
317 unsigned short *p01 = px + ox1;
318 unsigned short *p10 = px + oy1;
319 unsigned short *p11 = p10 + ox1;
320
321 //
322 // 2D wavelet decoding
323 //
324
325 if (w14)
326 {
327 wdec14 (*px, *p10, i00, i10);
328 wdec14 (*p01, *p11, i01, i11);
329 wdec14 (i00, i01, *px, *p01);
330 wdec14 (i10, i11, *p10, *p11);
331 }
332 else
333 {
334 wdec16 (*px, *p10, i00, i10);
335 wdec16 (*p01, *p11, i01, i11);
336 wdec16 (i00, i01, *px, *p01);
337 wdec16 (i10, i11, *p10, *p11);
338 }
339 }
340
341 //
342 // Decode (1D) odd column (still in Y loop)
343 //
344
345 if (nx & p)
346 {
347 unsigned short *p10 = px + oy1;
348
349 if (w14)
350 wdec14 (*px, *p10, i00, *p10);
351 else
352 wdec16 (*px, *p10, i00, *p10);
353
354 *px= i00;
355 }
356 }
357
358 //
359 // Decode (1D) odd line (must loop in X)
360 //
361
362 if (ny & p)
363 {
364 unsigned short *px = py;
365 unsigned short *ex = py + ox * (nx - p2);
366
367 for (; px <= ex; px += ox2)
368 {
369 unsigned short *p01 = px + ox1;
370
371 if (w14)
372 wdec14 (*px, *p01, i00, *p01);
373 else
374 wdec16 (*px, *p01, i00, *p01);
375
376 *px= i00;
377 }
378 }
379
380 //
381 // Next level
382 //
383
384 p2 = p;
385 p >>= 1;
386 }
387 }
388
389
390 } // namespace Imf
391