summaryrefslogtreecommitdiff
path: root/src/lib/openjp3d/dwt.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/lib/openjp3d/dwt.c')
-rw-r--r--src/lib/openjp3d/dwt.c1683
1 files changed, 910 insertions, 773 deletions
diff --git a/src/lib/openjp3d/dwt.c b/src/lib/openjp3d/dwt.c
index a1e4e104..90a807ed 100644
--- a/src/lib/openjp3d/dwt.c
+++ b/src/lib/openjp3d/dwt.c
@@ -1,6 +1,6 @@
/*
- * The copyright in this software is being made available under the 2-clauses
- * BSD License, included below. This software may be subject to other third
+ * The copyright in this software is being made available under the 2-clauses
+ * BSD License, included below. This software may be subject to other third
* party and contributor rights, including patent rights, and no such rights
* are granted under this license.
*
@@ -49,7 +49,7 @@
* the major bottleneck in the former version.
*
* I have also removed the "Add Patrick" part because it is not longer
- * needed.
+ * needed.
*
* 6/6/2005
* -Ive (aka Reiner Wahler)
@@ -101,11 +101,13 @@ static void dwt_decode_97(int *a, int dn, int sn, int cas);
/**
Computing of wavelet transform L2 norms for arbitrary transforms
*/
-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);
+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);
/**
Encoding of quantification stepsize
*/
-static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno_stepsize);
+static void dwt_encode_stepsize(int stepsize, int numbps,
+ opj_stepsize_t *bandno_stepsize);
/*@}*/
/*@}*/
@@ -125,209 +127,230 @@ static double dwt_norm[10][10][10][8];
static int flagnorm[10][10][10][8];
/*static const double dwt_norms[5][8][10] = {
- {//ResZ=1
- {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
- {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
- {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
- {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
- },{//ResZ=2
- {1.000, 1.8371, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
- {1.2717, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
- {1.2717, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
- {.8803, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
- {1.2717},
- {.8803},
- {.8803},
- {.6093},
- },{ //ResZ=3
- {1.000, 1.8371, 4.5604, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
- {1.2717, 2.6403, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
- {1.2717, 2.6403, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
- {.8803, 1.5286, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
- {1.2717, 2.6403},
- {.8803, 1.5286},
- {.8803, 1.5286},
- {.6093, 0.8850},
- },{ //ResZ=4
- {1.000, 1.8371, 4.5604, 12.4614, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
- {1.2717, 2.6403, 6.7691 , 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
- {1.2717, 2.6403, 6.7691 , 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
- {.8803, 1.5286, 3.6770 , 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
- {1.2717, 2.6403, 6.7691 },
- {.8803, 1.5286, 3.6770 },
- {.8803, 1.5286, 3.6770 },
- {.6093, 0.8850, 1.9974 },
- },{ //ResZ=5
- {1.000, 1.8371, 4.5604, 12.4614, 34.9025, 21.34, 42.67, 85.33, 170.7, 341.3},
- {1.2717, 2.6403, 6.7691 , 18.6304 , 11.33, 22.64, 45.25, 90.48, 180.9},
- {1.2717, 2.6403, 6.7691 , 18.6304, 11.33, 22.64, 45.25, 90.48, 180.9},
- {.8803, 1.5286, 3.6770 , 9.9446, 6.019, 12.01, 24.00, 47.97, 95.93},
- {1.2717, 2.6403, 6.7691, 18.6304},
- {.8803, 1.5286, 3.6770, 9.9446 },
- {.8803, 1.5286, 3.6770, 9.9446 },
- {.6093, 0.8850, 1.9974, 5.3083 },
- }
+ {//ResZ=1
+ {1.000, 1.500, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
+ {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
+ {1.038, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
+ {.7186, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93}
+ },{//ResZ=2
+ {1.000, 1.8371, 2.750, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
+ {1.2717, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
+ {1.2717, 1.592, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
+ {.8803, .9218, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
+ {1.2717},
+ {.8803},
+ {.8803},
+ {.6093},
+ },{ //ResZ=3
+ {1.000, 1.8371, 4.5604, 5.375, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
+ {1.2717, 2.6403, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
+ {1.2717, 2.6403, 2.919, 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
+ {.8803, 1.5286, 1.586, 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
+ {1.2717, 2.6403},
+ {.8803, 1.5286},
+ {.8803, 1.5286},
+ {.6093, 0.8850},
+ },{ //ResZ=4
+ {1.000, 1.8371, 4.5604, 12.4614, 10.68, 21.34, 42.67, 85.33, 170.7, 341.3},
+ {1.2717, 2.6403, 6.7691 , 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
+ {1.2717, 2.6403, 6.7691 , 5.703, 11.33, 22.64, 45.25, 90.48, 180.9},
+ {.8803, 1.5286, 3.6770 , 3.043, 6.019, 12.01, 24.00, 47.97, 95.93},
+ {1.2717, 2.6403, 6.7691 },
+ {.8803, 1.5286, 3.6770 },
+ {.8803, 1.5286, 3.6770 },
+ {.6093, 0.8850, 1.9974 },
+ },{ //ResZ=5
+ {1.000, 1.8371, 4.5604, 12.4614, 34.9025, 21.34, 42.67, 85.33, 170.7, 341.3},
+ {1.2717, 2.6403, 6.7691 , 18.6304 , 11.33, 22.64, 45.25, 90.48, 180.9},
+ {1.2717, 2.6403, 6.7691 , 18.6304, 11.33, 22.64, 45.25, 90.48, 180.9},
+ {.8803, 1.5286, 3.6770 , 9.9446, 6.019, 12.01, 24.00, 47.97, 95.93},
+ {1.2717, 2.6403, 6.7691, 18.6304},
+ {.8803, 1.5286, 3.6770, 9.9446 },
+ {.8803, 1.5286, 3.6770, 9.9446 },
+ {.6093, 0.8850, 1.9974, 5.3083 },
+ }
};*/
/* <summary> */
/* This table contains the norms of the 9-7 wavelets for different bands. */
/* </summary> */
/*static const double dwt_norms_real[5][8][10] = {
- {//ResZ==1
- {1.000, 1.9659, 4.1224, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
- {1.0113, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
- {1.0113, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
- {0.5202, 0.9672, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722}
- }, { //ResZ==2
- {1.000, 2.7564, 4.1224, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
- {1.4179, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
- {1.4179, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
- {0.7294, 0.9672, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
- {1.4179},
- {0.7294},
- {0.7294},
- {0.3752} //HHH
- },{ //ResZ==3
- {1.000, 2.7564, 8.3700, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
- {1.4179, 4.0543, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
- {1.4179, 4.0543, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
- {0.7294, 1.9638, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
- {1.4179, 4.0543},
- {0.7294, 1.9638},
- {0.7294, 1.9638},
- {0.3752, 0.9512} //HHH
- },{ //ResZ==4
- {1.000, 2.7564, 8.3700, 24.4183, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
- {1.4179, 4.0543, 12.1366, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
- {1.4179, 4.0543, 12.1366, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
- {0.7294, 1.9638, 6.0323, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
- {1.4179, 4.0543, 12.1366},
- {0.7294, 1.9638, 6.0323},
- {0.7294, 1.9638, 6.0323},
- {0.3752, 0.9512, 2.9982} //HHH
- },{ //ResZ==5
- {1.000, 2.7564, 8.3700, 24.4183, 69.6947, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
- {1.4179, 4.0543, 12.1366, 35.1203, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
- {1.4179, 4.0543, 12.1366, 35.1203, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
- {0.7294, 1.9638, 6.0323, 17.6977, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
- {1.4179, 4.0543, 12.1366, 35.1203},
- {0.7294, 1.9638, 6.0323, 17.6977},
- {0.7294, 1.9638, 6.0323, 17.6977},
- {0.3752, 0.9512, 2.9982, 8.9182} //HHH
- }
+ {//ResZ==1
+ {1.000, 1.9659, 4.1224, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
+ {1.0113, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
+ {1.0113, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
+ {0.5202, 0.9672, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722}
+ }, { //ResZ==2
+ {1.000, 2.7564, 4.1224, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
+ {1.4179, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
+ {1.4179, 1.9968, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
+ {0.7294, 0.9672, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
+ {1.4179},
+ {0.7294},
+ {0.7294},
+ {0.3752} //HHH
+ },{ //ResZ==3
+ {1.000, 2.7564, 8.3700, 8.4167, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
+ {1.4179, 4.0543, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
+ {1.4179, 4.0543, 4.1834, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
+ {0.7294, 1.9638, 2.0793, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
+ {1.4179, 4.0543},
+ {0.7294, 1.9638},
+ {0.7294, 1.9638},
+ {0.3752, 0.9512} //HHH
+ },{ //ResZ==4
+ {1.000, 2.7564, 8.3700, 24.4183, 16.9356, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
+ {1.4179, 4.0543, 12.1366, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
+ {1.4179, 4.0543, 12.1366, 8.5341, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
+ {0.7294, 1.9638, 6.0323, 4.3005, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
+ {1.4179, 4.0543, 12.1366},
+ {0.7294, 1.9638, 6.0323},
+ {0.7294, 1.9638, 6.0323},
+ {0.3752, 0.9512, 2.9982} //HHH
+ },{ //ResZ==5
+ {1.000, 2.7564, 8.3700, 24.4183, 69.6947, 33.9249, 67.8772, 135.7680, 271.5430, 543.0894},
+ {1.4179, 4.0543, 12.1366, 35.1203, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
+ {1.4179, 4.0543, 12.1366, 35.1203, 17.1667, 34.3852, 68.7967, 137.6065, 275.2196},
+ {0.7294, 1.9638, 6.0323, 17.6977, 8.6867, 17.4188, 34.8608, 69.7332, 139.4722},
+ {1.4179, 4.0543, 12.1366, 35.1203},
+ {0.7294, 1.9638, 6.0323, 17.6977},
+ {0.7294, 1.9638, 6.0323, 17.6977},
+ {0.3752, 0.9512, 2.9982, 8.9182} //HHH
+ }
};*/
static opj_atk_t atk_info_wt[] = {
- {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*/
- {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*/
- {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*/
- {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*/
- {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*/
- {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*/
- {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},
- {-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*/
- {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*/
- {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*/
+ {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*/
+ {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*/
+ {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*/
+ {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*/
+ {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*/
+ {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*/
+ {
+ 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},
+ {-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*/
+ {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*/
+ {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*/
};
-/*
+/*
==========================================================
local functions
==========================================================
*/
-/* <summary> */
+/* <summary> */
/* Forward lazy transform (horizontal). */
-/* </summary> */
-static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas) {
- int i;
- for (i=0; i<sn; i++) b[i]=a[2*i+cas];
- for (i=0; i<dn; i++) b[sn+i]=a[(2*i+1-cas)];
+/* </summary> */
+static void dwt_deinterleave_h(int *a, int *b, int dn, int sn, int cas)
+{
+ int i;
+ for (i = 0; i < sn; i++) {
+ b[i] = a[2 * i + cas];
+ }
+ for (i = 0; i < dn; i++) {
+ b[sn + i] = a[(2 * i + 1 - cas)];
+ }
}
-/* <summary> */
+/* <summary> */
/* Forward lazy transform (vertical). */
-/* </summary> */
-static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
+/* </summary> */
+static void dwt_deinterleave_v(int *a, int *b, int dn, int sn, int x, int cas)
+{
int i;
- for (i=0; i<sn; i++) b[i*x]=a[2*i+cas];
- for (i=0; i<dn; i++) b[(sn+i)*x]=a[(2*i+1-cas)];
+ for (i = 0; i < sn; i++) {
+ b[i * x] = a[2 * i + cas];
+ }
+ for (i = 0; i < dn; i++) {
+ b[(sn + i)*x] = a[(2 * i + 1 - cas)];
+ }
}
-/* <summary> */
+/* <summary> */
/* Forward lazy transform (axial). */
-/* </summary> */
-static void dwt_deinterleave_z(int *a, int *b, int dn, int sn, int xy, int cas) {
+/* </summary> */
+static void dwt_deinterleave_z(int *a, int *b, int dn, int sn, int xy, int cas)
+{
int i;
- for (i=0; i<sn; i++) b[i*xy]=a[2*i+cas];
- for (i=0; i<dn; i++) b[(sn+i)*xy]=a[(2*i+1-cas)];
+ for (i = 0; i < sn; i++) {
+ b[i * xy] = a[2 * i + cas];
+ }
+ for (i = 0; i < dn; i++) {
+ b[(sn + i)*xy] = a[(2 * i + 1 - cas)];
+ }
}
/* <summary> */
/* Inverse lazy transform (horizontal). */
/* </summary> */
-static void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas) {
+static void dwt_interleave_h(int *a, int *b, int dn, int sn, int cas)
+{
int i;
int *ai = NULL;
int *bi = NULL;
ai = a;
bi = b + cas;
for (i = 0; i < sn; i++) {
- *bi = *ai;
- bi += 2;
- ai++;
+ *bi = *ai;
+ bi += 2;
+ ai++;
}
ai = a + sn;
bi = b + 1 - cas;
for (i = 0; i < dn; i++) {
- *bi = *ai;
- bi += 2;
- ai++;
+ *bi = *ai;
+ bi += 2;
+ ai++;
}
}
-/* <summary> */
+/* <summary> */
/* Inverse lazy transform (vertical). */
-/* </summary> */
-static void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas) {
+/* </summary> */
+static void dwt_interleave_v(int *a, int *b, int dn, int sn, int x, int cas)
+{
int i;
int *ai = NULL;
int *bi = NULL;
ai = a;
bi = b + cas;
for (i = 0; i < sn; i++) {
- *bi = *ai;
- bi += 2;
- ai += x;
+ *bi = *ai;
+ bi += 2;
+ ai += x;
}
ai = a + (sn * x);
bi = b + 1 - cas;
for (i = 0; i < dn; i++) {
- *bi = *ai;
- bi += 2;
- ai += x;
+ *bi = *ai;
+ bi += 2;
+ ai += x;
}
}
/* <summary> */
/* Inverse lazy transform (axial). */
/* </summary> */
-static void dwt_interleave_z(int *a, int *b, int dn, int sn, int xy, int cas) {
+static void dwt_interleave_z(int *a, int *b, int dn, int sn, int xy, int cas)
+{
int i;
int *ai = NULL;
int *bi = NULL;
ai = a;
bi = b + cas;
for (i = 0; i < sn; i++) {
- *bi = *ai;
- bi += 2;
- ai += xy;
+ *bi = *ai;
+ bi += 2;
+ ai += xy;
}
ai = a + (sn * xy);
bi = b + 1 - cas;
for (i = 0; i < dn; i++) {
- *bi = *ai;
- bi += 2;
- ai += xy;
+ *bi = *ai;
+ bi += 2;
+ ai += xy;
}
}
@@ -335,327 +358,389 @@ static void dwt_interleave_z(int *a, int *b, int dn, int sn, int xy, int cas) {
/* <summary> */
/* Forward 5-3 or 9-7 wavelet transform in 1-D. */
/* </summary> */
-static void dwt_encode_53(int *a, int dn, int sn, int cas) {
- int i;
-
- if (!cas) {
- if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
- /*for (i = 0; i < dn; i++) D(i) -= (S_(i) + S_(i + 1)) >> 1;*/
- /*for (i = 0; i < sn; i++) S(i) += (D_(i - 1) + D_(i) + 2) >> 2;*/
- for (i = 0; i < dn; i++){
- D(i) -= (S_(i) + S_(i + 1)) >> 1;
- /*ops += 2;*/
- }
- for (i = 0; i < sn; i++){
- S(i) += (D_(i - 1) + D_(i) + 2) >> 2;
- /*ops += 3;*/
- }
- }
- } else {
- /*if (!sn && dn == 1)
- S(0) *= 2;
- else {
- for (i = 0; i < dn; i++) S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
- for (i = 0; i < sn; i++) D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
- }*/
- if (!sn && dn == 1){
- S(0) *= 2;
- /*ops++;*/
- } else {
- for (i = 0; i < dn; i++){
- S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
- /* ops += 2;*/
- }
- for (i = 0; i < sn; i++){
- D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
- /* ops += 3;*/
- }
- }
- }
+static void dwt_encode_53(int *a, int dn, int sn, int cas)
+{
+ int i;
+
+ if (!cas) {
+ if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
+ /*for (i = 0; i < dn; i++) D(i) -= (S_(i) + S_(i + 1)) >> 1;*/
+ /*for (i = 0; i < sn; i++) S(i) += (D_(i - 1) + D_(i) + 2) >> 2;*/
+ for (i = 0; i < dn; i++) {
+ D(i) -= (S_(i) + S_(i + 1)) >> 1;
+ /*ops += 2;*/
+ }
+ for (i = 0; i < sn; i++) {
+ S(i) += (D_(i - 1) + D_(i) + 2) >> 2;
+ /*ops += 3;*/
+ }
+ }
+ } else {
+ /*if (!sn && dn == 1)
+ S(0) *= 2;
+ else {
+ for (i = 0; i < dn; i++) S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
+ for (i = 0; i < sn; i++) D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
+ }*/
+ if (!sn && dn == 1) {
+ S(0) *= 2;
+ /*ops++;*/
+ } else {
+ for (i = 0; i < dn; i++) {
+ S(i) -= (DD_(i) + DD_(i - 1)) >> 1;
+ /* ops += 2;*/
+ }
+ for (i = 0; i < sn; i++) {
+ D(i) += (SS_(i) + SS_(i + 1) + 2) >> 2;
+ /* ops += 3;*/
+ }
+ }
+ }
}
-static void dwt_encode_97(int *a, int dn, int sn, int cas) {
- int i;
-
- if (!cas) {
- if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
- for (i = 0; i < dn; i++)
- D(i) -= fix_mul(S_(i) + S_(i + 1), 12993);
- for (i = 0; i < sn; i++)
- S(i) -= fix_mul(D_(i - 1) + D_(i), 434);
- for (i = 0; i < dn; i++)
- D(i) += fix_mul(S_(i) + S_(i + 1), 7233);
- for (i = 0; i < sn; i++)
- S(i) += fix_mul(D_(i - 1) + D_(i), 3633);
- for (i = 0; i < dn; i++)
- D(i) = fix_mul(D(i), 5038); /*5038 */
- for (i = 0; i < sn; i++)
- S(i) = fix_mul(S(i), 6659); /*6660 */
- }
- } else {
- if ((sn > 0) || (dn > 1)) { /* NEW : CASE ONE ELEMENT */
- for (i = 0; i < dn; i++)
- S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993);
- for (i = 0; i < sn; i++)
- D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434);
- for (i = 0; i < dn; i++)
- S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233);
- for (i = 0; i < sn; i++)
- D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633);
- for (i = 0; i < dn; i++)
- S(i) = fix_mul(S(i), 5038); /*5038 */
- for (i = 0; i < sn; i++)
- D(i) = fix_mul(D(i), 6659); /*6660 */
- }
- }
+static void dwt_encode_97(int *a, int dn, int sn, int cas)
+{
+ int i;
+
+ if (!cas) {
+ if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
+ for (i = 0; i < dn; i++) {
+ D(i) -= fix_mul(S_(i) + S_(i + 1), 12993);
+ }
+ for (i = 0; i < sn; i++) {
+ S(i) -= fix_mul(D_(i - 1) + D_(i), 434);
+ }
+ for (i = 0; i < dn; i++) {
+ D(i) += fix_mul(S_(i) + S_(i + 1), 7233);
+ }
+ for (i = 0; i < sn; i++) {
+ S(i) += fix_mul(D_(i - 1) + D_(i), 3633);
+ }
+ for (i = 0; i < dn; i++) {
+ D(i) = fix_mul(D(i), 5038); /*5038 */
+ }
+ for (i = 0; i < sn; i++) {
+ S(i) = fix_mul(S(i), 6659); /*6660 */
+ }
+ }
+ } else {
+ if ((sn > 0) || (dn > 1)) { /* NEW : CASE ONE ELEMENT */
+ for (i = 0; i < dn; i++) {
+ S(i) -= fix_mul(DD_(i) + DD_(i - 1), 12993);
+ }
+ for (i = 0; i < sn; i++) {
+ D(i) -= fix_mul(SS_(i) + SS_(i + 1), 434);
+ }
+ for (i = 0; i < dn; i++) {
+ S(i) += fix_mul(DD_(i) + DD_(i - 1), 7233);
+ }
+ for (i = 0; i < sn; i++) {
+ D(i) += fix_mul(SS_(i) + SS_(i + 1), 3633);
+ }
+ for (i = 0; i < dn; i++) {
+ S(i) = fix_mul(S(i), 5038); /*5038 */
+ }
+ for (i = 0; i < sn; i++) {
+ D(i) = fix_mul(D(i), 6659); /*6660 */
+ }
+ }
+ }
}
/* <summary> */
/* Inverse 5-3 or 9-7 wavelet transform in 1-D. */
-/* </summary> */
-static void dwt_decode_53(int *a, int dn, int sn, int cas) {
- int i;
- if (!cas) {
- if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
- for (i = 0; i < sn; i++) S(i) -= (D_(i - 1) + D_(i) + 2) >> 2;
- for (i = 0; i < dn; i++) D(i) += (S_(i) + S_(i + 1)) >> 1;
- }
- } else {
- if (!sn && dn == 1) /* NEW : CASE ONE ELEMENT */
- S(0) /= 2;
- else {
- for (i = 0; i < sn; i++) D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2;
- for (i = 0; i < dn; i++) S(i) += (DD_(i) + DD_(i - 1)) >> 1;
- }
- }
+/* </summary> */
+static void dwt_decode_53(int *a, int dn, int sn, int cas)
+{
+ int i;
+ if (!cas) {
+ if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
+ for (i = 0; i < sn; i++) {
+ S(i) -= (D_(i - 1) + D_(i) + 2) >> 2;
+ }
+ for (i = 0; i < dn; i++) {
+ D(i) += (S_(i) + S_(i + 1)) >> 1;
+ }
+ }
+ } else {
+ if (!sn && dn == 1) { /* NEW : CASE ONE ELEMENT */
+ S(0) /= 2;
+ } else {
+ for (i = 0; i < sn; i++) {
+ D(i) -= (SS_(i) + SS_(i + 1) + 2) >> 2;
+ }
+ for (i = 0; i < dn; i++) {
+ S(i) += (DD_(i) + DD_(i - 1)) >> 1;
+ }
+ }
+ }
}
-static void dwt_decode_97(int *a, int dn, int sn, int cas) {
- int i;
-
- if (!cas) {
- if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
- for (i = 0; i < sn; i++)
- S(i) = fix_mul(S(i), 10078); /* 10076 */
- for (i = 0; i < dn; i++)
- D(i) = fix_mul(D(i), 13318); /* 13320 */
- for (i = 0; i < sn; i++)
- S(i) -= fix_mul(D_(i - 1) + D_(i), 3633);
- for (i = 0; i < dn; i++)
- D(i) -= fix_mul(S_(i) + S_(i + 1), 7233);
- for (i = 0; i < sn; i++)
- S(i) += fix_mul(D_(i - 1) + D_(i), 434);
- for (i = 0; i < dn; i++)
- D(i) += fix_mul(S_(i) + S_(i + 1), 12994); /* 12993 */
- }
- } else {
- if ((sn > 0) || (dn > 1)) { /* NEW : CASE ONE ELEMENT */
- for (i = 0; i < sn; i++)
- D(i) = fix_mul(D(i), 10078); /* 10076 */
- for (i = 0; i < dn; i++)
- S(i) = fix_mul(S(i), 13318); /* 13320 */
- for (i = 0; i < sn; i++)
- D(i) -= fix_mul(SS_(i) + SS_(i + 1), 3633);
- for (i = 0; i < dn; i++)
- S(i) -= fix_mul(DD_(i) + DD_(i - 1), 7233);
- for (i = 0; i < sn; i++)
- D(i) += fix_mul(SS_(i) + SS_(i + 1), 434);
- for (i = 0; i < dn; i++)
- S(i) += fix_mul(DD_(i) + DD_(i - 1), 12994); /* 12993 */
- }
- }
+static void dwt_decode_97(int *a, int dn, int sn, int cas)
+{
+ int i;
+
+ if (!cas) {
+ if ((dn > 0) || (sn > 1)) { /* NEW : CASE ONE ELEMENT */
+ for (i = 0; i < sn; i++) {
+ S(i) = fix_mul(S(i), 10078); /* 10076 */
+ }
+ for (i = 0; i < dn; i++) {
+ D(i) = fix_mul(D(i), 13318); /* 13320 */
+ }
+ for (i = 0; i < sn; i++) {
+ S(i) -= fix_mul(D_(i - 1) + D_(i), 3633);
+ }
+ for (i = 0; i < dn; i++) {
+ D(i) -= fix_mul(S_(i) + S_(i + 1), 7233);
+ }
+ for (i = 0; i < sn; i++) {
+ S(i) += fix_mul(D_(i - 1) + D_(i), 434);
+ }
+ for (i = 0; i < dn; i++) {
+ D(i) += fix_mul(S_(i) + S_(i + 1), 12994); /* 12993 */
+ }
+ }
+ } else {
+ if ((sn > 0) || (dn > 1)) { /* NEW : CASE ONE ELEMENT */
+ for (i = 0; i < sn; i++) {
+ D(i) = fix_mul(D(i), 10078); /* 10076 */
+ }
+ for (i = 0; i < dn; i++) {
+ S(i) = fix_mul(S(i), 13318); /* 13320 */
+ }
+ for (i = 0; i < sn; i++) {
+ D(i) -= fix_mul(SS_(i) + SS_(i + 1), 3633);
+ }
+ for (i = 0; i < dn; i++) {
+ S(i) -= fix_mul(DD_(i) + DD_(i - 1), 7233);
+ }
+ for (i = 0; i < sn; i++) {
+ D(i) += fix_mul(SS_(i) + SS_(i + 1), 434);
+ }
+ for (i = 0; i < dn; i++) {
+ S(i) += fix_mul(DD_(i) + DD_(i - 1), 12994); /* 12993 */
+ }
+ }
+ }
}
/* <summary> */
/* Get norm of arbitrary wavelet transform. */
/* </summary> */
-static int upandconv(double *nXPS, double *LPS, int lenXPS, int lenLPS) {
- /* Perform the convolution of the vectors. */
- int i,j;
- double *tmp = (double *)opj_malloc(2*lenXPS * sizeof(double));
- /*Upsample*/
- memset(tmp, 0, 2*lenXPS*sizeof(double));
- for (i = 0; i < lenXPS; i++) {
- *(tmp + 2*i) = *(nXPS + i);
- *(nXPS + i) = 0;
- }
- /*Convolution*/
- for (i = 0; i < 2*lenXPS; i++) {
- for (j = 0; j < lenLPS; j++) {
- *(nXPS+i+j) = *(nXPS+i+j) + *(tmp + i) * *(LPS + j);
- /*fprintf(stdout,"*(tmp + %d) * *(LPS + %d) = %f * %f \n",i,j,*(tmp + i),*(LPS + j));*/
- }
- }
- free(tmp);
- return 2*lenXPS+lenLPS-1;
+static int upandconv(double *nXPS, double *LPS, int lenXPS, int lenLPS)
+{
+ /* Perform the convolution of the vectors. */
+ int i, j;
+ double *tmp = (double *)opj_malloc(2 * lenXPS * sizeof(double));
+ /*Upsample*/
+ memset(tmp, 0, 2 * lenXPS * sizeof(double));
+ for (i = 0; i < lenXPS; i++) {
+ *(tmp + 2 * i) = *(nXPS + i);
+ *(nXPS + i) = 0;
+ }
+ /*Convolution*/
+ for (i = 0; i < 2 * lenXPS; i++) {
+ for (j = 0; j < lenLPS; j++) {
+ *(nXPS + i + j) = *(nXPS + i + j) + *(tmp + i) * *(LPS + j);
+ /*fprintf(stdout,"*(tmp + %d) * *(LPS + %d) = %f * %f \n",i,j,*(tmp + i),*(LPS + j));*/
+ }
+ }
+ free(tmp);
+ return 2 * lenXPS + lenLPS - 1;
}
-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) {
- int i, lenLPS, lenHPS;
- double Lx = 0, Ly= 0, Hx= 0, Hy= 0, Lz= 0, Hz= 0;
- double *nLPSx, *nHPSx,*nLPSy, *nHPSy,*nLPSz, *nHPSz;
- int levelx, levely, levelz;
-
- levelx = (orient == 0) ? level[0]-1 : level[0];
- levely = (orient == 0) ? level[1]-1 : level[1];
- levelz = (orient == 0) ? level[2]-1 : level[2];
-
- /*X axis*/
- lenLPS = wtfiltX->lenLPS;
- lenHPS = wtfiltX->lenHPS;
- for (i = 0; i < levelx; i++) {
- lenLPS *= 2;
- lenHPS *= 2;
- lenLPS += wtfiltX->lenLPS - 1;
- lenHPS += wtfiltX->lenLPS - 1;
- }
- nLPSx = (double *)opj_malloc(lenLPS * sizeof(double));
- nHPSx = (double *)opj_malloc(lenHPS * sizeof(double));
-
- memcpy(nLPSx, wtfiltX->LPS, wtfiltX->lenLPS * sizeof(double));
- memcpy(nHPSx, wtfiltX->HPS, wtfiltX->lenHPS * sizeof(double));
- lenLPS = wtfiltX->lenLPS;
- lenHPS = wtfiltX->lenHPS;
- for (i = 0; i < levelx; i++) {
- lenLPS = upandconv(nLPSx, wtfiltX->LPS, lenLPS, wtfiltX->lenLPS);
- lenHPS = upandconv(nHPSx, wtfiltX->LPS, lenHPS, wtfiltX->lenLPS);
- }
- for (i = 0; i < lenLPS; i++)
- Lx += nLPSx[i] * nLPSx[i];
- for (i = 0; i < lenHPS; i++)
- Hx += nHPSx[i] * nHPSx[i];
- Lx = sqrt(Lx);
- Hx = sqrt(Hx);
- free(nLPSx);
- free(nHPSx);
-
- /*Y axis*/
- if (dwtid[0] != dwtid[1] || level[0] != level[1]){
- lenLPS = wtfiltY->lenLPS;
- lenHPS = wtfiltY->lenHPS;
- for (i = 0; i < levely; i++) {
- lenLPS *= 2;
- lenHPS *= 2;
- lenLPS += wtfiltY->lenLPS - 1;
- lenHPS += wtfiltY->lenLPS - 1;
- }
- nLPSy = (double *)opj_malloc(lenLPS * sizeof(double));
- nHPSy = (double *)opj_malloc(lenHPS * sizeof(double));
-
- memcpy(nLPSy, wtfiltY->LPS, wtfiltY->lenLPS * sizeof(double));
- memcpy(nHPSy, wtfiltY->HPS, wtfiltY->lenHPS * sizeof(double));
- lenLPS = wtfiltY->lenLPS;
- lenHPS = wtfiltY->lenHPS;
- for (i = 0; i < levely; i++) {
- lenLPS = upandconv(nLPSy, wtfiltY->LPS, lenLPS, wtfiltY->lenLPS);
- lenHPS = upandconv(nHPSy, wtfiltY->LPS, lenHPS, wtfiltY->lenLPS);
- }
- for (i = 0; i < lenLPS; i++)
- Ly += nLPSy[i] * nLPSy[i];
- for (i = 0; i < lenHPS; i++)
- Hy += nHPSy[i] * nHPSy[i];
- Ly = sqrt(Ly);
- Hy = sqrt(Hy);
- free(nLPSy);
- free(nHPSy);
- } else {
- Ly = Lx;
- Hy = Hx;
- }
- /*Z axis*/
- if (levelz >= 0) {
- lenLPS = wtfiltZ->lenLPS;
- lenHPS = wtfiltZ->lenHPS;
- for (i = 0; i < levelz; i++) {
- lenLPS *= 2;
- lenHPS *= 2;
- lenLPS += wtfiltZ->lenLPS - 1;
- lenHPS += wtfiltZ->lenLPS - 1;
- }
- nLPSz = (double *)opj_malloc(lenLPS * sizeof(double));
- nHPSz = (double *)opj_malloc(lenHPS * sizeof(double));
-
- memcpy(nLPSz, wtfiltZ->LPS, wtfiltZ->lenLPS * sizeof(double));
- memcpy(nHPSz, wtfiltZ->HPS, wtfiltZ->lenHPS * sizeof(double));
- lenLPS = wtfiltZ->lenLPS;
- lenHPS = wtfiltZ->lenHPS;
- for (i = 0; i < levelz; i++) {
- lenLPS = upandconv(nLPSz, wtfiltZ->LPS, lenLPS, wtfiltZ->lenLPS);
- lenHPS = upandconv(nHPSz, wtfiltZ->LPS, lenHPS, wtfiltZ->lenLPS);
- }
- for (i = 0; i < lenLPS; i++)
- Lz += nLPSz[i] * nLPSz[i];
- for (i = 0; i < lenHPS; i++)
- Hz += nHPSz[i] * nHPSz[i];
- Lz = sqrt(Lz);
- Hz = sqrt(Hz);
- free(nLPSz);
- free(nHPSz);
- } else {
- Lz = 1.0; Hz = 1.0;
- }
- switch (orient) {
- case 0:
- return Lx * Ly * Lz;
- case 1:
- return Lx * Hy * Lz;
- case 2:
- return Hx * Ly * Lz;
- case 3:
- return Hx * Hy * Lz;
- case 4:
- return Lx * Ly * Hz;
- case 5:
- return Lx * Hy * Hz;
- case 6:
- return Hx * Ly * Hz;
- case 7:
- return Hx * Hy * Hz;
- default:
- return -1;
- }
-
+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)
+{
+ int i, lenLPS, lenHPS;
+ double Lx = 0, Ly = 0, Hx = 0, Hy = 0, Lz = 0, Hz = 0;
+ double *nLPSx, *nHPSx, *nLPSy, *nHPSy, *nLPSz, *nHPSz;
+ int levelx, levely, levelz;
+
+ levelx = (orient == 0) ? level[0] - 1 : level[0];
+ levely = (orient == 0) ? level[1] - 1 : level[1];
+ levelz = (orient == 0) ? level[2] - 1 : level[2];
+
+ /*X axis*/
+ lenLPS = wtfiltX->lenLPS;
+ lenHPS = wtfiltX->lenHPS;
+ for (i = 0; i < levelx; i++) {
+ lenLPS *= 2;
+ lenHPS *= 2;
+ lenLPS += wtfiltX->lenLPS - 1;
+ lenHPS += wtfiltX->lenLPS - 1;
+ }
+ nLPSx = (double *)opj_malloc(lenLPS * sizeof(double));
+ nHPSx = (double *)opj_malloc(lenHPS * sizeof(double));
+
+ memcpy(nLPSx, wtfiltX->LPS, wtfiltX->lenLPS * sizeof(double));
+ memcpy(nHPSx, wtfiltX->HPS, wtfiltX->lenHPS * sizeof(double));
+ lenLPS = wtfiltX->lenLPS;
+ lenHPS = wtfiltX->lenHPS;
+ for (i = 0; i < levelx; i++) {
+ lenLPS = upandconv(nLPSx, wtfiltX->LPS, lenLPS, wtfiltX->lenLPS);
+ lenHPS = upandconv(nHPSx, wtfiltX->LPS, lenHPS, wtfiltX->lenLPS);
+ }
+ for (i = 0; i < lenLPS; i++) {
+ Lx += nLPSx[i] * nLPSx[i];
+ }
+ for (i = 0; i < lenHPS; i++) {
+ Hx += nHPSx[i] * nHPSx[i];
+ }
+ Lx = sqrt(Lx);
+ Hx = sqrt(Hx);
+ free(nLPSx);
+ free(nHPSx);
+
+ /*Y axis*/
+ if (dwtid[0] != dwtid[1] || level[0] != level[1]) {
+ lenLPS = wtfiltY->lenLPS;
+ lenHPS = wtfiltY->lenHPS;
+ for (i = 0; i < levely; i++) {
+ lenLPS *= 2;
+ lenHPS *= 2;
+ lenLPS += wtfiltY->lenLPS - 1;
+ lenHPS += wtfiltY->lenLPS - 1;
+ }
+ nLPSy = (double *)opj_malloc(lenLPS * sizeof(double));
+ nHPSy = (double *)opj_malloc(lenHPS * sizeof(double));
+
+ memcpy(nLPSy, wtfiltY->LPS, wtfiltY->lenLPS * sizeof(double));
+ memcpy(nHPSy, wtfiltY->HPS, wtfiltY->lenHPS * sizeof(double));
+ lenLPS = wtfiltY->lenLPS;
+ lenHPS = wtfiltY->lenHPS;
+ for (i = 0; i < levely; i++) {
+ lenLPS = upandconv(nLPSy, wtfiltY->LPS, lenLPS, wtfiltY->lenLPS);
+ lenHPS = upandconv(nHPSy, wtfiltY->LPS, lenHPS, wtfiltY->lenLPS);
+ }
+ for (i = 0; i < lenLPS; i++) {
+ Ly += nLPSy[i] * nLPSy[i];
+ }
+ for (i = 0; i < lenHPS; i++) {
+ Hy += nHPSy[i] * nHPSy[i];
+ }
+ Ly = sqrt(Ly);
+ Hy = sqrt(Hy);
+ free(nLPSy);
+ free(nHPSy);
+ } else {
+ Ly = Lx;
+ Hy = Hx;
+ }
+ /*Z axis*/
+ if (levelz >= 0) {
+ lenLPS = wtfiltZ->lenLPS;
+ lenHPS = wtfiltZ->lenHPS;
+ for (i = 0; i < levelz; i++) {
+ lenLPS *= 2;
+ lenHPS *= 2;
+ lenLPS += wtfiltZ->lenLPS - 1;
+ lenHPS += wtfiltZ->lenLPS - 1;
+ }
+ nLPSz = (double *)opj_malloc(lenLPS * sizeof(double));
+ nHPSz = (double *)opj_malloc(lenHPS * sizeof(double));
+
+ memcpy(nLPSz, wtfiltZ->LPS, wtfiltZ->lenLPS * sizeof(double));
+ memcpy(nHPSz, wtfiltZ->HPS, wtfiltZ->lenHPS * sizeof(double));
+ lenLPS = wtfiltZ->lenLPS;
+ lenHPS = wtfiltZ->lenHPS;
+ for (i = 0; i < levelz; i++) {
+ lenLPS = upandconv(nLPSz, wtfiltZ->LPS, lenLPS, wtfiltZ->lenLPS);
+ lenHPS = upandconv(nHPSz, wtfiltZ->LPS, lenHPS, wtfiltZ->lenLPS);
+ }
+ for (i = 0; i < lenLPS; i++) {
+ Lz += nLPSz[i] * nLPSz[i];
+ }
+ for (i = 0; i < lenHPS; i++) {
+ Hz += nHPSz[i] * nHPSz[i];
+ }
+ Lz = sqrt(Lz);
+ Hz = sqrt(Hz);
+ free(nLPSz);
+ free(nHPSz);
+ } else {
+ Lz = 1.0;
+ Hz = 1.0;
+ }
+ switch (orient) {
+ case 0:
+ return Lx * Ly * Lz;
+ case 1:
+ return Lx * Hy * Lz;
+ case 2:
+ return Hx * Ly * Lz;
+ case 3:
+ return Hx * Hy * Lz;
+ case 4:
+ return Lx * Ly * Hz;
+ case 5:
+ return Lx * Hy * Hz;
+ case 6:
+ return Hx * Ly * Hz;
+ case 7:
+ return Hx * Hy * Hz;
+ default:
+ return -1;
+ }
+
}
-static void dwt_getwtfilters(opj_wtfilt_t *wtfilt, int dwtid) {
- if (dwtid == 0) { /*DWT 9-7 */
- wtfilt->lenLPS = 7; wtfilt->lenHPS = 9;
- wtfilt->LPS = (double *)opj_malloc(wtfilt->lenLPS * sizeof(double));
- wtfilt->HPS = (double *)opj_malloc(wtfilt->lenHPS * sizeof(double));
- wtfilt->LPS[0] = -0.091271763114; wtfilt->HPS[0] = 0.026748757411;
- wtfilt->LPS[1] = -0.057543526228; wtfilt->HPS[1] = 0.016864118443;
- wtfilt->LPS[2] = 0.591271763114; wtfilt->HPS[2] = -0.078223266529;
- wtfilt->LPS[3] = 1.115087052457; wtfilt->HPS[3] = -0.266864118443;
- wtfilt->LPS[4] = 0.591271763114; wtfilt->HPS[4] = 0.602949018236;
- wtfilt->LPS[5] = -0.057543526228; wtfilt->HPS[5] = -0.266864118443;
- wtfilt->LPS[6] = -0.091271763114; wtfilt->HPS[6] = -0.078223266529;
- wtfilt->HPS[7] = 0.016864118443;
- wtfilt->HPS[8] = 0.026748757411;
- } else if (dwtid == 1) { /*DWT 5-3 */
- wtfilt->lenLPS = 3; wtfilt->lenHPS = 5;
- wtfilt->LPS = (double *)opj_malloc(wtfilt->lenLPS * sizeof(double));
- wtfilt->HPS = (double *)opj_malloc(wtfilt->lenHPS * sizeof(double));
- wtfilt->LPS[0] = 0.5; wtfilt->HPS[0] = -0.125;
- wtfilt->LPS[1] = 1; wtfilt->HPS[1] = -0.25;
- wtfilt->LPS[2] = 0.5; wtfilt->HPS[2] = 0.75;
- wtfilt->HPS[3] = -0.25;
- wtfilt->HPS[4] = -0.125;
- } else {
- fprintf(stdout,"[ERROR] Sorry, this wavelet hasn't been implemented so far ... Try another one :-)\n");
- exit(1);
- }
+static void dwt_getwtfilters(opj_wtfilt_t *wtfilt, int dwtid)
+{
+ if (dwtid == 0) { /*DWT 9-7 */
+ wtfilt->lenLPS = 7;
+ wtfilt->lenHPS = 9;
+ wtfilt->LPS = (double *)opj_malloc(wtfilt->lenLPS * sizeof(double));
+ wtfilt->HPS = (double *)opj_malloc(wtfilt->lenHPS * sizeof(double));
+ wtfilt->LPS[0] = -0.091271763114;
+ wtfilt->HPS[0] = 0.026748757411;
+ wtfilt->LPS[1] = -0.057543526228;
+ wtfilt->HPS[1] = 0.016864118443;
+ wtfilt->LPS[2] = 0.591271763114;
+ wtfilt->HPS[2] = -0.078223266529;
+ wtfilt->LPS[3] = 1.115087052457;
+ wtfilt->HPS[3] = -0.266864118443;
+ wtfilt->LPS[4] = 0.591271763114;
+ wtfilt->HPS[4] = 0.602949018236;
+ wtfilt->LPS[5] = -0.057543526228;
+ wtfilt->HPS[5] = -0.266864118443;
+ wtfilt->LPS[6] = -0.091271763114;
+ wtfilt->HPS[6] = -0.078223266529;
+ wtfilt->HPS[7] = 0.016864118443;
+ wtfilt->HPS[8] = 0.026748757411;
+ } else if (dwtid == 1) { /*DWT 5-3 */
+ wtfilt->lenLPS = 3;
+ wtfilt->lenHPS = 5;
+ wtfilt->LPS = (double *)opj_malloc(wtfilt->lenLPS * sizeof(double));
+ wtfilt->HPS = (double *)opj_malloc(wtfilt->lenHPS * sizeof(double));
+ wtfilt->LPS[0] = 0.5;
+ wtfilt->HPS[0] = -0.125;
+ wtfilt->LPS[1] = 1;
+ wtfilt->HPS[1] = -0.25;
+ wtfilt->LPS[2] = 0.5;
+ wtfilt->HPS[2] = 0.75;
+ wtfilt->HPS[3] = -0.25;
+ wtfilt->HPS[4] = -0.125;
+ } else {
+ fprintf(stdout,
+ "[ERROR] Sorry, this wavelet hasn't been implemented so far ... Try another one :-)\n");
+ exit(1);
+ }
}
/* <summary> */
/* Encoding of quantization stepsize for each subband. */
-/* </summary> */
-static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno_stepsize) {
- int p, n;
- p = int_floorlog2(stepsize) - 13;
- n = 11 - int_floorlog2(stepsize);
- bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
- bandno_stepsize->expn = numbps - p;
- /*if J3D_CCP_QNTSTY_NOQNT --> stepsize = 8192.0 --> p = 0, n = -2 --> mant = 0; expn = (prec+gain)*/
- /*else --> bandno_stepsize = (1<<(numbps - expn)) + (1<<(numbps - expn - 11)) * Ub*/
+/* </summary> */
+static void dwt_encode_stepsize(int stepsize, int numbps,
+ opj_stepsize_t *bandno_stepsize)
+{
+ int p, n;
+ p = int_floorlog2(stepsize) - 13;
+ n = 11 - int_floorlog2(stepsize);
+ bandno_stepsize->mant = (n < 0 ? stepsize >> -n : stepsize << n) & 0x7ff;
+ bandno_stepsize->expn = numbps - p;
+ /*if J3D_CCP_QNTSTY_NOQNT --> stepsize = 8192.0 --> p = 0, n = -2 --> mant = 0; expn = (prec+gain)*/
+ /*else --> bandno_stepsize = (1<<(numbps - expn)) + (1<<(numbps - expn - 11)) * Ub*/
}
-/*
+/*
==========================================================
DWT interface
==========================================================
@@ -663,261 +748,298 @@ static void dwt_encode_stepsize(int stepsize, int numbps, opj_stepsize_t *bandno
/* <summary> */
/* Forward 5-3 wavelet transform in 3-D. */
/* </summary> */
-void dwt_encode(opj_tcd_tilecomp_t * tilec, int dwtid[3]) {
- int i, j, k;
- int x, y, z;
- int w, h, wh, d;
- int level,levelx,levely,levelz,diff;
- int *a = NULL;
- int *aj = NULL;
- int *bj = NULL;
- int *cj = NULL;
-
- /*ops = 0;*/
-
- memset(flagnorm,0,8000*sizeof(int));
- w = tilec->x1-tilec->x0;
- h = tilec->y1-tilec->y0;
- d = tilec->z1-tilec->z0;
- wh = w * h;
- levelx = tilec->numresolution[0]-1;
- levely = tilec->numresolution[1]-1;
- levelz = tilec->numresolution[2]-1;
- level = int_max(levelx,int_max(levely,levelz));
- diff = tilec->numresolution[0] - tilec->numresolution[2];
-
- a = tilec->data;
-
- for (x = 0, y = 0, z = 0; (x < levelx) && (y < levely); x++, y++, z++) {
- int rw; /* width of the resolution level computed */
- int rh; /* heigth of the resolution level computed */
- int rd; /* depth of the resolution level computed */
- int rw1; /* width of the resolution level once lower than computed one */
- int rh1; /* height of the resolution level once lower than computed one */
- int rd1; /* depth of the resolution level once lower than computed one */
- int cas_col; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
- int cas_row; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */
- int cas_axl; /* 0 = non inversion on axial filtering 1 = inversion between low-pass and high-pass filtering */
- int dn, sn;
-
- rw = tilec->resolutions[level - x].x1 - tilec->resolutions[level - x].x0;
- rh = tilec->resolutions[level - y].y1 - tilec->resolutions[level - y].y0;
- rd = tilec->resolutions[level - z].z1 - tilec->resolutions[level - z].z0;
- rw1= tilec->resolutions[level - x - 1].x1 - tilec->resolutions[level - x - 1].x0;
- rh1= tilec->resolutions[level - y - 1].y1 - tilec->resolutions[level - y - 1].y0;
- rd1= tilec->resolutions[level - z - 1].z1 - tilec->resolutions[level - z - 1].z0;
-
- cas_col = tilec->resolutions[level - x].x0 % 2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
- cas_row = tilec->resolutions[level - y].y0 % 2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */
- cas_axl = tilec->resolutions[level - z].z0 % 2;
-
- /*fprintf(stdout," x %d y %d z %d \n",x,y,z);
- fprintf(stdout," levelx %d levely %d levelz %d \n",levelx,levely,levelz);
- fprintf(stdout," z1 %d z0 %d\n",tilec->resolutions[level - z].z1,tilec->resolutions[level - z].z0);
- fprintf(stdout," rw %d rh %d rd %d \n rw1 %d rh1 %d rd1 %d \n",rw,rh,rd,rw1,rh1,rd1);*/
-
- for (i = 0; i < rd; i++) {
-
- cj = a + (i * wh);
-
- /*Horizontal*/
- sn = rw1;
- dn = rw - rw1;
- bj = (int*)opj_malloc(rw * sizeof(int));
- if (dwtid[0] == 0) {
- for (j = 0; j < rh; j++) {
- aj = cj + j * w;
- for (k = 0; k < rw; k++) bj[k] = aj[k];
- dwt_encode_97(bj, dn, sn, cas_row);
- dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
- }
- } else if (dwtid[0] == 1) {
- for (j = 0; j < rh; j++) {
- aj = cj + j * w;
- for (k = 0; k < rw; k++) bj[k] = aj[k];
- dwt_encode_53(bj, dn, sn, cas_row);
- dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
- }
- }
- opj_free(bj);
-
- /*Vertical*/
- sn = rh1;
- dn = rh - rh1;
- bj = (int*)opj_malloc(rh * sizeof(int));
- if (dwtid[1] == 0) { /*DWT 9-7*/
- for (j = 0; j < rw; j++) {
- aj = cj + j;
- for (k = 0; k < rh; k++) bj[k] = aj[k*w];
- dwt_encode_97(bj, dn, sn, cas_col);
- dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
- }
+void dwt_encode(opj_tcd_tilecomp_t * tilec, int dwtid[3])
+{
+ int i, j, k;
+ int x, y, z;
+ int w, h, wh, d;
+ int level, levelx, levely, levelz, diff;
+ int *a = NULL;
+ int *aj = NULL;
+ int *bj = NULL;
+ int *cj = NULL;
+
+ /*ops = 0;*/
+
+ memset(flagnorm, 0, 8000 * sizeof(int));
+ w = tilec->x1 - tilec->x0;
+ h = tilec->y1 - tilec->y0;
+ d = tilec->z1 - tilec->z0;
+ wh = w * h;
+ levelx = tilec->numresolution[0] - 1;
+ levely = tilec->numresolution[1] - 1;
+ levelz = tilec->numresolution[2] - 1;
+ level = int_max(levelx, int_max(levely, levelz));
+ diff = tilec->numresolution[0] - tilec->numresolution[2];
+
+ a = tilec->data;
+
+ for (x = 0, y = 0, z = 0; (x < levelx) && (y < levely); x++, y++, z++) {
+ int rw; /* width of the resolution level computed */
+ int rh; /* heigth of the resolution level computed */
+ int rd; /* depth of the resolution level computed */
+ int rw1; /* width of the resolution level once lower than computed one */
+ int rh1; /* height of the resolution level once lower than computed one */
+ int rd1; /* depth of the resolution level once lower than computed one */
+ int cas_col; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
+ int cas_row; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */
+ int cas_axl; /* 0 = non inversion on axial filtering 1 = inversion between low-pass and high-pass filtering */
+ int dn, sn;
+
+ rw = tilec->resolutions[level - x].x1 - tilec->resolutions[level - x].x0;
+ rh = tilec->resolutions[level - y].y1 - tilec->resolutions[level - y].y0;
+ rd = tilec->resolutions[level - z].z1 - tilec->resolutions[level - z].z0;
+ rw1 = tilec->resolutions[level - x - 1].x1 - tilec->resolutions[level - x -
+ 1].x0;
+ rh1 = tilec->resolutions[level - y - 1].y1 - tilec->resolutions[level - y -
+ 1].y0;
+ rd1 = tilec->resolutions[level - z - 1].z1 - tilec->resolutions[level - z -
+ 1].z0;
+
+ cas_col = tilec->resolutions[level - x].x0 %
+ 2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
+ cas_row = tilec->resolutions[level - y].y0 %
+ 2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */
+ cas_axl = tilec->resolutions[level - z].z0 % 2;
+
+ /*fprintf(stdout," x %d y %d z %d \n",x,y,z);
+ fprintf(stdout," levelx %d levely %d levelz %d \n",levelx,levely,levelz);
+ fprintf(stdout," z1 %d z0 %d\n",tilec->resolutions[level - z].z1,tilec->resolutions[level - z].z0);
+ fprintf(stdout," rw %d rh %d rd %d \n rw1 %d rh1 %d rd1 %d \n",rw,rh,rd,rw1,rh1,rd1);*/
+
+ for (i = 0; i < rd; i++) {
+
+ cj = a + (i * wh);
+
+ /*Horizontal*/
+ sn = rw1;
+ dn = rw - rw1;
+ bj = (int*)opj_malloc(rw * sizeof(int));
+ if (dwtid[0] == 0) {
+ for (j = 0; j < rh; j++) {
+ aj = cj + j * w;
+ for (k = 0; k < rw; k++) {
+ bj[k] = aj[k];
+ }
+ dwt_encode_97(bj, dn, sn, cas_row);
+ dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
+ }
+ } else if (dwtid[0] == 1) {
+ for (j = 0; j < rh; j++) {
+ aj = cj + j * w;
+ for (k = 0; k < rw; k++) {
+ bj[k] = aj[k];
+ }
+ dwt_encode_53(bj, dn, sn, cas_row);
+ dwt_deinterleave_h(bj, aj, dn, sn, cas_row);
+ }
+ }
+ opj_free(bj);
+
+ /*Vertical*/
+ sn = rh1;
+ dn = rh - rh1;
+ bj = (int*)opj_malloc(rh * sizeof(int));
+ if (dwtid[1] == 0) { /*DWT 9-7*/
+ for (j = 0; j < rw; j++) {
+ aj = cj + j;
+ for (k = 0; k < rh; k++) {
+ bj[k] = aj[k * w];
+ }
+ dwt_encode_97(bj, dn, sn, cas_col);
+ dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
+ }
} else if (dwtid[1] == 1) { /*DWT 5-3*/
- for (j = 0; j < rw; j++) {
- aj = cj + j;
- for (k = 0; k < rh; k++) bj[k] = aj[k*w];
- dwt_encode_53(bj, dn, sn, cas_col);
- dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
- }
- }
- opj_free(bj);
- }
-
- if (z < levelz){
- /*Axial fprintf(stdout,"Axial DWT Transform %d %d %d\n",z,rd,rd1);*/
- sn = rd1;
- dn = rd - rd1;
- bj = (int*)opj_malloc(rd * sizeof(int));
- if (dwtid[2] == 0) {
- for (j = 0; j < (rw*rh); j++) {
- aj = a + j;
- for (k = 0; k < rd; k++) bj[k] = aj[k*wh];
- dwt_encode_97(bj, dn, sn, cas_axl);
- dwt_deinterleave_z(bj, aj, dn, sn, wh, cas_axl);
- }
- } else if (dwtid[2] == 1) {
- for (j = 0; j < (rw*rh); j++) {
- aj = a + j;
- for (k = 0; k < rd; k++) bj[k] = aj[k*wh];
- dwt_encode_53(bj, dn, sn, cas_axl);
- dwt_deinterleave_z(bj, aj, dn, sn, wh, cas_axl);
- }
- }
- opj_free(bj);
- }
- }
-
- /*fprintf(stdout,"[INFO] Ops: %d \n",ops);*/
+ for (j = 0; j < rw; j++) {
+ aj = cj + j;
+ for (k = 0; k < rh; k++) {
+ bj[k] = aj[k * w];
+ }
+ dwt_encode_53(bj, dn, sn, cas_col);
+ dwt_deinterleave_v(bj, aj, dn, sn, w, cas_col);
+ }
+ }
+ opj_free(bj);
+ }
+
+ if (z < levelz) {
+ /*Axial fprintf(stdout,"Axial DWT Transform %d %d %d\n",z,rd,rd1);*/
+ sn = rd1;
+ dn = rd - rd1;
+ bj = (int*)opj_malloc(rd * sizeof(int));
+ if (dwtid[2] == 0) {
+ for (j = 0; j < (rw * rh); j++) {
+ aj = a + j;
+ for (k = 0; k < rd; k++) {
+ bj[k] = aj[k * wh];
+ }
+ dwt_encode_97(bj, dn, sn, cas_axl);
+ dwt_deinterleave_z(bj, aj, dn, sn, wh, cas_axl);
+ }
+ } else if (dwtid[2] == 1) {
+ for (j = 0; j < (rw * rh); j++) {
+ aj = a + j;
+ for (k = 0; k < rd; k++) {
+ bj[k] = aj[k * wh];
+ }
+ dwt_encode_53(bj, dn, sn, cas_axl);
+ dwt_deinterleave_z(bj, aj, dn, sn, wh, cas_axl);
+ }
+ }
+ opj_free(bj);
+ }
+ }
+
+ /*fprintf(stdout,"[INFO] Ops: %d \n",ops);*/
}
/* <summary> */
/* Inverse 5-3 wavelet transform in 3-D. */
/* </summary> */
-void dwt_decode(opj_tcd_tilecomp_t * tilec, int stops[3], int dwtid[3]) {
- int i, j, k;
- int x, y, z;
- int w, h, wh, d;
- int level, levelx, levely, levelz, diff;
- int *a = NULL;
- int *aj = NULL;
- int *bj = NULL;
- int *cj = NULL;
-
- a = tilec->data;
-
- w = tilec->x1-tilec->x0;
- h = tilec->y1-tilec->y0;
- d = tilec->z1-tilec->z0;
- wh = w * h;
- levelx = tilec->numresolution[0]-1;
- levely = tilec->numresolution[1]-1;
- levelz = tilec->numresolution[2]-1;
- level = int_max(levelx,int_max(levely,levelz));
- diff = tilec->numresolution[0] - tilec->numresolution[2];
-
-/* General lifting framework -- DCCS-LIWT */
- for (x = level - 1, y = level - 1, z = level - 1; (x >= stops[0]) && (y >= stops[1]); x--, y--, z--) {
- int rw; /* width of the resolution level computed */
- int rh; /* heigth of the resolution level computed */
- int rd; /* depth of the resolution level computed */
- int rw1; /* width of the resolution level once lower than computed one */
- int rh1; /* height of the resolution level once lower than computed one */
- int rd1; /* depth of the resolution level once lower than computed one */
- int cas_col; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
- int cas_row; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */
- int cas_axl; /* 0 = non inversion on axial filtering 1 = inversion between low-pass and high-pass filtering */
- int dn, sn;
-
- rw = tilec->resolutions[level - x].x1 - tilec->resolutions[level - x].x0;
- rh = tilec->resolutions[level - y].y1 - tilec->resolutions[level - y].y0;
- rd = tilec->resolutions[level - z].z1 - tilec->resolutions[level - z].z0;
- rw1= tilec->resolutions[level - x - 1].x1 - tilec->resolutions[level - x - 1].x0;
- rh1= tilec->resolutions[level - y - 1].y1 - tilec->resolutions[level - y - 1].y0;
- rd1= tilec->resolutions[level - z - 1].z1 - tilec->resolutions[level - z - 1].z0;
-
- cas_col = tilec->resolutions[level - x].x0 % 2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
- cas_row = tilec->resolutions[level - y].y0 % 2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */
- cas_axl = tilec->resolutions[level - z].z0 % 2;
-
- /*fprintf(stdout," x %d y %d z %d \n",x,y,z);
- fprintf(stdout," levelx %d levely %d levelz %d \n",levelx,levely,levelz);
- fprintf(stdout," dwtid[0] %d [1] %d [2] %d \n",dwtid[0],dwtid[1],dwtid[2]);
- fprintf(stdout," rw %d rh %d rd %d \n rw1 %d rh1 %d rd1 %d \n",rw,rh,rd,rw1,rh1,rd1);
- fprintf(stdout,"IDWT Transform %d %d %d %d\n",level, z, rd,rd1);*/
-
- if (z >= stops[2] && rd != rd1) {
- /*fprintf(stdout,"Axial Transform %d %d %d %d\n",levelz, z, rd,rd1);*/
- sn = rd1;
- dn = rd - rd1;
- bj = (int*)opj_malloc(rd * sizeof(int));
- if (dwtid[2] == 0) {
- for (j = 0; j < (rw*rh); j++) {
- aj = a + j;
- dwt_interleave_z(aj, bj, dn, sn, wh, cas_axl);
- dwt_decode_97(bj, dn, sn, cas_axl);
- for (k = 0; k < rd; k++) aj[k * wh] = bj[k];
- }
- } else if (dwtid[2] == 1) {
- for (j = 0; j < (rw*rh); j++) {
- aj = a + j;
- dwt_interleave_z(aj, bj, dn, sn, wh, cas_axl);
- dwt_decode_53(bj, dn, sn, cas_axl);
- for (k = 0; k < rd; k++) aj[k * wh] = bj[k];
- }
- }
- opj_free(bj);
- }
-
- for (i = 0; i < rd; i++) {
- /*Fetch corresponding slice for doing DWT-2D*/
- cj = tilec->data + (i * wh);
-
- /*Vertical*/
- sn = rh1;
- dn = rh - rh1;
- bj = (int*)opj_malloc(rh * sizeof(int));
- if (dwtid[1] == 0) {
- for (j = 0; j < rw; j++) {
- aj = cj + j;
- dwt_interleave_v(aj, bj, dn, sn, w, cas_col);
- dwt_decode_97(bj, dn, sn, cas_col);
- for (k = 0; k < rh; k++) aj[k * w] = bj[k];
- }
- } else if (dwtid[1] == 1) {
- for (j = 0; j < rw; j++) {
- aj = cj + j;
- dwt_interleave_v(aj, bj, dn, sn, w, cas_col);
- dwt_decode_53(bj, dn, sn, cas_col);
- for (k = 0; k < rh; k++) aj[k * w] = bj[k];
- }
- }
- opj_free(bj);
-
- /*Horizontal*/
- sn = rw1;
- dn = rw - rw1;
- bj = (int*)opj_malloc(rw * sizeof(int));
- if (dwtid[0]==0) {
- for (j = 0; j < rh; j++) {
- aj = cj + j*w;
- dwt_interleave_h(aj, bj, dn, sn, cas_row);
- dwt_decode_97(bj, dn, sn, cas_row);
- for (k = 0; k < rw; k++) aj[k] = bj[k];
- }
- } else if (dwtid[0]==1) {
- for (j = 0; j < rh; j++) {
- aj = cj + j*w;
- dwt_interleave_h(aj, bj, dn, sn, cas_row);
- dwt_decode_53(bj, dn, sn, cas_row);
- for (k = 0; k < rw; k++) aj[k] = bj[k];
- }
- }
- opj_free(bj);
-
- }
-
- }
+void dwt_decode(opj_tcd_tilecomp_t * tilec, int stops[3], int dwtid[3])
+{
+ int i, j, k;
+ int x, y, z;
+ int w, h, wh, d;
+ int level, levelx, levely, levelz, diff;
+ int *a = NULL;
+ int *aj = NULL;
+ int *bj = NULL;
+ int *cj = NULL;
+
+ a = tilec->data;
+
+ w = tilec->x1 - tilec->x0;
+ h = tilec->y1 - tilec->y0;
+ d = tilec->z1 - tilec->z0;
+ wh = w * h;
+ levelx = tilec->numresolution[0] - 1;
+ levely = tilec->numresolution[1] - 1;
+ levelz = tilec->numresolution[2] - 1;
+ level = int_max(levelx, int_max(levely, levelz));
+ diff = tilec->numresolution[0] - tilec->numresolution[2];
+
+ /* General lifting framework -- DCCS-LIWT */
+ for (x = level - 1, y = level - 1, z = level - 1; (x >= stops[0]) &&
+ (y >= stops[1]); x--, y--, z--) {
+ int rw; /* width of the resolution level computed */
+ int rh; /* heigth of the resolution level computed */
+ int rd; /* depth of the resolution level computed */
+ int rw1; /* width of the resolution level once lower than computed one */
+ int rh1; /* height of the resolution level once lower than computed one */
+ int rd1; /* depth of the resolution level once lower than computed one */
+ int cas_col; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
+ int cas_row; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */
+ int cas_axl; /* 0 = non inversion on axial filtering 1 = inversion between low-pass and high-pass filtering */
+ int dn, sn;
+
+ rw = tilec->resolutions[level - x].x1 - tilec->resolutions[level - x].x0;
+ rh = tilec->resolutions[level - y].y1 - tilec->resolutions[level - y].y0;
+ rd = tilec->resolutions[level - z].z1 - tilec->resolutions[level - z].z0;
+ rw1 = tilec->resolutions[level - x - 1].x1 - tilec->resolutions[level - x -
+ 1].x0;
+ rh1 = tilec->resolutions[level - y - 1].y1 - tilec->resolutions[level - y -
+ 1].y0;
+ rd1 = tilec->resolutions[level - z - 1].z1 - tilec->resolutions[level - z -
+ 1].z0;
+
+ cas_col = tilec->resolutions[level - x].x0 %
+ 2; /* 0 = non inversion on horizontal filtering 1 = inversion between low-pass and high-pass filtering */
+ cas_row = tilec->resolutions[level - y].y0 %
+ 2; /* 0 = non inversion on vertical filtering 1 = inversion between low-pass and high-pass filtering */
+ cas_axl = tilec->resolutions[level - z].z0 % 2;
+
+ /*fprintf(stdout," x %d y %d z %d \n",x,y,z);
+ fprintf(stdout," levelx %d levely %d levelz %d \n",levelx,levely,levelz);
+ fprintf(stdout," dwtid[0] %d [1] %d [2] %d \n",dwtid[0],dwtid[1],dwtid[2]);
+ fprintf(stdout," rw %d rh %d rd %d \n rw1 %d rh1 %d rd1 %d \n",rw,rh,rd,rw1,rh1,rd1);
+ fprintf(stdout,"IDWT Transform %d %d %d %d\n",level, z, rd,rd1);*/
+
+ if (z >= stops[2] && rd != rd1) {
+ /*fprintf(stdout,"Axial Transform %d %d %d %d\n",levelz, z, rd,rd1);*/
+ sn = rd1;
+ dn = rd - rd1;
+ bj = (int*)opj_malloc(rd * sizeof(int));
+ if (dwtid[2] == 0) {
+ for (j = 0; j < (rw * rh); j++) {
+ aj = a + j;
+ dwt_interleave_z(aj, bj, dn, sn, wh, cas_axl);
+ dwt_decode_97(bj, dn, sn, cas_axl);
+ for (k = 0; k < rd; k++) {
+ aj[k * wh] = bj[k];
+ }
+ }
+ } else if (dwtid[2] == 1) {
+ for (j = 0; j < (rw * rh); j++) {
+ aj = a + j;
+ dwt_interleave_z(aj, bj, dn, sn, wh, cas_axl);
+ dwt_decode_53(bj, dn, sn, cas_axl);
+ for (k = 0; k < rd; k++) {
+ aj[k * wh] = bj[k];
+ }
+ }
+ }
+ opj_free(bj);
+ }
+
+ for (i = 0; i < rd; i++) {
+ /*Fetch corresponding slice for doing DWT-2D*/
+ cj = tilec->data + (i * wh);
+
+ /*Vertical*/
+ sn = rh1;
+ dn = rh - rh1;
+ bj = (int*)opj_malloc(rh * sizeof(int));
+ if (dwtid[1] == 0) {
+ for (j = 0; j < rw; j++) {
+ aj = cj + j;
+ dwt_interleave_v(aj, bj, dn, sn, w, cas_col);
+ dwt_decode_97(bj, dn, sn, cas_col);
+ for (k = 0; k < rh; k++) {
+ aj[k * w] = bj[k];
+ }
+ }
+ } else if (dwtid[1] == 1) {
+ for (j = 0; j < rw; j++) {
+ aj = cj + j;
+ dwt_interleave_v(aj, bj, dn, sn, w, cas_col);
+ dwt_decode_53(bj, dn, sn, cas_col);
+ for (k = 0; k < rh; k++) {
+ aj[k * w] = bj[k];
+ }
+ }
+ }
+ opj_free(bj);
+
+ /*Horizontal*/
+ sn = rw1;
+ dn = rw - rw1;
+ bj = (int*)opj_malloc(rw * sizeof(int));
+ if (dwtid[0] == 0) {
+ for (j = 0; j < rh; j++) {
+ aj = cj + j * w;
+ dwt_interleave_h(aj, bj, dn, sn, cas_row);
+ dwt_decode_97(bj, dn, sn, cas_row);
+ for (k = 0; k < rw; k++) {
+ aj[k] = bj[k];
+ }
+ }
+ } else if (dwtid[0] == 1) {
+ for (j = 0; j < rh; j++) {
+ aj = cj + j * w;
+ dwt_interleave_h(aj, bj, dn, sn, cas_row);
+ dwt_decode_53(bj, dn, sn, cas_row);
+ for (k = 0; k < rw; k++) {
+ aj[k] = bj[k];
+ }
+ }
+ }
+ opj_free(bj);
+
+ }
+
+ }
}
@@ -925,96 +1047,111 @@ void dwt_decode(opj_tcd_tilecomp_t * tilec, int stops[3], int dwtid[3]) {
/* <summary> */
/* Get gain of wavelet transform. */
/* </summary> */
-int dwt_getgain(int orient, int reversible) {
- if (reversible == 1) {
- if (orient == 0)
- return 0;
- else if (orient == 1 || orient == 2 || orient == 4 )
- return 1;
- else if (orient == 3 || orient == 5 || orient == 6 )
- return 2;
- else
- return 3;
- }
- /*else if (reversible == 0){*/
- return 0;
+int dwt_getgain(int orient, int reversible)
+{
+ if (reversible == 1) {
+ if (orient == 0) {
+ return 0;
+ } else if (orient == 1 || orient == 2 || orient == 4) {
+ return 1;
+ } else if (orient == 3 || orient == 5 || orient == 6) {
+ return 2;
+ } else {
+ return 3;
+ }
+ }
+ /*else if (reversible == 0){*/
+ return 0;
}
/* <summary> */
/* Get norm of wavelet transform. */
/* </summary> */
-double dwt_getnorm(int orient, int level[3], int dwtid[3]) {
- int levelx = level[0];
- int levely = level[1];
- int levelz = (level[2] < 0) ? 0 : level[2];
- double norm;
-
- if (flagnorm[levelx][levely][levelz][orient] == 1) {
- norm = dwt_norm[levelx][levely][levelz][orient];
- /*fprintf(stdout,"[INFO] Level: %d %d %d Orient %d Dwt_norm: %f \n",level[0],level[1],level[2],orient,norm);*/
- } else {
- opj_wtfilt_t *wtfiltx =(opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
- opj_wtfilt_t *wtfilty =(opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
- opj_wtfilt_t *wtfiltz =(opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
- /*Fetch equivalent filters for each dimension*/
- dwt_getwtfilters(wtfiltx, dwtid[0]);
- dwt_getwtfilters(wtfilty, dwtid[1]);
- dwt_getwtfilters(wtfiltz, dwtid[2]);
- /*Calculate the corresponding norm */
- norm = dwt_calc_wtnorms(orient, level, dwtid, wtfiltx, wtfilty, wtfiltz);
- /*Save norm in array (no recalculation)*/
- dwt_norm[levelx][levely][levelz][orient] = norm;
- flagnorm[levelx][levely][levelz][orient] = 1;
- /*Free reserved space*/
- opj_free(wtfiltx->LPS); opj_free(wtfilty->LPS); opj_free(wtfiltz->LPS);
- opj_free(wtfiltx->HPS); opj_free(wtfilty->HPS); opj_free(wtfiltz->HPS);
- opj_free(wtfiltx); opj_free(wtfilty); opj_free(wtfiltz);
- /*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);*/
- }
- return norm;
+double dwt_getnorm(int orient, int level[3], int dwtid[3])
+{
+ int levelx = level[0];
+ int levely = level[1];
+ int levelz = (level[2] < 0) ? 0 : level[2];
+ double norm;
+
+ if (flagnorm[levelx][levely][levelz][orient] == 1) {
+ norm = dwt_norm[levelx][levely][levelz][orient];
+ /*fprintf(stdout,"[INFO] Level: %d %d %d Orient %d Dwt_norm: %f \n",level[0],level[1],level[2],orient,norm);*/
+ } else {
+ opj_wtfilt_t *wtfiltx = (opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
+ opj_wtfilt_t *wtfilty = (opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
+ opj_wtfilt_t *wtfiltz = (opj_wtfilt_t *) opj_malloc(sizeof(opj_wtfilt_t));
+ /*Fetch equivalent filters for each dimension*/
+ dwt_getwtfilters(wtfiltx, dwtid[0]);
+ dwt_getwtfilters(wtfilty, dwtid[1]);
+ dwt_getwtfilters(wtfiltz, dwtid[2]);
+ /*Calculate the corresponding norm */
+ norm = dwt_calc_wtnorms(orient, level, dwtid, wtfiltx, wtfilty, wtfiltz);
+ /*Save norm in array (no recalculation)*/
+ dwt_norm[levelx][levely][levelz][orient] = norm;
+ flagnorm[levelx][levely][levelz][orient] = 1;
+ /*Free reserved space*/
+ opj_free(wtfiltx->LPS);
+ opj_free(wtfilty->LPS);
+ opj_free(wtfiltz->LPS);
+ opj_free(wtfiltx->HPS);
+ opj_free(wtfilty->HPS);
+ opj_free(wtfiltz->HPS);
+ opj_free(wtfiltx);
+ opj_free(wtfilty);
+ opj_free(wtfiltz);
+ /*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);*/
+ }
+ return norm;
}
-/* <summary> */
-/* Calculate explicit stepsizes for DWT. */
-/* </summary> */
-void dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, int prec) {
- int totnumbands, bandno, diff;
-
- assert(tccp->numresolution[0] >= tccp->numresolution[2]);
- diff = tccp->numresolution[0] - tccp->numresolution[2]; /*if RESx=RESy != RESz */
- totnumbands = (7 * tccp->numresolution[0] - 6) - 4 * diff; /* 3-D */
-
- for (bandno = 0; bandno < totnumbands; bandno++) {
- double stepsize;
- int resno, level[3], orient, gain;
-
- /* Bandno: 0 - LLL 1 - LHL
- 2 - HLL 3 - HHL
- 4 - LLH 5 - LHH
- 6 - HLH 7 - HHH */
-
- resno = (bandno == 0) ? 0 : ( (bandno <= 3 * diff) ? ((bandno - 1) / 3 + 1) : ((bandno + 4*diff - 1) / 7 + 1));
- orient = (bandno == 0) ? 0 : ( (bandno <= 3 * diff) ? ((bandno - 1) % 3 + 1) : ((bandno + 4*diff - 1) % 7 + 1));
- level[0] = tccp->numresolution[0] - 1 - resno;
- level[1] = tccp->numresolution[1] - 1 - resno;
- level[2] = tccp->numresolution[2] - 1 - resno;
-
- /* Gain: 0 - LLL 1 - LHL
- 1 - HLL 2 - HHL
- 1 - LLH 2 - LHH
- 2 - HLH 3 - HHH */
- gain = (tccp->reversible == 0) ? 0 : ( (orient == 0) ? 0 :
- ( ((orient == 1) || (orient == 2) || (orient == 4)) ? 1 :
- (((orient == 3) || (orient == 5) || (orient == 6)) ? 2 : 3)) );
-
- if (tccp->qntsty == J3D_CCP_QNTSTY_NOQNT) {
- stepsize = 1.0;
- } else {
- double norm = dwt_getnorm(orient,level,tccp->dwtid); /*Fetch norms if irreversible transform (by the moment only I9.7)*/
- stepsize = (1 << (gain + 1)) / norm;
- }
- /*fprintf(stdout,"[INFO] Bandno: %d Orient: %d Level: %d %d %d Stepsize: %f\n",bandno,orient,level[0],level[1],level[2],stepsize);*/
- dwt_encode_stepsize((int) floor(stepsize * 8192.0), prec + gain, &tccp->stepsizes[bandno]);
- }
+/* <summary> */
+/* Calculate explicit stepsizes for DWT. */
+/* </summary> */
+void dwt_calc_explicit_stepsizes(opj_tccp_t * tccp, int prec)
+{
+ int totnumbands, bandno, diff;
+
+ assert(tccp->numresolution[0] >= tccp->numresolution[2]);
+ diff = tccp->numresolution[0] -
+ tccp->numresolution[2]; /*if RESx=RESy != RESz */
+ totnumbands = (7 * tccp->numresolution[0] - 6) - 4 * diff; /* 3-D */
+
+ for (bandno = 0; bandno < totnumbands; bandno++) {
+ double stepsize;
+ int resno, level[3], orient, gain;
+
+ /* Bandno: 0 - LLL 1 - LHL
+ 2 - HLL 3 - HHL
+ 4 - LLH 5 - LHH
+ 6 - HLH 7 - HHH */
+
+ resno = (bandno == 0) ? 0 : ((bandno <= 3 * diff) ? ((bandno - 1) / 3 + 1) : ((
+ bandno + 4 * diff - 1) / 7 + 1));
+ orient = (bandno == 0) ? 0 : ((bandno <= 3 * diff) ? ((bandno - 1) % 3 + 1) : ((
+ bandno + 4 * diff - 1) % 7 + 1));
+ level[0] = tccp->numresolution[0] - 1 - resno;
+ level[1] = tccp->numresolution[1] - 1 - resno;
+ level[2] = tccp->numresolution[2] - 1 - resno;
+
+ /* Gain: 0 - LLL 1 - LHL
+ 1 - HLL 2 - HHL
+ 1 - LLH 2 - LHH
+ 2 - HLH 3 - HHH */
+ gain = (tccp->reversible == 0) ? 0 : ((orient == 0) ? 0 :
+ (((orient == 1) || (orient == 2) || (orient == 4)) ? 1 :
+ (((orient == 3) || (orient == 5) || (orient == 6)) ? 2 : 3)));
+
+ if (tccp->qntsty == J3D_CCP_QNTSTY_NOQNT) {
+ stepsize = 1.0;
+ } else {
+ double norm = dwt_getnorm(orient, level,
+ tccp->dwtid); /*Fetch norms if irreversible transform (by the moment only I9.7)*/
+ stepsize = (1 << (gain + 1)) / norm;
+ }
+ /*fprintf(stdout,"[INFO] Bandno: %d Orient: %d Level: %d %d %d Stepsize: %f\n",bandno,orient,level[0],level[1],level[2],stepsize);*/
+ dwt_encode_stepsize((int) floor(stepsize * 8192.0), prec + gain,
+ &tccp->stepsizes[bandno]);
+ }
}
generated by cgit v1.2.3 (git 2.46.0) at 2026-04-20 18:23:29 +0000