2 * The copyright in this software is being made available under the 2-clauses
3 * BSD License, included below. This software may be subject to other third
4 * party and contributor rights, including patent rights, and no such rights
5 * are granted under this license.
7 * Copyright (c) 2002-2014, Universite catholique de Louvain (UCL), Belgium
8 * Copyright (c) 2002-2014, Professor Benoit Macq
9 * Copyright (c) 2001-2003, David Janssens
10 * Copyright (c) 2002-2003, Yannick Verschueren
11 * Copyright (c) 2003-2007, Francois-Olivier Devaux
12 * Copyright (c) 2003-2014, Antonin Descampe
13 * Copyright (c) 2005, Herve Drolon, FreeImage Team
14 * Copyright (c) 2007, Jonathan Ballard <dzonatas@dzonux.net>
15 * Copyright (c) 2007, Callum Lerwick <seg@haxxed.com>
16 * Copyright (c) 2017, IntoPIX SA <support@intopix.com>
17 * All rights reserved.
19 * Redistribution and use in source and binary forms, with or without
20 * modification, are permitted provided that the following conditions
22 * 1. Redistributions of source code must retain the above copyright
23 * notice, this list of conditions and the following disclaimer.
24 * 2. Redistributions in binary form must reproduce the above copyright
25 * notice, this list of conditions and the following disclaimer in the
26 * documentation and/or other materials provided with the distribution.
28 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
29 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
30 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
31 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
32 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
33 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
34 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
35 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
36 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
37 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
38 * POSSIBILITY OF SUCH DAMAGE.
43 #define OPJ_SKIP_POISON
44 #include "opj_includes.h"
47 #include <xmmintrin.h>
50 #include <emmintrin.h>
53 #include <tmmintrin.h>
56 #include <immintrin.h>
60 #pragma GCC poison malloc calloc realloc free
63 /** @defgroup DWT DWT - Implementation of a discrete wavelet transform */
66 #define OPJ_WS(i) v->mem[(i)*2]
67 #define OPJ_WD(i) v->mem[(1+(i)*2)]
70 /** Number of int32 values in a AVX2 register */
71 #define VREG_INT_COUNT 8
73 /** Number of int32 values in a SSE2 register */
74 #define VREG_INT_COUNT 4
77 /** Number of columns that we can process in parallel in the vertical pass */
78 #define PARALLEL_COLS_53 (2*VREG_INT_COUNT)
80 /** @name Local data structures */
83 typedef struct dwt_local {
85 OPJ_INT32 dn; /* number of elements in high pass band */
86 OPJ_INT32 sn; /* number of elements in low pass band */
87 OPJ_INT32 cas; /* 0 = start on even coord, 1 = start on odd coord */
93 OPJ_FLOAT32 f[NB_ELTS_V8];
96 typedef struct v8dwt_local {
98 OPJ_INT32 dn ; /* number of elements in high pass band */
99 OPJ_INT32 sn ; /* number of elements in low pass band */
100 OPJ_INT32 cas ; /* 0 = start on even coord, 1 = start on odd coord */
101 OPJ_UINT32 win_l_x0; /* start coord in low pass band */
102 OPJ_UINT32 win_l_x1; /* end coord in low pass band */
103 OPJ_UINT32 win_h_x0; /* start coord in high pass band */
104 OPJ_UINT32 win_h_x1; /* end coord in high pass band */
107 /* From table F.4 from the standard */
108 static const OPJ_FLOAT32 opj_dwt_alpha = -1.586134342f;
109 static const OPJ_FLOAT32 opj_dwt_beta = -0.052980118f;
110 static const OPJ_FLOAT32 opj_dwt_gamma = 0.882911075f;
111 static const OPJ_FLOAT32 opj_dwt_delta = 0.443506852f;
113 static const OPJ_FLOAT32 opj_K = 1.230174105f;
114 static const OPJ_FLOAT32 opj_invK = (OPJ_FLOAT32)(1.0 / 1.230174105);
119 Virtual function type for wavelet transform in 1-D
121 typedef void (*DWT1DFN)(const opj_dwt_t* v);
123 /** @name Local static functions */
127 Forward lazy transform (horizontal)
129 static void opj_dwt_deinterleave_h(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
130 OPJ_INT32 sn, OPJ_INT32 cas);
132 Forward lazy transform (vertical)
134 static void opj_dwt_deinterleave_v(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
135 OPJ_INT32 sn, OPJ_UINT32 x, OPJ_INT32 cas);
137 Forward 5-3 wavelet transform in 1-D
139 static void opj_dwt_encode_1(void *a, OPJ_INT32 dn, OPJ_INT32 sn,
142 Forward 9-7 wavelet transform in 1-D
144 static void opj_dwt_encode_1_real(void *a, OPJ_INT32 dn, OPJ_INT32 sn,
147 Explicit calculation of the Quantization Stepsizes
149 static void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps,
150 opj_stepsize_t *bandno_stepsize);
152 Inverse wavelet transform in 2-D.
154 static OPJ_BOOL opj_dwt_decode_tile(opj_thread_pool_t* tp,
155 opj_tcd_tilecomp_t* tilec, OPJ_UINT32 i);
157 static OPJ_BOOL opj_dwt_decode_partial_tile(
158 opj_tcd_tilecomp_t* tilec,
161 /* Where void* is a OPJ_INT32* for 5x3 and OPJ_FLOAT32* for 9x7 */
162 typedef void (*opj_encode_one_row_fnptr_type)(void *, OPJ_INT32, OPJ_INT32,
165 static OPJ_BOOL opj_dwt_encode_procedure(opj_thread_pool_t* tp,
166 opj_tcd_tilecomp_t * tilec,
167 opj_encode_one_row_fnptr_type p_function);
169 static OPJ_UINT32 opj_dwt_max_resolution(opj_tcd_resolution_t* OPJ_RESTRICT r,
173 /* Inverse 9-7 wavelet transform in 1-D. */
180 #define OPJ_S(i) a[(i)*2]
181 #define OPJ_D(i) a[(1+(i)*2)]
182 #define OPJ_S_(i) ((i)<0?OPJ_S(0):((i)>=sn?OPJ_S(sn-1):OPJ_S(i)))
183 #define OPJ_D_(i) ((i)<0?OPJ_D(0):((i)>=dn?OPJ_D(dn-1):OPJ_D(i)))
185 #define OPJ_SS_(i) ((i)<0?OPJ_S(0):((i)>=dn?OPJ_S(dn-1):OPJ_S(i)))
186 #define OPJ_DD_(i) ((i)<0?OPJ_D(0):((i)>=sn?OPJ_D(sn-1):OPJ_D(i)))
189 /* This table contains the norms of the 5-3 wavelets for different bands. */
191 /* FIXME! the array should really be extended up to 33 resolution levels */
192 /* See https://github.com/uclouvain/openjpeg/issues/493 */
193 static const OPJ_FLOAT64 opj_dwt_norms[4][10] = {
194 {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
195 {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
196 {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
197 {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
201 /* This table contains the norms of the 9-7 wavelets for different bands. */
203 /* FIXME! the array should really be extended up to 33 resolution levels */
204 /* See https://github.com/uclouvain/openjpeg/issues/493 */
205 static const OPJ_FLOAT64 opj_dwt_norms_real[4][10] = {
206 {1.000, 1.965, 4.177, 8.403, 16.90, 33.84, 67.69, 135.3, 270.6, 540.9},
207 {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
208 {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
209 {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2}
213 ==========================================================
215 ==========================================================
219 /* Forward lazy transform (horizontal). */
221 static void opj_dwt_deinterleave_h(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
222 OPJ_INT32 sn, OPJ_INT32 cas)
225 OPJ_INT32 * l_dest = b;
226 OPJ_INT32 * l_src = a + cas;
228 for (i = 0; i < sn; ++i) {
236 for (i = 0; i < dn; ++i) {
243 /* Forward lazy transform (vertical). */
245 static void opj_dwt_deinterleave_v(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
246 OPJ_INT32 sn, OPJ_UINT32 x, OPJ_INT32 cas)
249 OPJ_INT32 * l_dest = b;
250 OPJ_INT32 * l_src = a + cas;
256 } /* b[i*x]=a[2*i+cas]; */
258 l_dest = b + (OPJ_SIZE_T)sn * (OPJ_SIZE_T)x;
266 } /*b[(sn+i)*x]=a[(2*i+1-cas)];*/
269 #ifdef STANDARD_SLOW_VERSION
271 /* Inverse lazy transform (horizontal). */
273 static void opj_dwt_interleave_h(const opj_dwt_t* h, OPJ_INT32 *a)
276 OPJ_INT32 *bi = h->mem + h->cas;
283 bi = h->mem + 1 - h->cas;
292 /* Inverse lazy transform (vertical). */
294 static void opj_dwt_interleave_v(const opj_dwt_t* v, OPJ_INT32 *a, OPJ_INT32 x)
297 OPJ_INT32 *bi = v->mem + v->cas;
304 ai = a + (v->sn * (OPJ_SIZE_T)x);
305 bi = v->mem + 1 - v->cas;
314 #endif /* STANDARD_SLOW_VERSION */
317 /* Forward 5-3 wavelet transform in 1-D. */
319 static void opj_dwt_encode_1(void *aIn, OPJ_INT32 dn, OPJ_INT32 sn,
323 OPJ_INT32* a = (OPJ_INT32*)aIn;
326 if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
327 for (i = 0; i < dn; i++) {
328 OPJ_D(i) -= (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
330 for (i = 0; i < sn; i++) {
331 OPJ_S(i) += (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
335 if (!sn && dn == 1) { /* NEW : CASE ONE ELEMENT */
338 for (i = 0; i < dn; i++) {
339 OPJ_S(i) -= (OPJ_DD_(i) + OPJ_DD_(i - 1)) >> 1;
341 for (i = 0; i < sn; i++) {
342 OPJ_D(i) += (OPJ_SS_(i) + OPJ_SS_(i + 1) + 2) >> 2;
348 #ifdef STANDARD_SLOW_VERSION
350 /* Inverse 5-3 wavelet transform in 1-D. */
352 static void opj_dwt_decode_1_(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
358 if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
359 for (i = 0; i < sn; i++) {
360 OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
362 for (i = 0; i < dn; i++) {
363 OPJ_D(i) += (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
367 if (!sn && dn == 1) { /* NEW : CASE ONE ELEMENT */
370 for (i = 0; i < sn; i++) {
371 OPJ_D(i) -= (OPJ_SS_(i) + OPJ_SS_(i + 1) + 2) >> 2;
373 for (i = 0; i < dn; i++) {
374 OPJ_S(i) += (OPJ_DD_(i) + OPJ_DD_(i - 1)) >> 1;
380 static void opj_dwt_decode_1(const opj_dwt_t *v)
382 opj_dwt_decode_1_(v->mem, v->dn, v->sn, v->cas);
385 #endif /* STANDARD_SLOW_VERSION */
387 #if !defined(STANDARD_SLOW_VERSION)
388 static void opj_idwt53_h_cas0(OPJ_INT32* tmp,
394 const OPJ_INT32* in_even = &tiledp[0];
395 const OPJ_INT32* in_odd = &tiledp[sn];
397 #ifdef TWO_PASS_VERSION
398 /* For documentation purpose: performs lifting in two iterations, */
399 /* but without explicit interleaving */
404 tmp[0] = in_even[0] - ((in_odd[0] + 1) >> 1);
405 for (i = 2, j = 0; i <= len - 2; i += 2, j++) {
406 tmp[i] = in_even[j + 1] - ((in_odd[j] + in_odd[j + 1] + 2) >> 2);
408 if (len & 1) { /* if len is odd */
409 tmp[len - 1] = in_even[(len - 1) / 2] - ((in_odd[(len - 2) / 2] + 1) >> 1);
413 for (i = 1, j = 0; i < len - 1; i += 2, j++) {
414 tmp[i] = in_odd[j] + ((tmp[i - 1] + tmp[i + 1]) >> 1);
416 if (!(len & 1)) { /* if len is even */
417 tmp[len - 1] = in_odd[(len - 1) / 2] + tmp[len - 2];
420 OPJ_INT32 d1c, d1n, s1n, s0c, s0n;
424 /* Improved version of the TWO_PASS_VERSION: */
425 /* Performs lifting in one single iteration. Saves memory */
426 /* accesses and explicit interleaving. */
429 s0n = s1n - ((d1n + 1) >> 1);
431 for (i = 0, j = 1; i < (len - 3); i += 2, j++) {
438 s0n = s1n - ((d1c + d1n + 2) >> 2);
441 tmp[i + 1] = d1c + ((s0c + s0n) >> 1);
447 tmp[len - 1] = in_even[(len - 1) / 2] - ((d1n + 1) >> 1);
448 tmp[len - 2] = d1n + ((s0n + tmp[len - 1]) >> 1);
450 tmp[len - 1] = d1n + s0n;
453 memcpy(tiledp, tmp, (OPJ_UINT32)len * sizeof(OPJ_INT32));
456 static void opj_idwt53_h_cas1(OPJ_INT32* tmp,
462 const OPJ_INT32* in_even = &tiledp[sn];
463 const OPJ_INT32* in_odd = &tiledp[0];
465 #ifdef TWO_PASS_VERSION
466 /* For documentation purpose: performs lifting in two iterations, */
467 /* but without explicit interleaving */
472 for (i = 1, j = 0; i < len - 1; i += 2, j++) {
473 tmp[i] = in_odd[j] - ((in_even[j] + in_even[j + 1] + 2) >> 2);
476 tmp[len - 1] = in_odd[len / 2 - 1] - ((in_even[len / 2 - 1] + 1) >> 1);
480 tmp[0] = in_even[0] + tmp[1];
481 for (i = 2, j = 1; i < len - 1; i += 2, j++) {
482 tmp[i] = in_even[j] + ((tmp[i + 1] + tmp[i - 1]) >> 1);
485 tmp[len - 1] = in_even[len / 2] + tmp[len - 2];
488 OPJ_INT32 s1, s2, dc, dn;
492 /* Improved version of the TWO_PASS_VERSION: */
493 /* Performs lifting in one single iteration. Saves memory */
494 /* accesses and explicit interleaving. */
497 dc = in_odd[0] - ((in_even[0] + s1 + 2) >> 2);
498 tmp[0] = in_even[0] + dc;
500 for (i = 1, j = 1; i < (len - 2 - !(len & 1)); i += 2, j++) {
504 dn = in_odd[j] - ((s1 + s2 + 2) >> 2);
506 tmp[i + 1] = s1 + ((dn + dc) >> 1);
515 dn = in_odd[len / 2 - 1] - ((s1 + 1) >> 1);
516 tmp[len - 2] = s1 + ((dn + dc) >> 1);
519 tmp[len - 1] = s1 + dc;
522 memcpy(tiledp, tmp, (OPJ_UINT32)len * sizeof(OPJ_INT32));
526 #endif /* !defined(STANDARD_SLOW_VERSION) */
529 /* Inverse 5-3 wavelet transform in 1-D for one row. */
531 /* Performs interleave, inverse wavelet transform and copy back to buffer */
532 static void opj_idwt53_h(const opj_dwt_t *dwt,
535 #ifdef STANDARD_SLOW_VERSION
536 /* For documentation purpose */
537 opj_dwt_interleave_h(dwt, tiledp);
538 opj_dwt_decode_1(dwt);
539 memcpy(tiledp, dwt->mem, (OPJ_UINT32)(dwt->sn + dwt->dn) * sizeof(OPJ_INT32));
541 const OPJ_INT32 sn = dwt->sn;
542 const OPJ_INT32 len = sn + dwt->dn;
543 if (dwt->cas == 0) { /* Left-most sample is on even coordinate */
545 opj_idwt53_h_cas0(dwt->mem, sn, len, tiledp);
547 /* Unmodified value */
549 } else { /* Left-most sample is on odd coordinate */
552 } else if (len == 2) {
553 OPJ_INT32* out = dwt->mem;
554 const OPJ_INT32* in_even = &tiledp[sn];
555 const OPJ_INT32* in_odd = &tiledp[0];
556 out[1] = in_odd[0] - ((in_even[0] + 1) >> 1);
557 out[0] = in_even[0] + out[1];
558 memcpy(tiledp, dwt->mem, (OPJ_UINT32)len * sizeof(OPJ_INT32));
559 } else if (len > 2) {
560 opj_idwt53_h_cas1(dwt->mem, sn, len, tiledp);
566 #if (defined(__SSE2__) || defined(__AVX2__)) && !defined(STANDARD_SLOW_VERSION)
568 /* Conveniency macros to improve the readabilty of the formulas */
571 #define LOAD_CST(x) _mm256_set1_epi32(x)
572 #define LOAD(x) _mm256_load_si256((const VREG*)(x))
573 #define LOADU(x) _mm256_loadu_si256((const VREG*)(x))
574 #define STORE(x,y) _mm256_store_si256((VREG*)(x),(y))
575 #define STOREU(x,y) _mm256_storeu_si256((VREG*)(x),(y))
576 #define ADD(x,y) _mm256_add_epi32((x),(y))
577 #define SUB(x,y) _mm256_sub_epi32((x),(y))
578 #define SAR(x,y) _mm256_srai_epi32((x),(y))
581 #define LOAD_CST(x) _mm_set1_epi32(x)
582 #define LOAD(x) _mm_load_si128((const VREG*)(x))
583 #define LOADU(x) _mm_loadu_si128((const VREG*)(x))
584 #define STORE(x,y) _mm_store_si128((VREG*)(x),(y))
585 #define STOREU(x,y) _mm_storeu_si128((VREG*)(x),(y))
586 #define ADD(x,y) _mm_add_epi32((x),(y))
587 #define SUB(x,y) _mm_sub_epi32((x),(y))
588 #define SAR(x,y) _mm_srai_epi32((x),(y))
590 #define ADD3(x,y,z) ADD(ADD(x,y),z)
593 void opj_idwt53_v_final_memcpy(OPJ_INT32* tiledp_col,
594 const OPJ_INT32* tmp,
599 for (i = 0; i < len; ++i) {
600 /* A memcpy(&tiledp_col[i * stride + 0],
601 &tmp[PARALLEL_COLS_53 * i + 0],
602 PARALLEL_COLS_53 * sizeof(OPJ_INT32))
603 would do but would be a tiny bit slower.
604 We can take here advantage of our knowledge of alignment */
605 STOREU(&tiledp_col[(OPJ_SIZE_T)i * stride + 0],
606 LOAD(&tmp[PARALLEL_COLS_53 * i + 0]));
607 STOREU(&tiledp_col[(OPJ_SIZE_T)i * stride + VREG_INT_COUNT],
608 LOAD(&tmp[PARALLEL_COLS_53 * i + VREG_INT_COUNT]));
612 /** Vertical inverse 5x3 wavelet transform for 8 columns in SSE2, or
613 * 16 in AVX2, when top-most pixel is on even coordinate */
614 static void opj_idwt53_v_cas0_mcols_SSE2_OR_AVX2(
618 OPJ_INT32* tiledp_col,
619 const OPJ_SIZE_T stride)
621 const OPJ_INT32* in_even = &tiledp_col[0];
622 const OPJ_INT32* in_odd = &tiledp_col[(OPJ_SIZE_T)sn * stride];
626 VREG d1c_0, d1n_0, s1n_0, s0c_0, s0n_0;
627 VREG d1c_1, d1n_1, s1n_1, s0c_1, s0n_1;
628 const VREG two = LOAD_CST(2);
632 assert(PARALLEL_COLS_53 == 16);
633 assert(VREG_INT_COUNT == 8);
635 assert(PARALLEL_COLS_53 == 8);
636 assert(VREG_INT_COUNT == 4);
639 /* Note: loads of input even/odd values must be done in a unaligned */
640 /* fashion. But stores in tmp can be done with aligned store, since */
641 /* the temporary buffer is properly aligned */
642 assert((OPJ_SIZE_T)tmp % (sizeof(OPJ_INT32) * VREG_INT_COUNT) == 0);
644 s1n_0 = LOADU(in_even + 0);
645 s1n_1 = LOADU(in_even + VREG_INT_COUNT);
646 d1n_0 = LOADU(in_odd);
647 d1n_1 = LOADU(in_odd + VREG_INT_COUNT);
649 /* s0n = s1n - ((d1n + 1) >> 1); <==> */
650 /* s0n = s1n - ((d1n + d1n + 2) >> 2); */
651 s0n_0 = SUB(s1n_0, SAR(ADD3(d1n_0, d1n_0, two), 2));
652 s0n_1 = SUB(s1n_1, SAR(ADD3(d1n_1, d1n_1, two), 2));
654 for (i = 0, j = 1; i < (len - 3); i += 2, j++) {
660 s1n_0 = LOADU(in_even + j * stride);
661 s1n_1 = LOADU(in_even + j * stride + VREG_INT_COUNT);
662 d1n_0 = LOADU(in_odd + j * stride);
663 d1n_1 = LOADU(in_odd + j * stride + VREG_INT_COUNT);
665 /*s0n = s1n - ((d1c + d1n + 2) >> 2);*/
666 s0n_0 = SUB(s1n_0, SAR(ADD3(d1c_0, d1n_0, two), 2));
667 s0n_1 = SUB(s1n_1, SAR(ADD3(d1c_1, d1n_1, two), 2));
669 STORE(tmp + PARALLEL_COLS_53 * (i + 0), s0c_0);
670 STORE(tmp + PARALLEL_COLS_53 * (i + 0) + VREG_INT_COUNT, s0c_1);
672 /* d1c + ((s0c + s0n) >> 1) */
673 STORE(tmp + PARALLEL_COLS_53 * (i + 1) + 0,
674 ADD(d1c_0, SAR(ADD(s0c_0, s0n_0), 1)));
675 STORE(tmp + PARALLEL_COLS_53 * (i + 1) + VREG_INT_COUNT,
676 ADD(d1c_1, SAR(ADD(s0c_1, s0n_1), 1)));
679 STORE(tmp + PARALLEL_COLS_53 * (i + 0) + 0, s0n_0);
680 STORE(tmp + PARALLEL_COLS_53 * (i + 0) + VREG_INT_COUNT, s0n_1);
683 VREG tmp_len_minus_1;
684 s1n_0 = LOADU(in_even + (OPJ_SIZE_T)((len - 1) / 2) * stride);
685 /* tmp_len_minus_1 = s1n - ((d1n + 1) >> 1); */
686 tmp_len_minus_1 = SUB(s1n_0, SAR(ADD3(d1n_0, d1n_0, two), 2));
687 STORE(tmp + PARALLEL_COLS_53 * (len - 1), tmp_len_minus_1);
688 /* d1n + ((s0n + tmp_len_minus_1) >> 1) */
689 STORE(tmp + PARALLEL_COLS_53 * (len - 2),
690 ADD(d1n_0, SAR(ADD(s0n_0, tmp_len_minus_1), 1)));
692 s1n_1 = LOADU(in_even + (OPJ_SIZE_T)((len - 1) / 2) * stride + VREG_INT_COUNT);
693 /* tmp_len_minus_1 = s1n - ((d1n + 1) >> 1); */
694 tmp_len_minus_1 = SUB(s1n_1, SAR(ADD3(d1n_1, d1n_1, two), 2));
695 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + VREG_INT_COUNT,
697 /* d1n + ((s0n + tmp_len_minus_1) >> 1) */
698 STORE(tmp + PARALLEL_COLS_53 * (len - 2) + VREG_INT_COUNT,
699 ADD(d1n_1, SAR(ADD(s0n_1, tmp_len_minus_1), 1)));
703 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 0,
705 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + VREG_INT_COUNT,
709 opj_idwt53_v_final_memcpy(tiledp_col, tmp, len, stride);
713 /** Vertical inverse 5x3 wavelet transform for 8 columns in SSE2, or
714 * 16 in AVX2, when top-most pixel is on odd coordinate */
715 static void opj_idwt53_v_cas1_mcols_SSE2_OR_AVX2(
719 OPJ_INT32* tiledp_col,
720 const OPJ_SIZE_T stride)
725 VREG s1_0, s2_0, dc_0, dn_0;
726 VREG s1_1, s2_1, dc_1, dn_1;
727 const VREG two = LOAD_CST(2);
729 const OPJ_INT32* in_even = &tiledp_col[(OPJ_SIZE_T)sn * stride];
730 const OPJ_INT32* in_odd = &tiledp_col[0];
734 assert(PARALLEL_COLS_53 == 16);
735 assert(VREG_INT_COUNT == 8);
737 assert(PARALLEL_COLS_53 == 8);
738 assert(VREG_INT_COUNT == 4);
741 /* Note: loads of input even/odd values must be done in a unaligned */
742 /* fashion. But stores in tmp can be done with aligned store, since */
743 /* the temporary buffer is properly aligned */
744 assert((OPJ_SIZE_T)tmp % (sizeof(OPJ_INT32) * VREG_INT_COUNT) == 0);
746 s1_0 = LOADU(in_even + stride);
747 /* in_odd[0] - ((in_even[0] + s1 + 2) >> 2); */
748 dc_0 = SUB(LOADU(in_odd + 0),
749 SAR(ADD3(LOADU(in_even + 0), s1_0, two), 2));
750 STORE(tmp + PARALLEL_COLS_53 * 0, ADD(LOADU(in_even + 0), dc_0));
752 s1_1 = LOADU(in_even + stride + VREG_INT_COUNT);
753 /* in_odd[0] - ((in_even[0] + s1 + 2) >> 2); */
754 dc_1 = SUB(LOADU(in_odd + VREG_INT_COUNT),
755 SAR(ADD3(LOADU(in_even + VREG_INT_COUNT), s1_1, two), 2));
756 STORE(tmp + PARALLEL_COLS_53 * 0 + VREG_INT_COUNT,
757 ADD(LOADU(in_even + VREG_INT_COUNT), dc_1));
759 for (i = 1, j = 1; i < (len - 2 - !(len & 1)); i += 2, j++) {
761 s2_0 = LOADU(in_even + (j + 1) * stride);
762 s2_1 = LOADU(in_even + (j + 1) * stride + VREG_INT_COUNT);
764 /* dn = in_odd[j * stride] - ((s1 + s2 + 2) >> 2); */
765 dn_0 = SUB(LOADU(in_odd + j * stride),
766 SAR(ADD3(s1_0, s2_0, two), 2));
767 dn_1 = SUB(LOADU(in_odd + j * stride + VREG_INT_COUNT),
768 SAR(ADD3(s1_1, s2_1, two), 2));
770 STORE(tmp + PARALLEL_COLS_53 * i, dc_0);
771 STORE(tmp + PARALLEL_COLS_53 * i + VREG_INT_COUNT, dc_1);
773 /* tmp[i + 1] = s1 + ((dn + dc) >> 1); */
774 STORE(tmp + PARALLEL_COLS_53 * (i + 1) + 0,
775 ADD(s1_0, SAR(ADD(dn_0, dc_0), 1)));
776 STORE(tmp + PARALLEL_COLS_53 * (i + 1) + VREG_INT_COUNT,
777 ADD(s1_1, SAR(ADD(dn_1, dc_1), 1)));
784 STORE(tmp + PARALLEL_COLS_53 * i, dc_0);
785 STORE(tmp + PARALLEL_COLS_53 * i + VREG_INT_COUNT, dc_1);
788 /*dn = in_odd[(len / 2 - 1) * stride] - ((s1 + 1) >> 1); */
789 dn_0 = SUB(LOADU(in_odd + (OPJ_SIZE_T)(len / 2 - 1) * stride),
790 SAR(ADD3(s1_0, s1_0, two), 2));
791 dn_1 = SUB(LOADU(in_odd + (OPJ_SIZE_T)(len / 2 - 1) * stride + VREG_INT_COUNT),
792 SAR(ADD3(s1_1, s1_1, two), 2));
794 /* tmp[len - 2] = s1 + ((dn + dc) >> 1); */
795 STORE(tmp + PARALLEL_COLS_53 * (len - 2) + 0,
796 ADD(s1_0, SAR(ADD(dn_0, dc_0), 1)));
797 STORE(tmp + PARALLEL_COLS_53 * (len - 2) + VREG_INT_COUNT,
798 ADD(s1_1, SAR(ADD(dn_1, dc_1), 1)));
800 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 0, dn_0);
801 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + VREG_INT_COUNT, dn_1);
803 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 0, ADD(s1_0, dc_0));
804 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + VREG_INT_COUNT,
808 opj_idwt53_v_final_memcpy(tiledp_col, tmp, len, stride);
822 #endif /* (defined(__SSE2__) || defined(__AVX2__)) && !defined(STANDARD_SLOW_VERSION) */
824 #if !defined(STANDARD_SLOW_VERSION)
825 /** Vertical inverse 5x3 wavelet transform for one column, when top-most
826 * pixel is on even coordinate */
827 static void opj_idwt3_v_cas0(OPJ_INT32* tmp,
830 OPJ_INT32* tiledp_col,
831 const OPJ_SIZE_T stride)
834 OPJ_INT32 d1c, d1n, s1n, s0c, s0n;
838 /* Performs lifting in one single iteration. Saves memory */
839 /* accesses and explicit interleaving. */
842 d1n = tiledp_col[(OPJ_SIZE_T)sn * stride];
843 s0n = s1n - ((d1n + 1) >> 1);
845 for (i = 0, j = 0; i < (len - 3); i += 2, j++) {
849 s1n = tiledp_col[(OPJ_SIZE_T)(j + 1) * stride];
850 d1n = tiledp_col[(OPJ_SIZE_T)(sn + j + 1) * stride];
852 s0n = s1n - ((d1c + d1n + 2) >> 2);
855 tmp[i + 1] = d1c + ((s0c + s0n) >> 1);
862 tiledp_col[(OPJ_SIZE_T)((len - 1) / 2) * stride] -
864 tmp[len - 2] = d1n + ((s0n + tmp[len - 1]) >> 1);
866 tmp[len - 1] = d1n + s0n;
869 for (i = 0; i < len; ++i) {
870 tiledp_col[(OPJ_SIZE_T)i * stride] = tmp[i];
875 /** Vertical inverse 5x3 wavelet transform for one column, when top-most
876 * pixel is on odd coordinate */
877 static void opj_idwt3_v_cas1(OPJ_INT32* tmp,
880 OPJ_INT32* tiledp_col,
881 const OPJ_SIZE_T stride)
884 OPJ_INT32 s1, s2, dc, dn;
885 const OPJ_INT32* in_even = &tiledp_col[(OPJ_SIZE_T)sn * stride];
886 const OPJ_INT32* in_odd = &tiledp_col[0];
890 /* Performs lifting in one single iteration. Saves memory */
891 /* accesses and explicit interleaving. */
893 s1 = in_even[stride];
894 dc = in_odd[0] - ((in_even[0] + s1 + 2) >> 2);
895 tmp[0] = in_even[0] + dc;
896 for (i = 1, j = 1; i < (len - 2 - !(len & 1)); i += 2, j++) {
898 s2 = in_even[(OPJ_SIZE_T)(j + 1) * stride];
900 dn = in_odd[(OPJ_SIZE_T)j * stride] - ((s1 + s2 + 2) >> 2);
902 tmp[i + 1] = s1 + ((dn + dc) >> 1);
909 dn = in_odd[(OPJ_SIZE_T)(len / 2 - 1) * stride] - ((s1 + 1) >> 1);
910 tmp[len - 2] = s1 + ((dn + dc) >> 1);
913 tmp[len - 1] = s1 + dc;
916 for (i = 0; i < len; ++i) {
917 tiledp_col[(OPJ_SIZE_T)i * stride] = tmp[i];
920 #endif /* !defined(STANDARD_SLOW_VERSION) */
923 /* Inverse vertical 5-3 wavelet transform in 1-D for several columns. */
925 /* Performs interleave, inverse wavelet transform and copy back to buffer */
926 static void opj_idwt53_v(const opj_dwt_t *dwt,
927 OPJ_INT32* tiledp_col,
931 #ifdef STANDARD_SLOW_VERSION
932 /* For documentation purpose */
934 for (c = 0; c < nb_cols; c ++) {
935 opj_dwt_interleave_v(dwt, tiledp_col + c, stride);
936 opj_dwt_decode_1(dwt);
937 for (k = 0; k < dwt->sn + dwt->dn; ++k) {
938 tiledp_col[c + k * stride] = dwt->mem[k];
942 const OPJ_INT32 sn = dwt->sn;
943 const OPJ_INT32 len = sn + dwt->dn;
945 /* If len == 1, unmodified value */
947 #if (defined(__SSE2__) || defined(__AVX2__))
948 if (len > 1 && nb_cols == PARALLEL_COLS_53) {
949 /* Same as below general case, except that thanks to SSE2/AVX2 */
950 /* we can efficiently process 8/16 columns in parallel */
951 opj_idwt53_v_cas0_mcols_SSE2_OR_AVX2(dwt->mem, sn, len, tiledp_col, stride);
957 for (c = 0; c < nb_cols; c++, tiledp_col++) {
958 opj_idwt3_v_cas0(dwt->mem, sn, len, tiledp_col, stride);
965 for (c = 0; c < nb_cols; c++, tiledp_col++) {
973 OPJ_INT32* out = dwt->mem;
974 for (c = 0; c < nb_cols; c++, tiledp_col++) {
976 const OPJ_INT32* in_even = &tiledp_col[(OPJ_SIZE_T)sn * stride];
977 const OPJ_INT32* in_odd = &tiledp_col[0];
979 out[1] = in_odd[0] - ((in_even[0] + 1) >> 1);
980 out[0] = in_even[0] + out[1];
982 for (i = 0; i < len; ++i) {
983 tiledp_col[(OPJ_SIZE_T)i * stride] = out[i];
990 #if (defined(__SSE2__) || defined(__AVX2__))
991 if (len > 2 && nb_cols == PARALLEL_COLS_53) {
992 /* Same as below general case, except that thanks to SSE2/AVX2 */
993 /* we can efficiently process 8/16 columns in parallel */
994 opj_idwt53_v_cas1_mcols_SSE2_OR_AVX2(dwt->mem, sn, len, tiledp_col, stride);
1000 for (c = 0; c < nb_cols; c++, tiledp_col++) {
1001 opj_idwt3_v_cas1(dwt->mem, sn, len, tiledp_col, stride);
1009 static void opj_dwt_encode_step1(OPJ_FLOAT32* fw,
1012 const OPJ_FLOAT32 c)
1015 for (i = start; i < end; ++i) {
1019 static void opj_dwt_encode_step2(OPJ_FLOAT32* fl, OPJ_FLOAT32* fw,
1026 OPJ_UINT32 imax = opj_uint_min(end, m);
1031 for (i = start; i < imax; ++i) {
1032 fw[-1] += (fl[0] + fw[0]) * c;
1037 assert(m + 1 == end);
1038 fw[-1] += (2 * fl[0]) * c;
1042 static void opj_dwt_encode_1_real(void *aIn, OPJ_INT32 dn, OPJ_INT32 sn,
1045 OPJ_FLOAT32* w = (OPJ_FLOAT32*)aIn;
1048 if (!((dn > 0) || (sn > 1))) {
1054 if (!((sn > 0) || (dn > 1))) {
1060 opj_dwt_encode_step2(w + a, w + b + 1,
1062 (OPJ_UINT32)opj_int_min(dn, sn - b),
1064 opj_dwt_encode_step2(w + b, w + a + 1,
1066 (OPJ_UINT32)opj_int_min(sn, dn - a),
1068 opj_dwt_encode_step2(w + a, w + b + 1,
1070 (OPJ_UINT32)opj_int_min(dn, sn - b),
1072 opj_dwt_encode_step2(w + b, w + a + 1,
1074 (OPJ_UINT32)opj_int_min(sn, dn - a),
1076 opj_dwt_encode_step1(w + b, 0, (OPJ_UINT32)dn,
1078 opj_dwt_encode_step1(w + a, 0, (OPJ_UINT32)sn,
1082 static void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps,
1083 opj_stepsize_t *bandno_stepsize)
1086 p = opj_int_floorlog2(stepsize) - 13;
1087 n = 11 - opj_int_floorlog2(stepsize);
1088 bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
1089 bandno_stepsize->expn = numbps - p;
1093 ==========================================================
1095 ==========================================================
1103 OPJ_INT32 * OPJ_RESTRICT tiledp;
1106 opj_encode_one_row_fnptr_type p_function;
1107 } opj_dwt_encode_h_job_t;
1109 static void opj_dwt_encode_h_func(void* user_data, opj_tls_t* tls)
1112 opj_dwt_encode_h_job_t* job;
1115 job = (opj_dwt_encode_h_job_t*)user_data;
1116 for (j = job->min_j; j < job->max_j; j++) {
1117 OPJ_INT32* OPJ_RESTRICT aj = job->tiledp + j * job->w;
1119 for (k = 0; k < job->rw; k++) {
1120 job->h.mem[k] = aj[k];
1122 (*job->p_function)(job->h.mem, job->h.dn, job->h.sn, job->h.cas);
1123 opj_dwt_deinterleave_h(job->h.mem, aj, job->h.dn, job->h.sn, job->h.cas);
1126 opj_aligned_free(job->h.mem);
1134 OPJ_INT32 * OPJ_RESTRICT tiledp;
1137 opj_encode_one_row_fnptr_type p_function;
1138 } opj_dwt_encode_v_job_t;
1140 static void opj_dwt_encode_v_func(void* user_data, opj_tls_t* tls)
1143 opj_dwt_encode_v_job_t* job;
1146 job = (opj_dwt_encode_v_job_t*)user_data;
1147 for (j = job->min_j; j < job->max_j; j++) {
1148 OPJ_INT32* OPJ_RESTRICT aj = job->tiledp + j;
1150 for (k = 0; k < job->rh; ++k) {
1151 job->v.mem[k] = aj[k * job->w];
1154 (*job->p_function)(job->v.mem, job->v.dn, job->v.sn, job->v.cas);
1156 opj_dwt_deinterleave_v(job->v.mem, aj, job->v.dn, job->v.sn, job->w,
1160 opj_aligned_free(job->v.mem);
1165 /* Forward 5-3 wavelet transform in 2-D. */
1167 static INLINE OPJ_BOOL opj_dwt_encode_procedure(opj_thread_pool_t* tp,
1168 opj_tcd_tilecomp_t * tilec,
1169 opj_encode_one_row_fnptr_type p_function)
1176 OPJ_SIZE_T l_data_size;
1178 opj_tcd_resolution_t * l_cur_res = 0;
1179 opj_tcd_resolution_t * l_last_res = 0;
1180 const int num_threads = opj_thread_pool_get_thread_count(tp);
1181 OPJ_INT32 * OPJ_RESTRICT tiledp = tilec->data;
1183 w = (OPJ_UINT32)(tilec->x1 - tilec->x0);
1184 l = (OPJ_INT32)tilec->numresolutions - 1;
1186 l_cur_res = tilec->resolutions + l;
1187 l_last_res = l_cur_res - 1;
1189 l_data_size = opj_dwt_max_resolution(tilec->resolutions, tilec->numresolutions);
1190 /* overflow check */
1191 if (l_data_size > (SIZE_MAX / sizeof(OPJ_INT32))) {
1192 /* FIXME event manager error callback */
1195 l_data_size *= sizeof(OPJ_INT32);
1196 bj = (OPJ_INT32*)opj_aligned_32_malloc(l_data_size);
1197 /* l_data_size is equal to 0 when numresolutions == 1 but bj is not used */
1198 /* in that case, so do not error out */
1199 if (l_data_size != 0 && ! bj) {
1206 OPJ_UINT32 rw; /* width of the resolution level computed */
1207 OPJ_UINT32 rh; /* height of the resolution level computed */
1209 rw1; /* width of the resolution level once lower than computed one */
1211 rh1; /* height of the resolution level once lower than computed one */
1212 OPJ_INT32 cas_col; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
1213 OPJ_INT32 cas_row; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */
1216 rw = (OPJ_UINT32)(l_cur_res->x1 - l_cur_res->x0);
1217 rh = (OPJ_UINT32)(l_cur_res->y1 - l_cur_res->y0);
1218 rw1 = (OPJ_UINT32)(l_last_res->x1 - l_last_res->x0);
1219 rh1 = (OPJ_UINT32)(l_last_res->y1 - l_last_res->y0);
1221 cas_row = l_cur_res->x0 & 1;
1222 cas_col = l_cur_res->y0 & 1;
1224 sn = (OPJ_INT32)rh1;
1225 dn = (OPJ_INT32)(rh - rh1);
1227 /* Perform vertical pass */
1228 if (num_threads <= 1 || rw <= 1) {
1229 for (j = 0; j < rw; ++j) {
1230 OPJ_INT32* OPJ_RESTRICT aj = tiledp + j;
1232 for (k = 0; k < rh; ++k) {
1236 (*p_function)(bj, dn, sn, cas_col);
1238 opj_dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
1241 OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
1244 if (rw < num_jobs) {
1247 step_j = (rw / num_jobs);
1249 for (j = 0; j < num_jobs; j++) {
1250 opj_dwt_encode_v_job_t* job;
1252 job = (opj_dwt_encode_v_job_t*) opj_malloc(sizeof(opj_dwt_encode_v_job_t));
1254 opj_thread_pool_wait_completion(tp, 0);
1255 opj_aligned_free(bj);
1258 job->v.mem = (OPJ_INT32*)opj_aligned_32_malloc(l_data_size);
1260 opj_thread_pool_wait_completion(tp, 0);
1262 opj_aligned_free(bj);
1267 job->v.cas = cas_col;
1270 job->tiledp = tiledp;
1271 job->min_j = j * step_j;
1272 job->max_j = (j + 1U) * step_j; /* this can overflow */
1273 if (j == (num_jobs - 1U)) { /* this will take care of the overflow */
1276 job->p_function = p_function;
1277 opj_thread_pool_submit_job(tp, opj_dwt_encode_v_func, job);
1279 opj_thread_pool_wait_completion(tp, 0);
1282 sn = (OPJ_INT32)rw1;
1283 dn = (OPJ_INT32)(rw - rw1);
1285 /* Perform horizontal pass */
1286 if (num_threads <= 1 || rh <= 1) {
1287 for (j = 0; j < rh; j++) {
1288 OPJ_INT32* OPJ_RESTRICT aj = tiledp + j * w;
1290 for (k = 0; k < rw; k++) {
1293 (*p_function)(bj, dn, sn, cas_row);
1294 opj_dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
1297 OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
1300 if (rh < num_jobs) {
1303 step_j = (rh / num_jobs);
1305 for (j = 0; j < num_jobs; j++) {
1306 opj_dwt_encode_h_job_t* job;
1308 job = (opj_dwt_encode_h_job_t*) opj_malloc(sizeof(opj_dwt_encode_h_job_t));
1310 opj_thread_pool_wait_completion(tp, 0);
1311 opj_aligned_free(bj);
1314 job->h.mem = (OPJ_INT32*)opj_aligned_32_malloc(l_data_size);
1316 opj_thread_pool_wait_completion(tp, 0);
1318 opj_aligned_free(bj);
1323 job->h.cas = cas_row;
1326 job->tiledp = tiledp;
1327 job->min_j = j * step_j;
1328 job->max_j = (j + 1U) * step_j; /* this can overflow */
1329 if (j == (num_jobs - 1U)) { /* this will take care of the overflow */
1332 job->p_function = p_function;
1333 opj_thread_pool_submit_job(tp, opj_dwt_encode_h_func, job);
1335 opj_thread_pool_wait_completion(tp, 0);
1338 l_cur_res = l_last_res;
1343 opj_aligned_free(bj);
1347 /* Forward 5-3 wavelet transform in 2-D. */
1349 OPJ_BOOL opj_dwt_encode(opj_tcd_t *p_tcd,
1350 opj_tcd_tilecomp_t * tilec)
1352 return opj_dwt_encode_procedure(p_tcd->thread_pool, tilec, opj_dwt_encode_1);
1356 /* Inverse 5-3 wavelet transform in 2-D. */
1358 OPJ_BOOL opj_dwt_decode(opj_tcd_t *p_tcd, opj_tcd_tilecomp_t* tilec,
1361 if (p_tcd->whole_tile_decoding) {
1362 return opj_dwt_decode_tile(p_tcd->thread_pool, tilec, numres);
1364 return opj_dwt_decode_partial_tile(tilec, numres);
1369 /* Get norm of 5-3 wavelet. */
1371 OPJ_FLOAT64 opj_dwt_getnorm(OPJ_UINT32 level, OPJ_UINT32 orient)
1373 /* FIXME ! This is just a band-aid to avoid a buffer overflow */
1374 /* but the array should really be extended up to 33 resolution levels */
1375 /* See https://github.com/uclouvain/openjpeg/issues/493 */
1376 if (orient == 0 && level >= 10) {
1378 } else if (orient > 0 && level >= 9) {
1381 return opj_dwt_norms[orient][level];
1385 /* Forward 9-7 wavelet transform in 2-D. */
1387 OPJ_BOOL opj_dwt_encode_real(opj_tcd_t *p_tcd,
1388 opj_tcd_tilecomp_t * tilec)
1390 return opj_dwt_encode_procedure(p_tcd->thread_pool, tilec,
1391 opj_dwt_encode_1_real);
1395 /* Get norm of 9-7 wavelet. */
1397 OPJ_FLOAT64 opj_dwt_getnorm_real(OPJ_UINT32 level, OPJ_UINT32 orient)
1399 /* FIXME ! This is just a band-aid to avoid a buffer overflow */
1400 /* but the array should really be extended up to 33 resolution levels */
1401 /* See https://github.com/uclouvain/openjpeg/issues/493 */
1402 if (orient == 0 && level >= 10) {
1404 } else if (orient > 0 && level >= 9) {
1407 return opj_dwt_norms_real[orient][level];
1410 void opj_dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, OPJ_UINT32 prec)
1412 OPJ_UINT32 numbands, bandno;
1413 numbands = 3 * tccp->numresolutions - 2;
1414 for (bandno = 0; bandno < numbands; bandno++) {
1415 OPJ_FLOAT64 stepsize;
1416 OPJ_UINT32 resno, level, orient, gain;
1418 resno = (bandno == 0) ? 0 : ((bandno - 1) / 3 + 1);
1419 orient = (bandno == 0) ? 0 : ((bandno - 1) % 3 + 1);
1420 level = tccp->numresolutions - 1 - resno;
1421 gain = (tccp->qmfbid == 0) ? 0 : ((orient == 0) ? 0 : (((orient == 1) ||
1422 (orient == 2)) ? 1 : 2));
1423 if (tccp->qntsty == J2K_CCP_QNTSTY_NOQNT) {
1426 OPJ_FLOAT64 norm = opj_dwt_norms_real[orient][level];
1427 stepsize = (1 << (gain)) / norm;
1429 opj_dwt_encode_stepsize((OPJ_INT32) floor(stepsize * 8192.0),
1430 (OPJ_INT32)(prec + gain), &tccp->stepsizes[bandno]);
1435 /* Determine maximum computed resolution level for inverse wavelet transform */
1437 static OPJ_UINT32 opj_dwt_max_resolution(opj_tcd_resolution_t* OPJ_RESTRICT r,
1444 if (mr < (w = (OPJ_UINT32)(r->x1 - r->x0))) {
1447 if (mr < (w = (OPJ_UINT32)(r->y1 - r->y0))) {
1458 OPJ_INT32 * OPJ_RESTRICT tiledp;
1461 } opj_dwt_decode_h_job_t;
1463 static void opj_dwt_decode_h_func(void* user_data, opj_tls_t* tls)
1466 opj_dwt_decode_h_job_t* job;
1469 job = (opj_dwt_decode_h_job_t*)user_data;
1470 for (j = job->min_j; j < job->max_j; j++) {
1471 opj_idwt53_h(&job->h, &job->tiledp[j * job->w]);
1474 opj_aligned_free(job->h.mem);
1482 OPJ_INT32 * OPJ_RESTRICT tiledp;
1485 } opj_dwt_decode_v_job_t;
1487 static void opj_dwt_decode_v_func(void* user_data, opj_tls_t* tls)
1490 opj_dwt_decode_v_job_t* job;
1493 job = (opj_dwt_decode_v_job_t*)user_data;
1494 for (j = job->min_j; j + PARALLEL_COLS_53 <= job->max_j;
1495 j += PARALLEL_COLS_53) {
1496 opj_idwt53_v(&job->v, &job->tiledp[j], (OPJ_SIZE_T)job->w,
1500 opj_idwt53_v(&job->v, &job->tiledp[j], (OPJ_SIZE_T)job->w,
1501 (OPJ_INT32)(job->max_j - j));
1503 opj_aligned_free(job->v.mem);
1509 /* Inverse wavelet transform in 2-D. */
1511 static OPJ_BOOL opj_dwt_decode_tile(opj_thread_pool_t* tp,
1512 opj_tcd_tilecomp_t* tilec, OPJ_UINT32 numres)
1517 opj_tcd_resolution_t* tr = tilec->resolutions;
1519 OPJ_UINT32 rw = (OPJ_UINT32)(tr->x1 -
1520 tr->x0); /* width of the resolution level computed */
1521 OPJ_UINT32 rh = (OPJ_UINT32)(tr->y1 -
1522 tr->y0); /* height of the resolution level computed */
1524 OPJ_UINT32 w = (OPJ_UINT32)(tilec->resolutions[tilec->minimum_num_resolutions -
1526 tilec->resolutions[tilec->minimum_num_resolutions - 1].x0);
1527 OPJ_SIZE_T h_mem_size;
1533 num_threads = opj_thread_pool_get_thread_count(tp);
1534 h_mem_size = opj_dwt_max_resolution(tr, numres);
1535 /* overflow check */
1536 if (h_mem_size > (SIZE_MAX / PARALLEL_COLS_53 / sizeof(OPJ_INT32))) {
1537 /* FIXME event manager error callback */
1540 /* We need PARALLEL_COLS_53 times the height of the array, */
1541 /* since for the vertical pass */
1542 /* we process PARALLEL_COLS_53 columns at a time */
1543 h_mem_size *= PARALLEL_COLS_53 * sizeof(OPJ_INT32);
1544 h.mem = (OPJ_INT32*)opj_aligned_32_malloc(h_mem_size);
1546 /* FIXME event manager error callback */
1553 OPJ_INT32 * OPJ_RESTRICT tiledp = tilec->data;
1557 h.sn = (OPJ_INT32)rw;
1558 v.sn = (OPJ_INT32)rh;
1560 rw = (OPJ_UINT32)(tr->x1 - tr->x0);
1561 rh = (OPJ_UINT32)(tr->y1 - tr->y0);
1563 h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
1566 if (num_threads <= 1 || rh <= 1) {
1567 for (j = 0; j < rh; ++j) {
1568 opj_idwt53_h(&h, &tiledp[(OPJ_SIZE_T)j * w]);
1571 OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
1574 if (rh < num_jobs) {
1577 step_j = (rh / num_jobs);
1579 for (j = 0; j < num_jobs; j++) {
1580 opj_dwt_decode_h_job_t* job;
1582 job = (opj_dwt_decode_h_job_t*) opj_malloc(sizeof(opj_dwt_decode_h_job_t));
1584 /* It would be nice to fallback to single thread case, but */
1585 /* unfortunately some jobs may be launched and have modified */
1586 /* tiledp, so it is not practical to recover from that error */
1587 /* FIXME event manager error callback */
1588 opj_thread_pool_wait_completion(tp, 0);
1589 opj_aligned_free(h.mem);
1595 job->tiledp = tiledp;
1596 job->min_j = j * step_j;
1597 job->max_j = (j + 1U) * step_j; /* this can overflow */
1598 if (j == (num_jobs - 1U)) { /* this will take care of the overflow */
1601 job->h.mem = (OPJ_INT32*)opj_aligned_32_malloc(h_mem_size);
1603 /* FIXME event manager error callback */
1604 opj_thread_pool_wait_completion(tp, 0);
1606 opj_aligned_free(h.mem);
1609 opj_thread_pool_submit_job(tp, opj_dwt_decode_h_func, job);
1611 opj_thread_pool_wait_completion(tp, 0);
1614 v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
1617 if (num_threads <= 1 || rw <= 1) {
1618 for (j = 0; j + PARALLEL_COLS_53 <= rw;
1619 j += PARALLEL_COLS_53) {
1620 opj_idwt53_v(&v, &tiledp[j], (OPJ_SIZE_T)w, PARALLEL_COLS_53);
1623 opj_idwt53_v(&v, &tiledp[j], (OPJ_SIZE_T)w, (OPJ_INT32)(rw - j));
1626 OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
1629 if (rw < num_jobs) {
1632 step_j = (rw / num_jobs);
1634 for (j = 0; j < num_jobs; j++) {
1635 opj_dwt_decode_v_job_t* job;
1637 job = (opj_dwt_decode_v_job_t*) opj_malloc(sizeof(opj_dwt_decode_v_job_t));
1639 /* It would be nice to fallback to single thread case, but */
1640 /* unfortunately some jobs may be launched and have modified */
1641 /* tiledp, so it is not practical to recover from that error */
1642 /* FIXME event manager error callback */
1643 opj_thread_pool_wait_completion(tp, 0);
1644 opj_aligned_free(v.mem);
1650 job->tiledp = tiledp;
1651 job->min_j = j * step_j;
1652 job->max_j = (j + 1U) * step_j; /* this can overflow */
1653 if (j == (num_jobs - 1U)) { /* this will take care of the overflow */
1656 job->v.mem = (OPJ_INT32*)opj_aligned_32_malloc(h_mem_size);
1658 /* FIXME event manager error callback */
1659 opj_thread_pool_wait_completion(tp, 0);
1661 opj_aligned_free(v.mem);
1664 opj_thread_pool_submit_job(tp, opj_dwt_decode_v_func, job);
1666 opj_thread_pool_wait_completion(tp, 0);
1669 opj_aligned_free(h.mem);
1673 static void opj_dwt_interleave_partial_h(OPJ_INT32 *dest,
1675 opj_sparse_array_int32_t* sa,
1678 OPJ_UINT32 win_l_x0,
1679 OPJ_UINT32 win_l_x1,
1680 OPJ_UINT32 win_h_x0,
1681 OPJ_UINT32 win_h_x1)
1684 ret = opj_sparse_array_int32_read(sa,
1686 win_l_x1, sa_line + 1,
1687 dest + cas + 2 * win_l_x0,
1690 ret = opj_sparse_array_int32_read(sa,
1691 sn + win_h_x0, sa_line,
1692 sn + win_h_x1, sa_line + 1,
1693 dest + 1 - cas + 2 * win_h_x0,
1700 static void opj_dwt_interleave_partial_v(OPJ_INT32 *dest,
1702 opj_sparse_array_int32_t* sa,
1706 OPJ_UINT32 win_l_y0,
1707 OPJ_UINT32 win_l_y1,
1708 OPJ_UINT32 win_h_y0,
1709 OPJ_UINT32 win_h_y1)
1712 ret = opj_sparse_array_int32_read(sa,
1714 sa_col + nb_cols, win_l_y1,
1715 dest + cas * 4 + 2 * 4 * win_l_y0,
1716 1, 2 * 4, OPJ_TRUE);
1718 ret = opj_sparse_array_int32_read(sa,
1719 sa_col, sn + win_h_y0,
1720 sa_col + nb_cols, sn + win_h_y1,
1721 dest + (1 - cas) * 4 + 2 * 4 * win_h_y0,
1722 1, 2 * 4, OPJ_TRUE);
1727 static void opj_dwt_decode_partial_1(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
1737 if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
1739 /* Naive version is :
1740 for (i = win_l_x0; i < i_max; i++) {
1741 OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
1743 for (i = win_h_x0; i < win_h_x1; i++) {
1744 OPJ_D(i) += (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
1746 but the compiler doesn't manage to unroll it to avoid bound
1747 checking in OPJ_S_ and OPJ_D_ macros
1754 /* Left-most case */
1755 OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
1762 for (; i < i_max; i++) {
1763 /* No bound checking */
1764 OPJ_S(i) -= (OPJ_D(i - 1) + OPJ_D(i) + 2) >> 2;
1766 for (; i < win_l_x1; i++) {
1767 /* Right-most case */
1768 OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
1774 OPJ_INT32 i_max = win_h_x1;
1778 for (; i < i_max; i++) {
1779 /* No bound checking */
1780 OPJ_D(i) += (OPJ_S(i) + OPJ_S(i + 1)) >> 1;
1782 for (; i < win_h_x1; i++) {
1783 /* Right-most case */
1784 OPJ_D(i) += (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
1789 if (!sn && dn == 1) { /* NEW : CASE ONE ELEMENT */
1792 for (i = win_l_x0; i < win_l_x1; i++) {
1793 OPJ_D(i) -= (OPJ_SS_(i) + OPJ_SS_(i + 1) + 2) >> 2;
1795 for (i = win_h_x0; i < win_h_x1; i++) {
1796 OPJ_S(i) += (OPJ_DD_(i) + OPJ_DD_(i - 1)) >> 1;
1802 #define OPJ_S_off(i,off) a[(OPJ_UINT32)(i)*2*4+off]
1803 #define OPJ_D_off(i,off) a[(1+(OPJ_UINT32)(i)*2)*4+off]
1804 #define OPJ_S__off(i,off) ((i)<0?OPJ_S_off(0,off):((i)>=sn?OPJ_S_off(sn-1,off):OPJ_S_off(i,off)))
1805 #define OPJ_D__off(i,off) ((i)<0?OPJ_D_off(0,off):((i)>=dn?OPJ_D_off(dn-1,off):OPJ_D_off(i,off)))
1806 #define OPJ_SS__off(i,off) ((i)<0?OPJ_S_off(0,off):((i)>=dn?OPJ_S_off(dn-1,off):OPJ_S_off(i,off)))
1807 #define OPJ_DD__off(i,off) ((i)<0?OPJ_D_off(0,off):((i)>=sn?OPJ_D_off(sn-1,off):OPJ_D_off(i,off)))
1809 static void opj_dwt_decode_partial_1_parallel(OPJ_INT32 *a,
1811 OPJ_INT32 dn, OPJ_INT32 sn,
1824 if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
1826 /* Naive version is :
1827 for (i = win_l_x0; i < i_max; i++) {
1828 OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
1830 for (i = win_h_x0; i < win_h_x1; i++) {
1831 OPJ_D(i) += (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
1833 but the compiler doesn't manage to unroll it to avoid bound
1834 checking in OPJ_S_ and OPJ_D_ macros
1841 /* Left-most case */
1842 for (off = 0; off < 4; off++) {
1843 OPJ_S_off(i, off) -= (OPJ_D__off(i - 1, off) + OPJ_D__off(i, off) + 2) >> 2;
1853 if (i + 1 < i_max) {
1854 const __m128i two = _mm_set1_epi32(2);
1855 __m128i Dm1 = _mm_load_si128((__m128i * const)(a + 4 + (i - 1) * 8));
1856 for (; i + 1 < i_max; i += 2) {
1857 /* No bound checking */
1858 __m128i S = _mm_load_si128((__m128i * const)(a + i * 8));
1859 __m128i D = _mm_load_si128((__m128i * const)(a + 4 + i * 8));
1860 __m128i S1 = _mm_load_si128((__m128i * const)(a + (i + 1) * 8));
1861 __m128i D1 = _mm_load_si128((__m128i * const)(a + 4 + (i + 1) * 8));
1862 S = _mm_sub_epi32(S,
1863 _mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(Dm1, D), two), 2));
1864 S1 = _mm_sub_epi32(S1,
1865 _mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(D, D1), two), 2));
1866 _mm_store_si128((__m128i*)(a + i * 8), S);
1867 _mm_store_si128((__m128i*)(a + (i + 1) * 8), S1);
1873 for (; i < i_max; i++) {
1874 /* No bound checking */
1875 for (off = 0; off < 4; off++) {
1876 OPJ_S_off(i, off) -= (OPJ_D_off(i - 1, off) + OPJ_D_off(i, off) + 2) >> 2;
1879 for (; i < win_l_x1; i++) {
1880 /* Right-most case */
1881 for (off = 0; off < 4; off++) {
1882 OPJ_S_off(i, off) -= (OPJ_D__off(i - 1, off) + OPJ_D__off(i, off) + 2) >> 2;
1889 OPJ_INT32 i_max = win_h_x1;
1895 if (i + 1 < i_max) {
1896 __m128i S = _mm_load_si128((__m128i * const)(a + i * 8));
1897 for (; i + 1 < i_max; i += 2) {
1898 /* No bound checking */
1899 __m128i D = _mm_load_si128((__m128i * const)(a + 4 + i * 8));
1900 __m128i S1 = _mm_load_si128((__m128i * const)(a + (i + 1) * 8));
1901 __m128i D1 = _mm_load_si128((__m128i * const)(a + 4 + (i + 1) * 8));
1902 __m128i S2 = _mm_load_si128((__m128i * const)(a + (i + 2) * 8));
1903 D = _mm_add_epi32(D, _mm_srai_epi32(_mm_add_epi32(S, S1), 1));
1904 D1 = _mm_add_epi32(D1, _mm_srai_epi32(_mm_add_epi32(S1, S2), 1));
1905 _mm_store_si128((__m128i*)(a + 4 + i * 8), D);
1906 _mm_store_si128((__m128i*)(a + 4 + (i + 1) * 8), D1);
1912 for (; i < i_max; i++) {
1913 /* No bound checking */
1914 for (off = 0; off < 4; off++) {
1915 OPJ_D_off(i, off) += (OPJ_S_off(i, off) + OPJ_S_off(i + 1, off)) >> 1;
1918 for (; i < win_h_x1; i++) {
1919 /* Right-most case */
1920 for (off = 0; off < 4; off++) {
1921 OPJ_D_off(i, off) += (OPJ_S__off(i, off) + OPJ_S__off(i + 1, off)) >> 1;
1927 if (!sn && dn == 1) { /* NEW : CASE ONE ELEMENT */
1928 for (off = 0; off < 4; off++) {
1929 OPJ_S_off(0, off) /= 2;
1932 for (i = win_l_x0; i < win_l_x1; i++) {
1933 for (off = 0; off < 4; off++) {
1934 OPJ_D_off(i, off) -= (OPJ_SS__off(i, off) + OPJ_SS__off(i + 1, off) + 2) >> 2;
1937 for (i = win_h_x0; i < win_h_x1; i++) {
1938 for (off = 0; off < 4; off++) {
1939 OPJ_S_off(i, off) += (OPJ_DD__off(i, off) + OPJ_DD__off(i - 1, off)) >> 1;
1946 static void opj_dwt_get_band_coordinates(opj_tcd_tilecomp_t* tilec,
1958 /* Compute number of decomposition for this band. See table F-1 */
1959 OPJ_UINT32 nb = (resno == 0) ?
1960 tilec->numresolutions - 1 :
1961 tilec->numresolutions - resno;
1962 /* Map above tile-based coordinates to sub-band-based coordinates per */
1963 /* equation B-15 of the standard */
1964 OPJ_UINT32 x0b = bandno & 1;
1965 OPJ_UINT32 y0b = bandno >> 1;
1967 *tbx0 = (nb == 0) ? tcx0 :
1968 (tcx0 <= (1U << (nb - 1)) * x0b) ? 0 :
1969 opj_uint_ceildivpow2(tcx0 - (1U << (nb - 1)) * x0b, nb);
1972 *tby0 = (nb == 0) ? tcy0 :
1973 (tcy0 <= (1U << (nb - 1)) * y0b) ? 0 :
1974 opj_uint_ceildivpow2(tcy0 - (1U << (nb - 1)) * y0b, nb);
1977 *tbx1 = (nb == 0) ? tcx1 :
1978 (tcx1 <= (1U << (nb - 1)) * x0b) ? 0 :
1979 opj_uint_ceildivpow2(tcx1 - (1U << (nb - 1)) * x0b, nb);
1982 *tby1 = (nb == 0) ? tcy1 :
1983 (tcy1 <= (1U << (nb - 1)) * y0b) ? 0 :
1984 opj_uint_ceildivpow2(tcy1 - (1U << (nb - 1)) * y0b, nb);
1988 static void opj_dwt_segment_grow(OPJ_UINT32 filter_width,
1989 OPJ_UINT32 max_size,
1993 *start = opj_uint_subs(*start, filter_width);
1994 *end = opj_uint_adds(*end, filter_width);
1995 *end = opj_uint_min(*end, max_size);
1999 static opj_sparse_array_int32_t* opj_dwt_init_sparse_array(
2000 opj_tcd_tilecomp_t* tilec,
2003 opj_tcd_resolution_t* tr_max = &(tilec->resolutions[numres - 1]);
2004 OPJ_UINT32 w = (OPJ_UINT32)(tr_max->x1 - tr_max->x0);
2005 OPJ_UINT32 h = (OPJ_UINT32)(tr_max->y1 - tr_max->y0);
2006 OPJ_UINT32 resno, bandno, precno, cblkno;
2007 opj_sparse_array_int32_t* sa = opj_sparse_array_int32_create(
2008 w, h, opj_uint_min(w, 64), opj_uint_min(h, 64));
2013 for (resno = 0; resno < numres; ++resno) {
2014 opj_tcd_resolution_t* res = &tilec->resolutions[resno];
2016 for (bandno = 0; bandno < res->numbands; ++bandno) {
2017 opj_tcd_band_t* band = &res->bands[bandno];
2019 for (precno = 0; precno < res->pw * res->ph; ++precno) {
2020 opj_tcd_precinct_t* precinct = &band->precincts[precno];
2021 for (cblkno = 0; cblkno < precinct->cw * precinct->ch; ++cblkno) {
2022 opj_tcd_cblk_dec_t* cblk = &precinct->cblks.dec[cblkno];
2023 if (cblk->decoded_data != NULL) {
2024 OPJ_UINT32 x = (OPJ_UINT32)(cblk->x0 - band->x0);
2025 OPJ_UINT32 y = (OPJ_UINT32)(cblk->y0 - band->y0);
2026 OPJ_UINT32 cblk_w = (OPJ_UINT32)(cblk->x1 - cblk->x0);
2027 OPJ_UINT32 cblk_h = (OPJ_UINT32)(cblk->y1 - cblk->y0);
2029 if (band->bandno & 1) {
2030 opj_tcd_resolution_t* pres = &tilec->resolutions[resno - 1];
2031 x += (OPJ_UINT32)(pres->x1 - pres->x0);
2033 if (band->bandno & 2) {
2034 opj_tcd_resolution_t* pres = &tilec->resolutions[resno - 1];
2035 y += (OPJ_UINT32)(pres->y1 - pres->y0);
2038 if (!opj_sparse_array_int32_write(sa, x, y,
2039 x + cblk_w, y + cblk_h,
2041 1, cblk_w, OPJ_TRUE)) {
2042 opj_sparse_array_int32_free(sa);
2055 static OPJ_BOOL opj_dwt_decode_partial_tile(
2056 opj_tcd_tilecomp_t* tilec,
2059 opj_sparse_array_int32_t* sa;
2063 /* This value matches the maximum left/right extension given in tables */
2064 /* F.2 and F.3 of the standard. */
2065 const OPJ_UINT32 filter_width = 2U;
2067 opj_tcd_resolution_t* tr = tilec->resolutions;
2068 opj_tcd_resolution_t* tr_max = &(tilec->resolutions[numres - 1]);
2070 OPJ_UINT32 rw = (OPJ_UINT32)(tr->x1 -
2071 tr->x0); /* width of the resolution level computed */
2072 OPJ_UINT32 rh = (OPJ_UINT32)(tr->y1 -
2073 tr->y0); /* height of the resolution level computed */
2075 OPJ_SIZE_T h_mem_size;
2077 /* Compute the intersection of the area of interest, expressed in tile coordinates */
2078 /* with the tile coordinates */
2079 OPJ_UINT32 win_tcx0 = tilec->win_x0;
2080 OPJ_UINT32 win_tcy0 = tilec->win_y0;
2081 OPJ_UINT32 win_tcx1 = tilec->win_x1;
2082 OPJ_UINT32 win_tcy1 = tilec->win_y1;
2084 if (tr_max->x0 == tr_max->x1 || tr_max->y0 == tr_max->y1) {
2088 sa = opj_dwt_init_sparse_array(tilec, numres);
2094 OPJ_BOOL ret = opj_sparse_array_int32_read(sa,
2095 tr_max->win_x0 - (OPJ_UINT32)tr_max->x0,
2096 tr_max->win_y0 - (OPJ_UINT32)tr_max->y0,
2097 tr_max->win_x1 - (OPJ_UINT32)tr_max->x0,
2098 tr_max->win_y1 - (OPJ_UINT32)tr_max->y0,
2100 1, tr_max->win_x1 - tr_max->win_x0,
2104 opj_sparse_array_int32_free(sa);
2107 h_mem_size = opj_dwt_max_resolution(tr, numres);
2108 /* overflow check */
2109 /* in vertical pass, we process 4 columns at a time */
2110 if (h_mem_size > (SIZE_MAX / (4 * sizeof(OPJ_INT32)))) {
2111 /* FIXME event manager error callback */
2112 opj_sparse_array_int32_free(sa);
2116 h_mem_size *= 4 * sizeof(OPJ_INT32);
2117 h.mem = (OPJ_INT32*)opj_aligned_32_malloc(h_mem_size);
2119 /* FIXME event manager error callback */
2120 opj_sparse_array_int32_free(sa);
2126 for (resno = 1; resno < numres; resno ++) {
2128 /* Window of interest subband-based coordinates */
2129 OPJ_UINT32 win_ll_x0, win_ll_y0, win_ll_x1, win_ll_y1;
2130 OPJ_UINT32 win_hl_x0, win_hl_x1;
2131 OPJ_UINT32 win_lh_y0, win_lh_y1;
2132 /* Window of interest tile-resolution-based coordinates */
2133 OPJ_UINT32 win_tr_x0, win_tr_x1, win_tr_y0, win_tr_y1;
2134 /* Tile-resolution subband-based coordinates */
2135 OPJ_UINT32 tr_ll_x0, tr_ll_y0, tr_hl_x0, tr_lh_y0;
2139 h.sn = (OPJ_INT32)rw;
2140 v.sn = (OPJ_INT32)rh;
2142 rw = (OPJ_UINT32)(tr->x1 - tr->x0);
2143 rh = (OPJ_UINT32)(tr->y1 - tr->y0);
2145 h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
2148 v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
2151 /* Get the subband coordinates for the window of interest */
2153 opj_dwt_get_band_coordinates(tilec, resno, 0,
2154 win_tcx0, win_tcy0, win_tcx1, win_tcy1,
2155 &win_ll_x0, &win_ll_y0,
2156 &win_ll_x1, &win_ll_y1);
2159 opj_dwt_get_band_coordinates(tilec, resno, 1,
2160 win_tcx0, win_tcy0, win_tcx1, win_tcy1,
2161 &win_hl_x0, NULL, &win_hl_x1, NULL);
2164 opj_dwt_get_band_coordinates(tilec, resno, 2,
2165 win_tcx0, win_tcy0, win_tcx1, win_tcy1,
2166 NULL, &win_lh_y0, NULL, &win_lh_y1);
2168 /* Beware: band index for non-LL0 resolution are 0=HL, 1=LH and 2=HH */
2169 tr_ll_x0 = (OPJ_UINT32)tr->bands[1].x0;
2170 tr_ll_y0 = (OPJ_UINT32)tr->bands[0].y0;
2171 tr_hl_x0 = (OPJ_UINT32)tr->bands[0].x0;
2172 tr_lh_y0 = (OPJ_UINT32)tr->bands[1].y0;
2174 /* Subtract the origin of the bands for this tile, to the subwindow */
2175 /* of interest band coordinates, so as to get them relative to the */
2177 win_ll_x0 = opj_uint_subs(win_ll_x0, tr_ll_x0);
2178 win_ll_y0 = opj_uint_subs(win_ll_y0, tr_ll_y0);
2179 win_ll_x1 = opj_uint_subs(win_ll_x1, tr_ll_x0);
2180 win_ll_y1 = opj_uint_subs(win_ll_y1, tr_ll_y0);
2181 win_hl_x0 = opj_uint_subs(win_hl_x0, tr_hl_x0);
2182 win_hl_x1 = opj_uint_subs(win_hl_x1, tr_hl_x0);
2183 win_lh_y0 = opj_uint_subs(win_lh_y0, tr_lh_y0);
2184 win_lh_y1 = opj_uint_subs(win_lh_y1, tr_lh_y0);
2186 opj_dwt_segment_grow(filter_width, (OPJ_UINT32)h.sn, &win_ll_x0, &win_ll_x1);
2187 opj_dwt_segment_grow(filter_width, (OPJ_UINT32)h.dn, &win_hl_x0, &win_hl_x1);
2189 opj_dwt_segment_grow(filter_width, (OPJ_UINT32)v.sn, &win_ll_y0, &win_ll_y1);
2190 opj_dwt_segment_grow(filter_width, (OPJ_UINT32)v.dn, &win_lh_y0, &win_lh_y1);
2192 /* Compute the tile-resolution-based coordinates for the window of interest */
2194 win_tr_x0 = opj_uint_min(2 * win_ll_x0, 2 * win_hl_x0 + 1);
2195 win_tr_x1 = opj_uint_min(opj_uint_max(2 * win_ll_x1, 2 * win_hl_x1 + 1), rw);
2197 win_tr_x0 = opj_uint_min(2 * win_hl_x0, 2 * win_ll_x0 + 1);
2198 win_tr_x1 = opj_uint_min(opj_uint_max(2 * win_hl_x1, 2 * win_ll_x1 + 1), rw);
2202 win_tr_y0 = opj_uint_min(2 * win_ll_y0, 2 * win_lh_y0 + 1);
2203 win_tr_y1 = opj_uint_min(opj_uint_max(2 * win_ll_y1, 2 * win_lh_y1 + 1), rh);
2205 win_tr_y0 = opj_uint_min(2 * win_lh_y0, 2 * win_ll_y0 + 1);
2206 win_tr_y1 = opj_uint_min(opj_uint_max(2 * win_lh_y1, 2 * win_ll_y1 + 1), rh);
2209 for (j = 0; j < rh; ++j) {
2210 if ((j >= win_ll_y0 && j < win_ll_y1) ||
2211 (j >= win_lh_y0 + (OPJ_UINT32)v.sn && j < win_lh_y1 + (OPJ_UINT32)v.sn)) {
2213 /* Avoids dwt.c:1584:44 (in opj_dwt_decode_partial_1): runtime error: */
2214 /* signed integer overflow: -1094795586 + -1094795586 cannot be represented in type 'int' */
2215 /* on opj_decompress -i ../../openjpeg/MAPA.jp2 -o out.tif -d 0,0,256,256 */
2216 /* This is less extreme than memsetting the whole buffer to 0 */
2217 /* although we could potentially do better with better handling of edge conditions */
2218 if (win_tr_x1 >= 1 && win_tr_x1 < rw) {
2219 h.mem[win_tr_x1 - 1] = 0;
2221 if (win_tr_x1 < rw) {
2222 h.mem[win_tr_x1] = 0;
2225 opj_dwt_interleave_partial_h(h.mem,
2234 opj_dwt_decode_partial_1(h.mem, h.dn, h.sn, h.cas,
2235 (OPJ_INT32)win_ll_x0,
2236 (OPJ_INT32)win_ll_x1,
2237 (OPJ_INT32)win_hl_x0,
2238 (OPJ_INT32)win_hl_x1);
2239 if (!opj_sparse_array_int32_write(sa,
2244 /* FIXME event manager error callback */
2245 opj_sparse_array_int32_free(sa);
2246 opj_aligned_free(h.mem);
2252 for (i = win_tr_x0; i < win_tr_x1;) {
2253 OPJ_UINT32 nb_cols = opj_uint_min(4U, win_tr_x1 - i);
2254 opj_dwt_interleave_partial_v(v.mem,
2264 opj_dwt_decode_partial_1_parallel(v.mem, nb_cols, v.dn, v.sn, v.cas,
2265 (OPJ_INT32)win_ll_y0,
2266 (OPJ_INT32)win_ll_y1,
2267 (OPJ_INT32)win_lh_y0,
2268 (OPJ_INT32)win_lh_y1);
2269 if (!opj_sparse_array_int32_write(sa,
2271 i + nb_cols, win_tr_y1,
2272 v.mem + 4 * win_tr_y0,
2274 /* FIXME event manager error callback */
2275 opj_sparse_array_int32_free(sa);
2276 opj_aligned_free(h.mem);
2283 opj_aligned_free(h.mem);
2286 OPJ_BOOL ret = opj_sparse_array_int32_read(sa,
2287 tr_max->win_x0 - (OPJ_UINT32)tr_max->x0,
2288 tr_max->win_y0 - (OPJ_UINT32)tr_max->y0,
2289 tr_max->win_x1 - (OPJ_UINT32)tr_max->x0,
2290 tr_max->win_y1 - (OPJ_UINT32)tr_max->y0,
2292 1, tr_max->win_x1 - tr_max->win_x0,
2297 opj_sparse_array_int32_free(sa);
2301 static void opj_v8dwt_interleave_h(opj_v8dwt_t* OPJ_RESTRICT dwt,
2302 OPJ_FLOAT32* OPJ_RESTRICT a,
2304 OPJ_UINT32 remaining_height)
2306 OPJ_FLOAT32* OPJ_RESTRICT bi = (OPJ_FLOAT32*)(dwt->wavelet + dwt->cas);
2308 OPJ_UINT32 x0 = dwt->win_l_x0;
2309 OPJ_UINT32 x1 = dwt->win_l_x1;
2311 for (k = 0; k < 2; ++k) {
2312 if (remaining_height >= NB_ELTS_V8 && ((OPJ_SIZE_T) a & 0x0f) == 0 &&
2313 ((OPJ_SIZE_T) bi & 0x0f) == 0) {
2314 /* Fast code path */
2315 for (i = x0; i < x1; ++i) {
2317 OPJ_FLOAT32* OPJ_RESTRICT dst = bi + i * 2 * NB_ELTS_V8;
2335 /* Slow code path */
2336 for (i = x0; i < x1; ++i) {
2338 OPJ_FLOAT32* OPJ_RESTRICT dst = bi + i * 2 * NB_ELTS_V8;
2341 if (remaining_height == 1) {
2346 if (remaining_height == 2) {
2351 if (remaining_height == 3) {
2356 if (remaining_height == 4) {
2361 if (remaining_height == 5) {
2366 if (remaining_height == 6) {
2371 if (remaining_height == 7) {
2378 bi = (OPJ_FLOAT32*)(dwt->wavelet + 1 - dwt->cas);
2385 static void opj_v8dwt_interleave_partial_h(opj_v8dwt_t* dwt,
2386 opj_sparse_array_int32_t* sa,
2388 OPJ_UINT32 remaining_height)
2391 for (i = 0; i < remaining_height; i++) {
2393 ret = opj_sparse_array_int32_read(sa,
2394 dwt->win_l_x0, sa_line + i,
2395 dwt->win_l_x1, sa_line + i + 1,
2396 /* Nasty cast from float* to int32* */
2397 (OPJ_INT32*)(dwt->wavelet + dwt->cas + 2 * dwt->win_l_x0) + i,
2398 2 * NB_ELTS_V8, 0, OPJ_TRUE);
2400 ret = opj_sparse_array_int32_read(sa,
2401 (OPJ_UINT32)dwt->sn + dwt->win_h_x0, sa_line + i,
2402 (OPJ_UINT32)dwt->sn + dwt->win_h_x1, sa_line + i + 1,
2403 /* Nasty cast from float* to int32* */
2404 (OPJ_INT32*)(dwt->wavelet + 1 - dwt->cas + 2 * dwt->win_h_x0) + i,
2405 2 * NB_ELTS_V8, 0, OPJ_TRUE);
2411 static INLINE void opj_v8dwt_interleave_v(opj_v8dwt_t* OPJ_RESTRICT dwt,
2412 OPJ_FLOAT32* OPJ_RESTRICT a,
2414 OPJ_UINT32 nb_elts_read)
2416 opj_v8_t* OPJ_RESTRICT bi = dwt->wavelet + dwt->cas;
2419 for (i = dwt->win_l_x0; i < dwt->win_l_x1; ++i) {
2420 memcpy(&bi[i * 2], &a[i * (OPJ_SIZE_T)width],
2421 (OPJ_SIZE_T)nb_elts_read * sizeof(OPJ_FLOAT32));
2424 a += (OPJ_UINT32)dwt->sn * (OPJ_SIZE_T)width;
2425 bi = dwt->wavelet + 1 - dwt->cas;
2427 for (i = dwt->win_h_x0; i < dwt->win_h_x1; ++i) {
2428 memcpy(&bi[i * 2], &a[i * (OPJ_SIZE_T)width],
2429 (OPJ_SIZE_T)nb_elts_read * sizeof(OPJ_FLOAT32));
2433 static void opj_v8dwt_interleave_partial_v(opj_v8dwt_t* OPJ_RESTRICT dwt,
2434 opj_sparse_array_int32_t* sa,
2436 OPJ_UINT32 nb_elts_read)
2439 ret = opj_sparse_array_int32_read(sa,
2440 sa_col, dwt->win_l_x0,
2441 sa_col + nb_elts_read, dwt->win_l_x1,
2442 (OPJ_INT32*)(dwt->wavelet + dwt->cas + 2 * dwt->win_l_x0),
2443 1, 2 * NB_ELTS_V8, OPJ_TRUE);
2445 ret = opj_sparse_array_int32_read(sa,
2446 sa_col, (OPJ_UINT32)dwt->sn + dwt->win_h_x0,
2447 sa_col + nb_elts_read, (OPJ_UINT32)dwt->sn + dwt->win_h_x1,
2448 (OPJ_INT32*)(dwt->wavelet + 1 - dwt->cas + 2 * dwt->win_h_x0),
2449 1, 2 * NB_ELTS_V8, OPJ_TRUE);
2456 static void opj_v8dwt_decode_step1_sse(opj_v8_t* w,
2461 __m128* OPJ_RESTRICT vw = (__m128*) w;
2462 OPJ_UINT32 i = start;
2463 /* To be adapted if NB_ELTS_V8 changes */
2465 /* Note: attempt at loop unrolling x2 doesn't help */
2466 for (; i < end; ++i, vw += 4) {
2467 vw[0] = _mm_mul_ps(vw[0], c);
2468 vw[1] = _mm_mul_ps(vw[1], c);
2472 static void opj_v8dwt_decode_step2_sse(opj_v8_t* l, opj_v8_t* w,
2478 __m128* OPJ_RESTRICT vl = (__m128*) l;
2479 __m128* OPJ_RESTRICT vw = (__m128*) w;
2480 /* To be adapted if NB_ELTS_V8 changes */
2482 OPJ_UINT32 imax = opj_uint_min(end, m);
2485 vw[-2] = _mm_add_ps(vw[-2], _mm_mul_ps(_mm_add_ps(vl[0], vw[0]), c));
2486 vw[-1] = _mm_add_ps(vw[-1], _mm_mul_ps(_mm_add_ps(vl[1], vw[1]), c));
2495 /* Note: attempt at loop unrolling x2 doesn't help */
2496 for (; i < imax; ++i) {
2497 vw[-2] = _mm_add_ps(vw[-2], _mm_mul_ps(_mm_add_ps(vw[-4], vw[0]), c));
2498 vw[-1] = _mm_add_ps(vw[-1], _mm_mul_ps(_mm_add_ps(vw[-3], vw[1]), c));
2502 assert(m + 1 == end);
2503 c = _mm_add_ps(c, c);
2504 vw[-2] = _mm_add_ps(vw[-2], _mm_mul_ps(c, vw[-4]));
2505 vw[-1] = _mm_add_ps(vw[-1], _mm_mul_ps(c, vw[-3]));
2511 static void opj_v8dwt_decode_step1(opj_v8_t* w,
2514 const OPJ_FLOAT32 c)
2516 OPJ_FLOAT32* OPJ_RESTRICT fw = (OPJ_FLOAT32*) w;
2518 /* To be adapted if NB_ELTS_V8 changes */
2519 for (i = start; i < end; ++i) {
2520 fw[i * 2 * 8 ] = fw[i * 2 * 8 ] * c;
2521 fw[i * 2 * 8 + 1] = fw[i * 2 * 8 + 1] * c;
2522 fw[i * 2 * 8 + 2] = fw[i * 2 * 8 + 2] * c;
2523 fw[i * 2 * 8 + 3] = fw[i * 2 * 8 + 3] * c;
2524 fw[i * 2 * 8 + 4] = fw[i * 2 * 8 + 4] * c;
2525 fw[i * 2 * 8 + 5] = fw[i * 2 * 8 + 5] * c;
2526 fw[i * 2 * 8 + 6] = fw[i * 2 * 8 + 6] * c;
2527 fw[i * 2 * 8 + 7] = fw[i * 2 * 8 + 7] * c;
2531 static void opj_v8dwt_decode_step2(opj_v8_t* l, opj_v8_t* w,
2537 OPJ_FLOAT32* fl = (OPJ_FLOAT32*) l;
2538 OPJ_FLOAT32* fw = (OPJ_FLOAT32*) w;
2540 OPJ_UINT32 imax = opj_uint_min(end, m);
2542 fw += 2 * NB_ELTS_V8 * start;
2543 fl = fw - 2 * NB_ELTS_V8;
2545 /* To be adapted if NB_ELTS_V8 changes */
2546 for (i = start; i < imax; ++i) {
2547 fw[-8] = fw[-8] + ((fl[0] + fw[0]) * c);
2548 fw[-7] = fw[-7] + ((fl[1] + fw[1]) * c);
2549 fw[-6] = fw[-6] + ((fl[2] + fw[2]) * c);
2550 fw[-5] = fw[-5] + ((fl[3] + fw[3]) * c);
2551 fw[-4] = fw[-4] + ((fl[4] + fw[4]) * c);
2552 fw[-3] = fw[-3] + ((fl[5] + fw[5]) * c);
2553 fw[-2] = fw[-2] + ((fl[6] + fw[6]) * c);
2554 fw[-1] = fw[-1] + ((fl[7] + fw[7]) * c);
2556 fw += 2 * NB_ELTS_V8;
2559 assert(m + 1 == end);
2561 fw[-8] = fw[-8] + fl[0] * c;
2562 fw[-7] = fw[-7] + fl[1] * c;
2563 fw[-6] = fw[-6] + fl[2] * c;
2564 fw[-5] = fw[-5] + fl[3] * c;
2565 fw[-4] = fw[-4] + fl[4] * c;
2566 fw[-3] = fw[-3] + fl[5] * c;
2567 fw[-2] = fw[-2] + fl[6] * c;
2568 fw[-1] = fw[-1] + fl[7] * c;
2575 /* Inverse 9-7 wavelet transform in 1-D. */
2577 static void opj_v8dwt_decode(opj_v8dwt_t* OPJ_RESTRICT dwt)
2580 /* BUG_WEIRD_TWO_INVK (look for this identifier in tcd.c) */
2581 /* Historic value for 2 / opj_invK */
2582 /* Normally, we should use invK, but if we do so, we have failures in the */
2583 /* conformance test, due to MSE and peak errors significantly higher than */
2584 /* accepted value */
2585 /* Due to using two_invK instead of invK, we have to compensate in tcd.c */
2586 /* the computation of the stepsize for the non LL subbands */
2587 const float two_invK = 1.625732422f;
2588 if (dwt->cas == 0) {
2589 if (!((dwt->dn > 0) || (dwt->sn > 1))) {
2595 if (!((dwt->sn > 0) || (dwt->dn > 1))) {
2602 opj_v8dwt_decode_step1_sse(dwt->wavelet + a, dwt->win_l_x0, dwt->win_l_x1,
2603 _mm_set1_ps(opj_K));
2604 opj_v8dwt_decode_step1_sse(dwt->wavelet + b, dwt->win_h_x0, dwt->win_h_x1,
2605 _mm_set1_ps(two_invK));
2606 opj_v8dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1,
2607 dwt->win_l_x0, dwt->win_l_x1,
2608 (OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a),
2609 _mm_set1_ps(-opj_dwt_delta));
2610 opj_v8dwt_decode_step2_sse(dwt->wavelet + a, dwt->wavelet + b + 1,
2611 dwt->win_h_x0, dwt->win_h_x1,
2612 (OPJ_UINT32)opj_int_min(dwt->dn, dwt->sn - b),
2613 _mm_set1_ps(-opj_dwt_gamma));
2614 opj_v8dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1,
2615 dwt->win_l_x0, dwt->win_l_x1,
2616 (OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a),
2617 _mm_set1_ps(-opj_dwt_beta));
2618 opj_v8dwt_decode_step2_sse(dwt->wavelet + a, dwt->wavelet + b + 1,
2619 dwt->win_h_x0, dwt->win_h_x1,
2620 (OPJ_UINT32)opj_int_min(dwt->dn, dwt->sn - b),
2621 _mm_set1_ps(-opj_dwt_alpha));
2623 opj_v8dwt_decode_step1(dwt->wavelet + a, dwt->win_l_x0, dwt->win_l_x1,
2625 opj_v8dwt_decode_step1(dwt->wavelet + b, dwt->win_h_x0, dwt->win_h_x1,
2627 opj_v8dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1,
2628 dwt->win_l_x0, dwt->win_l_x1,
2629 (OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a),
2631 opj_v8dwt_decode_step2(dwt->wavelet + a, dwt->wavelet + b + 1,
2632 dwt->win_h_x0, dwt->win_h_x1,
2633 (OPJ_UINT32)opj_int_min(dwt->dn, dwt->sn - b),
2635 opj_v8dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1,
2636 dwt->win_l_x0, dwt->win_l_x1,
2637 (OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a),
2639 opj_v8dwt_decode_step2(dwt->wavelet + a, dwt->wavelet + b + 1,
2640 dwt->win_h_x0, dwt->win_h_x1,
2641 (OPJ_UINT32)opj_int_min(dwt->dn, dwt->sn - b),
2650 OPJ_FLOAT32 * OPJ_RESTRICT aj;
2652 } opj_dwt97_decode_h_job_t;
2654 static void opj_dwt97_decode_h_func(void* user_data, opj_tls_t* tls)
2657 opj_dwt97_decode_h_job_t* job;
2658 OPJ_FLOAT32 * OPJ_RESTRICT aj;
2662 job = (opj_dwt97_decode_h_job_t*)user_data;
2665 assert((job->nb_rows % NB_ELTS_V8) == 0);
2668 for (j = 0; j + NB_ELTS_V8 <= job->nb_rows; j += NB_ELTS_V8) {
2670 opj_v8dwt_interleave_h(&job->h, aj, job->w, NB_ELTS_V8);
2671 opj_v8dwt_decode(&job->h);
2673 /* To be adapted if NB_ELTS_V8 changes */
2674 for (k = 0; k < job->rw; k++) {
2675 aj[k ] = job->h.wavelet[k].f[0];
2676 aj[k + (OPJ_SIZE_T)w ] = job->h.wavelet[k].f[1];
2677 aj[k + (OPJ_SIZE_T)w * 2] = job->h.wavelet[k].f[2];
2678 aj[k + (OPJ_SIZE_T)w * 3] = job->h.wavelet[k].f[3];
2680 for (k = 0; k < job->rw; k++) {
2681 aj[k + (OPJ_SIZE_T)w * 4] = job->h.wavelet[k].f[4];
2682 aj[k + (OPJ_SIZE_T)w * 5] = job->h.wavelet[k].f[5];
2683 aj[k + (OPJ_SIZE_T)w * 6] = job->h.wavelet[k].f[6];
2684 aj[k + (OPJ_SIZE_T)w * 7] = job->h.wavelet[k].f[7];
2687 aj += w * NB_ELTS_V8;
2690 opj_aligned_free(job->h.wavelet);
2699 OPJ_FLOAT32 * OPJ_RESTRICT aj;
2700 OPJ_UINT32 nb_columns;
2701 } opj_dwt97_decode_v_job_t;
2703 static void opj_dwt97_decode_v_func(void* user_data, opj_tls_t* tls)
2706 opj_dwt97_decode_v_job_t* job;
2707 OPJ_FLOAT32 * OPJ_RESTRICT aj;
2710 job = (opj_dwt97_decode_v_job_t*)user_data;
2712 assert((job->nb_columns % NB_ELTS_V8) == 0);
2715 for (j = 0; j + NB_ELTS_V8 <= job->nb_columns; j += NB_ELTS_V8) {
2718 opj_v8dwt_interleave_v(&job->v, aj, job->w, NB_ELTS_V8);
2719 opj_v8dwt_decode(&job->v);
2721 for (k = 0; k < job->rh; ++k) {
2722 memcpy(&aj[k * (OPJ_SIZE_T)job->w], &job->v.wavelet[k],
2723 NB_ELTS_V8 * sizeof(OPJ_FLOAT32));
2728 opj_aligned_free(job->v.wavelet);
2734 /* Inverse 9-7 wavelet transform in 2-D. */
2737 OPJ_BOOL opj_dwt_decode_tile_97(opj_thread_pool_t* tp,
2738 opj_tcd_tilecomp_t* OPJ_RESTRICT tilec,
2744 opj_tcd_resolution_t* res = tilec->resolutions;
2746 OPJ_UINT32 rw = (OPJ_UINT32)(res->x1 -
2747 res->x0); /* width of the resolution level computed */
2748 OPJ_UINT32 rh = (OPJ_UINT32)(res->y1 -
2749 res->y0); /* height of the resolution level computed */
2751 OPJ_UINT32 w = (OPJ_UINT32)(tilec->resolutions[tilec->minimum_num_resolutions -
2753 tilec->resolutions[tilec->minimum_num_resolutions - 1].x0);
2755 OPJ_SIZE_T l_data_size;
2756 const int num_threads = opj_thread_pool_get_thread_count(tp);
2762 l_data_size = opj_dwt_max_resolution(res, numres);
2763 /* overflow check */
2764 if (l_data_size > (SIZE_MAX / sizeof(opj_v8_t))) {
2765 /* FIXME event manager error callback */
2768 h.wavelet = (opj_v8_t*) opj_aligned_malloc(l_data_size * sizeof(opj_v8_t));
2770 /* FIXME event manager error callback */
2773 v.wavelet = h.wavelet;
2776 OPJ_FLOAT32 * OPJ_RESTRICT aj = (OPJ_FLOAT32*) tilec->data;
2779 h.sn = (OPJ_INT32)rw;
2780 v.sn = (OPJ_INT32)rh;
2784 rw = (OPJ_UINT32)(res->x1 -
2785 res->x0); /* width of the resolution level computed */
2786 rh = (OPJ_UINT32)(res->y1 -
2787 res->y0); /* height of the resolution level computed */
2789 h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
2790 h.cas = res->x0 % 2;
2793 h.win_l_x1 = (OPJ_UINT32)h.sn;
2795 h.win_h_x1 = (OPJ_UINT32)h.dn;
2797 if (num_threads <= 1 || rh < 2 * NB_ELTS_V8) {
2798 for (j = 0; j + (NB_ELTS_V8 - 1) < rh; j += NB_ELTS_V8) {
2800 opj_v8dwt_interleave_h(&h, aj, w, NB_ELTS_V8);
2801 opj_v8dwt_decode(&h);
2803 /* To be adapted if NB_ELTS_V8 changes */
2804 for (k = 0; k < rw; k++) {
2805 aj[k ] = h.wavelet[k].f[0];
2806 aj[k + (OPJ_SIZE_T)w ] = h.wavelet[k].f[1];
2807 aj[k + (OPJ_SIZE_T)w * 2] = h.wavelet[k].f[2];
2808 aj[k + (OPJ_SIZE_T)w * 3] = h.wavelet[k].f[3];
2810 for (k = 0; k < rw; k++) {
2811 aj[k + (OPJ_SIZE_T)w * 4] = h.wavelet[k].f[4];
2812 aj[k + (OPJ_SIZE_T)w * 5] = h.wavelet[k].f[5];
2813 aj[k + (OPJ_SIZE_T)w * 6] = h.wavelet[k].f[6];
2814 aj[k + (OPJ_SIZE_T)w * 7] = h.wavelet[k].f[7];
2817 aj += w * NB_ELTS_V8;
2820 OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
2823 if ((rh / NB_ELTS_V8) < num_jobs) {
2824 num_jobs = rh / NB_ELTS_V8;
2826 step_j = ((rh / num_jobs) / NB_ELTS_V8) * NB_ELTS_V8;
2827 for (j = 0; j < num_jobs; j++) {
2828 opj_dwt97_decode_h_job_t* job;
2830 job = (opj_dwt97_decode_h_job_t*) opj_malloc(sizeof(opj_dwt97_decode_h_job_t));
2832 opj_thread_pool_wait_completion(tp, 0);
2833 opj_aligned_free(h.wavelet);
2836 job->h.wavelet = (opj_v8_t*)opj_aligned_malloc(l_data_size * sizeof(opj_v8_t));
2837 if (!job->h.wavelet) {
2838 opj_thread_pool_wait_completion(tp, 0);
2840 opj_aligned_free(h.wavelet);
2846 job->h.win_l_x0 = h.win_l_x0;
2847 job->h.win_l_x1 = h.win_l_x1;
2848 job->h.win_h_x0 = h.win_h_x0;
2849 job->h.win_h_x1 = h.win_h_x1;
2853 job->nb_rows = (j + 1 == num_jobs) ? (rh & (OPJ_UINT32)~
2854 (NB_ELTS_V8 - 1)) - j * step_j : step_j;
2855 aj += w * job->nb_rows;
2856 opj_thread_pool_submit_job(tp, opj_dwt97_decode_h_func, job);
2858 opj_thread_pool_wait_completion(tp, 0);
2859 j = rh & (OPJ_UINT32)~(NB_ELTS_V8 - 1);
2864 opj_v8dwt_interleave_h(&h, aj, w, rh - j);
2865 opj_v8dwt_decode(&h);
2866 for (k = 0; k < rw; k++) {
2868 for (l = 0; l < rh - j; l++) {
2869 aj[k + (OPJ_SIZE_T)w * l ] = h.wavelet[k].f[l];
2874 v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
2875 v.cas = res->y0 % 2;
2877 v.win_l_x1 = (OPJ_UINT32)v.sn;
2879 v.win_h_x1 = (OPJ_UINT32)v.dn;
2881 aj = (OPJ_FLOAT32*) tilec->data;
2882 if (num_threads <= 1 || rw < 2 * NB_ELTS_V8) {
2883 for (j = rw; j > (NB_ELTS_V8 - 1); j -= NB_ELTS_V8) {
2886 opj_v8dwt_interleave_v(&v, aj, w, NB_ELTS_V8);
2887 opj_v8dwt_decode(&v);
2889 for (k = 0; k < rh; ++k) {
2890 memcpy(&aj[k * (OPJ_SIZE_T)w], &v.wavelet[k], NB_ELTS_V8 * sizeof(OPJ_FLOAT32));
2895 /* "bench_dwt -I" shows that scaling is poor, likely due to RAM
2896 transfer being the limiting factor. So limit the number of
2899 OPJ_UINT32 num_jobs = opj_uint_max((OPJ_UINT32)num_threads / 2, 2U);
2902 if ((rw / NB_ELTS_V8) < num_jobs) {
2903 num_jobs = rw / NB_ELTS_V8;
2905 step_j = ((rw / num_jobs) / NB_ELTS_V8) * NB_ELTS_V8;
2906 for (j = 0; j < num_jobs; j++) {
2907 opj_dwt97_decode_v_job_t* job;
2909 job = (opj_dwt97_decode_v_job_t*) opj_malloc(sizeof(opj_dwt97_decode_v_job_t));
2911 opj_thread_pool_wait_completion(tp, 0);
2912 opj_aligned_free(h.wavelet);
2915 job->v.wavelet = (opj_v8_t*)opj_aligned_malloc(l_data_size * sizeof(opj_v8_t));
2916 if (!job->v.wavelet) {
2917 opj_thread_pool_wait_completion(tp, 0);
2919 opj_aligned_free(h.wavelet);
2925 job->v.win_l_x0 = v.win_l_x0;
2926 job->v.win_l_x1 = v.win_l_x1;
2927 job->v.win_h_x0 = v.win_h_x0;
2928 job->v.win_h_x1 = v.win_h_x1;
2932 job->nb_columns = (j + 1 == num_jobs) ? (rw & (OPJ_UINT32)~
2933 (NB_ELTS_V8 - 1)) - j * step_j : step_j;
2934 aj += job->nb_columns;
2935 opj_thread_pool_submit_job(tp, opj_dwt97_decode_v_func, job);
2937 opj_thread_pool_wait_completion(tp, 0);
2940 if (rw & (NB_ELTS_V8 - 1)) {
2943 j = rw & (NB_ELTS_V8 - 1);
2945 opj_v8dwt_interleave_v(&v, aj, w, j);
2946 opj_v8dwt_decode(&v);
2948 for (k = 0; k < rh; ++k) {
2949 memcpy(&aj[k * (OPJ_SIZE_T)w], &v.wavelet[k],
2950 (OPJ_SIZE_T)j * sizeof(OPJ_FLOAT32));
2955 opj_aligned_free(h.wavelet);
2960 OPJ_BOOL opj_dwt_decode_partial_97(opj_tcd_tilecomp_t* OPJ_RESTRICT tilec,
2963 opj_sparse_array_int32_t* sa;
2967 /* This value matches the maximum left/right extension given in tables */
2968 /* F.2 and F.3 of the standard. Note: in opj_tcd_is_subband_area_of_interest() */
2969 /* we currently use 3. */
2970 const OPJ_UINT32 filter_width = 4U;
2972 opj_tcd_resolution_t* tr = tilec->resolutions;
2973 opj_tcd_resolution_t* tr_max = &(tilec->resolutions[numres - 1]);
2975 OPJ_UINT32 rw = (OPJ_UINT32)(tr->x1 -
2976 tr->x0); /* width of the resolution level computed */
2977 OPJ_UINT32 rh = (OPJ_UINT32)(tr->y1 -
2978 tr->y0); /* height of the resolution level computed */
2980 OPJ_SIZE_T l_data_size;
2982 /* Compute the intersection of the area of interest, expressed in tile coordinates */
2983 /* with the tile coordinates */
2984 OPJ_UINT32 win_tcx0 = tilec->win_x0;
2985 OPJ_UINT32 win_tcy0 = tilec->win_y0;
2986 OPJ_UINT32 win_tcx1 = tilec->win_x1;
2987 OPJ_UINT32 win_tcy1 = tilec->win_y1;
2989 if (tr_max->x0 == tr_max->x1 || tr_max->y0 == tr_max->y1) {
2993 sa = opj_dwt_init_sparse_array(tilec, numres);
2999 OPJ_BOOL ret = opj_sparse_array_int32_read(sa,
3000 tr_max->win_x0 - (OPJ_UINT32)tr_max->x0,
3001 tr_max->win_y0 - (OPJ_UINT32)tr_max->y0,
3002 tr_max->win_x1 - (OPJ_UINT32)tr_max->x0,
3003 tr_max->win_y1 - (OPJ_UINT32)tr_max->y0,
3005 1, tr_max->win_x1 - tr_max->win_x0,
3009 opj_sparse_array_int32_free(sa);
3013 l_data_size = opj_dwt_max_resolution(tr, numres);
3014 /* overflow check */
3015 if (l_data_size > (SIZE_MAX / sizeof(opj_v8_t))) {
3016 /* FIXME event manager error callback */
3017 opj_sparse_array_int32_free(sa);
3020 h.wavelet = (opj_v8_t*) opj_aligned_malloc(l_data_size * sizeof(opj_v8_t));
3022 /* FIXME event manager error callback */
3023 opj_sparse_array_int32_free(sa);
3026 v.wavelet = h.wavelet;
3028 for (resno = 1; resno < numres; resno ++) {
3030 /* Window of interest subband-based coordinates */
3031 OPJ_UINT32 win_ll_x0, win_ll_y0, win_ll_x1, win_ll_y1;
3032 OPJ_UINT32 win_hl_x0, win_hl_x1;
3033 OPJ_UINT32 win_lh_y0, win_lh_y1;
3034 /* Window of interest tile-resolution-based coordinates */
3035 OPJ_UINT32 win_tr_x0, win_tr_x1, win_tr_y0, win_tr_y1;
3036 /* Tile-resolution subband-based coordinates */
3037 OPJ_UINT32 tr_ll_x0, tr_ll_y0, tr_hl_x0, tr_lh_y0;
3041 h.sn = (OPJ_INT32)rw;
3042 v.sn = (OPJ_INT32)rh;
3044 rw = (OPJ_UINT32)(tr->x1 - tr->x0);
3045 rh = (OPJ_UINT32)(tr->y1 - tr->y0);
3047 h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
3050 v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
3053 /* Get the subband coordinates for the window of interest */
3055 opj_dwt_get_band_coordinates(tilec, resno, 0,
3056 win_tcx0, win_tcy0, win_tcx1, win_tcy1,
3057 &win_ll_x0, &win_ll_y0,
3058 &win_ll_x1, &win_ll_y1);
3061 opj_dwt_get_band_coordinates(tilec, resno, 1,
3062 win_tcx0, win_tcy0, win_tcx1, win_tcy1,
3063 &win_hl_x0, NULL, &win_hl_x1, NULL);
3066 opj_dwt_get_band_coordinates(tilec, resno, 2,
3067 win_tcx0, win_tcy0, win_tcx1, win_tcy1,
3068 NULL, &win_lh_y0, NULL, &win_lh_y1);
3070 /* Beware: band index for non-LL0 resolution are 0=HL, 1=LH and 2=HH */
3071 tr_ll_x0 = (OPJ_UINT32)tr->bands[1].x0;
3072 tr_ll_y0 = (OPJ_UINT32)tr->bands[0].y0;
3073 tr_hl_x0 = (OPJ_UINT32)tr->bands[0].x0;
3074 tr_lh_y0 = (OPJ_UINT32)tr->bands[1].y0;
3076 /* Subtract the origin of the bands for this tile, to the subwindow */
3077 /* of interest band coordinates, so as to get them relative to the */
3079 win_ll_x0 = opj_uint_subs(win_ll_x0, tr_ll_x0);
3080 win_ll_y0 = opj_uint_subs(win_ll_y0, tr_ll_y0);
3081 win_ll_x1 = opj_uint_subs(win_ll_x1, tr_ll_x0);
3082 win_ll_y1 = opj_uint_subs(win_ll_y1, tr_ll_y0);
3083 win_hl_x0 = opj_uint_subs(win_hl_x0, tr_hl_x0);
3084 win_hl_x1 = opj_uint_subs(win_hl_x1, tr_hl_x0);
3085 win_lh_y0 = opj_uint_subs(win_lh_y0, tr_lh_y0);
3086 win_lh_y1 = opj_uint_subs(win_lh_y1, tr_lh_y0);
3088 opj_dwt_segment_grow(filter_width, (OPJ_UINT32)h.sn, &win_ll_x0, &win_ll_x1);
3089 opj_dwt_segment_grow(filter_width, (OPJ_UINT32)h.dn, &win_hl_x0, &win_hl_x1);
3091 opj_dwt_segment_grow(filter_width, (OPJ_UINT32)v.sn, &win_ll_y0, &win_ll_y1);
3092 opj_dwt_segment_grow(filter_width, (OPJ_UINT32)v.dn, &win_lh_y0, &win_lh_y1);
3094 /* Compute the tile-resolution-based coordinates for the window of interest */
3096 win_tr_x0 = opj_uint_min(2 * win_ll_x0, 2 * win_hl_x0 + 1);
3097 win_tr_x1 = opj_uint_min(opj_uint_max(2 * win_ll_x1, 2 * win_hl_x1 + 1), rw);
3099 win_tr_x0 = opj_uint_min(2 * win_hl_x0, 2 * win_ll_x0 + 1);
3100 win_tr_x1 = opj_uint_min(opj_uint_max(2 * win_hl_x1, 2 * win_ll_x1 + 1), rw);
3104 win_tr_y0 = opj_uint_min(2 * win_ll_y0, 2 * win_lh_y0 + 1);
3105 win_tr_y1 = opj_uint_min(opj_uint_max(2 * win_ll_y1, 2 * win_lh_y1 + 1), rh);
3107 win_tr_y0 = opj_uint_min(2 * win_lh_y0, 2 * win_ll_y0 + 1);
3108 win_tr_y1 = opj_uint_min(opj_uint_max(2 * win_lh_y1, 2 * win_ll_y1 + 1), rh);
3111 h.win_l_x0 = win_ll_x0;
3112 h.win_l_x1 = win_ll_x1;
3113 h.win_h_x0 = win_hl_x0;
3114 h.win_h_x1 = win_hl_x1;
3115 for (j = 0; j + (NB_ELTS_V8 - 1) < rh; j += NB_ELTS_V8) {
3116 if ((j + (NB_ELTS_V8 - 1) >= win_ll_y0 && j < win_ll_y1) ||
3117 (j + (NB_ELTS_V8 - 1) >= win_lh_y0 + (OPJ_UINT32)v.sn &&
3118 j < win_lh_y1 + (OPJ_UINT32)v.sn)) {
3119 opj_v8dwt_interleave_partial_h(&h, sa, j, opj_uint_min(NB_ELTS_V8, rh - j));
3120 opj_v8dwt_decode(&h);
3121 if (!opj_sparse_array_int32_write(sa,
3123 win_tr_x1, j + NB_ELTS_V8,
3124 (OPJ_INT32*)&h.wavelet[win_tr_x0].f[0],
3125 NB_ELTS_V8, 1, OPJ_TRUE)) {
3126 /* FIXME event manager error callback */
3127 opj_sparse_array_int32_free(sa);
3128 opj_aligned_free(h.wavelet);
3135 ((j + (NB_ELTS_V8 - 1) >= win_ll_y0 && j < win_ll_y1) ||
3136 (j + (NB_ELTS_V8 - 1) >= win_lh_y0 + (OPJ_UINT32)v.sn &&
3137 j < win_lh_y1 + (OPJ_UINT32)v.sn))) {
3138 opj_v8dwt_interleave_partial_h(&h, sa, j, rh - j);
3139 opj_v8dwt_decode(&h);
3140 if (!opj_sparse_array_int32_write(sa,
3143 (OPJ_INT32*)&h.wavelet[win_tr_x0].f[0],
3144 NB_ELTS_V8, 1, OPJ_TRUE)) {
3145 /* FIXME event manager error callback */
3146 opj_sparse_array_int32_free(sa);
3147 opj_aligned_free(h.wavelet);
3152 v.win_l_x0 = win_ll_y0;
3153 v.win_l_x1 = win_ll_y1;
3154 v.win_h_x0 = win_lh_y0;
3155 v.win_h_x1 = win_lh_y1;
3156 for (j = win_tr_x0; j < win_tr_x1; j += NB_ELTS_V8) {
3157 OPJ_UINT32 nb_elts = opj_uint_min(NB_ELTS_V8, win_tr_x1 - j);
3159 opj_v8dwt_interleave_partial_v(&v, sa, j, nb_elts);
3160 opj_v8dwt_decode(&v);
3162 if (!opj_sparse_array_int32_write(sa,
3164 j + nb_elts, win_tr_y1,
3165 (OPJ_INT32*)&h.wavelet[win_tr_y0].f[0],
3166 1, NB_ELTS_V8, OPJ_TRUE)) {
3167 /* FIXME event manager error callback */
3168 opj_sparse_array_int32_free(sa);
3169 opj_aligned_free(h.wavelet);
3176 OPJ_BOOL ret = opj_sparse_array_int32_read(sa,
3177 tr_max->win_x0 - (OPJ_UINT32)tr_max->x0,
3178 tr_max->win_y0 - (OPJ_UINT32)tr_max->y0,
3179 tr_max->win_x1 - (OPJ_UINT32)tr_max->x0,
3180 tr_max->win_y1 - (OPJ_UINT32)tr_max->y0,
3182 1, tr_max->win_x1 - tr_max->win_x0,
3187 opj_sparse_array_int32_free(sa);
3189 opj_aligned_free(h.wavelet);
3194 OPJ_BOOL opj_dwt_decode_real(opj_tcd_t *p_tcd,
3195 opj_tcd_tilecomp_t* OPJ_RESTRICT tilec,
3198 if (p_tcd->whole_tile_decoding) {
3199 return opj_dwt_decode_tile_97(p_tcd->thread_pool, tilec, numres);
3201 return opj_dwt_decode_partial_97(tilec, numres);