~cosme/ubuntu/precise/freeimage/freeimage-3.15.1

« back to all changes in this revision

Viewing changes to Source/OpenEXR/IlmImf/ImfWav.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Cosme Domínguez Díaz
  • Date: 2010-07-20 13:42:15 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20100720134215-xt1454zaedv3b604
Tags: 3.13.1-0ubuntu1
* New upstream release. Closes: (LP: #607800)
 - Updated debian/freeimage-get-orig-source script.
 - Removing no longer necessary debian/patches/* and
   the patch system in debian/rules.
 - Updated debian/rules to work with the new Makefiles.
 - Drop from -O3 to -O2 and use lzma compression saves
   ~10 MB of free space. 
* lintian stuff
 - fixed debhelper-but-no-misc-depends
 - fixed ldconfig-symlink-missing-for-shlib

Show diffs side-by-side

added added

removed removed

Lines of Context:
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
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
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
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
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
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
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
 
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
 
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
 
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
 
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
 
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
 
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
 
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