src/gdcmopenjpeg/libopenjpeg/dwt.c

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