libopenjpeg/dwt.c

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