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