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, Herve 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
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions
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.
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.
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.
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.
46 * I have also removed the "Add Patrick" part because it is not longer
50 * -Ive (aka Reiner Wahler)
51 * mail: ive@lilysoft.com
54 #include "opj_includes.h"
56 /** @defgroup DWT DWT - Implementation of a discrete wavelet transform */
59 /** @name Local static functions */
63 Forward lazy transform (horizontal)
65 static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas);
67 Forward lazy transform (vertical)
69 static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas);
71 Forward lazy transform (axial)
73 static void dwt_deinterleave_z(int *a, int *b, int dn, int sn, int xy, int cas);
75 Inverse lazy transform (horizontal)
77 static void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas);
79 Inverse lazy transform (vertical)
81 static void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas);
83 Inverse lazy transform (axial)
85 static void dwt_interleave_z(int *a, int *b, int dn, int sn, int xy, int cas);
87 Forward 5-3 wavelet tranform in 1-D
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);
92 Inverse 5-3 wavelet tranform in 1-D
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);
97 Computing of wavelet transform L2 norms for arbitrary transforms
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);
101 Encoding of quantification stepsize
103 static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno_stepsize);
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)))
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)))
117 /* This table contains the norms of the 5-3 wavelets for different bands. */
119 static double dwt_norm[10][10][10][8];
120 static int flagnorm[10][10][10][8];
122 /*static const double dwt_norms[5][8][10] = {
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}
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},
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},
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 },
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 },
168 /* This table contains the norms of the 9-7 wavelets for different bands. */
170 /*static const double dwt_norms_real[5][8][10] = {
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}
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},
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},
193 {0.3752, 0.9512} //HHH
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
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
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*/
228 ==========================================================
230 ==========================================================
234 /* Forward lazy transform (horizontal). */
236 static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas) {
238 for (i=0; i<sn; i++) b[i]=a[2*i+cas];
239 for (i=0; i<dn; i++) b[sn+i]=a[(2*i+1-cas)];
243 /* Forward lazy transform (vertical). */
245 static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
247 for (i=0; i<sn; i++) b[i*x]=a[2*i+cas];
248 for (i=0; i<dn; i++) b[(sn+i)*x]=a[(2*i+1-cas)];
252 /* Forward lazy transform (axial). */
254 static void dwt_deinterleave_z(int *a, int *b, int dn, int sn, int xy, int cas) {
256 for (i=0; i<sn; i++) b[i*xy]=a[2*i+cas];
257 for (i=0; i<dn; i++) b[(sn+i)*xy]=a[(2*i+1-cas)];
261 /* Inverse lazy transform (horizontal). */
263 static void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas) {
269 for (i = 0; i < sn; i++) {
276 for (i = 0; i < dn; i++) {
284 /* Inverse lazy transform (vertical). */
286 static void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
292 for (i = 0; i < sn; i++) {
299 for (i = 0; i < dn; i++) {
307 /* Inverse lazy transform (axial). */
309 static void dwt_interleave_z(int *a, int *b, int dn, int sn, int xy, int cas) {
315 for (i = 0; i < sn; i++) {
322 for (i = 0; i < dn; i++) {
331 /* Forward 5-3 or 9-7 wavelet tranform in 1-D. */
333 static void dwt_encode_53(int *a, int dn, int sn, int 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;
344 for (i = 0; i < sn; i++){
345 S(i) += (D_(i - 1) + D_(i) + 2) >> 2;
350 /*if (!sn && dn == 1)
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;
360 for (i = 0; i < dn; i++){
361 S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
364 for (i = 0; i < sn; i++){
365 D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
371 static void dwt_encode_97(int *a, int dn, int sn, int 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 */
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 */
407 /* Inverse 5-3 or 9-7 wavelet tranform in 1-D. */
409 static void dwt_decode_53(int *a, int dn, int sn, int 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;
417 if (!sn && dn == 1) /* NEW : CASE ONE ELEMENT */
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;
425 static void dwt_decode_97(int *a, int dn, int sn, int 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 */
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 */
463 /* Get norm of arbitrary wavelet transform. */
465 static int upandconv(double *nXPS, double *LPS, int lenXPS, int lenLPS) {
466 /* Perform the convolution of the vectors. */
468 double *tmp = (double *)opj_malloc(2*lenXPS * sizeof(double));
470 memset(tmp, 0, 2*lenXPS*sizeof(double));
471 for (i = 0; i < lenXPS; i++) {
472 *(tmp + 2*i) = *(nXPS + i);
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));
483 return 2*lenXPS+lenLPS-1;
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;
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];
497 lenLPS = wtfiltX->lenLPS;
498 lenHPS = wtfiltX->lenHPS;
499 for (i = 0; i < levelx; i++) {
502 lenLPS += wtfiltX->lenLPS - 1;
503 lenHPS += wtfiltX->lenLPS - 1;
505 nLPSx = (double *)opj_malloc(lenLPS * sizeof(double));
506 nHPSx = (double *)opj_malloc(lenHPS * sizeof(double));
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);
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];
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++) {
532 lenLPS += wtfiltY->lenLPS - 1;
533 lenHPS += wtfiltY->lenLPS - 1;
535 nLPSy = (double *)opj_malloc(lenLPS * sizeof(double));
536 nHPSy = (double *)opj_malloc(lenHPS * sizeof(double));
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);
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];
560 lenLPS = wtfiltZ->lenLPS;
561 lenHPS = wtfiltZ->lenHPS;
562 for (i = 0; i < levelz; i++) {
565 lenLPS += wtfiltZ->lenLPS - 1;
566 lenHPS += wtfiltZ->lenLPS - 1;
568 nLPSz = (double *)opj_malloc(lenLPS * sizeof(double));
569 nHPSz = (double *)opj_malloc(lenHPS * sizeof(double));
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);
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];
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;
636 fprintf(stdout,"[ERROR] Sorry, this wavelet hasn't been implemented so far ... Try another one :-)\n");
641 /* Encoding of quantization stepsize for each subband. */
643 static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno_stepsize) {
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
654 ==========================================================
656 ==========================================================
659 /* Forward 5-3 wavelet tranform in 3-D. */
661 void dwt_encode(opj_tcd_tilecomp_t * tilec, int dwtid[3]) {
665 int level,levelx,levely,levelz,diff;
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;
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];
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 */
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;
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;
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);*/
714 for (i = 0; i < rd; i++) {
721 bj = (int*)opj_malloc(rw * sizeof(int));
723 for (j = 0; j < rh; j++) {
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);
729 } else if (dwtid[0] == 1) {
730 for (j = 0; j < rh; j++) {
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);
742 bj = (int*)opj_malloc(rh * sizeof(int));
743 if (dwtid[1] == 0) { /*DWT 9-7*/
744 for (j = 0; j < rw; 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);
750 } else if (dwtid[1] == 1) { /*DWT 5-3*/
751 for (j = 0; j < rw; 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);
762 //Axial fprintf(stdout,"Axial DWT Transform %d %d %d\n",z,rd,rd1);
765 bj = (int*)opj_malloc(rd * sizeof(int));
767 for (j = 0; j < (rw*rh); 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);
773 } else if (dwtid[2] == 1) {
774 for (j = 0; j < (rw*rh); 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);
785 //fprintf(stdout,"[INFO] Ops: %d \n",ops);
790 /* Inverse 5-3 wavelet tranform in 3-D. */
792 void dwt_decode(opj_tcd_tilecomp_t * tilec, int stops[3], int dwtid[3]) {
796 int level, levelx, levely, levelz, diff;
804 w = tilec->x1-tilec->x0;
805 h = tilec->y1-tilec->y0;
806 d = tilec->z1-tilec->z0;
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];
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 */
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;
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;
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);*/
844 if (z >= stops[2] && rd != rd1) {
845 //fprintf(stdout,"Axial Transform %d %d %d %d\n",levelz, z, rd,rd1);
848 bj = (int*)opj_malloc(rd * sizeof(int));
850 for (j = 0; j < (rw*rh); 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];
856 } else if (dwtid[2] == 1) {
857 for (j = 0; j < (rw*rh); 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];
867 for (i = 0; i < rd; i++) {
868 //Fetch corresponding slice for doing DWT-2D
869 cj = tilec->data + (i * wh);
874 bj = (int*)opj_malloc(rh * sizeof(int));
876 for (j = 0; j < rw; 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];
882 } else if (dwtid[1] == 1) {
883 for (j = 0; j < rw; 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];
895 bj = (int*)opj_malloc(rw * sizeof(int));
897 for (j = 0; j < rh; j++) {
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];
903 } else if (dwtid[0]==1) {
904 for (j = 0; j < rh; j++) {
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];
921 /* Get gain of wavelet transform. */
923 int dwt_getgain(int orient, int reversible) {
924 if (reversible == 1) {
927 else if (orient == 1 || orient == 2 || orient == 4 )
929 else if (orient == 3 || orient == 5 || orient == 6 )
934 //else if (reversible == 0){
939 /* Get norm of wavelet transform. */
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];
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);
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);
972 /* Calculate explicit stepsizes for DWT. */
974 void dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, int prec) {
975 int totnumbands, bandno, diff;
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 */
981 for (bandno = 0; bandno < totnumbands; bandno++) {
983 int resno, level[3], orient, gain;
985 /* Bandno: 0 - LLL 1 - LHL
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;
996 /* Gain: 0 - LLL 1 - LHL
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)) );
1004 if (tccp->qntsty == J3D_CCP_QNTSTY_NOQNT) {
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;
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]);