Subversion Repositories Kolibri OS

Rev

Go to most recent revision | Blame | Last modification | View Log | RSS feed

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