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 tranform 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 tranform 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], opj_wtfilt_t *wtfiltx, opj_wtfilt_t *wtfilty, opj_wtfilt_t *wtfiltz);
 105 /**
 106 Encoding of quantification stepsize
 107 */
 108 static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno_stepsize);
 109 /*@}*/
 110
 111 /*@}*/
 112
 113 #define S(i) a[(i)*2]
 114 #define D(i) a[(1+(i)*2)]
 115 #define S_(i) ((i)<0?S(0):((i)>=sn?S(sn-1):S(i)))
 116 #define D_(i) ((i)<0?D(0):((i)>=dn?D(dn-1):D(i)))
 117 /* new */
 118 #define SS_(i) ((i)<0?S(0):((i)>=dn?S(dn-1):S(i)))
 119 #define DD_(i) ((i)<0?D(0):((i)>=sn?D(sn-1):D(i)))
 120
 121 /* <summary>                                                              */
 122 /* This table contains the norms of the 5-3 wavelets for different bands. */
 123 /* </summary>                                                             */
 124 static double dwt_norm[10][10][10][8];
 125 static int flagnorm[10][10][10][8];
 126
 127 /*static const double dwt_norms[5][8][10] = {
 128         {//ResZ=1
 129                 {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 130                 {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 131                 {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 132                 {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
 133         },{//ResZ=2
 134                 {1.000, 1.8371, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 135                 {1.2717, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 136                 {1.2717, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 137                 {.8803, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
 138                 {1.2717},
 139                 {.8803},
 140                 {.8803},
 141                 {.6093},
 142         },{ //ResZ=3
 143                 {1.000, 1.8371, 4.5604, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 144                 {1.2717, 2.6403, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 145                 {1.2717, 2.6403, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 146                 {.8803, 1.5286, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
 147                 {1.2717, 2.6403},
 148                 {.8803, 1.5286},
 149                 {.8803, 1.5286},
 150                 {.6093, 0.8850},
 151         },{ //ResZ=4
 152                 {1.000, 1.8371, 4.5604, 12.4614, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
 153                 {1.2717, 2.6403, 6.7691 , 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 154                 {1.2717, 2.6403, 6.7691 , 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
 155                 {.8803, 1.5286, 3.6770 , 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
 156                 {1.2717, 2.6403, 6.7691 },
 157                 {.8803, 1.5286, 3.6770 },
 158                 {.8803, 1.5286, 3.6770 },
 159                 {.6093, 0.8850, 1.9974 },
 160         },{ //ResZ=5
 161                 {1.000, 1.8371, 4.5604, 12.4614, 34.9025, 21.34, 42.67, 85.33, 170.7, 341.3},
 162                 {1.2717, 2.6403, 6.7691 , 18.6304 , 11.33, 22.64, 45.25, 90.48, 180.9},
 163                 {1.2717, 2.6403, 6.7691 , 18.6304, 11.33, 22.64, 45.25, 90.48, 180.9},
 164                 {.8803, 1.5286, 3.6770 , 9.9446, 6.019, 12.01, 24.00, 47.97, 95.93},
 165                 {1.2717, 2.6403, 6.7691, 18.6304},
 166                 {.8803, 1.5286, 3.6770, 9.9446 },
 167                 {.8803, 1.5286, 3.6770, 9.9446 },
 168                 {.6093, 0.8850, 1.9974, 5.3083 },
 169         }
 170 };*/
 171
 172 /* <summary>                                                              */
 173 /* This table contains the norms of the 9-7 wavelets for different bands. */
 174 /* </summary>                                                             */
 175 /*static const double dwt_norms_real[5][8][10] = {
 176         {//ResZ==1
 177                 {1.000, 1.9659, 4.1224, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
 178                 {1.0113, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 179                 {1.0113, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 180                 {0.5202, 0.9672, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722}
 181         }, { //ResZ==2
 182                 {1.000, 2.7564, 4.1224, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
 183                 {1.4179, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 184                 {1.4179, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 185                 {0.7294, 0.9672, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
 186                 {1.4179},
 187                 {0.7294},
 188                 {0.7294},
 189                 {0.3752} //HHH
 190         },{ //ResZ==3
 191                 {1.000, 2.7564, 8.3700, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
 192                 {1.4179, 4.0543, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 193                 {1.4179, 4.0543, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 194                 {0.7294, 1.9638, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
 195                 {1.4179, 4.0543},
 196                 {0.7294, 1.9638},
 197                 {0.7294, 1.9638},
 198                 {0.3752, 0.9512} //HHH
 199         },{ //ResZ==4
 200                 {1.000, 2.7564, 8.3700, 24.4183, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
 201                 {1.4179, 4.0543, 12.1366, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 202                 {1.4179, 4.0543, 12.1366, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 203                 {0.7294, 1.9638, 6.0323, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
 204                 {1.4179, 4.0543, 12.1366},
 205                 {0.7294, 1.9638, 6.0323},
 206                 {0.7294, 1.9638, 6.0323},
 207                 {0.3752, 0.9512, 2.9982} //HHH
 208         },{ //ResZ==5
 209                 {1.000, 2.7564, 8.3700, 24.4183, 69.6947, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
 210                 {1.4179, 4.0543, 12.1366, 35.1203, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 211                 {1.4179, 4.0543, 12.1366, 35.1203, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
 212                 {0.7294, 1.9638, 6.0323, 17.6977, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
 213                 {1.4179, 4.0543, 12.1366, 35.1203},
 214                 {0.7294, 1.9638, 6.0323, 17.6977},
 215                 {0.7294, 1.9638, 6.0323, 17.6977},
 216                 {0.3752, 0.9512, 2.9982, 8.9182} //HHH
 217         }
 218 };*/
 219
 220 static opj_atk_t atk_info_wt[] = {
 221         {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*/
 222         {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*/
 223         {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*/
 224         {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*/
 225         {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*/
 226         {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*/
 227         {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},
 228                 {-1.62682532350, 1.62682532350},{0.52040357631},{0.60404664250, -0.60404664250},{-0.82775064841},{-0.06615812964, 0.29402137720, 0, -0.29402137720, 0.06615812964}}}, /* WT 10-18 IRR*/
 229         {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*/
 230         {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*/
 231 };
 232 /*
 233 ==========================================================
 234    local functions
 235 ==========================================================
 236 */
 237
 238 /* <summary>                                     */
 239 /* Forward lazy transform (horizontal).  */
 240 /* </summary>                            */
 241 static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas) {
 242         int i;
 243     for (i=0; i<sn; i++) b[i]=a[2*i+cas];
 244     for (i=0; i<dn; i++) b[sn+i]=a[(2*i+1-cas)];
 245 }
 246
 247 /* <summary>                             */
 248 /* Forward lazy transform (vertical).    */
 249 /* </summary>                            */
 250 static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
 251     int i;
 252     for (i=0; i<sn; i++) b[i*x]=a[2*i+cas];
 253     for (i=0; i<dn; i++) b[(sn+i)*x]=a[(2*i+1-cas)];
 254 }
 255
 256 /* <summary>                             */
 257 /* Forward lazy transform (axial).       */
 258 /* </summary>                            */
 259 static void dwt_deinterleave_z(int *a, int *b, int dn, int sn, int xy, int cas) {
 260     int i;
 261     for (i=0; i<sn; i++) b[i*xy]=a[2*i+cas];
 262     for (i=0; i<dn; i++) b[(sn+i)*xy]=a[(2*i+1-cas)];
 263 }
 264
 265 /* <summary>                             */
 266 /* Inverse lazy transform (horizontal).  */
 267 /* </summary>                            */
 268 static void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas) {
 269     int i;
 270     int *ai = NULL;
 271     int *bi = NULL;
 272     ai = a;
 273     bi = b + cas;
 274     for (i = 0; i < sn; i++) {
 275       *bi = *ai;
 276           bi += 2;
 277           ai++;
 278     }
 279     ai = a + sn;
 280     bi = b + 1 - cas;
 281     for (i = 0; i < dn; i++) {
 282       *bi = *ai;
 283           bi += 2;
 284           ai++;
 285     }
 286 }
 287
 288 /* <summary>                             */
 289 /* Inverse lazy transform (vertical).    */
 290 /* </summary>                            */
 291 static void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
 292     int i;
 293     int *ai = NULL;
 294     int *bi = NULL;
 295     ai = a;
 296     bi = b + cas;
 297     for (i = 0; i < sn; i++) {
 298       *bi = *ai;
 299           bi += 2;
 300           ai += x;
 301     }
 302     ai = a + (sn * x);
 303     bi = b + 1 - cas;
 304     for (i = 0; i < dn; i++) {
 305       *bi = *ai;
 306           bi += 2;
 307           ai += x;
 308     }
 309 }
 310
 311 /* <summary>                             */
 312 /* Inverse lazy transform (axial).  */
 313 /* </summary>                            */
 314 static void dwt_interleave_z(int *a, int *b, int dn, int sn, int xy, int cas) {
 315     int i;
 316     int *ai = NULL;
 317     int *bi = NULL;
 318     ai = a;
 319     bi = b + cas;
 320     for (i = 0; i < sn; i++) {
 321       *bi = *ai;
 322           bi += 2;
 323           ai += xy;
 324     }
 325     ai = a + (sn * xy);
 326     bi = b + 1 - cas;
 327     for (i = 0; i < dn; i++) {
 328       *bi = *ai;
 329           bi += 2;
 330           ai += xy;
 331     }
 332 }
 333
 334
 335 /* <summary>                            */
 336 /* Forward 5-3 or 9-7 wavelet tranform in 1-D. */
 337 /* </summary>                           */
 338 static void dwt_encode_53(int *a, int dn, int sn, int cas) {
 339         int i;
 340
 341         if (!cas) {
 342                 if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 343                         /*for (i = 0; i < dn; i++) D(i) -= (S_(i) + S_(i + 1)) >> 1;*/
 344                         /*for (i = 0; i < sn; i++) S(i) += (D_(i - 1) + D_(i) + 2) >> 2;*/
 345                         for (i = 0; i < dn; i++){
 346                                 D(i) -= (S_(i) + S_(i + 1)) >> 1;
 347                                 /*ops += 2;*/
 348                         }
 349                         for (i = 0; i < sn; i++){
 350                                 S(i) += (D_(i - 1) + D_(i) + 2) >> 2;
 351                                 /*ops += 3;*/
 352                         }
 353                 }
 354         } else {
 355                 /*if (!sn && dn == 1)
 356                         S(0) *= 2;
 357                 else {
 358                         for (i = 0; i < dn; i++) S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
 359                         for (i = 0; i < sn; i++) D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
 360                 }*/
 361                 if (!sn && dn == 1){
 362                         S(0) *= 2;
 363                         /*ops++;*/
 364                 } else {
 365                         for (i = 0; i < dn; i++){
 366                                 S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
 367                         /*      ops += 2;*/
 368                         }
 369                         for (i = 0; i < sn; i++){
 370                                 D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
 371                         /*      ops += 3;*/
 372                         }
 373                 }
 374         }
 375 }
 376 static void dwt_encode_97(int *a, int dn, int sn, int cas) {
 377         int i;
 378
 379         if (!cas) {
 380                         if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 381                                 for (i = 0; i < dn; i++)
 382                                         D(i) -= fix_mul(S_(i) + S_(i + 1), 12993);
 383                                 for (i = 0; i < sn; i++)
 384                                         S(i) -= fix_mul(D_(i - 1) + D_(i), 434);
 385                                 for (i = 0; i < dn; i++)
 386                                         D(i) += fix_mul(S_(i) + S_(i + 1), 7233);
 387                                 for (i = 0; i < sn; i++)
 388                                         S(i) += fix_mul(D_(i - 1) + D_(i), 3633);
 389                                 for (i = 0; i < dn; i++)
 390                                         D(i) = fix_mul(D(i), 5038);     /*5038 */
 391                                 for (i = 0; i < sn; i++)
 392                                         S(i) = fix_mul(S(i), 6659);     /*6660 */
 393                         }
 394                 } else {
 395                         if ((sn > 0) || (dn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 396                                 for (i = 0; i < dn; i++)
 397                                         S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993);
 398                                 for (i = 0; i < sn; i++)
 399                                         D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434);
 400                                 for (i = 0; i < dn; i++)
 401                                         S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233);
 402                                 for (i = 0; i < sn; i++)
 403                                         D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633);
 404                                 for (i = 0; i < dn; i++)
 405                                         S(i) = fix_mul(S(i), 5038);     /*5038 */
 406                                 for (i = 0; i < sn; i++)
 407                                         D(i) = fix_mul(D(i), 6659);     /*6660 */
 408                         }
 409                 }
 410 }
 411 /* <summary>                            */
 412 /* Inverse 5-3 or 9-7 wavelet tranform in 1-D. */
 413 /* </summary>                           */
 414 static void dwt_decode_53(int *a, int dn, int sn, int cas) {
 415         int i;
 416         if (!cas) {
 417                 if ((dn > 0) || (sn > 1)) { /* NEW :  CASE ONE ELEMENT */
 418                         for (i = 0; i < sn; i++) S(i) -= (D_(i - 1) + D_(i) + 2) >> 2;
 419                         for (i = 0; i < dn; i++) D(i) += (S_(i) + S_(i + 1)) >> 1;
 420                 }
 421         } else {
 422                 if (!sn  && dn == 1)          /* NEW :  CASE ONE ELEMENT */
 423                         S(0) /= 2;
 424                 else {
 425                         for (i = 0; i < sn; i++) D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2;
 426                         for (i = 0; i < dn; i++) S(i) += (DD_(i) + DD_(i - 1)) >> 1;
 427                 }
 428         }
 429 }
 430 static void dwt_decode_97(int *a, int dn, int sn, int cas) {
 431         int i;
 432
 433         if (!cas) {
 434                 if ((dn > 0) || (sn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 435                         for (i = 0; i < sn; i++)
 436                                 S(i) = fix_mul(S(i), 10078);    /* 10076 */
 437                         for (i = 0; i < dn; i++)
 438                                 D(i) = fix_mul(D(i), 13318);    /* 13320 */
 439                         for (i = 0; i < sn; i++)
 440                                 S(i) -= fix_mul(D_(i - 1) + D_(i), 3633);
 441                         for (i = 0; i < dn; i++)
 442                                 D(i) -= fix_mul(S_(i) + S_(i + 1), 7233);
 443                         for (i = 0; i < sn; i++)
 444                                 S(i) += fix_mul(D_(i - 1) + D_(i), 434);
 445                         for (i = 0; i < dn; i++)
 446                                 D(i) += fix_mul(S_(i) + S_(i + 1), 12994);      /* 12993 */
 447                 }
 448         } else {
 449                 if ((sn > 0) || (dn > 1)) {     /* NEW :  CASE ONE ELEMENT */
 450                         for (i = 0; i < sn; i++)
 451                                 D(i) = fix_mul(D(i), 10078);    /* 10076 */
 452                         for (i = 0; i < dn; i++)
 453                                 S(i) = fix_mul(S(i), 13318);    /* 13320 */
 454                         for (i = 0; i < sn; i++)
 455                                 D(i) -= fix_mul(SS_(i) + SS_(i + 1), 3633);
 456                         for (i = 0; i < dn; i++)
 457                                 S(i) -= fix_mul(DD_(i) + DD_(i - 1), 7233);
 458                         for (i = 0; i < sn; i++)
 459                                 D(i) += fix_mul(SS_(i) + SS_(i + 1), 434);
 460                         for (i = 0; i < dn; i++)
 461                                 S(i) += fix_mul(DD_(i) + DD_(i - 1), 12994);    /* 12993 */
 462                 }
 463         }
 464 }
 465
 466
 467 /* <summary>                */
 468 /* Get norm of arbitrary wavelet transform. */
 469 /* </summary>               */
 470 static int upandconv(double *nXPS, double *LPS, int lenXPS, int lenLPS) {
 471         /* Perform the convolution of the vectors. */
 472         int i,j;
 473         double *tmp = (double *)opj_malloc(2*lenXPS * sizeof(double));
 474         /*Upsample*/
 475         memset(tmp, 0, 2*lenXPS*sizeof(double));
 476         for (i = 0; i < lenXPS; i++) {
 477                 *(tmp + 2*i) = *(nXPS + i);
 478                 *(nXPS + i) = 0;
 479         }
 480         /*Convolution*/
 481         for (i = 0; i < 2*lenXPS; i++) {
 482                 for (j = 0; j < lenLPS; j++) {
 483                         *(nXPS+i+j) = *(nXPS+i+j) + *(tmp + i) * *(LPS + j);
 484                         /*fprintf(stdout,"*(tmp + %d) * *(LPS + %d) = %f * %f \n",i,j,*(tmp + i),*(LPS + j));*/
 485                 }
 486         }
 487         free(tmp);
 488         return 2*lenXPS+lenLPS-1;
 489 }
 490
 491 static double dwt_calc_wtnorms(int orient, int level[3], int dwtid[3], opj_wtfilt_t *wtfiltX,  opj_wtfilt_t *wtfiltY,  opj_wtfilt_t *wtfiltZ) {
 492         int i, lenLPS, lenHPS;
 493         double  Lx = 0, Ly= 0, Hx= 0, Hy= 0, Lz= 0, Hz= 0;
 494         double *nLPSx, *nHPSx,*nLPSy, *nHPSy,*nLPSz, *nHPSz;
 495         int levelx, levely, levelz;
 496
 497         levelx = (orient == 0) ? level[0]-1 : level[0];
 498         levely = (orient == 0) ? level[1]-1 : level[1];
 499         levelz = (orient == 0) ? level[2]-1 : level[2];
 500
 501         /*X axis*/
 502         lenLPS = wtfiltX->lenLPS;
 503         lenHPS = wtfiltX->lenHPS;
 504         for (i = 0; i < levelx; i++) {
 505                 lenLPS *= 2;
 506                 lenHPS *= 2;
 507                 lenLPS += wtfiltX->lenLPS - 1;
 508                 lenHPS += wtfiltX->lenLPS - 1;
 509         }
 510         nLPSx = (double *)opj_malloc(lenLPS * sizeof(double));
 511         nHPSx = (double *)opj_malloc(lenHPS * sizeof(double));
 512
 513         memcpy(nLPSx, wtfiltX->LPS, wtfiltX->lenLPS * sizeof(double));
 514         memcpy(nHPSx, wtfiltX->HPS, wtfiltX->lenHPS * sizeof(double));
 515         lenLPS = wtfiltX->lenLPS;
 516         lenHPS = wtfiltX->lenHPS;
 517         for (i = 0; i < levelx; i++) {
 518                 lenLPS = upandconv(nLPSx, wtfiltX->LPS, lenLPS, wtfiltX->lenLPS);
 519                 lenHPS = upandconv(nHPSx, wtfiltX->LPS, lenHPS, wtfiltX->lenLPS);
 520         }
 521         for (i = 0; i < lenLPS; i++)
 522                 Lx += nLPSx[i] * nLPSx[i];
 523         for (i = 0; i < lenHPS; i++)
 524                 Hx += nHPSx[i] * nHPSx[i];
 525         Lx = sqrt(Lx);
 526         Hx = sqrt(Hx);
 527         free(nLPSx);
 528         free(nHPSx);
 529
 530         /*Y axis*/
 531         if (dwtid[0] != dwtid[1] || level[0] != level[1]){
 532                 lenLPS = wtfiltY->lenLPS;
 533                 lenHPS = wtfiltY->lenHPS;
 534                 for (i = 0; i < levely; i++) {
 535                         lenLPS *= 2;
 536                         lenHPS *= 2;
 537                         lenLPS += wtfiltY->lenLPS - 1;
 538                         lenHPS += wtfiltY->lenLPS - 1;
 539                 }
 540                 nLPSy = (double *)opj_malloc(lenLPS * sizeof(double));
 541                 nHPSy = (double *)opj_malloc(lenHPS * sizeof(double));
 542
 543                 memcpy(nLPSy, wtfiltY->LPS, wtfiltY->lenLPS * sizeof(double));
 544                 memcpy(nHPSy, wtfiltY->HPS, wtfiltY->lenHPS * sizeof(double));
 545                 lenLPS = wtfiltY->lenLPS;
 546                 lenHPS = wtfiltY->lenHPS;
 547                 for (i = 0; i < levely; i++) {
 548                         lenLPS = upandconv(nLPSy, wtfiltY->LPS, lenLPS, wtfiltY->lenLPS);
 549                         lenHPS = upandconv(nHPSy, wtfiltY->LPS, lenHPS, wtfiltY->lenLPS);
 550                 }
 551                 for (i = 0; i < lenLPS; i++)
 552                         Ly += nLPSy[i] * nLPSy[i];
 553                 for (i = 0; i < lenHPS; i++)
 554                         Hy += nHPSy[i] * nHPSy[i];
 555                 Ly = sqrt(Ly);
 556                 Hy = sqrt(Hy);
 557                 free(nLPSy);
 558                 free(nHPSy);
 559         } else {
 560                 Ly = Lx;
 561                 Hy = Hx;
 562         }
 563         /*Z axis*/
 564         if (levelz >= 0) {
 565                 lenLPS = wtfiltZ->lenLPS;
 566                 lenHPS = wtfiltZ->lenHPS;
 567                 for (i = 0; i < levelz; i++) {
 568                         lenLPS *= 2;
 569                         lenHPS *= 2;
 570                         lenLPS += wtfiltZ->lenLPS - 1;
 571                         lenHPS += wtfiltZ->lenLPS - 1;
 572                 }
 573                 nLPSz = (double *)opj_malloc(lenLPS * sizeof(double));
 574                 nHPSz = (double *)opj_malloc(lenHPS * sizeof(double));
 575
 576                 memcpy(nLPSz, wtfiltZ->LPS, wtfiltZ->lenLPS * sizeof(double));
 577                 memcpy(nHPSz, wtfiltZ->HPS, wtfiltZ->lenHPS * sizeof(double));
 578                 lenLPS = wtfiltZ->lenLPS;
 579                 lenHPS = wtfiltZ->lenHPS;
 580                 for (i = 0; i < levelz; i++) {
 581                         lenLPS = upandconv(nLPSz, wtfiltZ->LPS, lenLPS, wtfiltZ->lenLPS);
 582                         lenHPS = upandconv(nHPSz, wtfiltZ->LPS, lenHPS, wtfiltZ->lenLPS);
 583                 }
 584                 for (i = 0; i < lenLPS; i++)
 585                         Lz += nLPSz[i] * nLPSz[i];
 586                 for (i = 0; i < lenHPS; i++)
 587                         Hz += nHPSz[i] * nHPSz[i];
 588                 Lz = sqrt(Lz);
 589                 Hz = sqrt(Hz);
 590                 free(nLPSz);
 591                 free(nHPSz);
 592         } else {
 593                 Lz = 1.0; Hz = 1.0;
 594         }
 595         switch (orient) {
 596                 case 0:
 597                         return Lx * Ly * Lz;
 598                 case 1:
 599                         return Lx * Hy * Lz;
 600                 case 2:
 601                         return Hx * Ly * Lz;
 602                 case 3:
 603                         return Hx * Hy * Lz;
 604                 case 4:
 605                         return Lx * Ly * Hz;
 606                 case 5:
 607                         return Lx * Hy * Hz;
 608                 case 6:
 609                         return Hx * Ly * Hz;
 610                 case 7:
 611                         return Hx * Hy * Hz;
 612                 default:
 613                         return -1;
 614         }
 615
 616 }
 617 static void dwt_getwtfilters(opj_wtfilt_t *wtfilt, int dwtid) {
 618         if (dwtid == 0) { /*DWT 9-7 */
 619                         wtfilt->lenLPS = 7;             wtfilt->lenHPS = 9;
 620                         wtfilt->LPS = (double *)opj_malloc(wtfilt->lenLPS * sizeof(double));
 621                         wtfilt->HPS = (double *)opj_malloc(wtfilt->lenHPS * sizeof(double));
 622                         wtfilt->LPS[0] = -0.091271763114;       wtfilt->HPS[0] = 0.026748757411;
 623                         wtfilt->LPS[1] = -0.057543526228;       wtfilt->HPS[1] = 0.016864118443;
 624                         wtfilt->LPS[2] = 0.591271763114;        wtfilt->HPS[2] = -0.078223266529;
 625                         wtfilt->LPS[3] = 1.115087052457;        wtfilt->HPS[3] = -0.266864118443;
 626                         wtfilt->LPS[4] = 0.591271763114;        wtfilt->HPS[4] = 0.602949018236;
 627                         wtfilt->LPS[5] = -0.057543526228;       wtfilt->HPS[5] = -0.266864118443;
 628                         wtfilt->LPS[6] = -0.091271763114;       wtfilt->HPS[6] = -0.078223266529;
 629                                                                                                 wtfilt->HPS[7] = 0.016864118443;
 630                                                                                                 wtfilt->HPS[8] = 0.026748757411;
 631         } else if (dwtid == 1) { /*DWT 5-3 */
 632                         wtfilt->lenLPS = 3;             wtfilt->lenHPS = 5;
 633                         wtfilt->LPS = (double *)opj_malloc(wtfilt->lenLPS * sizeof(double));
 634                         wtfilt->HPS = (double *)opj_malloc(wtfilt->lenHPS * sizeof(double));
 635                         wtfilt->LPS[0] = 0.5;   wtfilt->HPS[0] = -0.125;
 636                         wtfilt->LPS[1] = 1;             wtfilt->HPS[1] = -0.25;
 637                         wtfilt->LPS[2] = 0.5;   wtfilt->HPS[2] = 0.75;
 638                                                                         wtfilt->HPS[3] = -0.25;
 639                                                                         wtfilt->HPS[4] = -0.125;
 640         } else {
 641                 fprintf(stdout,"[ERROR] Sorry, this wavelet hasn't been implemented so far ... Try another one :-)\n");
 642                 exit(1);
 643         }
 644 }
 645 /* <summary>                            */
 646 /* Encoding of quantization stepsize for each subband. */
 647 /* </summary>                           */
 648 static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno_stepsize) {
 649         int p, n;
 650         p = int_floorlog2(stepsize) - 13;
 651         n = 11 - int_floorlog2(stepsize);
 652         bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
 653         bandno_stepsize->expn = numbps - p;
 654         /*if J3D_CCP_QNTSTY_NOQNT --> stepsize = 8192.0 --> p = 0, n = -2 --> mant = 0; expn = (prec+gain)*/
 655         /*else --> bandno_stepsize = (1<<(numbps - expn)) + (1<<(numbps - expn - 11)) * Ub*/
 656 }
 657
 658 /*
 659 ==========================================================
 660    DWT interface
 661 ==========================================================
 662 */
 663 /* <summary>                            */
 664 /* Forward 5-3 wavelet tranform in 3-D. */
 665 /* </summary>                           */
 666 void dwt_encode(opj_tcd_tilecomp_t * tilec, int dwtid[3]) {
 667         int i, j, k;
 668         int x, y, z;
 669         int w, h, wh, d;
 670         int level,levelx,levely,levelz,diff;
 671         int *a = NULL;
 672         int *aj = NULL;
 673         int *bj = NULL;
 674         int *cj = NULL;
 675
 676         /*ops = 0;*/
 677
 678         memset(flagnorm,0,8000*sizeof(int));
 679         w = tilec->x1-tilec->x0;
 680         h = tilec->y1-tilec->y0;
 681         d = tilec->z1-tilec->z0;
 682         wh = w * h;
 683         levelx = tilec->numresolution[0]-1;
 684         levely = tilec->numresolution[1]-1;
 685         levelz = tilec->numresolution[2]-1;
 686         level = int_max(levelx,int_max(levely,levelz));
 687         diff = tilec->numresolution[0] - tilec->numresolution[2];
 688
 689         a = tilec->data;
 690
 691         for (x = 0, y = 0, z = 0; (x < levelx) && (y < levely); x++, y++, z++) {
 692                 int rw;                 /* width of the resolution level computed                                                           */
 693                 int rh;                 /* heigth of the resolution level computed                                                          */
 694                 int rd;                 /* depth of the resolution level computed                                                          */
 695                 int rw1;                /* width of the resolution level once lower than computed one                                       */
 696                 int rh1;                /* height of the resolution level once lower than computed one                                      */
 697                 int rd1;                /* depth of the resolution level once lower than computed one                                      */
 698                 int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 699                 int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 700                 int cas_axl;    /* 0 = non inversion on axial filtering 1 = inversion between low-pass and high-pass filtering   */
 701                 int dn, sn;
 702
 703                 rw = tilec->resolutions[level - x].x1 - tilec->resolutions[level - x].x0;
 704                 rh = tilec->resolutions[level - y].y1 - tilec->resolutions[level - y].y0;
 705                 rd = tilec->resolutions[level - z].z1 - tilec->resolutions[level - z].z0;
 706                 rw1= tilec->resolutions[level - x - 1].x1 - tilec->resolutions[level - x - 1].x0;
 707                 rh1= tilec->resolutions[level - y - 1].y1 - tilec->resolutions[level - y - 1].y0;
 708                 rd1= tilec->resolutions[level - z - 1].z1 - tilec->resolutions[level - z - 1].z0;
 709
 710                 cas_col = tilec->resolutions[level - x].x0 % 2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 711                 cas_row = tilec->resolutions[level - y].y0 % 2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 712                 cas_axl = tilec->resolutions[level - z].z0 % 2;
 713
 714                 /*fprintf(stdout," x %d y %d z %d \n",x,y,z);
 715                 fprintf(stdout," levelx %d levely %d levelz %d \n",levelx,levely,levelz);
 716                 fprintf(stdout," z1 %d z0 %d\n",tilec->resolutions[level - z].z1,tilec->resolutions[level - z].z0);
 717                 fprintf(stdout," rw %d rh %d rd %d \n rw1 %d rh1 %d rd1 %d \n",rw,rh,rd,rw1,rh1,rd1);*/
 718
 719                 for (i = 0; i < rd; i++) {
 720
 721                         cj = a + (i * wh);
 722
 723                         /*Horizontal*/
 724                         sn = rw1;
 725                         dn = rw - rw1;
 726                         bj = (int*)opj_malloc(rw * sizeof(int));
 727                         if (dwtid[0] == 0) {
 728                                 for (j = 0; j < rh; j++) {
 729                                         aj = cj + j * w;
 730                                         for (k = 0; k < rw; k++)  bj[k] = aj[k];
 731                                         dwt_encode_97(bj, dn, sn, cas_row);
 732                                         dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 733                                 }
 734                         } else if (dwtid[0] == 1) {
 735                                 for (j = 0; j < rh; j++) {
 736                                         aj = cj + j * w;
 737                                         for (k = 0; k < rw; k++)  bj[k] = aj[k];
 738                                         dwt_encode_53(bj, dn, sn, cas_row);
 739                                         dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
 740                                 }
 741                         }
 742                         opj_free(bj);
 743
 744                         /*Vertical*/
 745                         sn = rh1;
 746                         dn = rh - rh1;
 747                         bj = (int*)opj_malloc(rh * sizeof(int));
 748                         if (dwtid[1] == 0) { /*DWT 9-7*/
 749                                 for (j = 0; j < rw; j++) {
 750                                         aj = cj + j;
 751                                         for (k = 0; k < rh; k++)  bj[k] = aj[k*w];
 752                                         dwt_encode_97(bj, dn, sn, cas_col);
 753                                         dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 754                                 }
 755             } else if (dwtid[1] == 1) { /*DWT 5-3*/
 756                                 for (j = 0; j < rw; j++) {
 757                                         aj = cj + j;
 758                                         for (k = 0; k < rh; k++)  bj[k] = aj[k*w];
 759                                         dwt_encode_53(bj, dn, sn, cas_col);
 760                                         dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
 761                                 }
 762                         }
 763                         opj_free(bj);
 764                 }
 765
 766                 if (z < levelz){
 767                         /*Axial fprintf(stdout,"Axial DWT Transform %d %d %d\n",z,rd,rd1);*/
 768                         sn = rd1;
 769                         dn = rd - rd1;
 770                         bj = (int*)opj_malloc(rd * sizeof(int));
 771                         if (dwtid[2] == 0) {
 772                 for (j = 0; j < (rw*rh); j++) {
 773                                         aj = a + j;
 774                                         for (k = 0; k < rd; k++)  bj[k] = aj[k*wh];
 775                                         dwt_encode_97(bj, dn, sn, cas_axl);
 776                                         dwt_deinterleave_z(bj, aj, dn, sn, wh, cas_axl);
 777                                 }
 778                         } else if (dwtid[2] == 1) {
 779                                 for (j = 0; j < (rw*rh); j++) {
 780                                         aj = a + j;
 781                                         for (k = 0; k < rd; k++)  bj[k] = aj[k*wh];
 782                                         dwt_encode_53(bj, dn, sn, cas_axl);
 783                                         dwt_deinterleave_z(bj, aj, dn, sn, wh, cas_axl);
 784                                 }
 785                         }
 786                         opj_free(bj);
 787                 }
 788         }
 789
 790         /*fprintf(stdout,"[INFO] Ops: %d \n",ops);*/
 791 }
 792
 793
 794 /* <summary>                            */
 795 /* Inverse 5-3 wavelet tranform in 3-D. */
 796 /* </summary>                           */
 797 void dwt_decode(opj_tcd_tilecomp_t * tilec, int stops[3], int dwtid[3]) {
 798         int i, j, k;
 799         int x, y, z;
 800         int w, h, wh, d;
 801         int level, levelx, levely, levelz, diff;
 802         int *a = NULL;
 803         int *aj = NULL;
 804         int *bj = NULL;
 805         int *cj = NULL;
 806
 807         a = tilec->data;
 808
 809         w = tilec->x1-tilec->x0;
 810         h = tilec->y1-tilec->y0;
 811         d = tilec->z1-tilec->z0;
 812         wh = w * h;
 813         levelx = tilec->numresolution[0]-1;
 814         levely = tilec->numresolution[1]-1;
 815         levelz = tilec->numresolution[2]-1;
 816         level = int_max(levelx,int_max(levely,levelz));
 817         diff = tilec->numresolution[0] - tilec->numresolution[2];
 818
 819 /* General lifting framework -- DCCS-LIWT */
 820         for (x = level - 1, y = level - 1, z = level - 1; (x >= stops[0]) && (y >= stops[1]); x--, y--, z--) {
 821                 int rw;                 /* width of the resolution level computed                                                           */
 822                 int rh;                 /* heigth of the resolution level computed                                                          */
 823                 int rd;                 /* depth of the resolution level computed                                                          */
 824                 int rw1;                /* width of the resolution level once lower than computed one                                       */
 825                 int rh1;                /* height of the resolution level once lower than computed one                                      */
 826                 int rd1;                /* depth of the resolution level once lower than computed one                                      */
 827                 int cas_col;    /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 828                 int cas_row;    /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 829                 int cas_axl;    /* 0 = non inversion on axial filtering 1 = inversion between low-pass and high-pass filtering   */
 830                 int dn, sn;
 831
 832                 rw = tilec->resolutions[level - x].x1 - tilec->resolutions[level - x].x0;
 833                 rh = tilec->resolutions[level - y].y1 - tilec->resolutions[level - y].y0;
 834                 rd = tilec->resolutions[level - z].z1 - tilec->resolutions[level - z].z0;
 835                 rw1= tilec->resolutions[level - x - 1].x1 - tilec->resolutions[level - x - 1].x0;
 836                 rh1= tilec->resolutions[level - y - 1].y1 - tilec->resolutions[level - y - 1].y0;
 837                 rd1= tilec->resolutions[level - z - 1].z1 - tilec->resolutions[level - z - 1].z0;
 838
 839                 cas_col = tilec->resolutions[level - x].x0 % 2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
 840                 cas_row = tilec->resolutions[level - y].y0 % 2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering   */
 841                 cas_axl = tilec->resolutions[level - z].z0 % 2;
 842
 843                 /*fprintf(stdout," x %d y %d z %d \n",x,y,z);
 844                 fprintf(stdout," levelx %d levely %d levelz %d \n",levelx,levely,levelz);
 845                 fprintf(stdout," dwtid[0] %d [1] %d [2] %d \n",dwtid[0],dwtid[1],dwtid[2]);
 846                 fprintf(stdout," rw %d rh %d rd %d \n rw1 %d rh1 %d rd1 %d \n",rw,rh,rd,rw1,rh1,rd1);
 847                 fprintf(stdout,"IDWT Transform %d %d %d %d\n",level, z, rd,rd1);*/
 848
 849                 if (z >= stops[2] && rd != rd1) {
 850                         /*fprintf(stdout,"Axial Transform %d %d %d %d\n",levelz, z, rd,rd1);*/
 851                         sn = rd1;
 852                         dn = rd - rd1;
 853                         bj = (int*)opj_malloc(rd * sizeof(int));
 854                         if (dwtid[2] == 0) {
 855                                 for (j = 0; j < (rw*rh); j++) {
 856                                         aj = a + j;
 857                                         dwt_interleave_z(aj, bj, dn, sn, wh, cas_axl);
 858                                         dwt_decode_97(bj, dn, sn, cas_axl);
 859                                         for (k = 0; k < rd; k++)  aj[k * wh] = bj[k];
 860                                 }
 861                         } else if (dwtid[2] == 1) {
 862                                 for (j = 0; j < (rw*rh); j++) {
 863                                         aj = a + j;
 864                                         dwt_interleave_z(aj, bj, dn, sn, wh, cas_axl);
 865                                         dwt_decode_53(bj, dn, sn, cas_axl);
 866                                         for (k = 0; k < rd; k++)  aj[k * wh] = bj[k];
 867                                 }
 868                         }
 869                         opj_free(bj);
 870                 }
 871
 872                 for (i = 0; i < rd; i++) {
 873                         /*Fetch corresponding slice for doing DWT-2D*/
 874                         cj = tilec->data + (i * wh);
 875
 876                         /*Vertical*/
 877                         sn = rh1;
 878                         dn = rh - rh1;
 879                         bj = (int*)opj_malloc(rh * sizeof(int));
 880                         if (dwtid[1] == 0) {
 881                                 for (j = 0; j < rw; j++) {
 882                                         aj = cj + j;
 883                                         dwt_interleave_v(aj, bj, dn, sn, w, cas_col);
 884                                         dwt_decode_97(bj, dn, sn, cas_col);
 885                                         for (k = 0; k < rh; k++)  aj[k * w] = bj[k];
 886                                 }
 887                         } else if (dwtid[1] == 1) {
 888                                 for (j = 0; j < rw; j++) {
 889                                         aj = cj + j;
 890                                         dwt_interleave_v(aj, bj, dn, sn, w, cas_col);
 891                                         dwt_decode_53(bj, dn, sn, cas_col);
 892                                         for (k = 0; k < rh; k++)  aj[k * w] = bj[k];
 893                                 }
 894                         }
 895                         opj_free(bj);
 896
 897                         /*Horizontal*/
 898                         sn = rw1;
 899                         dn = rw - rw1;
 900                         bj = (int*)opj_malloc(rw * sizeof(int));
 901                         if (dwtid[0]==0) {
 902                                 for (j = 0; j < rh; j++) {
 903                                         aj = cj + j*w;
 904                                         dwt_interleave_h(aj, bj, dn, sn, cas_row);
 905                                         dwt_decode_97(bj, dn, sn, cas_row);
 906                                         for (k = 0; k < rw; k++)  aj[k] = bj[k];
 907                                 }
 908                         } else if (dwtid[0]==1) {
 909                                 for (j = 0; j < rh; j++) {
 910                                         aj = cj + j*w;
 911                                         dwt_interleave_h(aj, bj, dn, sn, cas_row);
 912                                         dwt_decode_53(bj, dn, sn, cas_row);
 913                                         for (k = 0; k < rw; k++)  aj[k] = bj[k];
 914                                 }
 915                         }
 916                         opj_free(bj);
 917
 918                 }
 919
 920         }
 921
 922 }
 923
 924
 925 /* <summary>                          */
 926 /* Get gain of wavelet transform. */
 927 /* </summary>                         */
 928 int dwt_getgain(int orient, int reversible) {
 929         if (reversible == 1) {
 930                 if (orient == 0)
 931                         return 0;
 932                 else if (orient == 1 || orient == 2 || orient == 4 )
 933                         return 1;
 934                 else if (orient == 3 || orient == 5 || orient == 6 )
 935                         return 2;
 936                 else
 937                         return 3;
 938         }
 939         /*else if (reversible == 0){*/
 940         return 0;
 941 }
 942
 943 /* <summary>                */
 944 /* Get norm of wavelet transform. */
 945 /* </summary>               */
 946 double dwt_getnorm(int orient, int level[3], int dwtid[3]) {
 947         int levelx = level[0];
 948         int levely = level[1];
 949         int levelz = (level[2] < 0) ? 0 : level[2];
 950         double norm;
 951
 952         if (flagnorm[levelx][levely][levelz][orient] == 1) {
 953                 norm = dwt_norm[levelx][levely][levelz][orient];
 954                 /*fprintf(stdout,"[INFO] Level: %d %d %d Orient %d Dwt_norm: %f \n",level[0],level[1],level[2],orient,norm);*/
 955         } else {
 956                 opj_wtfilt_t *wtfiltx =(opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
 957                 opj_wtfilt_t *wtfilty =(opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
 958                 opj_wtfilt_t *wtfiltz =(opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
 959                 /*Fetch equivalent filters for each dimension*/
 960                 dwt_getwtfilters(wtfiltx, dwtid[0]);
 961                 dwt_getwtfilters(wtfilty, dwtid[1]);
 962                 dwt_getwtfilters(wtfiltz, dwtid[2]);
 963                 /*Calculate the corresponding norm */
 964                 norm = dwt_calc_wtnorms(orient, level, dwtid, wtfiltx, wtfilty, wtfiltz);
 965                 /*Save norm in array (no recalculation)*/
 966                 dwt_norm[levelx][levely][levelz][orient] = norm;
 967                 flagnorm[levelx][levely][levelz][orient] = 1;
 968                 /*Free reserved space*/
 969                 opj_free(wtfiltx->LPS); opj_free(wtfilty->LPS); opj_free(wtfiltz->LPS);
 970                 opj_free(wtfiltx->HPS); opj_free(wtfilty->HPS); opj_free(wtfiltz->HPS);
 971                 opj_free(wtfiltx);              opj_free(wtfilty);              opj_free(wtfiltz);
 972                 /*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);*/
 973         }
 974         return norm;
 975 }
 976 /* <summary>                                                            */
 977 /* Calculate explicit stepsizes for DWT.        */
 978 /* </summary>                                                           */
 979 void dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, int prec) {
 980         int totnumbands, bandno, diff;
 981
 982         assert(tccp->numresolution[0] >= tccp->numresolution[2]);
 983         diff = tccp->numresolution[0] - tccp->numresolution[2];         /*if RESx=RESy != RESz */
 984         totnumbands = (7 * tccp->numresolution[0] - 6) - 4 * diff; /* 3-D */
 985
 986         for (bandno = 0; bandno < totnumbands; bandno++) {
 987                 double stepsize;
 988                 int resno, level[3], orient, gain;
 989
 990                 /* Bandno:      0 - LLL         1 - LHL
 991                                         2 - HLL         3 - HHL
 992                                         4 - LLH         5 - LHH
 993                                         6 - HLH         7 - HHH */
 994
 995                 resno = (bandno == 0) ? 0 : ( (bandno <= 3 * diff) ? ((bandno - 1) / 3 + 1) : ((bandno + 4*diff - 1) / 7 + 1));
 996                 orient = (bandno == 0) ? 0 : ( (bandno <= 3 * diff) ? ((bandno - 1) % 3 + 1) : ((bandno + 4*diff - 1) % 7 + 1));
 997                 level[0] = tccp->numresolution[0] - 1 - resno;
 998                 level[1] = tccp->numresolution[1] - 1 - resno;
 999                 level[2] = tccp->numresolution[2] - 1 - resno;
1000
1001                 /* Gain:        0 - LLL         1 - LHL
1002                                         1 - HLL         2 - HHL
1003                                         1 - LLH         2 - LHH
1004                                         2 - HLH         3 - HHH         */
1005                 gain = (tccp->reversible == 0) ? 0 : ( (orient == 0) ? 0 :
1006                                 ( ((orient == 1) || (orient == 2) || (orient == 4)) ? 1 :
1007                                                 (((orient == 3) || (orient == 5) || (orient == 6)) ? 2 : 3)) );
1008
1009                 if (tccp->qntsty == J3D_CCP_QNTSTY_NOQNT) {
1010                         stepsize = 1.0;
1011                 } else {
1012                         double norm = dwt_getnorm(orient,level,tccp->dwtid); /*Fetch norms if irreversible transform (by the moment only I9.7)*/
1013                         stepsize = (1 << (gain + 1)) / norm;
1014                 }
1015                 /*fprintf(stdout,"[INFO] Bandno: %d Orient: %d Level: %d %d %d Stepsize: %f\n",bandno,orient,level[0],level[1],level[2],stepsize);*/
1016                 dwt_encode_stepsize((int) floor(stepsize * 8192.0), prec + gain, &tccp->stepsizes[bandno]);
1017         }
1018 }
1019
1020
1021