Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Newer
Older
100644 332 lines (268 sloc) 14.079 kB
6678052 @zmoratto Still working the repo organization
zmoratto authored
1 // __BEGIN_LICENSE__
e958375 @zmoratto Run relicense
zmoratto authored
2 // Copyright (C) 2006-2011 United States Government as represented by
19bd812 @zmoratto Update licensing
zmoratto authored
3 // the Administrator of the National Aeronautics and Space Administration.
4 // All Rights Reserved.
6678052 @zmoratto Still working the repo organization
zmoratto authored
5 // __END_LICENSE__
6
19bd812 @zmoratto Update licensing
zmoratto authored
7
6678052 @zmoratto Still working the repo organization
zmoratto authored
8 /// \file SurfaceNURBS.h
9 ///
10
11 #ifndef __SURFACENURBS_H__
12 #define __SURFACENURBS_H__
13
14 #include <asp_config.h>
15
16 #include <vw/Image/ImageView.h>
17 #include <vw/Image/PixelTypes.h>
18 #include <vw/Image/Interpolation.h>
19 #include <vw/Stereo/DisparityMap.h>
20 #include <vw/Core/ProgressCallback.h>
21 #include <vw/Math/Functions.h>
22
23 #if defined(ASP_HAVE_PKG_MBA) && ASP_HAVE_PKG_MBA==1
24 #include <MBA.h>
25 #include <UCButils.h>
26 #include <boost/shared_ptr.hpp>
27
28 // These type computation routines help us to decide which pixels are
29 // considered "valid" when computing the NURBS fit.
30 template <class PixelT> struct is_valid_pixel {
31 static bool value(PixelT const& pix) { return true; }
32 };
33
34
35 template <class ChannelT> struct is_valid_pixel<vw::PixelDisparity<ChannelT> > {
36 static bool value(vw::PixelDisparity<ChannelT> const& pix) {
37 return !pix.missing();
38 }
39 };
40
41 template <class ImageT>
42 static void init_mba_xy(vw::ImageViewBase<ImageT> const& image,
43 std::vector<double> &x_array,
44 std::vector<double> &y_array) {
45
46 int i = 0;
47 for (int32 row = 0; row < image.impl().rows(); row++) {
48 for (int32 col = 0; col < image.impl().cols(); col++) {
49 if (is_valid_pixel<typename ImageT::pixel_type>::value(image.impl()(col, row))) {
50 x_array.push_back(col);
51 y_array.push_back(row);
52 i++;
53 }
54 }
55 }
56
57 if ((image.impl().cols() * image.impl().rows() == 0) || (x_array.size() == 0))
58 throw vw::ArgumentErr() << "MBA NURBS Adaptation failed. Image does not contain any points or image has zero valid pixels.";
59 }
60
61
62 // Copy the data into STL vectors; this is the format expected by
63 // the MBA NURBS code. The View must have a PixelGray pixel type.
64 template <class ViewT>
65 static void init_mba_z(vw::ImageViewBase<ViewT> const& image,
66 std::vector<double> const& x_array,
67 std::vector<double> const& y_array,
68 std::vector<double> &z_array, int channel) {
69 z_array.resize(x_array.size());
70 for (unsigned i = 0; i < x_array.size(); ++i)
71 z_array[i] = vw::compound_select_channel<typename vw::CompoundChannelType<typename ViewT::pixel_type>::type>(image.impl()(int(x_array[i]),int(y_array[i])), channel);
72 }
73
74 template <class PixelT>
75 class SurfaceEvaluator {
76
77 std::vector<UCBspl::SplineSurface> m_fitted_surfaces;
78
79 public:
80 SurfaceEvaluator(std::vector<UCBspl::SplineSurface> const& fitted_surfaces) :
81 m_fitted_surfaces(fitted_surfaces) {}
82
83 PixelT operator()(int i, int j) const {
84 PixelT result;
85 if (i > m_fitted_surfaces[0].umin() && i < m_fitted_surfaces[0].umax() &&
86 j > m_fitted_surfaces[0].vmin() && j < m_fitted_surfaces[0].vmax()) {
87 for (int c = 0; c < vw::PixelNumChannels<PixelT>::value; ++c)
88 result[c] = m_fitted_surfaces[c].f((float)i,(float)j);
89 }
90 return result;
91 }
92 };
93
94 template <class ChannelT>
95 class SurfaceEvaluator<vw::PixelDisparity<ChannelT> > {
96
97 const std::vector<UCBspl::SplineSurface> &m_fitted_surfaces;
98
99 public:
100 SurfaceEvaluator(std::vector<UCBspl::SplineSurface> const& fitted_surfaces) :
101 m_fitted_surfaces(fitted_surfaces) {
102 VW_ASSERT(m_fitted_surfaces.size() == 2,
103 vw::ArgumentErr() << "SurfaceEvaluator: expecting two surfaces for pixels of PixelDisparty type.");
104 }
105
106 vw::PixelDisparity<ChannelT> operator()(int i, int j) const {
107 vw::PixelDisparity<ChannelT> result;
108 if (i > m_fitted_surfaces[0].umin() && i < m_fitted_surfaces[0].umax() &&
109 j > m_fitted_surfaces[0].vmin() && j < m_fitted_surfaces[0].vmax()) {
110 return vw::PixelDisparity<ChannelT> (m_fitted_surfaces[0].f((float)i,(float)j),
111 m_fitted_surfaces[1].f((float)i,(float)j));
112 } else {
113 return vw::PixelDisparity<ChannelT>();
114 }
115 }
116 };
117
118 // These two classes are just a quick fix for now since the
119 // PixelDisparity pixel type really only has two channels that we want
120 // to fit NURBS data on, but three channels total as far as the VW
121 // pixel type system is concerned.
122 template <class PixelT>
123 struct NURBSChannelCount { static int value() { return vw::PixelNumChannels<PixelT>::value; } };
124 template <class ChannelT>
125 struct NURBSChannelCount<vw::PixelDisparity<ChannelT> > { static int value() { return 2; } };
126
127 /// Fit a 2D spline surface to a given image of data. Each channel of
128 /// the image is handled as an independent 2D surface.
129 template <class ViewT>
130 vw::ImageView<typename ViewT::pixel_type> MBASurfaceNURBS(vw::ImageViewBase<ViewT> const& input,
131 int numIterations) {
132
133 typedef typename ViewT::pixel_type pixel_type;
134
135 // NOTE: we assume only one channel here!
136 VW_ASSERT(input.impl().planes() == 1, vw::NoImplErr() <<
137 "FindNURBSSurface() does not have support for images with more than one plane.");
138
139 boost::shared_ptr<std::vector<double> > x_arr(new std::vector<double>);
140 boost::shared_ptr<std::vector<double> > y_arr(new std::vector<double>);
141
142 // Set the size ratio of the spline space. This causes the
143 // algorithm to create more control points in one dimenion than in
144 // the other.
145 int m0 = 1, n0 = 1;
146 if (input.impl().cols() > input.impl().rows())
147 m0 = input.impl().cols() / input.impl().rows();
148 else
149 n0 = input.impl().rows() / input.impl().cols();
150
151 // Initialize the grid
152 init_mba_xy(input.impl(), *x_arr, *y_arr);
153 // Create a vector of spline surfaces
154 std::vector<UCBspl::SplineSurface> fitted_surfaces(NURBSChannelCount<pixel_type>::value());
155
156 // For each channel in the image, fit a spline surface
157 for (int channel = 0; channel < NURBSChannelCount<pixel_type>::value(); ++channel) {
158 boost::shared_ptr<std::vector<double> > z_arr(new std::vector<double>);
159 init_mba_z(input.impl(), *x_arr, *y_arr, *z_arr, channel);
160 std::cout << "\t Fitting spline surface for channel " << channel << ".\n";
161 MBA mba_z(x_arr, y_arr, z_arr); // Initialize with scattered data
162 mba_z.MBAalg(m0,n0,numIterations); // Create spline surface
163 fitted_surfaces[channel] = mba_z.getSplineSurface(); // Get the spline surface object
164 std::cout << "\t Range -- U: [" << fitted_surfaces[channel].umin() << " " << fitted_surfaces[channel].umax() << "]\n";
165 std::cout << "\t V: [" << fitted_surfaces[channel].vmin() << " " << fitted_surfaces[channel].vmax() << "]\n";
166 }
167
168 // Free up some memory.
169 (*x_arr).clear();
170 (*y_arr).clear();
171
172 // Rasterize the surface by iterating over all of the x and y values.
173 vw::ImageView<pixel_type> output(input.impl().cols(), input.impl().rows());
174 std::cout << "\t Evaluating surfaces.\n";
175 SurfaceEvaluator<pixel_type> evaluator(fitted_surfaces);
176 for (int j = 0; j < input.impl().rows(); ++j) {
177 for (int i = 0; i < input.impl().cols(); ++i) {
178 output(i,j) = evaluator(i,j);
179 }
180 }
181 // // Sample surface and print to VRML file.
182 // UCBspl::printVRMLgrid("qwe.wrl", surface, 50, 50, true);
183
184 return output;
185 }
186
187
188
189
190 // Rapid resample()
191 //
192 // Downsample a disparity map by averaging pixels. In this version of
193 // the algorithm, missing pixels are "consumed" by valid pixels.
194 template <class ViewT>
195 ImageView<typename ViewT::pixel_type > disparity_rapid_downsample(ImageViewBase<ViewT> const& view_,
196 int downsample_size) {
197
198 EdgeExtensionView<ViewT, ConstantEdgeExtension> view(view_.impl(), ConstantEdgeExtension());
199 typedef typename ViewT::pixel_type pixel_type;
200 typedef Vector2 accumulator_type;
201 typedef typename EdgeExtensionView<ViewT, ConstantEdgeExtension>::pixel_accessor pixel_accessor;
202
203 int32 resized_rows = view.impl().rows() / downsample_size;
204 int32 resized_cols = view.impl().cols() / downsample_size;
205 ImageView<pixel_type> result(resized_cols, resized_rows);
206
207 std::cout << "\t Subsampling disparity map to output size: " << resized_cols << " x " << resized_rows << "\n";;
208
209 TerminalProgressCallback progress_callback(ErrorMessage, "\t Subsampling: ");
210 progress_callback.report_progress(0);
211
212 for (int32 j = 0; j < view.impl().rows(); j+=downsample_size) {
213 progress_callback.report_progress((float)j / view.impl().rows());
214 for (int32 i = 0; i < view.impl().cols(); i+= downsample_size) {
215 accumulator_type sum = accumulator_type();
216 int num_good_pixels = 0;
217
218 pixel_accessor row_acc = view.origin().advance(i,j);
219 for (int32 v = 0; v < downsample_size; ++v) {
220 pixel_accessor col_acc = row_acc;
221 for (int32 u = 0; u < downsample_size; ++u) {
222 if ( !((*col_acc).missing()) ) {
223 sum[0] += (*col_acc).h();
224 sum[1] += (*col_acc).v();
225 ++num_good_pixels;
226 }
227 col_acc.next_col();
228 }
229 row_acc.next_row();
230 }
231
232 // If there were any valid pixels, we count them towards the
233 // average. If there are none, we mark this location as a
234 // missing pixel.
235 pixel_type avg;
236 if (num_good_pixels != 0) {
237 avg = pixel_type(sum[0] / num_good_pixels,
238 sum[1] / num_good_pixels);
239 } else {
240 avg = pixel_type();
241 }
242
243 // This bounds checking keeps us from straying outside of the
244 // result image.
245 if (i/downsample_size < resized_cols && j/downsample_size < resized_rows)
246 result(i/downsample_size, j/downsample_size) = avg;
247 }
248 }
249 progress_callback.report_finished();
250 return result;
251 }
252
253
254
255 class HoleFillView : public ImageViewBase<HoleFillView> {
256 ImageViewRef<PixelDisparity<float> > m_original_disparity_map;
257 ImageView<PixelDisparity<float> > m_lowres_disparity_map;
258 int m_subsample_factor;
259
260 public:
261 typedef PixelDisparity<float> pixel_type;
262 typedef const pixel_type result_type;
263 typedef ProceduralPixelAccessor<HoleFillView> pixel_accessor;
264
265 template <class DisparityViewT>
266 HoleFillView(DisparityViewT disparity_map, int subsample_factor = 16)
267 : m_original_disparity_map(disparity_map), m_subsample_factor(subsample_factor) {
268
269 // Compute the bounding box that encompasses all of the
270 // available points.
271 m_lowres_disparity_map = MBASurfaceNURBS(disparity_rapid_downsample(disparity_map, m_subsample_factor), 10);
272 }
273
274 inline int32 cols() const { return m_original_disparity_map.cols(); }
275 inline int32 rows() const { return m_original_disparity_map.rows(); }
276 inline int32 planes() const { return 1; }
277
278 inline pixel_accessor origin() const { return pixel_accessor(*this); }
279
280 inline result_type operator()( int i, int j, int p=0 ) const {
281
282 // If we can, we take the pixels from the original disparity
283 // image.
284 if ( !(m_original_disparity_map(i,j).missing()) ) {
285 return m_original_disparity_map(i,j);
286
287 // Otherwise, we resort to using the hole fill values. These
288 // values are the result of a low resolution NURBS fit followed by
289 // bicubic interpolation.
290 } else {
291 float ii = float(i)/m_subsample_factor;
292 float jj = float(j)/m_subsample_factor;
293
294 int32 x = vw::math::impl::_floor(ii), y = vw::math::impl::_floor(jj);
295 double normx = ii-x, normy = jj-y;
296
297 double s0 = ((2-normx)*normx-1)*normx; double t0 = ((2-normy)*normy-1)*normy;
298 double s1 = (3*normx-5)*normx*normx+2; double t1 = (3*normy-5)*normy*normy+2;
299 double s2 = ((4-3*normx)*normx+1)*normx; double t2 = ((4-3*normy)*normy+1)*normy;
300 double s3 = (normx-1)*normx*normx; double t3 = (normy-1)*normy*normy;
301
302 EdgeExtensionView<ImageView<PixelDisparity<float> >,ConstantEdgeExtension> extended_view = edge_extend(m_lowres_disparity_map,ConstantEdgeExtension());
303 if (extended_view(x-1,y-1).missing() || extended_view(x-1,y).missing() || extended_view(x-1,y+1).missing() ||
304 extended_view(x ,y-1).missing() || extended_view(x ,y).missing() || extended_view(x ,y+1).missing() ||
305 extended_view(x+1,y-1).missing() || extended_view(x+1,y).missing() || extended_view(x+1,y+1).missing() ) {
306 return PixelDisparity<float>();
307 } else {
308 float hdisp = ( ( s0*extended_view(x-1,y-1,p).h() + s1*extended_view(x+0,y-1,p).h() + s2*extended_view(x+1,y-1,p).h() + s3*extended_view(x+2,y-1,p).h() ) * t0 +
309 ( s0*extended_view(x-1,y+0,p).h() + s1*extended_view(x+0,y+0,p).h() + s2*extended_view(x+1,y+0,p).h() + s3*extended_view(x+2,y+0,p).h() ) * t1 +
310 ( s0*extended_view(x-1,y+1,p).h() + s1*extended_view(x+0,y+1,p).h() + s2*extended_view(x+1,y+1,p).h() + s3*extended_view(x+2,y+1,p).h() ) * t2 +
311 ( s0*extended_view(x-1,y+2,p).h() + s1*extended_view(x+0,y+2,p).h() + s2*extended_view(x+1,y+2,p).h() + s3*extended_view(x+2,y+2,p).h() ) * t3 ) * 0.25;
312 float vdisp = ( ( s0*extended_view(x-1,y-1,p).v() + s1*extended_view(x+0,y-1,p).v() + s2*extended_view(x+1,y-1,p).v() + s3*extended_view(x+2,y-1,p).v() ) * t0 +
313 ( s0*extended_view(x-1,y+0,p).v() + s1*extended_view(x+0,y+0,p).v() + s2*extended_view(x+1,y+0,p).v() + s3*extended_view(x+2,y+0,p).v() ) * t1 +
314 ( s0*extended_view(x-1,y+1,p).v() + s1*extended_view(x+0,y+1,p).v() + s2*extended_view(x+1,y+1,p).v() + s3*extended_view(x+2,y+1,p).v() ) * t2 +
315 ( s0*extended_view(x-1,y+2,p).v() + s1*extended_view(x+0,y+2,p).v() + s2*extended_view(x+1,y+2,p).v() + s3*extended_view(x+2,y+2,p).v() ) * t3 ) * 0.25;
316 return PixelDisparity<float>(hdisp,vdisp);
317 }
318 }
319 }
320
321 /// \cond INTERNAL
322 typedef HoleFillView prerasterize_type;
323 inline prerasterize_type prerasterize( BBox2i const& bbox ) const { return *this; }
324 template <class DestT> inline void rasterize( DestT const& dest, BBox2i const& bbox ) const { vw::rasterize( prerasterize(bbox), dest, bbox ); }
325 /// \endcond
326
327 };
328
329 #endif // have MBA nurbs
330
331 #endif // __SURFACENURBS_H__
Something went wrong with that request. Please try again.