src/gdcmopenjpeg/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 }