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.
42 #include <xmmintrin.h>
45 #include <emmintrin.h>
50 #include "opj_includes.h"
52 /** @defgroup DWT DWT - Implementation of a discrete wavelet transform */
55 #define OPJ_WS(i) v->mem[(i)*2]
56 #define OPJ_WD(i) v->mem[(1+(i)*2)]
58 #define PARALLEL_COLS_53 8
60 /** @name Local data structures */
63 typedef struct dwt_local {
74 typedef struct v4dwt_local {
81 static const OPJ_FLOAT32 opj_dwt_alpha = 1.586134342f; /* 12994 */
82 static const OPJ_FLOAT32 opj_dwt_beta = 0.052980118f; /* 434 */
83 static const OPJ_FLOAT32 opj_dwt_gamma = -0.882911075f; /* -7233 */
84 static const OPJ_FLOAT32 opj_dwt_delta = -0.443506852f; /* -3633 */
86 static const OPJ_FLOAT32 opj_K = 1.230174105f; /* 10078 */
87 static const OPJ_FLOAT32 opj_c13318 = 1.625732422f;
92 Virtual function type for wavelet transform in 1-D
94 typedef void (*DWT1DFN)(const opj_dwt_t* v);
96 /** @name Local static functions */
100 Forward lazy transform (horizontal)
102 static void opj_dwt_deinterleave_h(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
103 OPJ_INT32 sn, OPJ_INT32 cas);
105 Forward lazy transform (vertical)
107 static void opj_dwt_deinterleave_v(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
108 OPJ_INT32 sn, OPJ_INT32 x, OPJ_INT32 cas);
110 Forward 5-3 wavelet transform in 1-D
112 static void opj_dwt_encode_1(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
115 Forward 9-7 wavelet transform in 1-D
117 static void opj_dwt_encode_1_real(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
120 Explicit calculation of the Quantization Stepsizes
122 static void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps,
123 opj_stepsize_t *bandno_stepsize);
125 Inverse wavelet transform in 2-D.
127 static OPJ_BOOL opj_dwt_decode_tile(opj_thread_pool_t* tp,
128 opj_tcd_tilecomp_t* tilec, OPJ_UINT32 i);
130 static OPJ_BOOL opj_dwt_encode_procedure(opj_tcd_tilecomp_t * tilec,
131 void (*p_function)(OPJ_INT32 *, OPJ_INT32, OPJ_INT32, OPJ_INT32));
133 static OPJ_UINT32 opj_dwt_max_resolution(opj_tcd_resolution_t* OPJ_RESTRICT r,
137 /* Inverse 9-7 wavelet transform in 1-D. */
139 static void opj_v4dwt_decode(opj_v4dwt_t* OPJ_RESTRICT dwt);
141 static void opj_v4dwt_interleave_h(opj_v4dwt_t* OPJ_RESTRICT w,
142 OPJ_FLOAT32* OPJ_RESTRICT a, OPJ_INT32 x, OPJ_INT32 size);
144 static void opj_v4dwt_interleave_v(opj_v4dwt_t* OPJ_RESTRICT v,
145 OPJ_FLOAT32* OPJ_RESTRICT a, OPJ_INT32 x, OPJ_INT32 nb_elts_read);
148 static void opj_v4dwt_decode_step1_sse(opj_v4_t* w, OPJ_INT32 count,
151 static void opj_v4dwt_decode_step2_sse(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k,
152 OPJ_INT32 m, __m128 c);
155 static void opj_v4dwt_decode_step1(opj_v4_t* w, OPJ_INT32 count,
156 const OPJ_FLOAT32 c);
158 static void opj_v4dwt_decode_step2(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k,
159 OPJ_INT32 m, OPJ_FLOAT32 c);
167 #define OPJ_S(i) a[(i)*2]
168 #define OPJ_D(i) a[(1+(i)*2)]
169 #define OPJ_S_(i) ((i)<0?OPJ_S(0):((i)>=sn?OPJ_S(sn-1):OPJ_S(i)))
170 #define OPJ_D_(i) ((i)<0?OPJ_D(0):((i)>=dn?OPJ_D(dn-1):OPJ_D(i)))
172 #define OPJ_SS_(i) ((i)<0?OPJ_S(0):((i)>=dn?OPJ_S(dn-1):OPJ_S(i)))
173 #define OPJ_DD_(i) ((i)<0?OPJ_D(0):((i)>=sn?OPJ_D(sn-1):OPJ_D(i)))
176 /* This table contains the norms of the 5-3 wavelets for different bands. */
178 static const OPJ_FLOAT64 opj_dwt_norms[4][10] = {
179 {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
180 {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
181 {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
182 {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
186 /* This table contains the norms of the 9-7 wavelets for different bands. */
188 static const OPJ_FLOAT64 opj_dwt_norms_real[4][10] = {
189 {1.000, 1.965, 4.177, 8.403, 16.90, 33.84, 67.69, 135.3, 270.6, 540.9},
190 {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
191 {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
192 {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2}
196 ==========================================================
198 ==========================================================
202 /* Forward lazy transform (horizontal). */
204 static void opj_dwt_deinterleave_h(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
205 OPJ_INT32 sn, OPJ_INT32 cas)
208 OPJ_INT32 * l_dest = b;
209 OPJ_INT32 * l_src = a + cas;
211 for (i = 0; i < sn; ++i) {
219 for (i = 0; i < dn; ++i) {
226 /* Forward lazy transform (vertical). */
228 static void opj_dwt_deinterleave_v(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
229 OPJ_INT32 sn, OPJ_INT32 x, OPJ_INT32 cas)
232 OPJ_INT32 * l_dest = b;
233 OPJ_INT32 * l_src = a + cas;
239 } /* b[i*x]=a[2*i+cas]; */
249 } /*b[(sn+i)*x]=a[(2*i+1-cas)];*/
252 #ifdef STANDARD_SLOW_VERSION
254 /* Inverse lazy transform (horizontal). */
256 static void opj_dwt_interleave_h(const opj_dwt_t* h, OPJ_INT32 *a)
259 OPJ_INT32 *bi = h->mem + h->cas;
266 bi = h->mem + 1 - h->cas;
275 /* Inverse lazy transform (vertical). */
277 static void opj_dwt_interleave_v(const opj_dwt_t* v, OPJ_INT32 *a, OPJ_INT32 x)
280 OPJ_INT32 *bi = v->mem + v->cas;
287 ai = a + (v->sn * x);
288 bi = v->mem + 1 - v->cas;
297 #endif /* STANDARD_SLOW_VERSION */
300 /* Forward 5-3 wavelet transform in 1-D. */
302 static void opj_dwt_encode_1(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
308 if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
309 for (i = 0; i < dn; i++) {
310 OPJ_D(i) -= (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
312 for (i = 0; i < sn; i++) {
313 OPJ_S(i) += (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
317 if (!sn && dn == 1) { /* NEW : CASE ONE ELEMENT */
320 for (i = 0; i < dn; i++) {
321 OPJ_S(i) -= (OPJ_DD_(i) + OPJ_DD_(i - 1)) >> 1;
323 for (i = 0; i < sn; i++) {
324 OPJ_D(i) += (OPJ_SS_(i) + OPJ_SS_(i + 1) + 2) >> 2;
330 #ifdef STANDARD_SLOW_VERSION
332 /* Inverse 5-3 wavelet transform in 1-D. */
334 static void opj_dwt_decode_1_(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
340 if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
341 for (i = 0; i < sn; i++) {
342 OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
344 for (i = 0; i < dn; i++) {
345 OPJ_D(i) += (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
349 if (!sn && dn == 1) { /* NEW : CASE ONE ELEMENT */
352 for (i = 0; i < sn; i++) {
353 OPJ_D(i) -= (OPJ_SS_(i) + OPJ_SS_(i + 1) + 2) >> 2;
355 for (i = 0; i < dn; i++) {
356 OPJ_S(i) += (OPJ_DD_(i) + OPJ_DD_(i - 1)) >> 1;
362 static void opj_dwt_decode_1(const opj_dwt_t *v)
364 opj_dwt_decode_1_(v->mem, v->dn, v->sn, v->cas);
367 #endif /* STANDARD_SLOW_VERSION */
369 #if !defined(STANDARD_SLOW_VERSION)
370 static void opj_idwt53_h_cas0(OPJ_INT32* tmp,
376 const OPJ_INT32* in_even = &tiledp[0];
377 const OPJ_INT32* in_odd = &tiledp[sn];
379 #ifdef TWO_PASS_VERSION
380 /* For documentation purpose: performs lifting in two iterations, */
381 /* but withtmp explicit interleaving */
386 tmp[0] = in_even[0] - ((in_odd[0] + 1) >> 1);
387 for (i = 2, j = 0; i <= len - 2; i += 2, j++) {
388 tmp[i] = in_even[j + 1] - ((in_odd[j] + in_odd[j + 1] + 2) >> 2);
390 if (len & 1) { /* if len is odd */
391 tmp[len - 1] = in_even[(len - 1) / 2] - ((in_odd[(len - 2) / 2] + 1) >> 1);
395 for (i = 1, j = 0; i < len - 1; i += 2, j++) {
396 tmp[i] = in_odd[j] + ((tmp[i - 1] + tmp[i + 1]) >> 1);
398 if (!(len & 1)) { /* if len is even */
399 tmp[len - 1] = in_odd[(len - 1) / 2] + tmp[len - 2];
402 OPJ_INT32 d1c, d1n, s1n, s0c, s0n;
406 /* Improved version of the TWO_PASS_VERSION: */
407 /* Performs lifting in one single iteration. Saves memory */
408 /* accesses and explicit interleaving. */
411 s0n = s1n - ((d1n + 1) >> 1);
413 for (i = 0, j = 1; i < (len - 3); i += 2, j++) {
420 s0n = s1n - ((d1c + d1n + 2) >> 2);
423 tmp[i + 1] = d1c + ((s0c + s0n) >> 1);
429 tmp[len - 1] = in_even[(len - 1) / 2] - ((d1n + 1) >> 1);
430 tmp[len - 2] = d1n + ((s0n + tmp[len - 1]) >> 1);
432 tmp[len - 1] = d1n + s0n;
435 memcpy(tiledp, tmp, (OPJ_UINT32)len * sizeof(OPJ_INT32));
438 static void opj_idwt53_h_cas1(OPJ_INT32* tmp,
444 const OPJ_INT32* in_even = &tiledp[sn];
445 const OPJ_INT32* in_odd = &tiledp[0];
447 #ifdef TWO_PASS_VERSION
448 /* For documentation purpose: performs lifting in two iterations, */
449 /* but withtmp explicit interleaving */
454 for (i = 1, j = 0; i < len - 1; i += 2, j++) {
455 tmp[i] = in_odd[j] - ((in_even[j] + in_even[j + 1] + 2) >> 2);
458 tmp[len - 1] = in_odd[len / 2 - 1] - ((in_even[len / 2 - 1] + 1) >> 1);
462 tmp[0] = in_even[0] + tmp[1];
463 for (i = 2, j = 1; i < len - 1; i += 2, j++) {
464 tmp[i] = in_even[j] + ((tmp[i + 1] + tmp[i - 1]) >> 1);
467 tmp[len - 1] = in_even[len / 2] + tmp[len - 2];
470 OPJ_INT32 s1, s2, dc, dn;
474 /* Improved version of the TWO_PASS_VERSION: */
475 /* Performs lifting in one single iteration. Saves memory */
476 /* accesses and explicit interleaving. */
479 dc = in_odd[0] - ((in_even[0] + s1 + 2) >> 2);
480 tmp[0] = in_even[0] + dc;
482 for (i = 1, j = 1; i < (len - 2 - !(len & 1)); i += 2, j++) {
486 dn = in_odd[j] - ((s1 + s2 + 2) >> 2);
488 tmp[i + 1] = s1 + ((dn + dc) >> 1);
497 dn = in_odd[len / 2 - 1] - ((s1 + 1) >> 1);
498 tmp[len - 2] = s1 + ((dn + dc) >> 1);
501 tmp[len - 1] = s1 + dc;
504 memcpy(tiledp, tmp, (OPJ_UINT32)len * sizeof(OPJ_INT32));
508 #endif /* !defined(STANDARD_SLOW_VERSION) */
511 /* Inverse 5-3 wavelet transform in 1-D for one row. */
513 /* Performs interleave, inverse wavelet transform and copy back to buffer */
514 static void opj_idwt53_h(const opj_dwt_t *dwt,
517 #ifdef STANDARD_SLOW_VERSION
518 /* For documentation purpose */
519 opj_dwt_interleave_h(dwt, tiledp);
520 opj_dwt_decode_1(dwt);
521 memcpy(tiledp, dwt->mem, (OPJ_UINT32)(dwt->sn + dwt->dn) * sizeof(OPJ_INT32));
523 const OPJ_INT32 sn = dwt->sn;
524 const OPJ_INT32 len = sn + dwt->dn;
525 if (dwt->cas == 0) { /* Left-most sample is on even coordinate */
527 opj_idwt53_h_cas0(dwt->mem, sn, len, tiledp);
529 /* Unmodified value */
531 } else { /* Left-most sample is on odd coordinate */
534 } else if (len == 2) {
535 OPJ_INT32* out = dwt->mem;
536 const OPJ_INT32* in_even = &tiledp[sn];
537 const OPJ_INT32* in_odd = &tiledp[0];
538 out[1] = in_odd[0] - ((in_even[0] + 1) >> 1);
539 out[0] = in_even[0] + out[1];
540 memcpy(tiledp, dwt->mem, (OPJ_UINT32)len * sizeof(OPJ_INT32));
541 } else if (len > 2) {
542 opj_idwt53_h_cas1(dwt->mem, sn, len, tiledp);
548 #if defined(__SSE2__) && !defined(STANDARD_SLOW_VERSION)
550 /* Conveniency macros to improve the readabilty of the formulas */
551 #define LOADU(x) _mm_loadu_si128((const __m128i*)(x))
552 #define STORE(x,y) _mm_store_si128((__m128i*)(x),(y))
553 #define ADD(x,y) _mm_add_epi32((x),(y))
554 #define ADD3(x,y,z) ADD(ADD(x,y),z)
555 #define SUB(x,y) _mm_sub_epi32((x),(y))
556 #define SAR(x,y) _mm_srai_epi32((x),(y))
558 /** Vertical inverse 5x3 wavelet transform for 8 columns, when top-most
559 * pixel is on even coordinate */
560 static void opj_idwt53_v_cas0_8cols_SSE2(
564 OPJ_INT32* tiledp_col,
565 const OPJ_INT32 stride)
567 const OPJ_INT32* in_even = &tiledp_col[0];
568 const OPJ_INT32* in_odd = &tiledp_col[sn * stride];
571 __m128i d1c_0, d1n_0, s1n_0, s0c_0, s0n_0;
572 __m128i d1c_1, d1n_1, s1n_1, s0c_1, s0n_1;
573 const __m128i two = _mm_set1_epi32(2);
576 assert(PARALLEL_COLS_53 == 8);
578 s1n_0 = LOADU(in_even + 0);
579 s1n_1 = LOADU(in_even + 4);
580 d1n_0 = LOADU(in_odd);
581 d1n_1 = LOADU(in_odd + 4);
583 /* s0n = s1n - ((d1n + 1) >> 1); <==> */
584 /* s0n = s1n - ((d1n + d1n + 2) >> 2); */
585 s0n_0 = SUB(s1n_0, SAR(ADD3(d1n_0, d1n_0, two), 2));
586 s0n_1 = SUB(s1n_1, SAR(ADD3(d1n_1, d1n_1, two), 2));
588 for (i = 0, j = 1; i < (len - 3); i += 2, j++) {
594 s1n_0 = LOADU(in_even + j * stride);
595 s1n_1 = LOADU(in_even + j * stride + 4);
596 d1n_0 = LOADU(in_odd + j * stride);
597 d1n_1 = LOADU(in_odd + j * stride + 4);
599 /*s0n = s1n - ((d1c + d1n + 2) >> 2);*/
600 s0n_0 = SUB(s1n_0, SAR(ADD3(d1c_0, d1n_0, two), 2));
601 s0n_1 = SUB(s1n_1, SAR(ADD3(d1c_1, d1n_1, two), 2));
603 STORE(tmp + PARALLEL_COLS_53 * (i + 0), s0c_0);
604 STORE(tmp + PARALLEL_COLS_53 * (i + 0) + 4, s0c_1);
606 /* d1c + ((s0c + s0n) >> 1) */
607 STORE(tmp + PARALLEL_COLS_53 * (i + 1) + 0,
608 ADD(d1c_0, SAR(ADD(s0c_0, s0n_0), 1)));
609 STORE(tmp + PARALLEL_COLS_53 * (i + 1) + 4,
610 ADD(d1c_1, SAR(ADD(s0c_1, s0n_1), 1)));
613 STORE(tmp + PARALLEL_COLS_53 * (i + 0) + 0, s0n_0);
614 STORE(tmp + PARALLEL_COLS_53 * (i + 0) + 4, s0n_1);
617 __m128i tmp_len_minus_1;
618 s1n_0 = LOADU(in_even + ((len - 1) / 2) * stride);
619 /* tmp_len_minus_1 = s1n - ((d1n + 1) >> 1); */
620 tmp_len_minus_1 = SUB(s1n_0, SAR(ADD3(d1n_0, d1n_0, two), 2));
621 STORE(tmp + 8 * (len - 1), tmp_len_minus_1);
622 /* d1n + ((s0n + tmp_len_minus_1) >> 1) */
623 STORE(tmp + 8 * (len - 2),
624 ADD(d1n_0, SAR(ADD(s0n_0, tmp_len_minus_1), 1)));
626 s1n_1 = LOADU(in_even + ((len - 1) / 2) * stride + 4);
627 /* tmp_len_minus_1 = s1n - ((d1n + 1) >> 1); */
628 tmp_len_minus_1 = SUB(s1n_1, SAR(ADD3(d1n_1, d1n_1, two), 2));
629 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 4, tmp_len_minus_1);
630 /* d1n + ((s0n + tmp_len_minus_1) >> 1) */
631 STORE(tmp + PARALLEL_COLS_53 * (len - 2) + 4,
632 ADD(d1n_1, SAR(ADD(s0n_1, tmp_len_minus_1), 1)));
636 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 0, ADD(d1n_0, s0n_0));
637 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 4, ADD(d1n_1, s0n_1));
640 for (i = 0; i < len; ++i) {
641 memcpy(&tiledp_col[i * stride],
642 &tmp[PARALLEL_COLS_53 * i],
643 PARALLEL_COLS_53 * sizeof(OPJ_INT32));
648 /** Vertical inverse 5x3 wavelet transform for 8 columns, when top-most
649 * pixel is on odd coordinate */
650 static void opj_idwt53_v_cas1_8cols_SSE2(
654 OPJ_INT32* tiledp_col,
655 const OPJ_INT32 stride)
659 __m128i s1_0, s2_0, dc_0, dn_0;
660 __m128i s1_1, s2_1, dc_1, dn_1;
661 const __m128i two = _mm_set1_epi32(2);
663 const OPJ_INT32* in_even = &tiledp_col[sn * stride];
664 const OPJ_INT32* in_odd = &tiledp_col[0];
667 assert(PARALLEL_COLS_53 == 8);
669 s1_0 = LOADU(in_even + stride);
670 /* in_odd[0] - ((in_even[0] + s1 + 2) >> 2); */
671 dc_0 = _mm_sub_epi32(
673 SAR(ADD3(LOADU(in_even + 0), s1_0, two), 2));
674 STORE(tmp + PARALLEL_COLS_53 * 0, ADD(LOADU(in_even + 0), dc_0));
676 s1_1 = LOADU(in_even + stride + 4);
677 /* in_odd[0] - ((in_even[0] + s1 + 2) >> 2); */
678 dc_1 = _mm_sub_epi32(
680 SAR(ADD3(LOADU(in_even + 4), s1_1, two), 2));
681 STORE(tmp + PARALLEL_COLS_53 * 0 + 4, ADD(LOADU(in_even + 4), dc_1));
683 for (i = 1, j = 1; i < (len - 2 - !(len & 1)); i += 2, j++) {
685 s2_0 = LOADU(in_even + (j + 1) * stride);
686 s2_1 = LOADU(in_even + (j + 1) * stride + 4);
688 /* dn = in_odd[j * stride] - ((s1 + s2 + 2) >> 2); */
689 dn_0 = SUB(LOADU(in_odd + j * stride),
690 SAR(ADD3(s1_0, s2_0, two), 2));
691 dn_1 = SUB(LOADU(in_odd + j * stride + 4),
692 SAR(ADD3(s1_1, s2_1, two), 2));
694 STORE(tmp + PARALLEL_COLS_53 * i, dc_0);
695 STORE(tmp + PARALLEL_COLS_53 * i + 4, dc_1);
697 /* tmp[i + 1] = s1 + ((dn + dc) >> 1); */
698 STORE(tmp + PARALLEL_COLS_53 * (i + 1) + 0,
699 ADD(s1_0, SAR(ADD(dn_0, dc_0), 1)));
700 STORE(tmp + PARALLEL_COLS_53 * (i + 1) + 4,
701 ADD(s1_1, SAR(ADD(dn_1, dc_1), 1)));
708 STORE(tmp + PARALLEL_COLS_53 * i, dc_0);
709 STORE(tmp + PARALLEL_COLS_53 * i + 4, dc_1);
712 /*dn = in_odd[(len / 2 - 1) * stride] - ((s1 + 1) >> 1); */
713 dn_0 = SUB(LOADU(in_odd + (len / 2 - 1) * stride),
714 SAR(ADD3(s1_0, s1_0, two), 2));
715 dn_1 = SUB(LOADU(in_odd + (len / 2 - 1) * stride + 4),
716 SAR(ADD3(s1_1, s1_1, two), 2));
718 /* tmp[len - 2] = s1 + ((dn + dc) >> 1); */
719 STORE(tmp + PARALLEL_COLS_53 * (len - 2) + 0,
720 ADD(s1_0, SAR(ADD(dn_0, dc_0), 1)));
721 STORE(tmp + PARALLEL_COLS_53 * (len - 2) + 4,
722 ADD(s1_1, SAR(ADD(dn_1, dc_1), 1)));
724 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 0, dn_0);
725 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 4, dn_1);
727 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 0, ADD(s1_0, dc_0));
728 STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 4, ADD(s1_1, dc_1));
731 for (i = 0; i < len; ++i) {
732 memcpy(&tiledp_col[i * stride],
733 &tmp[PARALLEL_COLS_53 * i],
734 PARALLEL_COLS_53 * sizeof(OPJ_INT32));
745 #endif /* defined(__SSE2__) && !defined(STANDARD_SLOW_VERSION) */
747 #if !defined(STANDARD_SLOW_VERSION)
748 /** Vertical inverse 5x3 wavelet transform for one column, when top-most
749 * pixel is on even coordinate */
750 static void opj_idwt3_v_cas0(OPJ_INT32* tmp,
753 OPJ_INT32* tiledp_col,
754 const OPJ_INT32 stride)
757 OPJ_INT32 d1c, d1n, s1n, s0c, s0n;
761 /* Performs lifting in one single iteration. Saves memory */
762 /* accesses and explicit interleaving. */
765 d1n = tiledp_col[sn * stride];
766 s0n = s1n - ((d1n + 1) >> 1);
768 for (i = 0, j = 0; i < (len - 3); i += 2, j++) {
772 s1n = tiledp_col[(j + 1) * stride];
773 d1n = tiledp_col[(sn + j + 1) * stride];
775 s0n = s1n - ((d1c + d1n + 2) >> 2);
778 tmp[i + 1] = d1c + ((s0c + s0n) >> 1);
785 tiledp_col[((len - 1) / 2) * stride] -
787 tmp[len - 2] = d1n + ((s0n + tmp[len - 1]) >> 1);
789 tmp[len - 1] = d1n + s0n;
792 for (i = 0; i < len; ++i) {
793 tiledp_col[i * stride] = tmp[i];
798 /** Vertical inverse 5x3 wavelet transform for one column, when top-most
799 * pixel is on odd coordinate */
800 static void opj_idwt3_v_cas1(OPJ_INT32* tmp,
803 OPJ_INT32* tiledp_col,
804 const OPJ_INT32 stride)
807 OPJ_INT32 s1, s2, dc, dn;
808 const OPJ_INT32* in_even = &tiledp_col[sn * stride];
809 const OPJ_INT32* in_odd = &tiledp_col[0];
813 /* Performs lifting in one single iteration. Saves memory */
814 /* accesses and explicit interleaving. */
816 s1 = in_even[stride];
817 dc = in_odd[0] - ((in_even[0] + s1 + 2) >> 2);
818 tmp[0] = in_even[0] + dc;
819 for (i = 1, j = 1; i < (len - 2 - !(len & 1)); i += 2, j++) {
821 s2 = in_even[(j + 1) * stride];
823 dn = in_odd[j * stride] - ((s1 + s2 + 2) >> 2);
825 tmp[i + 1] = s1 + ((dn + dc) >> 1);
832 dn = in_odd[(len / 2 - 1) * stride] - ((s1 + 1) >> 1);
833 tmp[len - 2] = s1 + ((dn + dc) >> 1);
836 tmp[len - 1] = s1 + dc;
839 for (i = 0; i < len; ++i) {
840 tiledp_col[i * stride] = tmp[i];
843 #endif /* !defined(STANDARD_SLOW_VERSION) */
846 /* Inverse vertical 5-3 wavelet transform in 1-D for several columns. */
848 /* Performs interleave, inverse wavelet transform and copy back to buffer */
849 static void opj_idwt53_v(const opj_dwt_t *dwt,
850 OPJ_INT32* tiledp_col,
854 #ifdef STANDARD_SLOW_VERSION
855 /* For documentation purpose */
857 for (c = 0; c < nb_cols; c ++) {
858 opj_dwt_interleave_v(dwt, tiledp_col + c, stride);
859 opj_dwt_decode_1(dwt);
860 for (k = 0; k < dwt->sn + dwt->dn; ++k) {
861 tiledp_col[c + k * stride] = dwt->mem[k];
865 const OPJ_INT32 sn = dwt->sn;
866 const OPJ_INT32 len = sn + dwt->dn;
868 /* If len == 1, unmodified value */
871 if (len > 1 && nb_cols == PARALLEL_COLS_53) {
872 /* Same as below general case, except that thanks to SSE2 */
873 /* we can efficently process 8 columns in parallel */
874 opj_idwt53_v_cas0_8cols_SSE2(dwt->mem, sn, len, tiledp_col, stride);
880 for (c = 0; c < nb_cols; c++, tiledp_col++) {
881 opj_idwt3_v_cas0(dwt->mem, sn, len, tiledp_col, stride);
888 for (c = 0; c < nb_cols; c++, tiledp_col++) {
896 OPJ_INT32* out = dwt->mem;
897 for (c = 0; c < nb_cols; c++, tiledp_col++) {
899 const OPJ_INT32* in_even = &tiledp_col[sn * stride];
900 const OPJ_INT32* in_odd = &tiledp_col[0];
902 out[1] = in_odd[0] - ((in_even[0] + 1) >> 1);
903 out[0] = in_even[0] + out[1];
905 for (i = 0; i < len; ++i) {
906 tiledp_col[i * stride] = out[i];
914 if (len > 2 && nb_cols == PARALLEL_COLS_53) {
915 /* Same as below general case, except that thanks to SSE2 */
916 /* we can efficently process 8 columns in parallel */
917 opj_idwt53_v_cas1_8cols_SSE2(dwt->mem, sn, len, tiledp_col, stride);
923 for (c = 0; c < nb_cols; c++, tiledp_col++) {
924 opj_idwt3_v_cas1(dwt->mem, sn, len, tiledp_col, stride);
934 /* Forward 9-7 wavelet transform in 1-D. */
936 static void opj_dwt_encode_1_real(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
941 if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
942 for (i = 0; i < dn; i++) {
943 OPJ_D(i) -= opj_int_fix_mul(OPJ_S_(i) + OPJ_S_(i + 1), 12993);
945 for (i = 0; i < sn; i++) {
946 OPJ_S(i) -= opj_int_fix_mul(OPJ_D_(i - 1) + OPJ_D_(i), 434);
948 for (i = 0; i < dn; i++) {
949 OPJ_D(i) += opj_int_fix_mul(OPJ_S_(i) + OPJ_S_(i + 1), 7233);
951 for (i = 0; i < sn; i++) {
952 OPJ_S(i) += opj_int_fix_mul(OPJ_D_(i - 1) + OPJ_D_(i), 3633);
954 for (i = 0; i < dn; i++) {
955 OPJ_D(i) = opj_int_fix_mul(OPJ_D(i), 5038); /*5038 */
957 for (i = 0; i < sn; i++) {
958 OPJ_S(i) = opj_int_fix_mul(OPJ_S(i), 6659); /*6660 */
962 if ((sn > 0) || (dn > 1)) { /* NEW : CASE ONE ELEMENT */
963 for (i = 0; i < dn; i++) {
964 OPJ_S(i) -= opj_int_fix_mul(OPJ_DD_(i) + OPJ_DD_(i - 1), 12993);
966 for (i = 0; i < sn; i++) {
967 OPJ_D(i) -= opj_int_fix_mul(OPJ_SS_(i) + OPJ_SS_(i + 1), 434);
969 for (i = 0; i < dn; i++) {
970 OPJ_S(i) += opj_int_fix_mul(OPJ_DD_(i) + OPJ_DD_(i - 1), 7233);
972 for (i = 0; i < sn; i++) {
973 OPJ_D(i) += opj_int_fix_mul(OPJ_SS_(i) + OPJ_SS_(i + 1), 3633);
975 for (i = 0; i < dn; i++) {
976 OPJ_S(i) = opj_int_fix_mul(OPJ_S(i), 5038); /*5038 */
978 for (i = 0; i < sn; i++) {
979 OPJ_D(i) = opj_int_fix_mul(OPJ_D(i), 6659); /*6660 */
985 static void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps,
986 opj_stepsize_t *bandno_stepsize)
989 p = opj_int_floorlog2(stepsize) - 13;
990 n = 11 - opj_int_floorlog2(stepsize);
991 bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
992 bandno_stepsize->expn = numbps - p;
996 ==========================================================
998 ==========================================================
1003 /* Forward 5-3 wavelet transform in 2-D. */
1005 static INLINE OPJ_BOOL opj_dwt_encode_procedure(opj_tcd_tilecomp_t * tilec,
1006 void (*p_function)(OPJ_INT32 *, OPJ_INT32, OPJ_INT32, OPJ_INT32))
1014 OPJ_INT32 rw; /* width of the resolution level computed */
1015 OPJ_INT32 rh; /* height of the resolution level computed */
1018 opj_tcd_resolution_t * l_cur_res = 0;
1019 opj_tcd_resolution_t * l_last_res = 0;
1021 w = tilec->x1 - tilec->x0;
1022 l = (OPJ_INT32)tilec->numresolutions - 1;
1025 l_cur_res = tilec->resolutions + l;
1026 l_last_res = l_cur_res - 1;
1028 l_data_size = opj_dwt_max_resolution(tilec->resolutions, tilec->numresolutions);
1029 /* overflow check */
1030 if (l_data_size > (SIZE_MAX / sizeof(OPJ_INT32))) {
1031 /* FIXME event manager error callback */
1034 l_data_size *= sizeof(OPJ_INT32);
1035 bj = (OPJ_INT32*)opj_malloc(l_data_size);
1036 /* l_data_size is equal to 0 when numresolutions == 1 but bj is not used */
1037 /* in that case, so do not error out */
1038 if (l_data_size != 0 && ! bj) {
1044 OPJ_INT32 rw1; /* width of the resolution level once lower than computed one */
1045 OPJ_INT32 rh1; /* height of the resolution level once lower than computed one */
1046 OPJ_INT32 cas_col; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
1047 OPJ_INT32 cas_row; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */
1050 rw = l_cur_res->x1 - l_cur_res->x0;
1051 rh = l_cur_res->y1 - l_cur_res->y0;
1052 rw1 = l_last_res->x1 - l_last_res->x0;
1053 rh1 = l_last_res->y1 - l_last_res->y0;
1055 cas_row = l_cur_res->x0 & 1;
1056 cas_col = l_cur_res->y0 & 1;
1060 for (j = 0; j < rw; ++j) {
1062 for (k = 0; k < rh; ++k) {
1066 (*p_function)(bj, dn, sn, cas_col);
1068 opj_dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
1074 for (j = 0; j < rh; j++) {
1076 for (k = 0; k < rw; k++) {
1079 (*p_function)(bj, dn, sn, cas_row);
1080 opj_dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
1083 l_cur_res = l_last_res;
1092 /* Forward 5-3 wavelet transform in 2-D. */
1094 OPJ_BOOL opj_dwt_encode(opj_tcd_tilecomp_t * tilec)
1096 return opj_dwt_encode_procedure(tilec, opj_dwt_encode_1);
1100 /* Inverse 5-3 wavelet transform in 2-D. */
1102 OPJ_BOOL opj_dwt_decode(opj_thread_pool_t* tp, opj_tcd_tilecomp_t* tilec,
1105 return opj_dwt_decode_tile(tp, tilec, numres);
1110 /* Get gain of 5-3 wavelet transform. */
1112 OPJ_UINT32 opj_dwt_getgain(OPJ_UINT32 orient)
1117 if (orient == 1 || orient == 2) {
1124 /* Get norm of 5-3 wavelet. */
1126 OPJ_FLOAT64 opj_dwt_getnorm(OPJ_UINT32 level, OPJ_UINT32 orient)
1128 return opj_dwt_norms[orient][level];
1132 /* Forward 9-7 wavelet transform in 2-D. */
1134 OPJ_BOOL opj_dwt_encode_real(opj_tcd_tilecomp_t * tilec)
1136 return opj_dwt_encode_procedure(tilec, opj_dwt_encode_1_real);
1140 /* Get gain of 9-7 wavelet transform. */
1142 OPJ_UINT32 opj_dwt_getgain_real(OPJ_UINT32 orient)
1149 /* Get norm of 9-7 wavelet. */
1151 OPJ_FLOAT64 opj_dwt_getnorm_real(OPJ_UINT32 level, OPJ_UINT32 orient)
1153 return opj_dwt_norms_real[orient][level];
1156 void opj_dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, OPJ_UINT32 prec)
1158 OPJ_UINT32 numbands, bandno;
1159 numbands = 3 * tccp->numresolutions - 2;
1160 for (bandno = 0; bandno < numbands; bandno++) {
1161 OPJ_FLOAT64 stepsize;
1162 OPJ_UINT32 resno, level, orient, gain;
1164 resno = (bandno == 0) ? 0 : ((bandno - 1) / 3 + 1);
1165 orient = (bandno == 0) ? 0 : ((bandno - 1) % 3 + 1);
1166 level = tccp->numresolutions - 1 - resno;
1167 gain = (tccp->qmfbid == 0) ? 0 : ((orient == 0) ? 0 : (((orient == 1) ||
1168 (orient == 2)) ? 1 : 2));
1169 if (tccp->qntsty == J2K_CCP_QNTSTY_NOQNT) {
1172 OPJ_FLOAT64 norm = opj_dwt_norms_real[orient][level];
1173 stepsize = (1 << (gain)) / norm;
1175 opj_dwt_encode_stepsize((OPJ_INT32) floor(stepsize * 8192.0),
1176 (OPJ_INT32)(prec + gain), &tccp->stepsizes[bandno]);
1181 /* Determine maximum computed resolution level for inverse wavelet transform */
1183 static OPJ_UINT32 opj_dwt_max_resolution(opj_tcd_resolution_t* OPJ_RESTRICT r,
1190 if (mr < (w = (OPJ_UINT32)(r->x1 - r->x0))) {
1193 if (mr < (w = (OPJ_UINT32)(r->y1 - r->y0))) {
1204 OPJ_INT32 * OPJ_RESTRICT tiledp;
1207 } opj_dwd_decode_h_job_t;
1209 static void opj_dwt_decode_h_func(void* user_data, opj_tls_t* tls)
1212 opj_dwd_decode_h_job_t* job;
1215 job = (opj_dwd_decode_h_job_t*)user_data;
1216 for (j = job->min_j; j < job->max_j; j++) {
1217 opj_idwt53_h(&job->h, &job->tiledp[j * job->w]);
1220 opj_aligned_free(job->h.mem);
1228 OPJ_INT32 * OPJ_RESTRICT tiledp;
1231 } opj_dwd_decode_v_job_t;
1233 static void opj_dwt_decode_v_func(void* user_data, opj_tls_t* tls)
1236 opj_dwd_decode_v_job_t* job;
1239 job = (opj_dwd_decode_v_job_t*)user_data;
1240 for (j = job->min_j; j + PARALLEL_COLS_53 <= job->max_j;
1241 j += PARALLEL_COLS_53) {
1242 opj_idwt53_v(&job->v, &job->tiledp[j], (OPJ_INT32)job->w,
1246 opj_idwt53_v(&job->v, &job->tiledp[j], (OPJ_INT32)job->w,
1247 (OPJ_INT32)(job->max_j - j));
1249 opj_aligned_free(job->v.mem);
1255 /* Inverse wavelet transform in 2-D. */
1257 static OPJ_BOOL opj_dwt_decode_tile(opj_thread_pool_t* tp,
1258 opj_tcd_tilecomp_t* tilec, OPJ_UINT32 numres)
1263 opj_tcd_resolution_t* tr = tilec->resolutions;
1265 OPJ_UINT32 rw = (OPJ_UINT32)(tr->x1 -
1266 tr->x0); /* width of the resolution level computed */
1267 OPJ_UINT32 rh = (OPJ_UINT32)(tr->y1 -
1268 tr->y0); /* height of the resolution level computed */
1270 OPJ_UINT32 w = (OPJ_UINT32)(tilec->x1 - tilec->x0);
1277 num_threads = opj_thread_pool_get_thread_count(tp);
1278 h_mem_size = opj_dwt_max_resolution(tr, numres);
1279 /* overflow check */
1280 if (h_mem_size > (SIZE_MAX / PARALLEL_COLS_53 / sizeof(OPJ_INT32))) {
1281 /* FIXME event manager error callback */
1284 /* We need PARALLEL_COLS_53 times the height of the array, */
1285 /* since for the vertical pass */
1286 /* we process PARALLEL_COLS_53 columns at a time */
1287 h_mem_size *= PARALLEL_COLS_53 * sizeof(OPJ_INT32);
1288 h.mem = (OPJ_INT32*)opj_aligned_malloc(h_mem_size);
1290 /* FIXME event manager error callback */
1297 OPJ_INT32 * OPJ_RESTRICT tiledp = tilec->data;
1301 h.sn = (OPJ_INT32)rw;
1302 v.sn = (OPJ_INT32)rh;
1304 rw = (OPJ_UINT32)(tr->x1 - tr->x0);
1305 rh = (OPJ_UINT32)(tr->y1 - tr->y0);
1307 h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
1310 if (num_threads <= 1 || rh <= 1) {
1311 for (j = 0; j < rh; ++j) {
1312 opj_idwt53_h(&h, &tiledp[j * w]);
1315 OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
1318 if (rh < num_jobs) {
1321 step_j = (rh / num_jobs);
1323 for (j = 0; j < num_jobs; j++) {
1324 opj_dwd_decode_h_job_t* job;
1326 job = (opj_dwd_decode_h_job_t*) opj_malloc(sizeof(opj_dwd_decode_h_job_t));
1328 /* It would be nice to fallback to single thread case, but */
1329 /* unfortunately some jobs may be launched and have modified */
1330 /* tiledp, so it is not practical to recover from that error */
1331 /* FIXME event manager error callback */
1332 opj_thread_pool_wait_completion(tp, 0);
1333 opj_aligned_free(h.mem);
1339 job->tiledp = tiledp;
1340 job->min_j = j * step_j;
1341 job->max_j = (j + 1U) * step_j; /* this can overflow */
1342 if (j == (num_jobs - 1U)) { /* this will take care of the overflow */
1345 job->h.mem = (OPJ_INT32*)opj_aligned_malloc(h_mem_size);
1347 /* FIXME event manager error callback */
1348 opj_thread_pool_wait_completion(tp, 0);
1350 opj_aligned_free(h.mem);
1353 opj_thread_pool_submit_job(tp, opj_dwt_decode_h_func, job);
1355 opj_thread_pool_wait_completion(tp, 0);
1358 v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
1361 if (num_threads <= 1 || rw <= 1) {
1362 for (j = 0; j + PARALLEL_COLS_53 <= rw;
1363 j += PARALLEL_COLS_53) {
1364 opj_idwt53_v(&v, &tiledp[j], (OPJ_INT32)w, PARALLEL_COLS_53);
1367 opj_idwt53_v(&v, &tiledp[j], (OPJ_INT32)w, (OPJ_INT32)(rw - j));
1370 OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
1373 if (rw < num_jobs) {
1376 step_j = (rw / num_jobs);
1378 for (j = 0; j < num_jobs; j++) {
1379 opj_dwd_decode_v_job_t* job;
1381 job = (opj_dwd_decode_v_job_t*) opj_malloc(sizeof(opj_dwd_decode_v_job_t));
1383 /* It would be nice to fallback to single thread case, but */
1384 /* unfortunately some jobs may be launched and have modified */
1385 /* tiledp, so it is not practical to recover from that error */
1386 /* FIXME event manager error callback */
1387 opj_thread_pool_wait_completion(tp, 0);
1388 opj_aligned_free(v.mem);
1394 job->tiledp = tiledp;
1395 job->min_j = j * step_j;
1396 job->max_j = (j + 1U) * step_j; /* this can overflow */
1397 if (j == (num_jobs - 1U)) { /* this will take care of the overflow */
1400 job->v.mem = (OPJ_INT32*)opj_aligned_malloc(h_mem_size);
1402 /* FIXME event manager error callback */
1403 opj_thread_pool_wait_completion(tp, 0);
1405 opj_aligned_free(v.mem);
1408 opj_thread_pool_submit_job(tp, opj_dwt_decode_v_func, job);
1410 opj_thread_pool_wait_completion(tp, 0);
1413 opj_aligned_free(h.mem);
1417 static void opj_v4dwt_interleave_h(opj_v4dwt_t* OPJ_RESTRICT w,
1418 OPJ_FLOAT32* OPJ_RESTRICT a, OPJ_INT32 x, OPJ_INT32 size)
1420 OPJ_FLOAT32* OPJ_RESTRICT bi = (OPJ_FLOAT32*)(w->wavelet + w->cas);
1421 OPJ_INT32 count = w->sn;
1424 for (k = 0; k < 2; ++k) {
1425 if (count + 3 * x < size && ((size_t) a & 0x0f) == 0 &&
1426 ((size_t) bi & 0x0f) == 0 && (x & 0x0f) == 0) {
1427 /* Fast code path */
1428 for (i = 0; i < count; ++i) {
1432 bi[i * 8 + 1] = a[j];
1434 bi[i * 8 + 2] = a[j];
1436 bi[i * 8 + 3] = a[j];
1439 /* Slow code path */
1440 for (i = 0; i < count; ++i) {
1447 bi[i * 8 + 1] = a[j];
1452 bi[i * 8 + 2] = a[j];
1457 bi[i * 8 + 3] = a[j]; /* This one*/
1461 bi = (OPJ_FLOAT32*)(w->wavelet + 1 - w->cas);
1468 static void opj_v4dwt_interleave_v(opj_v4dwt_t* OPJ_RESTRICT v,
1469 OPJ_FLOAT32* OPJ_RESTRICT a, OPJ_INT32 x, OPJ_INT32 nb_elts_read)
1471 opj_v4_t* OPJ_RESTRICT bi = v->wavelet + v->cas;
1474 for (i = 0; i < v->sn; ++i) {
1475 memcpy(&bi[i * 2], &a[i * x], (size_t)nb_elts_read * sizeof(OPJ_FLOAT32));
1479 bi = v->wavelet + 1 - v->cas;
1481 for (i = 0; i < v->dn; ++i) {
1482 memcpy(&bi[i * 2], &a[i * x], (size_t)nb_elts_read * sizeof(OPJ_FLOAT32));
1488 static void opj_v4dwt_decode_step1_sse(opj_v4_t* w, OPJ_INT32 count,
1491 __m128* OPJ_RESTRICT vw = (__m128*) w;
1493 /* 4x unrolled loop */
1494 for (i = 0; i < count >> 2; ++i) {
1495 *vw = _mm_mul_ps(*vw, c);
1497 *vw = _mm_mul_ps(*vw, c);
1499 *vw = _mm_mul_ps(*vw, c);
1501 *vw = _mm_mul_ps(*vw, c);
1505 for (i = 0; i < count; ++i) {
1506 *vw = _mm_mul_ps(*vw, c);
1511 void opj_v4dwt_decode_step2_sse(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k,
1512 OPJ_INT32 m, __m128 c)
1514 __m128* OPJ_RESTRICT vl = (__m128*) l;
1515 __m128* OPJ_RESTRICT vw = (__m128*) w;
1517 __m128 tmp1, tmp2, tmp3;
1519 for (i = 0; i < m; ++i) {
1522 vw[-1] = _mm_add_ps(tmp2, _mm_mul_ps(_mm_add_ps(tmp1, tmp3), c));
1530 c = _mm_add_ps(c, c);
1531 c = _mm_mul_ps(c, vl[0]);
1532 for (; m < k; ++m) {
1533 __m128 tmp = vw[-1];
1534 vw[-1] = _mm_add_ps(tmp, c);
1541 static void opj_v4dwt_decode_step1(opj_v4_t* w, OPJ_INT32 count,
1542 const OPJ_FLOAT32 c)
1544 OPJ_FLOAT32* OPJ_RESTRICT fw = (OPJ_FLOAT32*) w;
1546 for (i = 0; i < count; ++i) {
1547 OPJ_FLOAT32 tmp1 = fw[i * 8 ];
1548 OPJ_FLOAT32 tmp2 = fw[i * 8 + 1];
1549 OPJ_FLOAT32 tmp3 = fw[i * 8 + 2];
1550 OPJ_FLOAT32 tmp4 = fw[i * 8 + 3];
1551 fw[i * 8 ] = tmp1 * c;
1552 fw[i * 8 + 1] = tmp2 * c;
1553 fw[i * 8 + 2] = tmp3 * c;
1554 fw[i * 8 + 3] = tmp4 * c;
1558 static void opj_v4dwt_decode_step2(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k,
1559 OPJ_INT32 m, OPJ_FLOAT32 c)
1561 OPJ_FLOAT32* fl = (OPJ_FLOAT32*) l;
1562 OPJ_FLOAT32* fw = (OPJ_FLOAT32*) w;
1564 for (i = 0; i < m; ++i) {
1565 OPJ_FLOAT32 tmp1_1 = fl[0];
1566 OPJ_FLOAT32 tmp1_2 = fl[1];
1567 OPJ_FLOAT32 tmp1_3 = fl[2];
1568 OPJ_FLOAT32 tmp1_4 = fl[3];
1569 OPJ_FLOAT32 tmp2_1 = fw[-4];
1570 OPJ_FLOAT32 tmp2_2 = fw[-3];
1571 OPJ_FLOAT32 tmp2_3 = fw[-2];
1572 OPJ_FLOAT32 tmp2_4 = fw[-1];
1573 OPJ_FLOAT32 tmp3_1 = fw[0];
1574 OPJ_FLOAT32 tmp3_2 = fw[1];
1575 OPJ_FLOAT32 tmp3_3 = fw[2];
1576 OPJ_FLOAT32 tmp3_4 = fw[3];
1577 fw[-4] = tmp2_1 + ((tmp1_1 + tmp3_1) * c);
1578 fw[-3] = tmp2_2 + ((tmp1_2 + tmp3_2) * c);
1579 fw[-2] = tmp2_3 + ((tmp1_3 + tmp3_3) * c);
1580 fw[-1] = tmp2_4 + ((tmp1_4 + tmp3_4) * c);
1594 for (; m < k; ++m) {
1595 OPJ_FLOAT32 tmp1 = fw[-4];
1596 OPJ_FLOAT32 tmp2 = fw[-3];
1597 OPJ_FLOAT32 tmp3 = fw[-2];
1598 OPJ_FLOAT32 tmp4 = fw[-1];
1611 /* Inverse 9-7 wavelet transform in 1-D. */
1613 static void opj_v4dwt_decode(opj_v4dwt_t* OPJ_RESTRICT dwt)
1616 if (dwt->cas == 0) {
1617 if (!((dwt->dn > 0) || (dwt->sn > 1))) {
1623 if (!((dwt->sn > 0) || (dwt->dn > 1))) {
1630 opj_v4dwt_decode_step1_sse(dwt->wavelet + a, dwt->sn, _mm_set1_ps(opj_K));
1631 opj_v4dwt_decode_step1_sse(dwt->wavelet + b, dwt->dn, _mm_set1_ps(opj_c13318));
1632 opj_v4dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
1633 opj_int_min(dwt->sn, dwt->dn - a), _mm_set1_ps(opj_dwt_delta));
1634 opj_v4dwt_decode_step2_sse(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
1635 opj_int_min(dwt->dn, dwt->sn - b), _mm_set1_ps(opj_dwt_gamma));
1636 opj_v4dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
1637 opj_int_min(dwt->sn, dwt->dn - a), _mm_set1_ps(opj_dwt_beta));
1638 opj_v4dwt_decode_step2_sse(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
1639 opj_int_min(dwt->dn, dwt->sn - b), _mm_set1_ps(opj_dwt_alpha));
1641 opj_v4dwt_decode_step1(dwt->wavelet + a, dwt->sn, opj_K);
1642 opj_v4dwt_decode_step1(dwt->wavelet + b, dwt->dn, opj_c13318);
1643 opj_v4dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
1644 opj_int_min(dwt->sn, dwt->dn - a), opj_dwt_delta);
1645 opj_v4dwt_decode_step2(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
1646 opj_int_min(dwt->dn, dwt->sn - b), opj_dwt_gamma);
1647 opj_v4dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
1648 opj_int_min(dwt->sn, dwt->dn - a), opj_dwt_beta);
1649 opj_v4dwt_decode_step2(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
1650 opj_int_min(dwt->dn, dwt->sn - b), opj_dwt_alpha);
1656 /* Inverse 9-7 wavelet transform in 2-D. */
1658 OPJ_BOOL opj_dwt_decode_real(opj_tcd_tilecomp_t* OPJ_RESTRICT tilec,
1664 opj_tcd_resolution_t* res = tilec->resolutions;
1666 OPJ_UINT32 rw = (OPJ_UINT32)(res->x1 -
1667 res->x0); /* width of the resolution level computed */
1668 OPJ_UINT32 rh = (OPJ_UINT32)(res->y1 -
1669 res->y0); /* height of the resolution level computed */
1671 OPJ_UINT32 w = (OPJ_UINT32)(tilec->x1 - tilec->x0);
1675 l_data_size = opj_dwt_max_resolution(res, numres);
1676 /* overflow check */
1677 if (l_data_size > (SIZE_MAX - 5U)) {
1678 /* FIXME event manager error callback */
1682 /* overflow check */
1683 if (l_data_size > (SIZE_MAX / sizeof(opj_v4_t))) {
1684 /* FIXME event manager error callback */
1687 h.wavelet = (opj_v4_t*) opj_aligned_malloc(l_data_size * sizeof(opj_v4_t));
1689 /* FIXME event manager error callback */
1692 v.wavelet = h.wavelet;
1695 OPJ_FLOAT32 * OPJ_RESTRICT aj = (OPJ_FLOAT32*) tilec->data;
1696 OPJ_UINT32 bufsize = (OPJ_UINT32)((tilec->x1 - tilec->x0) *
1697 (tilec->y1 - tilec->y0));
1700 h.sn = (OPJ_INT32)rw;
1701 v.sn = (OPJ_INT32)rh;
1705 rw = (OPJ_UINT32)(res->x1 -
1706 res->x0); /* width of the resolution level computed */
1707 rh = (OPJ_UINT32)(res->y1 -
1708 res->y0); /* height of the resolution level computed */
1710 h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
1711 h.cas = res->x0 % 2;
1713 for (j = (OPJ_INT32)rh; j > 3; j -= 4) {
1715 opj_v4dwt_interleave_h(&h, aj, (OPJ_INT32)w, (OPJ_INT32)bufsize);
1716 opj_v4dwt_decode(&h);
1718 for (k = (OPJ_INT32)rw; --k >= 0;) {
1719 aj[k ] = h.wavelet[k].f[0];
1720 aj[k + (OPJ_INT32)w ] = h.wavelet[k].f[1];
1721 aj[k + (OPJ_INT32)w * 2] = h.wavelet[k].f[2];
1722 aj[k + (OPJ_INT32)w * 3] = h.wavelet[k].f[3];
1732 opj_v4dwt_interleave_h(&h, aj, (OPJ_INT32)w, (OPJ_INT32)bufsize);
1733 opj_v4dwt_decode(&h);
1734 for (k = (OPJ_INT32)rw; --k >= 0;) {
1737 aj[k + (OPJ_INT32)w * 2] = h.wavelet[k].f[2];
1740 aj[k + (OPJ_INT32)w ] = h.wavelet[k].f[1];
1743 aj[k ] = h.wavelet[k].f[0];
1748 v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
1749 v.cas = res->y0 % 2;
1751 aj = (OPJ_FLOAT32*) tilec->data;
1752 for (j = (OPJ_INT32)rw; j > 3; j -= 4) {
1755 opj_v4dwt_interleave_v(&v, aj, (OPJ_INT32)w, 4);
1756 opj_v4dwt_decode(&v);
1758 for (k = 0; k < rh; ++k) {
1759 memcpy(&aj[k * w], &v.wavelet[k], 4 * sizeof(OPJ_FLOAT32));
1769 opj_v4dwt_interleave_v(&v, aj, (OPJ_INT32)w, j);
1770 opj_v4dwt_decode(&v);
1772 for (k = 0; k < rh; ++k) {
1773 memcpy(&aj[k * w], &v.wavelet[k], (size_t)j * sizeof(OPJ_FLOAT32));
1778 opj_aligned_free(h.wavelet);