← Back to branch summary

~ubuntu-branches/ubuntu/saucy/emscripten/saucy-proposed

« back to all changes in this revision

Viewing changes to tests/openjpeg/libopenjpeg/dwt.c

  • Committer: Package Import Robot
  • Author(s): Sylvestre Ledru
  • Date: 2013-05-02 13:11:51 UTC
  • Revision ID: package-import@ubuntu.com-20130502131151-q8dvteqr1ef2x7xz
Tags: upstream-1.4.1~20130504~adb56cb
Import upstream version 1.4.1~20130504~adb56cb

Show diffs side-by-side

added added

removed removed

Lines of Context:
tests/openjpeg/libopenjpeg/dwt.c
 
1
/*
 
2
 * Copyright (c) 2002-2007, Communications and Remote Sensing Laboratory, Universite catholique de Louvain (UCL), Belgium
 
3
 * Copyright (c) 2002-2007, Professor Benoit Macq
 
4
 * Copyright (c) 2001-2003, David Janssens
 
5
 * Copyright (c) 2002-2003, Yannick Verschueren
 
6
 * Copyright (c) 2003-2007, Francois-Olivier Devaux and Antonin Descampe
 
7
 * Copyright (c) 2005, Herve Drolon, FreeImage Team
 
8
 * Copyright (c) 2007, Jonathan Ballard <dzonatas@dzonux.net>
 
9
 * Copyright (c) 2007, Callum Lerwick <seg@haxxed.com>
 
10
 * All rights reserved.
 
11
 *
 
12
 * Redistribution and use in source and binary forms, with or without
 
13
 * modification, are permitted provided that the following conditions
 
14
 * are met:
 
15
 * 1. Redistributions of source code must retain the above copyright
 
16
 *    notice, this list of conditions and the following disclaimer.
 
17
 * 2. Redistributions in binary form must reproduce the above copyright
 
18
 *    notice, this list of conditions and the following disclaimer in the
 
19
 *    documentation and/or other materials provided with the distribution.
 
20
 *
 
21
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
 
22
 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 
23
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 
24
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
 
25
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 
26
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 
27
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 
28
 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
 
29
 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
 
30
 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
 
31
 * POSSIBILITY OF SUCH DAMAGE.
 
32
 */
 
33
 
 
34
#ifdef __SSE__
 
35
#include <xmmintrin.h>
 
36
#endif
 
37
 
 
38
#include "opj_includes.h"
 
39
 
 
40
/** @defgroup DWT DWT - Implementation of a discrete wavelet transform */
 
41
/*@{*/
 
42
 
 
43
#define WS(i) v->mem[(i)*2]
 
44
#define WD(i) v->mem[(1+(i)*2)]
 
45
 
 
46
/** @name Local data structures */
 
47
/*@{*/
 
48
 
 
49
typedef struct dwt_local {
 
50
        int* mem;
 
51
        int dn;
 
52
        int sn;
 
53
        int cas;
 
54
} dwt_t;
 
55
 
 
56
typedef union {
 
57
        float   f[4];
 
58
} v4;
 
59
 
 
60
typedef struct v4dwt_local {
 
61
        v4*     wavelet ;
 
62
        int             dn ;
 
63
        int             sn ;
 
64
        int             cas ;
 
65
} v4dwt_t ;
 
66
 
 
67
static const float dwt_alpha =  1.586134342f; //  12994
 
68
static const float dwt_beta  =  0.052980118f; //    434
 
69
static const float dwt_gamma = -0.882911075f; //  -7233
 
70
static const float dwt_delta = -0.443506852f; //  -3633
 
71
 
 
72
static const float K      = 1.230174105f; //  10078
 
73
/* FIXME: What is this constant? */
 
74
static const float c13318 = 1.625732422f;
 
75
 
 
76
/*@}*/
 
77
 
 
78
/**
 
79
Virtual function type for wavelet transform in 1-D 
 
80
*/
 
81
typedef void (*DWT1DFN)(dwt_t* v);
 
82
 
 
83
/** @name Local static functions */
 
84
/*@{*/
 
85
 
 
86
/**
 
87
Forward lazy transform (horizontal)
 
88
*/
 
89
static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas);
 
90
/**
 
91
Forward lazy transform (vertical)
 
92
*/
 
93
static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas);
 
94
/**
 
95
Inverse lazy transform (horizontal)
 
96
*/
 
97
static void dwt_interleave_h(dwt_t* h, int *a);
 
98
/**
 
99
Inverse lazy transform (vertical)
 
100
*/
 
101
static void dwt_interleave_v(dwt_t* v, int *a, int x);
 
102
/**
 
103
Forward 5-3 wavelet transform in 1-D
 
104
*/
 
105
static void dwt_encode_1(int *a, int dn, int sn, int cas);
 
106
/**
 
107
Inverse 5-3 wavelet transform in 1-D
 
108
*/
 
