src/lib/openmj2/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 WS(i) v->mem[(i)*2]
  50 #define WD(i) v->mem[(1+(i)*2)]
  51
  52 /** @name Local data structures */
  53 /*@{*/
  54
  55 typedef struct dwt_local {
  56     int* mem;
  57     int dn;
  58     int sn;
  59     int cas;
  60 } dwt_t;
  61
  62 typedef union {
  63     float   f[4];
  64 } v4;
  65
  66 typedef struct v4dwt_local {
  67     v4* wavelet ;
  68     int     dn ;
  69     int     sn ;
  70     int     cas ;
  71 } v4dwt_t ;
  72
  73 static const float dwt_alpha =  1.586134342f; /*  12994 */
  74 static const float dwt_beta  =  0.052980118f; /*    434 */
  75 static const float dwt_gamma = -0.882911075f; /*  -7233 */
  76 static const float dwt_delta = -0.443506852f; /*  -3633 */
  77
  78 static const float K      = 1.230174105f; /*  10078 */
  79 /* FIXME: What is this constant? */
  80 static const float c13318 = 1.625732422f;
  81
  82 /*@}*/
  83
  84 /**
  85 Virtual function type for wavelet transform in 1-D
  86 */
  87 typedef void (*DWT1DFN)(dwt_t* v);
  88
  89 /** @name Local static functions */
  90 /*@{*/
  91
  92 /**
  93 Forward lazy transform (horizontal)
  94 */
  95 static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas);
  96 /**
  97 Forward lazy transform (vertical)
  98 */
  99 static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas);
 100 /**
 101 Inverse lazy transform (horizontal)
 102 */
 103 static void dwt_interleave_h(dwt_t* h, int *a);
 104 /**
 105 Inverse lazy transform (vertical)
 106 */
 107 static void dwt_interleave_v(dwt_t* v, int *a, int x);
 108 /**
 109 Forward 5-3 wavelet transform in 1-D
 110 */
 111 static void dwt_encode_1(int *a, int dn, int sn, int cas);
 112 /**
 113 Inverse 5-3 wavelet transform in 1-D
 114 */
 115 static void dwt_decode_1(dwt_t *v);
 116 /**
 117 Forward 9-7 wavelet transform in 1-D
 118 */
 119 static void dwt_encode_1_real(int *a, int dn, int sn, int cas);
 120 /**
 121 Explicit calculation of the Quantization Stepsizes
 122 */
 123 static void dwt_encode_stepsize(int stepsize, int numbps,
 124                                 opj_stepsize_t *bandno_stepsize);
 125 /**
 126 Inverse wavelet transform in 2-D.
 127 */
 128 static void dwt_decode_tile(opj_tcd_tilecomp_t* tilec, int i, DWT1DFN fn);
 129
 130 /*@}*/
 131
 132 /*@}*/
 133
 134 #define S(i) a[(i)*2]
 135 #define D(i) a[(1+(i)*2)]
 136 #define S_(i) ((i)<0?S(0):((i)>=sn?S(sn-1):S(i)))
 137 #define D_(i) ((i)<0?D(0):((i)>=dn?D(dn-1):D(i)))
 138 /* new */
 139 #define SS_(i) ((i)<0?S(0):((i)>=dn?S(dn-1):S(i)))
 140 #define DD_(i) ((i)<0?D(0):((i)>=sn?D(sn-1):D(i)))
 141
 142 /* <summary>                                                              */
 143 /* This table contains the norms of the 5-3 wavelets for different bands. */
 144 /* </summary>                                                             */
 145 static const double dwt_norms[4][10] = {
 146     {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 147     {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 148     {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 149     {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
 150 };
 151
 152 /* <summary>                                                              */
 153 /* This table contains the norms of the 9-7 wavelets for different bands. */
 154 /* </summary>                                                             */
 155 static const double dwt_norms_real[4][10] = {
 156     {1.000, 1.965, 4.177, 8.403, 16.90, 33.84, 67.69, 135.3, 270.6, 540.9},
 157     {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
 158     {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
 159     {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2}
 160 };
 161
 162 /*
 163 ==========================================================
 164    local functions
 165 ==========================================================
 166 */
 167
 168 /* <summary>                             */
 169 /* Forward lazy transform (horizontal).  */
 170 /* </summary>                            */
 171 static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas)
 172 {
 173     int i;
 174     for (i = 0; i < sn; i++) {
 175         b[i] = a[2 * i + cas];
 176     }
 177     for (i = 0; i < dn; i++) {
 178         b[sn + i] = a[(2 * i + 1 - cas)];
 179     }
 180 }
 181
 182 /* <summary>                             */
 183 /* Forward lazy transform (vertical).    */
 184 /* </summary>                            */
 185 static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas)
 186 {
 187     int i;
 188     for (i = 0; i < sn; i++) {
 189         b[i * x] = a[2 * i + cas];
 190     }
 191     for (i = 0; i < dn; i++) {
 192         b[(sn + i)*x] = a[(2 * i + 1 - cas)];
 193     }
 194 }
 195
 196 /* <summary>                             */
 197 /* Inverse lazy transform (horizontal).  */
 198 /* </summary>                            */
 199 static void dwt_interleave_h(dwt_t* h, int *a)
 200 {
 201     int *ai = a;
 202     int *bi = h->mem + h->cas;
 203     int  i  = h->sn;
 204     while (i--) {
 205         *bi = *(ai++);
 206         bi += 2;
 207     }
 208     ai  = a + h->sn;
 209     bi  = h->mem + 1 - h->cas;
 210     i   = h->dn ;
 211     while (i--) {
 212         *bi = *(ai++);
 213         bi += 2;
 214     }
 215 }
 216
 217 /* <summary>                             */
 218 /* Inverse lazy transform (vertical).    */
 219 /* </summary>                            */
 220 static void dwt_interleave_v(dwt_t* v, int *a, int x)
 221 {
 222     int *ai = a;
 223     int *bi = v->mem + v->cas;
 224     int  i = v->sn;
 225     while (i--) {
 226         *bi = *ai;
 227         bi += 2;
 228         ai += x;
 229     }
 230     ai = a + (v->sn * x);
 231     bi = v->mem + 1 - v->cas;
 232     i = v->dn ;
 233     while (i--) {
 234         *bi = *ai;
 235         bi += 2;
 236         ai += x;
 237     }
 238 }
 239
 240
 241 /* <summary>                            */
 242 /* Forward 5-3 wavelet transform in 1-D. */
 243 /* </summary>                           */
 244 static void dwt_encode_1(int *a, int dn, int sn, int cas)
 245 {
 246     int i;
 247
 248     if (!cas) {
 249         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 250             for (i = 0; i < dn; i++) {
 251                 D(i) -= (S_(i) + S_(i + 1)) >> 1;
 252             }
 253             for (i = 0; i < sn; i++) {
 254                 S(i) += (D_(i - 1) + D_(i) + 2) >> 2;
 255             }
 256         }
 257     } else {
 258         if (!sn && dn == 1) {       /* NEW :  CASE ONE ELEMENT */
 259             S(0) *= 2;
 260         } else {
 261             for (i = 0; i < dn; i++) {
 262                 S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
 263             }
 264             for (i = 0; i < sn; i++) {
 265                 D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
 266             }
 267         }
 268     }
 269 }
 270
 271 /* <summary>                            */
 272 /* Inverse 5-3 wavelet transform in 1-D. */
 273 /* </summary>                           */
 274 static void dwt_decode_1_(int *a, int dn, int sn, int cas)
 275 {
 276     int i;
 277
 278     if (!cas) {
 279         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 280             for (i = 0; i < sn; i++) {
 281                 S(i) -= (D_(i - 1) + D_(i) + 2) >> 2;
 282             }
 283             for (i = 0; i < dn; i++) {
 284                 D(i) += (S_(i) + S_(i + 1)) >> 1;
 285             }
 286         }
 287     } else {
 288         if (!sn  && dn == 1) {        /* NEW :  CASE ONE ELEMENT */
 289             S(0) /= 2;
 290         } else {
 291             for (i = 0; i < sn; i++) {
 292                 D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2;
 293             }
 294             for (i = 0; i < dn; i++) {
 295                 S(i) += (DD_(i) + DD_(i - 1)) >> 1;
 296             }
 297         }
 298     }
 299 }
 300
 301 /* <summary>                            */
 302 /* Inverse 5-3 wavelet transform in 1-D. */
 303 /* </summary>                           */
 304 static void dwt_decode_1(dwt_t *v)
 305 {
 306     dwt_decode_1_(v->mem, v->dn, v->sn, v->cas);
 307 }
 308
 309 /* <summary>                             */
 310 /* Forward 9-7 wavelet transform in 1-D. */
 311 /* </summary>                            */
 312 static void dwt_encode_1_real(int *a, int dn, int sn, int cas)
 313 {
 314     int i;
 315     if (!cas) {
 316         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 317             for (i = 0; i < dn; i++) {
 318                 D(i) -= fix_mul(S_(i) + S_(i + 1), 12993);
 319             }
 320             for (i = 0; i < sn; i++) {
 321                 S(i) -= fix_mul(D_(i - 1) + D_(i), 434);
 322             }
 323             for (i = 0; i < dn; i++) {
 324                 D(i) += fix_mul(S_(i) + S_(i + 1), 7233);
 325             }
 326             for (i = 0; i < sn; i++) {
 327                 S(i) += fix_mul(D_(i - 1) + D_(i), 3633);
 328             }
 329             for (i = 0; i < dn; i++) {
 330                 D(i) = fix_mul(D(i), 5038);    /*5038 */
 331             }
 332             for (i = 0; i < sn; i++) {
 333                 S(i) = fix_mul(S(i), 6659);    /*6660 */
 334             }
 335         }
 336     } else {
 337         if ((sn > 0) || (dn > 1)) { /* NEW :  CASE ONE ELEMENT */
 338             for (i = 0; i < dn; i++) {
 339                 S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993);
 340             }
 341             for (i = 0; i < sn; i++) {
 342                 D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434);
 343             }
 344             for (i = 0; i < dn; i++) {
 345                 S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233);
 346             }
 347             for (i = 0; i < sn; i++) {
 348                 D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633);
 349             }
 350             for (i = 0; i < dn; i++) {
 351                 S(i) = fix_mul(S(i), 5038);    /*5038 */
 352             }
 353             for (i = 0; i < sn; i++) {
 354                 D(i) = fix_mul(D(i), 6659);    /*6660 */
 355             }
 356         }
 357     }
 358 }
 359
 360 static void dwt_encode_stepsize(int stepsize, int numbps,
 361                                 opj_stepsize_t *bandno_stepsize)
 362 {
 363     int p, n;
 364     p = int_floorlog2(stepsize) - 13;
 365     n = 11 - int_floorlog2(stepsize);
 366     bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
 367     bandno_stepsize->expn = numbps - p;
 368 }
 369
 370 /*
 371 ==========================================================
 372    DWT interface
 373 ==========================================================
 374 */
 375
 376 /* <summary>                            */
 377 /* Forward 5-3 wavelet transform in 2-D. */
 378 /* </summary>                           */
 379 void dwt_encode(opj_tcd_tilecomp_t * tilec)
 380 {
 381     int i, j, k;
 382     int *a = NULL;
 383     int *aj = NULL;
 384     int *bj = NULL;
 385     int w, l;
 386
 387     w = tilec->x1 - tilec->x0;
 388     l = tilec->numresolutions - 1;
 389     a = tilec->data;
 390
 391     for (i = 0; i < l; i++) {
 392         int rw;         /* width of the resolution level computed                                                           */
 393         int rh;         /* height of the resolution level computed                                                          */
 394         int rw1;        /* width of the resolution level once lower than computed one                                       */
 395         int rh1;        /* height of the resolution level once lower than computed one                                      */
 396         int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 397         int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 398         int dn, sn;
 399
 400         rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 401         rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 402         rw1 = tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 403         rh1 = tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 404
 405         cas_row = tilec->resolutions[l - i].x0 % 2;
 406         cas_col = tilec->resolutions[l - i].y0 % 2;
 407
 408         sn = rh1;
 409         dn = rh - rh1;
 410         bj = (int*)opj_malloc(rh * sizeof(int));
 411         for (j = 0; j < rw; j++) {
 412             aj = a + j;
 413             for (k = 0; k < rh; k++) {
 414                 bj[k] = aj[k * w];
 415             }
 416             dwt_encode_1(bj, dn, sn, cas_col);
 417             dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 418         }
 419         opj_free(bj);
 420
 421         sn = rw1;
 422         dn = rw - rw1;
 423         bj = (int*)opj_malloc(rw * sizeof(int));
 424         for (j = 0; j < rh; j++) {
 425             aj = a + j * w;
 426             for (k = 0; k < rw; k++) {
 427                 bj[k] = aj[k];
 428             }
 429             dwt_encode_1(bj, dn, sn, cas_row);
 430             dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 431         }
 432         opj_free(bj);
 433     }
 434 }
 435
 436
 437 /* <summary>                            */
 438 /* Inverse 5-3 wavelet transform in 2-D. */
 439 /* </summary>                           */
 440 void dwt_decode(opj_tcd_tilecomp_t* tilec, int numres)
 441 {
 442     dwt_decode_tile(tilec, numres, &dwt_decode_1);
 443 }
 444
 445
 446 /* <summary>                          */
 447 /* Get gain of 5-3 wavelet transform. */
 448 /* </summary>                         */
 449 int dwt_getgain(int orient)
 450 {
 451     if (orient == 0) {
 452         return 0;
 453     }
 454     if (orient == 1 || orient == 2) {
 455         return 1;
 456     }
 457     return 2;
 458 }
 459
 460 /* <summary>                */
 461 /* Get norm of 5-3 wavelet. */
 462 /* </summary>               */
 463 double dwt_getnorm(int level, int orient)
 464 {
 465     return dwt_norms[orient][level];
 466 }
 467
 468 /* <summary>                             */
 469 /* Forward 9-7 wavelet transform in 2-D. */
 470 /* </summary>                            */
 471
 472 void dwt_encode_real(opj_tcd_tilecomp_t * tilec)
 473 {
 474     int i, j, k;
 475     int *a = NULL;
 476     int *aj = NULL;
 477     int *bj = NULL;
 478     int w, l;
 479
 480     w = tilec->x1 - tilec->x0;
 481     l = tilec->numresolutions - 1;
 482     a = tilec->data;
 483
 484     for (i = 0; i < l; i++) {
 485         int rw;         /* width of the resolution level computed                                                     */
 486         int rh;         /* height of the resolution level computed                                                    */
 487         int rw1;        /* width of the resolution level once lower than computed one                                 */
 488         int rh1;        /* height of the resolution level once lower than computed one                                */
 489         int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 490         int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 491         int dn, sn;
 492
 493         rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 494         rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 495         rw1 = tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 496         rh1 = tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 497
 498         cas_row = tilec->resolutions[l - i].x0 % 2;
 499         cas_col = tilec->resolutions[l - i].y0 % 2;
 500
 501         sn = rh1;
 502         dn = rh - rh1;
 503         bj = (int*)opj_malloc(rh * sizeof(int));
 504         for (j = 0; j < rw; j++) {
 505             aj = a + j;
 506             for (k = 0; k < rh; k++) {
 507                 bj[k] = aj[k * w];
 508             }
 509             dwt_encode_1_real(bj, dn, sn, cas_col);
 510             dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 511         }
 512         opj_free(bj);
 513
 514         sn = rw1;
 515         dn = rw - rw1;
 516         bj = (int*)opj_malloc(rw * sizeof(int));
 517         for (j = 0; j < rh; j++) {
 518             aj = a + j * w;
 519             for (k = 0; k < rw; k++) {
 520                 bj[k] = aj[k];
 521             }
 522             dwt_encode_1_real(bj, dn, sn, cas_row);
 523             dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 524         }
 525         opj_free(bj);
 526     }
 527 }
 528
 529
 530 /* <summary>                          */
 531 /* Get gain of 9-7 wavelet transform. */
 532 /* </summary>                         */
 533 int dwt_getgain_real(int orient)
 534 {
 535     (void)orient;
 536     return 0;
 537 }
 538
 539 /* <summary>                */
 540 /* Get norm of 9-7 wavelet. */
 541 /* </summary>               */
 542 double dwt_getnorm_real(int level, int orient)
 543 {
 544     return dwt_norms_real[orient][level];
 545 }
 546
 547 void dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, int prec)
 548 {
 549     int numbands, bandno;
 550     numbands = 3 * tccp->numresolutions - 2;
 551     for (bandno = 0; bandno < numbands; bandno++) {
 552         double stepsize;
 553         int resno, level, orient, gain;
 554
 555         resno = (bandno == 0) ? 0 : ((bandno - 1) / 3 + 1);
 556         orient = (bandno == 0) ? 0 : ((bandno - 1) % 3 + 1);
 557         level = tccp->numresolutions - 1 - resno;
 558         gain = (tccp->qmfbid == 0) ? 0 : ((orient == 0) ? 0 : (((orient == 1) ||
 559                                           (orient == 2)) ? 1 : 2));
 560         if (tccp->qntsty == J2K_CCP_QNTSTY_NOQNT) {
 561             stepsize = 1.0;
 562         } else {
 563             double norm = dwt_norms_real[orient][level];
 564             stepsize = (1 << (gain)) / norm;
 565         }
 566         dwt_encode_stepsize((int) floor(stepsize * 8192.0), prec + gain,
 567                             &tccp->stepsizes[bandno]);
 568     }
 569 }
 570
 571
 572 /* <summary>                             */
 573 /* Determine maximum computed resolution level for inverse wavelet transform */
 574 /* </summary>                            */
 575 static int dwt_decode_max_resolution(opj_tcd_resolution_t* restrict r, int i)
 576 {
 577     int mr  = 1;
 578     int w;
 579     while (--i) {
 580         r++;
 581         if (mr < (w = r->x1 - r->x0)) {
 582             mr = w ;
 583         }
 584         if (mr < (w = r->y1 - r->y0)) {
 585             mr = w ;
 586         }
 587     }
 588     return mr ;
 589 }
 590
 591
 592 /* <summary>                            */
 593 /* Inverse wavelet transform in 2-D.     */
 594 /* </summary>                           */
 595 static void dwt_decode_tile(opj_tcd_tilecomp_t* tilec, int numres,
 596                             DWT1DFN dwt_1D)
 597 {
 598     dwt_t h;
 599     dwt_t v;
 600
 601     opj_tcd_resolution_t* tr = tilec->resolutions;
 602
 603     int rw = tr->x1 - tr->x0;   /* width of the resolution level computed */
 604     int rh = tr->y1 - tr->y0;   /* height of the resolution level computed */
 605
 606     int w = tilec->x1 - tilec->x0;
 607
 608     h.mem = (int*)opj_aligned_malloc(dwt_decode_max_resolution(tr,
 609                                      numres) * sizeof(int));
 610     v.mem = h.mem;
 611
 612     while (--numres) {
 613         int * restrict tiledp = tilec->data;
 614         int j;
 615
 616         ++tr;
 617         h.sn = rw;
 618         v.sn = rh;
 619
 620         rw = tr->x1 - tr->x0;
 621         rh = tr->y1 - tr->y0;
 622
 623         h.dn = rw - h.sn;
 624         h.cas = tr->x0 % 2;
 625
 626         for (j = 0; j < rh; ++j) {
 627             dwt_interleave_h(&h, &tiledp[j * w]);
 628             (dwt_1D)(&h);
 629             memcpy(&tiledp[j * w], h.mem, rw * sizeof(int));
 630         }
 631
 632         v.dn = rh - v.sn;
 633         v.cas = tr->y0 % 2;
 634
 635         for (j = 0; j < rw; ++j) {
 636             int k;
 637             dwt_interleave_v(&v, &tiledp[j], w);
 638             (dwt_1D)(&v);
 639             for (k = 0; k < rh; ++k) {
 640                 tiledp[k * w + j] = v.mem[k];
 641             }
 642         }
 643     }
 644     opj_aligned_free(h.mem);
 645 }
 646
 647 static void v4dwt_interleave_h(v4dwt_t* restrict w, float* restrict a, int x,
 648                                int size)
 649 {
 650     float* restrict bi = (float*)(w->wavelet + w->cas);
 651     int count = w->sn;
 652     int i, k;
 653     for (k = 0; k < 2; ++k) {
 654         if (count + 3 * x < size && ((size_t) a & 0x0f) == 0 &&
 655                 ((size_t) bi & 0x0f) == 0 && (x & 0x0f) == 0) {
 656             /* Fast code path */
 657             for (i = 0; i < count; ++i) {
 658                 int j = i;
 659                 bi[i * 8    ] = a[j];
 660                 j += x;
 661                 bi[i * 8 + 1] = a[j];
 662                 j += x;
 663                 bi[i * 8 + 2] = a[j];
 664                 j += x;
 665                 bi[i * 8 + 3] = a[j];
 666             }
 667         } else {
 668             /* Slow code path */
 669             for (i = 0; i < count; ++i) {
 670                 int j = i;
 671                 bi[i * 8    ] = a[j];
 672                 j += x;
 673                 if (j > size) {
 674                     continue;
 675                 }
 676                 bi[i * 8 + 1] = a[j];
 677                 j += x;
 678                 if (j > size) {
 679                     continue;
 680                 }
 681                 bi[i * 8 + 2] = a[j];
 682                 j += x;
 683                 if (j > size) {
 684                     continue;
 685                 }
 686                 bi[i * 8 + 3] = a[j];
 687             }
 688         }
 689         bi = (float*)(w->wavelet + 1 - w->cas);
 690         a += w->sn;
 691         size -= w->sn;
 692         count = w->dn;
 693     }
 694 }
 695
 696 static void v4dwt_interleave_v(v4dwt_t* restrict v, float* restrict a, int x)
 697 {
 698     v4* restrict bi = v->wavelet + v->cas;
 699     int i;
 700     for (i = 0; i < v->sn; ++i) {
 701         memcpy(&bi[i * 2], &a[i * x], 4 * sizeof(float));
 702     }
 703     a += v->sn * x;
 704     bi = v->wavelet + 1 - v->cas;
 705     for (i = 0; i < v->dn; ++i) {
 706         memcpy(&bi[i * 2], &a[i * x], 4 * sizeof(float));
 707     }
 708 }
 709
 710 #ifdef __SSE__
 711
 712 static void v4dwt_decode_step1_sse(v4* w, int count, const __m128 c)
 713 {
 714     __m128* restrict vw = (__m128*) w;
 715     int i;
 716     /* 4x unrolled loop */
 717     for (i = 0; i < count >> 2; ++i) {
 718         *vw = _mm_mul_ps(*vw, c);
 719         vw += 2;
 720         *vw = _mm_mul_ps(*vw, c);
 721         vw += 2;
 722         *vw = _mm_mul_ps(*vw, c);
 723         vw += 2;
 724         *vw = _mm_mul_ps(*vw, c);
 725         vw += 2;
 726     }
 727     count &= 3;
 728     for (i = 0; i < count; ++i) {
 729         *vw = _mm_mul_ps(*vw, c);
 730         vw += 2;
 731     }
 732 }
 733
 734 static void v4dwt_decode_step2_sse(v4* l, v4* w, int k, int m, __m128 c)
 735 {
 736     __m128* restrict vl = (__m128*) l;
 737     __m128* restrict vw = (__m128*) w;
 738     int i;
 739     __m128 tmp1, tmp2, tmp3;
 740     tmp1 = vl[0];
 741     for (i = 0; i < m; ++i) {
 742         tmp2 = vw[-1];
 743         tmp3 = vw[ 0];
 744         vw[-1] = _mm_add_ps(tmp2, _mm_mul_ps(_mm_add_ps(tmp1, tmp3), c));
 745         tmp1 = tmp3;
 746         vw += 2;
 747     }
 748     vl = vw - 2;
 749     if (m >= k) {
 750         return;
 751     }
 752     c = _mm_add_ps(c, c);
 753     c = _mm_mul_ps(c, vl[0]);
 754     for (; m < k; ++m) {
 755         __m128 tmp = vw[-1];
 756         vw[-1] = _mm_add_ps(tmp, c);
 757         vw += 2;
 758     }
 759 }
 760
 761 #else
 762
 763 static void v4dwt_decode_step1(v4* w, int count, const float c)
 764 {
 765     float* restrict fw = (float*) w;
 766     int i;
 767     for (i = 0; i < count; ++i) {
 768         float tmp1 = fw[i * 8    ];
 769         float tmp2 = fw[i * 8 + 1];
 770         float tmp3 = fw[i * 8 + 2];
 771         float tmp4 = fw[i * 8 + 3];
 772         fw[i * 8    ] = tmp1 * c;
 773         fw[i * 8 + 1] = tmp2 * c;
 774         fw[i * 8 + 2] = tmp3 * c;
 775         fw[i * 8 + 3] = tmp4 * c;
 776     }
 777 }
 778
 779 static void v4dwt_decode_step2(v4* l, v4* w, int k, int m, float c)
 780 {
 781     float* restrict fl = (float*) l;
 782     float* restrict fw = (float*) w;
 783     int i;
 784     for (i = 0; i < m; ++i) {
 785         float tmp1_1 = fl[0];
 786         float tmp1_2 = fl[1];
 787         float tmp1_3 = fl[2];
 788         float tmp1_4 = fl[3];
 789         float tmp2_1 = fw[-4];
 790         float tmp2_2 = fw[-3];
 791         float tmp2_3 = fw[-2];
 792         float tmp2_4 = fw[-1];
 793         float tmp3_1 = fw[0];
 794         float tmp3_2 = fw[1];
 795         float tmp3_3 = fw[2];
 796         float tmp3_4 = fw[3];
 797         fw[-4] = tmp2_1 + ((tmp1_1 + tmp3_1) * c);
 798         fw[-3] = tmp2_2 + ((tmp1_2 + tmp3_2) * c);
 799         fw[-2] = tmp2_3 + ((tmp1_3 + tmp3_3) * c);
 800         fw[-1] = tmp2_4 + ((tmp1_4 + tmp3_4) * c);
 801         fl = fw;
 802         fw += 8;
 803     }
 804     if (m < k) {
 805         float c1;
 806         float c2;
 807         float c3;
 808         float c4;
 809         c += c;
 810         c1 = fl[0] * c;
 811         c2 = fl[1] * c;
 812         c3 = fl[2] * c;
 813         c4 = fl[3] * c;
 814         for (; m < k; ++m) {
 815             float tmp1 = fw[-4];
 816             float tmp2 = fw[-3];
 817             float tmp3 = fw[-2];
 818             float tmp4 = fw[-1];
 819             fw[-4] = tmp1 + c1;
 820             fw[-3] = tmp2 + c2;
 821             fw[-2] = tmp3 + c3;
 822             fw[-1] = tmp4 + c4;
 823             fw += 8;
 824         }
 825     }
 826 }
 827
 828 #endif
 829
 830 /* <summary>                             */
 831 /* Inverse 9-7 wavelet transform in 1-D. */
 832 /* </summary>                            */
 833 static void v4dwt_decode(v4dwt_t* restrict dwt)
 834 {
 835     int a, b;
 836     if (dwt->cas == 0) {
 837         if (!((dwt->dn > 0) || (dwt->sn > 1))) {
 838             return;
 839         }
 840         a = 0;
 841         b = 1;
 842     } else {
 843         if (!((dwt->sn > 0) || (dwt->dn > 1))) {
 844             return;
 845         }
 846         a = 1;
 847         b = 0;
 848     }
 849 #ifdef __SSE__
 850     v4dwt_decode_step1_sse(dwt->wavelet + a, dwt->sn, _mm_set1_ps(K));
 851     v4dwt_decode_step1_sse(dwt->wavelet + b, dwt->dn, _mm_set1_ps(c13318));
 852     v4dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
 853                            int_min(dwt->sn, dwt->dn - a), _mm_set1_ps(dwt_delta));
 854     v4dwt_decode_step2_sse(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
 855                            int_min(dwt->dn, dwt->sn - b), _mm_set1_ps(dwt_gamma));
 856     v4dwt_decode_step2_sse(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
 857                            int_min(dwt->sn, dwt->dn - a), _mm_set1_ps(dwt_beta));
 858     v4dwt_decode_step2_sse(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
 859                            int_min(dwt->dn, dwt->sn - b), _mm_set1_ps(dwt_alpha));
 860 #else
 861     v4dwt_decode_step1(dwt->wavelet + a, dwt->sn, K);
 862     v4dwt_decode_step1(dwt->wavelet + b, dwt->dn, c13318);
 863     v4dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
 864                        int_min(dwt->sn, dwt->dn - a), dwt_delta);
 865     v4dwt_decode_step2(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
 866                        int_min(dwt->dn, dwt->sn - b), dwt_gamma);
 867     v4dwt_decode_step2(dwt->wavelet + b, dwt->wavelet + a + 1, dwt->sn,
 868                        int_min(dwt->sn, dwt->dn - a), dwt_beta);
 869     v4dwt_decode_step2(dwt->wavelet + a, dwt->wavelet + b + 1, dwt->dn,
 870                        int_min(dwt->dn, dwt->sn - b), dwt_alpha);
 871 #endif
 872 }
 873
 874 /* <summary>                             */
 875 /* Inverse 9-7 wavelet transform in 2-D. */
 876 /* </summary>                            */
 877 void dwt_decode_real(opj_tcd_tilecomp_t* restrict tilec, int numres)
 878 {
 879     v4dwt_t h;
 880     v4dwt_t v;
 881
 882     opj_tcd_resolution_t* res = tilec->resolutions;
 883
 884     int rw = res->x1 - res->x0; /* width of the resolution level computed */
 885     int rh = res->y1 - res->y0; /* height of the resolution level computed */
 886
 887     int w = tilec->x1 - tilec->x0;
 888
 889     h.wavelet = (v4*) opj_aligned_malloc((dwt_decode_max_resolution(res,
 890                                           numres) + 5) * sizeof(v4));
 891     v.wavelet = h.wavelet;
 892
 893     while (--numres) {
 894         float * restrict aj = (float*) tilec->data;
 895         int bufsize = (tilec->x1 - tilec->x0) * (tilec->y1 - tilec->y0);
 896         int j;
 897
 898         h.sn = rw;
 899         v.sn = rh;
 900
 901         ++res;
 902
 903         rw = res->x1 - res->x0; /* width of the resolution level computed */
 904         rh = res->y1 - res->y0; /* height of the resolution level computed */
 905
 906         h.dn = rw - h.sn;
 907         h.cas = res->x0 % 2;
 908
 909         for (j = rh; j > 3; j -= 4) {
 910             int k;
 911             v4dwt_interleave_h(&h, aj, w, bufsize);
 912             v4dwt_decode(&h);
 913             for (k = rw; --k >= 0;) {
 914                 aj[k    ] = h.wavelet[k].f[0];
 915                 aj[k + w  ] = h.wavelet[k].f[1];
 916                 aj[k + w * 2] = h.wavelet[k].f[2];
 917                 aj[k + w * 3] = h.wavelet[k].f[3];
 918             }
 919             aj += w * 4;
 920             bufsize -= w * 4;
 921         }
 922         if (rh & 0x03) {
 923             int k;
 924             j = rh & 0x03;
 925             v4dwt_interleave_h(&h, aj, w, bufsize);
 926             v4dwt_decode(&h);
 927             for (k = rw; --k >= 0;) {
 928                 switch (j) {
 929                 case 3:
 930                     aj[k + w * 2] = h.wavelet[k].f[2];
 931                 case 2:
 932                     aj[k + w  ] = h.wavelet[k].f[1];
 933                 case 1:
 934                     aj[k    ] = h.wavelet[k].f[0];
 935                 }
 936             }
 937         }
 938
 939         v.dn = rh - v.sn;
 940         v.cas = res->y0 % 2;
 941
 942         aj = (float*) tilec->data;
 943         for (j = rw; j > 3; j -= 4) {
 944             int k;
 945             v4dwt_interleave_v(&v, aj, w);
 946             v4dwt_decode(&v);
 947             for (k = 0; k < rh; ++k) {
 948                 memcpy(&aj[k * w], &v.wavelet[k], 4 * sizeof(float));
 949             }
 950             aj += 4;
 951         }
 952         if (rw & 0x03) {
 953             int k;
 954             j = rw & 0x03;
 955             v4dwt_interleave_v(&v, aj, w);
 956             v4dwt_decode(&v);
 957             for (k = 0; k < rh; ++k) {
 958                 memcpy(&aj[k * w], &v.wavelet[k], j * sizeof(float));
 959             }
 960         }
 961     }
 962
 963     opj_aligned_free(h.wavelet);
 964 }
 965