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