src/lib/openjp2/dwt.c

   1 /*
   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.
   6  *
   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.
  18  *
  19  * Redistribution and use in source and binary forms, with or without
  20  * modification, are permitted provided that the following conditions
  21  * are met:
  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.
  27  *
  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.
  39  */
  40
  41 #include <assert.h>
  42
  43 #define OPJ_SKIP_POISON
  44 #include "opj_includes.h"
  45
  46 #ifdef __SSE__
  47 #include <xmmintrin.h>
  48 #endif
  49 #ifdef __SSE2__
  50 #include <emmintrin.h>
  51 #endif
  52 #ifdef __SSSE3__
  53 #include <tmmintrin.h>
  54 #endif
  55 #ifdef __AVX2__
  56 #include <immintrin.h>
  57 #endif
  58
  59 #if defined(__GNUC__)
  60 #pragma GCC poison malloc calloc realloc free
  61 #endif
  62
  63 /** @defgroup DWT DWT - Implementation of a discrete wavelet transform */
  64 /*@{*/
  65
  66 #define OPJ_WS(i) v->mem[(i)*2]
  67 #define OPJ_WD(i) v->mem[(1+(i)*2)]
  68
  69 #ifdef __AVX2__
  70 /** Number of int32 values in a AVX2 register */
  71 #define VREG_INT_COUNT       8
  72 #else
  73 /** Number of int32 values in a SSE2 register */
  74 #define VREG_INT_COUNT       4
  75 #endif
  76
  77 /** Number of columns that we can process in parallel in the vertical pass */
  78 #define PARALLEL_COLS_53     (2*VREG_INT_COUNT)
  79
  80 /** @name Local data structures */
  81 /*@{*/
  82
  83 typedef struct dwt_local {
  84     OPJ_INT32* mem;
  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 */
  88 } opj_dwt_t;
  89
  90 #define NB_ELTS_V8  8
  91
  92 typedef union {
  93     OPJ_FLOAT32 f[NB_ELTS_V8];
  94 } opj_v8_t;
  95
  96 typedef struct v8dwt_local {
  97     opj_v8_t*   wavelet ;
  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 */
 105 } opj_v8dwt_t ;
 106
 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;
 112
 113 static const OPJ_FLOAT32 opj_K      = 1.230174105f;
 114 static const OPJ_FLOAT32 opj_invK   = (OPJ_FLOAT32)(1.0 / 1.230174105);
 115
 116 /*@}*/
 117
 118 /**
 119 Virtual function type for wavelet transform in 1-D
 120 */
 121 typedef void (*DWT1DFN)(const opj_dwt_t* v);
 122
 123 /** @name Local static functions */
 124 /*@{*/
 125
 126 /**
 127 Forward lazy transform (horizontal)
 128 */
 129 static void opj_dwt_deinterleave_h(const OPJ_INT32 * OPJ_RESTRICT a,
 130                                    OPJ_INT32 * OPJ_RESTRICT b,
 131                                    OPJ_INT32 dn,
 132                                    OPJ_INT32 sn, OPJ_INT32 cas);
 133 /**
 134 Forward lazy transform (vertical)
 135 */
 136 static void opj_dwt_deinterleave_v(const OPJ_INT32 * OPJ_RESTRICT a,
 137                                    OPJ_INT32 * OPJ_RESTRICT b,
 138                                    OPJ_INT32 dn,
 139                                    OPJ_INT32 sn, OPJ_UINT32 x, OPJ_INT32 cas);
 140 /**
 141 Forward 5-3 wavelet transform in 1-D
 142 */
 143 static void opj_dwt_encode_1(void *a, OPJ_INT32 dn, OPJ_INT32 sn,
 144                              OPJ_INT32 cas);
 145 /**
 146 Forward 9-7 wavelet transform in 1-D
 147 */
 148 static void opj_dwt_encode_1_real(void *a, OPJ_INT32 dn, OPJ_INT32 sn,
 149                                   OPJ_INT32 cas);
 150 /**
 151 Explicit calculation of the Quantization Stepsizes
 152 */
 153 static void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps,
 154                                     opj_stepsize_t *bandno_stepsize);
 155 /**
 156 Inverse wavelet transform in 2-D.
 157 */
 158 static OPJ_BOOL opj_dwt_decode_tile(opj_thread_pool_t* tp,
 159                                     opj_tcd_tilecomp_t* tilec, OPJ_UINT32 i);
 160
 161 static OPJ_BOOL opj_dwt_decode_partial_tile(
 162     opj_tcd_tilecomp_t* tilec,
 163     OPJ_UINT32 numres);
 164
 165 /* Where void* is a OPJ_INT32* for 5x3 and OPJ_FLOAT32* for 9x7 */
 166 typedef void (*opj_encode_one_row_fnptr_type)(void *, OPJ_INT32, OPJ_INT32,
 167         OPJ_INT32);
 168
 169 typedef void (*opj_encode_and_deinterleave_h_one_row_fnptr_type)(
 170     void *row,
 171     void *tmp,
 172     OPJ_UINT32 width,
 173     OPJ_BOOL even);
 174
 175 static OPJ_BOOL opj_dwt_encode_procedure(opj_thread_pool_t* tp,
 176         opj_tcd_tilecomp_t * tilec,
 177         opj_encode_one_row_fnptr_type p_function,
 178         opj_encode_and_deinterleave_h_one_row_fnptr_type
 179         p_encode_and_deinterleave_h_one_row);
 180
 181 static OPJ_UINT32 opj_dwt_max_resolution(opj_tcd_resolution_t* OPJ_RESTRICT r,
 182         OPJ_UINT32 i);
 183
 184 /* <summary>                             */
 185 /* Inverse 9-7 wavelet transform in 1-D. */
 186 /* </summary>                            */
 187
 188 /*@}*/
 189
 190 /*@}*/
 191
 192 #define OPJ_S(i) a[(i)*2]
 193 #define OPJ_D(i) a[(1+(i)*2)]
 194 #define OPJ_S_(i) ((i)<0?OPJ_S(0):((i)>=sn?OPJ_S(sn-1):OPJ_S(i)))
 195 #define OPJ_D_(i) ((i)<0?OPJ_D(0):((i)>=dn?OPJ_D(dn-1):OPJ_D(i)))
 196 /* new */
 197 #define OPJ_SS_(i) ((i)<0?OPJ_S(0):((i)>=dn?OPJ_S(dn-1):OPJ_S(i)))
 198 #define OPJ_DD_(i) ((i)<0?OPJ_D(0):((i)>=sn?OPJ_D(sn-1):OPJ_D(i)))
 199
 200 /* <summary>                                                              */
 201 /* This table contains the norms of the 5-3 wavelets for different bands. */
 202 /* </summary>                                                             */
 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[4][10] = {
 206     {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 207     {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 208     {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 209     {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
 210 };
 211
 212 /* <summary>                                                              */
 213 /* This table contains the norms of the 9-7 wavelets for different bands. */
 214 /* </summary>                                                             */
 215 /* FIXME! the array should really be extended up to 33 resolution levels */
 216 /* See https://github.com/uclouvain/openjpeg/issues/493 */
 217 static const OPJ_FLOAT64 opj_dwt_norms_real[4][10] = {
 218     {1.000, 1.965, 4.177, 8.403, 16.90, 33.84, 67.69, 135.3, 270.6, 540.9},
 219     {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
 220     {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
 221     {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2}
 222 };
 223
 224 /*
 225 ==========================================================
 226    local functions
 227 ==========================================================
 228 */
 229
 230 /* <summary>                             */
 231 /* Forward lazy transform (horizontal).  */
 232 /* </summary>                            */
 233 static void opj_dwt_deinterleave_h(const OPJ_INT32 * OPJ_RESTRICT a,
 234                                    OPJ_INT32 * OPJ_RESTRICT b,
 235                                    OPJ_INT32 dn,
 236                                    OPJ_INT32 sn, OPJ_INT32 cas)
 237 {
 238     OPJ_INT32 i;
 239     OPJ_INT32 * OPJ_RESTRICT l_dest = b;
 240     const OPJ_INT32 * OPJ_RESTRICT l_src = a + cas;
 241
 242     for (i = 0; i < sn; ++i) {
 243         *l_dest++ = *l_src;
 244         l_src += 2;
 245     }
 246
 247     l_dest = b + sn;
 248     l_src = a + 1 - cas;
 249
 250     for (i = 0; i < dn; ++i)  {
 251         *l_dest++ = *l_src;
 252         l_src += 2;
 253     }
 254 }
 255
 256 /* <summary>                             */
 257 /* Forward lazy transform (vertical).    */
 258 /* </summary>                            */
 259 static void opj_dwt_deinterleave_v(const OPJ_INT32 * OPJ_RESTRICT a,
 260                                    OPJ_INT32 * OPJ_RESTRICT b,
 261                                    OPJ_INT32 dn,
 262                                    OPJ_INT32 sn, OPJ_UINT32 x, OPJ_INT32 cas)
 263 {
 264     OPJ_INT32 i = sn;
 265     OPJ_INT32 * OPJ_RESTRICT l_dest = b;
 266     const OPJ_INT32 * OPJ_RESTRICT l_src = a + cas;
 267
 268     while (i--) {
 269         *l_dest = *l_src;
 270         l_dest += x;
 271         l_src += 2;
 272     } /* b[i*x]=a[2*i+cas]; */
 273
 274     l_dest = b + (OPJ_SIZE_T)sn * (OPJ_SIZE_T)x;
 275     l_src = a + 1 - cas;
 276
 277     i = dn;
 278     while (i--) {
 279         *l_dest = *l_src;
 280         l_dest += x;
 281         l_src += 2;
 282     } /*b[(sn+i)*x]=a[(2*i+1-cas)];*/
 283 }
 284
 285 #ifdef STANDARD_SLOW_VERSION
 286 /* <summary>                             */
 287 /* Inverse lazy transform (horizontal).  */
 288 /* </summary>                            */
 289 static void opj_dwt_interleave_h(const opj_dwt_t* h, OPJ_INT32 *a)
 290 {
 291     const OPJ_INT32 *ai = a;
 292     OPJ_INT32 *bi = h->mem + h->cas;
 293     OPJ_INT32  i    = h->sn;
 294     while (i--) {
 295         *bi = *(ai++);
 296         bi += 2;
 297     }
 298     ai  = a + h->sn;
 299     bi  = h->mem + 1 - h->cas;
 300     i   = h->dn ;
 301     while (i--) {
 302         *bi = *(ai++);
 303         bi += 2;
 304     }
 305 }
 306
 307 /* <summary>                             */
 308 /* Inverse lazy transform (vertical).    */
 309 /* </summary>                            */
 310 static void opj_dwt_interleave_v(const opj_dwt_t* v, OPJ_INT32 *a, OPJ_INT32 x)
 311 {
 312     const OPJ_INT32 *ai = a;
 313     OPJ_INT32 *bi = v->mem + v->cas;
 314     OPJ_INT32  i = v->sn;
 315     while (i--) {
 316         *bi = *ai;
 317         bi += 2;
 318         ai += x;
 319     }
 320     ai = a + (v->sn * (OPJ_SIZE_T)x);
 321     bi = v->mem + 1 - v->cas;
 322     i = v->dn ;
 323     while (i--) {
 324         *bi = *ai;
 325         bi += 2;
 326         ai += x;
 327     }
 328 }
 329
 330 #endif /* STANDARD_SLOW_VERSION */
 331
 332 /* <summary>                            */
 333 /* Forward 5-3 wavelet transform in 1-D. */
 334 /* </summary>                           */
 335 static void opj_dwt_encode_1(void *aIn, OPJ_INT32 dn, OPJ_INT32 sn,
 336                              OPJ_INT32 cas)
 337 {
 338     OPJ_INT32 i;
 339     OPJ_INT32* a = (OPJ_INT32*)aIn;
 340
 341     if (!cas) {
 342         if (sn + dn > 1) {
 343             for (i = 0; i < sn - 1; i++) {
 344                 OPJ_D(i) -= (OPJ_S(i) + OPJ_S(i + 1)) >> 1;
 345             }
 346             if (((sn + dn) % 2) == 0) {
 347                 OPJ_D(i) -= OPJ_S(i);
 348             }
 349             OPJ_S(0) += (OPJ_D(0) + OPJ_D(0) + 2) >> 2;
 350             for (i = 1; i < dn; i++) {
 351                 OPJ_S(i) += (OPJ_D(i - 1) + OPJ_D(i) + 2) >> 2;
 352             }
 353             if (((sn + dn) % 2) == 1) {
 354                 OPJ_S(i) += (OPJ_D(i - 1) + OPJ_D(i - 1) + 2) >> 2;
 355             }
 356         }
 357     } else {
 358         if (sn + dn == 1) {
 359             a[0] *= 2;
 360         } else {
 361             OPJ_S(0) -= OPJ_D(0);
 362             for (i = 1; i < sn; i++) {
 363                 OPJ_S(i) -= (OPJ_D(i) + OPJ_D(i - 1)) >> 1;
 364             }
 365             if (((sn + dn) % 2) == 1) {
 366                 OPJ_S(i) -= OPJ_D(i - 1);
 367             }
 368             for (i = 0; i < dn - 1; i++) {
 369                 OPJ_D(i) += (OPJ_S(i) + OPJ_S(i + 1) + 2) >> 2;
 370             }
 371             if (((sn + dn) % 2) == 0) {
 372                 OPJ_D(i) += (OPJ_S(i) + OPJ_S(i) + 2) >> 2;
 373             }
 374         }
 375     }
 376 }
 377
 378 #ifdef STANDARD_SLOW_VERSION
 379 /* <summary>                            */
 380 /* Inverse 5-3 wavelet transform in 1-D. */
 381 /* </summary>                           */
 382 static void opj_dwt_decode_1_(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
 383                               OPJ_INT32 cas)
 384 {
 385     OPJ_INT32 i;
 386
 387     if (!cas) {
 388         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 389             for (i = 0; i < sn; i++) {
 390                 OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
 391             }
 392             for (i = 0; i < dn; i++) {
 393                 OPJ_D(i) += (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
 394             }
 395         }
 396     } else {
 397         if (!sn  && dn == 1) {        /* NEW :  CASE ONE ELEMENT */
 398             OPJ_S(0) /= 2;
 399         } else {
 400             for (i = 0; i < sn; i++) {
 401                 OPJ_D(i) -= (OPJ_SS_(i) + OPJ_SS_(i + 1) + 2) >> 2;
 402             }
 403             for (i = 0; i < dn; i++) {
 404                 OPJ_S(i) += (OPJ_DD_(i) + OPJ_DD_(i - 1)) >> 1;
 405             }
 406         }
 407     }
 408 }
 409
 410 static void opj_dwt_decode_1(const opj_dwt_t *v)
 411 {
 412     opj_dwt_decode_1_(v->mem, v->dn, v->sn, v->cas);
 413 }
 414
 415 #endif /* STANDARD_SLOW_VERSION */
 416
 417 #if !defined(STANDARD_SLOW_VERSION)
 418 static void  opj_idwt53_h_cas0(OPJ_INT32* tmp,
 419                                const OPJ_INT32 sn,
 420                                const OPJ_INT32 len,
 421                                OPJ_INT32* tiledp)
 422 {
 423     OPJ_INT32 i, j;
 424     const OPJ_INT32* in_even = &tiledp[0];
 425     const OPJ_INT32* in_odd = &tiledp[sn];
 426
 427 #ifdef TWO_PASS_VERSION
 428     /* For documentation purpose: performs lifting in two iterations, */
 429     /* but without explicit interleaving */
 430
 431     assert(len > 1);
 432
 433     /* Even */
 434     tmp[0] = in_even[0] - ((in_odd[0] + 1) >> 1);
 435     for (i = 2, j = 0; i <= len - 2; i += 2, j++) {
 436         tmp[i] = in_even[j + 1] - ((in_odd[j] + in_odd[j + 1] + 2) >> 2);
 437     }
 438     if (len & 1) { /* if len is odd */
 439         tmp[len - 1] = in_even[(len - 1) / 2] - ((in_odd[(len - 2) / 2] + 1) >> 1);
 440     }
 441
 442     /* Odd */
 443     for (i = 1, j = 0; i < len - 1; i += 2, j++) {
 444         tmp[i] = in_odd[j] + ((tmp[i - 1] + tmp[i + 1]) >> 1);
 445     }
 446     if (!(len & 1)) { /* if len is even */
 447         tmp[len - 1] = in_odd[(len - 1) / 2] + tmp[len - 2];
 448     }
 449 #else
 450     OPJ_INT32 d1c, d1n, s1n, s0c, s0n;
 451
 452     assert(len > 1);
 453
 454     /* Improved version of the TWO_PASS_VERSION: */
 455     /* Performs lifting in one single iteration. Saves memory */
 456     /* accesses and explicit interleaving. */
 457     s1n = in_even[0];
 458     d1n = in_odd[0];
 459     s0n = s1n - ((d1n + 1) >> 1);
 460
 461     for (i = 0, j = 1; i < (len - 3); i += 2, j++) {
 462         d1c = d1n;
 463         s0c = s0n;
 464
 465         s1n = in_even[j];
 466         d1n = in_odd[j];
 467
 468         s0n = s1n - ((d1c + d1n + 2) >> 2);
 469
 470         tmp[i  ] = s0c;
 471         tmp[i + 1] = d1c + ((s0c + s0n) >> 1);
 472     }
 473
 474     tmp[i] = s0n;
 475
 476     if (len & 1) {
 477         tmp[len - 1] = in_even[(len - 1) / 2] - ((d1n + 1) >> 1);
 478         tmp[len - 2] = d1n + ((s0n + tmp[len - 1]) >> 1);
 479     } else {
 480         tmp[len - 1] = d1n + s0n;
 481     }
 482 #endif
 483     memcpy(tiledp, tmp, (OPJ_UINT32)len * sizeof(OPJ_INT32));
 484 }
 485
 486 static void  opj_idwt53_h_cas1(OPJ_INT32* tmp,
 487                                const OPJ_INT32 sn,
 488                                const OPJ_INT32 len,
 489                                OPJ_INT32* tiledp)
 490 {
 491     OPJ_INT32 i, j;
 492     const OPJ_INT32* in_even = &tiledp[sn];
 493     const OPJ_INT32* in_odd = &tiledp[0];
 494
 495 #ifdef TWO_PASS_VERSION
 496     /* For documentation purpose: performs lifting in two iterations, */
 497     /* but without explicit interleaving */
 498
 499     assert(len > 2);
 500
 501     /* Odd */
 502     for (i = 1, j = 0; i < len - 1; i += 2, j++) {
 503         tmp[i] = in_odd[j] - ((in_even[j] + in_even[j + 1] + 2) >> 2);
 504     }
 505     if (!(len & 1)) {
 506         tmp[len - 1] = in_odd[len / 2 - 1] - ((in_even[len / 2 - 1] + 1) >> 1);
 507     }
 508
 509     /* Even */
 510     tmp[0] = in_even[0] + tmp[1];
 511     for (i = 2, j = 1; i < len - 1; i += 2, j++) {
 512         tmp[i] = in_even[j] + ((tmp[i + 1] + tmp[i - 1]) >> 1);
 513     }
 514     if (len & 1) {
 515         tmp[len - 1] = in_even[len / 2] + tmp[len - 2];
 516     }
 517 #else
 518     OPJ_INT32 s1, s2, dc, dn;
 519
 520     assert(len > 2);
 521
 522     /* Improved version of the TWO_PASS_VERSION: */
 523     /* Performs lifting in one single iteration. Saves memory */
 524     /* accesses and explicit interleaving. */
 525
 526     s1 = in_even[1];
 527     dc = in_odd[0] - ((in_even[0] + s1 + 2) >> 2);
 528     tmp[0] = in_even[0] + dc;
 529
 530     for (i = 1, j = 1; i < (len - 2 - !(len & 1)); i += 2, j++) {
 531
 532         s2 = in_even[j + 1];
 533
 534         dn = in_odd[j] - ((s1 + s2 + 2) >> 2);
 535         tmp[i  ] = dc;
 536         tmp[i + 1] = s1 + ((dn + dc) >> 1);
 537
 538         dc = dn;
 539         s1 = s2;
 540     }
 541
 542     tmp[i] = dc;
 543
 544     if (!(len & 1)) {
 545         dn = in_odd[len / 2 - 1] - ((s1 + 1) >> 1);
 546         tmp[len - 2] = s1 + ((dn + dc) >> 1);
 547         tmp[len - 1] = dn;
 548     } else {
 549         tmp[len - 1] = s1 + dc;
 550     }
 551 #endif
 552     memcpy(tiledp, tmp, (OPJ_UINT32)len * sizeof(OPJ_INT32));
 553 }
 554
 555
 556 #endif /* !defined(STANDARD_SLOW_VERSION) */
 557
 558 /* <summary>                            */
 559 /* Inverse 5-3 wavelet transform in 1-D for one row. */
 560 /* </summary>                           */
 561 /* Performs interleave, inverse wavelet transform and copy back to buffer */
 562 static void opj_idwt53_h(const opj_dwt_t *dwt,
 563                          OPJ_INT32* tiledp)
 564 {
 565 #ifdef STANDARD_SLOW_VERSION
 566     /* For documentation purpose */
 567     opj_dwt_interleave_h(dwt, tiledp);
 568     opj_dwt_decode_1(dwt);
 569     memcpy(tiledp, dwt->mem, (OPJ_UINT32)(dwt->sn + dwt->dn) * sizeof(OPJ_INT32));
 570 #else
 571     const OPJ_INT32 sn = dwt->sn;
 572     const OPJ_INT32 len = sn + dwt->dn;
 573     if (dwt->cas == 0) { /* Left-most sample is on even coordinate */
 574         if (len > 1) {
 575             opj_idwt53_h_cas0(dwt->mem, sn, len, tiledp);
 576         } else {
 577             /* Unmodified value */
 578         }
 579     } else { /* Left-most sample is on odd coordinate */
 580         if (len == 1) {
 581             tiledp[0] /= 2;
 582         } else if (len == 2) {
 583             OPJ_INT32* out = dwt->mem;
 584             const OPJ_INT32* in_even = &tiledp[sn];
 585             const OPJ_INT32* in_odd = &tiledp[0];
 586             out[1] = in_odd[0] - ((in_even[0] + 1) >> 1);
 587             out[0] = in_even[0] + out[1];
 588             memcpy(tiledp, dwt->mem, (OPJ_UINT32)len * sizeof(OPJ_INT32));
 589         } else if (len > 2) {
 590             opj_idwt53_h_cas1(dwt->mem, sn, len, tiledp);
 591         }
 592     }
 593 #endif
 594 }
 595
 596 #if (defined(__SSE2__) || defined(__AVX2__)) && !defined(STANDARD_SLOW_VERSION)
 597
 598 /* Conveniency macros to improve the readabilty of the formulas */
 599 #if __AVX2__
 600 #define VREG        __m256i
 601 #define LOAD_CST(x) _mm256_set1_epi32(x)
 602 #define LOAD(x)     _mm256_load_si256((const VREG*)(x))
 603 #define LOADU(x)    _mm256_loadu_si256((const VREG*)(x))
 604 #define STORE(x,y)  _mm256_store_si256((VREG*)(x),(y))
 605 #define STOREU(x,y) _mm256_storeu_si256((VREG*)(x),(y))
 606 #define ADD(x,y)    _mm256_add_epi32((x),(y))
 607 #define SUB(x,y)    _mm256_sub_epi32((x),(y))
 608 #define SAR(x,y)    _mm256_srai_epi32((x),(y))
 609 #else
 610 #define VREG        __m128i
 611 #define LOAD_CST(x) _mm_set1_epi32(x)
 612 #define LOAD(x)     _mm_load_si128((const VREG*)(x))
 613 #define LOADU(x)    _mm_loadu_si128((const VREG*)(x))
 614 #define STORE(x,y)  _mm_store_si128((VREG*)(x),(y))
 615 #define STOREU(x,y) _mm_storeu_si128((VREG*)(x),(y))
 616 #define ADD(x,y)    _mm_add_epi32((x),(y))
 617 #define SUB(x,y)    _mm_sub_epi32((x),(y))
 618 #define SAR(x,y)    _mm_srai_epi32((x),(y))
 619 #endif
 620 #define ADD3(x,y,z) ADD(ADD(x,y),z)
 621
 622 static
 623 void opj_idwt53_v_final_memcpy(OPJ_INT32* tiledp_col,
 624                                const OPJ_INT32* tmp,
 625                                OPJ_INT32 len,
 626                                OPJ_SIZE_T stride)
 627 {
 628     OPJ_INT32 i;
 629     for (i = 0; i < len; ++i) {
 630         /* A memcpy(&tiledp_col[i * stride + 0],
 631                     &tmp[PARALLEL_COLS_53 * i + 0],
 632                     PARALLEL_COLS_53 * sizeof(OPJ_INT32))
 633            would do but would be a tiny bit slower.
 634            We can take here advantage of our knowledge of alignment */
 635         STOREU(&tiledp_col[(OPJ_SIZE_T)i * stride + 0],
 636                LOAD(&tmp[PARALLEL_COLS_53 * i + 0]));
 637         STOREU(&tiledp_col[(OPJ_SIZE_T)i * stride + VREG_INT_COUNT],
 638                LOAD(&tmp[PARALLEL_COLS_53 * i + VREG_INT_COUNT]));
 639     }
 640 }
 641
 642 /** Vertical inverse 5x3 wavelet transform for 8 columns in SSE2, or
 643  * 16 in AVX2, when top-most pixel is on even coordinate */
 644 static void opj_idwt53_v_cas0_mcols_SSE2_OR_AVX2(
 645     OPJ_INT32* tmp,
 646     const OPJ_INT32 sn,
 647     const OPJ_INT32 len,
 648     OPJ_INT32* tiledp_col,
 649     const OPJ_SIZE_T stride)
 650 {
 651     const OPJ_INT32* in_even = &tiledp_col[0];
 652     const OPJ_INT32* in_odd = &tiledp_col[(OPJ_SIZE_T)sn * stride];
 653
 654     OPJ_INT32 i;
 655     OPJ_SIZE_T j;
 656     VREG d1c_0, d1n_0, s1n_0, s0c_0, s0n_0;
 657     VREG d1c_1, d1n_1, s1n_1, s0c_1, s0n_1;
 658     const VREG two = LOAD_CST(2);
 659
 660     assert(len > 1);
 661 #if __AVX2__
 662     assert(PARALLEL_COLS_53 == 16);
 663     assert(VREG_INT_COUNT == 8);
 664 #else
 665     assert(PARALLEL_COLS_53 == 8);
 666     assert(VREG_INT_COUNT == 4);
 667 #endif
 668
 669     /* Note: loads of input even/odd values must be done in a unaligned */
 670     /* fashion. But stores in tmp can be done with aligned store, since */
 671     /* the temporary buffer is properly aligned */
 672     assert((OPJ_SIZE_T)tmp % (sizeof(OPJ_INT32) * VREG_INT_COUNT) == 0);
 673
 674     s1n_0 = LOADU(in_even + 0);
 675     s1n_1 = LOADU(in_even + VREG_INT_COUNT);
 676     d1n_0 = LOADU(in_odd);
 677     d1n_1 = LOADU(in_odd + VREG_INT_COUNT);
 678
 679     /* s0n = s1n - ((d1n + 1) >> 1); <==> */
 680     /* s0n = s1n - ((d1n + d1n + 2) >> 2); */
 681     s0n_0 = SUB(s1n_0, SAR(ADD3(d1n_0, d1n_0, two), 2));
 682     s0n_1 = SUB(s1n_1, SAR(ADD3(d1n_1, d1n_1, two), 2));
 683
 684     for (i = 0, j = 1; i < (len - 3); i += 2, j++) {
 685         d1c_0 = d1n_0;
 686         s0c_0 = s0n_0;
 687         d1c_1 = d1n_1;
 688         s0c_1 = s0n_1;
 689
 690         s1n_0 = LOADU(in_even + j * stride);
 691         s1n_1 = LOADU(in_even + j * stride + VREG_INT_COUNT);
 692         d1n_0 = LOADU(in_odd + j * stride);
 693         d1n_1 = LOADU(in_odd + j * stride + VREG_INT_COUNT);
 694
 695         /*s0n = s1n - ((d1c + d1n + 2) >> 2);*/
 696         s0n_0 = SUB(s1n_0, SAR(ADD3(d1c_0, d1n_0, two), 2));
 697         s0n_1 = SUB(s1n_1, SAR(ADD3(d1c_1, d1n_1, two), 2));
 698
 699         STORE(tmp + PARALLEL_COLS_53 * (i + 0), s0c_0);
 700         STORE(tmp + PARALLEL_COLS_53 * (i + 0) + VREG_INT_COUNT, s0c_1);
 701
 702         /* d1c + ((s0c + s0n) >> 1) */
 703         STORE(tmp + PARALLEL_COLS_53 * (i + 1) + 0,
 704               ADD(d1c_0, SAR(ADD(s0c_0, s0n_0), 1)));
 705         STORE(tmp + PARALLEL_COLS_53 * (i + 1) + VREG_INT_COUNT,
 706               ADD(d1c_1, SAR(ADD(s0c_1, s0n_1), 1)));
 707     }
 708
 709     STORE(tmp + PARALLEL_COLS_53 * (i + 0) + 0, s0n_0);
 710     STORE(tmp + PARALLEL_COLS_53 * (i + 0) + VREG_INT_COUNT, s0n_1);
 711
 712     if (len & 1) {
 713         VREG tmp_len_minus_1;
 714         s1n_0 = LOADU(in_even + (OPJ_SIZE_T)((len - 1) / 2) * stride);
 715         /* tmp_len_minus_1 = s1n - ((d1n + 1) >> 1); */
 716         tmp_len_minus_1 = SUB(s1n_0, SAR(ADD3(d1n_0, d1n_0, two), 2));
 717         STORE(tmp + PARALLEL_COLS_53 * (len - 1), tmp_len_minus_1);
 718         /* d1n + ((s0n + tmp_len_minus_1) >> 1) */
 719         STORE(tmp + PARALLEL_COLS_53 * (len - 2),
 720               ADD(d1n_0, SAR(ADD(s0n_0, tmp_len_minus_1), 1)));
 721
 722         s1n_1 = LOADU(in_even + (OPJ_SIZE_T)((len - 1) / 2) * stride + VREG_INT_COUNT);
 723         /* tmp_len_minus_1 = s1n - ((d1n + 1) >> 1); */
 724         tmp_len_minus_1 = SUB(s1n_1, SAR(ADD3(d1n_1, d1n_1, two), 2));
 725         STORE(tmp + PARALLEL_COLS_53 * (len - 1) + VREG_INT_COUNT,
 726               tmp_len_minus_1);
 727         /* d1n + ((s0n + tmp_len_minus_1) >> 1) */
 728         STORE(tmp + PARALLEL_COLS_53 * (len - 2) + VREG_INT_COUNT,
 729               ADD(d1n_1, SAR(ADD(s0n_1, tmp_len_minus_1), 1)));
 730
 731
 732     } else {
 733         STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 0,
 734               ADD(d1n_0, s0n_0));
 735         STORE(tmp + PARALLEL_COLS_53 * (len - 1) + VREG_INT_COUNT,
 736               ADD(d1n_1, s0n_1));
 737     }
 738
 739     opj_idwt53_v_final_memcpy(tiledp_col, tmp, len, stride);
 740 }
 741
 742
 743 /** Vertical inverse 5x3 wavelet transform for 8 columns in SSE2, or
 744  * 16 in AVX2, when top-most pixel is on odd coordinate */
 745 static void opj_idwt53_v_cas1_mcols_SSE2_OR_AVX2(
 746     OPJ_INT32* tmp,
 747     const OPJ_INT32 sn,
 748     const OPJ_INT32 len,
 749     OPJ_INT32* tiledp_col,
 750     const OPJ_SIZE_T stride)
 751 {
 752     OPJ_INT32 i;
 753     OPJ_SIZE_T j;
 754
 755     VREG s1_0, s2_0, dc_0, dn_0;
 756     VREG s1_1, s2_1, dc_1, dn_1;
 757     const VREG two = LOAD_CST(2);
 758
 759     const OPJ_INT32* in_even = &tiledp_col[(OPJ_SIZE_T)sn * stride];
 760     const OPJ_INT32* in_odd = &tiledp_col[0];
 761
 762     assert(len > 2);
 763 #if __AVX2__
 764     assert(PARALLEL_COLS_53 == 16);
 765     assert(VREG_INT_COUNT == 8);
 766 #else
 767     assert(PARALLEL_COLS_53 == 8);
 768     assert(VREG_INT_COUNT == 4);
 769 #endif
 770
 771     /* Note: loads of input even/odd values must be done in a unaligned */
 772     /* fashion. But stores in tmp can be done with aligned store, since */
 773     /* the temporary buffer is properly aligned */
 774     assert((OPJ_SIZE_T)tmp % (sizeof(OPJ_INT32) * VREG_INT_COUNT) == 0);
 775
 776     s1_0 = LOADU(in_even + stride);
 777     /* in_odd[0] - ((in_even[0] + s1 + 2) >> 2); */
 778     dc_0 = SUB(LOADU(in_odd + 0),
 779                SAR(ADD3(LOADU(in_even + 0), s1_0, two), 2));
 780     STORE(tmp + PARALLEL_COLS_53 * 0, ADD(LOADU(in_even + 0), dc_0));
 781
 782     s1_1 = LOADU(in_even + stride + VREG_INT_COUNT);
 783     /* in_odd[0] - ((in_even[0] + s1 + 2) >> 2); */
 784     dc_1 = SUB(LOADU(in_odd + VREG_INT_COUNT),
 785                SAR(ADD3(LOADU(in_even + VREG_INT_COUNT), s1_1, two), 2));
 786     STORE(tmp + PARALLEL_COLS_53 * 0 + VREG_INT_COUNT,
 787           ADD(LOADU(in_even + VREG_INT_COUNT), dc_1));
 788
 789     for (i = 1, j = 1; i < (len - 2 - !(len & 1)); i += 2, j++) {
 790
 791         s2_0 = LOADU(in_even + (j + 1) * stride);
 792         s2_1 = LOADU(in_even + (j + 1) * stride + VREG_INT_COUNT);
 793
 794         /* dn = in_odd[j * stride] - ((s1 + s2 + 2) >> 2); */
 795         dn_0 = SUB(LOADU(in_odd + j * stride),
 796                    SAR(ADD3(s1_0, s2_0, two), 2));
 797         dn_1 = SUB(LOADU(in_odd + j * stride + VREG_INT_COUNT),
 798                    SAR(ADD3(s1_1, s2_1, two), 2));
 799
 800         STORE(tmp + PARALLEL_COLS_53 * i, dc_0);
 801         STORE(tmp + PARALLEL_COLS_53 * i + VREG_INT_COUNT, dc_1);
 802
 803         /* tmp[i + 1] = s1 + ((dn + dc) >> 1); */
 804         STORE(tmp + PARALLEL_COLS_53 * (i + 1) + 0,
 805               ADD(s1_0, SAR(ADD(dn_0, dc_0), 1)));
 806         STORE(tmp + PARALLEL_COLS_53 * (i + 1) + VREG_INT_COUNT,
 807               ADD(s1_1, SAR(ADD(dn_1, dc_1), 1)));
 808
 809         dc_0 = dn_0;
 810         s1_0 = s2_0;
 811         dc_1 = dn_1;
 812         s1_1 = s2_1;
 813     }
 814     STORE(tmp + PARALLEL_COLS_53 * i, dc_0);
 815     STORE(tmp + PARALLEL_COLS_53 * i + VREG_INT_COUNT, dc_1);
 816
 817     if (!(len & 1)) {
 818         /*dn = in_odd[(len / 2 - 1) * stride] - ((s1 + 1) >> 1); */
 819         dn_0 = SUB(LOADU(in_odd + (OPJ_SIZE_T)(len / 2 - 1) * stride),
 820                    SAR(ADD3(s1_0, s1_0, two), 2));
 821         dn_1 = SUB(LOADU(in_odd + (OPJ_SIZE_T)(len / 2 - 1) * stride + VREG_INT_COUNT),
 822                    SAR(ADD3(s1_1, s1_1, two), 2));
 823
 824         /* tmp[len - 2] = s1 + ((dn + dc) >> 1); */
 825         STORE(tmp + PARALLEL_COLS_53 * (len - 2) + 0,
 826               ADD(s1_0, SAR(ADD(dn_0, dc_0), 1)));
 827         STORE(tmp + PARALLEL_COLS_53 * (len - 2) + VREG_INT_COUNT,
 828               ADD(s1_1, SAR(ADD(dn_1, dc_1), 1)));
 829
 830         STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 0, dn_0);
 831         STORE(tmp + PARALLEL_COLS_53 * (len - 1) + VREG_INT_COUNT, dn_1);
 832     } else {
 833         STORE(tmp + PARALLEL_COLS_53 * (len - 1) + 0, ADD(s1_0, dc_0));
 834         STORE(tmp + PARALLEL_COLS_53 * (len - 1) + VREG_INT_COUNT,
 835               ADD(s1_1, dc_1));
 836     }
 837
 838     opj_idwt53_v_final_memcpy(tiledp_col, tmp, len, stride);
 839 }
 840
 841 #undef VREG
 842 #undef LOAD_CST
 843 #undef LOADU
 844 #undef LOAD
 845 #undef STORE
 846 #undef STOREU
 847 #undef ADD
 848 #undef ADD3
 849 #undef SUB
 850 #undef SAR
 851
 852 #endif /* (defined(__SSE2__) || defined(__AVX2__)) && !defined(STANDARD_SLOW_VERSION) */
 853
 854 #if !defined(STANDARD_SLOW_VERSION)
 855 /** Vertical inverse 5x3 wavelet transform for one column, when top-most
 856  * pixel is on even coordinate */
 857 static void opj_idwt3_v_cas0(OPJ_INT32* tmp,
 858                              const OPJ_INT32 sn,
 859                              const OPJ_INT32 len,
 860                              OPJ_INT32* tiledp_col,
 861                              const OPJ_SIZE_T stride)
 862 {
 863     OPJ_INT32 i, j;
 864     OPJ_INT32 d1c, d1n, s1n, s0c, s0n;
 865
 866     assert(len > 1);
 867
 868     /* Performs lifting in one single iteration. Saves memory */
 869     /* accesses and explicit interleaving. */
 870
 871     s1n = tiledp_col[0];
 872     d1n = tiledp_col[(OPJ_SIZE_T)sn * stride];
 873     s0n = s1n - ((d1n + 1) >> 1);
 874
 875     for (i = 0, j = 0; i < (len - 3); i += 2, j++) {
 876         d1c = d1n;
 877         s0c = s0n;
 878
 879         s1n = tiledp_col[(OPJ_SIZE_T)(j + 1) * stride];
 880         d1n = tiledp_col[(OPJ_SIZE_T)(sn + j + 1) * stride];
 881
 882         s0n = s1n - ((d1c + d1n + 2) >> 2);
 883
 884         tmp[i  ] = s0c;
 885         tmp[i + 1] = d1c + ((s0c + s0n) >> 1);
 886     }
 887
 888     tmp[i] = s0n;
 889
 890     if (len & 1) {
 891         tmp[len - 1] =
 892             tiledp_col[(OPJ_SIZE_T)((len - 1) / 2) * stride] -
 893             ((d1n + 1) >> 1);
 894         tmp[len - 2] = d1n + ((s0n + tmp[len - 1]) >> 1);
 895     } else {
 896         tmp[len - 1] = d1n + s0n;
 897     }
 898
 899     for (i = 0; i < len; ++i) {
 900         tiledp_col[(OPJ_SIZE_T)i * stride] = tmp[i];
 901     }
 902 }
 903
 904
 905 /** Vertical inverse 5x3 wavelet transform for one column, when top-most
 906  * pixel is on odd coordinate */
 907 static void opj_idwt3_v_cas1(OPJ_INT32* tmp,
 908                              const OPJ_INT32 sn,
 909                              const OPJ_INT32 len,
 910                              OPJ_INT32* tiledp_col,
 911                              const OPJ_SIZE_T stride)
 912 {
 913     OPJ_INT32 i, j;
 914     OPJ_INT32 s1, s2, dc, dn;
 915     const OPJ_INT32* in_even = &tiledp_col[(OPJ_SIZE_T)sn * stride];
 916     const OPJ_INT32* in_odd = &tiledp_col[0];
 917
 918     assert(len > 2);
 919
 920     /* Performs lifting in one single iteration. Saves memory */
 921     /* accesses and explicit interleaving. */
 922
 923     s1 = in_even[stride];
 924     dc = in_odd[0] - ((in_even[0] + s1 + 2) >> 2);
 925     tmp[0] = in_even[0] + dc;
 926     for (i = 1, j = 1; i < (len - 2 - !(len & 1)); i += 2, j++) {
 927
 928         s2 = in_even[(OPJ_SIZE_T)(j + 1) * stride];
 929
 930         dn = in_odd[(OPJ_SIZE_T)j * stride] - ((s1 + s2 + 2) >> 2);
 931         tmp[i  ] = dc;
 932         tmp[i + 1] = s1 + ((dn + dc) >> 1);
 933
 934         dc = dn;
 935         s1 = s2;
 936     }
 937     tmp[i] = dc;
 938     if (!(len & 1)) {
 939         dn = in_odd[(OPJ_SIZE_T)(len / 2 - 1) * stride] - ((s1 + 1) >> 1);
 940         tmp[len - 2] = s1 + ((dn + dc) >> 1);
 941         tmp[len - 1] = dn;
 942     } else {
 943         tmp[len - 1] = s1 + dc;
 944     }
 945
 946     for (i = 0; i < len; ++i) {
 947         tiledp_col[(OPJ_SIZE_T)i * stride] = tmp[i];
 948     }
 949 }
 950 #endif /* !defined(STANDARD_SLOW_VERSION) */
 951
 952 /* <summary>                            */
 953 /* Inverse vertical 5-3 wavelet transform in 1-D for several columns. */
 954 /* </summary>                           */
 955 /* Performs interleave, inverse wavelet transform and copy back to buffer */
 956 static void opj_idwt53_v(const opj_dwt_t *dwt,
 957                          OPJ_INT32* tiledp_col,
 958                          OPJ_SIZE_T stride,
 959                          OPJ_INT32 nb_cols)
 960 {
 961 #ifdef STANDARD_SLOW_VERSION
 962     /* For documentation purpose */
 963     OPJ_INT32 k, c;
 964     for (c = 0; c < nb_cols; c ++) {
 965         opj_dwt_interleave_v(dwt, tiledp_col + c, stride);
 966         opj_dwt_decode_1(dwt);
 967         for (k = 0; k < dwt->sn + dwt->dn; ++k) {
 968             tiledp_col[c + k * stride] = dwt->mem[k];
 969         }
 970     }
 971 #else
 972     const OPJ_INT32 sn = dwt->sn;
 973     const OPJ_INT32 len = sn + dwt->dn;
 974     if (dwt->cas == 0) {
 975         /* If len == 1, unmodified value */
 976
 977 #if (defined(__SSE2__) || defined(__AVX2__))
 978         if (len > 1 && nb_cols == PARALLEL_COLS_53) {
 979             /* Same as below general case, except that thanks to SSE2/AVX2 */
 980             /* we can efficiently process 8/16 columns in parallel */
 981             opj_idwt53_v_cas0_mcols_SSE2_OR_AVX2(dwt->mem, sn, len, tiledp_col, stride);
 982             return;
 983         }
 984 #endif
 985         if (len > 1) {
 986             OPJ_INT32 c;
 987             for (c = 0; c < nb_cols; c++, tiledp_col++) {
 988                 opj_idwt3_v_cas0(dwt->mem, sn, len, tiledp_col, stride);
 989             }
 990             return;
 991         }
 992     } else {
 993         if (len == 1) {
 994             OPJ_INT32 c;
 995             for (c = 0; c < nb_cols; c++, tiledp_col++) {
 996                 tiledp_col[0] /= 2;
 997             }
 998             return;
 999         }
1000
1001         if (len == 2) {
1002             OPJ_INT32 c;
1003             OPJ_INT32* out = dwt->mem;
1004             for (c = 0; c < nb_cols; c++, tiledp_col++) {
1005                 OPJ_INT32 i;
1006                 const OPJ_INT32* in_even = &tiledp_col[(OPJ_SIZE_T)sn * stride];
1007                 const OPJ_INT32* in_odd = &tiledp_col[0];
1008
1009                 out[1] = in_odd[0] - ((in_even[0] + 1) >> 1);
1010                 out[0] = in_even[0] + out[1];
1011
1012                 for (i = 0; i < len; ++i) {
1013                     tiledp_col[(OPJ_SIZE_T)i * stride] = out[i];
1014                 }
1015             }
1016
1017             return;
1018         }
1019
1020 #if (defined(__SSE2__) || defined(__AVX2__))
1021         if (len > 2 && nb_cols == PARALLEL_COLS_53) {
1022             /* Same as below general case, except that thanks to SSE2/AVX2 */
1023             /* we can efficiently process 8/16 columns in parallel */
1024             opj_idwt53_v_cas1_mcols_SSE2_OR_AVX2(dwt->mem, sn, len, tiledp_col, stride);
1025             return;
1026         }
1027 #endif
1028         if (len > 2) {
1029             OPJ_INT32 c;
1030             for (c = 0; c < nb_cols; c++, tiledp_col++) {
1031                 opj_idwt3_v_cas1(dwt->mem, sn, len, tiledp_col, stride);
1032             }
1033             return;
1034         }
1035     }
1036 #endif
1037 }
1038
1039 static void opj_dwt_encode_step1(OPJ_FLOAT32* fw,
1040                                  OPJ_UINT32 start,
1041                                  OPJ_UINT32 end,
1042                                  const OPJ_FLOAT32 c)
1043 {
1044     OPJ_UINT32 i;
1045     for (i = start; i < end; ++i) {
1046         fw[i * 2] *= c;
1047     }
1048 }
1049 static void opj_dwt_encode_step2(OPJ_FLOAT32* fl, OPJ_FLOAT32* fw,
1050                                  OPJ_UINT32 start,
1051                                  OPJ_UINT32 end,
1052                                  OPJ_UINT32 m,
1053                                  OPJ_FLOAT32 c)
1054 {
1055     OPJ_UINT32 i;
1056     OPJ_UINT32 imax = opj_uint_min(end, m);
1057     if (start > 0) {
1058         fw += 2 * start;
1059         fl = fw - 2;
1060     }
1061     for (i = start; i < imax; ++i) {
1062         fw[-1] += (fl[0] + fw[0]) * c;
1063         fl = fw;
1064         fw += 2;
1065     }
1066     if (m < end) {
1067         assert(m + 1 == end);
1068         fw[-1] += (2 * fl[0]) * c;
1069     }
1070 }
1071
1072 static void opj_dwt_encode_1_real(void *aIn, OPJ_INT32 dn, OPJ_INT32 sn,
1073                                   OPJ_INT32 cas)
1074 {
1075     OPJ_FLOAT32* w = (OPJ_FLOAT32*)aIn;
1076     OPJ_INT32 a, b;
1077     if (cas == 0) {
1078         if (!((dn > 0) || (sn > 1))) {
1079             return;
1080         }
1081         a = 0;
1082         b = 1;
1083     } else {
1084         if (!((sn > 0) || (dn > 1))) {
1085             return;
1086         }
1087         a = 1;
1088         b = 0;
1089     }
1090     opj_dwt_encode_step2(w + a, w + b + 1,
1091                          0, (OPJ_UINT32)dn,
1092                          (OPJ_UINT32)opj_int_min(dn, sn - b),
1093                          opj_dwt_alpha);
1094     opj_dwt_encode_step2(w + b, w + a + 1,
1095                          0, (OPJ_UINT32)sn,
1096                          (OPJ_UINT32)opj_int_min(sn, dn - a),
1097                          opj_dwt_beta);
1098     opj_dwt_encode_step2(w + a, w + b + 1,
1099                          0, (OPJ_UINT32)dn,
1100                          (OPJ_UINT32)opj_int_min(dn, sn - b),
1101                          opj_dwt_gamma);
1102     opj_dwt_encode_step2(w + b, w + a + 1,
1103                          0, (OPJ_UINT32)sn,
1104                          (OPJ_UINT32)opj_int_min(sn, dn - a),
1105                          opj_dwt_delta);
1106     opj_dwt_encode_step1(w + b, 0, (OPJ_UINT32)dn,
1107                          opj_K);
1108     opj_dwt_encode_step1(w + a, 0, (OPJ_UINT32)sn,
1109                          opj_invK);
1110 }
1111
1112 static void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps,
1113                                     opj_stepsize_t *bandno_stepsize)
1114 {
1115     OPJ_INT32 p, n;
1116     p = opj_int_floorlog2(stepsize) - 13;
1117     n = 11 - opj_int_floorlog2(stepsize);
1118     bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
1119     bandno_stepsize->expn = numbps - p;
1120 }
1121
1122 /*
1123 ==========================================================
1124    DWT interface
1125 ==========================================================
1126 */
1127
1128 /** Process one line for the horizontal pass of the 5x3 forward transform */
1129 static
1130 void opj_dwt_encode_and_deinterleave_h_one_row(void* rowIn,
1131         void* tmpIn,
1132         OPJ_UINT32 width,
1133         OPJ_BOOL even)
1134 {
1135     OPJ_INT32* OPJ_RESTRICT row = (OPJ_INT32*)rowIn;
1136     OPJ_INT32* OPJ_RESTRICT tmp = (OPJ_INT32*)tmpIn;
1137     const OPJ_INT32 sn = (OPJ_INT32)((width + (even ? 1 : 0)) >> 1);
1138     const OPJ_INT32 dn = (OPJ_INT32)(width - (OPJ_UINT32)sn);
1139
1140     if (even) {
1141         if (width > 1) {
1142             OPJ_INT32 i;
1143             for (i = 0; i < sn - 1; i++) {
1144                 tmp[sn + i] = row[2 * i + 1] - ((row[(i) * 2] + row[(i + 1) * 2]) >> 1);
1145             }
1146             if ((width % 2) == 0) {
1147                 tmp[sn + i] = row[2 * i + 1] - row[(i) * 2];
1148             }
1149             row[0] += (tmp[sn] + tmp[sn] + 2) >> 2;
1150             for (i = 1; i < dn; i++) {
1151                 row[i] = row[2 * i] + ((tmp[sn + (i - 1)] + tmp[sn + i] + 2) >> 2);
1152             }
1153             if ((width % 2) == 1) {
1154                 row[i] = row[2 * i] + ((tmp[sn + (i - 1)] + tmp[sn + (i - 1)] + 2) >> 2);
1155             }
1156             memcpy(row + sn, tmp + sn, (OPJ_SIZE_T)dn * sizeof(OPJ_INT32));
1157         }
1158     } else {
1159         if (width == 1) {
1160             row[0] *= 2;
1161         } else {
1162             OPJ_INT32 i;
1163             tmp[sn + 0] = row[0] - row[1];
1164             for (i = 1; i < sn; i++) {
1165                 tmp[sn + i] = row[2 * i] - ((row[2 * i + 1] + row[2 * (i - 1) + 1]) >> 1);
1166             }
1167             if ((width % 2) == 1) {
1168                 tmp[sn + i] = row[2 * i] - row[2 * (i - 1) + 1];
1169             }
1170
1171             for (i = 0; i < dn - 1; i++) {
1172                 row[i] = row[2 * i + 1] + ((tmp[sn + i] + tmp[sn + i + 1] + 2) >> 2);
1173             }
1174             if ((width % 2) == 0) {
1175                 row[i] = row[2 * i + 1] + ((tmp[sn + i] + tmp[sn + i] + 2) >> 2);
1176             }
1177             memcpy(row + sn, tmp + sn, (OPJ_SIZE_T)dn * sizeof(OPJ_INT32));
1178         }
1179     }
1180 }
1181
1182 /** Process one line for the horizontal pass of the 9x7 forward transform */
1183 static
1184 void opj_dwt_encode_and_deinterleave_h_one_row_real(void* rowIn,
1185         void* tmpIn,
1186         OPJ_UINT32 width,
1187         OPJ_BOOL even)
1188 {
1189     OPJ_FLOAT32* OPJ_RESTRICT row = (OPJ_FLOAT32*)rowIn;
1190     OPJ_FLOAT32* OPJ_RESTRICT tmp = (OPJ_FLOAT32*)tmpIn;
1191     const OPJ_INT32 sn = (OPJ_INT32)((width + (even ? 1 : 0)) >> 1);
1192     const OPJ_INT32 dn = (OPJ_INT32)(width - (OPJ_UINT32)sn);
1193     memcpy(tmp, row, width * sizeof(OPJ_FLOAT32));
1194     opj_dwt_encode_1_real(tmp, dn, sn, even ? 0 : 1);
1195     opj_dwt_deinterleave_h((OPJ_INT32 * OPJ_RESTRICT)tmp,
1196                            (OPJ_INT32 * OPJ_RESTRICT)row,
1197                            dn, sn, even ? 0 : 1);
1198 }
1199
1200 typedef struct {
1201     opj_dwt_t h;
1202     OPJ_UINT32 rw; /* Width of the resolution to process */
1203     OPJ_UINT32 w; /* Width of tiledp */
1204     OPJ_INT32 * OPJ_RESTRICT tiledp;
1205     OPJ_UINT32 min_j;
1206     OPJ_UINT32 max_j;
1207     opj_encode_and_deinterleave_h_one_row_fnptr_type p_function;
1208 } opj_dwt_encode_h_job_t;
1209
1210 static void opj_dwt_encode_h_func(void* user_data, opj_tls_t* tls)
1211 {
1212     OPJ_UINT32 j;
1213     opj_dwt_encode_h_job_t* job;
1214     (void)tls;
1215
1216     job = (opj_dwt_encode_h_job_t*)user_data;
1217     for (j = job->min_j; j < job->max_j; j++) {
1218         OPJ_INT32* OPJ_RESTRICT aj = job->tiledp + j * job->w;
1219         (*job->p_function)(aj, job->h.mem, job->rw,
1220                            job->h.cas == 0 ? OPJ_TRUE : OPJ_FALSE);
1221     }
1222
1223     opj_aligned_free(job->h.mem);
1224     opj_free(job);
1225 }
1226
1227 typedef struct {
1228     opj_dwt_t v;
1229     OPJ_UINT32 rh;
1230     OPJ_UINT32 w;
1231     OPJ_INT32 * OPJ_RESTRICT tiledp;
1232     OPJ_UINT32 min_j;
1233     OPJ_UINT32 max_j;
1234     opj_encode_one_row_fnptr_type p_function;
1235 } opj_dwt_encode_v_job_t;
1236
1237 static void opj_dwt_encode_v_func(void* user_data, opj_tls_t* tls)
1238 {
1239     OPJ_UINT32 j;
1240     opj_dwt_encode_v_job_t* job;
1241     (void)tls;
1242
1243     job = (opj_dwt_encode_v_job_t*)user_data;
1244     for (j = job->min_j; j < job->max_j; j++) {
1245         OPJ_INT32* OPJ_RESTRICT aj = job->tiledp + j;
1246         OPJ_UINT32 k;
1247         for (k = 0; k < job->rh; ++k) {
1248             job->v.mem[k] = aj[k * job->w];
1249         }
1250
1251         (*job->p_function)(job->v.mem, job->v.dn, job->v.sn, job->v.cas);
1252
1253         opj_dwt_deinterleave_v(job->v.mem, aj, job->v.dn, job->v.sn, job->w,
1254                                job->v.cas);
1255     }
1256
1257     opj_aligned_free(job->v.mem);
1258     opj_free(job);
1259 }
1260
1261 /* <summary>                            */
1262 /* Forward 5-3 wavelet transform in 2-D. */
1263 /* </summary>                           */
1264 static INLINE OPJ_BOOL opj_dwt_encode_procedure(opj_thread_pool_t* tp,
1265         opj_tcd_tilecomp_t * tilec,
1266         opj_encode_one_row_fnptr_type p_function,
1267         opj_encode_and_deinterleave_h_one_row_fnptr_type
1268         p_encode_and_deinterleave_h_one_row)
1269 {
1270     OPJ_INT32 i;
1271     OPJ_INT32 *bj = 00;
1272     OPJ_UINT32 w;
1273     OPJ_INT32 l;
1274
1275     OPJ_SIZE_T l_data_size;
1276
1277     opj_tcd_resolution_t * l_cur_res = 0;
1278     opj_tcd_resolution_t * l_last_res = 0;
1279     const int num_threads = opj_thread_pool_get_thread_count(tp);
1280     OPJ_INT32 * OPJ_RESTRICT tiledp = tilec->data;
1281
1282     w = (OPJ_UINT32)(tilec->x1 - tilec->x0);
1283     l = (OPJ_INT32)tilec->numresolutions - 1;
1284
1285     l_cur_res = tilec->resolutions + l;
1286     l_last_res = l_cur_res - 1;
1287
1288     l_data_size = opj_dwt_max_resolution(tilec->resolutions, tilec->numresolutions);
1289     /* overflow check */
1290     if (l_data_size > (SIZE_MAX / sizeof(OPJ_INT32))) {
1291         /* FIXME event manager error callback */
1292         return OPJ_FALSE;
1293     }
1294     l_data_size *= sizeof(OPJ_INT32);
1295     bj = (OPJ_INT32*)opj_aligned_32_malloc(l_data_size);
1296     /* l_data_size is equal to 0 when numresolutions == 1 but bj is not used */
1297     /* in that case, so do not error out */
1298     if (l_data_size != 0 && ! bj) {
1299         return OPJ_FALSE;
1300     }
1301     i = l;
1302
1303     while (i--) {
1304         OPJ_UINT32 j;
1305         OPJ_UINT32 rw;           /* width of the resolution level computed   */
1306         OPJ_UINT32 rh;           /* height of the resolution level computed  */
1307         OPJ_UINT32
1308         rw1;      /* width of the resolution level once lower than computed one                                       */
1309         OPJ_UINT32
1310         rh1;      /* height of the resolution level once lower than computed one                                      */
1311         OPJ_INT32 cas_col;  /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
1312         OPJ_INT32 cas_row;  /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
1313         OPJ_INT32 dn, sn;
1314
1315         rw  = (OPJ_UINT32)(l_cur_res->x1 - l_cur_res->x0);
1316         rh  = (OPJ_UINT32)(l_cur_res->y1 - l_cur_res->y0);
1317         rw1 = (OPJ_UINT32)(l_last_res->x1 - l_last_res->x0);
1318         rh1 = (OPJ_UINT32)(l_last_res->y1 - l_last_res->y0);
1319
1320         cas_row = l_cur_res->x0 & 1;
1321         cas_col = l_cur_res->y0 & 1;
1322
1323         sn = (OPJ_INT32)rh1;
1324         dn = (OPJ_INT32)(rh - rh1);
1325
1326         /* Perform vertical pass */
1327         if (num_threads <= 1 || rw <= 1) {
1328             for (j = 0; j < rw; ++j) {
1329                 OPJ_INT32* OPJ_RESTRICT aj = tiledp + j;
1330                 OPJ_UINT32 k;
1331                 for (k = 0; k < rh; ++k) {
1332                     bj[k] = aj[k * w];
1333                 }
1334
1335                 (*p_function)(bj, dn, sn, cas_col);
1336
1337                 opj_dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
1338             }
1339         }  else {
1340             OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
1341             OPJ_UINT32 step_j;
1342
1343             if (rw < num_jobs) {
1344                 num_jobs = rw;
1345             }
1346             step_j = (rw / num_jobs);
1347
1348             for (j = 0; j < num_jobs; j++) {
1349                 opj_dwt_encode_v_job_t* job;
1350
1351                 job = (opj_dwt_encode_v_job_t*) opj_malloc(sizeof(opj_dwt_encode_v_job_t));
1352                 if (!job) {
1353                     opj_thread_pool_wait_completion(tp, 0);
1354                     opj_aligned_free(bj);
1355                     return OPJ_FALSE;
1356                 }
1357                 job->v.mem = (OPJ_INT32*)opj_aligned_32_malloc(l_data_size);
1358                 if (!job->v.mem) {
1359                     opj_thread_pool_wait_completion(tp, 0);
1360                     opj_free(job);
1361                     opj_aligned_free(bj);
1362                     return OPJ_FALSE;
1363                 }
1364                 job->v.dn = dn;
1365                 job->v.sn = sn;
1366                 job->v.cas = cas_col;
1367                 job->rh = rh;
1368                 job->w = w;
1369                 job->tiledp = tiledp;
1370                 job->min_j = j * step_j;
1371                 job->max_j = (j + 1U) * step_j; /* this can overflow */
1372                 if (j == (num_jobs - 1U)) {  /* this will take care of the overflow */
1373                     job->max_j = rw;
1374                 }
1375                 job->p_function = p_function;
1376                 opj_thread_pool_submit_job(tp, opj_dwt_encode_v_func, job);
1377             }
1378             opj_thread_pool_wait_completion(tp, 0);
1379         }
1380
1381         sn = (OPJ_INT32)rw1;
1382         dn = (OPJ_INT32)(rw - rw1);
1383
1384         /* Perform horizontal pass */
1385         if (num_threads <= 1 || rh <= 1) {
1386             for (j = 0; j < rh; j++) {
1387                 OPJ_INT32* OPJ_RESTRICT aj = tiledp + j * w;
1388                 (*p_encode_and_deinterleave_h_one_row)(aj, bj, rw,
1389                                                        cas_row == 0 ? OPJ_TRUE : OPJ_FALSE);
1390             }
1391         }  else {
1392             OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
1393             OPJ_UINT32 step_j;
1394
1395             if (rh < num_jobs) {
1396                 num_jobs = rh;
1397             }
1398             step_j = (rh / num_jobs);
1399
1400             for (j = 0; j < num_jobs; j++) {
1401                 opj_dwt_encode_h_job_t* job;
1402
1403                 job = (opj_dwt_encode_h_job_t*) opj_malloc(sizeof(opj_dwt_encode_h_job_t));
1404                 if (!job) {
1405                     opj_thread_pool_wait_completion(tp, 0);
1406                     opj_aligned_free(bj);
1407                     return OPJ_FALSE;
1408                 }
1409                 job->h.mem = (OPJ_INT32*)opj_aligned_32_malloc(l_data_size);
1410                 if (!job->h.mem) {
1411                     opj_thread_pool_wait_completion(tp, 0);
1412                     opj_free(job);
1413                     opj_aligned_free(bj);
1414                     return OPJ_FALSE;
1415                 }
1416                 job->h.dn = dn;
1417                 job->h.sn = sn;
1418                 job->h.cas = cas_row;
1419                 job->rw = rw;
1420                 job->w = w;
1421                 job->tiledp = tiledp;
1422                 job->min_j = j * step_j;
1423                 job->max_j = (j + 1U) * step_j; /* this can overflow */
1424                 if (j == (num_jobs - 1U)) {  /* this will take care of the overflow */
1425                     job->max_j = rh;
1426                 }
1427                 job->p_function = p_encode_and_deinterleave_h_one_row;
1428                 opj_thread_pool_submit_job(tp, opj_dwt_encode_h_func, job);
1429             }
1430             opj_thread_pool_wait_completion(tp, 0);
1431         }
1432
1433         l_cur_res = l_last_res;
1434
1435         --l_last_res;
1436     }
1437
1438     opj_aligned_free(bj);
1439     return OPJ_TRUE;
1440 }
1441
1442 /* Forward 5-3 wavelet transform in 2-D. */
1443 /* </summary>                           */
1444 OPJ_BOOL opj_dwt_encode(opj_tcd_t *p_tcd,
1445                         opj_tcd_tilecomp_t * tilec)
1446 {
1447     return opj_dwt_encode_procedure(p_tcd->thread_pool, tilec,
1448                                     opj_dwt_encode_1,
1449                                     opj_dwt_encode_and_deinterleave_h_one_row);
1450 }
1451
1452 /* <summary>                            */
1453 /* Inverse 5-3 wavelet transform in 2-D. */
1454 /* </summary>                           */
1455 OPJ_BOOL opj_dwt_decode(opj_tcd_t *p_tcd, opj_tcd_tilecomp_t* tilec,
1456                         OPJ_UINT32 numres)
1457 {
1458     if (p_tcd->whole_tile_decoding) {
1459         return opj_dwt_decode_tile(p_tcd->thread_pool, tilec, numres);
1460     } else {
1461         return opj_dwt_decode_partial_tile(tilec, numres);
1462     }
1463 }
1464
1465 /* <summary>                */
1466 /* Get norm of 5-3 wavelet. */
1467 /* </summary>               */
1468 OPJ_FLOAT64 opj_dwt_getnorm(OPJ_UINT32 level, OPJ_UINT32 orient)
1469 {
1470     /* FIXME ! This is just a band-aid to avoid a buffer overflow */
1471     /* but the array should really be extended up to 33 resolution levels */
1472     /* See https://github.com/uclouvain/openjpeg/issues/493 */
1473     if (orient == 0 && level >= 10) {
1474         level = 9;
1475     } else if (orient > 0 && level >= 9) {
1476         level = 8;
1477     }
1478     return opj_dwt_norms[orient][level];
1479 }
1480
1481 /* <summary>                             */
1482 /* Forward 9-7 wavelet transform in 2-D. */
1483 /* </summary>                            */
1484 OPJ_BOOL opj_dwt_encode_real(opj_tcd_t *p_tcd,
1485                              opj_tcd_tilecomp_t * tilec)
1486 {
1487     return opj_dwt_encode_procedure(p_tcd->thread_pool, tilec,
1488                                     opj_dwt_encode_1_real,
1489                                     opj_dwt_encode_and_deinterleave_h_one_row_real);
1490 }
1491
1492 /* <summary>                */
1493 /* Get norm of 9-7 wavelet. */
1494 /* </summary>               */
1495 OPJ_FLOAT64 opj_dwt_getnorm_real(OPJ_UINT32 level, OPJ_UINT32 orient)
1496 {
1497     /* FIXME ! This is just a band-aid to avoid a buffer overflow */
1498     /* but the array should really be extended up to 33 resolution levels */
1499     /* See https://github.com/uclouvain/openjpeg/issues/493 */
1500     if (orient == 0 && level >= 10) {
1501         level = 9;
1502     } else if (orient > 0 && level >= 9) {
1503         level = 8;
1504     }
1505     return opj_dwt_norms_real[orient][level];
1506 }
1507
1508 void opj_dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, OPJ_UINT32 prec)
1509 {
1510     OPJ_UINT32 numbands, bandno;
1511     numbands = 3 * tccp->numresolutions - 2;
1512     for (bandno = 0; bandno < numbands; bandno++) {
1513         OPJ_FLOAT64 stepsize;
1514         OPJ_UINT32 resno, level, orient, gain;
1515
1516         resno = (bandno == 0) ? 0 : ((bandno - 1) / 3 + 1);
1517         orient = (bandno == 0) ? 0 : ((bandno - 1) % 3 + 1);
1518         level = tccp->numresolutions - 1 - resno;
1519         gain = (tccp->qmfbid == 0) ? 0 : ((orient == 0) ? 0 : (((orient == 1) ||
1520                                           (orient == 2)) ? 1 : 2));
1521         if (tccp->qntsty == J2K_CCP_QNTSTY_NOQNT) {
1522             stepsize = 1.0;
1523         } else {
1524             OPJ_FLOAT64 norm = opj_dwt_norms_real[orient][level];
1525             stepsize = (1 << (gain)) / norm;
1526         }
1527         opj_dwt_encode_stepsize((OPJ_INT32) floor(stepsize * 8192.0),
1528                                 (OPJ_INT32)(prec + gain), &tccp->stepsizes[bandno]);
1529     }
1530 }
1531
1532 /* <summary>                             */
1533 /* Determine maximum computed resolution level for inverse wavelet transform */
1534 /* </summary>                            */
1535 static OPJ_UINT32 opj_dwt_max_resolution(opj_tcd_resolution_t* OPJ_RESTRICT r,
1536         OPJ_UINT32 i)
1537 {
1538     OPJ_UINT32 mr   = 0;
1539     OPJ_UINT32 w;
1540     while (--i) {
1541         ++r;
1542         if (mr < (w = (OPJ_UINT32)(r->x1 - r->x0))) {
1543             mr = w ;
1544         }
1545         if (mr < (w = (OPJ_UINT32)(r->y1 - r->y0))) {
1546             mr = w ;
1547         }
1548     }
1549     return mr ;
1550 }
1551
1552 typedef struct {
1553     opj_dwt_t h;
1554     OPJ_UINT32 rw;
1555     OPJ_UINT32 w;
1556     OPJ_INT32 * OPJ_RESTRICT tiledp;
1557     OPJ_UINT32 min_j;
1558     OPJ_UINT32 max_j;
1559 } opj_dwt_decode_h_job_t;
1560
1561 static void opj_dwt_decode_h_func(void* user_data, opj_tls_t* tls)
1562 {
1563     OPJ_UINT32 j;
1564     opj_dwt_decode_h_job_t* job;
1565     (void)tls;
1566
1567     job = (opj_dwt_decode_h_job_t*)user_data;
1568     for (j = job->min_j; j < job->max_j; j++) {
1569         opj_idwt53_h(&job->h, &job->tiledp[j * job->w]);
1570     }
1571
1572     opj_aligned_free(job->h.mem);
1573     opj_free(job);
1574 }
1575
1576 typedef struct {
1577     opj_dwt_t v;
1578     OPJ_UINT32 rh;
1579     OPJ_UINT32 w;
1580     OPJ_INT32 * OPJ_RESTRICT tiledp;
1581     OPJ_UINT32 min_j;
1582     OPJ_UINT32 max_j;
1583 } opj_dwt_decode_v_job_t;
1584
1585 static void opj_dwt_decode_v_func(void* user_data, opj_tls_t* tls)
1586 {
1587     OPJ_UINT32 j;
1588     opj_dwt_decode_v_job_t* job;
1589     (void)tls;
1590
1591     job = (opj_dwt_decode_v_job_t*)user_data;
1592     for (j = job->min_j; j + PARALLEL_COLS_53 <= job->max_j;
1593             j += PARALLEL_COLS_53) {
1594         opj_idwt53_v(&job->v, &job->tiledp[j], (OPJ_SIZE_T)job->w,
1595                      PARALLEL_COLS_53);
1596     }
1597     if (j < job->max_j)
1598         opj_idwt53_v(&job->v, &job->tiledp[j], (OPJ_SIZE_T)job->w,
1599                      (OPJ_INT32)(job->max_j - j));
1600
1601     opj_aligned_free(job->v.mem);
1602     opj_free(job);
1603 }
1604
1605
1606 /* <summary>                            */
1607 /* Inverse wavelet transform in 2-D.    */
1608 /* </summary>                           */
1609 static OPJ_BOOL opj_dwt_decode_tile(opj_thread_pool_t* tp,
1610                                     opj_tcd_tilecomp_t* tilec, OPJ_UINT32 numres)
1611 {
1612     opj_dwt_t h;
1613     opj_dwt_t v;
1614
1615     opj_tcd_resolution_t* tr = tilec->resolutions;
1616
1617     OPJ_UINT32 rw = (OPJ_UINT32)(tr->x1 -
1618                                  tr->x0);  /* width of the resolution level computed */
1619     OPJ_UINT32 rh = (OPJ_UINT32)(tr->y1 -
1620                                  tr->y0);  /* height of the resolution level computed */
1621
1622     OPJ_UINT32 w = (OPJ_UINT32)(tilec->resolutions[tilec->minimum_num_resolutions -
1623                                                                1].x1 -
1624                                 tilec->resolutions[tilec->minimum_num_resolutions - 1].x0);
1625     OPJ_SIZE_T h_mem_size;
1626     int num_threads;
1627
1628     if (numres == 1U) {
1629         return OPJ_TRUE;
1630     }
1631     num_threads = opj_thread_pool_get_thread_count(tp);
1632     h_mem_size = opj_dwt_max_resolution(tr, numres);
1633     /* overflow check */
1634     if (h_mem_size > (SIZE_MAX / PARALLEL_COLS_53 / sizeof(OPJ_INT32))) {
1635         /* FIXME event manager error callback */
1636         return OPJ_FALSE;
1637     }
1638     /* We need PARALLEL_COLS_53 times the height of the array, */
1639     /* since for the vertical pass */
1640     /* we process PARALLEL_COLS_53 columns at a time */
1641     h_mem_size *= PARALLEL_COLS_53 * sizeof(OPJ_INT32);
1642     h.mem = (OPJ_INT32*)opj_aligned_32_malloc(h_mem_size);
1643     if (! h.mem) {
1644         /* FIXME event manager error callback */
1645         return OPJ_FALSE;
1646     }
1647
1648     v.mem = h.mem;
1649
1650     while (--numres) {
1651         OPJ_INT32 * OPJ_RESTRICT tiledp = tilec->data;
1652         OPJ_UINT32 j;
1653
1654         ++tr;
1655         h.sn = (OPJ_INT32)rw;
1656         v.sn = (OPJ_INT32)rh;
1657
1658         rw = (OPJ_UINT32)(tr->x1 - tr->x0);
1659         rh = (OPJ_UINT32)(tr->y1 - tr->y0);
1660
1661         h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
1662         h.cas = tr->x0 % 2;
1663
1664         if (num_threads <= 1 || rh <= 1) {
1665             for (j = 0; j < rh; ++j) {
1666                 opj_idwt53_h(&h, &tiledp[(OPJ_SIZE_T)j * w]);
1667             }
1668         } else {
1669             OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
1670             OPJ_UINT32 step_j;
1671
1672             if (rh < num_jobs) {
1673                 num_jobs = rh;
1674             }
1675             step_j = (rh / num_jobs);
1676
1677             for (j = 0; j < num_jobs; j++) {
1678                 opj_dwt_decode_h_job_t* job;
1679
1680                 job = (opj_dwt_decode_h_job_t*) opj_malloc(sizeof(opj_dwt_decode_h_job_t));
1681                 if (!job) {
1682                     /* It would be nice to fallback to single thread case, but */
1683                     /* unfortunately some jobs may be launched and have modified */
1684                     /* tiledp, so it is not practical to recover from that error */
1685                     /* FIXME event manager error callback */
1686                     opj_thread_pool_wait_completion(tp, 0);
1687                     opj_aligned_free(h.mem);
1688                     return OPJ_FALSE;
1689                 }
1690                 job->h = h;
1691                 job->rw = rw;
1692                 job->w = w;
1693                 job->tiledp = tiledp;
1694                 job->min_j = j * step_j;
1695                 job->max_j = (j + 1U) * step_j; /* this can overflow */
1696                 if (j == (num_jobs - 1U)) {  /* this will take care of the overflow */
1697                     job->max_j = rh;
1698                 }
1699                 job->h.mem = (OPJ_INT32*)opj_aligned_32_malloc(h_mem_size);
1700                 if (!job->h.mem) {
1701                     /* FIXME event manager error callback */
1702                     opj_thread_pool_wait_completion(tp, 0);
1703                     opj_free(job);
1704                     opj_aligned_free(h.mem);
1705                     return OPJ_FALSE;
1706                 }
1707                 opj_thread_pool_submit_job(tp, opj_dwt_decode_h_func, job);
1708             }
1709             opj_thread_pool_wait_completion(tp, 0);
1710         }
1711
1712         v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
1713         v.cas = tr->y0 % 2;
1714
1715         if (num_threads <= 1 || rw <= 1) {
1716             for (j = 0; j + PARALLEL_COLS_53 <= rw;
1717                     j += PARALLEL_COLS_53) {
1718                 opj_idwt53_v(&v, &tiledp[j], (OPJ_SIZE_T)w, PARALLEL_COLS_53);
1719             }
1720             if (j < rw) {
1721                 opj_idwt53_v(&v, &tiledp[j], (OPJ_SIZE_T)w, (OPJ_INT32)(rw - j));
1722             }
1723         } else {
1724             OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
1725             OPJ_UINT32 step_j;
1726
1727             if (rw < num_jobs) {
1728                 num_jobs = rw;
1729             }
1730             step_j = (rw / num_jobs);
1731
1732             for (j = 0; j < num_jobs; j++) {
1733                 opj_dwt_decode_v_job_t* job;
1734
1735                 job = (opj_dwt_decode_v_job_t*) opj_malloc(sizeof(opj_dwt_decode_v_job_t));
1736                 if (!job) {
1737                     /* It would be nice to fallback to single thread case, but */
1738                     /* unfortunately some jobs may be launched and have modified */
1739                     /* tiledp, so it is not practical to recover from that error */
1740                     /* FIXME event manager error callback */
1741                     opj_thread_pool_wait_completion(tp, 0);
1742                     opj_aligned_free(v.mem);
1743                     return OPJ_FALSE;
1744                 }
1745                 job->v = v;
1746                 job->rh = rh;
1747                 job->w = w;
1748                 job->tiledp = tiledp;
1749                 job->min_j = j * step_j;
1750                 job->max_j = (j + 1U) * step_j; /* this can overflow */
1751                 if (j == (num_jobs - 1U)) {  /* this will take care of the overflow */
1752                     job->max_j = rw;
1753                 }
1754                 job->v.mem = (OPJ_INT32*)opj_aligned_32_malloc(h_mem_size);
1755                 if (!job->v.mem) {
1756                     /* FIXME event manager error callback */
1757                     opj_thread_pool_wait_completion(tp, 0);
1758                     opj_free(job);
1759                     opj_aligned_free(v.mem);
1760                     return OPJ_FALSE;
1761                 }
1762                 opj_thread_pool_submit_job(tp, opj_dwt_decode_v_func, job);
1763             }
1764             opj_thread_pool_wait_completion(tp, 0);
1765         }
1766     }
1767     opj_aligned_free(h.mem);
1768     return OPJ_TRUE;
1769 }
1770
1771 static void opj_dwt_interleave_partial_h(OPJ_INT32 *dest,
1772         OPJ_INT32 cas,
1773         opj_sparse_array_int32_t* sa,
1774         OPJ_UINT32 sa_line,
1775         OPJ_UINT32 sn,
1776         OPJ_UINT32 win_l_x0,
1777         OPJ_UINT32 win_l_x1,
1778         OPJ_UINT32 win_h_x0,
1779         OPJ_UINT32 win_h_x1)
1780 {
1781     OPJ_BOOL ret;
1782     ret = opj_sparse_array_int32_read(sa,
1783                                       win_l_x0, sa_line,
1784                                       win_l_x1, sa_line + 1,
1785                                       dest + cas + 2 * win_l_x0,
1786                                       2, 0, OPJ_TRUE);
1787     assert(ret);
1788     ret = opj_sparse_array_int32_read(sa,
1789                                       sn + win_h_x0, sa_line,
1790                                       sn + win_h_x1, sa_line + 1,
1791                                       dest + 1 - cas + 2 * win_h_x0,
1792                                       2, 0, OPJ_TRUE);
1793     assert(ret);
1794     OPJ_UNUSED(ret);
1795 }
1796
1797
1798 static void opj_dwt_interleave_partial_v(OPJ_INT32 *dest,
1799         OPJ_INT32 cas,
1800         opj_sparse_array_int32_t* sa,
1801         OPJ_UINT32 sa_col,
1802         OPJ_UINT32 nb_cols,
1803         OPJ_UINT32 sn,
1804         OPJ_UINT32 win_l_y0,
1805         OPJ_UINT32 win_l_y1,
1806         OPJ_UINT32 win_h_y0,
1807         OPJ_UINT32 win_h_y1)
1808 {
1809     OPJ_BOOL ret;
1810     ret  = opj_sparse_array_int32_read(sa,
1811                                        sa_col, win_l_y0,
1812                                        sa_col + nb_cols, win_l_y1,
1813                                        dest + cas * 4 + 2 * 4 * win_l_y0,
1814                                        1, 2 * 4, OPJ_TRUE);
1815     assert(ret);
1816     ret = opj_sparse_array_int32_read(sa,
1817                                       sa_col, sn + win_h_y0,
1818                                       sa_col + nb_cols, sn + win_h_y1,
1819                                       dest + (1 - cas) * 4 + 2 * 4 * win_h_y0,
1820                                       1, 2 * 4, OPJ_TRUE);
1821     assert(ret);
1822     OPJ_UNUSED(ret);
1823 }
1824
1825 static void opj_dwt_decode_partial_1(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
1826                                      OPJ_INT32 cas,
1827                                      OPJ_INT32 win_l_x0,
1828                                      OPJ_INT32 win_l_x1,
1829                                      OPJ_INT32 win_h_x0,
1830                                      OPJ_INT32 win_h_x1)
1831 {
1832     OPJ_INT32 i;
1833
1834     if (!cas) {
1835         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
1836
1837             /* Naive version is :
1838             for (i = win_l_x0; i < i_max; i++) {
1839                 OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
1840             }
1841             for (i = win_h_x0; i < win_h_x1; i++) {
1842                 OPJ_D(i) += (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
1843             }
1844             but the compiler doesn't manage to unroll it to avoid bound
1845             checking in OPJ_S_ and OPJ_D_ macros
1846             */
1847
1848             i = win_l_x0;
1849             if (i < win_l_x1) {
1850                 OPJ_INT32 i_max;
1851
1852                 /* Left-most case */
1853                 OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
1854                 i ++;
1855
1856                 i_max = win_l_x1;
1857                 if (i_max > dn) {
1858                     i_max = dn;
1859                 }
1860                 for (; i < i_max; i++) {
1861                     /* No bound checking */
1862                     OPJ_S(i) -= (OPJ_D(i - 1) + OPJ_D(i) + 2) >> 2;
1863                 }
1864                 for (; i < win_l_x1; i++) {
1865                     /* Right-most case */
1866                     OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
1867                 }
1868             }
1869
1870             i = win_h_x0;
1871             if (i < win_h_x1) {
1872                 OPJ_INT32 i_max = win_h_x1;
1873                 if (i_max >= sn) {
1874                     i_max = sn - 1;
1875                 }
1876                 for (; i < i_max; i++) {
1877                     /* No bound checking */
1878                     OPJ_D(i) += (OPJ_S(i) + OPJ_S(i + 1)) >> 1;
1879                 }
1880                 for (; i < win_h_x1; i++) {
1881                     /* Right-most case */
1882                     OPJ_D(i) += (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
1883                 }
1884             }
1885         }
1886     } else {
1887         if (!sn  && dn == 1) {        /* NEW :  CASE ONE ELEMENT */
1888             OPJ_S(0) /= 2;
1889         } else {
1890             for (i = win_l_x0; i < win_l_x1; i++) {
1891                 OPJ_D(i) -= (OPJ_SS_(i) + OPJ_SS_(i + 1) + 2) >> 2;
1892             }
1893             for (i = win_h_x0; i < win_h_x1; i++) {
1894                 OPJ_S(i) += (OPJ_DD_(i) + OPJ_DD_(i - 1)) >> 1;
1895             }
1896         }
1897     }
1898 }
1899
1900 #define OPJ_S_off(i,off) a[(OPJ_UINT32)(i)*2*4+off]
1901 #define OPJ_D_off(i,off) a[(1+(OPJ_UINT32)(i)*2)*4+off]
1902 #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)))
1903 #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)))
1904 #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)))
1905 #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)))
1906
1907 static void opj_dwt_decode_partial_1_parallel(OPJ_INT32 *a,
1908         OPJ_UINT32 nb_cols,
1909         OPJ_INT32 dn, OPJ_INT32 sn,
1910         OPJ_INT32 cas,
1911         OPJ_INT32 win_l_x0,
1912         OPJ_INT32 win_l_x1,
1913         OPJ_INT32 win_h_x0,
1914         OPJ_INT32 win_h_x1)
1915 {
1916     OPJ_INT32 i;
1917     OPJ_UINT32 off;
1918
1919     (void)nb_cols;
1920
1921     if (!cas) {
1922         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
1923
1924             /* Naive version is :
1925             for (i = win_l_x0; i < i_max; i++) {
1926                 OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
1927             }
1928             for (i = win_h_x0; i < win_h_x1; i++) {
1929                 OPJ_D(i) += (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
1930             }
1931             but the compiler doesn't manage to unroll it to avoid bound
1932             checking in OPJ_S_ and OPJ_D_ macros
1933             */
1934
1935             i = win_l_x0;
1936             if (i < win_l_x1) {
1937                 OPJ_INT32 i_max;
1938
1939                 /* Left-most case */
1940                 for (off = 0; off < 4; off++) {
1941                     OPJ_S_off(i, off) -= (OPJ_D__off(i - 1, off) + OPJ_D__off(i, off) + 2) >> 2;
1942                 }
1943                 i ++;
1944
1945                 i_max = win_l_x1;
1946                 if (i_max > dn) {
1947                     i_max = dn;
1948                 }
1949
1950 #ifdef __SSE2__
1951                 if (i + 1 < i_max) {
1952                     const __m128i two = _mm_set1_epi32(2);
1953                     __m128i Dm1 = _mm_load_si128((__m128i * const)(a + 4 + (i - 1) * 8));
1954                     for (; i + 1 < i_max; i += 2) {
1955                         /* No bound checking */
1956                         __m128i S = _mm_load_si128((__m128i * const)(a + i * 8));
1957                         __m128i D = _mm_load_si128((__m128i * const)(a + 4 + i * 8));
1958                         __m128i S1 = _mm_load_si128((__m128i * const)(a + (i + 1) * 8));
1959                         __m128i D1 = _mm_load_si128((__m128i * const)(a + 4 + (i + 1) * 8));
1960                         S = _mm_sub_epi32(S,
1961                                           _mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(Dm1, D), two), 2));
1962                         S1 = _mm_sub_epi32(S1,
1963                                            _mm_srai_epi32(_mm_add_epi32(_mm_add_epi32(D, D1), two), 2));
1964                         _mm_store_si128((__m128i*)(a + i * 8), S);
1965                         _mm_store_si128((__m128i*)(a + (i + 1) * 8), S1);
1966                         Dm1 = D1;
1967                     }
1968                 }
1969 #endif
1970
1971                 for (; i < i_max; i++) {
1972                     /* No bound checking */
1973                     for (off = 0; off < 4; off++) {
1974                         OPJ_S_off(i, off) -= (OPJ_D_off(i - 1, off) + OPJ_D_off(i, off) + 2) >> 2;
1975                     }
1976                 }
1977                 for (; i < win_l_x1; i++) {
1978                     /* Right-most case */
1979                     for (off = 0; off < 4; off++) {
1980                         OPJ_S_off(i, off) -= (OPJ_D__off(i - 1, off) + OPJ_D__off(i, off) + 2) >> 2;
1981                     }
1982                 }
1983             }
1984
1985             i = win_h_x0;
1986             if (i < win_h_x1) {
1987                 OPJ_INT32 i_max = win_h_x1;
1988                 if (i_max >= sn) {
1989                     i_max = sn - 1;
1990                 }
1991
1992 #ifdef __SSE2__
1993                 if (i + 1 < i_max) {
1994                     __m128i S =  _mm_load_si128((__m128i * const)(a + i * 8));
1995                     for (; i + 1 < i_max; i += 2) {
1996                         /* No bound checking */
1997                         __m128i D = _mm_load_si128((__m128i * const)(a + 4 + i * 8));
1998                         __m128i S1 = _mm_load_si128((__m128i * const)(a + (i + 1) * 8));
1999                         __m128i D1 = _mm_load_si128((__m128i * const)(a + 4 + (i + 1) * 8));
2000                         __m128i S2 = _mm_load_si128((__m128i * const)(a + (i + 2) * 8));
2001                         D = _mm_add_epi32(D, _mm_srai_epi32(_mm_add_epi32(S, S1), 1));
2002                         D1 = _mm_add_epi32(D1, _mm_srai_epi32(_mm_add_epi32(S1, S2), 1));
2003                         _mm_store_si128((__m128i*)(a + 4 + i * 8), D);
2004                         _mm_store_si128((__m128i*)(a + 4 + (i + 1) * 8), D1);
2005                         S = S2;
2006                     }
2007                 }
2008 #endif
2009
2010                 for (; i < i_max; i++) {
2011                     /* No bound checking */
2012                     for (off = 0; off < 4; off++) {
2013                         OPJ_D_off(i, off) += (OPJ_S_off(i, off) + OPJ_S_off(i + 1, off)) >> 1;
2014                     }
2015                 }
2016                 for (; i < win_h_x1; i++) {
2017                     /* Right-most case */
2018                     for (off = 0; off < 4; off++) {
2019                         OPJ_D_off(i, off) += (OPJ_S__off(i, off) + OPJ_S__off(i + 1, off)) >> 1;
2020                     }
2021                 }
2022             }
2023         }
2024     } else {
2025         if (!sn  && dn == 1) {        /* NEW :  CASE ONE ELEMENT */
2026             for (off = 0; off < 4; off++) {
2027                 OPJ_S_off(0, off) /= 2;
2028             }
2029         } else {
2030             for (i = win_l_x0; i < win_l_x1; i++) {
2031                 for (off = 0; off < 4; off++) {
2032                     OPJ_D_off(i, off) -= (OPJ_SS__off(i, off) + OPJ_SS__off(i + 1, off) + 2) >> 2;
2033                 }
2034             }
2035             for (i = win_h_x0; i < win_h_x1; i++) {
2036                 for (off = 0; off < 4; off++) {
2037                     OPJ_S_off(i, off) += (OPJ_DD__off(i, off) + OPJ_DD__off(i - 1, off)) >> 1;
2038                 }
2039             }
2040         }
2041     }
2042 }
2043
2044 static void opj_dwt_get_band_coordinates(opj_tcd_tilecomp_t* tilec,
2045         OPJ_UINT32 resno,
2046         OPJ_UINT32 bandno,
2047         OPJ_UINT32 tcx0,
2048         OPJ_UINT32 tcy0,
2049         OPJ_UINT32 tcx1,
2050         OPJ_UINT32 tcy1,
2051         OPJ_UINT32* tbx0,
2052         OPJ_UINT32* tby0,
2053         OPJ_UINT32* tbx1,
2054         OPJ_UINT32* tby1)
2055 {
2056     /* Compute number of decomposition for this band. See table F-1 */
2057     OPJ_UINT32 nb = (resno == 0) ?
2058                     tilec->numresolutions - 1 :
2059                     tilec->numresolutions - resno;
2060     /* Map above tile-based coordinates to sub-band-based coordinates per */
2061     /* equation B-15 of the standard */
2062     OPJ_UINT32 x0b = bandno & 1;
2063     OPJ_UINT32 y0b = bandno >> 1;
2064     if (tbx0) {
2065         *tbx0 = (nb == 0) ? tcx0 :
2066                 (tcx0 <= (1U << (nb - 1)) * x0b) ? 0 :
2067                 opj_uint_ceildivpow2(tcx0 - (1U << (nb - 1)) * x0b, nb);
2068     }
2069     if (tby0) {
2070         *tby0 = (nb == 0) ? tcy0 :
2071                 (tcy0 <= (1U << (nb - 1)) * y0b) ? 0 :
2072                 opj_uint_ceildivpow2(tcy0 - (1U << (nb - 1)) * y0b, nb);
2073     }
2074     if (tbx1) {
2075         *tbx1 = (nb == 0) ? tcx1 :
2076                 (tcx1 <= (1U << (nb - 1)) * x0b) ? 0 :
2077                 opj_uint_ceildivpow2(tcx1 - (1U << (nb - 1)) * x0b, nb);
2078     }
2079     if (tby1) {
2080         *tby1 = (nb == 0) ? tcy1 :
2081                 (tcy1 <= (1U << (nb - 1)) * y0b) ? 0 :
2082                 opj_uint_ceildivpow2(tcy1 - (1U << (nb - 1)) * y0b, nb);
2083     }
2084 }
2085
2086 static void opj_dwt_segment_grow(OPJ_UINT32 filter_width,
2087                                  OPJ_UINT32 max_size,
2088                                  OPJ_UINT32* start,
2089                                  OPJ_UINT32* end)
2090 {
2091     *start = opj_uint_subs(*start, filter_width);
2092     *end = opj_uint_adds(*end, filter_width);
2093     *end = opj_uint_min(*end, max_size);
2094 }
2095
2096
2097 static opj_sparse_array_int32_t* opj_dwt_init_sparse_array(
2098     opj_tcd_tilecomp_t* tilec,
2099     OPJ_UINT32 numres)
2100 {
2101     opj_tcd_resolution_t* tr_max = &(tilec->resolutions[numres - 1]);
2102     OPJ_UINT32 w = (OPJ_UINT32)(tr_max->x1 - tr_max->x0);
2103     OPJ_UINT32 h = (OPJ_UINT32)(tr_max->y1 - tr_max->y0);
2104     OPJ_UINT32 resno, bandno, precno, cblkno;
2105     opj_sparse_array_int32_t* sa = opj_sparse_array_int32_create(
2106                                        w, h, opj_uint_min(w, 64), opj_uint_min(h, 64));
2107     if (sa == NULL) {
2108         return NULL;
2109     }
2110
2111     for (resno = 0; resno < numres; ++resno) {
2112         opj_tcd_resolution_t* res = &tilec->resolutions[resno];
2113
2114         for (bandno = 0; bandno < res->numbands; ++bandno) {
2115             opj_tcd_band_t* band = &res->bands[bandno];
2116
2117             for (precno = 0; precno < res->pw * res->ph; ++precno) {
2118                 opj_tcd_precinct_t* precinct = &band->precincts[precno];
2119                 for (cblkno = 0; cblkno < precinct->cw * precinct->ch; ++cblkno) {
2120                     opj_tcd_cblk_dec_t* cblk = &precinct->cblks.dec[cblkno];
2121                     if (cblk->decoded_data != NULL) {
2122                         OPJ_UINT32 x = (OPJ_UINT32)(cblk->x0 - band->x0);
2123                         OPJ_UINT32 y = (OPJ_UINT32)(cblk->y0 - band->y0);
2124                         OPJ_UINT32 cblk_w = (OPJ_UINT32)(cblk->x1 - cblk->x0);
2125                         OPJ_UINT32 cblk_h = (OPJ_UINT32)(cblk->y1 - cblk->y0);
2126
2127                         if (band->bandno & 1) {
2128                             opj_tcd_resolution_t* pres = &tilec->resolutions[resno - 1];
2129                             x += (OPJ_UINT32)(pres->x1 - pres->x0);
2130                         }
2131                         if (band->bandno & 2) {
2132                             opj_tcd_resolution_t* pres = &tilec->resolutions[resno - 1];
2133                             y += (OPJ_UINT32)(pres->y1 - pres->y0);
2134                         }
2135
2136                         if (!opj_sparse_array_int32_write(sa, x, y,
2137                                                           x + cblk_w, y + cblk_h,
2138                                                           cblk->decoded_data,
2139                                                           1, cblk_w, OPJ_TRUE)) {
2140                             opj_sparse_array_int32_free(sa);
2141                             return NULL;
2142                         }
2143                     }
2144                 }
2145             }
2146         }
2147     }
2148
2149     return sa;
2150 }
2151
2152
2153 static OPJ_BOOL opj_dwt_decode_partial_tile(
2154     opj_tcd_tilecomp_t* tilec,
2155     OPJ_UINT32 numres)
2156 {
2157     opj_sparse_array_int32_t* sa;
2158     opj_dwt_t h;
2159     opj_dwt_t v;
2160     OPJ_UINT32 resno;
2161     /* This value matches the maximum left/right extension given in tables */
2162     /* F.2 and F.3 of the standard. */
2163     const OPJ_UINT32 filter_width = 2U;
2164
2165     opj_tcd_resolution_t* tr = tilec->resolutions;
2166     opj_tcd_resolution_t* tr_max = &(tilec->resolutions[numres - 1]);
2167
2168     OPJ_UINT32 rw = (OPJ_UINT32)(tr->x1 -
2169                                  tr->x0);  /* width of the resolution level computed */
2170     OPJ_UINT32 rh = (OPJ_UINT32)(tr->y1 -
2171                                  tr->y0);  /* height of the resolution level computed */
2172
2173     OPJ_SIZE_T h_mem_size;
2174
2175     /* Compute the intersection of the area of interest, expressed in tile coordinates */
2176     /* with the tile coordinates */
2177     OPJ_UINT32 win_tcx0 = tilec->win_x0;
2178     OPJ_UINT32 win_tcy0 = tilec->win_y0;
2179     OPJ_UINT32 win_tcx1 = tilec->win_x1;
2180     OPJ_UINT32 win_tcy1 = tilec->win_y1;
2181
2182     if (tr_max->x0 == tr_max->x1 || tr_max->y0 == tr_max->y1) {
2183         return OPJ_TRUE;
2184     }
2185
2186     sa = opj_dwt_init_sparse_array(tilec, numres);
2187     if (sa == NULL) {
2188         return OPJ_FALSE;
2189     }
2190
2191     if (numres == 1U) {
2192         OPJ_BOOL ret = opj_sparse_array_int32_read(sa,
2193                        tr_max->win_x0 - (OPJ_UINT32)tr_max->x0,
2194                        tr_max->win_y0 - (OPJ_UINT32)tr_max->y0,
2195                        tr_max->win_x1 - (OPJ_UINT32)tr_max->x0,
2196                        tr_max->win_y1 - (OPJ_UINT32)tr_max->y0,
2197                        tilec->data_win,
2198                        1, tr_max->win_x1 - tr_max->win_x0,
2199                        OPJ_TRUE);
2200         assert(ret);
2201         OPJ_UNUSED(ret);
2202         opj_sparse_array_int32_free(sa);
2203         return OPJ_TRUE;
2204     }
2205     h_mem_size = opj_dwt_max_resolution(tr, numres);
2206     /* overflow check */
2207     /* in vertical pass, we process 4 columns at a time */
2208     if (h_mem_size > (SIZE_MAX / (4 * sizeof(OPJ_INT32)))) {
2209         /* FIXME event manager error callback */
2210         opj_sparse_array_int32_free(sa);
2211         return OPJ_FALSE;
2212     }
2213
2214     h_mem_size *= 4 * sizeof(OPJ_INT32);
2215     h.mem = (OPJ_INT32*)opj_aligned_32_malloc(h_mem_size);
2216     if (! h.mem) {
2217         /* FIXME event manager error callback */
2218         opj_sparse_array_int32_free(sa);
2219         return OPJ_FALSE;
2220     }
2221
2222     v.mem = h.mem;
2223
2224     for (resno = 1; resno < numres; resno ++) {
2225         OPJ_UINT32 i, j;
2226         /* Window of interest subband-based coordinates */
2227         OPJ_UINT32 win_ll_x0, win_ll_y0, win_ll_x1, win_ll_y1;
2228         OPJ_UINT32 win_hl_x0, win_hl_x1;
2229         OPJ_UINT32 win_lh_y0, win_lh_y1;
2230         /* Window of interest tile-resolution-based coordinates */
2231         OPJ_UINT32 win_tr_x0, win_tr_x1, win_tr_y0, win_tr_y1;
2232         /* Tile-resolution subband-based coordinates */
2233         OPJ_UINT32 tr_ll_x0, tr_ll_y0, tr_hl_x0, tr_lh_y0;
2234
2235         ++tr;
2236
2237         h.sn = (OPJ_INT32)rw;
2238         v.sn = (OPJ_INT32)rh;
2239
2240         rw = (OPJ_UINT32)(tr->x1 - tr->x0);
2241         rh = (OPJ_UINT32)(tr->y1 - tr->y0);
2242
2243         h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
2244         h.cas = tr->x0 % 2;
2245
2246         v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
2247         v.cas = tr->y0 % 2;
2248
2249         /* Get the subband coordinates for the window of interest */
2250         /* LL band */
2251         opj_dwt_get_band_coordinates(tilec, resno, 0,
2252                                      win_tcx0, win_tcy0, win_tcx1, win_tcy1,
2253                                      &win_ll_x0, &win_ll_y0,
2254                                      &win_ll_x1, &win_ll_y1);
2255
2256         /* HL band */
2257         opj_dwt_get_band_coordinates(tilec, resno, 1,
2258                                      win_tcx0, win_tcy0, win_tcx1, win_tcy1,
2259                                      &win_hl_x0, NULL, &win_hl_x1, NULL);
2260
2261         /* LH band */
2262         opj_dwt_get_band_coordinates(tilec, resno, 2,
2263                                      win_tcx0, win_tcy0, win_tcx1, win_tcy1,
2264                                      NULL, &win_lh_y0, NULL, &win_lh_y1);
2265
2266         /* Beware: band index for non-LL0 resolution are 0=HL, 1=LH and 2=HH */
2267         tr_ll_x0 = (OPJ_UINT32)tr->bands[1].x0;
2268         tr_ll_y0 = (OPJ_UINT32)tr->bands[0].y0;
2269         tr_hl_x0 = (OPJ_UINT32)tr->bands[0].x0;
2270         tr_lh_y0 = (OPJ_UINT32)tr->bands[1].y0;
2271
2272         /* Subtract the origin of the bands for this tile, to the subwindow */
2273         /* of interest band coordinates, so as to get them relative to the */
2274         /* tile */
2275         win_ll_x0 = opj_uint_subs(win_ll_x0, tr_ll_x0);
2276         win_ll_y0 = opj_uint_subs(win_ll_y0, tr_ll_y0);
2277         win_ll_x1 = opj_uint_subs(win_ll_x1, tr_ll_x0);
2278         win_ll_y1 = opj_uint_subs(win_ll_y1, tr_ll_y0);
2279         win_hl_x0 = opj_uint_subs(win_hl_x0, tr_hl_x0);
2280         win_hl_x1 = opj_uint_subs(win_hl_x1, tr_hl_x0);
2281         win_lh_y0 = opj_uint_subs(win_lh_y0, tr_lh_y0);
2282         win_lh_y1 = opj_uint_subs(win_lh_y1, tr_lh_y0);
2283
2284         opj_dwt_segment_grow(filter_width, (OPJ_UINT32)h.sn, &win_ll_x0, &win_ll_x1);
2285         opj_dwt_segment_grow(filter_width, (OPJ_UINT32)h.dn, &win_hl_x0, &win_hl_x1);
2286
2287         opj_dwt_segment_grow(filter_width, (OPJ_UINT32)v.sn, &win_ll_y0, &win_ll_y1);
2288         opj_dwt_segment_grow(filter_width, (OPJ_UINT32)v.dn, &win_lh_y0, &win_lh_y1);
2289
2290         /* Compute the tile-resolution-based coordinates for the window of interest */
2291         if (h.cas == 0) {
2292             win_tr_x0 = opj_uint_min(2 * win_ll_x0, 2 * win_hl_x0 + 1);
2293             win_tr_x1 = opj_uint_min(opj_uint_max(2 * win_ll_x1, 2 * win_hl_x1 + 1), rw);
2294         } else {
2295             win_tr_x0 = opj_uint_min(2 * win_hl_x0, 2 * win_ll_x0 + 1);
2296             win_tr_x1 = opj_uint_min(opj_uint_max(2 * win_hl_x1, 2 * win_ll_x1 + 1), rw);
2297         }
2298
2299         if (v.cas == 0) {
2300             win_tr_y0 = opj_uint_min(2 * win_ll_y0, 2 * win_lh_y0 + 1);
2301             win_tr_y1 = opj_uint_min(opj_uint_max(2 * win_ll_y1, 2 * win_lh_y1 + 1), rh);
2302         } else {
2303             win_tr_y0 = opj_uint_min(2 * win_lh_y0, 2 * win_ll_y0 + 1);
2304             win_tr_y1 = opj_uint_min(opj_uint_max(2 * win_lh_y1, 2 * win_ll_y1 + 1), rh);
2305         }
2306
2307         for (j = 0; j < rh; ++j) {
2308             if ((j >= win_ll_y0 && j < win_ll_y1) ||
2309                     (j >= win_lh_y0 + (OPJ_UINT32)v.sn && j < win_lh_y1 + (OPJ_UINT32)v.sn)) {
2310
2311                 /* Avoids dwt.c:1584:44 (in opj_dwt_decode_partial_1): runtime error: */
2312                 /* signed integer overflow: -1094795586 + -1094795586 cannot be represented in type 'int' */
2313                 /* on opj_decompress -i  ../../openjpeg/MAPA.jp2 -o out.tif -d 0,0,256,256 */
2314                 /* This is less extreme than memsetting the whole buffer to 0 */
2315                 /* although we could potentially do better with better handling of edge conditions */
2316                 if (win_tr_x1 >= 1 && win_tr_x1 < rw) {
2317                     h.mem[win_tr_x1 - 1] = 0;
2318                 }
2319                 if (win_tr_x1 < rw) {
2320                     h.mem[win_tr_x1] = 0;
2321                 }
2322
2323                 opj_dwt_interleave_partial_h(h.mem,
2324                                              h.cas,
2325                                              sa,
2326                                              j,
2327                                              (OPJ_UINT32)h.sn,
2328                                              win_ll_x0,
2329                                              win_ll_x1,
2330                                              win_hl_x0,
2331                                              win_hl_x1);
2332                 opj_dwt_decode_partial_1(h.mem, h.dn, h.sn, h.cas,
2333                                          (OPJ_INT32)win_ll_x0,
2334                                          (OPJ_INT32)win_ll_x1,
2335                                          (OPJ_INT32)win_hl_x0,
2336                                          (OPJ_INT32)win_hl_x1);
2337                 if (!opj_sparse_array_int32_write(sa,
2338                                                   win_tr_x0, j,
2339                                                   win_tr_x1, j + 1,
2340                                                   h.mem + win_tr_x0,
2341                                                   1, 0, OPJ_TRUE)) {
2342                     /* FIXME event manager error callback */
2343                     opj_sparse_array_int32_free(sa);
2344                     opj_aligned_free(h.mem);
2345                     return OPJ_FALSE;
2346                 }
2347             }
2348         }
2349
2350         for (i = win_tr_x0; i < win_tr_x1;) {
2351             OPJ_UINT32 nb_cols = opj_uint_min(4U, win_tr_x1 - i);
2352             opj_dwt_interleave_partial_v(v.mem,
2353                                          v.cas,
2354                                          sa,
2355                                          i,
2356                                          nb_cols,
2357                                          (OPJ_UINT32)v.sn,
2358                                          win_ll_y0,
2359                                          win_ll_y1,
2360                                          win_lh_y0,
2361                                          win_lh_y1);
2362             opj_dwt_decode_partial_1_parallel(v.mem, nb_cols, v.dn, v.sn, v.cas,
2363                                               (OPJ_INT32)win_ll_y0,
2364                                               (OPJ_INT32)win_ll_y1,
2365                                               (OPJ_INT32)win_lh_y0,
2366                                               (OPJ_INT32)win_lh_y1);
2367             if (!opj_sparse_array_int32_write(sa,
2368                                               i, win_tr_y0,
2369                                               i + nb_cols, win_tr_y1,
2370                                               v.mem + 4 * win_tr_y0,
2371                                               1, 4, OPJ_TRUE)) {
2372                 /* FIXME event manager error callback */
2373                 opj_sparse_array_int32_free(sa);
2374                 opj_aligned_free(h.mem);
2375                 return OPJ_FALSE;
2376             }
2377
2378             i += nb_cols;
2379         }
2380     }
2381     opj_aligned_free(h.mem);
2382
2383     {
2384         OPJ_BOOL ret = opj_sparse_array_int32_read(sa,
2385                        tr_max->win_x0 - (OPJ_UINT32)tr_max->x0,
2386                        tr_max->win_y0 - (OPJ_UINT32)tr_max->y0,
2387                        tr_max->win_x1 - (OPJ_UINT32)tr_max->x0,
2388                        tr_max->win_y1 - (OPJ_UINT32)tr_max->y0,
2389                        tilec->data_win,
2390                        1, tr_max->win_x1 - tr_max->win_x0,
2391                        OPJ_TRUE);
2392         assert(ret);
2393         OPJ_UNUSED(ret);
2394     }
2395     opj_sparse_array_int32_free(sa);
2396     return OPJ_TRUE;
2397 }
2398
2399 static void opj_v8dwt_interleave_h(opj_v8dwt_t* OPJ_RESTRICT dwt,
2400                                    OPJ_FLOAT32* OPJ_RESTRICT a,
2401                                    OPJ_UINT32 width,
2402                                    OPJ_UINT32 remaining_height)
2403 {
2404     OPJ_FLOAT32* OPJ_RESTRICT bi = (OPJ_FLOAT32*)(dwt->wavelet + dwt->cas);
2405     OPJ_UINT32 i, k;
2406     OPJ_UINT32 x0 = dwt->win_l_x0;
2407     OPJ_UINT32 x1 = dwt->win_l_x1;
2408
2409     for (k = 0; k < 2; ++k) {
2410         if (remaining_height >= NB_ELTS_V8 && ((OPJ_SIZE_T) a & 0x0f) == 0 &&
2411                 ((OPJ_SIZE_T) bi & 0x0f) == 0) {
2412             /* Fast code path */
2413             for (i = x0; i < x1; ++i) {
2414                 OPJ_UINT32 j = i;
2415                 OPJ_FLOAT32* OPJ_RESTRICT dst = bi + i * 2 * NB_ELTS_V8;
2416                 dst[0] = a[j];
2417                 j += width;
2418                 dst[1] = a[j];
2419                 j += width;
2420                 dst[2] = a[j];
2421                 j += width;
2422                 dst[3] = a[j];
2423                 j += width;
2424                 dst[4] = a[j];
2425                 j += width;
2426                 dst[5] = a[j];
2427                 j += width;
2428                 dst[6] = a[j];
2429                 j += width;
2430                 dst[7] = a[j];
2431             }
2432         } else {
2433             /* Slow code path */
2434             for (i = x0; i < x1; ++i) {
2435                 OPJ_UINT32 j = i;
2436                 OPJ_FLOAT32* OPJ_RESTRICT dst = bi + i * 2 * NB_ELTS_V8;
2437                 dst[0] = a[j];
2438                 j += width;
2439                 if (remaining_height == 1) {
2440                     continue;
2441                 }
2442                 dst[1] = a[j];
2443                 j += width;
2444                 if (remaining_height == 2) {
2445                     continue;
2446                 }
2447                 dst[2] = a[j];
2448                 j += width;
2449                 if (remaining_height == 3) {
2450                     continue;
2451                 }
2452                 dst[3] = a[j];
2453                 j += width;
2454                 if (remaining_height == 4) {
2455                     continue;
2456                 }
2457                 dst[4] = a[j];
2458                 j += width;
2459                 if (remaining_height == 5) {
2460                     continue;
2461                 }
2462                 dst[5] = a[j];
2463                 j += width;
2464                 if (remaining_height == 6) {
2465                     continue;
2466                 }
2467                 dst[6] = a[j];
2468                 j += width;
2469                 if (remaining_height == 7) {
2470                     continue;
2471                 }
2472                 dst[7] = a[j];
2473             }
2474         }
2475
2476         bi = (OPJ_FLOAT32*)(dwt->wavelet + 1 - dwt->cas);
2477         a += dwt->sn;
2478         x0 = dwt->win_h_x0;
2479         x1 = dwt->win_h_x1;
2480     }
2481 }
2482
2483 static void opj_v8dwt_interleave_partial_h(opj_v8dwt_t* dwt,
2484         opj_sparse_array_int32_t* sa,
2485         OPJ_UINT32 sa_line,
2486         OPJ_UINT32 remaining_height)
2487 {
2488     OPJ_UINT32 i;
2489     for (i = 0; i < remaining_height; i++) {
2490         OPJ_BOOL ret;
2491         ret = opj_sparse_array_int32_read(sa,
2492                                           dwt->win_l_x0, sa_line + i,
2493                                           dwt->win_l_x1, sa_line + i + 1,
2494                                           /* Nasty cast from float* to int32* */
2495                                           (OPJ_INT32*)(dwt->wavelet + dwt->cas + 2 * dwt->win_l_x0) + i,
2496                                           2 * NB_ELTS_V8, 0, OPJ_TRUE);
2497         assert(ret);
2498         ret = opj_sparse_array_int32_read(sa,
2499                                           (OPJ_UINT32)dwt->sn + dwt->win_h_x0, sa_line + i,
2500                                           (OPJ_UINT32)dwt->sn + dwt->win_h_x1, sa_line + i + 1,
2501                                           /* Nasty cast from float* to int32* */
2502                                           (OPJ_INT32*)(dwt->wavelet + 1 - dwt->cas + 2 * dwt->win_h_x0) + i,
2503                                           2 * NB_ELTS_V8, 0, OPJ_TRUE);
2504         assert(ret);
2505         OPJ_UNUSED(ret);
2506     }
2507 }
2508
2509 static INLINE void opj_v8dwt_interleave_v(opj_v8dwt_t* OPJ_RESTRICT dwt,
2510         OPJ_FLOAT32* OPJ_RESTRICT a,
2511         OPJ_UINT32 width,
2512         OPJ_UINT32 nb_elts_read)
2513 {
2514     opj_v8_t* OPJ_RESTRICT bi = dwt->wavelet + dwt->cas;
2515     OPJ_UINT32 i;
2516
2517     for (i = dwt->win_l_x0; i < dwt->win_l_x1; ++i) {
2518         memcpy(&bi[i * 2], &a[i * (OPJ_SIZE_T)width],
2519                (OPJ_SIZE_T)nb_elts_read * sizeof(OPJ_FLOAT32));
2520     }
2521
2522     a += (OPJ_UINT32)dwt->sn * (OPJ_SIZE_T)width;
2523     bi = dwt->wavelet + 1 - dwt->cas;
2524
2525     for (i = dwt->win_h_x0; i < dwt->win_h_x1; ++i) {
2526         memcpy(&bi[i * 2], &a[i * (OPJ_SIZE_T)width],
2527                (OPJ_SIZE_T)nb_elts_read * sizeof(OPJ_FLOAT32));
2528     }
2529 }
2530
2531 static void opj_v8dwt_interleave_partial_v(opj_v8dwt_t* OPJ_RESTRICT dwt,
2532         opj_sparse_array_int32_t* sa,
2533         OPJ_UINT32 sa_col,
2534         OPJ_UINT32 nb_elts_read)
2535 {
2536     OPJ_BOOL ret;
2537     ret = opj_sparse_array_int32_read(sa,
2538                                       sa_col, dwt->win_l_x0,
2539                                       sa_col + nb_elts_read, dwt->win_l_x1,
2540                                       (OPJ_INT32*)(dwt->wavelet + dwt->cas + 2 * dwt->win_l_x0),
2541                                       1, 2 * NB_ELTS_V8, OPJ_TRUE);
2542     assert(ret);
2543     ret = opj_sparse_array_int32_read(sa,
2544                                       sa_col, (OPJ_UINT32)dwt->sn + dwt->win_h_x0,
2545                                       sa_col + nb_elts_read, (OPJ_UINT32)dwt->sn + dwt->win_h_x1,
2546                                       (OPJ_INT32*)(dwt->wavelet + 1 - dwt->cas + 2 * dwt->win_h_x0),
2547                                       1, 2 * NB_ELTS_V8, OPJ_TRUE);
2548     assert(ret);
2549     OPJ_UNUSED(ret);
2550 }
2551
2552 #ifdef __SSE__
2553
2554 static void opj_v8dwt_decode_step1_sse(opj_v8_t* w,
2555                                        OPJ_UINT32 start,
2556                                        OPJ_UINT32 end,
2557                                        const __m128 c)
2558 {
2559     __m128* OPJ_RESTRICT vw = (__m128*) w;
2560     OPJ_UINT32 i = start;
2561     /* To be adapted if NB_ELTS_V8 changes */
2562     vw += 4 * start;
2563     /* Note: attempt at loop unrolling x2 doesn't help */
2564     for (; i < end; ++i, vw += 4) {
2565         vw[0] = _mm_mul_ps(vw[0], c);
2566         vw[1] = _mm_mul_ps(vw[1], c);
2567     }
2568 }
2569
2570 static void opj_v8dwt_decode_step2_sse(opj_v8_t* l, opj_v8_t* w,
2571                                        OPJ_UINT32 start,
2572                                        OPJ_UINT32 end,
2573                                        OPJ_UINT32 m,
2574                                        __m128 c)
2575 {
2576     __m128* OPJ_RESTRICT vl = (__m128*) l;
2577     __m128* OPJ_RESTRICT vw = (__m128*) w;
2578     /* To be adapted if NB_ELTS_V8 changes */
2579     OPJ_UINT32 i;
2580     OPJ_UINT32 imax = opj_uint_min(end, m);
2581     if (start == 0) {
2582         if (imax >= 1) {
2583             vw[-2] = _mm_add_ps(vw[-2], _mm_mul_ps(_mm_add_ps(vl[0], vw[0]), c));
2584             vw[-1] = _mm_add_ps(vw[-1], _mm_mul_ps(_mm_add_ps(vl[1], vw[1]), c));
2585             vw += 4;
2586             start = 1;
2587         }
2588     } else {
2589         vw += start * 4;
2590     }
2591
2592     i = start;
2593     /* Note: attempt at loop unrolling x2 doesn't help */
2594     for (; i < imax; ++i) {
2595         vw[-2] = _mm_add_ps(vw[-2], _mm_mul_ps(_mm_add_ps(vw[-4], vw[0]), c));
2596         vw[-1] = _mm_add_ps(vw[-1], _mm_mul_ps(_mm_add_ps(vw[-3], vw[1]), c));
2597         vw += 4;
2598     }
2599     if (m < end) {
2600         assert(m + 1 == end);
2601         c = _mm_add_ps(c, c);
2602         vw[-2] = _mm_add_ps(vw[-2], _mm_mul_ps(c, vw[-4]));
2603         vw[-1] = _mm_add_ps(vw[-1], _mm_mul_ps(c, vw[-3]));
2604     }
2605 }
2606
2607 #else
2608
2609 static void opj_v8dwt_decode_step1(opj_v8_t* w,
2610                                    OPJ_UINT32 start,
2611                                    OPJ_UINT32 end,
2612                                    const OPJ_FLOAT32 c)
2613 {
2614     OPJ_FLOAT32* OPJ_RESTRICT fw = (OPJ_FLOAT32*) w;
2615     OPJ_UINT32 i;
2616     /* To be adapted if NB_ELTS_V8 changes */
2617     for (i = start; i < end; ++i) {
2618         fw[i * 2 * 8    ] = fw[i * 2 * 8    ] * c;
2619         fw[i * 2 * 8 + 1] = fw[i * 2 * 8 + 1] * c;
2620         fw[i * 2 * 8 + 2] = fw[i * 2 * 8 + 2] * c;
2621         fw[i * 2 * 8 + 3] = fw[i * 2 * 8 + 3] * c;
2622         fw[i * 2 * 8 + 4] = fw[i * 2 * 8 + 4] * c;
2623         fw[i * 2 * 8 + 5] = fw[i * 2 * 8 + 5] * c;
2624         fw[i * 2 * 8 + 6] = fw[i * 2 * 8 + 6] * c;
2625         fw[i * 2 * 8 + 7] = fw[i * 2 * 8 + 7] * c;
2626     }
2627 }
2628
2629 static void opj_v8dwt_decode_step2(opj_v8_t* l, opj_v8_t* w,
2630                                    OPJ_UINT32 start,
2631                                    OPJ_UINT32 end,
2632                                    OPJ_UINT32 m,
2633                                    OPJ_FLOAT32 c)
2634 {
2635     OPJ_FLOAT32* fl = (OPJ_FLOAT32*) l;
2636     OPJ_FLOAT32* fw = (OPJ_FLOAT32*) w;
2637     OPJ_UINT32 i;
2638     OPJ_UINT32 imax = opj_uint_min(end, m);
2639     if (start > 0) {
2640         fw += 2 * NB_ELTS_V8 * start;
2641         fl = fw - 2 * NB_ELTS_V8;
2642     }
2643     /* To be adapted if NB_ELTS_V8 changes */
2644     for (i = start; i < imax; ++i) {
2645         fw[-8] = fw[-8] + ((fl[0] + fw[0]) * c);
2646         fw[-7] = fw[-7] + ((fl[1] + fw[1]) * c);
2647         fw[-6] = fw[-6] + ((fl[2] + fw[2]) * c);
2648         fw[-5] = fw[-5] + ((fl[3] + fw[3]) * c);
2649         fw[-4] = fw[-4] + ((fl[4] + fw[4]) * c);
2650         fw[-3] = fw[-3] + ((fl[5] + fw[5]) * c);
2651         fw[-2] = fw[-2] + ((fl[6] + fw[6]) * c);
2652         fw[-1] = fw[-1] + ((fl[7] + fw[7]) * c);
2653         fl = fw;
2654         fw += 2 * NB_ELTS_V8;
2655     }
2656     if (m < end) {
2657         assert(m + 1 == end);
2658         c += c;
2659         fw[-8] = fw[-8] + fl[0] * c;
2660         fw[-7] = fw[-7] + fl[1] * c;
2661         fw[-6] = fw[-6] + fl[2] * c;
2662         fw[-5] = fw[-5] + fl[3] * c;
2663         fw[-4] = fw[-4] + fl[4] * c;
2664         fw[-3] = fw[-3] + fl[5] * c;
2665         fw[-2] = fw[-2] + fl[6] * c;
2666         fw[-1] = fw[-1] + fl[7] * c;
2667     }
2668 }
2669
2670 #endif
2671
2672 /* <summary>                             */
2673 /* Inverse 9-7 wavelet transform in 1-D. */
2674 /* </summary>                            */
2675 static void opj_v8dwt_decode(opj_v8dwt_t* OPJ_RESTRICT dwt)
2676 {
2677     OPJ_INT32 a, b;
2678     /* BUG_WEIRD_TWO_INVK (look for this identifier in tcd.c) */
2679     /* Historic value for 2 / opj_invK */
2680     /* Normally, we should use invK, but if we do so, we have failures in the */
2681     /* conformance test, due to MSE and peak errors significantly higher than */
2682     /* accepted value */
2683     /* Due to using two_invK instead of invK, we have to compensate in tcd.c */
2684     /* the computation of the stepsize for the non LL subbands */
2685     const float two_invK = 1.625732422f;
2686     if (dwt->cas == 0) {
2687         if (!((dwt->dn > 0) || (dwt->sn > 1))) {
2688             return;
2689         }
2690         a = 0;
2691         b = 1;
2692     } else {
2693         if (!((dwt->sn > 0) || (dwt->dn > 1))) {
2694             return;
2695         }
2696         a = 1;
2697         b = 0;
2698     }
2699 #ifdef __SSE__
2700     opj_v8dwt_decode_step1_sse(dwt->wavelet + a, dwt->win_l_x0, dwt->win_l_x1,
2701                                _mm_set1_ps(opj_K));
2702     opj_v8dwt_decode_step1_sse(dwt->wavelet + b, dwt->win_h_x0, dwt->win_h_x1,
2703                                _mm_set1_ps(two_invK));
2704     opj_v8dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1,
2705                                dwt->win_l_x0, dwt->win_l_x1,
2706                                (OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a),
2707                                _mm_set1_ps(-opj_dwt_delta));
2708     opj_v8dwt_decode_step2_sse(dwt->wavelet + a, dwt->wavelet + b + 1,
2709                                dwt->win_h_x0, dwt->win_h_x1,
2710                                (OPJ_UINT32)opj_int_min(dwt->dn, dwt->sn - b),
2711                                _mm_set1_ps(-opj_dwt_gamma));
2712     opj_v8dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1,
2713                                dwt->win_l_x0, dwt->win_l_x1,
2714                                (OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a),
2715                                _mm_set1_ps(-opj_dwt_beta));
2716     opj_v8dwt_decode_step2_sse(dwt->wavelet + a, dwt->wavelet + b + 1,
2717                                dwt->win_h_x0, dwt->win_h_x1,
2718                                (OPJ_UINT32)opj_int_min(dwt->dn, dwt->sn - b),
2719                                _mm_set1_ps(-opj_dwt_alpha));
2720 #else
2721     opj_v8dwt_decode_step1(dwt->wavelet + a, dwt->win_l_x0, dwt->win_l_x1,
2722                            opj_K);
2723     opj_v8dwt_decode_step1(dwt->wavelet + b, dwt->win_h_x0, dwt->win_h_x1,
2724                            two_invK);
2725     opj_v8dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1,
2726                            dwt->win_l_x0, dwt->win_l_x1,
2727                            (OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a),
2728                            -opj_dwt_delta);
2729     opj_v8dwt_decode_step2(dwt->wavelet + a, dwt->wavelet + b + 1,
2730                            dwt->win_h_x0, dwt->win_h_x1,
2731                            (OPJ_UINT32)opj_int_min(dwt->dn, dwt->sn - b),
2732                            -opj_dwt_gamma);
2733     opj_v8dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1,
2734                            dwt->win_l_x0, dwt->win_l_x1,
2735                            (OPJ_UINT32)opj_int_min(dwt->sn, dwt->dn - a),
2736                            -opj_dwt_beta);
2737     opj_v8dwt_decode_step2(dwt->wavelet + a, dwt->wavelet + b + 1,
2738                            dwt->win_h_x0, dwt->win_h_x1,
2739                            (OPJ_UINT32)opj_int_min(dwt->dn, dwt->sn - b),
2740                            -opj_dwt_alpha);
2741 #endif
2742 }
2743
2744 typedef struct {
2745     opj_v8dwt_t h;
2746     OPJ_UINT32 rw;
2747     OPJ_UINT32 w;
2748     OPJ_FLOAT32 * OPJ_RESTRICT aj;
2749     OPJ_UINT32 nb_rows;
2750 } opj_dwt97_decode_h_job_t;
2751
2752 static void opj_dwt97_decode_h_func(void* user_data, opj_tls_t* tls)
2753 {
2754     OPJ_UINT32 j;
2755     opj_dwt97_decode_h_job_t* job;
2756     OPJ_FLOAT32 * OPJ_RESTRICT aj;
2757     OPJ_UINT32 w;
2758     (void)tls;
2759
2760     job = (opj_dwt97_decode_h_job_t*)user_data;
2761     w = job->w;
2762
2763     assert((job->nb_rows % NB_ELTS_V8) == 0);
2764
2765     aj = job->aj;
2766     for (j = 0; j + NB_ELTS_V8 <= job->nb_rows; j += NB_ELTS_V8) {
2767         OPJ_UINT32 k;
2768         opj_v8dwt_interleave_h(&job->h, aj, job->w, NB_ELTS_V8);
2769         opj_v8dwt_decode(&job->h);
2770
2771         /* To be adapted if NB_ELTS_V8 changes */
2772         for (k = 0; k < job->rw; k++) {
2773             aj[k      ] = job->h.wavelet[k].f[0];
2774             aj[k + (OPJ_SIZE_T)w  ] = job->h.wavelet[k].f[1];
2775             aj[k + (OPJ_SIZE_T)w * 2] = job->h.wavelet[k].f[2];
2776             aj[k + (OPJ_SIZE_T)w * 3] = job->h.wavelet[k].f[3];
2777         }
2778         for (k = 0; k < job->rw; k++) {
2779             aj[k + (OPJ_SIZE_T)w * 4] = job->h.wavelet[k].f[4];
2780             aj[k + (OPJ_SIZE_T)w * 5] = job->h.wavelet[k].f[5];
2781             aj[k + (OPJ_SIZE_T)w * 6] = job->h.wavelet[k].f[6];
2782             aj[k + (OPJ_SIZE_T)w * 7] = job->h.wavelet[k].f[7];
2783         }
2784
2785         aj += w * NB_ELTS_V8;
2786     }
2787
2788     opj_aligned_free(job->h.wavelet);
2789     opj_free(job);
2790 }
2791
2792
2793 typedef struct {
2794     opj_v8dwt_t v;
2795     OPJ_UINT32 rh;
2796     OPJ_UINT32 w;
2797     OPJ_FLOAT32 * OPJ_RESTRICT aj;
2798     OPJ_UINT32 nb_columns;
2799 } opj_dwt97_decode_v_job_t;
2800
2801 static void opj_dwt97_decode_v_func(void* user_data, opj_tls_t* tls)
2802 {
2803     OPJ_UINT32 j;
2804     opj_dwt97_decode_v_job_t* job;
2805     OPJ_FLOAT32 * OPJ_RESTRICT aj;
2806     (void)tls;
2807
2808     job = (opj_dwt97_decode_v_job_t*)user_data;
2809
2810     assert((job->nb_columns % NB_ELTS_V8) == 0);
2811
2812     aj = job->aj;
2813     for (j = 0; j + NB_ELTS_V8 <= job->nb_columns; j += NB_ELTS_V8) {
2814         OPJ_UINT32 k;
2815
2816         opj_v8dwt_interleave_v(&job->v, aj, job->w, NB_ELTS_V8);
2817         opj_v8dwt_decode(&job->v);
2818
2819         for (k = 0; k < job->rh; ++k) {
2820             memcpy(&aj[k * (OPJ_SIZE_T)job->w], &job->v.wavelet[k],
2821                    NB_ELTS_V8 * sizeof(OPJ_FLOAT32));
2822         }
2823         aj += NB_ELTS_V8;
2824     }
2825
2826     opj_aligned_free(job->v.wavelet);
2827     opj_free(job);
2828 }
2829
2830
2831 /* <summary>                             */
2832 /* Inverse 9-7 wavelet transform in 2-D. */
2833 /* </summary>                            */
2834 static
2835 OPJ_BOOL opj_dwt_decode_tile_97(opj_thread_pool_t* tp,
2836                                 opj_tcd_tilecomp_t* OPJ_RESTRICT tilec,
2837                                 OPJ_UINT32 numres)
2838 {
2839     opj_v8dwt_t h;
2840     opj_v8dwt_t v;
2841
2842     opj_tcd_resolution_t* res = tilec->resolutions;
2843
2844     OPJ_UINT32 rw = (OPJ_UINT32)(res->x1 -
2845                                  res->x0);    /* width of the resolution level computed */
2846     OPJ_UINT32 rh = (OPJ_UINT32)(res->y1 -
2847                                  res->y0);    /* height of the resolution level computed */
2848
2849     OPJ_UINT32 w = (OPJ_UINT32)(tilec->resolutions[tilec->minimum_num_resolutions -
2850                                                                1].x1 -
2851                                 tilec->resolutions[tilec->minimum_num_resolutions - 1].x0);
2852
2853     OPJ_SIZE_T l_data_size;
2854     const int num_threads = opj_thread_pool_get_thread_count(tp);
2855
2856     if (numres == 1) {
2857         return OPJ_TRUE;
2858     }
2859
2860     l_data_size = opj_dwt_max_resolution(res, numres);
2861     /* overflow check */
2862     if (l_data_size > (SIZE_MAX / sizeof(opj_v8_t))) {
2863         /* FIXME event manager error callback */
2864         return OPJ_FALSE;
2865     }
2866     h.wavelet = (opj_v8_t*) opj_aligned_malloc(l_data_size * sizeof(opj_v8_t));
2867     if (!h.wavelet) {
2868         /* FIXME event manager error callback */
2869         return OPJ_FALSE;
2870     }
2871     v.wavelet = h.wavelet;
2872
2873     while (--numres) {
2874         OPJ_FLOAT32 * OPJ_RESTRICT aj = (OPJ_FLOAT32*) tilec->data;
2875         OPJ_UINT32 j;
2876
2877         h.sn = (OPJ_INT32)rw;
2878         v.sn = (OPJ_INT32)rh;
2879
2880         ++res;
2881
2882         rw = (OPJ_UINT32)(res->x1 -
2883                           res->x0);   /* width of the resolution level computed */
2884         rh = (OPJ_UINT32)(res->y1 -
2885                           res->y0);   /* height of the resolution level computed */
2886
2887         h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
2888         h.cas = res->x0 % 2;
2889
2890         h.win_l_x0 = 0;
2891         h.win_l_x1 = (OPJ_UINT32)h.sn;
2892         h.win_h_x0 = 0;
2893         h.win_h_x1 = (OPJ_UINT32)h.dn;
2894
2895         if (num_threads <= 1 || rh < 2 * NB_ELTS_V8) {
2896             for (j = 0; j + (NB_ELTS_V8 - 1) < rh; j += NB_ELTS_V8) {
2897                 OPJ_UINT32 k;
2898                 opj_v8dwt_interleave_h(&h, aj, w, NB_ELTS_V8);
2899                 opj_v8dwt_decode(&h);
2900
2901                 /* To be adapted if NB_ELTS_V8 changes */
2902                 for (k = 0; k < rw; k++) {
2903                     aj[k      ] = h.wavelet[k].f[0];
2904                     aj[k + (OPJ_SIZE_T)w  ] = h.wavelet[k].f[1];
2905                     aj[k + (OPJ_SIZE_T)w * 2] = h.wavelet[k].f[2];
2906                     aj[k + (OPJ_SIZE_T)w * 3] = h.wavelet[k].f[3];
2907                 }
2908                 for (k = 0; k < rw; k++) {
2909                     aj[k + (OPJ_SIZE_T)w * 4] = h.wavelet[k].f[4];
2910                     aj[k + (OPJ_SIZE_T)w * 5] = h.wavelet[k].f[5];
2911                     aj[k + (OPJ_SIZE_T)w * 6] = h.wavelet[k].f[6];
2912                     aj[k + (OPJ_SIZE_T)w * 7] = h.wavelet[k].f[7];
2913                 }
2914
2915                 aj += w * NB_ELTS_V8;
2916             }
2917         } else {
2918             OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
2919             OPJ_UINT32 step_j;
2920
2921             if ((rh / NB_ELTS_V8) < num_jobs) {
2922                 num_jobs = rh / NB_ELTS_V8;
2923             }
2924             step_j = ((rh / num_jobs) / NB_ELTS_V8) * NB_ELTS_V8;
2925             for (j = 0; j < num_jobs; j++) {
2926                 opj_dwt97_decode_h_job_t* job;
2927
2928                 job = (opj_dwt97_decode_h_job_t*) opj_malloc(sizeof(opj_dwt97_decode_h_job_t));
2929                 if (!job) {
2930                     opj_thread_pool_wait_completion(tp, 0);
2931                     opj_aligned_free(h.wavelet);
2932                     return OPJ_FALSE;
2933                 }
2934                 job->h.wavelet = (opj_v8_t*)opj_aligned_malloc(l_data_size * sizeof(opj_v8_t));
2935                 if (!job->h.wavelet) {
2936                     opj_thread_pool_wait_completion(tp, 0);
2937                     opj_free(job);
2938                     opj_aligned_free(h.wavelet);
2939                     return OPJ_FALSE;
2940                 }
2941                 job->h.dn = h.dn;
2942                 job->h.sn = h.sn;
2943                 job->h.cas = h.cas;
2944                 job->h.win_l_x0 = h.win_l_x0;
2945                 job->h.win_l_x1 = h.win_l_x1;
2946                 job->h.win_h_x0 = h.win_h_x0;
2947                 job->h.win_h_x1 = h.win_h_x1;
2948                 job->rw = rw;
2949                 job->w = w;
2950                 job->aj = aj;
2951                 job->nb_rows = (j + 1 == num_jobs) ? (rh & (OPJ_UINT32)~
2952                                                       (NB_ELTS_V8 - 1)) - j * step_j : step_j;
2953                 aj += w * job->nb_rows;
2954                 opj_thread_pool_submit_job(tp, opj_dwt97_decode_h_func, job);
2955             }
2956             opj_thread_pool_wait_completion(tp, 0);
2957             j = rh & (OPJ_UINT32)~(NB_ELTS_V8 - 1);
2958         }
2959
2960         if (j < rh) {
2961             OPJ_UINT32 k;
2962             opj_v8dwt_interleave_h(&h, aj, w, rh - j);
2963             opj_v8dwt_decode(&h);
2964             for (k = 0; k < rw; k++) {
2965                 OPJ_UINT32 l;
2966                 for (l = 0; l < rh - j; l++) {
2967                     aj[k + (OPJ_SIZE_T)w  * l ] = h.wavelet[k].f[l];
2968                 }
2969             }
2970         }
2971
2972         v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
2973         v.cas = res->y0 % 2;
2974         v.win_l_x0 = 0;
2975         v.win_l_x1 = (OPJ_UINT32)v.sn;
2976         v.win_h_x0 = 0;
2977         v.win_h_x1 = (OPJ_UINT32)v.dn;
2978
2979         aj = (OPJ_FLOAT32*) tilec->data;
2980         if (num_threads <= 1 || rw < 2 * NB_ELTS_V8) {
2981             for (j = rw; j > (NB_ELTS_V8 - 1); j -= NB_ELTS_V8) {
2982                 OPJ_UINT32 k;
2983
2984                 opj_v8dwt_interleave_v(&v, aj, w, NB_ELTS_V8);
2985                 opj_v8dwt_decode(&v);
2986
2987                 for (k = 0; k < rh; ++k) {
2988                     memcpy(&aj[k * (OPJ_SIZE_T)w], &v.wavelet[k], NB_ELTS_V8 * sizeof(OPJ_FLOAT32));
2989                 }
2990                 aj += NB_ELTS_V8;
2991             }
2992         } else {
2993             /* "bench_dwt -I" shows that scaling is poor, likely due to RAM
2994                 transfer being the limiting factor. So limit the number of
2995                 threads.
2996              */
2997             OPJ_UINT32 num_jobs = opj_uint_max((OPJ_UINT32)num_threads / 2, 2U);
2998             OPJ_UINT32 step_j;
2999
3000             if ((rw / NB_ELTS_V8) < num_jobs) {
3001                 num_jobs = rw / NB_ELTS_V8;
3002             }
3003             step_j = ((rw / num_jobs) / NB_ELTS_V8) * NB_ELTS_V8;
3004             for (j = 0; j < num_jobs; j++) {
3005                 opj_dwt97_decode_v_job_t* job;
3006
3007                 job = (opj_dwt97_decode_v_job_t*) opj_malloc(sizeof(opj_dwt97_decode_v_job_t));
3008                 if (!job) {
3009                     opj_thread_pool_wait_completion(tp, 0);
3010                     opj_aligned_free(h.wavelet);
3011                     return OPJ_FALSE;
3012                 }
3013                 job->v.wavelet = (opj_v8_t*)opj_aligned_malloc(l_data_size * sizeof(opj_v8_t));
3014                 if (!job->v.wavelet) {
3015                     opj_thread_pool_wait_completion(tp, 0);
3016                     opj_free(job);
3017                     opj_aligned_free(h.wavelet);
3018                     return OPJ_FALSE;
3019                 }
3020                 job->v.dn = v.dn;
3021                 job->v.sn = v.sn;
3022                 job->v.cas = v.cas;
3023                 job->v.win_l_x0 = v.win_l_x0;
3024                 job->v.win_l_x1 = v.win_l_x1;
3025                 job->v.win_h_x0 = v.win_h_x0;
3026                 job->v.win_h_x1 = v.win_h_x1;
3027                 job->rh = rh;
3028                 job->w = w;
3029                 job->aj = aj;
3030                 job->nb_columns = (j + 1 == num_jobs) ? (rw & (OPJ_UINT32)~
3031                                   (NB_ELTS_V8 - 1)) - j * step_j : step_j;
3032                 aj += job->nb_columns;
3033                 opj_thread_pool_submit_job(tp, opj_dwt97_decode_v_func, job);
3034             }
3035             opj_thread_pool_wait_completion(tp, 0);
3036         }
3037
3038         if (rw & (NB_ELTS_V8 - 1)) {
3039             OPJ_UINT32 k;
3040
3041             j = rw & (NB_ELTS_V8 - 1);
3042
3043             opj_v8dwt_interleave_v(&v, aj, w, j);
3044             opj_v8dwt_decode(&v);
3045
3046             for (k = 0; k < rh; ++k) {
3047                 memcpy(&aj[k * (OPJ_SIZE_T)w], &v.wavelet[k],
3048                        (OPJ_SIZE_T)j * sizeof(OPJ_FLOAT32));
3049             }
3050         }
3051     }
3052
3053     opj_aligned_free(h.wavelet);
3054     return OPJ_TRUE;
3055 }
3056
3057 static
3058 OPJ_BOOL opj_dwt_decode_partial_97(opj_tcd_tilecomp_t* OPJ_RESTRICT tilec,
3059                                    OPJ_UINT32 numres)
3060 {
3061     opj_sparse_array_int32_t* sa;
3062     opj_v8dwt_t h;
3063     opj_v8dwt_t v;
3064     OPJ_UINT32 resno;
3065     /* This value matches the maximum left/right extension given in tables */
3066     /* F.2 and F.3 of the standard. Note: in opj_tcd_is_subband_area_of_interest() */
3067     /* we currently use 3. */
3068     const OPJ_UINT32 filter_width = 4U;
3069
3070     opj_tcd_resolution_t* tr = tilec->resolutions;
3071     opj_tcd_resolution_t* tr_max = &(tilec->resolutions[numres - 1]);
3072
3073     OPJ_UINT32 rw = (OPJ_UINT32)(tr->x1 -
3074                                  tr->x0);    /* width of the resolution level computed */
3075     OPJ_UINT32 rh = (OPJ_UINT32)(tr->y1 -
3076                                  tr->y0);    /* height of the resolution level computed */
3077
3078     OPJ_SIZE_T l_data_size;
3079
3080     /* Compute the intersection of the area of interest, expressed in tile coordinates */
3081     /* with the tile coordinates */
3082     OPJ_UINT32 win_tcx0 = tilec->win_x0;
3083     OPJ_UINT32 win_tcy0 = tilec->win_y0;
3084     OPJ_UINT32 win_tcx1 = tilec->win_x1;
3085     OPJ_UINT32 win_tcy1 = tilec->win_y1;
3086
3087     if (tr_max->x0 == tr_max->x1 || tr_max->y0 == tr_max->y1) {
3088         return OPJ_TRUE;
3089     }
3090
3091     sa = opj_dwt_init_sparse_array(tilec, numres);
3092     if (sa == NULL) {
3093         return OPJ_FALSE;
3094     }
3095
3096     if (numres == 1U) {
3097         OPJ_BOOL ret = opj_sparse_array_int32_read(sa,
3098                        tr_max->win_x0 - (OPJ_UINT32)tr_max->x0,
3099                        tr_max->win_y0 - (OPJ_UINT32)tr_max->y0,
3100                        tr_max->win_x1 - (OPJ_UINT32)tr_max->x0,
3101                        tr_max->win_y1 - (OPJ_UINT32)tr_max->y0,
3102                        tilec->data_win,
3103                        1, tr_max->win_x1 - tr_max->win_x0,
3104                        OPJ_TRUE);
3105         assert(ret);
3106         OPJ_UNUSED(ret);
3107         opj_sparse_array_int32_free(sa);
3108         return OPJ_TRUE;
3109     }
3110
3111     l_data_size = opj_dwt_max_resolution(tr, numres);
3112     /* overflow check */
3113     if (l_data_size > (SIZE_MAX / sizeof(opj_v8_t))) {
3114         /* FIXME event manager error callback */
3115         opj_sparse_array_int32_free(sa);
3116         return OPJ_FALSE;
3117     }
3118     h.wavelet = (opj_v8_t*) opj_aligned_malloc(l_data_size * sizeof(opj_v8_t));
3119     if (!h.wavelet) {
3120         /* FIXME event manager error callback */
3121         opj_sparse_array_int32_free(sa);
3122         return OPJ_FALSE;
3123     }
3124     v.wavelet = h.wavelet;
3125
3126     for (resno = 1; resno < numres; resno ++) {
3127         OPJ_UINT32 j;
3128         /* Window of interest subband-based coordinates */
3129         OPJ_UINT32 win_ll_x0, win_ll_y0, win_ll_x1, win_ll_y1;
3130         OPJ_UINT32 win_hl_x0, win_hl_x1;
3131         OPJ_UINT32 win_lh_y0, win_lh_y1;
3132         /* Window of interest tile-resolution-based coordinates */
3133         OPJ_UINT32 win_tr_x0, win_tr_x1, win_tr_y0, win_tr_y1;
3134         /* Tile-resolution subband-based coordinates */
3135         OPJ_UINT32 tr_ll_x0, tr_ll_y0, tr_hl_x0, tr_lh_y0;
3136
3137         ++tr;
3138
3139         h.sn = (OPJ_INT32)rw;
3140         v.sn = (OPJ_INT32)rh;
3141
3142         rw = (OPJ_UINT32)(tr->x1 - tr->x0);
3143         rh = (OPJ_UINT32)(tr->y1 - tr->y0);
3144
3145         h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
3146         h.cas = tr->x0 % 2;
3147
3148         v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
3149         v.cas = tr->y0 % 2;
3150
3151         /* Get the subband coordinates for the window of interest */
3152         /* LL band */
3153         opj_dwt_get_band_coordinates(tilec, resno, 0,
3154                                      win_tcx0, win_tcy0, win_tcx1, win_tcy1,
3155                                      &win_ll_x0, &win_ll_y0,
3156                                      &win_ll_x1, &win_ll_y1);
3157
3158         /* HL band */
3159         opj_dwt_get_band_coordinates(tilec, resno, 1,
3160                                      win_tcx0, win_tcy0, win_tcx1, win_tcy1,
3161                                      &win_hl_x0, NULL, &win_hl_x1, NULL);
3162
3163         /* LH band */
3164         opj_dwt_get_band_coordinates(tilec, resno, 2,
3165                                      win_tcx0, win_tcy0, win_tcx1, win_tcy1,
3166                                      NULL, &win_lh_y0, NULL, &win_lh_y1);
3167
3168         /* Beware: band index for non-LL0 resolution are 0=HL, 1=LH and 2=HH */
3169         tr_ll_x0 = (OPJ_UINT32)tr->bands[1].x0;
3170         tr_ll_y0 = (OPJ_UINT32)tr->bands[0].y0;
3171         tr_hl_x0 = (OPJ_UINT32)tr->bands[0].x0;
3172         tr_lh_y0 = (OPJ_UINT32)tr->bands[1].y0;
3173
3174         /* Subtract the origin of the bands for this tile, to the subwindow */
3175         /* of interest band coordinates, so as to get them relative to the */
3176         /* tile */
3177         win_ll_x0 = opj_uint_subs(win_ll_x0, tr_ll_x0);
3178         win_ll_y0 = opj_uint_subs(win_ll_y0, tr_ll_y0);
3179         win_ll_x1 = opj_uint_subs(win_ll_x1, tr_ll_x0);
3180         win_ll_y1 = opj_uint_subs(win_ll_y1, tr_ll_y0);
3181         win_hl_x0 = opj_uint_subs(win_hl_x0, tr_hl_x0);
3182         win_hl_x1 = opj_uint_subs(win_hl_x1, tr_hl_x0);
3183         win_lh_y0 = opj_uint_subs(win_lh_y0, tr_lh_y0);
3184         win_lh_y1 = opj_uint_subs(win_lh_y1, tr_lh_y0);
3185
3186         opj_dwt_segment_grow(filter_width, (OPJ_UINT32)h.sn, &win_ll_x0, &win_ll_x1);
3187         opj_dwt_segment_grow(filter_width, (OPJ_UINT32)h.dn, &win_hl_x0, &win_hl_x1);
3188
3189         opj_dwt_segment_grow(filter_width, (OPJ_UINT32)v.sn, &win_ll_y0, &win_ll_y1);
3190         opj_dwt_segment_grow(filter_width, (OPJ_UINT32)v.dn, &win_lh_y0, &win_lh_y1);
3191
3192         /* Compute the tile-resolution-based coordinates for the window of interest */
3193         if (h.cas == 0) {
3194             win_tr_x0 = opj_uint_min(2 * win_ll_x0, 2 * win_hl_x0 + 1);
3195             win_tr_x1 = opj_uint_min(opj_uint_max(2 * win_ll_x1, 2 * win_hl_x1 + 1), rw);
3196         } else {
3197             win_tr_x0 = opj_uint_min(2 * win_hl_x0, 2 * win_ll_x0 + 1);
3198             win_tr_x1 = opj_uint_min(opj_uint_max(2 * win_hl_x1, 2 * win_ll_x1 + 1), rw);
3199         }
3200
3201         if (v.cas == 0) {
3202             win_tr_y0 = opj_uint_min(2 * win_ll_y0, 2 * win_lh_y0 + 1);
3203             win_tr_y1 = opj_uint_min(opj_uint_max(2 * win_ll_y1, 2 * win_lh_y1 + 1), rh);
3204         } else {
3205             win_tr_y0 = opj_uint_min(2 * win_lh_y0, 2 * win_ll_y0 + 1);
3206             win_tr_y1 = opj_uint_min(opj_uint_max(2 * win_lh_y1, 2 * win_ll_y1 + 1), rh);
3207         }
3208
3209         h.win_l_x0 = win_ll_x0;
3210         h.win_l_x1 = win_ll_x1;
3211         h.win_h_x0 = win_hl_x0;
3212         h.win_h_x1 = win_hl_x1;
3213         for (j = 0; j + (NB_ELTS_V8 - 1) < rh; j += NB_ELTS_V8) {
3214             if ((j + (NB_ELTS_V8 - 1) >= win_ll_y0 && j < win_ll_y1) ||
3215                     (j + (NB_ELTS_V8 - 1) >= win_lh_y0 + (OPJ_UINT32)v.sn &&
3216                      j < win_lh_y1 + (OPJ_UINT32)v.sn)) {
3217                 opj_v8dwt_interleave_partial_h(&h, sa, j, opj_uint_min(NB_ELTS_V8, rh - j));
3218                 opj_v8dwt_decode(&h);
3219                 if (!opj_sparse_array_int32_write(sa,
3220                                                   win_tr_x0, j,
3221                                                   win_tr_x1, j + NB_ELTS_V8,
3222                                                   (OPJ_INT32*)&h.wavelet[win_tr_x0].f[0],
3223                                                   NB_ELTS_V8, 1, OPJ_TRUE)) {
3224                     /* FIXME event manager error callback */
3225                     opj_sparse_array_int32_free(sa);
3226                     opj_aligned_free(h.wavelet);
3227                     return OPJ_FALSE;
3228                 }
3229             }
3230         }
3231
3232         if (j < rh &&
3233                 ((j + (NB_ELTS_V8 - 1) >= win_ll_y0 && j < win_ll_y1) ||
3234                  (j + (NB_ELTS_V8 - 1) >= win_lh_y0 + (OPJ_UINT32)v.sn &&
3235                   j < win_lh_y1 + (OPJ_UINT32)v.sn))) {
3236             opj_v8dwt_interleave_partial_h(&h, sa, j, rh - j);
3237             opj_v8dwt_decode(&h);
3238             if (!opj_sparse_array_int32_write(sa,
3239                                               win_tr_x0, j,
3240                                               win_tr_x1, rh,
3241                                               (OPJ_INT32*)&h.wavelet[win_tr_x0].f[0],
3242                                               NB_ELTS_V8, 1, OPJ_TRUE)) {
3243                 /* FIXME event manager error callback */
3244                 opj_sparse_array_int32_free(sa);
3245                 opj_aligned_free(h.wavelet);
3246                 return OPJ_FALSE;
3247             }
3248         }
3249
3250         v.win_l_x0 = win_ll_y0;
3251         v.win_l_x1 = win_ll_y1;
3252         v.win_h_x0 = win_lh_y0;
3253         v.win_h_x1 = win_lh_y1;
3254         for (j = win_tr_x0; j < win_tr_x1; j += NB_ELTS_V8) {
3255             OPJ_UINT32 nb_elts = opj_uint_min(NB_ELTS_V8, win_tr_x1 - j);
3256
3257             opj_v8dwt_interleave_partial_v(&v, sa, j, nb_elts);
3258             opj_v8dwt_decode(&v);
3259
3260             if (!opj_sparse_array_int32_write(sa,
3261                                               j, win_tr_y0,
3262                                               j + nb_elts, win_tr_y1,
3263                                               (OPJ_INT32*)&h.wavelet[win_tr_y0].f[0],
3264                                               1, NB_ELTS_V8, OPJ_TRUE)) {
3265                 /* FIXME event manager error callback */
3266                 opj_sparse_array_int32_free(sa);
3267                 opj_aligned_free(h.wavelet);
3268                 return OPJ_FALSE;
3269             }
3270         }
3271     }
3272
3273     {
3274         OPJ_BOOL ret = opj_sparse_array_int32_read(sa,
3275                        tr_max->win_x0 - (OPJ_UINT32)tr_max->x0,
3276                        tr_max->win_y0 - (OPJ_UINT32)tr_max->y0,
3277                        tr_max->win_x1 - (OPJ_UINT32)tr_max->x0,
3278                        tr_max->win_y1 - (OPJ_UINT32)tr_max->y0,
3279                        tilec->data_win,
3280                        1, tr_max->win_x1 - tr_max->win_x0,
3281                        OPJ_TRUE);
3282         assert(ret);
3283         OPJ_UNUSED(ret);
3284     }
3285     opj_sparse_array_int32_free(sa);
3286
3287     opj_aligned_free(h.wavelet);
3288     return OPJ_TRUE;
3289 }
3290
3291
3292 OPJ_BOOL opj_dwt_decode_real(opj_tcd_t *p_tcd,
3293                              opj_tcd_tilecomp_t* OPJ_RESTRICT tilec,
3294                              OPJ_UINT32 numres)
3295 {
3296     if (p_tcd->whole_tile_decoding) {
3297         return opj_dwt_decode_tile_97(p_tcd->thread_pool, tilec, numres);
3298     } else {
3299         return opj_dwt_decode_partial_97(tilec, numres);
3300     }
3301 }