Browse files

Implement threshold for masking dark pixels

  • Loading branch information...
1 parent c9899c6 commit 4de161763a58e7bfcfa3fa8bc4070b671291a861 @oleg-alexandrov oleg-alexandrov committed Dec 11, 2012
View
2 src/asp/Core.h
@@ -16,7 +16,7 @@
// __END_LICENSE__
-/// \file stereo.cc
+/// \file Core.h
///
#include <asp/asp_config.h>
View
2 src/asp/Core/ErodeView.h
@@ -93,4 +93,4 @@ class ErodeView : public vw::ImageViewBase<ErodeView<ViewT> > {
};
-#endif//__INPAINT_H__
+#endif//__ERODEVIEW_H__
View
123 src/asp/Core/InpaintView.h
@@ -16,11 +16,11 @@
// __END_LICENSE__
-/// \file Inpaint.h
+/// \file InpaintView.h
///
-#ifndef __INPAINT_H__
-#define __INPAINT_H__
+#ifndef __INPAINTVIEW_H__
+#define __INPAINTVIEW_H__
// Standard
#include <vector>
@@ -48,6 +48,8 @@ namespace asp {
vw::ImageViewBase<SourceT> const& m_view;
blob::BlobCompressed m_c_blob;
+ bool m_use_grassfire;
+ typename SourceT::pixel_type m_default_inpaint_val;
SparseView<SparseT> & m_patches;
int m_id;
boost::shared_ptr<vw::Mutex> m_crop;
@@ -56,21 +58,31 @@ namespace asp {
public:
InpaintTask( vw::ImageViewBase<SourceT> const& view,
blob::BlobCompressed const& c_blob,
+ bool use_grassfire,
+ typename SourceT::pixel_type default_inpaint_val,
SparseView<SparseT> & sparse,
int const& id,
boost::shared_ptr<vw::Mutex> crop,
boost::shared_ptr<vw::Mutex> insert ) :
- m_view(view), m_c_blob(c_blob), m_patches(sparse), m_id(id), m_crop(crop), m_insert(insert) {}
+ m_view(view), m_c_blob(c_blob),
+ m_use_grassfire(use_grassfire), m_default_inpaint_val(default_inpaint_val),
+ m_patches(sparse), m_id(id), m_crop(crop), m_insert(insert) {}
void operator()() {
vw_out(vw::VerboseDebugMessage,"inpaint") << "Task " << m_id << ": started\n";
// Gathering information about blob
vw::BBox2i bbox = m_c_blob.bounding_box();
- bbox.expand(10);
+ if (m_use_grassfire){
+ // I don't understand why this is necessary
+ bbox.expand(10);
+ }else{
+ bbox.expand(1);
+ }
+
// How do we want to handle spots on the edges?
if ( bbox.min().x() < 0 || bbox.min().y() < 0 ||
- bbox.max().x() > m_view.impl().cols() || bbox.max().y() > m_view.impl().rows() ) {
+ bbox.max().x() >= m_view.impl().cols() || bbox.max().y() >= m_view.impl().rows() ) {
vw_out(vw::VerboseDebugMessage,"inpaint") << "Task " << m_id << ": early exiting\n";
return;
}
@@ -88,42 +100,50 @@ namespace asp {
cropped_copy = crop( m_view, bbox );
}
- // Creating binary image to highlight hole (then make grassfire)
+ // Creating binary image to highlight hole
vw::ImageView<vw::uint8> mask( bbox.width(), bbox.height() );
fill( mask, 0 );
for ( std::list<vw::Vector2i>::const_iterator iter = blob.begin();
iter != blob.end(); iter++ )
mask( iter->x(), iter->y() ) = 255;
- vw::ImageView<vw::int32> distance = grassfire(mask);
- int max_distance = max_pixel_value( distance );
-
- // Working out order of convolution
- std::list<vw::Vector2i> processing_order;
- for ( int d = 1; d < max_distance+1; d++ )
- for ( int i = 0; i < bbox.width(); i++ )
- for ( int j = 0; j < bbox.height(); j++ )
- if ( distance(i,j) == d ) {
- processing_order.push_back( vw::Vector2i(i,j) );
+
+ if (m_use_grassfire){
+
+ vw::ImageView<vw::int32> distance = grassfire(mask);
+ int max_distance = max_pixel_value( distance );
+
+ // Working out order of convolution
+ std::list<vw::Vector2i> processing_order;
+ for ( int d = 1; d < max_distance+1; d++ )
+ for ( int i = 0; i < bbox.width(); i++ )
+ for ( int j = 0; j < bbox.height(); j++ )
+ if ( distance(i,j) == d ) {
+ processing_order.push_back( vw::Vector2i(i,j) );
+ }
+
+ // Iterate and apply convolution seperately to each channel
+ for ( vw::uint32 c = 0; c < vw::PixelNumChannels<typename SourceT::pixel_type>::value;
+ c++ )
+ for ( int d = 0; d < 10*max_distance*max_distance; d++ )
+ for ( std::list<vw::Vector2i>::const_iterator iter = processing_order.begin();
+ iter != processing_order.end(); iter++ ) {
+ float sum = 0;
+ sum += .176765*cropped_copy(iter->x()-1,iter->y()-1)[c];
+ sum += .176765*cropped_copy(iter->x()-1,iter->y()+1)[c];
+ sum += .176765*cropped_copy(iter->x()+1,iter->y()-1)[c];
+ sum += .176765*cropped_copy(iter->x()+1,iter->y()+1)[c];
+ sum += .073235*cropped_copy(iter->x()+1,iter->y())[c];
+ sum += .073235*cropped_copy(iter->x()-1,iter->y())[c];
+ sum += .073235*cropped_copy(iter->x(),iter->y()+1)[c];
+ sum += .073235*cropped_copy(iter->x(),iter->y()-1)[c];
+ cropped_copy(iter->x(),iter->y())[c] = sum;
}
-
- // Iterate and apply convolution seperately to each channel
- for ( vw::uint32 c = 0; c < vw::PixelNumChannels<typename SourceT::pixel_type>::value;
- c++ )
- for ( int d = 0; d < 10*max_distance*max_distance; d++ )
- for ( std::list<vw::Vector2i>::const_iterator iter = processing_order.begin();
- iter != processing_order.end(); iter++ ) {
- float sum = 0;
- sum += .176765*cropped_copy(iter->x()-1,iter->y()-1)[c];
- sum += .176765*cropped_copy(iter->x()-1,iter->y()+1)[c];
- sum += .176765*cropped_copy(iter->x()+1,iter->y()-1)[c];
- sum += .176765*cropped_copy(iter->x()+1,iter->y()+1)[c];
- sum += .073235*cropped_copy(iter->x()+1,iter->y())[c];
- sum += .073235*cropped_copy(iter->x()-1,iter->y())[c];
- sum += .073235*cropped_copy(iter->x(),iter->y()+1)[c];
- sum += .073235*cropped_copy(iter->x(),iter->y()-1)[c];
- cropped_copy(iter->x(),iter->y())[c] = sum;
- }
-
+ }else{
+ for ( std::list<vw::Vector2i>::const_iterator iter = blob.begin();
+ iter != blob.end(); iter++ )
+ cropped_copy( iter->x(), iter->y() ) = m_default_inpaint_val;
+ }
+
{ // Insert results into sparse view
vw::Mutex::Lock lock( *m_insert );
m_patches.absorb(bbox.min(),copy_mask(cropped_copy,create_mask( mask, 0 )));
@@ -136,7 +156,7 @@ namespace asp {
} // end namespace inpaint_p
- /// InpaintView (feed all blobs before hand )
+ /// InpaintView (feed all blobs before hand)
/// Prerasterize -> do nothing
/// Constructor -> Perform all processing spawn own threads
/// Rasterize -> See if in blob area, then return pix in location,
@@ -147,12 +167,17 @@ namespace asp {
ViewT m_child;
SparseView<typename vw::UnmaskedPixelType<typename ViewT::pixel_type>::type> m_patches;
-
+ bool m_use_grassfire;
+ typename ViewT::pixel_type m_default_inpaint_val;
+
// A special copy constructor for prerasterization
template <class OViewT>
InpaintView( vw::ImageViewBase<ViewT> const& image,
- InpaintView<OViewT> const& other ) :
- m_child(image.impl()), m_patches(other.m_patches) {}
+ InpaintView<OViewT> const& other,
+ bool use_grassfire,
+ typename ViewT::pixel_type default_inpaint_val) :
+ m_child(image.impl()), m_patches(other.m_patches),
+ m_use_grassfire(use_grassfire), m_default_inpaint_val(default_inpaint_val) {}
public:
typedef typename vw::UnmaskedPixelType<typename ViewT::pixel_type>::type sparse_type;
@@ -161,7 +186,11 @@ namespace asp {
typedef vw::ProceduralPixelAccessor<InpaintView<ViewT> > pixel_accessor;
InpaintView( vw::ImageViewBase<ViewT> const& image,
- BlobIndexThreaded const& bindex ) : m_child(image.impl()) {
+ BlobIndexThreaded const& bindex,
+ bool use_grassfire,
+ typename ViewT::pixel_type default_inpaint_val): m_child(image.impl()),
+ m_use_grassfire(use_grassfire),
+ m_default_inpaint_val(default_inpaint_val){
{
vw::Stopwatch sw;
sw.start();
@@ -174,10 +203,12 @@ namespace asp {
for ( unsigned i = 0; i < bindex.num_blobs(); i++ ) {
boost::shared_ptr<task_type> task(new task_type(image, bindex.compressed_blob(i),
+ m_use_grassfire, m_default_inpaint_val,
m_patches, i, crop_mutex, insert_mutex ));
queue.add_task( task );
}
queue.join_all();
+ sw.stop();
vw_out(vw::VerboseDebugMessage,"inpaint") << "Time used in inpaint threads: " << sw.elapsed_seconds() << "s\n";
}
}
@@ -198,7 +229,7 @@ namespace asp {
typedef InpaintView<typename ViewT::prerasterize_type> prerasterize_type;
inline prerasterize_type prerasterize( vw::BBox2i const& bbox ) const {
typename ViewT::prerasterize_type preraster = m_child.prerasterize(bbox);
- return prerasterize_type( preraster, *this );
+ return prerasterize_type( preraster, *this, m_use_grassfire, m_default_inpaint_val );
}
template <class DestT>
inline void rasterize( DestT const& dest, vw::BBox2i const& bbox ) const {
@@ -212,10 +243,12 @@ namespace asp {
template <class SourceT>
inline InpaintView<SourceT> inpaint( vw::ImageViewBase<SourceT> const& src,
- BlobIndexThreaded const& bindex) {
- return InpaintView<SourceT>(src,bindex);
+ BlobIndexThreaded const& bindex,
+ bool use_grassfire,
+ typename SourceT::pixel_type default_inpaint_val) {
+ return InpaintView<SourceT>(src, bindex, use_grassfire, default_inpaint_val);
}
} //end namespace asp
-#endif//__INPAINT_H__
+#endif//__INPAINTVIEW_H__
View
15 src/asp/Core/StereoSettings.cc
@@ -48,14 +48,27 @@ namespace asp {
// Define our options that are available
// ----------------------------------------------------------
PreProcessingDescription::PreProcessingDescription() : po::options_description("Preprocessing Options") {
+
StereoSettings& global = stereo_settings();
+
+ double nan = std::numeric_limits<double>::quiet_NaN();
+ global.nodata_threshold = nan;
+ global.nodata_percentage = nan;
+ global.nodata_optimal_threshold_factor = nan;
+
(*this).add_options()
("alignment-method", po::value(&global.alignment_method),
"Rough alignment for input images. [Homography, Epipolar, None]")
("individually-normalize", po::bool_switch(&global.individually_normalize),
"Individually normalize the input images between 0.0-1.0 using +- 2.5 sigmas about their mean values.")
("force-use-entire-range", po::bool_switch(&global.force_max_min),
- "Normalize images based on the global min and max values from both images. Don't both using this option if you are using normalized cross correlation.");
+ "Normalize images based on the global min and max values from both images. Don't use this option if you are using normalized cross correlation.")
+ ("nodata-threshold", po::value(&global.nodata_threshold),
+ "Pixels with values less than this number times the maximum image value are treated as no-data.")
+ ("nodata-percentage", po::value(&global.nodata_percentage),
+ "The percentage of (low-value) pixels treated as no-data.")
+ ("nodata-optimal-threshold-factor", po::value(&global.nodata_optimal_threshold_factor),
+ "Pixels with values less than this factor times the optimal Otsu threshold are treated as no-data. Suggested value: 0.1 to 0.2.");
}
CorrelationDescription::CorrelationDescription() : po::options_description("Correlation Options") {
View
6 src/asp/Core/StereoSettings.h
@@ -70,7 +70,11 @@ namespace asp {
// individually with their
// own hi's and lo's
bool force_max_min; // Use entire dynamic range of image.
-
+ double nodata_threshold; // Pixels with value less than this are treated as no-data
+ double nodata_percentage; // the percentage of low-value pixels treated as no-data
+ double nodata_optimal_threshold_factor; // Pixels with values less than this factor times the optimal Otsu threshold
+ // are treated as no-data
+
// Correlation Options
float slogW; // Preprocessing filter width
vw::uint16 pre_filter_mode; // 0 = None
View
5 src/asp/Tools/stereo.cc
@@ -33,6 +33,11 @@ namespace asp {
boost::program_options::options_description const&
additional_options ) {
+ // Print the command being run in debug mode.
+ std::string run_cmd = "";
+ for (int s = 0; s < argc; s++) run_cmd += std::string(argv[s]) + " ";
+ VW_OUT(DebugMessage, "stereo") << "\n\n" << run_cmd << "\n";
+
po::options_description general_options_sub("");
general_options_sub.add_options()
("session-type,t", po::value(&opt.stereo_session_string), "Select the stereo session type to use for processing. [options: pinhole isis dg rpc]")
View
6 src/asp/Tools/stereo_corr.cc
@@ -354,8 +354,8 @@ void produce_lowres_disparity( int32 cols, int32 rows, Options const& opt ) {
Vector2f down_sample_scale( float(left_sub.cols()) / float(cols),
float(left_sub.rows()) / float(rows) );
- DiskImageView<uint8> left_mask( opt.out_prefix+"-lMask_sub.tif" ),
- right_mask( opt.out_prefix+"-rMask_sub.tif" );
+ DiskImageView<uint8> left_mask_sub( opt.out_prefix+"-lMask_sub.tif" ),
+ right_mask_sub( opt.out_prefix+"-rMask_sub.tif" );
BBox2i search_range( floor(elem_prod(down_sample_scale,stereo_settings().search_range.min())),
ceil(elem_prod(down_sample_scale,stereo_settings().search_range.max())) );
@@ -372,7 +372,7 @@ void produce_lowres_disparity( int32 cols, int32 rows, Options const& opt ) {
asp::block_write_gdal_image( opt.out_prefix + "-D_sub.tif",
remove_outliers(
stereo::pyramid_correlate( left_sub, right_sub,
- left_mask, right_mask,
+ left_mask_sub, right_mask_sub,
stereo::LaplacianOfGaussian(1.4),
search_range,
stereo_settings().kernel,
View
7 src/asp/Tools/stereo_fltr.cc
@@ -93,14 +93,17 @@ void write_good_pixel_and_filtered( ImageViewBase<ImageT> const& inputview,
opt, TerminalProgressCallback("asp", "\t--> Good Pxl Map: ") );
}
- // Fill Holes
+ // Fill holes
if(!stereo_settings().disable_fill_holes) {
vw_out() << "\t--> Filling holes with Inpainting method.\n";
BlobIndexThreaded bindex( invert_mask( inputview.impl() ),
stereo_settings().fill_hole_max_size );
vw_out() << "\t * Identified " << bindex.num_blobs() << " holes\n";
+ bool use_grassfire = true;
+ typename ImageT::pixel_type default_inpaint_val;
asp::block_write_gdal_image( opt.out_prefix + "-F.tif",
- asp::InpaintView<ImageT >(inputview.impl(), bindex ),
+ inpaint(inputview.impl(), bindex,
+ use_grassfire, default_inpaint_val),
opt, TerminalProgressCallback("asp","\t--> Filtering: ") );
} else {
View
71 src/asp/Tools/stereo_pprc.cc
@@ -22,7 +22,9 @@
#include <asp/Tools/stereo.h>
#include <asp/Core/ThreadedEdgeMask.h>
+#include <asp/Core/InpaintView.h>
#include <vw/Cartography/GeoTransform.h>
+#include <vw/Math/Functors.h>
using namespace vw;
using namespace asp;
@@ -80,6 +82,36 @@ resample_aa( ImageViewBase<ViewT> const& input, double factor ) {
int32(.5+(input.impl().rows()*factor)) );
}
+struct MaskAboveThreshold: public ReturnFixedType< PixelMask<uint8> > {
+ double m_threshold;
+
+ MaskAboveThreshold(double threshold): m_threshold(threshold){}
+
+ PixelMask<uint8> operator() (PixelGray<float> const& pix) const {
+ if (pix >= m_threshold)
+ return PixelMask<uint8>(255);
+ else
+ return PixelMask<uint8>();
+ }
+};
+template <class ImageT>
+UnaryPerPixelView<ImageT, MaskAboveThreshold>
+inline mask_above_threshold( ImageViewBase<ImageT> const& image, double threshold ) {
+ return UnaryPerPixelView<ImageT, MaskAboveThreshold>( image.impl(),
+ MaskAboveThreshold(threshold) );
+}
+
+ImageViewRef< PixelMask<uint8> > mask_and_fill_holes( ImageViewRef< PixelGray<float> > const& img,
+ double threshold ){
+ // Create the mask of pixels above threshold. Fix any holes in it.
+ ImageViewRef< PixelMask<uint8> > thresh_mask = mask_above_threshold(img, threshold);
+ int max_area = 0; // fill arbitrarily big holes
+ bool use_grassfire = false; // fill with default value
+ PixelMask<uint8> default_inpaint_val = uint8(255);
+ BlobIndexThreaded bindex( invert_mask( thresh_mask.impl() ), max_area );
+ return inpaint(thresh_mask.impl(), bindex, use_grassfire, default_inpaint_val);
+}
+
void stereo_preprocessing( Options& opt ) {
vw_out() << "\n[ " << current_posix_time_string()
@@ -118,8 +150,47 @@ void stereo_preprocessing( Options& opt ) {
ImageViewRef< PixelMask<uint8> > right_mask = copy_mask(constant_view(uint8(255),
right_image.cols(), right_image.rows() ),
asp::threaded_edge_mask(right_image,0,0,1024));
+
+ double left_threshold = stereo_settings().nodata_threshold;
+ double right_threshold = stereo_settings().nodata_threshold;
+ double nodata_fraction = stereo_settings().nodata_percentage/100.0;
+ double nodata_factor = stereo_settings().nodata_optimal_threshold_factor;
+ if (int(!std::isnan(left_threshold)) +
+ int(!std::isnan(nodata_fraction)) +
+ int(!std::isnan(nodata_factor)) >= 2
+ ){
+ vw_throw( ArgumentErr()
+ << "\nAt most one of the no-data settings "
+ << "(threshold, percentage, or optimal threshold factor) must be set.\n");
+ }
+
+ if ( !std::isnan(nodata_factor) ){
+ // Find the black pixels threshold using Otsu's optimal threshold method.
+ left_threshold = nodata_factor*optimal_threshold(left_image);
+ right_threshold = nodata_factor*optimal_threshold(right_image);
+ }
+ if ( !std::isnan(nodata_fraction) ){
+ // Declare a fixed proportion of pixels to be black.
+ math::CDFAccumulator< PixelGray<float> > left_cdf(1024, 1024), right_cdf(1024, 1024);
+ for_each_pixel( left_image, left_cdf );
+ for_each_pixel( right_image, right_cdf );
+
+ left_threshold = left_cdf.quantile(nodata_fraction);
+ right_threshold = right_cdf.quantile(nodata_fraction);
+ }
+
+ if ( !std::isnan(left_threshold) && !std::isnan(right_threshold) ){
+ // Mask pixels below threshold.
+
+ ImageViewRef< PixelMask<uint8> > left_thresh_mask = mask_and_fill_holes(left_image, left_threshold);
+ left_mask = intersect_mask(left_mask, left_thresh_mask);
+
+ ImageViewRef< PixelMask<uint8> > right_thresh_mask = mask_and_fill_holes(right_image, right_threshold);
+ right_mask = intersect_mask(right_mask, right_thresh_mask);
+ }
+
bool has_left_georef = read_georeference(left_georef, opt.in_file1);
bool has_right_georef = read_georeference(right_georef, opt.in_file2);
if (has_left_georef && has_right_georef){

0 comments on commit 4de1617

Please sign in to comment.