1 // This is mul/vimt/vimt_dog_pyramid_builder_2d.hxx
2 #ifndef vimt_dog_pyramid_builder_2d_hxx_
3 #define vimt_dog_pyramid_builder_2d_hxx_
4 //:
5 // \file
6 // \brief Build difference of gaussian pyramids of vimt_image_2d_of<T>
7 // \author Tim Cootes
8
9 #include <iostream>
10 #include <cstdlib>
11 #include <string>
12 #include "vimt_dog_pyramid_builder_2d.h"
13 #include "vimt_image_pyramid.h"
14
15 #ifdef _MSC_VER
16 # include <vcl_msvc_warnings.h>
17 #endif
18
19 #include <cassert>
20 #include <vnl/vnl_math.h> // for sqrt2
21 #include <vgl/vgl_point_2d.h>
22 #include <vgl/vgl_vector_2d.h>
23 #include <vil/algo/vil_gauss_filter.h>
24 #include <vil/vil_resample_bilin.h>
25 #include <vil/vil_math.h>
26
27 //=======================================================================
28
29 template <class T>
vimt_dog_pyramid_builder_2d()30 vimt_dog_pyramid_builder_2d<T>::vimt_dog_pyramid_builder_2d()
31
32 {
33 set_min_size(5, 5);
34 }
35
36 //=======================================================================
37
38 template<class T>
39 vimt_dog_pyramid_builder_2d<T>::~vimt_dog_pyramid_builder_2d() = default;
40
41 //=======================================================================
42 //: Define maximum number of levels to build
43 // Limits levels built in subsequent calls to build()
44 template<class T>
set_max_levels(int max_l)45 void vimt_dog_pyramid_builder_2d<T>::set_max_levels(int max_l)
46 {
47 if (max_l<1)
48 {
49 std::cerr<<"vimt_dog_pyramid_builder_2d<T>::set_max_levels() argument must be >=1\n";
50 std::abort();
51 }
52 max_levels_ = max_l;
53 }
54
55 //: Get the current maximum number levels allowed
56 template<class T>
max_levels() const57 int vimt_dog_pyramid_builder_2d<T>::max_levels() const
58 {
59 return max_levels_;
60 }
61
62 //=======================================================================
63 //: Create new (empty) pyramid on heap
64 // Caller responsible for its deletion
65 template<class T>
new_image_pyramid() const66 vimt_image_pyramid* vimt_dog_pyramid_builder_2d<T>::new_image_pyramid() const
67 {
68 return new vimt_image_pyramid;
69 }
70
71 //=======================================================================
72 //: Scale step between levels
73 template<class T>
scale_step() const74 double vimt_dog_pyramid_builder_2d<T>::scale_step() const
75 {
76 return 1.5;
77 }
78
79
80 //=======================================================================
81 //: Deletes all data in im_pyr
82 template<class T>
empty_pyr(vimt_image_pyramid & im_pyr) const83 void vimt_dog_pyramid_builder_2d<T>::empty_pyr(vimt_image_pyramid& im_pyr) const
84 {
85 for (int i=0; i<im_pyr.n_levels();++i)
86 delete im_pyr.data()[i];
87 }
88
89 //=======================================================================
90 //: Checks pyramid has at least n levels
91 template<class T>
check_pyr(vimt_image_pyramid & im_pyr,int n_levels) const92 void vimt_dog_pyramid_builder_2d<T>::check_pyr(vimt_image_pyramid& im_pyr, int n_levels) const
93 {
94 const int got_levels = im_pyr.n_levels();
95 if (got_levels >= n_levels && im_pyr(0).is_class(work_im_.is_a()))
96 {
97 if (im_pyr.n_levels()==n_levels) return;
98 else
99 {
100 for (int i=n_levels;i<got_levels;++i)
101 delete im_pyr.data()[i];
102 }
103 im_pyr.data().resize(n_levels);
104 return;
105 }
106
107 im_pyr.data().resize(n_levels);
108 empty_pyr(im_pyr);
109
110 for (int i=0;i<n_levels;++i)
111 im_pyr.data()[i] = new vimt_image_2d_of<T>;
112 }
113
114 //: Build pyramid
115 template<class T>
build(vimt_image_pyramid & dog_pyr,const vimt_image & im) const116 void vimt_dog_pyramid_builder_2d<T>::build(vimt_image_pyramid& dog_pyr,
117 const vimt_image& im) const
118 {
119 vimt_image_pyramid smooth_pyr;
120 build_dog(dog_pyr,smooth_pyr,im);
121 }
122
123 //=======================================================================
124 //: Build pyramid
125 template<class T>
build_dog(vimt_image_pyramid & dog_pyr,vimt_image_pyramid & smooth_pyr,const vimt_image & im,bool abs_diff) const126 void vimt_dog_pyramid_builder_2d<T>::build_dog(vimt_image_pyramid& dog_pyr,
127 vimt_image_pyramid& smooth_pyr,
128 const vimt_image& im, bool abs_diff) const
129 {
130 // Require image vimt_image_2d_of<T>
131 assert(im.is_class(work_im_.is_a()));
132
133 const vimt_image_2d_of<T>& base_image = static_cast<const vimt_image_2d_of<T>&>(im);
134
135 unsigned ni = base_image.image().ni();
136 unsigned nj = base_image.image().nj();
137
138 // Compute number of levels to pyramid so that top is no less
139 // than min_x_size_ x min_y_size_
140 int max_levels = 1;
141 while (ni>min_x_size_ && nj>min_y_size_)
142 {
143 max_levels++;
144 ni = 2*ni/3;
145 nj = 2*nj/3;
146 }
147
148 if (max_levels>max_levels_)
149 max_levels=max_levels_;
150
151 // Set up image pyramid
152 check_pyr(dog_pyr,max_levels);
153 check_pyr(smooth_pyr,max_levels);
154
155 vimt_image_2d_of<T>& smooth0 = static_cast<vimt_image_2d_of<T>&>( smooth_pyr(0));
156 vimt_image_2d_of<T>& dog0 = static_cast<vimt_image_2d_of<T>&>( dog_pyr(0));
157
158 vil_gauss_filter_5tap_params smooth_params(0.75);
159
160 vil_gauss_filter_5tap(base_image.image(),smooth0.image(),smooth_params,
161 dog0.image()); // Workspace
162
163 if (abs_diff)
164 vil_math_image_abs_difference(base_image.image(),smooth0.image(),dog0.image());
165 else
166 vil_math_image_difference(base_image.image(),smooth0.image(),dog0.image());
167
168 smooth0.set_world2im(base_image.world2im());
169 dog0.set_world2im(base_image.world2im());
170
171 unsigned n_planes = base_image.image().nplanes();
172
173 // Sort out world to image transformation for destination image
174 vimt_transform_2d scaling_trans;
175 scaling_trans.set_zoom_only(2.0/3.0,0,0);
176
177 // Workspace
178 vil_image_view<T> sub_sampled_image;
179
180
181 // Subsequent levels
182 for (int i=1;i<max_levels;++i)
183 {
184 vimt_image_2d_of<T>& smooth0 = static_cast<vimt_image_2d_of<T>&>( smooth_pyr(i-1));
185 vimt_image_2d_of<T>& smooth1 = static_cast<vimt_image_2d_of<T>&>( smooth_pyr(i));
186 vimt_image_2d_of<T>& dog1 = static_cast<vimt_image_2d_of<T>&>( dog_pyr(i));
187
188 // Subsample by a factor of 2/3
189 // Note - this could be implemented more efficiently
190 // since bilinear is sampling at pixel positions
191 // and on edges.
192 unsigned ni = smooth0.image().ni();
193 unsigned nj = smooth0.image().nj();
194 ni = 2*ni/3;
195 nj = 2*nj/3;
196 sub_sampled_image.set_size(ni,nj,n_planes);
197
198 vil_resample_bilin(smooth0.image(),sub_sampled_image,
199 0.0,0.0, 1.5,0.0, 0.0,1.5, ni,nj);
200
201 vil_gauss_filter_5tap(sub_sampled_image,smooth1.image(),
202 smooth_params,
203 dog1.image()); // Workspace
204 if (abs_diff)
205 vil_math_image_abs_difference(sub_sampled_image,smooth1.image(),
206 dog1.image());
207 else
208 vil_math_image_difference(sub_sampled_image,smooth1.image(),
209 dog1.image());
210
211 smooth1.set_world2im(scaling_trans*smooth0.world2im());
212 dog1.set_world2im(smooth1.world2im());
213 }
214
215 // Estimate width of pixels in base image
216 vgl_point_2d<double> c0(0.5*(ni-1),0.5*(nj-1));
217 vgl_point_2d<double> c1 = c0 + vgl_vector_2d<double> (1,1);
218 vimt_transform_2d im2world = base_image.world2im().inverse();
219 vgl_vector_2d<double> dw = im2world(c1) - im2world(c0);
220
221 double base_pixel_width = dw.length()/vnl_math::sqrt2;
222 double scale_step = 1.5;
223
224 dog_pyr.set_widths(base_pixel_width,scale_step);
225 smooth_pyr.set_widths(base_pixel_width,scale_step);
226 }
227
228 //=======================================================================
229 //: Extend pyramid (not implemented)
230 template<class T>
extend(vimt_image_pyramid &) const231 void vimt_dog_pyramid_builder_2d<T>::extend(vimt_image_pyramid& ) const
232 {
233 std::cerr << "vimt_dog_pyramid_builder_2d<T>::extend(vimt_image_pyramid&) NYI\n";
234 }
235
236
237 //=======================================================================
238
239 template<class T>
is_class(std::string const & s) const240 bool vimt_dog_pyramid_builder_2d<T>::is_class(std::string const& s) const
241 {
242 return s==vimt_dog_pyramid_builder_2d<T>::is_a() || vimt_image_pyramid_builder::is_class(s);
243 }
244
245 //=======================================================================
246
247 template<class T>
version_no() const248 short vimt_dog_pyramid_builder_2d<T>::version_no() const
249 {
250 return 1;
251 }
252
253 //=======================================================================
254
255 template<class T>
clone() const256 vimt_image_pyramid_builder* vimt_dog_pyramid_builder_2d<T>::clone() const
257 {
258 return new vimt_dog_pyramid_builder_2d<T>(*this);
259 }
260
261 //=======================================================================
262
263 template<class T>
print_summary(std::ostream &) const264 void vimt_dog_pyramid_builder_2d<T>::print_summary(std::ostream&) const
265 {
266 }
267
268 //=======================================================================
269
270 template<class T>
b_write(vsl_b_ostream & bfs) const271 void vimt_dog_pyramid_builder_2d<T>::b_write(vsl_b_ostream& bfs) const
272 {
273 vsl_b_write(bfs,version_no());
274 vsl_b_write(bfs,max_levels_);
275 vsl_b_write(bfs,min_x_size_);
276 vsl_b_write(bfs,min_y_size_);
277 }
278
279 //=======================================================================
280
281 template<class T>
b_read(vsl_b_istream & bfs)282 void vimt_dog_pyramid_builder_2d<T>::b_read(vsl_b_istream& bfs)
283 {
284 if (!bfs) return;
285
286 short version;
287 vsl_b_read(bfs,version);
288 switch (version)
289 {
290 // version number starts at 2 to follow on from the old mil stuff
291 case (1):
292 vsl_b_read(bfs,max_levels_);
293 vsl_b_read(bfs,min_x_size_);
294 vsl_b_read(bfs,min_y_size_);
295 break;
296 default:
297 std::cerr << "I/O ERROR: vimt_dog_pyramid_builder_2d<T>::b_read(vsl_b_istream&)\n"
298 << " Unknown version number "<< version << '\n';
299 bfs.is().clear(std::ios::badbit); // Set an unrecoverable IO error on stream
300 return;
301 }
302 }
303
304 template<class T>
gauss_reduce(const vimt_image_2d_of<T> &,vimt_image_2d_of<T> &) const305 void vimt_dog_pyramid_builder_2d<T>::gauss_reduce(const vimt_image_2d_of<T>& /*src_im*/,
306 vimt_image_2d_of<T>& /*dest_im*/) const
307 {
308 std::cerr << "ERROR: vimt_dog_pyramid_builder_2d<T>::gauss_reduce() not yet implemented\n";
309 }
310
311 #define VIMT_DOG_PYRAMID_BUILDER_2D_INSTANTIATE(T) \
312 template <> std::string vimt_dog_pyramid_builder_2d<T >::is_a() const \
313 { return std::string("vimt_dog_pyramid_builder_2d<" #T ">"); }\
314 template class vimt_dog_pyramid_builder_2d<T >
315
316 #endif // vimt_dog_pyramid_builder_2d_hxx_
317