109
static void dwt_decode_1(dwt_t *v);
 
110
/**
 
111
Forward 9-7 wavelet transform in 1-D
 
112
*/
 
113
static void dwt_encode_1_real(int *a, int dn, int sn, int cas);
 
114
/**
 
115
Explicit calculation of the Quantization Stepsizes 
 
116
*/
 
117
static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno_stepsize);
 
118
/**
 
119
Inverse wavelet transform in 2-D.
 
120
*/
 
121
static void dwt_decode_tile(opj_tcd_tilecomp_t* tilec, int i, DWT1DFN fn);
 
122
 
 
123
/*@}*/
 
124
 
 
125
/*@}*/
 
126
 
 
127
#define S(i) a[(i)*2]
 
128
#define D(i) a[(1+(i)*2)]
 
129
#define S_(i) ((i)<0?S(0):((i)>=sn?S(sn-1):S(i)))
 
130
#define D_(i) ((i)<0?D(0):((i)>=dn?D(dn-1):D(i)))
 
131
/* new */
 
132
#define SS_(i) ((i)<0?S(0):((i)>=dn?S(dn-1):S(i)))
 
133
#define DD_(i) ((i)<0?D(0):((i)>=sn?D(sn-1):D(i)))
 
134
 
 
135
/* <summary>                                                              */
 
136
/* This table contains the norms of the 5-3 wavelets for different bands. */
 
137
/* </summary>                                                             */
 
