Permalink
Browse files

Local homography

  • Loading branch information...
1 parent 1250176 commit 65f5d6e74746f4b6d2541dfa590e96b8fdde05bb @oleg-alexandrov oleg-alexandrov committed Apr 3, 2013
Showing with 305 additions and 48 deletions.
  1. +124 −3 src/asp/Core/LocalDisparity.cc
  2. +77 −1 src/asp/Core/LocalDisparity.h
  3. +2 −2 src/asp/Core/Makefile.am
  4. +23 −28 src/asp/Tools/stereo_corr.cc
  5. +79 −14 src/asp/Tools/stereo_rfne.cc
@@ -23,13 +23,134 @@
#include <vw/Image/Transform.h>
#include <vw/FileIO/DiskImageView.h>
#include <vw/Stereo/DisparityMap.h>
-#include <asp/Core/StereoSettings.h>
#include <asp/Core/LocalDisparity.h>
+#include <asp/Core/StereoSettings.h>
+#include <asp/Core/InterestPointMatching.h>
using namespace vw;
-using namespace vw::cartography;
namespace asp {
+ // We would like to split the numbers 0, ..., n - 1 into k buckets
+ // of approximately equal size. For example, for n = 8 and k = 3,
+ // we will have the split {0, 1, 2}, {3, 4, 5}, {6, 7}.
+
+ void split_n_into_k(int n, int k, std::vector<int> & partition){
+
+ VW_ASSERT(n >= k && k > 0,
+ ArgumentErr() << "split_n_into_k: Must have n >= k && k > 0.\n");
+ int rem = n % k;
+ int dx0 = n / k;
+
+ partition.clear();
+ int start = 0;
+ for (int i = 0; i < k; i++){
+ int dx = dx0;
+ if (rem > 0){
+ dx++;
+ rem--;
+ }
+ partition.push_back(start);
+ start += dx;
+ }
+ partition.push_back(start);
+
+ }
+
+ // Create a local homography for each correlation tile
+ void create_local_homographies(Options const& opt){
+
+ DiskImageView< PixelGray<float> > left_sub (opt.out_prefix + "-L_sub.tif");
+ DiskImageView< PixelGray<float> > left_img (opt.out_prefix + "-L.tif");
+ DiskImageView< PixelMask<Vector2i> >
+ sub_disparity(opt.out_prefix + "-D_sub.tif");
+
+ Vector2 upscale_factor( double(left_img.cols()) / double(left_sub.cols()),
+ double(left_img.rows()) / double(left_sub.rows()) );
+
+ int ts = Options::corr_tile_size();
+ int cols = (int)ceil(left_img.cols()/double(ts));
+ int rows = (int)ceil(left_img.rows()/double(ts));
+ ImageView<Matrix3x3> local_hom(cols, rows);
+
+ for (int col = 0; col < cols; col++){
+ for (int row = 0; row < rows; row++){
+
+ BBox2i bbox(col*ts, row*ts, ts, ts);
+ bbox.crop(bounding_box(left_img));
+
+ // The low-res version of bbox
+ BBox2i sub_bbox( elem_quot(bbox.min(), upscale_factor),
+ elem_quot(bbox.max(), upscale_factor) );
+
+ // Expand the box until square to make sure the local homography
+ // calculation does not fail.
+ int len = std::max(sub_bbox.width(), sub_bbox.height());
+ sub_bbox = BBox2i(sub_bbox.max() - Vector2(len, len), sub_bbox.max());
+ sub_bbox.expand(1);
+ sub_bbox.crop( bounding_box(sub_disparity) );
+
+ local_hom(col, row)
+ = homography_for_disparity(sub_bbox, crop(sub_disparity, sub_bbox) );
+ }
+ }
+
+ std::string local_hom_file = opt.out_prefix + "-local_hom.txt";
+ vw_out() << "Writing: " << local_hom_file << "\n";
+ write_local_homographies(local_hom_file, local_hom);
+
+ return;
+ }
+
+ void write_local_homographies(std::string const& local_hom_file,
+ ImageView<Matrix3x3> const& local_hom){
+
+ std::ofstream fh(local_hom_file.c_str());
+ fh.precision(18);
+ fh << local_hom.cols() << " " << local_hom.rows() << std::endl;
+
+ for (int col = 0; col < local_hom.cols(); col++){
+ for (int row = 0; row < local_hom.rows(); row++){
+ Vector<double> V = matrix_to_vector(local_hom(col, row));
+ for (int t = 0; t < int(V.size())-1; t++) fh << V[t] << " ";
+ if (V.size() > 0) fh << V[V.size()-1] << std::endl;
+ }
+ }
+ fh.close();
+
+ return;
+ }
+
+ void read_local_homographies(std::string const& local_hom_file,
+ ImageView<Matrix3x3> & local_hom){
+
+ std::ifstream fh(local_hom_file.c_str());
+ if (!fh.good())
+ vw_throw( IOErr() << "read_local_homographies: File does not exist: "
+ << local_hom_file << ".\n" );
+
+ int cols, rows;
+ if ( !(fh >> cols >> rows) )
+ vw_throw( IOErr() << "read_local_homographies: Invalid file: "
+ << local_hom_file << ".\n" );
+
+ local_hom.set_size(cols, rows);
+ for (int col = 0; col < local_hom.cols(); col++){
+ for (int row = 0; row < local_hom.rows(); row++){
+
+ Vector<double, 9> V;
+ for (int t = 0; t < int(V.size()); t++){
+ if (! (fh >> V[t]) )
+ vw_throw( IOErr() << "read_local_homographies: Invalid file: "
+ << local_hom_file << ".\n" );
+ }
+
+ local_hom(col, row) = vector_to_matrix(V);
+ }
+ }
+ fh.close();
+
+ return;
+ }
-}
+} // namespace asp
@@ -22,14 +22,90 @@
#ifndef __LOCAL_DISPARITY_H__
#define __LOCAL_DISPARITY_H__
+#include <vw/Image/ImageView.h>
+#include <vw/InterestPoint/InterestData.h>
+#include <vector>
+
// Forward declaration
namespace asp {
class Options;
}
namespace asp {
+ void split_n_into_k(int n, int k, std::vector<int> & partition);
-}
+ void create_local_homographies(Options const& opt);
+
+ void write_local_homographies(std::string const& local_hom_file,
+ vw::ImageView<vw::Matrix3x3> const& local_hom);
+ void read_local_homographies(std::string const& local_hom_file,
+ vw::ImageView<vw::Matrix3x3> & local_hom);
+
+ // Given a disparity map restricted to a subregion, find the homography
+ // transform which aligns best the two images based on this disparity.
+ template<class SeedDispT>
+ vw::math::Matrix<double> homography_for_disparity(vw::BBox2i subregion,
+ SeedDispT const& disparity){
+
+ // To do: I've never seen this, but it is possible that this
+ // homography computation may fail if the disparity in subregion
+ // has too few points or not scattered well. This function may
+ // need to be made more robust by taking in the entire disparity,
+ // not just cropped to subregion, as done now. And, we would crop
+ // the disparity to increasingly larger boxes containing
+ // subregion until the computation would succeed.
+
+ VW_ASSERT(subregion.width() == disparity.cols() &&
+ subregion.height() == disparity.rows(),
+ vw::ArgumentErr() << "homography_for_disparity: "
+ << "The sizes of subregion and disparity don't match.\n");
+
+ // We will split the subregion into N x N boxes, and average the
+ // disparity in each box, to reduce the run-time.
+ int N = 10;
+
+ std::vector<int> partitionx, partitiony;
+ split_n_into_k(disparity.cols(), std::min(disparity.cols(), N), partitionx);
+ split_n_into_k(disparity.rows(), std::min(disparity.rows(), N), partitiony);
+
+ std::vector<vw::ip::InterestPoint> left_ip, right_ip;
+ for (int ix = 0; ix < (int)partitionx.size()-1; ix++){
+ for (int iy = 0; iy < (int)partitiony.size()-1; iy++){
+
+ // First sum up the disparities in each subbox.
+ double lx = 0, ly = 0, rx = 0, ry = 0, count = 0; // int may cause overflow
+ for (int x = partitionx[ix]; x < partitionx[ix+1]; x++){
+ for (int y = partitiony[iy]; y < partitiony[iy+1]; y++){
+
+ typename SeedDispT::pixel_type disp = disparity(x, y);
+ if (!is_valid(disp)) continue;
+ lx += x; rx += (x + disp.child().x());
+ ly += y; ry += (y + disp.child().y());
+ count++;
+ }
+ }
+ if (count == 0) continue; // no valid points
+
+ // Do the averaging. We must add the box corner to the left and
+ // right interest points.
+ vw::ip::InterestPoint l, r;
+ l.x = subregion.min().x() + lx/count; r.x = subregion.min().x() + rx/count;
+ l.y = subregion.min().y() + ly/count; r.y = subregion.min().y() + ry/count;
+ left_ip.push_back(l);
+ right_ip.push_back(r);
+ }
+ }
+
+ try {
+ return homography_fit(right_ip, left_ip, bounding_box(disparity));
+ }catch ( const vw::ArgumentErr& e ){
+ // Will return the identity matrix.
+ }
+ return vw::math::identity_matrix<3>();
+
+ }
+
+} // namespace asp
#endif
View
@@ -34,11 +34,11 @@ include_HEADERS = BlobIndexThreaded.h StereoSettings.h SparseView.h \
SoftwareRenderer.h ErodeView.h $(ba_headers) Macros.h \
Common.h ThreadedEdgeMask.h GaussianClustering.h \
IntegralAutoGainDetector.h InterestPointMatching.h \
- DemDisparity.h
+ DemDisparity.h LocalDisparity.h
libaspCore_la_SOURCES = BlobIndexThreaded.cc Common.cc MedianFilter.cc \
SoftwareRenderer.cc StereoSettings.cc $(ba_sources) \
- InterestPointMatching.cc DemDisparity.cc
+ InterestPointMatching.cc DemDisparity.cc LocalDisparity.cc
libaspCore_la_LIBADD = @MODULE_CORE_LIBS@
@@ -46,8 +46,8 @@ void produce_lowres_disparity( Options & opt ) {
DiskImageView<PixelGray<float> > left_sub( opt.out_prefix+"-L_sub.tif" ),
right_sub( opt.out_prefix+"-R_sub.tif" );
- Vector2f downsample_scale( float(left_sub.cols()) / float(Lmask.cols()),
- float(left_sub.rows()) / float(Lmask.rows()) );
+ Vector2 downsample_scale( double(left_sub.cols()) / double(Lmask.cols()),
+ double(left_sub.rows()) / double(Lmask.rows()) );
DiskImageView<uint8> left_mask_sub( opt.out_prefix+"-lMask_sub.tif" ),
right_mask_sub( opt.out_prefix+"-rMask_sub.tif" );
@@ -125,11 +125,11 @@ void lowres_correlation( Options & opt ) {
// Pinhole + None
// DG + None
// Everything else should gather IP's all the time.
- float sub_scale =
- sum(elem_quot( Vector2f(file_image_size( opt.out_prefix+"-L_sub.tif" )),
- Vector2f(file_image_size( opt.out_prefix+"-L.tif" ) ) )) +
- sum(elem_quot( Vector2f(file_image_size( opt.out_prefix+"-R_sub.tif" )),
- Vector2f(file_image_size( opt.out_prefix+"-R.tif" ) ) ));
+ double sub_scale =
+ sum(elem_quot( Vector2(file_image_size( opt.out_prefix+"-L_sub.tif" )),
+ Vector2(file_image_size( opt.out_prefix+"-L.tif" ) ) )) +
+ sum(elem_quot( Vector2(file_image_size( opt.out_prefix+"-R_sub.tif" )),
+ Vector2(file_image_size( opt.out_prefix+"-R.tif" ) ) ));
sub_scale /= 4.0f;
stereo_settings().search_range =
@@ -183,7 +183,7 @@ void lowres_correlation( Options & opt ) {
}
// Create the local homographies based on D_sub
- if (stereo_settings().use_local_homography){
+ if (stereo_settings().seed_mode > 0 && stereo_settings().use_local_homography){
std::string local_hom_file = opt.out_prefix + "-local_hom.txt";
try {
ImageView<Matrix3x3> local_hom;
@@ -212,7 +212,7 @@ class SeededCorrelatorView : public ImageViewBase<SeededCorrelatorView<Image1T,
PProcT m_preproc_func;
// Settings
- Vector2f m_upscale_factor;
+ Vector2 m_upscale_factor;
BBox2i m_seed_bbox;
BBox2i m_left_image_crop_win;
stereo::CostFunctionType m_cost_mode;
@@ -234,8 +234,8 @@ class SeededCorrelatorView : public ImageViewBase<SeededCorrelatorView<Image1T,
m_sub_disparity_spread( sub_disparity_spread.impl() ),
m_local_hom(local_hom), m_preproc_func( filter.impl() ),
m_left_image_crop_win(left_image_crop_win), m_cost_mode(cost_mode) {
- m_upscale_factor[0] = float(m_left_image.cols()) / float(m_sub_disparity.cols());
- m_upscale_factor[1] = float(m_left_image.rows()) / float(m_sub_disparity.rows());
+ m_upscale_factor[0] = double(m_left_image.cols()) / m_sub_disparity.cols();
+ m_upscale_factor[1] = double(m_left_image.rows()) / m_sub_disparity.rows();
m_seed_bbox = bounding_box( m_sub_disparity );
}
@@ -351,8 +351,7 @@ class SeededCorrelatorView : public ImageViewBase<SeededCorrelatorView<Image1T,
= transform (copy_mask( m_right_image.impl(),
create_mask(m_right_mask.impl()) ),
HomographyTransform(fullres_hom),
- m_left_image.impl().cols(),
- m_left_image.impl().rows());
+ m_left_image.impl().cols(), m_left_image.impl().rows());
right_trans_img = apply_mask(right_trans_masked_img);
right_trans_mask
= channel_cast_rescale<uint8>(select_channel(right_trans_masked_img, 1));
@@ -391,11 +390,7 @@ class SeededCorrelatorView : public ImageViewBase<SeededCorrelatorView<Image1T,
stereo_settings().corr_kernel, m_cost_mode,
stereo_settings().xcorr_threshold,
stereo_settings().corr_max_levels );
- return prerasterize_type
- (transform_disparities(do_round, bbox, inverse(fullres_hom),
- crop(corr_view.prerasterize(bbox), bbox)),
- -bbox.min().x(), -bbox.min().y(),
- cols(), rows() );
+ return corr_view.prerasterize(bbox);
}else{
typedef stereo::PyramidCorrelationView<Image1T, Image2T, Mask1T, Mask2T, PProcT> CorrView;
CorrView corr_view( m_left_image, m_right_image,
@@ -460,6 +455,8 @@ void stereo_correlation( Options& opt ) {
// Load up for the actual native resolution processing
DiskImageView<PixelGray<float> > left_disk_image(opt.out_prefix+"-L.tif"),
right_disk_image(opt.out_prefix+"-R.tif");
+ DiskImageView<vw::uint8> Lmask(opt.out_prefix + "-lMask.tif"),
+ Rmask(opt.out_prefix + "-rMask.tif");
ImageViewRef<PixelMask<Vector2i> > sub_disparity;
if ( stereo_settings().seed_mode > 0 )
sub_disparity =
@@ -468,39 +465,37 @@ void stereo_correlation( Options& opt ) {
if ( stereo_settings().seed_mode == 2 )
sub_disparity_spread =
DiskImageView<PixelMask<Vector2i> >(opt.out_prefix+"-D_sub_spread.tif");
- ImageViewRef<PixelMask<Vector2i> > fullres_disparity;
ImageView<Matrix3x3> local_hom;
- if ( stereo_settings().use_local_homography ){
+ if ( stereo_settings().seed_mode > 0 && stereo_settings().use_local_homography ){
std::string local_hom_file = opt.out_prefix + "-local_hom.txt";
read_local_homographies(local_hom_file, local_hom);
}
- DiskImageView<vw::uint8> Lmask(opt.out_prefix + "-lMask.tif"),
- Rmask(opt.out_prefix + "-rMask.tif");
-
stereo::CostFunctionType cost_mode;
if (stereo_settings().cost_mode == 0) cost_mode = stereo::ABSOLUTE_DIFFERENCE;
else if (stereo_settings().cost_mode == 1) cost_mode = stereo::SQUARED_DIFFERENCE;
else if (stereo_settings().cost_mode == 2) cost_mode = stereo::CROSS_CORRELATION;
else
- vw_throw( ArgumentErr() << "Unknown value " << stereo_settings().cost_mode << " for cost-mode.\n" );
+ vw_throw( ArgumentErr() << "Unknown value " << stereo_settings().cost_mode
+ << " for cost-mode.\n" );
+ ImageViewRef<PixelMask<Vector2i> > fullres_disparity;
if ( stereo_settings().pre_filter_mode == 2 ) {
vw_out() << "\t--> Using LOG pre-processing filter with "
<< stereo_settings().slogW << " sigma blur.\n";
fullres_disparity =
seeded_correlation( left_disk_image, right_disk_image, Lmask, Rmask,
sub_disparity, sub_disparity_spread, local_hom,
- stereo::LaplacianOfGaussian(stereo_settings().slogW), opt.left_image_crop_win,
- cost_mode );
+ stereo::LaplacianOfGaussian(stereo_settings().slogW),
+ opt.left_image_crop_win, cost_mode );
} else if ( stereo_settings().pre_filter_mode == 1 ) {
vw_out() << "\t--> Using Subtracted Mean pre-processing filter with "
<< stereo_settings().slogW << " sigma blur.\n";
fullres_disparity =
seeded_correlation( left_disk_image, right_disk_image, Lmask, Rmask,
sub_disparity, sub_disparity_spread, local_hom,
- stereo::SubtractedMean(stereo_settings().slogW), opt.left_image_crop_win,
- cost_mode );
+ stereo::SubtractedMean(stereo_settings().slogW),
+ opt.left_image_crop_win, cost_mode );
} else {
vw_out() << "\t--> Using NO pre-processing filter." << std::endl;
fullres_disparity =
Oops, something went wrong.

0 comments on commit 65f5d6e

Please sign in to comment.