Skip to content
Browse files

wv_correct: Make multi-threaded

  • Loading branch information...
1 parent b40d42e commit 4816a53fbd29d64f0108c8cc054fe140245cef9f @oleg-alexandrov oleg-alexandrov committed Dec 17, 2013
Showing with 148 additions and 42 deletions.
  1. +148 −42 src/asp/Tools/wv_correct.cc
View
190 src/asp/Tools/wv_correct.cc
@@ -20,11 +20,34 @@
///
// Correct CCD artifacts in WorldView2 images with TDI 16.
-
// To do:
-// 1. Rm unnecessary digits from x and y offset values.
-// 2. Make the code multi-threaded.
-// 3. Add documentation.
+// Verify that the input image satisfies the above assumptions.
+// Add documentation.
+// Print default offsets.
+
+// The problem: A WV2 image is obtained by mosaicking from left to
+// right image blocks which are as tall is the entire image (each
+// block comes from an individual CCD image sensor). Blocks are
+// slightly misplaced in respect to each other by some unknown
+// subpixel offset. The goal of this tool is to locate these CCD
+// artifact boundary locations, and undo the offsets.
+//
+
+// Observations:
+// 1. The CCD artifact locations are periodic, but the starting offset
+// is not positioned at exactly one period. It is less by one fixed
+// value which we name 'shift'.
+// 2. There are CCD offsets in both x and y at each location.
+// 3. The magnitude of all CCD offsets in x is the same, but their sign
+// alternates. The same is true in y.
+// 4. The period of CCD offsets is inversely proportional to the detector
+// pitch.
+// 5. The CCD offsets are pretty consistent among all images of given
+// scan direction (forward or reverse).
+// We used all these and a lot of images to heuristically find the
+// period and offset of the artifacts, and the 'shift' value. We
+// correct these below. We allow the user to override the value of CCD
+// offsets if desired.
#include <vw/FileIO.h>
#include <vw/Image.h>
@@ -43,6 +66,7 @@ using namespace vw;
using namespace asp;
using namespace vw::cartography;
using namespace xercesc;
+using namespace std;
// Allows FileIO to correctly read/write these pixel types
namespace vw {
@@ -54,17 +78,23 @@ namespace vw {
}
struct Options : asp::BaseOptions {
- // Input
double xoffset, yoffset;
double xoffset_forward, yoffset_forward;
double xoffset_reverse, yoffset_reverse;
-
std::string camera_image_file, camera_model_file, output_image;
Options(){}
};
void handle_arguments( int argc, char *argv[], Options& opt ) {
+ // These quantities were heuristically obtained by averaging
+ // over a large set of runs while removing outliers and giving
+ // more weight to more reliable data.
+ double default_xoffset_forward = 0.2842;
+ double default_yoffset_forward = 0.2369;
+ double default_xoffset_reverse = 0.3396;
+ double default_yoffset_reverse = 0.3725;
+
po::options_description general_options("");
general_options.add_options()
("xoffset", po::value(&opt.xoffset), "Specify the CCD offset correction to apply in the x direction.")
@@ -92,11 +122,6 @@ void handle_arguments( int argc, char *argv[], Options& opt ) {
<< "in order to proceed.\n\n"
<< usage << general_options );
-
- double default_xoffset_forward = 0.284177215039376;
- double default_yoffset_forward = 0.23692707872382;
- double default_xoffset_reverse = 0.339580131804236;
- double default_yoffset_reverse = 0.372509828912177;
if (vm.count("xoffset")){
opt.xoffset_forward = opt.xoffset;
opt.xoffset_reverse = opt.xoffset;
@@ -114,14 +139,91 @@ void handle_arguments( int argc, char *argv[], Options& opt ) {
}
+template <class ImageT>
+class WVCorrectView: public ImageViewBase< WVCorrectView<ImageT> >{
+ ImageT m_img;
+ double m_shift, m_period, m_xoffset, m_yoffset;
+ typedef typename ImageT::pixel_type PixelT;
+
+public:
+ WVCorrectView( ImageT const& img,
+ double shift, double period, double xoffset, double yoffset):
+ m_img(img), m_shift(shift), m_period(period),
+ m_xoffset(xoffset), m_yoffset(yoffset){}
+
+ typedef PixelT pixel_type;
+ typedef PixelT result_type;
+ typedef ProceduralPixelAccessor<WVCorrectView> pixel_accessor;
+
+ inline int32 cols() const { return m_img.cols(); }
+ inline int32 rows() const { return m_img.rows(); }
+ inline int32 planes() const { return 1; }
+
+ inline pixel_accessor origin() const { return pixel_accessor( *this, 0, 0 ); }
+
+ inline pixel_type operator()( double/*i*/, double/*j*/, int32/*p*/ = 0 ) const {
+ vw_throw(NoImplErr() << "WVCorrectView::operator()(...) is not implemented");
+ return pixel_type();
+ }
+
+ typedef CropView<ImageView<pixel_type> > prerasterize_type;
+ inline prerasterize_type prerasterize(BBox2i const& bbox) const {
+
+ // Need to see a bit more of the input image for the purpose
+ // of interpolation.
+ int bias = (int)ceil(std::max(std::abs(m_xoffset), std::abs(m_yoffset)))
+ + BilinearInterpolation::pixel_buffer + 1;
+ BBox2i biased_box = bbox;
+ biased_box.expand(bias);
+ biased_box.crop(bounding_box(m_img));
+
+ ImageView<result_type> cropped_img = crop(m_img, biased_box);
+ InterpolationView<EdgeExtensionView< ImageView<result_type>, ValueEdgeExtension<result_type> >, BilinearInterpolation> interp_img
+ = interpolate(cropped_img, BilinearInterpolation(),
+ ValueEdgeExtension<result_type>(result_type()));
+
+ ImageView<result_type> tile(bbox.width(), bbox.height());
+ for (int col = bbox.min().x(); col < bbox.max().x(); col++){
+
+ // The sign of CCD offsets alternates as one moves along the image
+ // columns. As such, at "even" blocks, the offsets accumulated so
+ // far cancel each other, so we need to correct the "odd" blocks
+ // only.
+ int block_index = (int)floor((col - m_shift)/m_period);
+ double valx = 0, valy = 0;
+ if (block_index % 2 ){
+ valx = -m_xoffset;
+ valy = -m_yoffset;
+ }
+
+ for (int row = bbox.min().y(); row < bbox.max().y(); row++){
+ tile(col - bbox.min().x(), row - bbox.min().y() )
+ = interp_img(col - biased_box.min().x() + valx,
+ row - biased_box.min().y() + valy);
+ }
+ }
+
+ return prerasterize_type(tile, -bbox.min().x(), -bbox.min().y(),
+ cols(), rows() );
+ }
+
+ template <class DestT>
+ inline void rasterize(DestT const& dest, BBox2i bbox) const {
+ vw::rasterize(prerasterize(bbox), dest, bbox);
+ }
+};
+template <class ImageT>
+WVCorrectView<ImageT> wv_correct(ImageT const& img, double shift,
+ double period, double xoffset, double yoffset){
+ return WVCorrectView<ImageT>(img, shift, period, xoffset, yoffset);
+}
+
int main( int argc, char *argv[] ) {
Options opt;
try {
handle_arguments( argc, argv, opt );
- std::cout << "data is " << opt.camera_image_file << ' ' << opt.camera_model_file << ' ' << opt.output_image << std::endl;
-
GeometricXML geo;
AttitudeXML att;
EphemerisXML eph;
@@ -145,9 +247,14 @@ int main( int argc, char *argv[] ) {
vw_throw( ArgumentErr() << "XML file \"" << opt.camera_model_file << "\" is invalid.\n" );
}
- double period = 5.64/geo.detector_pixel_pitch;
+ // The first CCD artifact is at column period + shift,
+ // then they repeat with given period.
double shift = -30.0;
+ double period = 5.64/geo.detector_pixel_pitch;
+ // The offsets at the first CCD artifact location.
+ // The offsets keep the same magnitude but their sign
+ // alternates as one moves along image columns.
double xoffset, yoffset;
if (scan_dir == "forward"){
xoffset = opt.xoffset_forward;
@@ -157,44 +264,43 @@ int main( int argc, char *argv[] ) {
yoffset = opt.yoffset_reverse;
}
- vw_out() << "Using xoffset: " << xoffset << std::endl;
- vw_out() << "Using yoffset: " << yoffset << std::endl;
+ vw_out() << "Using x offset: " << xoffset << std::endl;
+ vw_out() << "Using y offset: " << yoffset << std::endl;
// Internal sign adjustments
if (scan_dir == "forward"){
yoffset = -yoffset;
}else{
xoffset = -xoffset;
}
- std::cout << "final values: " << xoffset << ' ' << yoffset << std::endl;
- std::cout << "period is " << period << std::endl;
- ImageView<float> D = copy(DiskImageView<float>(opt.camera_image_file));
- ImageView<float> E = copy(DiskImageView<float>(opt.camera_image_file));
-
- InterpolationView<EdgeExtensionView< ImageView<float>, ValueEdgeExtension<float> >, BilinearInterpolation> interp_E
- = interpolate(E, BilinearInterpolation(),
- ValueEdgeExtension<float>(0));
-
- int num_cols = D.cols();
- int num_rows = D.rows();
-
- for (int col = 0; col < num_cols; col++){
- int n = (int)floor((col - shift)/period);
- double valx = 0, valy = 0;
- if (n%2 == 1){
- valx = -xoffset;
- valy = -yoffset;
- }
-
- for (int row = 0; row < num_rows; row++){
- D(col, row) = interp_E(col + valx, row + valy);
- }
+ DiskImageView<float> input_img(opt.camera_image_file);
+ bool has_nodata = false;
+ double nodata = numeric_limits<double>::quiet_NaN();
+ boost::shared_ptr<DiskImageResource> img_rsrc
+ ( new DiskImageResourceGDAL(opt.camera_image_file) );
+ if (img_rsrc->has_nodata_read()){
+ has_nodata = true;
+ nodata = img_rsrc->nodata_read();
}
-
+
vw_out() << "Writing: " << opt.output_image << std::endl;
- asp::block_write_gdal_image(opt.output_image, D, opt,
- TerminalProgressCallback("asp", "\t-->: "));
+ if (has_nodata){
+ asp::block_write_gdal_image(opt.output_image,
+ apply_mask
+ (wv_correct(create_mask(input_img,
+ nodata),
+ shift, period, xoffset, yoffset),
+ nodata),
+ nodata, opt,
+ TerminalProgressCallback("asp", "\t-->: "));
+ }else{
+ asp::block_write_gdal_image(opt.output_image,
+ wv_correct(input_img, shift, period,
+ xoffset, yoffset),
+ opt,
+ TerminalProgressCallback("asp", "\t-->: "));
+ }
} ASP_STANDARD_CATCHES;

0 comments on commit 4816a53

Please sign in to comment.
Something went wrong with that request. Please try again.