138
static const double dwt_norms[4][10] = {
 
139
        {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 
140
        {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 
141
        {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 
142
        {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
 
143
};
 
144
 
 
145
/* <summary>                                                              */
 
146
/* This table contains the norms of the 9-7 wavelets for different bands. */
 
147
/* </summary>                                                             */
 
148
static const double dwt_norms_real[4][10] = {
 
149
        {1.000, 1.965, 4.177, 8.403, 16.90, 33.84, 67.69, 135.3, 270.6, 540.9},
 
150
        {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
 
151
        {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
 
152
        {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2}
 
153
};
 
154
 
 
155
/* 
 
156
==========================================================
 
157
   local functions
 
158
==========================================================
 
159
*/
 
160
 
 
161
/* <summary>                                     */
 
162
/* Forward lazy transform (horizontal).  */
 
163
/* </summary>                            */ 
 
164
static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas) {
 
165
        int i;
 
166
    for (i=0; i<sn; i++) b[i]=a[2*i+cas];
 
167
    for (i=0; i<dn; i++) b[sn+i]=a[(2*i+1-cas)];
 
168
}
 
169
 
 
170
/* <summary>                             */  
 
171
/* Forward lazy transform (vertical).    */
 
172
/* </summary>                            */ 
 
173
static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
 
174
    int i;
 
175
    for (i=0; i<sn; i++) b[i*x]=a[2*i+cas];
 
176
    for (i=0; i<dn; i++) b[(sn+i)*x]=a[(2*i+1-cas)];
 
177
}
 
178
 
 
179
/* <summary>                             */
 
180
/* Inverse lazy transform (horizontal).  */
 
181
/* </summary>                            */
 
182
static void dwt_interleave_h(dwt_t* h, int *a) {
 
183
    int *ai = a;
 
184
    int *bi = h->mem + h->cas;
 
185
    int  i      = h->sn;
 
186
    while( i-- ) {
 
187
      *bi = *(ai++);
 
188
          bi += 2;
 
189
    }
 
190
    ai  = a + h->sn;
 
191
    bi  = h->mem + 1 - h->cas;
 
192
    i   = h->dn ;
 
193
    while( i-- ) {
 
194
      *bi = *(ai++);
 
195
          bi += 2;
 
196
    }
 
197
}
 
198
 
 
199
/* <summary>                             */  
 
200
/* Inverse lazy transform (vertical).    */
 
201
/* </summary>                            */ 
 
202
static void dwt_interleave_v(dwt_t* v, int *a, int x) {
 
203
    int *ai = a;
 
204
    int *bi = v->mem + v->cas;
 
205
    int  i = v->sn;
 
206
    while( i-- ) {
 
207
      *bi = *ai;
 
208
          bi += 2;
 
209
          ai += x;
 
210
    }
 
211
    ai = a + (v->sn * x);
 
212
    bi = v->mem + 1 - v->cas;
 
213
    i = v->dn ;
 
214
    while( i-- ) {
 
215
      *bi = *ai;
 
216
          bi += 2;  
 
217
          ai += x;
 
218
    }
 
219
}
 
220
 
 
221
 
 
222
/* <summary>                            */
 
223
/* Forward 5-3 wavelet transform in 1-D. */
 
224
/* </summary>                           */
 
225
static void dwt_encode_1(int *a, int dn, int sn, int cas) {
 
226
        int i;
 
227
        
 
228
        if (!cas) {
 
229
                if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 
230
                        for (i = 0; i < dn; i++) D(i) -= (S_(i) + S_(i + 1)) >> 1;
 
231
                        for (i = 0; i < sn; i++) S(i) += (D_(i - 1) + D_(i) + 2) >> 2;
 
232
                }
 
233
        } else {
 
234
                if (!sn && dn == 1)                 /* NEW :  CASE ONE ELEMENT */
 
235
                        S(0) *= 2;
 
236
                else {
 
237
                        for (i = 0; i < dn; i++) S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
 
238
                        for (i = 0; i < sn; i++) D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
 
239
                }
 
240
        }
 
241
}
 
242
 
 
243
/* <summary>                            */
 
244
/* Inverse 5-3 wavelet transform in 1-D. */
 
245
/* </summary>                           */ 
 
246
static void dwt_decode_1_(int *a, int dn, int sn, int cas) {
 
247
        int i;
 
248
        
 
249
        if (!cas) {
 
250
                if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 
251
                        for (i = 0; i < sn; i++) S(i) -= (D_(i - 1) + D_(i) + 2) >> 2;
 
252
                        for (i = 0; i < dn; i++) D(i) += (S_(i) + S_(i + 1)) >> 1;
 
253
                }
 
254
        } else {
 
255
                if (!sn  && dn == 1)          /* NEW :  CASE ONE ELEMENT */
 
256
                        S(0) /= 2;
 
257
                else {
 
258
                        for (i = 0; i < sn; i++) D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2;
 
259
                        for (i = 0; i < dn; i++) S(i) += (DD_(i) + DD_(i - 1)) >> 1;
 
260
                }
 
261
        }
 
262
}
 
263
 
 
264
/* <summary>                            */
 
265
/* Inverse 5-3 wavelet transform in 1-D. */
 
266
/* </summary>                           */ 
 
267
static void dwt_decode_1(dwt_t *v) {
 
268
        dwt_decode_1_(v->mem, v->dn, v->sn, v->cas);
 
269
}
 
270
 
 
271
/* <summary>                             */
 
272
/* Forward 9-7 wavelet transform in 1-D. */
 
273
/* </summary>                            */
 
274
static void dwt_encode_1_real(int *a, int dn, int sn, int cas) {
 
275
        int i;
 
276
        if (!cas) {
 
277
                if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 
278
                        for (i = 0; i < dn; i++)
 
279
                                D(i) -= fix_mul(S_(i) + S_(i + 1), 12993);
 
280
                        for (i = 0; i < sn; i++)
 
281
                                S(i) -= fix_mul(D_(i - 1) + D_(i), 434);
 
282
                        for (i = 0; i < dn; i++)
 
283
                                D(i) += fix_mul(S_(i) + S_(i + 1), 7233);
 
284
                        for (i = 0; i < sn; i++)
 
285
                                S(i) += fix_mul(D_(i - 1) + D_(i), 3633);
 
286
                        for (i = 0; i < dn; i++)
 
287
                                D(i) = fix_mul(D(i), 5038);     /*5038 */
 
288
                        for (i = 0; i < sn; i++)
 
289
                                S(i) = fix_mul(S(i), 6659);     /*6660 */
 
290
                }
 
291
        } else {
 
292
                if ((sn > 0) || (dn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 
293
                        for (i = 0; i < dn; i++)
 
294
                                S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993);
 
295
                        for (i = 0; i < sn; i++)
 
296
                                D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434);
 
297
                        for (i = 0; i < dn; i++)
 
298
                                S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233);
 
299
                        for (i = 0; i < sn; i++)
 
300
                                D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633);
 
301
                        for (i = 0; i < dn; i++)
 
302
                                S(i) = fix_mul(S(i), 5038);     /*5038 */
 
303
                        for (i = 0; i < sn; i++)
 
304
                                D(i) = fix_mul(D(i), 6659);     /*6660 */
 
305
                }
 
306
        }
 
307
}
 
308
 
 
309
static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno_stepsize) {
 
310
        int p, n;
 
311
        p = int_floorlog2(stepsize) - 13;
 
312
        n = 11 - int_floorlog2(stepsize);
 
313
        bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
 
314
        bandno_stepsize->expn = numbps - p;
 
315
}
 
316
 
 
317
/* 
 
318
==========================================================
 
319
   DWT interface
 
320
==========================================================
 
321
*/
 
322
 
 
323
/* <summary>                            */
 
324
/* Forward 5-3 wavelet transform in 2-D. */
 
325
/* </summary>                           */
 
