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  * All rights reserved.
  17  *
  18  * Redistribution and use in source and binary forms, with or without
  19  * modification, are permitted provided that the following conditions
  20  * are met:
  21  * 1. Redistributions of source code must retain the above copyright
  22  *    notice, this list of conditions and the following disclaimer.
  23  * 2. Redistributions in binary form must reproduce the above copyright
  24  *    notice, this list of conditions and the following disclaimer in the
  25  *    documentation and/or other materials provided with the distribution.
  26  *
  27  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
  28  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  29  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  30  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  31  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  32  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  33  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  34  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  35  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  36  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  37  * POSSIBILITY OF SUCH DAMAGE.
  38  */
  39
  40 #ifdef __SSE__
  41 #include <xmmintrin.h>
  42 #endif
  43
  44 #include "opj_includes.h"
  45
  46 /** @defgroup DWT DWT - Implementation of a discrete wavelet transform */
  47 /*@{*/
  48
  49 #define OPJ_WS(i) v->mem[(i)*2]
  50 #define OPJ_WD(i) v->mem[(1+(i)*2)]
  51
  52 /** @name Local data structures */
  53 /*@{*/
  54
  55 typedef struct dwt_local {
  56     OPJ_INT32* mem;
  57     OPJ_INT32 dn;
  58     OPJ_INT32 sn;
  59     OPJ_INT32 cas;
  60 } opj_dwt_t;
  61
  62 typedef union {
  63     OPJ_FLOAT32 f[4];
  64 } opj_v4_t;
  65
  66 typedef struct v4dwt_local {
  67     opj_v4_t*   wavelet ;
  68     OPJ_INT32       dn ;
  69     OPJ_INT32       sn ;
  70     OPJ_INT32       cas ;
  71 } opj_v4dwt_t ;
  72
  73 static const OPJ_FLOAT32 opj_dwt_alpha =  1.586134342f; /*  12994 */
  74 static const OPJ_FLOAT32 opj_dwt_beta  =  0.052980118f; /*    434 */
  75 static const OPJ_FLOAT32 opj_dwt_gamma = -0.882911075f; /*  -7233 */
  76 static const OPJ_FLOAT32 opj_dwt_delta = -0.443506852f; /*  -3633 */
  77
  78 static const OPJ_FLOAT32 opj_K      = 1.230174105f; /*  10078 */
  79 static const OPJ_FLOAT32 opj_c13318 = 1.625732422f;
  80
  81 /*@}*/
  82
  83 /**
  84 Virtual function type for wavelet transform in 1-D
  85 */
  86 typedef void (*DWT1DFN)(opj_dwt_t* v);
  87
  88 /** @name Local static functions */
  89 /*@{*/
  90
  91 /**
  92 Forward lazy transform (horizontal)
  93 */
  94 static void opj_dwt_deinterleave_h(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
  95                                    OPJ_INT32 sn, OPJ_INT32 cas);
  96 /**
  97 Forward lazy transform (vertical)
  98 */
  99 static void opj_dwt_deinterleave_v(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
 100                                    OPJ_INT32 sn, OPJ_INT32 x, OPJ_INT32 cas);
 101 /**
 102 Inverse lazy transform (horizontal)
 103 */
 104 static void opj_dwt_interleave_h(opj_dwt_t* h, OPJ_INT32 *a);
 105 /**
 106 Inverse lazy transform (vertical)
 107 */
 108 static void opj_dwt_interleave_v(opj_dwt_t* v, OPJ_INT32 *a, OPJ_INT32 x);
 109 /**
 110 Forward 5-3 wavelet transform in 1-D
 111 */
 112 static void opj_dwt_encode_1(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
 113                              OPJ_INT32 cas);
 114 /**
 115 Inverse 5-3 wavelet transform in 1-D
 116 */
 117 static void opj_dwt_decode_1(opj_dwt_t *v);
 118 static void opj_dwt_decode_1_(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
 119                               OPJ_INT32 cas);
 120 /**
 121 Forward 9-7 wavelet transform in 1-D
 122 */
 123 static void opj_dwt_encode_1_real(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
 124                                   OPJ_INT32 cas);
 125 /**
 126 Explicit calculation of the Quantization Stepsizes
 127 */
 128 static void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps,
 129                                     opj_stepsize_t *bandno_stepsize);
 130 /**
 131 Inverse wavelet transform in 2-D.
 132 */
 133 static OPJ_BOOL opj_dwt_decode_tile(opj_thread_pool_t* tp,
 134                                     opj_tcd_tilecomp_t* tilec, OPJ_UINT32 i, DWT1DFN fn);
 135
 136 static OPJ_BOOL opj_dwt_encode_procedure(opj_tcd_tilecomp_t * tilec,
 137         void (*p_function)(OPJ_INT32 *, OPJ_INT32, OPJ_INT32, OPJ_INT32));
 138
 139 static OPJ_UINT32 opj_dwt_max_resolution(opj_tcd_resolution_t* OPJ_RESTRICT r,
 140         OPJ_UINT32 i);
 141
 142 /* <summary>                             */
 143 /* Inverse 9-7 wavelet transform in 1-D. */
 144 /* </summary>                            */
 145 static void opj_v4dwt_decode(opj_v4dwt_t* OPJ_RESTRICT dwt);
 146
 147 static void opj_v4dwt_interleave_h(opj_v4dwt_t* OPJ_RESTRICT w,
 148                                    OPJ_FLOAT32* OPJ_RESTRICT a, OPJ_INT32 x, OPJ_INT32 size);
 149
 150 static void opj_v4dwt_interleave_v(opj_v4dwt_t* OPJ_RESTRICT v,
 151                                    OPJ_FLOAT32* OPJ_RESTRICT a, OPJ_INT32 x, OPJ_INT32 nb_elts_read);
 152
 153 #ifdef __SSE__
 154 static void opj_v4dwt_decode_step1_sse(opj_v4_t* w, OPJ_INT32 count,
 155                                        const __m128 c);
 156
 157 static void opj_v4dwt_decode_step2_sse(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k,
 158                                        OPJ_INT32 m, __m128 c);
 159
 160 #else
 161 static void opj_v4dwt_decode_step1(opj_v4_t* w, OPJ_INT32 count,
 162                                    const OPJ_FLOAT32 c);
 163
 164 static void opj_v4dwt_decode_step2(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k,
 165                                    OPJ_INT32 m, OPJ_FLOAT32 c);
 166
 167 #endif
 168
 169 /*@}*/
 170
 171 /*@}*/
 172
 173 #define OPJ_S(i) a[(i)*2]
 174 #define OPJ_D(i) a[(1+(i)*2)]
 175 #define OPJ_S_(i) ((i)<0?OPJ_S(0):((i)>=sn?OPJ_S(sn-1):OPJ_S(i)))
 176 #define OPJ_D_(i) ((i)<0?OPJ_D(0):((i)>=dn?OPJ_D(dn-1):OPJ_D(i)))
 177 /* new */
 178 #define OPJ_SS_(i) ((i)<0?OPJ_S(0):((i)>=dn?OPJ_S(dn-1):OPJ_S(i)))
 179 #define OPJ_DD_(i) ((i)<0?OPJ_D(0):((i)>=sn?OPJ_D(sn-1):OPJ_D(i)))
 180
 181 /* <summary>                                                              */
 182 /* This table contains the norms of the 5-3 wavelets for different bands. */
 183 /* </summary>                                                             */
 184 static const OPJ_FLOAT64 opj_dwt_norms[4][10] = {
 185     {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 186     {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 187     {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 188     {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
 189 };
 190
 191 /* <summary>                                                              */
 192 /* This table contains the norms of the 9-7 wavelets for different bands. */
 193 /* </summary>                                                             */
 194 static const OPJ_FLOAT64 opj_dwt_norms_real[4][10] = {
 195     {1.000, 1.965, 4.177, 8.403, 16.90, 33.84, 67.69, 135.3, 270.6, 540.9},
 196     {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
 197     {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
 198     {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2}
 199 };
 200
 201 /*
 202 ==========================================================
 203    local functions
 204 ==========================================================
 205 */
 206
 207 /* <summary>                             */
 208 /* Forward lazy transform (horizontal).  */
 209 /* </summary>                            */
 210 static void opj_dwt_deinterleave_h(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
 211                                    OPJ_INT32 sn, OPJ_INT32 cas)
 212 {
 213     OPJ_INT32 i;
 214     OPJ_INT32 * l_dest = b;
 215     OPJ_INT32 * l_src = a + cas;
 216
 217     for (i = 0; i < sn; ++i) {
 218         *l_dest++ = *l_src;
 219         l_src += 2;
 220     }
 221
 222     l_dest = b + sn;
 223     l_src = a + 1 - cas;
 224
 225     for (i = 0; i < dn; ++i)  {
 226         *l_dest++ = *l_src;
 227         l_src += 2;
 228     }
 229 }
 230
 231 /* <summary>                             */
 232 /* Forward lazy transform (vertical).    */
 233 /* </summary>                            */
 234 static void opj_dwt_deinterleave_v(OPJ_INT32 *a, OPJ_INT32 *b, OPJ_INT32 dn,
 235                                    OPJ_INT32 sn, OPJ_INT32 x, OPJ_INT32 cas)
 236 {
 237     OPJ_INT32 i = sn;
 238     OPJ_INT32 * l_dest = b;
 239     OPJ_INT32 * l_src = a + cas;
 240
 241     while (i--) {
 242         *l_dest = *l_src;
 243         l_dest += x;
 244         l_src += 2;
 245     } /* b[i*x]=a[2*i+cas]; */
 246
 247     l_dest = b + sn * x;
 248     l_src = a + 1 - cas;
 249
 250     i = dn;
 251     while (i--) {
 252         *l_dest = *l_src;
 253         l_dest += x;
 254         l_src += 2;
 255     } /*b[(sn+i)*x]=a[(2*i+1-cas)];*/
 256 }
 257
 258 /* <summary>                             */
 259 /* Inverse lazy transform (horizontal).  */
 260 /* </summary>                            */
 261 static void opj_dwt_interleave_h(opj_dwt_t* h, OPJ_INT32 *a)
 262 {
 263     OPJ_INT32 *ai = a;
 264     OPJ_INT32 *bi = h->mem + h->cas;
 265     OPJ_INT32  i    = h->sn;
 266     while (i--) {
 267         *bi = *(ai++);
 268         bi += 2;
 269     }
 270     ai  = a + h->sn;
 271     bi  = h->mem + 1 - h->cas;
 272     i   = h->dn ;
 273     while (i--) {
 274         *bi = *(ai++);
 275         bi += 2;
 276     }
 277 }
 278
 279 /* <summary>                             */
 280 /* Inverse lazy transform (vertical).    */
 281 /* </summary>                            */
 282 static void opj_dwt_interleave_v(opj_dwt_t* v, OPJ_INT32 *a, OPJ_INT32 x)
 283 {
 284     OPJ_INT32 *ai = a;
 285     OPJ_INT32 *bi = v->mem + v->cas;
 286     OPJ_INT32  i = v->sn;
 287     while (i--) {
 288         *bi = *ai;
 289         bi += 2;
 290         ai += x;
 291     }
 292     ai = a + (v->sn * x);
 293     bi = v->mem + 1 - v->cas;
 294     i = v->dn ;
 295     while (i--) {
 296         *bi = *ai;
 297         bi += 2;
 298         ai += x;
 299     }
 300 }
 301
 302
 303 /* <summary>                            */
 304 /* Forward 5-3 wavelet transform in 1-D. */
 305 /* </summary>                           */
 306 static void opj_dwt_encode_1(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
 307                              OPJ_INT32 cas)
 308 {
 309     OPJ_INT32 i;
 310
 311     if (!cas) {
 312         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 313             for (i = 0; i < dn; i++) {
 314                 OPJ_D(i) -= (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
 315             }
 316             for (i = 0; i < sn; i++) {
 317                 OPJ_S(i) += (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
 318             }
 319         }
 320     } else {
 321         if (!sn && dn == 1) {       /* NEW :  CASE ONE ELEMENT */
 322             OPJ_S(0) *= 2;
 323         } else {
 324             for (i = 0; i < dn; i++) {
 325                 OPJ_S(i) -= (OPJ_DD_(i) + OPJ_DD_(i - 1)) >> 1;
 326             }
 327             for (i = 0; i < sn; i++) {
 328                 OPJ_D(i) += (OPJ_SS_(i) + OPJ_SS_(i + 1) + 2) >> 2;
 329             }
 330         }
 331     }
 332 }
 333
 334 /* <summary>                            */
 335 /* Inverse 5-3 wavelet transform in 1-D. */
 336 /* </summary>                           */
 337 static void opj_dwt_decode_1_(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
 338                               OPJ_INT32 cas)
 339 {
 340     OPJ_INT32 i;
 341
 342     if (!cas) {
 343         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 344             for (i = 0; i < sn; i++) {
 345                 OPJ_S(i) -= (OPJ_D_(i - 1) + OPJ_D_(i) + 2) >> 2;
 346             }
 347             for (i = 0; i < dn; i++) {
 348                 OPJ_D(i) += (OPJ_S_(i) + OPJ_S_(i + 1)) >> 1;
 349             }
 350         }
 351     } else {
 352         if (!sn  && dn == 1) {        /* NEW :  CASE ONE ELEMENT */
 353             OPJ_S(0) /= 2;
 354         } else {
 355             for (i = 0; i < sn; i++) {
 356                 OPJ_D(i) -= (OPJ_SS_(i) + OPJ_SS_(i + 1) + 2) >> 2;
 357             }
 358             for (i = 0; i < dn; i++) {
 359                 OPJ_S(i) += (OPJ_DD_(i) + OPJ_DD_(i - 1)) >> 1;
 360             }
 361         }
 362     }
 363 }
 364
 365 /* <summary>                            */
 366 /* Inverse 5-3 wavelet transform in 1-D. */
 367 /* </summary>                           */
 368 static void opj_dwt_decode_1(opj_dwt_t *v)
 369 {
 370     opj_dwt_decode_1_(v->mem, v->dn, v->sn, v->cas);
 371 }
 372
 373 /* <summary>                             */
 374 /* Forward 9-7 wavelet transform in 1-D. */
 375 /* </summary>                            */
 376 static void opj_dwt_encode_1_real(OPJ_INT32 *a, OPJ_INT32 dn, OPJ_INT32 sn,
 377                                   OPJ_INT32 cas)
 378 {
 379     OPJ_INT32 i;
 380     if (!cas) {
 381         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 382             for (i = 0; i < dn; i++) {
 383                 OPJ_D(i) -= opj_int_fix_mul(OPJ_S_(i) + OPJ_S_(i + 1), 12993);
 384             }
 385             for (i = 0; i < sn; i++) {
 386                 OPJ_S(i) -= opj_int_fix_mul(OPJ_D_(i - 1) + OPJ_D_(i), 434);
 387             }
 388             for (i = 0; i < dn; i++) {
 389                 OPJ_D(i) += opj_int_fix_mul(OPJ_S_(i) + OPJ_S_(i + 1), 7233);
 390             }
 391             for (i = 0; i < sn; i++) {
 392                 OPJ_S(i) += opj_int_fix_mul(OPJ_D_(i - 1) + OPJ_D_(i), 3633);
 393             }
 394             for (i = 0; i < dn; i++) {
 395                 OPJ_D(i) = opj_int_fix_mul(OPJ_D(i), 5038);    /*5038 */
 396             }
 397             for (i = 0; i < sn; i++) {
 398                 OPJ_S(i) = opj_int_fix_mul(OPJ_S(i), 6659);    /*6660 */
 399             }
 400         }
 401     } else {
 402         if ((sn > 0) || (dn > 1)) { /* NEW :  CASE ONE ELEMENT */
 403             for (i = 0; i < dn; i++) {
 404                 OPJ_S(i) -= opj_int_fix_mul(OPJ_DD_(i) + OPJ_DD_(i - 1), 12993);
 405             }
 406             for (i = 0; i < sn; i++) {
 407                 OPJ_D(i) -= opj_int_fix_mul(OPJ_SS_(i) + OPJ_SS_(i + 1), 434);
 408             }
 409             for (i = 0; i < dn; i++) {
 410                 OPJ_S(i) += opj_int_fix_mul(OPJ_DD_(i) + OPJ_DD_(i - 1), 7233);
 411             }
 412             for (i = 0; i < sn; i++) {
 413                 OPJ_D(i) += opj_int_fix_mul(OPJ_SS_(i) + OPJ_SS_(i + 1), 3633);
 414             }
 415             for (i = 0; i < dn; i++) {
 416                 OPJ_S(i) = opj_int_fix_mul(OPJ_S(i), 5038);    /*5038 */
 417             }
 418             for (i = 0; i < sn; i++) {
 419                 OPJ_D(i) = opj_int_fix_mul(OPJ_D(i), 6659);    /*6660 */
 420             }
 421         }
 422     }
 423 }
 424
 425 static void opj_dwt_encode_stepsize(OPJ_INT32 stepsize, OPJ_INT32 numbps,
 426                                     opj_stepsize_t *bandno_stepsize)
 427 {
 428     OPJ_INT32 p, n;
 429     p = opj_int_floorlog2(stepsize) - 13;
 430     n = 11 - opj_int_floorlog2(stepsize);
 431     bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
 432     bandno_stepsize->expn = numbps - p;
 433 }
 434
 435 /*
 436 ==========================================================
 437    DWT interface
 438 ==========================================================
 439 */
 440
 441
 442 /* <summary>                            */
 443 /* Forward 5-3 wavelet transform in 2-D. */
 444 /* </summary>                           */
 445 static INLINE OPJ_BOOL opj_dwt_encode_procedure(opj_tcd_tilecomp_t * tilec,
 446         void (*p_function)(OPJ_INT32 *, OPJ_INT32, OPJ_INT32, OPJ_INT32))
 447 {
 448     OPJ_INT32 i, j, k;
 449     OPJ_INT32 *a = 00;
 450     OPJ_INT32 *aj = 00;
 451     OPJ_INT32 *bj = 00;
 452     OPJ_INT32 w, l;
 453
 454     OPJ_INT32 rw;           /* width of the resolution level computed   */
 455     OPJ_INT32 rh;           /* height of the resolution level computed  */
 456     size_t l_data_size;
 457
 458     opj_tcd_resolution_t * l_cur_res = 0;
 459     opj_tcd_resolution_t * l_last_res = 0;
 460
 461     w = tilec->x1 - tilec->x0;
 462     l = (OPJ_INT32)tilec->numresolutions - 1;
 463     a = tilec->data;
 464
 465     l_cur_res = tilec->resolutions + l;
 466     l_last_res = l_cur_res - 1;
 467
 468     l_data_size = opj_dwt_max_resolution(tilec->resolutions, tilec->numresolutions);
 469     /* overflow check */
 470     if (l_data_size > (SIZE_MAX / sizeof(OPJ_INT32))) {
 471         /* FIXME event manager error callback */
 472         return OPJ_FALSE;
 473     }
 474     l_data_size *= sizeof(OPJ_INT32);
 475     bj = (OPJ_INT32*)opj_malloc(l_data_size);
 476     /* l_data_size is equal to 0 when numresolutions == 1 but bj is not used */
 477     /* in that case, so do not error out */
 478     if (l_data_size != 0 && ! bj) {
 479         return OPJ_FALSE;
 480     }
 481     i = l;
 482
 483     while (i--) {
 484         OPJ_INT32 rw1;      /* width of the resolution level once lower than computed one                                       */
 485         OPJ_INT32 rh1;      /* height of the resolution level once lower than computed one                                      */
 486         OPJ_INT32 cas_col;  /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 487         OPJ_INT32 cas_row;  /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 488         OPJ_INT32 dn, sn;
 489
 490         rw  = l_cur_res->x1 - l_cur_res->x0;
 491         rh  = l_cur_res->y1 - l_cur_res->y0;
 492         rw1 = l_last_res->x1 - l_last_res->x0;
 493         rh1 = l_last_res->y1 - l_last_res->y0;
 494
 495         cas_row = l_cur_res->x0 & 1;
 496         cas_col = l_cur_res->y0 & 1;
 497
 498         sn = rh1;
 499         dn = rh - rh1;
 500         for (j = 0; j < rw; ++j) {
 501             aj = a + j;
 502             for (k = 0; k < rh; ++k) {
 503                 bj[k] = aj[k * w];
 504             }
 505
 506             (*p_function)(bj, dn, sn, cas_col);
 507
 508             opj_dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 509         }
 510
 511         sn = rw1;
 512         dn = rw - rw1;
 513
 514         for (j = 0; j < rh; j++) {
 515             aj = a + j * w;
 516             for (k = 0; k < rw; k++) {
 517                 bj[k] = aj[k];
 518             }
 519             (*p_function)(bj, dn, sn, cas_row);
 520             opj_dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 521         }
 522
 523         l_cur_res = l_last_res;
 524
 525         --l_last_res;
 526     }
 527
 528     opj_free(bj);
 529     return OPJ_TRUE;
 530 }
 531
 532 /* Forward 5-3 wavelet transform in 2-D. */
 533 /* </summary>                           */
 534 OPJ_BOOL opj_dwt_encode(opj_tcd_tilecomp_t * tilec)
 535 {
 536     return opj_dwt_encode_procedure(tilec, opj_dwt_encode_1);
 537 }
 538
 539 /* <summary>                            */
 540 /* Inverse 5-3 wavelet transform in 2-D. */
 541 /* </summary>                           */
 542 OPJ_BOOL opj_dwt_decode(opj_thread_pool_t* tp, opj_tcd_tilecomp_t* tilec,
 543                         OPJ_UINT32 numres)
 544 {
 545     return opj_dwt_decode_tile(tp, tilec, numres, &opj_dwt_decode_1);
 546 }
 547
 548
 549 /* <summary>                          */
 550 /* Get gain of 5-3 wavelet transform. */
 551 /* </summary>                         */
 552 OPJ_UINT32 opj_dwt_getgain(OPJ_UINT32 orient)
 553 {
 554     if (orient == 0) {
 555         return 0;
 556     }
 557     if (orient == 1 || orient == 2) {
 558         return 1;
 559     }
 560     return 2;
 561 }
 562
 563 /* <summary>                */
 564 /* Get norm of 5-3 wavelet. */
 565 /* </summary>               */
 566 OPJ_FLOAT64 opj_dwt_getnorm(OPJ_UINT32 level, OPJ_UINT32 orient)
 567 {
 568     return opj_dwt_norms[orient][level];
 569 }
 570
 571 /* <summary>                             */
 572 /* Forward 9-7 wavelet transform in 2-D. */
 573 /* </summary>                            */
 574 OPJ_BOOL opj_dwt_encode_real(opj_tcd_tilecomp_t * tilec)
 575 {
 576     return opj_dwt_encode_procedure(tilec, opj_dwt_encode_1_real);
 577 }
 578
 579 /* <summary>                          */
 580 /* Get gain of 9-7 wavelet transform. */
 581 /* </summary>                         */
 582 OPJ_UINT32 opj_dwt_getgain_real(OPJ_UINT32 orient)
 583 {
 584     (void)orient;
 585     return 0;
 586 }
 587
 588 /* <summary>                */
 589 /* Get norm of 9-7 wavelet. */
 590 /* </summary>               */
 591 OPJ_FLOAT64 opj_dwt_getnorm_real(OPJ_UINT32 level, OPJ_UINT32 orient)
 592 {
 593     return opj_dwt_norms_real[orient][level];
 594 }
 595
 596 void opj_dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, OPJ_UINT32 prec)
 597 {
 598     OPJ_UINT32 numbands, bandno;
 599     numbands = 3 * tccp->numresolutions - 2;
 600     for (bandno = 0; bandno < numbands; bandno++) {
 601         OPJ_FLOAT64 stepsize;
 602         OPJ_UINT32 resno, level, orient, gain;
 603
 604         resno = (bandno == 0) ? 0 : ((bandno - 1) / 3 + 1);
 605         orient = (bandno == 0) ? 0 : ((bandno - 1) % 3 + 1);
 606         level = tccp->numresolutions - 1 - resno;
 607         gain = (tccp->qmfbid == 0) ? 0 : ((orient == 0) ? 0 : (((orient == 1) ||
 608                                           (orient == 2)) ? 1 : 2));
 609         if (tccp->qntsty == J2K_CCP_QNTSTY_NOQNT) {
 610             stepsize = 1.0;
 611         } else {
 612             OPJ_FLOAT64 norm = opj_dwt_norms_real[orient][level];
 613             stepsize = (1 << (gain)) / norm;
 614         }
 615         opj_dwt_encode_stepsize((OPJ_INT32) floor(stepsize * 8192.0),
 616                                 (OPJ_INT32)(prec + gain), &tccp->stepsizes[bandno]);
 617     }
 618 }
 619
 620 /* <summary>                             */
 621 /* Determine maximum computed resolution level for inverse wavelet transform */
 622 /* </summary>                            */
 623 static OPJ_UINT32 opj_dwt_max_resolution(opj_tcd_resolution_t* OPJ_RESTRICT r,
 624         OPJ_UINT32 i)
 625 {
 626     OPJ_UINT32 mr   = 0;
 627     OPJ_UINT32 w;
 628     while (--i) {
 629         ++r;
 630         if (mr < (w = (OPJ_UINT32)(r->x1 - r->x0))) {
 631             mr = w ;
 632         }
 633         if (mr < (w = (OPJ_UINT32)(r->y1 - r->y0))) {
 634             mr = w ;
 635         }
 636     }
 637     return mr ;
 638 }
 639
 640 typedef struct {
 641     opj_dwt_t h;
 642     DWT1DFN dwt_1D;
 643     OPJ_UINT32 rw;
 644     OPJ_UINT32 w;
 645     OPJ_INT32 * OPJ_RESTRICT tiledp;
 646     OPJ_UINT32 min_j;
 647     OPJ_UINT32 max_j;
 648 } opj_dwd_decode_h_job_t;
 649
 650 static void opj_dwt_decode_h_func(void* user_data, opj_tls_t* tls)
 651 {
 652     OPJ_UINT32 j;
 653     opj_dwd_decode_h_job_t* job;
 654     (void)tls;
 655
 656     job = (opj_dwd_decode_h_job_t*)user_data;
 657     for (j = job->min_j; j < job->max_j; j++) {
 658         opj_dwt_interleave_h(&job->h, &job->tiledp[j * job->w]);
 659         (job->dwt_1D)(&job->h);
 660         memcpy(&job->tiledp[j * job->w], job->h.mem, job->rw * sizeof(OPJ_INT32));
 661     }
 662
 663     opj_aligned_free(job->h.mem);
 664     opj_free(job);
 665 }
 666
 667 typedef struct {
 668     opj_dwt_t v;
 669     DWT1DFN dwt_1D;
 670     OPJ_UINT32 rh;
 671     OPJ_UINT32 w;
 672     OPJ_INT32 * OPJ_RESTRICT tiledp;
 673     OPJ_UINT32 min_j;
 674     OPJ_UINT32 max_j;
 675 } opj_dwd_decode_v_job_t;
 676
 677 static void opj_dwt_decode_v_func(void* user_data, opj_tls_t* tls)
 678 {
 679     OPJ_UINT32 j;
 680     opj_dwd_decode_v_job_t* job;
 681     (void)tls;
 682
 683     job = (opj_dwd_decode_v_job_t*)user_data;
 684     for (j = job->min_j; j < job->max_j; j++) {
 685         OPJ_UINT32 k;
 686         opj_dwt_interleave_v(&job->v, &job->tiledp[j], (OPJ_INT32)job->w);
 687         (job->dwt_1D)(&job->v);
 688         for (k = 0; k < job->rh; ++k) {
 689             job->tiledp[k * job->w + j] = job->v.mem[k];
 690         }
 691     }
 692
 693     opj_aligned_free(job->v.mem);
 694     opj_free(job);
 695 }
 696
 697
 698 /* <summary>                            */
 699 /* Inverse wavelet transform in 2-D.    */
 700 /* </summary>                           */
 701 static OPJ_BOOL opj_dwt_decode_tile(opj_thread_pool_t* tp,
 702                                     opj_tcd_tilecomp_t* tilec, OPJ_UINT32 numres, DWT1DFN dwt_1D)
 703 {
 704     opj_dwt_t h;
 705     opj_dwt_t v;
 706
 707     opj_tcd_resolution_t* tr = tilec->resolutions;
 708
 709     OPJ_UINT32 rw = (OPJ_UINT32)(tr->x1 -
 710                                  tr->x0);  /* width of the resolution level computed */
 711     OPJ_UINT32 rh = (OPJ_UINT32)(tr->y1 -
 712                                  tr->y0);  /* height of the resolution level computed */
 713
 714     OPJ_UINT32 w = (OPJ_UINT32)(tilec->x1 - tilec->x0);
 715     size_t h_mem_size;
 716     int num_threads;
 717
 718     if (numres == 1U) {
 719         return OPJ_TRUE;
 720     }
 721     num_threads = opj_thread_pool_get_thread_count(tp);
 722     h_mem_size = opj_dwt_max_resolution(tr, numres);
 723     /* overflow check */
 724     if (h_mem_size > (SIZE_MAX / sizeof(OPJ_INT32))) {
 725         /* FIXME event manager error callback */
 726         return OPJ_FALSE;
 727     }
 728     h_mem_size *= sizeof(OPJ_INT32);
 729     h.mem = (OPJ_INT32*)opj_aligned_malloc(h_mem_size);
 730     if (! h.mem) {
 731         /* FIXME event manager error callback */
 732         return OPJ_FALSE;
 733     }
 734
 735     v.mem = h.mem;
 736
 737     while (--numres) {
 738         OPJ_INT32 * OPJ_RESTRICT tiledp = tilec->data;
 739         OPJ_UINT32 j;
 740
 741         ++tr;
 742         h.sn = (OPJ_INT32)rw;
 743         v.sn = (OPJ_INT32)rh;
 744
 745         rw = (OPJ_UINT32)(tr->x1 - tr->x0);
 746         rh = (OPJ_UINT32)(tr->y1 - tr->y0);
 747
 748         h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
 749         h.cas = tr->x0 % 2;
 750
 751         if (num_threads <= 1 || rh <= 1) {
 752             for (j = 0; j < rh; ++j) {
 753                 opj_dwt_interleave_h(&h, &tiledp[j * w]);
 754                 (dwt_1D)(&h);
 755                 memcpy(&tiledp[j * w], h.mem, rw * sizeof(OPJ_INT32));
 756             }
 757         } else {
 758             OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
 759             OPJ_UINT32 step_j;
 760
 761             if (rh < num_jobs) {
 762                 num_jobs = rh;
 763             }
 764             step_j = (rh / num_jobs);
 765
 766             for (j = 0; j < num_jobs; j++) {
 767                 opj_dwd_decode_h_job_t* job;
 768
 769                 job = (opj_dwd_decode_h_job_t*) opj_malloc(sizeof(opj_dwd_decode_h_job_t));
 770                 if (!job) {
 771                     /* It would be nice to fallback to single thread case, but */
 772                     /* unfortunately some jobs may be launched and have modified */
 773                     /* tiledp, so it is not practical to recover from that error */
 774                     /* FIXME event manager error callback */
 775                     opj_thread_pool_wait_completion(tp, 0);
 776                     opj_aligned_free(h.mem);
 777                     return OPJ_FALSE;
 778                 }
 779                 job->h = h;
 780                 job->dwt_1D = dwt_1D;
 781                 job->rw = rw;
 782                 job->w = w;
 783                 job->tiledp = tiledp;
 784                 job->min_j = j * step_j;
 785                 job->max_j = (j + 1U) * step_j; /* this can overflow */
 786                 if (j == (num_jobs - 1U)) {  /* this will take care of the overflow */
 787                     job->max_j = rh;
 788                 }
 789                 job->h.mem = (OPJ_INT32*)opj_aligned_malloc(h_mem_size);
 790                 if (!job->h.mem) {
 791                     /* FIXME event manager error callback */
 792                     opj_thread_pool_wait_completion(tp, 0);
 793                     opj_free(job);
 794                     opj_aligned_free(h.mem);
 795                     return OPJ_FALSE;
 796                 }
 797                 opj_thread_pool_submit_job(tp, opj_dwt_decode_h_func, job);
 798             }
 799             opj_thread_pool_wait_completion(tp, 0);
 800         }
 801
 802         v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
 803         v.cas = tr->y0 % 2;
 804
 805         if (num_threads <= 1 || rw <= 1) {
 806             for (j = 0; j < rw; ++j) {
 807                 OPJ_UINT32 k;
 808
 809                 opj_dwt_interleave_v(&v, &tiledp[j], (OPJ_INT32)w);
 810                 (dwt_1D)(&v);
 811                 for (k = 0; k < rh; ++k) {
 812                     tiledp[k * w + j] = v.mem[k];
 813                 }
 814             }
 815         } else {
 816             OPJ_UINT32 num_jobs = (OPJ_UINT32)num_threads;
 817             OPJ_UINT32 step_j;
 818
 819             if (rw < num_jobs) {
 820                 num_jobs = rw;
 821             }
 822             step_j = (rw / num_jobs);
 823
 824             for (j = 0; j < num_jobs; j++) {
 825                 opj_dwd_decode_v_job_t* job;
 826
 827                 job = (opj_dwd_decode_v_job_t*) opj_malloc(sizeof(opj_dwd_decode_v_job_t));
 828                 if (!job) {
 829                     /* It would be nice to fallback to single thread case, but */
 830                     /* unfortunately some jobs may be launched and have modified */
 831                     /* tiledp, so it is not practical to recover from that error */
 832                     /* FIXME event manager error callback */
 833                     opj_thread_pool_wait_completion(tp, 0);
 834                     opj_aligned_free(v.mem);
 835                     return OPJ_FALSE;
 836                 }
 837                 job->v = v;
 838                 job->dwt_1D = dwt_1D;
 839                 job->rh = rh;
 840                 job->w = w;
 841                 job->tiledp = tiledp;
 842                 job->min_j = j * step_j;
 843                 job->max_j = (j + 1U) * step_j; /* this can overflow */
 844                 if (j == (num_jobs - 1U)) {  /* this will take care of the overflow */
 845                     job->max_j = rw;
 846                 }
 847                 job->v.mem = (OPJ_INT32*)opj_aligned_malloc(h_mem_size);
 848                 if (!job->v.mem) {
 849                     /* FIXME event manager error callback */
 850                     opj_thread_pool_wait_completion(tp, 0);
 851                     opj_free(job);
 852                     opj_aligned_free(v.mem);
 853                     return OPJ_FALSE;
 854                 }
 855                 opj_thread_pool_submit_job(tp, opj_dwt_decode_v_func, job);
 856             }
 857             opj_thread_pool_wait_completion(tp, 0);
 858         }
 859     }
 860     opj_aligned_free(h.mem);
 861     return OPJ_TRUE;
 862 }
 863
 864 static void opj_v4dwt_interleave_h(opj_v4dwt_t* OPJ_RESTRICT w,
 865                                    OPJ_FLOAT32* OPJ_RESTRICT a, OPJ_INT32 x, OPJ_INT32 size)
 866 {
 867     OPJ_FLOAT32* OPJ_RESTRICT bi = (OPJ_FLOAT32*)(w->wavelet + w->cas);
 868     OPJ_INT32 count = w->sn;
 869     OPJ_INT32 i, k;
 870
 871     for (k = 0; k < 2; ++k) {
 872         if (count + 3 * x < size && ((size_t) a & 0x0f) == 0 &&
 873                 ((size_t) bi & 0x0f) == 0 && (x & 0x0f) == 0) {
 874             /* Fast code path */
 875             for (i = 0; i < count; ++i) {
 876                 OPJ_INT32 j = i;
 877                 bi[i * 8    ] = a[j];
 878                 j += x;
 879                 bi[i * 8 + 1] = a[j];
 880                 j += x;
 881                 bi[i * 8 + 2] = a[j];
 882                 j += x;
 883                 bi[i * 8 + 3] = a[j];
 884             }
 885         } else {
 886             /* Slow code path */
 887             for (i = 0; i < count; ++i) {
 888                 OPJ_INT32 j = i;
 889                 bi[i * 8    ] = a[j];
 890                 j += x;
 891                 if (j >= size) {
 892                     continue;
 893                 }
 894                 bi[i * 8 + 1] = a[j];
 895                 j += x;
 896                 if (j >= size) {
 897                     continue;
 898                 }
 899                 bi[i * 8 + 2] = a[j];
 900                 j += x;
 901                 if (j >= size) {
 902                     continue;
 903                 }
 904                 bi[i * 8 + 3] = a[j]; /* This one*/
 905             }
 906         }
 907
 908         bi = (OPJ_FLOAT32*)(w->wavelet + 1 - w->cas);
 909         a += w->sn;
 910         size -= w->sn;
 911         count = w->dn;
 912     }
 913 }
 914
 915 static void opj_v4dwt_interleave_v(opj_v4dwt_t* OPJ_RESTRICT v,
 916                                    OPJ_FLOAT32* OPJ_RESTRICT a, OPJ_INT32 x, OPJ_INT32 nb_elts_read)
 917 {
 918     opj_v4_t* OPJ_RESTRICT bi = v->wavelet + v->cas;
 919     OPJ_INT32 i;
 920
 921     for (i = 0; i < v->sn; ++i) {
 922         memcpy(&bi[i * 2], &a[i * x], (size_t)nb_elts_read * sizeof(OPJ_FLOAT32));
 923     }
 924
 925     a += v->sn * x;
 926     bi = v->wavelet + 1 - v->cas;
 927
 928     for (i = 0; i < v->dn; ++i) {
 929         memcpy(&bi[i * 2], &a[i * x], (size_t)nb_elts_read * sizeof(OPJ_FLOAT32));
 930     }
 931 }
 932
 933 #ifdef __SSE__
 934
 935 static void opj_v4dwt_decode_step1_sse(opj_v4_t* w, OPJ_INT32 count,
 936                                        const __m128 c)
 937 {
 938     __m128* OPJ_RESTRICT vw = (__m128*) w;
 939     OPJ_INT32 i;
 940     /* 4x unrolled loop */
 941     for (i = 0; i < count >> 2; ++i) {
 942         *vw = _mm_mul_ps(*vw, c);
 943         vw += 2;
 944         *vw = _mm_mul_ps(*vw, c);
 945         vw += 2;
 946         *vw = _mm_mul_ps(*vw, c);
 947         vw += 2;
 948         *vw = _mm_mul_ps(*vw, c);
 949         vw += 2;
 950     }
 951     count &= 3;
 952     for (i = 0; i < count; ++i) {
 953         *vw = _mm_mul_ps(*vw, c);
 954         vw += 2;
 955     }
 956 }
 957
 958 void opj_v4dwt_decode_step2_sse(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k,
 959                                 OPJ_INT32 m, __m128 c)
 960 {
 961     __m128* OPJ_RESTRICT vl = (__m128*) l;
 962     __m128* OPJ_RESTRICT vw = (__m128*) w;
 963     OPJ_INT32 i;
 964     __m128 tmp1, tmp2, tmp3;
 965     tmp1 = vl[0];
 966     for (i = 0; i < m; ++i) {
 967         tmp2 = vw[-1];
 968         tmp3 = vw[ 0];
 969         vw[-1] = _mm_add_ps(tmp2, _mm_mul_ps(_mm_add_ps(tmp1, tmp3), c));
 970         tmp1 = tmp3;
 971         vw += 2;
 972     }
 973     vl = vw - 2;
 974     if (m >= k) {
 975         return;
 976     }
 977     c = _mm_add_ps(c, c);
 978     c = _mm_mul_ps(c, vl[0]);
 979     for (; m < k; ++m) {
 980         __m128 tmp = vw[-1];
 981         vw[-1] = _mm_add_ps(tmp, c);
 982         vw += 2;
 983     }
 984 }
 985
 986 #else
 987
 988 static void opj_v4dwt_decode_step1(opj_v4_t* w, OPJ_INT32 count,
 989                                    const OPJ_FLOAT32 c)
 990 {
 991     OPJ_FLOAT32* OPJ_RESTRICT fw = (OPJ_FLOAT32*) w;
 992     OPJ_INT32 i;
 993     for (i = 0; i < count; ++i) {
 994         OPJ_FLOAT32 tmp1 = fw[i * 8    ];
 995         OPJ_FLOAT32 tmp2 = fw[i * 8 + 1];
 996         OPJ_FLOAT32 tmp3 = fw[i * 8 + 2];
 997         OPJ_FLOAT32 tmp4 = fw[i * 8 + 3];
 998         fw[i * 8    ] = tmp1 * c;
 999         fw[i * 8 + 1] = tmp2 * c;
1000         fw[i * 8 + 2] = tmp3 * c;
1001         fw[i * 8 + 3] = tmp4 * c;
1002     }
1003 }
1004
1005 static void opj_v4dwt_decode_step2(opj_v4_t* l, opj_v4_t* w, OPJ_INT32 k,
1006                                    OPJ_INT32 m, OPJ_FLOAT32 c)
1007 {
1008     OPJ_FLOAT32* fl = (OPJ_FLOAT32*) l;
1009     OPJ_FLOAT32* fw = (OPJ_FLOAT32*) w;
1010     OPJ_INT32 i;
1011     for (i = 0; i < m; ++i) {
1012         OPJ_FLOAT32 tmp1_1 = fl[0];
1013         OPJ_FLOAT32 tmp1_2 = fl[1];
1014         OPJ_FLOAT32 tmp1_3 = fl[2];
1015         OPJ_FLOAT32 tmp1_4 = fl[3];
1016         OPJ_FLOAT32 tmp2_1 = fw[-4];
1017         OPJ_FLOAT32 tmp2_2 = fw[-3];
1018         OPJ_FLOAT32 tmp2_3 = fw[-2];
1019         OPJ_FLOAT32 tmp2_4 = fw[-1];
1020         OPJ_FLOAT32 tmp3_1 = fw[0];
1021         OPJ_FLOAT32 tmp3_2 = fw[1];
1022         OPJ_FLOAT32 tmp3_3 = fw[2];
1023         OPJ_FLOAT32 tmp3_4 = fw[3];
1024         fw[-4] = tmp2_1 + ((tmp1_1 + tmp3_1) * c);
1025         fw[-3] = tmp2_2 + ((tmp1_2 + tmp3_2) * c);
1026         fw[-2] = tmp2_3 + ((tmp1_3 + tmp3_3) * c);
1027         fw[-1] = tmp2_4 + ((tmp1_4 + tmp3_4) * c);
1028         fl = fw;
1029         fw += 8;
1030     }
1031     if (m < k) {
1032         OPJ_FLOAT32 c1;
1033         OPJ_FLOAT32 c2;
1034         OPJ_FLOAT32 c3;
1035         OPJ_FLOAT32 c4;
1036         c += c;
1037         c1 = fl[0] * c;
1038         c2 = fl[1] * c;
1039         c3 = fl[2] * c;
1040         c4 = fl[3] * c;
1041         for (; m < k; ++m) {
1042             OPJ_FLOAT32 tmp1 = fw[-4];
1043             OPJ_FLOAT32 tmp2 = fw[-3];
1044             OPJ_FLOAT32 tmp3 = fw[-2];
1045             OPJ_FLOAT32 tmp4 = fw[-1];
1046             fw[-4] = tmp1 + c1;
1047             fw[-3] = tmp2 + c2;
1048             fw[-2] = tmp3 + c3;
1049             fw[-1] = tmp4 + c4;
1050             fw += 8;
1051         }
1052     }
1053 }
1054
1055 #endif
1056
1057 /* <summary>                             */
1058 /* Inverse 9-7 wavelet transform in 1-D. */
1059 /* </summary>                            */
1060 static void opj_v4dwt_decode(opj_v4dwt_t* OPJ_RESTRICT dwt)
1061 {
1062     OPJ_INT32 a, b;
1063     if (dwt->cas == 0) {
1064         if (!((dwt->dn > 0) || (dwt->sn > 1))) {
1065             return;
1066         }
1067         a = 0;
1068         b = 1;
1069     } else {
1070         if (!((dwt->sn > 0) || (dwt->dn > 1))) {
1071             return;
1072         }
1073         a = 1;
1074         b = 0;
1075     }
1076 #ifdef __SSE__
1077     opj_v4dwt_decode_step1_sse(dwt->wavelet + a, dwt->sn, _mm_set1_ps(opj_K));
1078     opj_v4dwt_decode_step1_sse(dwt->wavelet + b, dwt->dn, _mm_set1_ps(opj_c13318));
1079     opj_v4dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
1080                                opj_int_min(dwt->sn, dwt->dn - a), _mm_set1_ps(opj_dwt_delta));
1081     opj_v4dwt_decode_step2_sse(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
1082                                opj_int_min(dwt->dn, dwt->sn - b), _mm_set1_ps(opj_dwt_gamma));
1083     opj_v4dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
1084                                opj_int_min(dwt->sn, dwt->dn - a), _mm_set1_ps(opj_dwt_beta));
1085     opj_v4dwt_decode_step2_sse(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
1086                                opj_int_min(dwt->dn, dwt->sn - b), _mm_set1_ps(opj_dwt_alpha));
1087 #else
1088     opj_v4dwt_decode_step1(dwt->wavelet + a, dwt->sn, opj_K);
1089     opj_v4dwt_decode_step1(dwt->wavelet + b, dwt->dn, opj_c13318);
1090     opj_v4dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
1091                            opj_int_min(dwt->sn, dwt->dn - a), opj_dwt_delta);
1092     opj_v4dwt_decode_step2(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
1093                            opj_int_min(dwt->dn, dwt->sn - b), opj_dwt_gamma);
1094     opj_v4dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
1095                            opj_int_min(dwt->sn, dwt->dn - a), opj_dwt_beta);
1096     opj_v4dwt_decode_step2(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
1097                            opj_int_min(dwt->dn, dwt->sn - b), opj_dwt_alpha);
1098 #endif
1099 }
1100
1101
1102 /* <summary>                             */
1103 /* Inverse 9-7 wavelet transform in 2-D. */
1104 /* </summary>                            */
1105 OPJ_BOOL opj_dwt_decode_real(opj_tcd_tilecomp_t* OPJ_RESTRICT tilec,
1106                              OPJ_UINT32 numres)
1107 {
1108     opj_v4dwt_t h;
1109     opj_v4dwt_t v;
1110
1111     opj_tcd_resolution_t* res = tilec->resolutions;
1112
1113     OPJ_UINT32 rw = (OPJ_UINT32)(res->x1 -
1114                                  res->x0);    /* width of the resolution level computed */
1115     OPJ_UINT32 rh = (OPJ_UINT32)(res->y1 -
1116                                  res->y0);    /* height of the resolution level computed */
1117
1118     OPJ_UINT32 w = (OPJ_UINT32)(tilec->x1 - tilec->x0);
1119
1120     size_t l_data_size;
1121
1122     l_data_size = opj_dwt_max_resolution(res, numres);
1123     /* overflow check */
1124     if (l_data_size > (SIZE_MAX - 5U)) {
1125         /* FIXME event manager error callback */
1126         return OPJ_FALSE;
1127     }
1128     l_data_size += 5U;
1129     /* overflow check */
1130     if (l_data_size > (SIZE_MAX / sizeof(opj_v4_t))) {
1131         /* FIXME event manager error callback */
1132         return OPJ_FALSE;
1133     }
1134     h.wavelet = (opj_v4_t*) opj_aligned_malloc(l_data_size * sizeof(opj_v4_t));
1135     if (!h.wavelet) {
1136         /* FIXME event manager error callback */
1137         return OPJ_FALSE;
1138     }
1139     v.wavelet = h.wavelet;
1140
1141     while (--numres) {
1142         OPJ_FLOAT32 * OPJ_RESTRICT aj = (OPJ_FLOAT32*) tilec->data;
1143         OPJ_UINT32 bufsize = (OPJ_UINT32)((tilec->x1 - tilec->x0) *
1144                                           (tilec->y1 - tilec->y0));
1145         OPJ_INT32 j;
1146
1147         h.sn = (OPJ_INT32)rw;
1148         v.sn = (OPJ_INT32)rh;
1149
1150         ++res;
1151
1152         rw = (OPJ_UINT32)(res->x1 -
1153                           res->x0);   /* width of the resolution level computed */
1154         rh = (OPJ_UINT32)(res->y1 -
1155                           res->y0);   /* height of the resolution level computed */
1156
1157         h.dn = (OPJ_INT32)(rw - (OPJ_UINT32)h.sn);
1158         h.cas = res->x0 % 2;
1159
1160         for (j = (OPJ_INT32)rh; j > 3; j -= 4) {
1161             OPJ_INT32 k;
1162             opj_v4dwt_interleave_h(&h, aj, (OPJ_INT32)w, (OPJ_INT32)bufsize);
1163             opj_v4dwt_decode(&h);
1164
1165             for (k = (OPJ_INT32)rw; --k >= 0;) {
1166                 aj[k               ] = h.wavelet[k].f[0];
1167                 aj[k + (OPJ_INT32)w  ] = h.wavelet[k].f[1];
1168                 aj[k + (OPJ_INT32)w * 2] = h.wavelet[k].f[2];
1169                 aj[k + (OPJ_INT32)w * 3] = h.wavelet[k].f[3];
1170             }
1171
1172             aj += w * 4;
1173             bufsize -= w * 4;
1174         }
1175
1176         if (rh & 0x03) {
1177             OPJ_INT32 k;
1178             j = rh & 0x03;
1179             opj_v4dwt_interleave_h(&h, aj, (OPJ_INT32)w, (OPJ_INT32)bufsize);
1180             opj_v4dwt_decode(&h);
1181             for (k = (OPJ_INT32)rw; --k >= 0;) {
1182                 switch (j) {
1183                 case 3:
1184                     aj[k + (OPJ_INT32)w * 2] = h.wavelet[k].f[2];
1185                 /* FALLTHRU */
1186                 case 2:
1187                     aj[k + (OPJ_INT32)w  ] = h.wavelet[k].f[1];
1188                 /* FALLTHRU */
1189                 case 1:
1190                     aj[k               ] = h.wavelet[k].f[0];
1191                 }
1192             }
1193         }
1194
1195         v.dn = (OPJ_INT32)(rh - (OPJ_UINT32)v.sn);
1196         v.cas = res->y0 % 2;
1197
1198         aj = (OPJ_FLOAT32*) tilec->data;
1199         for (j = (OPJ_INT32)rw; j > 3; j -= 4) {
1200             OPJ_UINT32 k;
1201
1202             opj_v4dwt_interleave_v(&v, aj, (OPJ_INT32)w, 4);
1203             opj_v4dwt_decode(&v);
1204
1205             for (k = 0; k < rh; ++k) {
1206                 memcpy(&aj[k * w], &v.wavelet[k], 4 * sizeof(OPJ_FLOAT32));
1207             }
1208             aj += 4;
1209         }
1210
1211         if (rw & 0x03) {
1212             OPJ_UINT32 k;
1213
1214             j = rw & 0x03;
1215
1216             opj_v4dwt_interleave_v(&v, aj, (OPJ_INT32)w, j);
1217             opj_v4dwt_decode(&v);
1218
1219             for (k = 0; k < rh; ++k) {
1220                 memcpy(&aj[k * w], &v.wavelet[k], (size_t)j * sizeof(OPJ_FLOAT32));
1221             }
1222         }
1223     }
1224
1225     opj_aligned_free(h.wavelet);
1226     return OPJ_TRUE;
1227 }