diff --git a/src/Makefile.am b/src/Makefile.am index 16c8f4d4c..45cbd0a6b 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -57,6 +57,8 @@ stereo_code = \ SurfaceNURBS.h \ TabulatedDataReader.cc \ TabulatedDataReader.h \ + Outliers.cc \ + Outliers.h \ $(rmax_code) \ $(bundle_code) diff --git a/src/Outliers.cc b/src/Outliers.cc new file mode 100644 index 000000000..1a96be480 --- /dev/null +++ b/src/Outliers.cc @@ -0,0 +1,282 @@ +#include +#include +#include + +using namespace vw; + + +// Pixel accessor function +class MaxFilter : public ReturnFixedType > { + int m_kern_size; + int m_threshold; + +public: + MaxFilter(int kern_size, int threshold) : m_kern_size(kern_size), m_threshold(threshold) {} + + BBox2i work_area() const { return BBox2i(-m_kern_size,-m_kern_size, + m_kern_size*2,m_kern_size*2); } + + template + typename boost::remove_reference::type + operator()(PixelAccessorT acc) const { + + float max_val = ScalarTypeLimits::lowest(); + + PixelAccessorT row_acc = acc; + row_acc.advance(-m_kern_size, -m_kern_size/2); + + for (int j = -m_kern_size/2; j < -m_kern_size/2+m_kern_size; ++j) { + PixelAccessorT col_acc = row_acc; + for (int i = -m_kern_size; i < -m_kern_size/2+m_kern_size; ++i) { + if (*col_acc > max_val && i*i+j*j <= pow(m_kern_size/2,2)) + max_val = *col_acc; + col_acc.next_col(); + } + row_acc.next_row(); + } + + if (max_val > m_threshold) + return max_val; + else + return *acc; + } +}; + + +// ------------------------------------------------------------------ + +class Histogram { + std::vector m_bins; + int m_nbins; + long int m_nvalid; + +public: + template + Histogram(ImageViewBase const& view, int nbins = 256) { + + TerminalProgressCallback callback(InfoMessage, "\t--> Building Histogram: "); + m_nbins = nbins; + m_nvalid = 0; + m_bins.resize(nbins); + for (unsigned i = 0; i < m_nbins; ++i) + m_bins[i] = 0.0; + + // Determine the standard range of pixel values. FIXME: this is + // hardwired for now for floating point images between [0, 1]. + float lo = 0.0; + float hi = 1.0; + + // Fill the bins using the pixel data in the image. + // + // FIXME: This is hardwired for PixelGray images for now! + for (unsigned j = 0; j < view.impl().rows(); ++j) { + callback.report_progress(float(j) / view.impl().rows()); + for (unsigned i = 0; i < view.impl().cols(); ++i) { + if ( is_valid( view.impl()(i,j) ) ) { + float val = view.impl()(i,j)[0]; + int idx = floor(val * (m_nbins-1)); + if (idx < 0) idx = 0; + if (idx > 255) idx = 255; + + m_bins[idx]++; + m_nvalid++; + } + } + } + + // Normalize the bins + for (unsigned i = 0; i < m_nbins; ++i) { + m_bins[i] /= (double)m_nvalid; + // std::cout << (float(i)/m_nbins) << " : " << m_bins[i] << "\n"; + } + callback.report_finished(); + + } + + float find_threshold(float percentile) { + float accum = 0; + for (unsigned i = 0; i < m_nbins; ++i) { + accum += m_bins[i]; + if (accum >= percentile) + return float(i) / m_nbins; + } + + // Should not be reached + return 1.0; + } + + + +}; + +// ------------------------------------------------------------------ + +/// DisparityTransform image transform functor +/// +/// Transform points in an image based on offset values in a disparity +/// map. +class DisparityTransform : public TransformBase { + ImageViewRef > m_offset_image; +public: + template + DisparityTransform(ImageViewBase const& offset_image) + : m_offset_image(offset_image.impl()) {} + + inline Vector2 reverse(const Vector2 &p) const { + int32 x = (int32) p.x(), y = (int32) p.y(); + // VW_DEBUG_ASSERT(x>=0 && y>=0 && x offset = m_offset_image(x,y); + if (offset.missing()) + return Vector2(-1, p.y()); // Missing pixels return a value + // outside of the bounds of the + // original image. Therefore, edge + // extension will determine what value + // missing pixels take on. + else + return Vector2(p.x() + offset.h(), p.y() + offset.v()); + } +}; + +// ------------------------------------------------------------------ + +void flood(ImageView > &in, int i, int j, int &num_pixels) { + + if (i >=0 && j >= 0 && i < in.cols() && j < in.rows() && in(i,j) == 1 ) { + num_pixels++; + in (i,j) = -1; // Mark as visited + + // Explore the 8 connected neighborhood + if (num_pixels < 500) { + flood(in, i+1, j-1, num_pixels); + flood(in, i , j-1, num_pixels); + flood(in, i-1, j-1, num_pixels); + flood(in, i+1, j , num_pixels); + flood(in, i-1, j , num_pixels); + flood(in, i+1, j+1, num_pixels); + flood(in, i , j+1, num_pixels); + flood(in, i-1, j+1, num_pixels); + } + } +} + +void fill(ImageView > const&in, int i, int j, int value) { + if (i >=0 && j >= 0 && i < in.cols() && j < in.rows() && in(i,j) == -1 ) { + in(i,j) = value; + + // Explore the 8 connected neighborhood + fill(in, i+1, j-1, value); + fill(in, i , j-1, value); + fill(in, i-1, j-1, value); + fill(in, i+1, j , value); + fill(in, i-1, j , value); + fill(in, i+1, j+1, value); + fill(in, i , j+1, value); + fill(in, i-1, j+1, value); + } +} + +// Region size determination (implemented as a recursive flood fill) +template +ImageView > region_size_image(ImageViewBase const& in) { + + ImageView > out = in.impl(); + + TerminalProgressCallback callback(InfoMessage, "\t--> Computing region sizes: "); + for (unsigned j=0; j < out.rows(); ++j) { + callback.report_progress(float(j) / out.rows()); + for (unsigned i=0; i < out.cols(); ++i) { + int num_pixels = 0; + flood(out, i, j, num_pixels); + fill(out, i, j, num_pixels); + } + } + callback.report_finished(); + + return out; +} + + + +// ------------------------------------------------------------------ + + +void photometric_outlier_rejection(std::string out_prefix, int kernel_size) { + + // DiskImageView > dust_mask(out_prefix + "-DustMask.tif"); + // ImageView > out = region_size_image(dust_mask); + // std::cout << out.cols() << " " << out.rows() << " " << out.planes() << " " << out.channels() << "\n"; + // write_image(out_prefix + "-DustMaskFixed.tif", + // out, + // TerminalProgressCallback(InfoMessage, "\t--> Writing region size image: ")); + + // ImageView > padded_dust_mask = per_pixel_accessor_filter(out, + // MaxFilter(kernel_size, 100)); + // ImageView > padded_dust_mask2 = per_pixel_accessor_filter(padded_dust_mask, + // MaxFilter(kernel_size, 10)); + + // write_image(out_prefix + "-PaddedDustMaskFixed.tif", + // padded_dust_mask2, + // TerminalProgressCallback(InfoMessage, "\t--> Writing padded dust mask: ")); + + // exit(0); + + + // Open the necessary files + DiskImageView > right_disk_image(out_prefix + "-R.tif"); + DiskImageView > disparity_disk_image(out_prefix + "-F.exr"); + + std::cout << "Dust-based outlier detection (warning: tested with Apollo imagery only!!)\n"; + DisparityTransform trans(disparity_disk_image); + write_image(out_prefix + "-R-reproj.tif", + transform(right_disk_image, trans, ZeroEdgeExtension()), + TerminalProgressCallback(InfoMessage, "\t--> Reprojecting Right Image: ")); + + DiskImageView > left_disk_image(out_prefix + "-L.tif"); + DiskImageView > right_reproj_image(out_prefix + "-R-reproj.tif"); + ImageViewRef > > left = create_mask(left_disk_image); + ImageViewRef > > right = create_mask(right_reproj_image, 0); + + ImageView > diff = abs(left_disk_image-right_reproj_image); + ImageView > > masked_diff = copy_mask(diff, right); + Histogram hist(masked_diff); + float thresh = hist.find_threshold(0.9999); + std::cout << "\t--> Using Threshold: " << thresh << "\n"; + + ImageViewRef > dust_mask = threshold(apply_mask(masked_diff,0),thresh,0.0,1.0); + + // Compute region sizes + ImageView > out = region_size_image(dust_mask); + write_image(out_prefix + "-RegionSize.tif", + out, + TerminalProgressCallback(InfoMessage, "\t--> Writing region size image: ")); + + // The first padding takes care of the large dust monsters, adding extra padding around them. + ImageView > padded_dust_mask = per_pixel_accessor_filter(out, + MaxFilter(kernel_size*1.5, 50)); + // The first padding takes care of the smaller dust spots. + ImageView > padded_dust_mask2 = 1.0 - clamp(per_pixel_accessor_filter(padded_dust_mask, + MaxFilter(kernel_size*1.5, 0))); + + DiskImageView > input_mask_image(out_prefix + "-NurbsMask.tif"); + ImageViewRef > final_mask = apply_mask(copy_mask(input_mask_image, + create_mask(padded_dust_mask2, 0))); + + + write_image(out_prefix + "-DustDiff.tif", + apply_mask(masked_diff, 0), + TerminalProgressCallback(InfoMessage, "\t--> Writing dust diff: ")); + + write_image(out_prefix + "-DustMask.tif", + dust_mask, + TerminalProgressCallback(InfoMessage, "\t--> Writing dust mask: ")); + + write_image(out_prefix + "-PaddedDustMask.tif", + padded_dust_mask2, + TerminalProgressCallback(InfoMessage, "\t--> Writing padded dust mask: ")); + + write_image(out_prefix + "-NurbsDustMask.tif", + final_mask, + TerminalProgressCallback(InfoMessage, "\t--> Writing dust NURBS mask: ")); + +} diff --git a/src/Outliers.h b/src/Outliers.h new file mode 100644 index 000000000..fbe68b587 --- /dev/null +++ b/src/Outliers.h @@ -0,0 +1,3 @@ +#include + +void photometric_outlier_rejection(std::string out_prefix, int kernel_size); diff --git a/src/stereo.cc b/src/stereo.cc index 08da0c451..1db06cd18 100644 --- a/src/stereo.cc +++ b/src/stereo.cc @@ -61,6 +61,7 @@ using namespace vw::cartography; #include "SurfaceNURBS.h" #include "BundleAdjustUtils.h" #include "MedianFilter.h" +#include "Outliers.h" // Support for Malin DDD image files #include "MRO/DiskImageResourceDDD.h" @@ -96,6 +97,7 @@ enum { PREPROCESSING = 0, CORRELATION, REFINEMENT, FILTERING, + DUST_REJECTION, POINT_CLOUD, WIRE_MESH, NUM_STAGES}; @@ -264,7 +266,7 @@ int main(int argc, char* argv[]) { // Set the Vision Workbench debug level set_debug_level(debug_level); - vw_system_cache().resize( cache_size*1024*1024 ); // Set cache to 1Gb + vw_system_cache().resize( cache_size*1024*1024 ); if ( num_threads != 0 ) { std::cout << "\t--> Setting number of processing threads to: " << num_threads << "\n"; vw_settings().set_default_num_threads(num_threads); @@ -673,6 +675,18 @@ int main(int argc, char* argv[]) { write_image(out_prefix + "-GoodPixelMap.tif", disparity::missing_pixel_image(filtered_disparity_map), TerminalProgressCallback(ErrorMessage, "\t Writing: ")); DiskImageView > good_pixel_image(out_prefix + "-GoodPixelMap.tif"); + + // Okay. This is ugly. Sorry. There is some sort of annoying + // edge artifact that creeps in when we apply the NURBS to an + // edge masked image, so we create two different edge masks + // here. The first buffers by one kernel size (to prevent the + // NURBS from matching against the real bad data along the + // edges), and the second edge mask adds some addional buffering + // to get rid of the annoying edge artifact. We should really + // only need the first edge mask here. + write_image(out_prefix + "-NurbsMask.tif", select_channel(edge_mask(good_pixel_image, + PixelRGB(255,0,0), + stereo_settings().subpixel_h_kern*1.0),3)); write_image(out_prefix + "-dMask.tif", select_channel(edge_mask(good_pixel_image, PixelRGB(255,0,0), stereo_settings().subpixel_h_kern*2.0),3)); @@ -688,7 +702,9 @@ int main(int argc, char* argv[]) { DiskImageView Dmask(out_prefix + "-dMask.tif"); if(stereo_settings().fill_holes_NURBS) { std::cout << "\t--> Filling holes with bicubicly interpolated B-SPLINE surface.\n"; - hole_filled_disp_map = HoleFillView(disparity::mask(filtered_disparity_map,Dmask,Rmask), 2); + hole_filled_disp_map = HoleFillView(disparity::mask(filtered_disparity_map, + DiskImageView(out_prefix + "-NurbsMask.tif"), + Rmask), 2); } DiskImageResourceOpenEXR disparity_map_rsrc(out_prefix + "-F.exr", hole_filled_disp_map.format() ); @@ -700,6 +716,52 @@ int main(int argc, char* argv[]) { } } + /******************************************************************************/ + /* DUST_REJECTION */ + /******************************************************************************/ + if (entry_point <= DUST_REJECTION) { + if (entry_point == DUST_REJECTION) + std::cout << "\nStarting at the DUST_REJECTION stage.\n"; + vw_out(0) << "\n[ " << current_posix_time_string() << " ] : Stage 4 --> DUST_REJECTION \n"; + + try { + + // NOTE: UNCOMMENT FOR APOLLO IMAGERY ONLY!! + // + // Compute the "dust mask" for the two images. + // + // + photometric_outlier_rejection(out_prefix, stereo_settings().subpixel_h_kern); + + DiskImageView Dmask(out_prefix + "-dMask.tif"); + DiskImageView Rmask(out_prefix + "-rMask.tif"); + + DiskImageView > disparity_map(out_prefix + "-F.exr"); + ImageView > masked_disp_map = disparity::mask(disparity_map, + DiskImageView(out_prefix + "-NurbsDustMask.tif"), + Rmask); + ImageViewRef > hole_filled_disp_map = masked_disp_map; + + // Disabled at Matt's request for now.... + if(0) { + // if(stereo_settings().fill_holes_NURBS) { + std::cout << "\t--> Re-filling holes with bicubicly interpolated B-SPLINE surface.\n"; + hole_filled_disp_map = HoleFillView(disparity::mask(disparity_map, + DiskImageView(out_prefix + "-NurbsDustMask.tif"), + Rmask), 2); + } + + DiskImageResourceOpenEXR disparity_map_rsrc(out_prefix + "-FD.exr", hole_filled_disp_map.format() ); + disparity_map_rsrc.set_tiled_write(std::min(1024,hole_filled_disp_map.cols()),std::min(1024, hole_filled_disp_map.rows())); + block_write_image(disparity_map_rsrc, disparity::mask(hole_filled_disp_map,Dmask,Rmask), TerminalProgressCallback(InfoMessage, "\t--> Filtering: ") ); + + } catch (IOErr &e) { + cout << "\nUnable to start at filtering stage -- could not read input files.\n" << e.what() << "\nExiting.\n\n"; + exit(0); + } + + + } /******************************************************************************/ /* TRIANGULATION */ @@ -709,6 +771,7 @@ int main(int argc, char* argv[]) { std::cout << "\nStarting at the TRIANGULATION stage.\n"; vw_out(0) << "\n[ " << current_posix_time_string() << " ] : Stage 4 --> TRIANGULATION \n"; + try { boost::shared_ptr camera_model1, camera_model2; session->camera_models(camera_model1, camera_model2); @@ -732,7 +795,7 @@ int main(int argc, char* argv[]) { } std::string prehook_filename; - session->pre_pointcloud_hook(out_prefix+"-F.exr", prehook_filename); + session->pre_pointcloud_hook(out_prefix+"-FD.exr", prehook_filename); DiskImageView > disparity_map(prehook_filename); // Apply the stereo model. This yields a image of 3D points in