326
void dwt_encode(opj_tcd_tilecomp_t * tilec) {
 
327
        int i, j, k;
 
328
        int *a = NULL;
 
329
        int *aj = NULL;
 
330
        int *bj = NULL;
 
331
        int w, l;
 
332
        
 
333
        w = tilec->x1-tilec->x0;
 
334
        l = tilec->numresolutions-1;
 
335
        a = tilec->data;
 
336
        
 
337
        for (i = 0; i < l; i++) {
 
338
                int rw;                 /* width of the resolution level computed                                                           */
 
339
                int rh;                 /* height of the resolution level computed                                                          */
 
340
                int rw1;                /* width of the resolution level once lower than computed one                                       */
 
341
                int rh1;                /* height of the resolution level once lower than computed one                                      */
 
342
                int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 
343
                int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 
344
                int dn, sn;
 
345
                
 
346
                rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 
347
                rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 
348
                rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 
349
                rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 
350
                
 
351
                cas_row = tilec->resolutions[l - i].x0 % 2;
 
352
                cas_col = tilec->resolutions[l - i].y0 % 2;
 
353
        
 
354
                sn = rh1;
 
355
                dn = rh - rh1;
 
356
                bj = (int*)opj_malloc(rh * sizeof(int));
 
357
                for (j = 0; j < rw; j++) {
 
358
                        aj = a + j;
 
359
                        for (k = 0; k < rh; k++)  bj[k] = aj[k*w];
 
360
                        dwt_encode_1(bj, dn, sn, cas_col);
 
361
                        dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 
362
                }
 
363
                opj_free(bj);
 
364
                
 
365
                sn = rw1;
 
366
                dn = rw - rw1;
 
367
                bj = (int*)opj_malloc(rw * sizeof(int));
 
368
                for (j = 0; j < rh; j++) {
 
369
                        aj = a + j * w;
 
370
                        for (k = 0; k < rw; k++)  bj[k] = aj[k];
 
371
                        dwt_encode_1(bj, dn, sn, cas_row);
 
372
                        dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 
373
                }
 
374
                opj_free(bj);
 
375
        }
 
376
}
 
377
 
 
378
 
 
379
/* <summary>                            */
 
380
/* Inverse 5-3 wavelet transform in 2-D. */
 
381
/* </summary>                           */
 
382
void dwt_decode(opj_tcd_tilecomp_t* tilec, int numres) {
 
383
        dwt_decode_tile(tilec, numres, &dwt_decode_1);
 
384
}
 
385
 
 
386
 
 
387
/* <summary>                          */
 
388
/* Get gain of 5-3 wavelet transform. */
 
389
/* </summary>                         */
 
390
int dwt_getgain(int orient) {
 
391
        if (orient == 0)
 
392
                return 0;
 
393
        if (orient == 1 || orient == 2)
 
394
                return 1;
 
395
        return 2;
 
396
}
 
397
 
 
398
/* <summary>                */
 
399
/* Get norm of 5-3 wavelet. */
 
400
/* </summary>               */
 
401
double dwt_getnorm(int level, int orient) {
 
402
        return dwt_norms[orient][level];
 
403
}
 
404
 
 
405
/* <summary>                             */
 
406
/* Forward 9-7 wavelet transform in 2-D. */
 
407
/* </summary>                            */
 
408
 
 
409
void dwt_encode_real(opj_tcd_tilecomp_t * tilec) {
 
410
        int i, j, k;
 
411
        int *a = NULL;
 
412
        int *aj = NULL;
 
413
        int *bj = NULL;
 
414
        int w, l;
 
415
        
 
416
        w = tilec->x1-tilec->x0;
 
417
        l = tilec->numresolutions-1;
 
418
        a = tilec->data;
 
419
        
 
420
        for (i = 0; i < l; i++) {
 
421
                int rw;                 /* width of the resolution level computed                                                     */
 
422
                int rh;                 /* height of the resolution level computed                                                    */
 
423
                int rw1;                /* width of the resolution level once lower than computed one                                 */
 
424
                int rh1;                /* height of the resolution level once lower than computed one                                */
 
425
                int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 
426
                int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 
427
                int dn, sn;
 
428
                
 
429
                rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 
430
                rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 
431
                rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 
432
                rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 
433
                
 
434
                cas_row = tilec->resolutions[l - i].x0 % 2;
 
435
                cas_col = tilec->resolutions[l - i].y0 % 2;
 
436
                
 
437
                sn = rh1;
 
438
                dn = rh - rh1;
 
439
                bj = (int*)opj_malloc(rh * sizeof(int));
 
440
                for (j = 0; j < rw; j++) {
 
441
                        aj = a + j;
 
442
                        for (k = 0; k < rh; k++)  bj[k] = aj[k*w];
 
443
                        dwt_encode_1_real(bj, dn, sn, cas_col);
 
444
                        dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 
445
                }
 
446
                opj_free(bj);
 
447
                
 
448
                sn = rw1;
 
449
                dn = rw - rw1;
 
450
                bj = (int*)opj_malloc(rw * sizeof(int));
 
451
                for (j = 0; j < rh; j++) {
 
452
                        aj = a + j * w;
 
453
                        for (k = 0; k < rw; k++)  bj[k] = aj[k];
 
454
                        dwt_encode_1_real(bj, dn, sn, cas_row);
 
455
                        dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 
456
                }
 
457
                opj_free(bj);
 
458
        }
 
459
}
 
