mj2/libopenjpeg_097/dwt.c

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