1 // This is brl/bbas/bil/algo/bil_harr_wavelet_transform.hxx
2 #ifndef bil_harr_wavelet_transform_hxx_
3 #define bil_harr_wavelet_transform_hxx_
4 //:
5 // \file
6
7 #include "bil_harr_wavelet_transform.h"
8
9 template< class T>
10 static inline
bil_harr_helper(const T * src,T * dest,std::ptrdiff_t s_step,std::ptrdiff_t d_step,unsigned int length)11 void bil_harr_helper(const T* src, T* dest,
12 std::ptrdiff_t s_step, std::ptrdiff_t d_step,
13 unsigned int length)
14 {
15 for (unsigned int i=0; i<length; ++i){
16 *dest = T(*src + *(src+s_step))/2;
17 *(dest+length*d_step) = *dest - *src;
18 src += 2*s_step;
19 dest += d_step;
20 }
21 }
22
23 //: Compute the Harr wavelet transform on each plane
24 template< class T >
bil_harr_wavelet_transform(const vil_image_view<T> & src,vil_image_view<T> & dest)25 void bil_harr_wavelet_transform(const vil_image_view<T>& src,
26 vil_image_view<T>& dest)
27 {
28 unsigned ni = src.ni();
29 unsigned nj = src.nj();
30 unsigned np = src.nplanes();
31 dest.set_size(ni,nj,np);
32
33 std::ptrdiff_t istepS=src.istep(),jstepS=src.jstep(),pstepS=src.planestep();
34 std::ptrdiff_t istepD=dest.istep(),jstepD=dest.jstep(),pstepD=dest.planestep();
35 const T* planeS = src.top_left_ptr();
36 T* planeD = dest.top_left_ptr();
37
38 T* temp = new T[ni>nj?ni:nj];
39 for (unsigned p=0;p<np;++p,planeS+=pstepS,planeD+=pstepD)
40 {
41 const T* rowS = planeS;
42 T* rowD = planeD;
43 for (unsigned j=0;j<nj;++j,rowS+=jstepS,rowD+=jstepD)
44 {
45 T* temp_ptr = temp;
46 const T* temp_src = rowS;
47 for (unsigned int i=0;i<ni;++i, ++temp_ptr, temp_src+=istepS)
48 *temp_ptr = *temp_src;
49 for (unsigned int cni = ni; cni%2 == 0; cni /= 2){
50 bil_harr_helper(temp,rowD,1,istepD,cni/2);
51 temp_ptr = temp;
52 temp_src = rowD;
53 for (unsigned int i=0;i<cni/2;++i, ++temp_ptr, temp_src+=istepS)
54 *temp_ptr = *temp_src;
55 }
56 }
57
58 T* colD = planeD;
59 for (unsigned i=0;i<ni;++i,colD+=istepD)
60 {
61 for (unsigned int cnj = nj; cnj%2 == 0; cnj /= 2){
62 T* temp_ptr = temp;
63 const T* temp_src = colD;
64 for (unsigned int j=0;j<cnj;++j, ++temp_ptr, temp_src+=jstepD)
65 *temp_ptr = *temp_src;
66 bil_harr_helper(temp,colD,1,jstepD,cnj/2);
67 }
68 }
69 }
70 delete [] temp;
71 }
72
73
74 template< class T>
75 static inline
bil_harr_inv_helper(const T * src,T * dest,std::ptrdiff_t s_step,std::ptrdiff_t d_step,unsigned int length)76 void bil_harr_inv_helper(const T* src, T* dest,
77 std::ptrdiff_t s_step, std::ptrdiff_t d_step,
78 unsigned int length)
79 {
80 for (unsigned int i=0; i<length; ++i){
81 const T& diff = *(src+length*s_step);
82 *dest = *src - diff;
83 *(dest+d_step) = *src + diff;
84 src += s_step;
85 dest += 2*d_step;
86 }
87 }
88
89
90 //: Computes the inverse of the Harr wavelet transform on each plane
91 template< class T >
bil_harr_wavelet_inverse(const vil_image_view<T> & src,vil_image_view<T> & dest)92 void bil_harr_wavelet_inverse(const vil_image_view<T>& src,
93 vil_image_view<T>& dest)
94 {
95 unsigned ni = src.ni();
96 unsigned nj = src.nj();
97 unsigned np = src.nplanes();
98 dest.set_size(ni,nj,np);
99 if (src != dest)
100 dest.deep_copy(src);
101
102 //std::ptrdiff_t istepS=src.istep(),jstepS=src.jstep(),pstepS=src.planestep();
103 std::ptrdiff_t istepD=dest.istep(),jstepD=dest.jstep(),pstepD=dest.planestep();
104 //const T* planeS = src.top_left_ptr();
105 T* planeD = dest.top_left_ptr();
106
107 unsigned int min_ni = ni, min_nj = nj;
108 while (min_ni%2==0) min_ni /= 2;
109 while (min_nj%2==0) min_nj /= 2;
110 T* temp = new T[ni>nj?ni:nj];
111 for (unsigned p=0;p<np;++p,planeD+=pstepD)
112 {
113 //const T* rowS = planeS;
114 T* rowD = planeD;
115 for (unsigned j=0;j<nj;++j,rowD+=jstepD)
116 {
117 for (unsigned int cni = min_ni; cni < ni; cni *= 2){
118 T* temp_ptr = temp;
119 const T* temp_src = rowD;
120 for (unsigned int i=0;i<2*cni;++i, ++temp_ptr, temp_src+=istepD)
121 *temp_ptr = *temp_src;
122
123 bil_harr_inv_helper(temp,rowD,1,istepD,cni);
124 }
125 }
126
127 //const T* colS = planeD;
128 T* colD = planeD;
129 for (unsigned i=0;i<ni;++i,colD+=istepD)
130 {
131 for (unsigned int cnj = min_nj; cnj < nj; cnj *= 2){
132 T* temp_ptr = temp;
133 const T* temp_src = colD;
134 for (unsigned int j=0;j<2*cnj;++j, ++temp_ptr, temp_src+=jstepD)
135 *temp_ptr = *temp_src;
136 bil_harr_inv_helper(temp,colD,1,jstepD,cnj);
137 }
138 }
139 }
140 delete [] temp;
141 }
142
143
144 #define BIL_HARR_WAVELET_TRANSFORM_INSTANTIATE(T) \
145 template void \
146 bil_harr_wavelet_transform< T >(const vil_image_view< T >& src,\
147 vil_image_view< T >& dest); \
148 template void \
149 bil_harr_wavelet_inverse< T >(const vil_image_view< T >& src,\
150 vil_image_view< T >& dest)
151
152 #endif // bil_harr_wavelet_transform_hxx_
153