460
 
 
461
 
 
462
/* <summary>                          */
 
463
/* Get gain of 9-7 wavelet transform. */
 
464
/* </summary>                         */
 
465
int dwt_getgain_real(int orient) {
 
466
        (void)orient;
 
467
        return 0;
 
468
}
 
469
 
 
470
/* <summary>                */
 
471
/* Get norm of 9-7 wavelet. */
 
472
/* </summary>               */
 
473
double dwt_getnorm_real(int level, int orient) {
 
474
        return dwt_norms_real[orient][level];
 
475
}
 
476
 
 
477
void dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, int prec) {
 
478
        int numbands, bandno;
 
479
        numbands = 3 * tccp->numresolutions - 2;
 
480
        for (bandno = 0; bandno < numbands; bandno++) {
 
481
                double stepsize;
 
482
                int resno, level, orient, gain;
 
483
 
 
484
                resno = (bandno == 0) ? 0 : ((bandno - 1) / 3 + 1);
 
485
                orient = (bandno == 0) ? 0 : ((bandno - 1) % 3 + 1);
 
486
                level = tccp->numresolutions - 1 - resno;
 
487
                gain = (tccp->qmfbid == 0) ? 0 : ((orient == 0) ? 0 : (((orient == 1) || (orient == 2)) ? 1 : 2));
 
488
                if (tccp->qntsty == J2K_CCP_QNTSTY_NOQNT) {
 
489
                        stepsize = 1.0;
 
490
                } else {
 
491
                        double norm = dwt_norms_real[orient][level];
 
492
                        stepsize = (1 << (gain)) / norm;
 
493
                }
 
494
                dwt_encode_stepsize((int) floor(stepsize * 8192.0), prec + gain, &tccp->stepsizes[bandno]);
 
495
        }
 
496
}
 
497
 
 
498
 
 
499
/* <summary>                             */
 
500
/* Determine maximum computed resolution level for inverse wavelet transform */
 
501
/* </summary>                            */
 
502
static int dwt_decode_max_resolution(opj_tcd_resolution_t* restrict r, int i) {
 
503
        int mr  = 1;
 
504
        int w;
 
505
        while( --i ) {
 
506
                r++;
 
507
                if( mr < ( w = r->x1 - r->x0 ) )
 
508
                        mr = w ;
 
509
                if( mr < ( w = r->y1 - r->y0 ) )
 
510
                        mr = w ;
 
511
        }
 
512
        return mr ;
 
513
}
 
514
 
 
515
 
 
516
/* <summary>                            */
 
517
/* Inverse wavelet transform in 2-D.     */
 
518
/* </summary>                           */
 
519
static void dwt_decode_tile(opj_tcd_tilecomp_t* tilec, int numres, DWT1DFN dwt_1D) {
 
520
        dwt_t h;
 
521
        dwt_t v;
 
522
 
 
523
        opj_tcd_resolution_t* tr = tilec->resolutions;
 
524
 
 
525
        int rw = tr->x1 - tr->x0;       /* width of the resolution level computed */
 
526
        int rh = tr->y1 - tr->y0;       /* height of the resolution level computed */
 
527
 
 
528
        int w = tilec->x1 - tilec->x0;
 
529
 
 
530
        h.mem = opj_aligned_malloc(dwt_decode_max_resolution(tr, numres) * sizeof(int));
 
531
        v.mem = h.mem;
 
532
 
 
533
        while( --numres) {
 
534
                int * restrict tiledp = tilec->data;
 
535
                int j;
 
536
 
 
537
                ++tr;
 
538
                h.sn = rw;
 
539
                v.sn = rh;
 
540
 
 
541
                rw = tr->x1 - tr->x0;
 
542
                rh = tr->y1 - tr->y0;
 
543
 
 
544
                h.dn = rw - h.sn;
 
545
                h.cas = tr->x0 % 2;
 
546
 
 
547
                for(j = 0; j < rh; ++j) {
 
548
                        dwt_interleave_h(&h, &tiledp[j*w]);
 
549
                        (dwt_1D)(&h);
 
550
                        memcpy(&tiledp[j*w], h.mem, rw * sizeof(int));
 
551
                }
 
552
 
 
553
                v.dn = rh - v.sn;
 
554
                v.cas = tr->y0 % 2;
 
555
 
 
556
                for(j = 0; j < rw; ++j){
 
557
                        int k;
 
558
                        dwt_interleave_v(&v, &tiledp[j], w);
 
559
                        (dwt_1D)(&v);
 
560
                        for(k = 0; k < rh; ++k) {
 
561
                                tiledp[k * w + j] = v.mem[k];
 
562
                        }
 
563
                }
 
564
        }
 
565
        opj_aligned_free(h.mem);
 
566
}
 
