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