libopenjpeg/dwt.c

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