567
 
 
568
static void v4dwt_interleave_h(v4dwt_t* restrict w, float* restrict a, int x, int size){
 
569
        float* restrict bi = (float*) (w->wavelet + w->cas);
 
570
        int count = w->sn;
 
571
        int i, k;
 
572
        for(k = 0; k < 2; ++k){
 
573
                if (count + 3 * x < size && ((int) a & 0x0f) == 0 && ((int) bi & 0x0f) == 0 && (x & 0x0f) == 0) {
 
574
                        /* Fast code path */
 
575
                        for(i = 0; i < count; ++i){
 
576
                                int j = i;
 
577
                                bi[i*8    ] = a[j];
 
578
                                j += x;
 
579
                                bi[i*8 + 1] = a[j];
 
580
                                j += x;
 
581
                                bi[i*8 + 2] = a[j];
 
582
                                j += x;
 
583
                                bi[i*8 + 3] = a[j];
 
584
                        }
 
585
                } else {
 
586
                        /* Slow code path */
 
587
                for(i = 0; i < count; ++i){
 
588
                        int j = i;
 
589
                        bi[i*8    ] = a[j];
 
590
                        j += x;
 
591
                        if(j > size) continue;
 
592
                        bi[i*8 + 1] = a[j];
 
593
                        j += x;
 
594
                        if(j > size) continue;
 
595
                        bi[i*8 + 2] = a[j];
 
596
                        j += x;
 
597
                        if(j > size) continue;
 
598
                        bi[i*8 + 3] = a[j];
 
599
                }
 
600
                }
 
601
                bi = (float*) (w->wavelet + 1 - w->cas);
 
602
                a += w->sn;
 
603
                size -= w->sn;
 
604
                count = w->dn;
 
605
        }
 
606
}
 
607
 
 
608
static void v4dwt_interleave_v(v4dwt_t* restrict v , float* restrict a , int x){
 
609
        v4* restrict bi = v->wavelet + v->cas;
 
610
        int i;
 
611
        for(i = 0; i < v->sn; ++i){
 
612
                memcpy(&bi[i*2], &a[i*x], 4 * sizeof(float));
 
613
        }
 
614
        a += v->sn * x;
 
615
        bi = v->wavelet + 1 - v->cas;
 
616
        for(i = 0; i < v->dn; ++i){
 
617
                memcpy(&bi[i*2], &a[i*x], 4 * sizeof(float));
 
618
        }
 
619
}
 
620
 
 
621
#ifdef __SSE__
 
622
 
 
623
static void v4dwt_decode_step1_sse(v4* w, int count, const __m128 c){
 
624
        __m128* restrict vw = (__m128*) w;
 
625
        int i;
 
626
        /* 4x unrolled loop */
 
627
        for(i = 0; i < count >> 2; ++i){
 
628
                *vw = _mm_mul_ps(*vw, c);
 
629
                vw += 2;
 
630
                *vw = _mm_mul_ps(*vw, c);
 
631
                vw += 2;
 
632
                *vw = _mm_mul_ps(*vw, c);
 
633
                vw += 2;
 
634
                *vw = _mm_mul_ps(*vw, c);
 
635
                vw += 2;
 
636
        }
 
637
        count &= 3;
 
638
        for(i = 0; i < count; ++i){
 
639
                *vw = _mm_mul_ps(*vw, c);
 
640
                vw += 2;
 
641
        }
 
642
}
 
643
 
 
644
static void v4dwt_decode_step2_sse(v4* l, v4* w, int k, int m, __m128 c){
 
645
        __m128* restrict vl = (__m128*) l;
 
646
        __m128* restrict vw = (__m128*) w;
 
647
        int i;
 
648
        __m128 tmp1, tmp2, tmp3;
 
649
        tmp1 = vl[0];
 
650
        for(i = 0; i < m; ++i){
 
651
                tmp2 = vw[-1];
 
652
                tmp3 = vw[ 0];
 
653
                vw[-1] = _mm_add_ps(tmp2, _mm_mul_ps(_mm_add_ps(tmp1, tmp3), c));
 
654
                tmp1 = tmp3;
 
655
                vw += 2;
 
656
        }
 
657
        vl = vw - 2;
 
658
        if(m >= k){
 
659
                return;
 
660
        }
 
661
        c = _mm_add_ps(c, c);
 
662
        c = _mm_mul_ps(c, vl[0]);
 
663
        for(; m < k; ++m){
 
664
                __m128 tmp = vw[-1];
 
665
                vw[-1] = _mm_add_ps(tmp, c);
 
666
                vw += 2;
 
667
        }
 
668
}
 
669
 
 
670
#else
 
