src/lib/openjp3d/dwt.c

   1 /*
   2  * The copyright in this software is being made available under the 2-clauses
   3  * BSD License, included below. This software may be subject to other third
   4  * party and contributor rights, including patent rights, and no such rights
   5  * are granted under this license.
   6  *
   7  * Copyright (c) 2001-2003, David Janssens
   8  * Copyright (c) 2002-2003, Yannick Verschueren
   9  * Copyright (c) 2003-2005, Francois Devaux and Antonin Descampe
  10  * Copyright (c) 2005, Herve Drolon, FreeImage Team
  11  * Copyright (c) 2002-2005, Communications and remote sensing Laboratory, Universite catholique de Louvain, Belgium
  12  * Copyrigth (c) 2006, Mónica Díez, LPI-UVA, Spain
  13  * All rights reserved.
  14  *
  15  * Redistribution and use in source and binary forms, with or without
  16  * modification, are permitted provided that the following conditions
  17  * are met:
  18  * 1. Redistributions of source code must retain the above copyright
  19  *    notice, this list of conditions and the following disclaimer.
  20  * 2. Redistributions in binary form must reproduce the above copyright
  21  *    notice, this list of conditions and the following disclaimer in the
  22  *    documentation and/or other materials provided with the distribution.
  23  *
  24  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS `AS IS'
  25  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  26  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  27  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  28  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  29  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  30  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  31  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  32  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  33  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  34  * POSSIBILITY OF SUCH DAMAGE.
  35  */
  36
  37 /*
  38  *  NOTE:
  39  *  This is a modified version of the openjpeg dwt.c file.
  40  *  Average speed improvement compared to the original file (measured on
  41  *  my own machine, a P4 running at 3.0 GHz):
  42  *  5x3 wavelets about 2 times faster
  43  *  9x7 wavelets about 3 times faster
  44  *  for both, encoding and decoding.
  45  *
  46  *  The better performance is caused by doing the 1-dimensional DWT
  47  *  within a temporary buffer where the data can be accessed sequential
  48  *  for both directions, horizontal and vertical. The 2d vertical DWT was
  49  *  the major bottleneck in the former version.
  50  *
  51  *  I have also removed the "Add Patrick" part because it is not longer
  52  *  needed.
  53  *
  54  *  6/6/2005
  55  *  -Ive (aka Reiner Wahler)
  56  *  mail: ive@lilysoft.com
  57  */
  58
  59 #include "opj_includes.h"
  60
  61 /** @defgroup DWT DWT - Implementation of a discrete wavelet transform */
  62 /*@{*/
  63
  64 /** @name Local static functions */
  65 /*@{*/
  66 /*unsigned int ops;*/
  67 /**
  68 Forward lazy transform (horizontal)
  69 */
  70 static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas);
  71 /**
  72 Forward lazy transform (vertical)
  73 */
  74 static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas);
  75 /**
  76 Forward lazy transform (axial)
  77 */
  78 static void dwt_deinterleave_z(int *a, int *b, int dn, int sn, int xy, int cas);
  79 /**
  80 Inverse lazy transform (horizontal)
  81 */
  82 static void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas);
  83 /**
  84 Inverse lazy transform (vertical)
  85 */
  86 static void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas);
  87 /**
  88 Inverse lazy transform (axial)
  89 */
  90 static void dwt_interleave_z(int *a, int *b, int dn, int sn, int xy, int cas);
  91 /**
  92 Forward 5-3 wavelet transform in 1-D
  93 */
  94 static void dwt_encode_53(int *a, int dn, int sn, int cas);
  95 static void dwt_encode_97(int *a, int dn, int sn, int cas);
  96 /**
  97 Inverse 5-3 wavelet transform in 1-D
  98 */
  99 static void dwt_decode_53(int *a, int dn, int sn, int cas);
 100 static void dwt_decode_97(int *a, int dn, int sn, int cas);
 101 /**
 102 Computing of wavelet transform L2 norms for arbitrary transforms
 103 */
 104 static double dwt_calc_wtnorms(int orient, int level[3], int dwtid[3],
 105                                opj_wtfilt_t *wtfiltx, opj_wtfilt_t *wtfilty, opj_wtfilt_t *wtfiltz);
 106 /**
 107 Encoding of quantification stepsize
 108 */
 109 static void dwt_encode_stepsize(int stepsize, int numbps,
 110                                 opj_stepsize_t *bandno_stepsize);
 111 /*@}*/
 112
 113 /*@}*/
 114
 115 #define S(i) a[(i)*2]
 116 #define D(i) a[(1+(i)*2)]
 117 #define S_(i) ((i)<0?S(0):((i)>=sn?S(sn-1):S(i)))
 118 #define D_(i) ((i)<0?D(0):((i)>=dn?D(dn-1):D(i)))
 119 /* new */
 120 #define SS_(i) ((i)<0?S(0):((i)>=dn?S(dn-1):S(i)))
 121 #define DD_(i) ((i)<0?D(0):((i)>=sn?D(sn-1):D(i)))
 122
 123 /* <summary>                                                              */
 124 /* This table contains the norms of the 5-3 wavelets for different bands. */
 125 /* </summary>                                                             */
 126 static double dwt_norm[10][10][10][8];
 127 static int flagnorm[10][10][10][8];
 128
 129 /*static const double dwt_norms[5][8][10] = {
 130     {//ResZ=1
 131         {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 132         {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 133         {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 134         {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
 135     },{//ResZ=2
 136         {1.000, 1.8371, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 137         {1.2717, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 138         {1.2717, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 139         {.8803, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
 140         {1.2717},
 141         {.8803},
 142         {.8803},
 143         {.6093},
 144     },{ //ResZ=3
 145         {1.000, 1.8371, 4.5604, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 146         {1.2717, 2.6403, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 147         {1.2717, 2.6403, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 148         {.8803, 1.5286, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
 149         {1.2717, 2.6403},
 150         {.8803, 1.5286},
 151         {.8803, 1.5286},
 152         {.6093, 0.8850},
 153     },{ //ResZ=4
 154         {1.000, 1.8371, 4.5604, 12.4614, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 155         {1.2717, 2.6403, 6.7691 , 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 156         {1.2717, 2.6403, 6.7691 , 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 157         {.8803, 1.5286, 3.6770 , 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
 158         {1.2717, 2.6403, 6.7691 },
 159         {.8803, 1.5286, 3.6770 },
 160         {.8803, 1.5286, 3.6770 },
 161         {.6093, 0.8850, 1.9974 },
 162     },{ //ResZ=5
 163         {1.000, 1.8371, 4.5604, 12.4614, 34.9025, 21.34, 42.67, 85.33, 170.7, 341.3},
 164         {1.2717, 2.6403, 6.7691 , 18.6304 , 11.33, 22.64, 45.25, 90.48, 180.9},
 165         {1.2717, 2.6403, 6.7691 , 18.6304, 11.33, 22.64, 45.25, 90.48, 180.9},
 166         {.8803, 1.5286, 3.6770 , 9.9446, 6.019, 12.01, 24.00, 47.97, 95.93},
 167         {1.2717, 2.6403, 6.7691, 18.6304},
 168         {.8803, 1.5286, 3.6770, 9.9446 },
 169         {.8803, 1.5286, 3.6770, 9.9446 },
 170         {.6093, 0.8850, 1.9974, 5.3083 },
 171     }
 172 };*/
 173
 174 /* <summary>                                                              */
 175 /* This table contains the norms of the 9-7 wavelets for different bands. */
 176 /* </summary>                                                             */
 177 /*static const double dwt_norms_real[5][8][10] = {
 178     {//ResZ==1
 179         {1.000, 1.9659, 4.1224, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
 180         {1.0113, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 181         {1.0113, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 182         {0.5202, 0.9672, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722}
 183     }, { //ResZ==2
 184         {1.000, 2.7564, 4.1224, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
 185         {1.4179, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 186         {1.4179, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 187         {0.7294, 0.9672, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
 188         {1.4179},
 189         {0.7294},
 190         {0.7294},
 191         {0.3752} //HHH
 192     },{ //ResZ==3
 193         {1.000, 2.7564, 8.3700, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
 194         {1.4179, 4.0543, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 195         {1.4179, 4.0543, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 196         {0.7294, 1.9638, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
 197         {1.4179, 4.0543},
 198         {0.7294, 1.9638},
 199         {0.7294, 1.9638},
 200         {0.3752, 0.9512} //HHH
 201     },{ //ResZ==4
 202         {1.000, 2.7564, 8.3700, 24.4183, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
 203         {1.4179, 4.0543, 12.1366, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 204         {1.4179, 4.0543, 12.1366, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 205         {0.7294, 1.9638, 6.0323, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
 206         {1.4179, 4.0543, 12.1366},
 207         {0.7294, 1.9638, 6.0323},
 208         {0.7294, 1.9638, 6.0323},
 209         {0.3752, 0.9512, 2.9982} //HHH
 210     },{ //ResZ==5
 211         {1.000, 2.7564, 8.3700, 24.4183, 69.6947, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
 212         {1.4179, 4.0543, 12.1366, 35.1203, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 213         {1.4179, 4.0543, 12.1366, 35.1203, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 214         {0.7294, 1.9638, 6.0323, 17.6977, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
 215         {1.4179, 4.0543, 12.1366, 35.1203},
 216         {0.7294, 1.9638, 6.0323, 17.6977},
 217         {0.7294, 1.9638, 6.0323, 17.6977},
 218         {0.3752, 0.9512, 2.9982, 8.9182} //HHH
 219     }
 220 };*/
 221
 222 static opj_atk_t atk_info_wt[] = {
 223     {0, 1, J3D_ATK_WS, J3D_ATK_IRR, 0, J3D_ATK_WS, 1.230174104, 4, {0}, {0}, {0}, {1, 1, 1, 1}, {-1.586134342059924, -0.052980118572961, 0.882911075530934, 0.443506852043971}}, /* WT 9-7 IRR*/
 224     {1, 0, J3D_ATK_WS, J3D_ATK_REV, 0, J3D_ATK_WS, 0, 2, {0}, {1, 2}, {1, 2}, {1, 1}, {-1, 1}}, /* WT 5-3 REV*/
 225     {2, 0, J3D_ATK_ARB, J3D_ATK_REV, 0, J3D_ATK_CON, 0, 2, {0, 0}, {0, 1}, {0, 1}, {1, 1}, {{-1}, {1}}}, /* WT 2-2 REV*/
 226     {3, 0, J3D_ATK_ARB, J3D_ATK_REV, 1, J3D_ATK_CON, 0, 3, {0, 0, -1}, {0, 1, 2}, {0, 1, 2}, {1, 1, 3}, {{-1}, {1}, {1, 0, -1}}}, /* WT 2-6 REV*/
 227     {4, 0, J3D_ATK_ARB, J3D_ATK_REV, 1, J3D_ATK_CON, 0, 3, {0, 0, -2}, {0, 1, 6}, {0, 1, 32}, {1, 1, 5}, {{-1}, {1}, {-3, 22, 0, -22, 3}}}, /* WT 2-10 REV*/
 228     {5, 1, J3D_ATK_ARB, J3D_ATK_IRR, 1, J3D_ATK_WS, 1, 7, {0}, {0}, {0}, {1, 1, 2, 1, 2, 1, 3}, {{-1}, {1.58613434206}, {-0.460348209828, 0.460348209828}, {0.25}, {0.374213867768, -0.374213867768}, {-1.33613434206}, {0.29306717103, 0, -0.29306717103}}}, /* WT 6-10 IRR*/
 229     {
 230         6, 1, J3D_ATK_ARB, J3D_ATK_IRR, 0, J3D_ATK_WS, 1, 11, {0}, {0}, {0}, {1, 1, 2, 1, 2, 1, 2, 1, 2, 1, 5}, {{-1}, {0, 99715069105}, {-1.00573127827, 1.00573127827}, {-0.27040357631}, {2.20509972343, -2.20509972343}, {0.08059995736},
 231             {-1.62682532350, 1.62682532350}, {0.52040357631}, {0.60404664250, -0.60404664250}, {-0.82775064841}, {-0.06615812964, 0.29402137720, 0, -0.29402137720, 0.06615812964}
 232         }
 233     }, /* WT 10-18 IRR*/
 234     {7, 1, J3D_ATK_WS, J3D_ATK_IRR, 0, J3D_ATK_WS, 1, 2, {0}, {0}, {0}, {1, 1}, {-0.5, 0.25}},  /* WT 5-3 IRR*/
 235     {8, 0, J3D_ATK_WS, J3D_ATK_REV, 0, J3D_ATK_WS, 0, 2, {0}, {4, 4}, {8, 8}, {2, 2}, {{-9, 1}, {5, -1}}} /* WT 13-7 REV*/
 236 };
 237 /*
 238 ==========================================================
 239    local functions
 240 ==========================================================
 241 */
 242
 243 /* <summary>                             */
 244 /* Forward lazy transform (horizontal).  */
 245 /* </summary>                            */
 246 static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas)
 247 {
 248     int i;
 249     for (i = 0; i < sn; i++) {
 250         b[i] = a[2 * i + cas];
 251     }
 252     for (i = 0; i < dn; i++) {
 253         b[sn + i] = a[(2 * i + 1 - cas)];
 254     }
 255 }
 256
 257 /* <summary>                             */
 258 /* Forward lazy transform (vertical).    */
 259 /* </summary>                            */
 260 static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas)
 261 {
 262     int i;
 263     for (i = 0; i < sn; i++) {
 264         b[i * x] = a[2 * i + cas];
 265     }
 266     for (i = 0; i < dn; i++) {
 267         b[(sn + i)*x] = a[(2 * i + 1 - cas)];
 268     }
 269 }
 270
 271 /* <summary>                             */
 272 /* Forward lazy transform (axial).       */
 273 /* </summary>                            */
 274 static void dwt_deinterleave_z(int *a, int *b, int dn, int sn, int xy, int cas)
 275 {
 276     int i;
 277     for (i = 0; i < sn; i++) {
 278         b[i * xy] = a[2 * i + cas];
 279     }
 280     for (i = 0; i < dn; i++) {
 281         b[(sn + i)*xy] = a[(2 * i + 1 - cas)];
 282     }
 283 }
 284
 285 /* <summary>                             */
 286 /* Inverse lazy transform (horizontal).  */
 287 /* </summary>                            */
 288 static void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas)
 289 {
 290     int i;
 291     int *ai = NULL;
 292     int *bi = NULL;
 293     ai = a;
 294     bi = b + cas;
 295     for (i = 0; i < sn; i++) {
 296         *bi = *ai;
 297         bi += 2;
 298         ai++;
 299     }
 300     ai = a + sn;
 301     bi = b + 1 - cas;
 302     for (i = 0; i < dn; i++) {
 303         *bi = *ai;
 304         bi += 2;
 305         ai++;
 306     }
 307 }
 308
 309 /* <summary>                             */
 310 /* Inverse lazy transform (vertical).    */
 311 /* </summary>                            */
 312 static void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas)
 313 {
 314     int i;
 315     int *ai = NULL;
 316     int *bi = NULL;
 317     ai = a;
 318     bi = b + cas;
 319     for (i = 0; i < sn; i++) {
 320         *bi = *ai;
 321         bi += 2;
 322         ai += x;
 323     }
 324     ai = a + (sn * x);
 325     bi = b + 1 - cas;
 326     for (i = 0; i < dn; i++) {
 327         *bi = *ai;
 328         bi += 2;
 329         ai += x;
 330     }
 331 }
 332
 333 /* <summary>                             */
 334 /* Inverse lazy transform (axial).  */
 335 /* </summary>                            */
 336 static void dwt_interleave_z(int *a, int *b, int dn, int sn, int xy, int cas)
 337 {
 338     int i;
 339     int *ai = NULL;
 340     int *bi = NULL;
 341     ai = a;
 342     bi = b + cas;
 343     for (i = 0; i < sn; i++) {
 344         *bi = *ai;
 345         bi += 2;
 346         ai += xy;
 347     }
 348     ai = a + (sn * xy);
 349     bi = b + 1 - cas;
 350     for (i = 0; i < dn; i++) {
 351         *bi = *ai;
 352         bi += 2;
 353         ai += xy;
 354     }
 355 }
 356
 357
 358 /* <summary>                            */
 359 /* Forward 5-3 or 9-7 wavelet transform in 1-D. */
 360 /* </summary>                           */
 361 static void dwt_encode_53(int *a, int dn, int sn, int cas)
 362 {
 363     int i;
 364
 365     if (!cas) {
 366         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 367             /*for (i = 0; i < dn; i++) D(i) -= (S_(i) + S_(i + 1)) >> 1;*/
 368             /*for (i = 0; i < sn; i++) S(i) += (D_(i - 1) + D_(i) + 2) >> 2;*/
 369             for (i = 0; i < dn; i++) {
 370                 D(i) -= (S_(i) + S_(i + 1)) >> 1;
 371                 /*ops += 2;*/
 372             }
 373             for (i = 0; i < sn; i++) {
 374                 S(i) += (D_(i - 1) + D_(i) + 2) >> 2;
 375                 /*ops += 3;*/
 376             }
 377         }
 378     } else {
 379         /*if (!sn && dn == 1)
 380             S(0) *= 2;
 381         else {
 382             for (i = 0; i < dn; i++) S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
 383             for (i = 0; i < sn; i++) D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
 384         }*/
 385         if (!sn && dn == 1) {
 386             S(0) *= 2;
 387             /*ops++;*/
 388         } else {
 389             for (i = 0; i < dn; i++) {
 390                 S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
 391                 /*  ops += 2;*/
 392             }
 393             for (i = 0; i < sn; i++) {
 394                 D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
 395                 /*  ops += 3;*/
 396             }
 397         }
 398     }
 399 }
 400 static void dwt_encode_97(int *a, int dn, int sn, int cas)
 401 {
 402     int i;
 403
 404     if (!cas) {
 405         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 406             for (i = 0; i < dn; i++) {
 407                 D(i) -= fix_mul(S_(i) + S_(i + 1), 12993);
 408             }
 409             for (i = 0; i < sn; i++) {
 410                 S(i) -= fix_mul(D_(i - 1) + D_(i), 434);
 411             }
 412             for (i = 0; i < dn; i++) {
 413                 D(i) += fix_mul(S_(i) + S_(i + 1), 7233);
 414             }
 415             for (i = 0; i < sn; i++) {
 416                 S(i) += fix_mul(D_(i - 1) + D_(i), 3633);
 417             }
 418             for (i = 0; i < dn; i++) {
 419                 D(i) = fix_mul(D(i), 5038);    /*5038 */
 420             }
 421             for (i = 0; i < sn; i++) {
 422                 S(i) = fix_mul(S(i), 6659);    /*6660 */
 423             }
 424         }
 425     } else {
 426         if ((sn > 0) || (dn > 1)) { /* NEW :  CASE ONE ELEMENT */
 427             for (i = 0; i < dn; i++) {
 428                 S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993);
 429             }
 430             for (i = 0; i < sn; i++) {
 431                 D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434);
 432             }
 433             for (i = 0; i < dn; i++) {
 434                 S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233);
 435             }
 436             for (i = 0; i < sn; i++) {
 437                 D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633);
 438             }
 439             for (i = 0; i < dn; i++) {
 440                 S(i) = fix_mul(S(i), 5038);    /*5038 */
 441             }
 442             for (i = 0; i < sn; i++) {
 443                 D(i) = fix_mul(D(i), 6659);    /*6660 */
 444             }
 445         }
 446     }
 447 }
 448 /* <summary>                            */
 449 /* Inverse 5-3 or 9-7 wavelet transform in 1-D. */
 450 /* </summary>                           */
 451 static void dwt_decode_53(int *a, int dn, int sn, int cas)
 452 {
 453     int i;
 454     if (!cas) {
 455         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 456             for (i = 0; i < sn; i++) {
 457                 S(i) -= (D_(i - 1) + D_(i) + 2) >> 2;
 458             }
 459             for (i = 0; i < dn; i++) {
 460                 D(i) += (S_(i) + S_(i + 1)) >> 1;
 461             }
 462         }
 463     } else {
 464         if (!sn  && dn == 1) {        /* NEW :  CASE ONE ELEMENT */
 465             S(0) /= 2;
 466         } else {
 467             for (i = 0; i < sn; i++) {
 468                 D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2;
 469             }
 470             for (i = 0; i < dn; i++) {
 471                 S(i) += (DD_(i) + DD_(i - 1)) >> 1;
 472             }
 473         }
 474     }
 475 }
 476 static void dwt_decode_97(int *a, int dn, int sn, int cas)
 477 {
 478     int i;
 479
 480     if (!cas) {
 481         if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 482             for (i = 0; i < sn; i++) {
 483                 S(i) = fix_mul(S(i), 10078);    /* 10076 */
 484             }
 485             for (i = 0; i < dn; i++) {
 486                 D(i) = fix_mul(D(i), 13318);    /* 13320 */
 487             }
 488             for (i = 0; i < sn; i++) {
 489                 S(i) -= fix_mul(D_(i - 1) + D_(i), 3633);
 490             }
 491             for (i = 0; i < dn; i++) {
 492                 D(i) -= fix_mul(S_(i) + S_(i + 1), 7233);
 493             }
 494             for (i = 0; i < sn; i++) {
 495                 S(i) += fix_mul(D_(i - 1) + D_(i), 434);
 496             }
 497             for (i = 0; i < dn; i++) {
 498                 D(i) += fix_mul(S_(i) + S_(i + 1), 12994);    /* 12993 */
 499             }
 500         }
 501     } else {
 502         if ((sn > 0) || (dn > 1)) { /* NEW :  CASE ONE ELEMENT */
 503             for (i = 0; i < sn; i++) {
 504                 D(i) = fix_mul(D(i), 10078);    /* 10076 */
 505             }
 506             for (i = 0; i < dn; i++) {
 507                 S(i) = fix_mul(S(i), 13318);    /* 13320 */
 508             }
 509             for (i = 0; i < sn; i++) {
 510                 D(i) -= fix_mul(SS_(i) + SS_(i + 1), 3633);
 511             }
 512             for (i = 0; i < dn; i++) {
 513                 S(i) -= fix_mul(DD_(i) + DD_(i - 1), 7233);
 514             }
 515             for (i = 0; i < sn; i++) {
 516                 D(i) += fix_mul(SS_(i) + SS_(i + 1), 434);
 517             }
 518             for (i = 0; i < dn; i++) {
 519                 S(i) += fix_mul(DD_(i) + DD_(i - 1), 12994);    /* 12993 */
 520             }
 521         }
 522     }
 523 }
 524
 525
 526 /* <summary>                */
 527 /* Get norm of arbitrary wavelet transform. */
 528 /* </summary>               */
 529 static int upandconv(double *nXPS, double *LPS, int lenXPS, int lenLPS)
 530 {
 531     /* Perform the convolution of the vectors. */
 532     int i, j;
 533     double *tmp = (double *)opj_malloc(2 * lenXPS * sizeof(double));
 534     /*Upsample*/
 535     memset(tmp, 0, 2 * lenXPS * sizeof(double));
 536     for (i = 0; i < lenXPS; i++) {
 537         *(tmp + 2 * i) = *(nXPS + i);
 538         *(nXPS + i) = 0;
 539     }
 540     /*Convolution*/
 541     for (i = 0; i < 2 * lenXPS; i++) {
 542         for (j = 0; j < lenLPS; j++) {
 543             *(nXPS + i + j) = *(nXPS + i + j) + *(tmp + i) * *(LPS + j);
 544             /*fprintf(stdout,"*(tmp + %d) * *(LPS + %d) = %f * %f \n",i,j,*(tmp + i),*(LPS + j));*/
 545         }
 546     }
 547     free(tmp);
 548     return 2 * lenXPS + lenLPS - 1;
 549 }
 550
 551 static double dwt_calc_wtnorms(int orient, int level[3], int dwtid[3],
 552                                opj_wtfilt_t *wtfiltX,  opj_wtfilt_t *wtfiltY,  opj_wtfilt_t *wtfiltZ)
 553 {
 554     int i, lenLPS, lenHPS;
 555     double  Lx = 0, Ly = 0, Hx = 0, Hy = 0, Lz = 0, Hz = 0;
 556     double *nLPSx, *nHPSx, *nLPSy, *nHPSy, *nLPSz, *nHPSz;
 557     int levelx, levely, levelz;
 558
 559     levelx = (orient == 0) ? level[0] - 1 : level[0];
 560     levely = (orient == 0) ? level[1] - 1 : level[1];
 561     levelz = (orient == 0) ? level[2] - 1 : level[2];
 562
 563     /*X axis*/
 564     lenLPS = wtfiltX->lenLPS;
 565     lenHPS = wtfiltX->lenHPS;
 566     for (i = 0; i < levelx; i++) {
 567         lenLPS *= 2;
 568         lenHPS *= 2;
 569         lenLPS += wtfiltX->lenLPS - 1;
 570         lenHPS += wtfiltX->lenLPS - 1;
 571     }
 572     nLPSx = (double *)opj_malloc(lenLPS * sizeof(double));
 573     nHPSx = (double *)opj_malloc(lenHPS * sizeof(double));
 574
 575     memcpy(nLPSx, wtfiltX->LPS, wtfiltX->lenLPS * sizeof(double));
 576     memcpy(nHPSx, wtfiltX->HPS, wtfiltX->lenHPS * sizeof(double));
 577     lenLPS = wtfiltX->lenLPS;
 578     lenHPS = wtfiltX->lenHPS;
 579     for (i = 0; i < levelx; i++) {
 580         lenLPS = upandconv(nLPSx, wtfiltX->LPS, lenLPS, wtfiltX->lenLPS);
 581         lenHPS = upandconv(nHPSx, wtfiltX->LPS, lenHPS, wtfiltX->lenLPS);
 582     }
 583     for (i = 0; i < lenLPS; i++) {
 584         Lx += nLPSx[i] * nLPSx[i];
 585     }
 586     for (i = 0; i < lenHPS; i++) {
 587         Hx += nHPSx[i] * nHPSx[i];
 588     }
 589     Lx = sqrt(Lx);
 590     Hx = sqrt(Hx);
 591     free(nLPSx);
 592     free(nHPSx);
 593
 594     /*Y axis*/
 595     if (dwtid[0] != dwtid[1] || level[0] != level[1]) {
 596         lenLPS = wtfiltY->lenLPS;
 597         lenHPS = wtfiltY->lenHPS;
 598         for (i = 0; i < levely; i++) {
 599             lenLPS *= 2;
 600             lenHPS *= 2;
 601             lenLPS += wtfiltY->lenLPS - 1;
 602             lenHPS += wtfiltY->lenLPS - 1;
 603         }
 604         nLPSy = (double *)opj_malloc(lenLPS * sizeof(double));
 605         nHPSy = (double *)opj_malloc(lenHPS * sizeof(double));
 606
 607         memcpy(nLPSy, wtfiltY->LPS, wtfiltY->lenLPS * sizeof(double));
 608         memcpy(nHPSy, wtfiltY->HPS, wtfiltY->lenHPS * sizeof(double));
 609         lenLPS = wtfiltY->lenLPS;
 610         lenHPS = wtfiltY->lenHPS;
 611         for (i = 0; i < levely; i++) {
 612             lenLPS = upandconv(nLPSy, wtfiltY->LPS, lenLPS, wtfiltY->lenLPS);
 613             lenHPS = upandconv(nHPSy, wtfiltY->LPS, lenHPS, wtfiltY->lenLPS);
 614         }
 615         for (i = 0; i < lenLPS; i++) {
 616             Ly += nLPSy[i] * nLPSy[i];
 617         }
 618         for (i = 0; i < lenHPS; i++) {
 619             Hy += nHPSy[i] * nHPSy[i];
 620         }
 621         Ly = sqrt(Ly);
 622         Hy = sqrt(Hy);
 623         free(nLPSy);
 624         free(nHPSy);
 625     } else {
 626         Ly = Lx;
 627         Hy = Hx;
 628     }
 629     /*Z axis*/
 630     if (levelz >= 0) {
 631         lenLPS = wtfiltZ->lenLPS;
 632         lenHPS = wtfiltZ->lenHPS;
 633         for (i = 0; i < levelz; i++) {
 634             lenLPS *= 2;
 635             lenHPS *= 2;
 636             lenLPS += wtfiltZ->lenLPS - 1;
 637             lenHPS += wtfiltZ->lenLPS - 1;
 638         }
 639         nLPSz = (double *)opj_malloc(lenLPS * sizeof(double));
 640         nHPSz = (double *)opj_malloc(lenHPS * sizeof(double));
 641
 642         memcpy(nLPSz, wtfiltZ->LPS, wtfiltZ->lenLPS * sizeof(double));
 643         memcpy(nHPSz, wtfiltZ->HPS, wtfiltZ->lenHPS * sizeof(double));
 644         lenLPS = wtfiltZ->lenLPS;
 645         lenHPS = wtfiltZ->lenHPS;
 646         for (i = 0; i < levelz; i++) {
 647             lenLPS = upandconv(nLPSz, wtfiltZ->LPS, lenLPS, wtfiltZ->lenLPS);
 648             lenHPS = upandconv(nHPSz, wtfiltZ->LPS, lenHPS, wtfiltZ->lenLPS);
 649         }
 650         for (i = 0; i < lenLPS; i++) {
 651             Lz += nLPSz[i] * nLPSz[i];
 652         }
 653         for (i = 0; i < lenHPS; i++) {
 654             Hz += nHPSz[i] * nHPSz[i];
 655         }
 656         Lz = sqrt(Lz);
 657         Hz = sqrt(Hz);
 658         free(nLPSz);
 659         free(nHPSz);
 660     } else {
 661         Lz = 1.0;
 662         Hz = 1.0;
 663     }
 664     switch (orient) {
 665     case 0:
 666         return Lx * Ly * Lz;
 667     case 1:
 668         return Lx * Hy * Lz;
 669     case 2:
 670         return Hx * Ly * Lz;
 671     case 3:
 672         return Hx * Hy * Lz;
 673     case 4:
 674         return Lx * Ly * Hz;
 675     case 5:
 676         return Lx * Hy * Hz;
 677     case 6:
 678         return Hx * Ly * Hz;
 679     case 7:
 680         return Hx * Hy * Hz;
 681     default:
 682         return -1;
 683     }
 684
 685 }
 686 static void dwt_getwtfilters(opj_wtfilt_t *wtfilt, int dwtid)
 687 {
 688     if (dwtid == 0) { /*DWT 9-7 */
 689         wtfilt->lenLPS = 7;
 690         wtfilt->lenHPS = 9;
 691         wtfilt->LPS = (double *)opj_malloc(wtfilt->lenLPS * sizeof(double));
 692         wtfilt->HPS = (double *)opj_malloc(wtfilt->lenHPS * sizeof(double));
 693         wtfilt->LPS[0] = -0.091271763114;
 694         wtfilt->HPS[0] = 0.026748757411;
 695         wtfilt->LPS[1] = -0.057543526228;
 696         wtfilt->HPS[1] = 0.016864118443;
 697         wtfilt->LPS[2] = 0.591271763114;
 698         wtfilt->HPS[2] = -0.078223266529;
 699         wtfilt->LPS[3] = 1.115087052457;
 700         wtfilt->HPS[3] = -0.266864118443;
 701         wtfilt->LPS[4] = 0.591271763114;
 702         wtfilt->HPS[4] = 0.602949018236;
 703         wtfilt->LPS[5] = -0.057543526228;
 704         wtfilt->HPS[5] = -0.266864118443;
 705         wtfilt->LPS[6] = -0.091271763114;
 706         wtfilt->HPS[6] = -0.078223266529;
 707         wtfilt->HPS[7] = 0.016864118443;
 708         wtfilt->HPS[8] = 0.026748757411;
 709     } else if (dwtid == 1) { /*DWT 5-3 */
 710         wtfilt->lenLPS = 3;
 711         wtfilt->lenHPS = 5;
 712         wtfilt->LPS = (double *)opj_malloc(wtfilt->lenLPS * sizeof(double));
 713         wtfilt->HPS = (double *)opj_malloc(wtfilt->lenHPS * sizeof(double));
 714         wtfilt->LPS[0] = 0.5;
 715         wtfilt->HPS[0] = -0.125;
 716         wtfilt->LPS[1] = 1;
 717         wtfilt->HPS[1] = -0.25;
 718         wtfilt->LPS[2] = 0.5;
 719         wtfilt->HPS[2] = 0.75;
 720         wtfilt->HPS[3] = -0.25;
 721         wtfilt->HPS[4] = -0.125;
 722     } else {
 723         fprintf(stdout,
 724                 "[ERROR] Sorry, this wavelet hasn't been implemented so far ... Try another one :-)\n");
 725         exit(1);
 726     }
 727 }
 728 /* <summary>                            */
 729 /* Encoding of quantization stepsize for each subband. */
 730 /* </summary>                           */
 731 static void dwt_encode_stepsize(int stepsize, int numbps,
 732                                 opj_stepsize_t *bandno_stepsize)
 733 {
 734     int p, n;
 735     p = int_floorlog2(stepsize) - 13;
 736     n = 11 - int_floorlog2(stepsize);
 737     bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
 738     bandno_stepsize->expn = numbps - p;
 739     /*if J3D_CCP_QNTSTY_NOQNT --> stepsize = 8192.0 --> p = 0, n = -2 --> mant = 0; expn = (prec+gain)*/
 740     /*else --> bandno_stepsize = (1<<(numbps - expn)) + (1<<(numbps - expn - 11)) * Ub*/
 741 }
 742
 743 /*
 744 ==========================================================
 745    DWT interface
 746 ==========================================================
 747 */
 748 /* <summary>                            */
 749 /* Forward 5-3 wavelet transform in 3-D. */
 750 /* </summary>                           */
 751 void dwt_encode(opj_tcd_tilecomp_t * tilec, int dwtid[3])
 752 {
 753     int i, j, k;
 754     int x, y, z;
 755     int w, h, wh, d;
 756     int level, levelx, levely, levelz, diff;
 757     int *a = NULL;
 758     int *aj = NULL;
 759     int *bj = NULL;
 760     int *cj = NULL;
 761
 762     /*ops = 0;*/
 763
 764     memset(flagnorm, 0, 8000 * sizeof(int));
 765     w = tilec->x1 - tilec->x0;
 766     h = tilec->y1 - tilec->y0;
 767     d = tilec->z1 - tilec->z0;
 768     wh = w * h;
 769     levelx = tilec->numresolution[0] - 1;
 770     levely = tilec->numresolution[1] - 1;
 771     levelz = tilec->numresolution[2] - 1;
 772     level = int_max(levelx, int_max(levely, levelz));
 773     diff = tilec->numresolution[0] - tilec->numresolution[2];
 774
 775     a = tilec->data;
 776
 777     for (x = 0, y = 0, z = 0; (x < levelx) && (y < levely); x++, y++, z++) {
 778         int rw;         /* width of the resolution level computed                                                           */
 779         int rh;         /* heigth of the resolution level computed                                                          */
 780         int rd;         /* depth of the resolution level computed                                                          */
 781         int rw1;        /* width of the resolution level once lower than computed one                                       */
 782         int rh1;        /* height of the resolution level once lower than computed one                                      */
 783         int rd1;        /* depth of the resolution level once lower than computed one                                      */
 784         int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 785         int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 786         int cas_axl;    /* 0 = non inversion on axial filtering 1 = inversion between low-pass and high-pass filtering   */
 787         int dn, sn;
 788
 789         rw = tilec->resolutions[level - x].x1 - tilec->resolutions[level - x].x0;
 790         rh = tilec->resolutions[level - y].y1 - tilec->resolutions[level - y].y0;
 791         rd = tilec->resolutions[level - z].z1 - tilec->resolutions[level - z].z0;
 792         rw1 = tilec->resolutions[level - x - 1].x1 - tilec->resolutions[level - x -
 793                 1].x0;
 794         rh1 = tilec->resolutions[level - y - 1].y1 - tilec->resolutions[level - y -
 795                 1].y0;
 796         rd1 = tilec->resolutions[level - z - 1].z1 - tilec->resolutions[level - z -
 797                 1].z0;
 798
 799         cas_col = tilec->resolutions[level - x].x0 %
 800                   2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 801         cas_row = tilec->resolutions[level - y].y0 %
 802                   2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 803         cas_axl = tilec->resolutions[level - z].z0 % 2;
 804
 805         /*fprintf(stdout," x %d y %d z %d \n",x,y,z);
 806         fprintf(stdout," levelx %d levely %d levelz %d \n",levelx,levely,levelz);
 807         fprintf(stdout," z1 %d z0 %d\n",tilec->resolutions[level - z].z1,tilec->resolutions[level - z].z0);
 808         fprintf(stdout," rw %d rh %d rd %d \n rw1 %d rh1 %d rd1 %d \n",rw,rh,rd,rw1,rh1,rd1);*/
 809
 810         for (i = 0; i < rd; i++) {
 811
 812             cj = a + (i * wh);
 813
 814             /*Horizontal*/
 815             sn = rw1;
 816             dn = rw - rw1;
 817             bj = (int*)opj_malloc(rw * sizeof(int));
 818             if (dwtid[0] == 0) {
 819                 for (j = 0; j < rh; j++) {
 820                     aj = cj + j * w;
 821                     for (k = 0; k < rw; k++) {
 822                         bj[k] = aj[k];
 823                     }
 824                     dwt_encode_97(bj, dn, sn, cas_row);
 825                     dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 826                 }
 827             } else if (dwtid[0] == 1) {
 828                 for (j = 0; j < rh; j++) {
 829                     aj = cj + j * w;
 830                     for (k = 0; k < rw; k++) {
 831                         bj[k] = aj[k];
 832                     }
 833                     dwt_encode_53(bj, dn, sn, cas_row);
 834                     dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 835                 }
 836             }
 837             opj_free(bj);
 838
 839             /*Vertical*/
 840             sn = rh1;
 841             dn = rh - rh1;
 842             bj = (int*)opj_malloc(rh * sizeof(int));
 843             if (dwtid[1] == 0) { /*DWT 9-7*/
 844                 for (j = 0; j < rw; j++) {
 845                     aj = cj + j;
 846                     for (k = 0; k < rh; k++) {
 847                         bj[k] = aj[k * w];
 848                     }
 849                     dwt_encode_97(bj, dn, sn, cas_col);
 850                     dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 851                 }
 852             } else if (dwtid[1] == 1) { /*DWT 5-3*/
 853                 for (j = 0; j < rw; j++) {
 854                     aj = cj + j;
 855                     for (k = 0; k < rh; k++) {
 856                         bj[k] = aj[k * w];
 857                     }
 858                     dwt_encode_53(bj, dn, sn, cas_col);
 859                     dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 860                 }
 861             }
 862             opj_free(bj);
 863         }
 864
 865         if (z < levelz) {
 866             /*Axial fprintf(stdout,"Axial DWT Transform %d %d %d\n",z,rd,rd1);*/
 867             sn = rd1;
 868             dn = rd - rd1;
 869             bj = (int*)opj_malloc(rd * sizeof(int));
 870             if (dwtid[2] == 0) {
 871                 for (j = 0; j < (rw * rh); j++) {
 872                     aj = a + j;
 873                     for (k = 0; k < rd; k++) {
 874                         bj[k] = aj[k * wh];
 875                     }
 876                     dwt_encode_97(bj, dn, sn, cas_axl);
 877                     dwt_deinterleave_z(bj, aj, dn, sn, wh, cas_axl);
 878                 }
 879             } else if (dwtid[2] == 1) {
 880                 for (j = 0; j < (rw * rh); j++) {
 881                     aj = a + j;
 882                     for (k = 0; k < rd; k++) {
 883                         bj[k] = aj[k * wh];
 884                     }
 885                     dwt_encode_53(bj, dn, sn, cas_axl);
 886                     dwt_deinterleave_z(bj, aj, dn, sn, wh, cas_axl);
 887                 }
 888             }
 889             opj_free(bj);
 890         }
 891     }
 892
 893     /*fprintf(stdout,"[INFO] Ops: %d \n",ops);*/
 894 }
 895
 896
 897 /* <summary>                            */
 898 /* Inverse 5-3 wavelet transform in 3-D. */
 899 /* </summary>                           */
 900 void dwt_decode(opj_tcd_tilecomp_t * tilec, int stops[3], int dwtid[3])
 901 {
 902     int i, j, k;
 903     int x, y, z;
 904     int w, h, wh, d;
 905     int level, levelx, levely, levelz, diff;
 906     int *a = NULL;
 907     int *aj = NULL;
 908     int *bj = NULL;
 909     int *cj = NULL;
 910
 911     a = tilec->data;
 912
 913     w = tilec->x1 - tilec->x0;
 914     h = tilec->y1 - tilec->y0;
 915     d = tilec->z1 - tilec->z0;
 916     wh = w * h;
 917     levelx = tilec->numresolution[0] - 1;
 918     levely = tilec->numresolution[1] - 1;
 919     levelz = tilec->numresolution[2] - 1;
 920     level = int_max(levelx, int_max(levely, levelz));
 921     diff = tilec->numresolution[0] - tilec->numresolution[2];
 922
 923     /* General lifting framework -- DCCS-LIWT */
 924     for (x = level - 1, y = level - 1, z = level - 1; (x >= stops[0]) &&
 925             (y >= stops[1]); x--, y--, z--) {
 926         int rw;         /* width of the resolution level computed                                                           */
 927         int rh;         /* heigth of the resolution level computed                                                          */
 928         int rd;         /* depth of the resolution level computed                                                          */
 929         int rw1;        /* width of the resolution level once lower than computed one                                       */
 930         int rh1;        /* height of the resolution level once lower than computed one                                      */
 931         int rd1;        /* depth of the resolution level once lower than computed one                                      */
 932         int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 933         int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 934         int cas_axl;    /* 0 = non inversion on axial filtering 1 = inversion between low-pass and high-pass filtering   */
 935         int dn, sn;
 936
 937         rw = tilec->resolutions[level - x].x1 - tilec->resolutions[level - x].x0;
 938         rh = tilec->resolutions[level - y].y1 - tilec->resolutions[level - y].y0;
 939         rd = tilec->resolutions[level - z].z1 - tilec->resolutions[level - z].z0;
 940         rw1 = tilec->resolutions[level - x - 1].x1 - tilec->resolutions[level - x -
 941                 1].x0;
 942         rh1 = tilec->resolutions[level - y - 1].y1 - tilec->resolutions[level - y -
 943                 1].y0;
 944         rd1 = tilec->resolutions[level - z - 1].z1 - tilec->resolutions[level - z -
 945                 1].z0;
 946
 947         cas_col = tilec->resolutions[level - x].x0 %
 948                   2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 949         cas_row = tilec->resolutions[level - y].y0 %
 950                   2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 951         cas_axl = tilec->resolutions[level - z].z0 % 2;
 952
 953         /*fprintf(stdout," x %d y %d z %d \n",x,y,z);
 954         fprintf(stdout," levelx %d levely %d levelz %d \n",levelx,levely,levelz);
 955         fprintf(stdout," dwtid[0] %d [1] %d [2] %d \n",dwtid[0],dwtid[1],dwtid[2]);
 956         fprintf(stdout," rw %d rh %d rd %d \n rw1 %d rh1 %d rd1 %d \n",rw,rh,rd,rw1,rh1,rd1);
 957         fprintf(stdout,"IDWT Transform %d %d %d %d\n",level, z, rd,rd1);*/
 958
 959         if (z >= stops[2] && rd != rd1) {
 960             /*fprintf(stdout,"Axial Transform %d %d %d %d\n",levelz, z, rd,rd1);*/
 961             sn = rd1;
 962             dn = rd - rd1;
 963             bj = (int*)opj_malloc(rd * sizeof(int));
 964             if (dwtid[2] == 0) {
 965                 for (j = 0; j < (rw * rh); j++) {
 966                     aj = a + j;
 967                     dwt_interleave_z(aj, bj, dn, sn, wh, cas_axl);
 968                     dwt_decode_97(bj, dn, sn, cas_axl);
 969                     for (k = 0; k < rd; k++) {
 970                         aj[k * wh] = bj[k];
 971                     }
 972                 }
 973             } else if (dwtid[2] == 1) {
 974                 for (j = 0; j < (rw * rh); j++) {
 975                     aj = a + j;
 976                     dwt_interleave_z(aj, bj, dn, sn, wh, cas_axl);
 977                     dwt_decode_53(bj, dn, sn, cas_axl);
 978                     for (k = 0; k < rd; k++) {
 979                         aj[k * wh] = bj[k];
 980                     }
 981                 }
 982             }
 983             opj_free(bj);
 984         }
 985
 986         for (i = 0; i < rd; i++) {
 987             /*Fetch corresponding slice for doing DWT-2D*/
 988             cj = tilec->data + (i * wh);
 989
 990             /*Vertical*/
 991             sn = rh1;
 992             dn = rh - rh1;
 993             bj = (int*)opj_malloc(rh * sizeof(int));
 994             if (dwtid[1] == 0) {
 995                 for (j = 0; j < rw; j++) {
 996                     aj = cj + j;
 997                     dwt_interleave_v(aj, bj, dn, sn, w, cas_col);
 998                     dwt_decode_97(bj, dn, sn, cas_col);
 999                     for (k = 0; k < rh; k++) {
1000                         aj[k * w] = bj[k];
1001                     }
1002                 }
1003             } else if (dwtid[1] == 1) {
1004                 for (j = 0; j < rw; j++) {
1005                     aj = cj + j;
1006                     dwt_interleave_v(aj, bj, dn, sn, w, cas_col);
1007                     dwt_decode_53(bj, dn, sn, cas_col);
1008                     for (k = 0; k < rh; k++) {
1009                         aj[k * w] = bj[k];
1010                     }
1011                 }
1012             }
1013             opj_free(bj);
1014
1015             /*Horizontal*/
1016             sn = rw1;
1017             dn = rw - rw1;
1018             bj = (int*)opj_malloc(rw * sizeof(int));
1019             if (dwtid[0] == 0) {
1020                 for (j = 0; j < rh; j++) {
1021                     aj = cj + j * w;
1022                     dwt_interleave_h(aj, bj, dn, sn, cas_row);
1023                     dwt_decode_97(bj, dn, sn, cas_row);
1024                     for (k = 0; k < rw; k++) {
1025                         aj[k] = bj[k];
1026                     }
1027                 }
1028             } else if (dwtid[0] == 1) {
1029                 for (j = 0; j < rh; j++) {
1030                     aj = cj + j * w;
1031                     dwt_interleave_h(aj, bj, dn, sn, cas_row);
1032                     dwt_decode_53(bj, dn, sn, cas_row);
1033                     for (k = 0; k < rw; k++) {
1034                         aj[k] = bj[k];
1035                     }
1036                 }
1037             }
1038             opj_free(bj);
1039
1040         }
1041
1042     }
1043
1044 }
1045
1046
1047 /* <summary>                          */
1048 /* Get gain of wavelet transform. */
1049 /* </summary>                         */
1050 int dwt_getgain(int orient, int reversible)
1051 {
1052     if (reversible == 1) {
1053         if (orient == 0) {
1054             return 0;
1055         } else if (orient == 1 || orient == 2 || orient == 4) {
1056             return 1;
1057         } else if (orient == 3 || orient == 5 || orient == 6) {
1058             return 2;
1059         } else {
1060             return 3;
1061         }
1062     }
1063     /*else if (reversible == 0){*/
1064     return 0;
1065 }
1066
1067 /* <summary>                */
1068 /* Get norm of wavelet transform. */
1069 /* </summary>               */
1070 double dwt_getnorm(int orient, int level[3], int dwtid[3])
1071 {
1072     int levelx = level[0];
1073     int levely = level[1];
1074     int levelz = (level[2] < 0) ? 0 : level[2];
1075     double norm;
1076
1077     if (flagnorm[levelx][levely][levelz][orient] == 1) {
1078         norm = dwt_norm[levelx][levely][levelz][orient];
1079         /*fprintf(stdout,"[INFO] Level: %d %d %d Orient %d Dwt_norm: %f \n",level[0],level[1],level[2],orient,norm);*/
1080     } else {
1081         opj_wtfilt_t *wtfiltx = (opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
1082         opj_wtfilt_t *wtfilty = (opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
1083         opj_wtfilt_t *wtfiltz = (opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
1084         /*Fetch equivalent filters for each dimension*/
1085         dwt_getwtfilters(wtfiltx, dwtid[0]);
1086         dwt_getwtfilters(wtfilty, dwtid[1]);
1087         dwt_getwtfilters(wtfiltz, dwtid[2]);
1088         /*Calculate the corresponding norm */
1089         norm = dwt_calc_wtnorms(orient, level, dwtid, wtfiltx, wtfilty, wtfiltz);
1090         /*Save norm in array (no recalculation)*/
1091         dwt_norm[levelx][levely][levelz][orient] = norm;
1092         flagnorm[levelx][levely][levelz][orient] = 1;
1093         /*Free reserved space*/
1094         opj_free(wtfiltx->LPS);
1095         opj_free(wtfilty->LPS);
1096         opj_free(wtfiltz->LPS);
1097         opj_free(wtfiltx->HPS);
1098         opj_free(wtfilty->HPS);
1099         opj_free(wtfiltz->HPS);
1100         opj_free(wtfiltx);
1101         opj_free(wtfilty);
1102         opj_free(wtfiltz);
1103         /*fprintf(stdout,"[INFO] Dwtid: %d %d %d Level: %d %d %d Orient %d Norm: %f \n",dwtid[0],dwtid[1],dwtid[2],level[0],level[1],level[2],orient,norm);*/
1104     }
1105     return norm;
1106 }
1107 /* <summary>                                */
1108 /* Calculate explicit stepsizes for DWT.    */
1109 /* </summary>                               */
1110 void dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, int prec)
1111 {
1112     int totnumbands, bandno, diff;
1113
1114     assert(tccp->numresolution[0] >= tccp->numresolution[2]);
1115     diff = tccp->numresolution[0] -
1116            tccp->numresolution[2];     /*if RESx=RESy != RESz */
1117     totnumbands = (7 * tccp->numresolution[0] - 6) - 4 * diff; /* 3-D */
1118
1119     for (bandno = 0; bandno < totnumbands; bandno++) {
1120         double stepsize;
1121         int resno, level[3], orient, gain;
1122
1123         /* Bandno:  0 - LLL     1 - LHL
1124                     2 - HLL     3 - HHL
1125                     4 - LLH     5 - LHH
1126                     6 - HLH     7 - HHH */
1127
1128         resno = (bandno == 0) ? 0 : ((bandno <= 3 * diff) ? ((bandno - 1) / 3 + 1) : ((
1129                                          bandno + 4 * diff - 1) / 7 + 1));
1130         orient = (bandno == 0) ? 0 : ((bandno <= 3 * diff) ? ((bandno - 1) % 3 + 1) : ((
1131                                           bandno + 4 * diff - 1) % 7 + 1));
1132         level[0] = tccp->numresolution[0] - 1 - resno;
1133         level[1] = tccp->numresolution[1] - 1 - resno;
1134         level[2] = tccp->numresolution[2] - 1 - resno;
1135
1136         /* Gain:    0 - LLL     1 - LHL
1137                     1 - HLL     2 - HHL
1138                     1 - LLH     2 - LHH
1139                     2 - HLH     3 - HHH     */
1140         gain = (tccp->reversible == 0) ? 0 : ((orient == 0) ? 0 :
1141                                               (((orient == 1) || (orient == 2) || (orient == 4)) ? 1 :
1142                                                (((orient == 3) || (orient == 5) || (orient == 6)) ? 2 : 3)));
1143
1144         if (tccp->qntsty == J3D_CCP_QNTSTY_NOQNT) {
1145             stepsize = 1.0;
1146         } else {
1147             double norm = dwt_getnorm(orient, level,
1148                                       tccp->dwtid); /*Fetch norms if irreversible transform (by the moment only I9.7)*/
1149             stepsize = (1 << (gain + 1)) / norm;
1150         }
1151         /*fprintf(stdout,"[INFO] Bandno: %d Orient: %d Level: %d %d %d Stepsize: %f\n",bandno,orient,level[0],level[1],level[2],stepsize);*/
1152         dwt_encode_stepsize((int) floor(stepsize * 8192.0), prec + gain,
1153                             &tccp->stepsizes[bandno]);
1154     }
1155 }
1156
1157
1158