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