src/gdcmopenjpeg/libopenjpeg/dwt.c

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