671
 
 
672
static void v4dwt_decode_step1(v4* w, int count, const float c){
 
673
        float* restrict fw = (float*) w;
 
674
        int i;
 
675
        for(i = 0; i < count; ++i){
 
676
                float tmp1 = fw[i*8    ];
 
677
                float tmp2 = fw[i*8 + 1];
 
678
                float tmp3 = fw[i*8 + 2];
 
679
                float tmp4 = fw[i*8 + 3];
 
680
                fw[i*8    ] = tmp1 * c;
 
681
                fw[i*8 + 1] = tmp2 * c;
 
682
                fw[i*8 + 2] = tmp3 * c;
 
683
                fw[i*8 + 3] = tmp4 * c;
 
684
        }
 
685
}
 
686
 
 
687
static void v4dwt_decode_step2(v4* l, v4* w, int k, int m, float c){
 
688
        float* restrict fl = (float*) l;
 
689
        float* restrict fw = (float*) w;
 
690
        int i;
 
691
        for(i = 0; i < m; ++i){
 
692
                float tmp1_1 = fl[0];
 
693
                float tmp1_2 = fl[1];
 
694
                float tmp1_3 = fl[2];
 
695
                float tmp1_4 = fl[3];
 
696
                float tmp2_1 = fw[-4];
 
697
                float tmp2_2 = fw[-3];
 
698
                float tmp2_3 = fw[-2];
 
699
                float tmp2_4 = fw[-1];
 
700
                float tmp3_1 = fw[0];
 
701
                float tmp3_2 = fw[1];
 
702
                float tmp3_3 = fw[2];
 
703
                float tmp3_4 = fw[3];
 
704
                fw[-4] = tmp2_1 + ((tmp1_1 + tmp3_1) * c);
 
705
                fw[-3] = tmp2_2 + ((tmp1_2 + tmp3_2) * c);
 
706
                fw[-2] = tmp2_3 + ((tmp1_3 + tmp3_3) * c);
 
707
                fw[-1] = tmp2_4 + ((tmp1_4 + tmp3_4) * c);
 
708
                fl = fw;
 
709
                fw += 8;
 
710
        }
 
711
        if(m < k){
 
712
                float c1;
 
713
                float c2;
 
714
                float c3;
 
715
                float c4;
 
716
                c += c;
 
717
                c1 = fl[0] * c;
 
718
                c2 = fl[1] * c;
 
719
                c3 = fl[2] * c;
 
720
                c4 = fl[3] * c;
 
721
                for(; m < k; ++m){
 
722
                        float tmp1 = fw[-4];
 
723
                        float tmp2 = fw[-3];
 
724
                        float tmp3 = fw[-2];
 
725
                        float tmp4 = fw[-1];
 
726
                        fw[-4] = tmp1 + c1;
 
727
                        fw[-3] = tmp2 + c2;
 
728
                        fw[-2] = tmp3 + c3;
 
729
                        fw[-1] = tmp4 + c4;
 
730
                        fw += 8;
 
731
                }
 
732
        }
 
733
}
 
734
 
 
735
#endif
 
736
 
 
737
/* <summary>                             */
 
738
/* Inverse 9-7 wavelet transform in 1-D. */
 
739
/* </summary>                            */
 
740
static void v4dwt_decode(v4dwt_t* restrict dwt){
 
741
        int a, b;
 
742
        if(dwt->cas == 0) {
 
743
                if(!((dwt->dn > 0) || (dwt->sn > 1))){
 
744
                        return;
 
745
                }
 
746
                a = 0;
 
747
                b = 1;
 
748
        }else{
 
749
                if(!((dwt->sn > 0) || (dwt->dn > 1))) {
 
750
                        return;
 
751
                }
 
752
                a = 1;
 
753
                b = 0;
 
754
        }
 
755
#ifdef __SSE__
 
756
        v4dwt_decode_step1_sse(dwt->wavelet+a, dwt->sn, _mm_set1_ps(K));
 
757
        v4dwt_decode_step1_sse(dwt->wavelet+b, dwt->dn, _mm_set1_ps(c13318));
 
758
        v4dwt_decode_step2_sse(dwt->wavelet+b, dwt->wavelet+a+1, dwt->sn, int_min(dwt->sn, dwt->dn-a), _mm_set1_ps(dwt_delta));
 
759
        v4dwt_decode_step2_sse(dwt->wavelet+a, dwt->wavelet+b+1, dwt->dn, int_min(dwt->dn, dwt->sn-b), _mm_set1_ps(dwt_gamma));
 
760
        v4dwt_decode_step2_sse(dwt->wavelet+b, dwt->wavelet+a+1, dwt->sn, int_min(dwt->sn, dwt->dn-a), _mm_set1_ps(dwt_beta));
 
761
        v4dwt_decode_step2_sse(dwt->wavelet+a, dwt->wavelet+b+1, dwt->dn, int_min(dwt->dn, dwt->sn-b), _mm_set1_ps(dwt_alpha));
 
762
#else
 
763
        v4dwt_decode_step1(dwt->wavelet+a, dwt->sn, K);
 
764
        v4dwt_decode_step1(dwt->wavelet+b, dwt->dn, c13318);
 
765
        v4dwt_decode_step2(dwt->wavelet+b, dwt->wavelet+a+1, dwt->sn, int_min(dwt->sn, dwt->dn-a), dwt_delta);
 
766
        v4dwt_decode_step2(dwt->wavelet+a, dwt->wavelet+b+1, dwt->dn, int_min(dwt->dn, dwt->sn-b), dwt_gamma);
 
767
        v4dwt_decode_step2(dwt->wavelet+b, dwt->wavelet+a+1, dwt->sn, int_min(dwt->sn, dwt->dn-a), dwt_beta);
 
768
        v4dwt_decode_step2(dwt->wavelet+a, dwt->wavelet+b+1, dwt->dn, int_min(dwt->dn, dwt->sn-b), dwt_alpha);
 
769
#endif
 
770
}
 
