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
  54 #include "opj_includes.h"
  55
  56 #define S(i) a[(i)*2]
  57 #define D(i) a[(1+(i)*2)]
  58 #define S_(i) ((i)<0?S(0):((i)>=sn?S(sn-1):S(i)))
  59 #define D_(i) ((i)<0?D(0):((i)>=dn?D(dn-1):D(i)))
  60 /* new */
  61 #define SS_(i) ((i)<0?S(0):((i)>=dn?S(dn-1):S(i)))
  62 #define DD_(i) ((i)<0?D(0):((i)>=sn?D(sn-1):D(i)))
  63
  64 /* <summary>                                                              */
  65 /* This table contains the norms of the 5-3 wavelets for different bands. */
  66 /* </summary>                                                             */
  67 static const double dwt_norms[4][10] = {
  68   {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
  69   {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
  70   {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
  71   {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
  72 };
  73
  74 /* <summary>                                                              */
  75 /* This table contains the norms of the 9-7 wavelets for different bands. */
  76 /* </summary>                                                             */
  77 static const double dwt_norms_real[4][10] = {
  78   {1.000, 1.965, 4.177, 8.403, 16.90, 33.84, 67.69, 135.3, 270.6, 540.9},
  79   {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
  80   {2.022, 3.989, 8.355, 17.04, 34.27, 68.63, 137.3, 274.6, 549.0},
  81   {2.080, 3.865, 8.307, 17.18, 34.71, 69.59, 139.3, 278.6, 557.2}
  82 };
  83
  84 /*
  85 ==========================================================
  86    local functions
  87 ==========================================================
  88 */
  89
  90 /* <summary>                       */
  91 /* Forward lazy transform (horizontal).  */
  92 /* </summary>                            */
  93 static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas) {
  94   int i;
  95     for (i=0; i<sn; i++) b[i]=a[2*i+cas];
  96     for (i=0; i<dn; i++) b[sn+i]=a[(2*i+1-cas)];
  97 }
  98
  99 /* <summary>                             */
 100 /* Forward lazy transform (vertical).    */
 101 /* </summary>                            */
 102 static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
 103     int i;
 104     for (i=0; i<sn; i++) b[i*x]=a[2*i+cas];
 105     for (i=0; i<dn; i++) b[(sn+i)*x]=a[(2*i+1-cas)];
 106 }
 107
 108 /* <summary>                             */
 109 /* Inverse lazy transform (horizontal).  */
 110 /* </summary>                            */
 111 static void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas) {
 112     int i;
 113     int *ai = NULL;
 114     int *bi = NULL;
 115     ai = a;
 116     bi = b + cas;
 117     for (i = 0; i < sn; i++) {
 118       *bi = *ai;
 119     bi += 2;
 120     ai++;
 121     }
 122     ai = a + sn;
 123     bi = b + 1 - cas;
 124     for (i = 0; i < dn; i++) {
 125       *bi = *ai;
 126     bi += 2;
 127     ai++;
 128     }
 129 }
 130
 131 /* <summary>                             */
 132 /* Inverse lazy transform (vertical).    */
 133 /* </summary>                            */
 134 static void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
 135     int i;
 136     int *ai = NULL;
 137     int *bi = NULL;
 138     ai = a;
 139     bi = b + cas;
 140     for (i = 0; i < sn; i++) {
 141       *bi = *ai;
 142     bi += 2;
 143     ai += x;
 144     }
 145     ai = a + (sn * x);
 146     bi = b + 1 - cas;
 147     for (i = 0; i < dn; i++) {
 148       *bi = *ai;
 149     bi += 2;
 150     ai += x;
 151     }
 152 }
 153
 154
 155 /* <summary>                            */
 156 /* Forward 5-3 wavelet tranform in 1-D. */
 157 /* </summary>                           */
 158 static void dwt_encode_1(int *a, int dn, int sn, int cas) {
 159   int i;
 160
 161   if (!cas) {
 162     if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 163       for (i = 0; i < dn; i++) D(i) -= (S_(i) + S_(i + 1)) >> 1;
 164       for (i = 0; i < sn; i++) S(i) += (D_(i - 1) + D_(i) + 2) >> 2;
 165     }
 166   } else {
 167     if (!sn && dn == 1)       /* NEW :  CASE ONE ELEMENT */
 168       S(0) *= 2;
 169     else {
 170       for (i = 0; i < dn; i++) S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
 171       for (i = 0; i < sn; i++) D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
 172     }
 173   }
 174 }
 175
 176 /* <summary>                            */
 177 /* Inverse 5-3 wavelet tranform in 1-D. */
 178 /* </summary>                           */
 179 static void dwt_decode_1(int *a, int dn, int sn, int cas) {
 180   int i;
 181
 182   if (!cas) {
 183     if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 184       for (i = 0; i < sn; i++) S(i) -= (D_(i - 1) + D_(i) + 2) >> 2;
 185       for (i = 0; i < dn; i++) D(i) += (S_(i) + S_(i + 1)) >> 1;
 186     }
 187   } else {
 188     if (!sn  && dn == 1)          /* NEW :  CASE ONE ELEMENT */
 189       S(0) /= 2;
 190     else {
 191       for (i = 0; i < sn; i++) D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2;
 192       for (i = 0; i < dn; i++) S(i) += (DD_(i) + DD_(i - 1)) >> 1;
 193     }
 194   }
 195 }
 196
 197 /* <summary>                             */
 198 /* Forward 9-7 wavelet transform in 1-D. */
 199 /* </summary>                            */
 200 static void dwt_encode_1_real(int *a, int dn, int sn, int cas) {
 201   int i;
 202   if (!cas) {
 203     if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 204       for (i = 0; i < dn; i++)
 205         D(i) -= fix_mul(S_(i) + S_(i + 1), 12993);
 206       for (i = 0; i < sn; i++)
 207         S(i) -= fix_mul(D_(i - 1) + D_(i), 434);
 208       for (i = 0; i < dn; i++)
 209         D(i) += fix_mul(S_(i) + S_(i + 1), 7233);
 210       for (i = 0; i < sn; i++)
 211         S(i) += fix_mul(D_(i - 1) + D_(i), 3633);
 212       for (i = 0; i < dn; i++)
 213         D(i) = fix_mul(D(i), 5038); /*5038 */
 214       for (i = 0; i < sn; i++)
 215         S(i) = fix_mul(S(i), 6659); /*6660 */
 216     }
 217   } else {
 218     if ((sn > 0) || (dn > 1)) { /* NEW :  CASE ONE ELEMENT */
 219       for (i = 0; i < dn; i++)
 220         S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993);
 221       for (i = 0; i < sn; i++)
 222         D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434);
 223       for (i = 0; i < dn; i++)
 224         S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233);
 225       for (i = 0; i < sn; i++)
 226         D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633);
 227       for (i = 0; i < dn; i++)
 228         S(i) = fix_mul(S(i), 5038); /*5038 */
 229       for (i = 0; i < sn; i++)
 230         D(i) = fix_mul(D(i), 6659); /*6660 */
 231     }
 232   }
 233 }
 234
 235 /* <summary>                             */
 236 /* Inverse 9-7 wavelet transform in 1-D. */
 237 /* </summary>                            */
 238 static void dwt_decode_1_real(int *a, int dn, int sn, int cas) {
 239   int i;
 240   if (!cas) {
 241     if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 242       for (i = 0; i < sn; i++)
 243         S(i) = fix_mul(S(i), 10078);  /* 10076 */
 244       for (i = 0; i < dn; i++)
 245         D(i) = fix_mul(D(i), 13318);  /* 13320 */
 246       for (i = 0; i < sn; i++)
 247         S(i) -= fix_mul(D_(i - 1) + D_(i), 3633);
 248       for (i = 0; i < dn; i++)
 249         D(i) -= fix_mul(S_(i) + S_(i + 1), 7233);
 250       for (i = 0; i < sn; i++)
 251         S(i) += fix_mul(D_(i - 1) + D_(i), 434);
 252       for (i = 0; i < dn; i++)
 253         D(i) += fix_mul(S_(i) + S_(i + 1), 12994);  /* 12993 */
 254     }
 255   } else {
 256     if ((sn > 0) || (dn > 1)) { /* NEW :  CASE ONE ELEMENT */
 257       for (i = 0; i < sn; i++)
 258         D(i) = fix_mul(D(i), 10078);  /* 10076 */
 259       for (i = 0; i < dn; i++)
 260         S(i) = fix_mul(S(i), 13318);  /* 13320 */
 261       for (i = 0; i < sn; i++)
 262         D(i) -= fix_mul(SS_(i) + SS_(i + 1), 3633);
 263       for (i = 0; i < dn; i++)
 264         S(i) -= fix_mul(DD_(i) + DD_(i - 1), 7233);
 265       for (i = 0; i < sn; i++)
 266         D(i) += fix_mul(SS_(i) + SS_(i + 1), 434);
 267       for (i = 0; i < dn; i++)
 268         S(i) += fix_mul(DD_(i) + DD_(i - 1), 12994);  /* 12993 */
 269     }
 270   }
 271 }
 272
 273 static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno_stepsize) {
 274   int p, n;
 275   p = int_floorlog2(stepsize) - 13;
 276   n = 11 - int_floorlog2(stepsize);
 277   bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
 278   bandno_stepsize->expn = numbps - p;
 279 }
 280
 281 /*
 282 ==========================================================
 283    DWT interface
 284 ==========================================================
 285 */
 286
 287 /* <summary>                            */
 288 /* Forward 5-3 wavelet tranform in 2-D. */
 289 /* </summary>                           */
 290 void dwt_encode(opj_tcd_tilecomp_t * tilec) {
 291   int i, j, k;
 292   int *a = NULL;
 293   int *aj = NULL;
 294   int *bj = NULL;
 295   int w, l;
 296
 297   w = tilec->x1-tilec->x0;
 298   l = tilec->numresolutions-1;
 299   a = tilec->data;
 300
 301   for (i = 0; i < l; i++) {
 302     int rw;     /* width of the resolution level computed                                                           */
 303     int rh;     /* heigth of the resolution level computed                                                          */
 304     int rw1;    /* width of the resolution level once lower than computed one                                       */
 305     int rh1;    /* height of the resolution level once lower than computed one                                      */
 306     int cas_col;  /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 307     int cas_row;  /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 308     int dn, sn;
 309
 310     rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 311     rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 312     rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 313     rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 314
 315     cas_row = tilec->resolutions[l - i].x0 % 2;
 316     cas_col = tilec->resolutions[l - i].y0 % 2;
 317
 318     sn = rh1;
 319     dn = rh - rh1;
 320     bj = (int*)opj_malloc(rh * sizeof(int));
 321     for (j = 0; j < rw; j++) {
 322       aj = a + j;
 323       for (k = 0; k < rh; k++)  bj[k] = aj[k*w];
 324       dwt_encode_1(bj, dn, sn, cas_col);
 325       dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 326     }
 327     opj_free(bj);
 328
 329     sn = rw1;
 330     dn = rw - rw1;
 331     bj = (int*)opj_malloc(rw * sizeof(int));
 332     for (j = 0; j < rh; j++) {
 333       aj = a + j * w;
 334       for (k = 0; k < rw; k++)  bj[k] = aj[k];
 335       dwt_encode_1(bj, dn, sn, cas_row);
 336       dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 337     }
 338     opj_free(bj);
 339   }
 340 }
 341
 342
 343 /* <summary>                            */
 344 /* Inverse 5-3 wavelet tranform in 2-D. */
 345 /* </summary>                           */
 346 void dwt_decode(opj_tcd_tilecomp_t * tilec, int stop) {
 347   int i, j, k;
 348   int *a = NULL;
 349   int *aj = NULL;
 350   int *bj = NULL;
 351   int w, l;
 352
 353   w = tilec->x1-tilec->x0;
 354   l = tilec->numresolutions-1;
 355   a = tilec->data;
 356
 357   for (i = l - 1; i >= stop; i--) {
 358     int rw;     /* width of the resolution level computed                                                           */
 359     int rh;     /* heigth of the resolution level computed                                                          */
 360     int rw1;    /* width of the resolution level once lower than computed one                                       */
 361     int rh1;    /* height of the resolution level once lower than computed one                                      */
 362     int cas_col;  /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 363     int cas_row;  /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 364     int dn, sn;
 365
 366     rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 367     rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 368     rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 369     rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 370
 371     cas_row = tilec->resolutions[l - i].x0 % 2;
 372     cas_col = tilec->resolutions[l - i].y0 % 2;
 373
 374     sn = rw1;
 375     dn = rw - rw1;
 376     bj = (int*)opj_malloc(rw * sizeof(int));
 377     for (j = 0; j < rh; j++) {
 378       aj = a + j*w;
 379       dwt_interleave_h(aj, bj, dn, sn, cas_row);
 380       dwt_decode_1(bj, dn, sn, cas_row);
 381       for (k = 0; k < rw; k++)  aj[k] = bj[k];
 382     }
 383     opj_free(bj);
 384
 385     sn = rh1;
 386     dn = rh - rh1;
 387     bj = (int*)opj_malloc(rh * sizeof(int));
 388     for (j = 0; j < rw; j++) {
 389       aj = a + j;
 390       dwt_interleave_v(aj, bj, dn, sn, w, cas_col);
 391       dwt_decode_1(bj, dn, sn, cas_col);
 392       for (k = 0; k < rh; k++)  aj[k * w] = bj[k];
 393     }
 394     opj_free(bj);
 395   }
 396 }
 397
 398
 399 /* <summary>                          */
 400 /* Get gain of 5-3 wavelet transform. */
 401 /* </summary>                         */
 402 int dwt_getgain(int orient) {
 403   if (orient == 0)
 404     return 0;
 405   if (orient == 1 || orient == 2)
 406     return 1;
 407   return 2;
 408 }
 409
 410 /* <summary>                */
 411 /* Get norm of 5-3 wavelet. */
 412 /* </summary>               */
 413 double dwt_getnorm(int level, int orient) {
 414   return dwt_norms[orient][level];
 415 }
 416
 417 /* <summary>                             */
 418 /* Forward 9-7 wavelet transform in 2-D. */
 419 /* </summary>                            */
 420
 421 void dwt_encode_real(opj_tcd_tilecomp_t * tilec) {
 422   int i, j, k;
 423   int *a = NULL;
 424   int *aj = NULL;
 425   int *bj = NULL;
 426   int w, l;
 427
 428   w = tilec->x1-tilec->x0;
 429   l = tilec->numresolutions-1;
 430   a = tilec->data;
 431
 432   for (i = 0; i < l; i++) {
 433     int rw;     /* width of the resolution level computed                                                     */
 434     int rh;     /* heigth of the resolution level computed                                                    */
 435     int rw1;    /* width of the resolution level once lower than computed one                                 */
 436     int rh1;    /* height of the resolution level once lower than computed one                                */
 437     int cas_col;  /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 438     int cas_row;  /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 439     int dn, sn;
 440
 441     rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 442     rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 443     rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 444     rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 445
 446     cas_row = tilec->resolutions[l - i].x0 % 2;
 447     cas_col = tilec->resolutions[l - i].y0 % 2;
 448
 449     sn = rh1;
 450     dn = rh - rh1;
 451     bj = (int*)opj_malloc(rh * sizeof(int));
 452     for (j = 0; j < rw; j++) {
 453       aj = a + j;
 454       for (k = 0; k < rh; k++)  bj[k] = aj[k*w];
 455       dwt_encode_1_real(bj, dn, sn, cas_col);
 456       dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 457     }
 458     opj_free(bj);
 459
 460     sn = rw1;
 461     dn = rw - rw1;
 462     bj = (int*)opj_malloc(rw * sizeof(int));
 463     for (j = 0; j < rh; j++) {
 464       aj = a + j * w;
 465       for (k = 0; k < rw; k++)  bj[k] = aj[k];
 466       dwt_encode_1_real(bj, dn, sn, cas_row);
 467       dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 468     }
 469     opj_free(bj);
 470   }
 471 }
 472
 473
 474 /* <summary>                             */
 475 /* Inverse 9-7 wavelet transform in 2-D. */
 476 /* </summary>                            */
 477 void dwt_decode_real(opj_tcd_tilecomp_t * tilec, int stop) {
 478   int i, j, k;
 479   int *a = NULL;
 480   int *aj = NULL;
 481   int *bj = NULL;
 482   int w, l;
 483
 484   w = tilec->x1-tilec->x0;
 485   l = tilec->numresolutions-1;
 486   a = tilec->data;
 487
 488   for (i = l-1; i >= stop; i--) {
 489     int rw;     /* width of the resolution level computed                       */
 490     int rh;     /* heigth of the resolution level computed                      */
 491     int rw1;    /* width of the resolution level once lower than computed one   */
 492     int rh1;    /* height of the resolution level once lower than computed one  */
 493     int cas_col;  /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 494     int cas_row;  /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 495     int dn, sn;
 496
 497     rw = tilec->resolutions[l - i].x1 - tilec->resolutions[l - i].x0;
 498     rh = tilec->resolutions[l - i].y1 - tilec->resolutions[l - i].y0;
 499     rw1= tilec->resolutions[l - i - 1].x1 - tilec->resolutions[l - i - 1].x0;
 500     rh1= tilec->resolutions[l - i - 1].y1 - tilec->resolutions[l - i - 1].y0;
 501
 502     cas_col = tilec->resolutions[l - i].x0 % 2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 503     cas_row = tilec->resolutions[l - i].y0 % 2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 504
 505     sn = rw1;
 506     dn = rw-rw1;
 507     bj = (int*)opj_malloc(rw * sizeof(int));
 508     for (j = 0; j < rh; j++) {
 509       aj = a + j * w;
 510       dwt_interleave_h(aj, bj, dn, sn, cas_col);
 511       dwt_decode_1_real(bj, dn, sn, cas_col);
 512       for (k = 0; k < rw; k++)  aj[k] = bj[k];
 513     }
 514     opj_free(bj);
 515
 516     sn = rh1;
 517     dn = rh-rh1;
 518     bj = (int*)opj_malloc(rh * sizeof(int));
 519     for (j = 0; j < rw; j++) {
 520       aj = a + j;
 521       dwt_interleave_v(aj, bj, dn, sn, w, cas_row);
 522       dwt_decode_1_real(bj, dn, sn, cas_row);
 523       for (k = 0; k < rh; k++)  aj[k * w] = bj[k];
 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   return 0;
 535 }
 536
 537 /* <summary>                */
 538 /* Get norm of 9-7 wavelet. */
 539 /* </summary>               */
 540 double dwt_getnorm_real(int level, int orient) {
 541   return dwt_norms_real[orient][level];
 542 }
 543
 544 void dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, int prec) {
 545   int numbands, bandno;
 546   numbands = 3 * tccp->numresolutions - 2;
 547   for (bandno = 0; bandno < numbands; bandno++) {
 548     double stepsize;
 549     int resno, level, orient, gain;
 550
 551     resno = (bandno == 0) ? 0 : ((bandno - 1) / 3 + 1);
 552     orient = (bandno == 0) ? 0 : ((bandno - 1) % 3 + 1);
 553     level = tccp->numresolutions - 1 - resno;
 554     gain = (tccp->qmfbid == 0) ? 0 : ((orient == 0) ? 0 : (((orient == 1) || (orient == 2)) ? 1 : 2));
 555     if (tccp->qntsty == J2K_CCP_QNTSTY_NOQNT) {
 556       stepsize = 1.0;
 557     } else {
 558       double norm = dwt_norms_real[orient][level];
 559       stepsize = (1 << (gain + 1)) / norm;
 560     }
 561     dwt_encode_stepsize((int) floor(stepsize * 8192.0), prec + gain, &tccp->stepsizes[bandno]);
 562   }
 563 }