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  * All rights reserved.
   9  *
  10  * Redistribution and use in source and binary forms, with or without
  11  * modification, are permitted provided that the following conditions
  12  * are met:
  13  * 1. Redistributions of source code must retain the above copyright
  14  *    notice, this list of conditions and the following disclaimer.
  15  * 2. Redistributions in binary form must reproduce the above copyright
  16  *    notice, this list of conditions and the following disclaimer in the
  17  *    documentation and/or other materials provided with the distribution.
  18  *
  19  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
  20  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  21  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  22  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  23  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  24  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  25  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  26  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  27  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  28  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  29  * POSSIBILITY OF SUCH DAMAGE.
  30  */
  31
  32 /*
  33  *  NOTE:
  34  *  This is a modified version of the openjpeg dwt.c file.
  35  *  Average speed improvement compared to the original file (measured on
  36  *  my own machine, a P4 running at 3.0 GHz):
  37  *  5x3 wavelets about 2 times faster
  38  *  9x7 wavelets about 3 times faster
  39  *  for both, encoding and decoding.
  40  *
  41  *  The better performance is caused by doing the 1-dimensional DWT
  42  *  within a temporary buffer where the data can be accessed sequential
  43  *  for both directions, horizontal and vertical. The 2d vertical DWT was
  44  *  the major bottleneck in the former version.
  45  *
  46  *  I have also removed the "Add Patrick" part because it is not longer
  47  *  needed.
  48  *
  49  *  6/6/2005
  50  *  -Ive (aka Reiner Wahler)
  51  *  mail: ive@lilysoft.com
  52  */
  53
  54 #include "opj_includes.h"
  55
  56 /** @defgroup DWT DWT - Implementation of a discrete wavelet transform */
  57 /*@{*/
  58
  59 /** @name Local static functions */
  60 /*@{*/
  61
  62 /**
  63 Forward lazy transform (horizontal)
  64 */
  65 static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas);
  66 /**
  67 Forward lazy transform (vertical)
  68 */
  69 static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas);
  70 /**
  71 Inverse lazy transform (horizontal)
  72 */
  73 static void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas);
  74 /**
  75 Inverse lazy transform (vertical)
  76 */
  77 static void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas);
  78 /**
  79 Forward 5-3 wavelet tranform in 1-D
  80 */
  81 static void dwt_encode_1(int *a, int dn, int sn, int cas);
  82 /**
  83 Inverse 5-3 wavelet tranform in 1-D
  84 */
  85 static void dwt_decode_1(int *a, int dn, int sn, int cas);
  86 /**
  87 Forward 9-7 wavelet transform in 1-D
  88 */
  89 static void dwt_encode_1_real(int *a, int dn, int sn, int cas);
  90 /**
  91 Inverse 9-7 wavelet transform in 1-D
  92 */
  93 static void dwt_decode_1_real(int *a, int dn, int sn, int cas);
  94 /**
  95 FIXME : comment ???
  96 */
  97 static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno_stepsize);
  98
  99 /*@}*/
 100
 101 /*@}*/
 102
 103 #define S(i) a[(i)*2]
 104 #define D(i) a[(1+(i)*2)]
 105 #define S_(i) ((i)<0?S(0):((i)>=sn?S(sn-1):S(i)))
 106 #define D_(i) ((i)<0?D(0):((i)>=dn?D(dn-1):D(i)))
 107 /* new */
 108 #define SS_(i) ((i)<0?S(0):((i)>=dn?S(dn-1):S(i)))
 109 #define DD_(i) ((i)<0?D(0):((i)>=sn?D(sn-1):D(i)))
 110
 111 /* <summary>                                                              */
 112 /* This table contains the norms of the 5-3 wavelets for different bands. */
 113 /* </summary>                                                             */
 114 static const double dwt_norms[4][10] = {
 115         {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 116         {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 117         {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 118         {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
 119 };
 120
 121 /* <summary>                                                              */
 122 /* This table contains the norms of the 9-7 wavelets for different bands. */
 123 /* </summary>                                                             */
 124 static const double dwt_norms_real[4][10] = {
 125         {1.000, 1.965, 4.177, 8.403, 16.90, 33.84, 67.69, 135.3, 270.6, 540.9},
 126         {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
 127         {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
 128         {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2}
 129 };
 130
 131 /*
 132 ==========================================================
 133    local functions
 134 ==========================================================
 135 */
 136
 137 /* <summary>                                     */
 138 /* Forward lazy transform (horizontal).  */
 139 /* </summary>                            */
 140 static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas) {
 141         int i;
 142     for (i=0; i<sn; i++) b[i]=a[2*i+cas];
 143     for (i=0; i<dn; i++) b[sn+i]=a[(2*i+1-cas)];
 144 }
 145
 146 /* <summary>                             */
 147 /* Forward lazy transform (vertical).    */
 148 /* </summary>                            */
 149 static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
 150     int i;
 151     for (i=0; i<sn; i++) b[i*x]=a[2*i+cas];
 152     for (i=0; i<dn; i++) b[(sn+i)*x]=a[(2*i+1-cas)];
 153 }
 154
 155 /* <summary>                             */
 156 /* Inverse lazy transform (horizontal).  */
 157 /* </summary>                            */
 158 static void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas) {
 159     int i;
 160     int *ai = NULL;
 161     int *bi = NULL;
 162     ai = a;
 163     bi = b + cas;
 164     for (i = 0; i < sn; i++) {
 165       *bi = *ai;
 166           bi += 2;
 167           ai++;
 168     }
 169     ai = a + sn;
 170     bi = b + 1 - cas;
 171     for (i = 0; i < dn; i++) {
 172       *bi = *ai;
 173           bi += 2;
 174           ai++;
 175     }
 176 }
 177
 178 /* <summary>                             */
 179 /* Inverse lazy transform (vertical).    */
 180 /* </summary>                            */
 181 static void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
 182     int i;
 183     int *ai = NULL;
 184     int *bi = NULL;
 185     ai = a;
 186     bi = b + cas;
 187     for (i = 0; i < sn; i++) {
 188       *bi = *ai;
 189           bi += 2;
 190           ai += x;
 191     }
 192     ai = a + (sn * x);
 193     bi = b + 1 - cas;
 194     for (i = 0; i < dn; i++) {
 195       *bi = *ai;
 196           bi += 2;
 197           ai += x;
 198     }
 199 }
 200
 201
 202 /* <summary>                            */
 203 /* Forward 5-3 wavelet tranform in 1-D. */
 204 /* </summary>                           */
 205 static void dwt_encode_1(int *a, int dn, int sn, int cas) {
 206         int i;
 207
 208         if (!cas) {
 209                 if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 210                         for (i = 0; i < dn; i++) D(i) -= (S_(i) + S_(i + 1)) >> 1;
 211                         for (i = 0; i < sn; i++) S(i) += (D_(i - 1) + D_(i) + 2) >> 2;
 212                 }
 213         } else {
 214                 if (!sn && dn == 1)                 /* NEW :  CASE ONE ELEMENT */
 215                         S(0) *= 2;
 216                 else {
 217                         for (i = 0; i < dn; i++) S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
 218                         for (i = 0; i < sn; i++) D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
 219                 }
 220         }
 221 }
 222
 223 /* <summary>                            */
 224 /* Inverse 5-3 wavelet tranform in 1-D. */
 225 /* </summary>                           */
 226 static void dwt_decode_1(int *a, int dn, int sn, int cas) {
 227         int i;
 228
 229         if (!cas) {
 230                 if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 231                         for (i = 0; i < sn; i++) S(i) -= (D_(i - 1) + D_(i) + 2) >> 2;
 232                         for (i = 0; i < dn; i++) D(i) += (S_(i) + S_(i + 1)) >> 1;
 233                 }
 234         } else {
 235                 if (!sn  && dn == 1)          /* NEW :  CASE ONE ELEMENT */
 236                         S(0) /= 2;
 237                 else {
 238                         for (i = 0; i < sn; i++) D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2;
 239                         for (i = 0; i < dn; i++) S(i) += (DD_(i) + DD_(i - 1)) >> 1;
 240                 }
 241         }
 242 }
 243
 244 /* <summary>                             */
 245 /* Forward 9-7 wavelet transform in 1-D. */
 246 /* </summary>                            */
 247 static void dwt_encode_1_real(int *a, int dn, int sn, int cas) {
 248         int i;
 249         if (!cas) {
 250                 if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 251                         for (i = 0; i < dn; i++)
 252                                 D(i) -= fix_mul(S_(i) + S_(i + 1), 12993);
 253                         for (i = 0; i < sn; i++)
 254                                 S(i) -= fix_mul(D_(i - 1) + D_(i), 434);
 255                         for (i = 0; i < dn; i++)
 256                                 D(i) += fix_mul(S_(i) + S_(i + 1), 7233);
 257                         for (i = 0; i < sn; i++)
 258                                 S(i) += fix_mul(D_(i - 1) + D_(i), 3633);
 259                         for (i = 0; i < dn; i++)
 260                                 D(i) = fix_mul(D(i), 5038);     /*5038 */
 261                         for (i = 0; i < sn; i++)
 262                                 S(i) = fix_mul(S(i), 6659);     /*6660 */
 263                 }
 264         } else {
 265                 if ((sn > 0) || (dn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 266                         for (i = 0; i < dn; i++)
 267                                 S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993);
 268                         for (i = 0; i < sn; i++)
 269                                 D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434);
 270                         for (i = 0; i < dn; i++)
 271                                 S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233);
 272                         for (i = 0; i < sn; i++)
 273                                 D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633);
 274                         for (i = 0; i < dn; i++)
 275                                 S(i) = fix_mul(S(i), 5038);     /*5038 */
 276                         for (i = 0; i < sn; i++)
 277                                 D(i) = fix_mul(D(i), 6659);     /*6660 */
 278                 }
 279         }
 280 }
 281
 282 /* <summary>                             */
 283 /* Inverse 9-7 wavelet transform in 1-D. */
 284 /* </summary>                            */
 285 static void dwt_decode_1_real(int *a, int dn, int sn, int cas) {
 286         int i;
 287         if (!cas) {
 288                 if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 289                         for (i = 0; i < sn; i++)
 290                                 S(i) = fix_mul(S(i), 10078);    /* 10076 */
 291                         for (i = 0; i < dn; i++)
 292                                 D(i) = fix_mul(D(i), 13318);    /* 13320 */
 293                         for (i = 0; i < sn; i++)
 294                                 S(i) -= fix_mul(D_(i - 1) + D_(i), 3633);
 295                         for (i = 0; i < dn; i++)
 296                                 D(i) -= fix_mul(S_(i) + S_(i + 1), 7233);
 297                         for (i = 0; i < sn; i++)
 298                                 S(i) += fix_mul(D_(i - 1) + D_(i), 434);
 299                         for (i = 0; i < dn; i++)
 300                                 D(i) += fix_mul(S_(i) + S_(i + 1), 12994);      /* 12993 */
 301                 }
 302         } else {
 303                 if ((sn > 0) || (dn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 304                         for (i = 0; i < sn; i++)
 305                                 D(i) = fix_mul(D(i), 10078);    /* 10076 */
 306                         for (i = 0; i < dn; i++)
 307                                 S(i) = fix_mul(S(i), 13318);    /* 13320 */
 308                         for (i = 0; i < sn; i++)
 309                                 D(i) -= fix_mul(SS_(i) + SS_(i + 1), 3633);
 310                         for (i = 0; i < dn; i++)
 311                                 S(i) -= fix_mul(DD_(i) + DD_(i - 1), 7233);
 312                         for (i = 0; i < sn; i++)
 313                                 D(i) += fix_mul(SS_(i) + SS_(i + 1), 434);
 314                         for (i = 0; i < dn; i++)
 315                                 S(i) += fix_mul(DD_(i) + DD_(i - 1), 12994);    /* 12993 */
 316                 }
 317         }
 318 }
 319
 320 static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno_stepsize) {
 321         int p, n;
 322         p = int_floorlog2(stepsize) - 13;
 323         n = 11 - int_floorlog2(stepsize);
 324         bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
 325         bandno_stepsize->expn = numbps - p;
 326 }
 327
 328 /*
 329 ==========================================================
 330    DWT interface
 331 ==========================================================
 332 */
 333
 334 /* <summary>                            */
 335 /* Forward 5-3 wavelet tranform in 2-D. */
 336 /* </summary>                           */
 337 void dwt_encode(opj_tcd_tilecomp_t * tilec) {
 338         int i, j, k;
 339         int *a = NULL;
 340         int *aj = NULL;
 341         int *bj = NULL;
 342         int w, l;
 343
 344         w = tilec->x1-tilec->x0;
 345         l = tilec->numresolutions-1;
 346         a = tilec->data;
 347
 348         for (i = 0; i < l; i++) {
 349                 int rw;                 /* width of the resolution level computed                                                           */
 350                 int rh;                 /* heigth of the resolution level computed                                                          */
 351                 int rw1;                /* width of the resolution level once lower than computed one                                       */
 352                 int rh1;                /* height of the resolution level once lower than computed one                                      */
 353                 int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 354                 int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 355                 int dn, sn;
 356
 357                 rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 358                 rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 359                 rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 360                 rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 361
 362                 cas_row = tilec->resolutions[l - i].x0 % 2;
 363                 cas_col = tilec->resolutions[l - i].y0 % 2;
 364
 365                 sn = rh1;
 366                 dn = rh - rh1;
 367                 bj = (int*)opj_malloc(rh * sizeof(int));
 368                 for (j = 0; j < rw; j++) {
 369                         aj = a + j;
 370                         for (k = 0; k < rh; k++)  bj[k] = aj[k*w];
 371                         dwt_encode_1(bj, dn, sn, cas_col);
 372                         dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 373                 }
 374                 opj_free(bj);
 375
 376                 sn = rw1;
 377                 dn = rw - rw1;
 378                 bj = (int*)opj_malloc(rw * sizeof(int));
 379                 for (j = 0; j < rh; j++) {
 380                         aj = a + j * w;
 381                         for (k = 0; k < rw; k++)  bj[k] = aj[k];
 382                         dwt_encode_1(bj, dn, sn, cas_row);
 383                         dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 384                 }
 385                 opj_free(bj);
 386         }
 387 }
 388
 389
 390 /* <summary>                            */
 391 /* Inverse 5-3 wavelet tranform in 2-D. */
 392 /* </summary>                           */
 393 void dwt_decode(opj_tcd_tilecomp_t * tilec, int stop) {
 394         int i, j, k;
 395         int *a = NULL;
 396         int *aj = NULL;
 397         int *bj = NULL;
 398         int w, l;
 399
 400         w = tilec->x1-tilec->x0;
 401         l = tilec->numresolutions-1;
 402         a = tilec->data;
 403
 404         for (i = l - 1; i >= stop; i--) {
 405                 int rw;                 /* width of the resolution level computed                                                           */
 406                 int rh;                 /* heigth of the resolution level computed                                                          */
 407                 int rw1;                /* width of the resolution level once lower than computed one                                       */
 408                 int rh1;                /* height of the resolution level once lower than computed one                                      */
 409                 int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 410                 int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 411                 int dn, sn;
 412
 413                 rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 414                 rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 415                 rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 416                 rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 417
 418                 cas_row = tilec->resolutions[l - i].x0 % 2;
 419                 cas_col = tilec->resolutions[l - i].y0 % 2;
 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                         dwt_interleave_h(aj, bj, dn, sn, cas_row);
 427                         dwt_decode_1(bj, dn, sn, cas_row);
 428                         for (k = 0; k < rw; k++)  aj[k] = bj[k];
 429                 }
 430                 opj_free(bj);
 431
 432                 sn = rh1;
 433                 dn = rh - rh1;
 434                 bj = (int*)opj_malloc(rh * sizeof(int));
 435                 for (j = 0; j < rw; j++) {
 436                         aj = a + j;
 437                         dwt_interleave_v(aj, bj, dn, sn, w, cas_col);
 438                         dwt_decode_1(bj, dn, sn, cas_col);
 439                         for (k = 0; k < rh; k++)  aj[k * w] = bj[k];
 440                 }
 441                 opj_free(bj);
 442         }
 443 }
 444
 445
 446 /* <summary>                          */
 447 /* Get gain of 5-3 wavelet transform. */
 448 /* </summary>                         */
 449 int dwt_getgain(int orient) {
 450         if (orient == 0)
 451                 return 0;
 452         if (orient == 1 || orient == 2)
 453                 return 1;
 454         return 2;
 455 }
 456
 457 /* <summary>                */
 458 /* Get norm of 5-3 wavelet. */
 459 /* </summary>               */
 460 double dwt_getnorm(int level, int orient) {
 461         return dwt_norms[orient][level];
 462 }
 463
 464 /* <summary>                             */
 465 /* Forward 9-7 wavelet transform in 2-D. */
 466 /* </summary>                            */
 467
 468 void dwt_encode_real(opj_tcd_tilecomp_t * tilec) {
 469         int i, j, k;
 470         int *a = NULL;
 471         int *aj = NULL;
 472         int *bj = NULL;
 473         int w, l;
 474
 475         w = tilec->x1-tilec->x0;
 476         l = tilec->numresolutions-1;
 477         a = tilec->data;
 478
 479         for (i = 0; i < l; i++) {
 480                 int rw;                 /* width of the resolution level computed                                                     */
 481                 int rh;                 /* heigth of the resolution level computed                                                    */
 482                 int rw1;                /* width of the resolution level once lower than computed one                                 */
 483                 int rh1;                /* height of the resolution level once lower than computed one                                */
 484                 int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 485                 int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 486                 int dn, sn;
 487
 488                 rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 489                 rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 490                 rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 491                 rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 492
 493                 cas_row = tilec->resolutions[l - i].x0 % 2;
 494                 cas_col = tilec->resolutions[l - i].y0 % 2;
 495
 496                 sn = rh1;
 497                 dn = rh - rh1;
 498                 bj = (int*)opj_malloc(rh * sizeof(int));
 499                 for (j = 0; j < rw; j++) {
 500                         aj = a + j;
 501                         for (k = 0; k < rh; k++)  bj[k] = aj[k*w];
 502                         dwt_encode_1_real(bj, dn, sn, cas_col);
 503                         dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 504                 }
 505                 opj_free(bj);
 506
 507                 sn = rw1;
 508                 dn = rw - rw1;
 509                 bj = (int*)opj_malloc(rw * sizeof(int));
 510                 for (j = 0; j < rh; j++) {
 511                         aj = a + j * w;
 512                         for (k = 0; k < rw; k++)  bj[k] = aj[k];
 513                         dwt_encode_1_real(bj, dn, sn, cas_row);
 514                         dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 515                 }
 516                 opj_free(bj);
 517         }
 518 }
 519
 520
 521 /* <summary>                             */
 522 /* Inverse 9-7 wavelet transform in 2-D. */
 523 /* </summary>                            */
 524 void dwt_decode_real(opj_tcd_tilecomp_t * tilec, int stop) {
 525         int i, j, k;
 526         int *a = NULL;
 527         int *aj = NULL;
 528         int *bj = NULL;
 529         int w, l;
 530
 531         w = tilec->x1-tilec->x0;
 532         l = tilec->numresolutions-1;
 533         a = tilec->data;
 534
 535         for (i = l-1; i >= stop; i--) {
 536                 int rw;                 /* width of the resolution level computed                       */
 537                 int rh;                 /* heigth of the resolution level computed                      */
 538                 int rw1;                /* width of the resolution level once lower than computed one   */
 539                 int rh1;                /* height of the resolution level once lower than computed one  */
 540                 int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 541                 int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 542                 int dn, sn;
 543
 544                 rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 545                 rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 546                 rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 547                 rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 548
 549                 cas_col = tilec->resolutions[l - i].x0 % 2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 550                 cas_row = tilec->resolutions[l - i].y0 % 2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 551
 552                 sn = rw1;
 553                 dn = rw-rw1;
 554                 bj = (int*)opj_malloc(rw * sizeof(int));
 555                 for (j = 0; j < rh; j++) {
 556                         aj = a + j * w;
 557                         dwt_interleave_h(aj, bj, dn, sn, cas_col);
 558                         dwt_decode_1_real(bj, dn, sn, cas_col);
 559                         for (k = 0; k < rw; k++)  aj[k] = bj[k];
 560                 }
 561                 opj_free(bj);
 562
 563                 sn = rh1;
 564                 dn = rh-rh1;
 565                 bj = (int*)opj_malloc(rh * sizeof(int));
 566                 for (j = 0; j < rw; j++) {
 567                         aj = a + j;
 568                         dwt_interleave_v(aj, bj, dn, sn, w, cas_row);
 569                         dwt_decode_1_real(bj, dn, sn, cas_row);
 570                         for (k = 0; k < rh; k++)  aj[k * w] = bj[k];
 571                 }
 572                 opj_free(bj);
 573         }
 574 }
 575
 576
 577 /* <summary>                          */
 578 /* Get gain of 9-7 wavelet transform. */
 579 /* </summary>                         */
 580 int dwt_getgain_real(int orient) {
 581         (void)orient;
 582         return 0;
 583 }
 584
 585 /* <summary>                */
 586 /* Get norm of 9-7 wavelet. */
 587 /* </summary>               */
 588 double dwt_getnorm_real(int level, int orient) {
 589         return dwt_norms_real[orient][level];
 590 }
 591
 592 void dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, int prec) {
 593         int numbands, bandno;
 594         numbands = 3 * tccp->numresolutions - 2;
 595         for (bandno = 0; bandno < numbands; bandno++) {
 596                 double stepsize;
 597                 int resno, level, orient, gain;
 598
 599                 resno = (bandno == 0) ? 0 : ((bandno - 1) / 3 + 1);
 600                 orient = (bandno == 0) ? 0 : ((bandno - 1) % 3 + 1);
 601                 level = tccp->numresolutions - 1 - resno;
 602                 gain = (tccp->qmfbid == 0) ? 0 : ((orient == 0) ? 0 : (((orient == 1) || (orient == 2)) ? 1 : 2));
 603                 if (tccp->qntsty == J2K_CCP_QNTSTY_NOQNT) {
 604                         stepsize = 1.0;
 605                 } else {
 606                         double norm = dwt_norms_real[orient][level];
 607                         stepsize = (1 << (gain)) / norm;
 608                 }
 609                 dwt_encode_stepsize((int) floor(stepsize * 8192.0), prec + gain, &tccp->stepsizes[bandno]);
 610         }
 611 }