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