1235368Sgnn///////////////////////////////////////////////////////////////////////////
2235368Sgnn//
3235368Sgnn// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
4235368Sgnn// Digital Ltd. LLC
5235368Sgnn//
6235368Sgnn// All rights reserved.
7235368Sgnn//
8235368Sgnn// Redistribution and use in source and binary forms, with or without
9235368Sgnn// modification, are permitted provided that the following conditions are
10235368Sgnn// met:
11235368Sgnn// *       Redistributions of source code must retain the above copyright
12235368Sgnn// notice, this list of conditions and the following disclaimer.
13235368Sgnn// *       Redistributions in binary form must reproduce the above
14235368Sgnn// copyright notice, this list of conditions and the following disclaimer
15235368Sgnn// in the documentation and/or other materials provided with the
16235368Sgnn// distribution.
17235368Sgnn// *       Neither the name of Industrial Light & Magic nor the names of
18235368Sgnn// its contributors may be used to endorse or promote products derived
19235368Sgnn// from this software without specific prior written permission.
20235368Sgnn//
21235368Sgnn// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
22235368Sgnn// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
23235368Sgnn// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
24235368Sgnn// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
25235368Sgnn// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
26235368Sgnn// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
27235368Sgnn// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
28235368Sgnn// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
29235368Sgnn// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
30235368Sgnn// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
31235368Sgnn// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32235368Sgnn//
33235368Sgnn///////////////////////////////////////////////////////////////////////////
34235368Sgnn
35235368Sgnn
36235368Sgnn
37235368Sgnn//-----------------------------------------------------------------------------
38235368Sgnn//
39235368Sgnn//	16-bit Haar Wavelet encoding and decoding
40235368Sgnn//
41235368Sgnn//	The source code in this file is derived from the encoding
42235368Sgnn//	and decoding routines written by Christian Rouet for his
43235368Sgnn//	PIZ image file format.
44235368Sgnn//
45235368Sgnn//-----------------------------------------------------------------------------
46235368Sgnn
47235368Sgnn
48235368Sgnn#include <ImfWav.h>
49235368Sgnn
50235368Sgnnnamespace Imf {
51235368Sgnnnamespace {
52235368Sgnn
53235368Sgnn
54235368Sgnn//
55235368Sgnn// Wavelet basis functions without modulo arithmetic; they produce
56235368Sgnn// the best compression ratios when the wavelet-transformed data are
57235368Sgnn// Huffman-encoded, but the wavelet transform works only for 14-bit
58235368Sgnn// data (untransformed data values must be less than (1 << 14)).
59235368Sgnn//
60235368Sgnn
61235368Sgnninline void
62235368Sgnnwenc14 (unsigned short  a, unsigned short  b,
63235368Sgnn        unsigned short &l, unsigned short &h)
64235368Sgnn{
65235368Sgnn    short as = a;
66235368Sgnn    short bs = b;
67235368Sgnn
68235368Sgnn    short ms = (as + bs) >> 1;
69235368Sgnn    short ds = as - bs;
70235368Sgnn
71235368Sgnn    l = ms;
72235368Sgnn    h = ds;
73235368Sgnn}
74235368Sgnn
75235368Sgnn
76235368Sgnninline void
77235368Sgnnwdec14 (unsigned short  l, unsigned short  h,
78235368Sgnn        unsigned short &a, unsigned short &b)
79235368Sgnn{
80235368Sgnn    short ls = l;
81235368Sgnn    short hs = h;
82235368Sgnn
83235368Sgnn    int hi = hs;
84235368Sgnn    int ai = ls + (hi & 1) + (hi >> 1);
85235368Sgnn
86235368Sgnn    short as = ai;
87235368Sgnn    short bs = ai - hi;
88235368Sgnn
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