771
 
 
772
/* <summary>                             */
 
773
/* Inverse 9-7 wavelet transform in 2-D. */
 
774
/* </summary>                            */
 
775
void dwt_decode_real(opj_tcd_tilecomp_t* restrict tilec, int numres){
 
776
        v4dwt_t h;
 
777
        v4dwt_t v;
 
778
 
 
779
        opj_tcd_resolution_t* res = tilec->resolutions;
 
780
 
 
781
        int rw = res->x1 - res->x0;     /* width of the resolution level computed */
 
782
        int rh = res->y1 - res->y0;     /* height of the resolution level computed */
 
783
 
 
784
        int w = tilec->x1 - tilec->x0;
 
785
 
 
786
        h.wavelet = (v4*) opj_aligned_malloc((dwt_decode_max_resolution(res, numres)+5) * sizeof(v4));
 
787
        v.wavelet = h.wavelet;
 
788
 
 
789
        while( --numres) {
 
790
                float * restrict aj = (float*) tilec->data;
 
791
                int bufsize = (tilec->x1 - tilec->x0) * (tilec->y1 - tilec->y0);
 
792
                int j;
 
793
 
 
794
                h.sn = rw;
 
795
                v.sn = rh;
 
796
 
 
797
                ++res;
 
798
 
 
799
                rw = res->x1 - res->x0; /* width of the resolution level computed */
 
800
                rh = res->y1 - res->y0; /* height of the resolution level computed */
 
801
 
 
802
                h.dn = rw - h.sn;
 
803
                h.cas = res->x0 % 2;
 
804
 
 
805
                for(j = rh; j > 3; j -= 4){
 
806
                        int k;
 
807
                        v4dwt_interleave_h(&h, aj, w, bufsize);
 
808
                        v4dwt_decode(&h);
 
809
                                for(k = rw; --k >= 0;){
 
810
                                        aj[k    ] = h.wavelet[k].f[0];
 
811
                                        aj[k+w  ] = h.wavelet[k].f[1];
 
812
                                        aj[k+w*2] = h.wavelet[k].f[2];
 
813
                                        aj[k+w*3] = h.wavelet[k].f[3];
 
814
                                }
 
815
                        aj += w*4;
 
816
                        bufsize -= w*4;
 
817
                }
 
818
                if (rh & 0x03) {
 
819
                                int k;
 
820
                        j = rh & 0x03;
 
821
                        v4dwt_interleave_h(&h, aj, w, bufsize);
 
822
                        v4dwt_decode(&h);
 
823
                                for(k = rw; --k >= 0;){
 
824
                                        switch(j) {
 
825
                                                case 3: aj[k+w*2] = h.wavelet[k].f[2];
 
826
                                                case 2: aj[k+w  ] = h.wavelet[k].f[1];
 
827
                                                case 1: aj[k    ] = h.wavelet[k].f[0];
 
828
                                        }
 
829
                                }
 
830
                        }
 
831
 
 
832
                v.dn = rh - v.sn;
 
833
                v.cas = res->y0 % 2;
 
834
 
 
835
                aj = (float*) tilec->data;
 
836
                for(j = rw; j > 3; j -= 4){
 
837
                        int k;
 
838
                        v4dwt_interleave_v(&v, aj, w);
 
839
                        v4dwt_decode(&v);
 
840
                                for(k = 0; k < rh; ++k){
 
841
                                        memcpy(&aj[k*w], &v.wavelet[k], 4 * sizeof(float));
 
842
                                }
 
843
                        aj += 4;
 
844
                }
 
845
                if (rw & 0x03){
 
846
                                int k;
 
847
                        j = rw & 0x03;
 
848
                        v4dwt_interleave_v(&v, aj, w);
 
849
                        v4dwt_decode(&v);
 
850
                                for(k = 0; k < rh; ++k){
 
851
                                        memcpy(&aj[k*w], &v.wavelet[k], j * sizeof(float));
 
852
                                }
 
853
                        }
 
854
        }
 
855
 
 
856
        opj_aligned_free(h.wavelet);
 
857
}
 
858