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
50namespace Imf {
51namespace {
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
61inline void
62wenc14 (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
76inline void
77wdec14 (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
100const int NBITS = 16;
101const int A_OFFSET =  1 << (NBITS  - 1);
102const int M_OFFSET =  1 << (NBITS  - 1);
103const int MOD_MASK = (1 <<  NBITS) - 1;
104
105
106inline void
107wenc16 (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
124inline void
125wdec16 (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
143void
144wav2Encode
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
263void
264wav2Decode
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