1 // This is mul/vil3d/algo/vil3d_gauss_reduce.hxx
2 #ifndef vil3d_gauss_reduce_hxx_
3 #define vil3d_gauss_reduce_hxx_
4 //:
5 // \file
6 // \brief Functions to smooth and sub-sample 3D images in one direction
7 //
8 //  These are not templated because
9 //  - Each type tends to need a slightly different implementation
10 //  - Let's not have too many templates.
11 // \author Tim Cootes
12 
13 #include "vil3d_gauss_reduce.h"
14 //
15 #include <vil/algo/vil_gauss_reduce.h>
16 
17 //: Smooth and subsample single plane src_im in i to produce dest_im
18 //  Applies 1-5-8-5-1 filter in i, then samples
19 //  every other pixel.  Fills [0,(ni+1)/2-1][0,nj-1][0,nk-1] elements of dest
20 //  Assumes dest_im has sufficient data allocated.
21 //
22 //  By applying three times we can obtain a full gaussian smoothed and
23 //  sub-sampled 3D image
24 template<class T>
vil3d_gauss_reduce_i(const T * src_im,unsigned src_ni,unsigned src_nj,unsigned src_nk,std::ptrdiff_t s_i_step,std::ptrdiff_t s_j_step,std::ptrdiff_t s_k_step,T * dest_im,std::ptrdiff_t d_i_step,std::ptrdiff_t d_j_step,std::ptrdiff_t d_k_step)25 void vil3d_gauss_reduce_i(const T* src_im,
26                           unsigned src_ni, unsigned src_nj, unsigned src_nk,
27                           std::ptrdiff_t s_i_step, std::ptrdiff_t s_j_step,
28                           std::ptrdiff_t s_k_step,
29                           T* dest_im,
30                           std::ptrdiff_t d_i_step, std::ptrdiff_t d_j_step,
31                           std::ptrdiff_t d_k_step)
32 {
33   for (unsigned k=0;k<src_nk;++k)
34   {
35     vil_gauss_reduce_1plane(src_im, src_ni,src_nj, s_i_step,s_j_step,
36                             dest_im,d_i_step, d_j_step);
37     dest_im += d_k_step;
38     src_im  += s_k_step;
39   }
40 }
41 
42 
43 //: Smooth and subsample src_im to produce dest_im
44 //  Applies filter in i,j and k directions, then samples every other pixel.
45 //  Resulting image is (ni+1)/2 x (nj+1)/2 x (nk+1)/2
46 template<class T>
vil3d_gauss_reduce(const vil3d_image_view<T> & src_im,vil3d_image_view<T> & dest_im,vil3d_image_view<T> & work_im1,vil3d_image_view<T> & work_im2)47 void vil3d_gauss_reduce(const vil3d_image_view<T>& src_im,
48                               vil3d_image_view<T>& dest_im,
49                               vil3d_image_view<T>& work_im1,
50                               vil3d_image_view<T>& work_im2)
51 {
52   unsigned ni = src_im.ni();
53   unsigned nj = src_im.nj();
54   unsigned nk = src_im.nk();
55   unsigned n_planes = src_im.nplanes();
56 
57   // Output image size
58   unsigned ni2 = (ni+1)/2;
59   unsigned nj2 = (nj+1)/2;
60   unsigned nk2 = (nk+1)/2;
61 
62   if (work_im1.ni()<ni2 || work_im1.nj()<nj || work_im1.nk()<nk)
63     work_im1.set_size(ni2, nj, nk, 1 ); // Only need one plane for this image (the larger of the work images)
64 
65   if (work_im2.ni()<ni2 || work_im2.nj()<nj2 || work_im2.nk()<nk || work_im2.nplanes()<n_planes)
66     work_im2.set_size(ni2, nj2, nk, n_planes);
67 
68 
69   // Smooth and subsample in i, result in work_im1
70   for (unsigned p=0;p<n_planes;++p)
71   {
72     vil3d_gauss_reduce_i(
73       src_im.origin_ptr()+p*src_im.planestep(), ni, nj, nk,
74       src_im.istep(), src_im.jstep(), src_im.kstep(),
75       work_im1.origin_ptr(), work_im1.istep(), work_im1.jstep(), work_im1.kstep());
76 
77   // Smooth and subsample in j (by implicitly transposing), result in work_im2
78     vil3d_gauss_reduce_i(
79       work_im1.origin_ptr() ,nj, ni2, nk,
80       work_im1.jstep(), work_im1.istep(), work_im1.kstep(),
81       work_im2.origin_ptr()+p*work_im2.planestep(),
82       work_im2.jstep() ,work_im2.istep(), work_im2.kstep());
83   }
84   // Can resize output now, in case it is the same as the input.
85   dest_im.set_size(ni2, nj2, nk2, n_planes);
86 
87   // Smooth and subsample in k (by implicitly transposing)
88   for (unsigned p=0; p<n_planes; ++p)
89     vil3d_gauss_reduce_i(
90       work_im2.origin_ptr()+p*work_im2.planestep(), nk, ni2, nj2,
91       work_im2.kstep(), work_im2.istep() ,work_im2.jstep(),
92       dest_im.origin_ptr()+p*dest_im.planestep(),
93       dest_im.kstep(), dest_im.istep(), dest_im.jstep());
94 }
95 
96 
97 //: Smooth and subsample src_im along i and j to produce dest_im
98 //  Applies filter in i,j directions, then samples every other pixel.
99 //  Resulting image is (ni+1)/2 x (nj+1)/2 x nk
100 template<class T>
vil3d_gauss_reduce_ij(const vil3d_image_view<T> & src_im,vil3d_image_view<T> & dest_im,vil3d_image_view<T> & work_im1)101 void vil3d_gauss_reduce_ij(const vil3d_image_view<T>& src_im,
102                                  vil3d_image_view<T>& dest_im,
103                                  vil3d_image_view<T>& work_im1)
104 {
105   unsigned ni = src_im.ni();
106   unsigned nj = src_im.nj();
107   unsigned nk = src_im.nk();
108   unsigned n_planes = src_im.nplanes();
109 
110   // Output image size
111   unsigned ni2 = (ni+1)/2;
112   unsigned nj2 = (nj+1)/2;
113   unsigned nk2 = nk;
114 
115   if (work_im1.ni()<ni2 || work_im1.nj()<nj || work_im1.nk()<nk)
116     work_im1.set_size(ni2,nj,nk);
117   dest_im.set_size(ni2,nj2,nk2,n_planes);
118 
119   for (unsigned p=0;p<n_planes;++p)
120   {
121     // Smooth and subsample in i, result in work_im1
122     vil3d_gauss_reduce_i(
123       src_im.origin_ptr()+p*src_im.planestep(),ni,nj,nk,
124       src_im.istep(),src_im.jstep(),src_im.kstep(),
125       work_im1.origin_ptr(),work_im1.istep(),work_im1.jstep(),work_im1.kstep());
126 
127     // Smooth and subsample in j (by implicitly transposing), result in dest_im
128     vil3d_gauss_reduce_i(
129       work_im1.origin_ptr(),nj,ni2,nk,
130       work_im1.jstep(),work_im1.istep(),work_im1.kstep(),
131       dest_im.origin_ptr()+p*dest_im.planestep(),
132       dest_im.jstep(),dest_im.istep(),dest_im.kstep());
133   }
134 }
135 
136 
137 //: Smooth and subsample src_im along i and k to produce dest_im
138 //  Applies filter in i,k directions, then samples every other pixel.
139 //  Resulting image is (ni+1)/2 x nj x (nk+1)/2
140 template<class T>
vil3d_gauss_reduce_ik(const vil3d_image_view<T> & src_im,vil3d_image_view<T> & dest_im,vil3d_image_view<T> & work_im1)141 void vil3d_gauss_reduce_ik(const vil3d_image_view<T>& src_im,
142                                  vil3d_image_view<T>& dest_im,
143                                  vil3d_image_view<T>& work_im1)
144 {
145   unsigned ni = src_im.ni();
146   unsigned nj = src_im.nj();
147   unsigned nk = src_im.nk();
148   unsigned n_planes = src_im.nplanes();
149 
150   // Output image size
151   unsigned ni2 = (ni+1)/2;
152   unsigned nj2 = nj;
153   unsigned nk2 = (nk+1)/2;
154 
155   if (work_im1.ni()<ni2 || work_im1.nj()<nj || work_im1.nk()<nk)
156     work_im1.set_size(ni2,nj,nk);
157   dest_im.set_size(ni2,nj2,nk2,n_planes);
158 
159   for (unsigned p=0;p<n_planes;++p)
160   {
161     // Smooth and subsample in i, result in work_im1
162     vil3d_gauss_reduce_i(
163       src_im.origin_ptr()+p*src_im.planestep(),ni,nj,nk,
164       src_im.istep(),src_im.jstep(),src_im.kstep(),
165       work_im1.origin_ptr(),work_im1.istep(),work_im1.jstep(),work_im1.kstep());
166 
167     // Smooth and subsample in k (by implicitly transposing), result in dest_im
168     vil3d_gauss_reduce_i(
169         work_im1.origin_ptr(),nk,ni2,nj,
170         work_im1.kstep(),work_im1.istep(),work_im1.jstep(),
171         dest_im.origin_ptr()+p*dest_im.planestep(),
172         dest_im.kstep(),dest_im.istep(),dest_im.jstep());
173   }
174 }
175 
176 
177 //: Smooth and subsample src_im along j and k to produce dest_im
178 //  Applies filter in j,k directions, then samples every other pixel.
179 //  Resulting image is ni x (nj+1)/2 x (nk+1)/2
180 template<class T>
vil3d_gauss_reduce_jk(const vil3d_image_view<T> & src_im,vil3d_image_view<T> & dest_im,vil3d_image_view<T> & work_im1)181 void vil3d_gauss_reduce_jk(const vil3d_image_view<T>& src_im,
182                                  vil3d_image_view<T>& dest_im,
183                                  vil3d_image_view<T>& work_im1)
184 {
185   unsigned ni = src_im.ni();
186   unsigned nj = src_im.nj();
187   unsigned nk = src_im.nk();
188   unsigned n_planes = src_im.nplanes();
189 
190   // Output image size
191   unsigned ni2 = ni;
192   unsigned nj2 = (nj+1)/2;
193   unsigned nk2 = (nk+1)/2;
194 
195   if (work_im1.ni()<ni || work_im1.nj()<nj2 || work_im1.nk()<nk)
196     work_im1.set_size(ni,nj2,nk);
197 
198   dest_im.set_size(ni2,nj2,nk2,n_planes);
199 
200   for (unsigned p=0;p<n_planes;++p)
201   {
202     // Smooth and subsample in j (by implicitly transposing), result in work_im1
203     vil3d_gauss_reduce_i(
204       src_im.origin_ptr()+p*src_im.planestep(),nj,ni,nk,
205       src_im.jstep(),src_im.istep(),src_im.kstep(),
206       work_im1.origin_ptr(),work_im1.jstep(),work_im1.istep(),work_im1.kstep());
207 
208     // Smooth and subsample in k (by implicitly transposing), result in dest_im
209     vil3d_gauss_reduce_i(
210       work_im1.origin_ptr(),nk,ni,nj2,
211       work_im1.kstep(),work_im1.istep(),work_im1.jstep(),
212       dest_im.origin_ptr()+p*dest_im.planestep(),
213       dest_im.kstep(),dest_im.istep(),dest_im.jstep());
214   }
215 }
216 
217 
218 #undef VIL3D_GAUSS_REDUCE_INSTANTIATE
219 #define VIL3D_GAUSS_REDUCE_INSTANTIATE(T) \
220 template void vil3d_gauss_reduce_i(const T* src_im,   \
221                                    unsigned src_ni, unsigned src_nj, unsigned src_nk, \
222                                    std::ptrdiff_t s_i_step, std::ptrdiff_t s_j_step,  \
223                                    std::ptrdiff_t s_k_step,  \
224                                    T* dest_im,  \
225                                    std::ptrdiff_t d_i_step,  \
226                                    std::ptrdiff_t d_j_step, std::ptrdiff_t d_k_step); \
227 template void vil3d_gauss_reduce(const vil3d_image_view<T >& src_im, \
228                                  vil3d_image_view<T >& dest_im,  \
229                                  vil3d_image_view<T >& work_im1, \
230                                  vil3d_image_view<T >& work_im2);  \
231 template void vil3d_gauss_reduce_ij(const vil3d_image_view<T >& src_im,  \
232                                     vil3d_image_view<T >& dest_im, \
233                                     vil3d_image_view<T >& work_im1); \
234 template void vil3d_gauss_reduce_ik(const vil3d_image_view<T >& src_im,  \
235                                     vil3d_image_view<T >& dest_im, \
236                                     vil3d_image_view<T >& work_im1); \
237 template void vil3d_gauss_reduce_jk(const vil3d_image_view<T >& src_im,  \
238                                     vil3d_image_view<T >& dest_im, \
239                                     vil3d_image_view<T >& work_im1)
240 
241 #endif // vil3d_gauss_reduce_hxx_
242