Permalink
Browse files

point2las and point2mesh: Support point clouds of size bw 3 and 6

  • Loading branch information...
1 parent bc1ded2 commit 0ea574cfc92b377bed71e47c7d8dc1b1719b7e3f @oleg-alexandrov oleg-alexandrov committed Mar 1, 2013
Showing with 112 additions and 93 deletions.
  1. +4 −4 graveyard/point2mesh.cc
  2. +7 −20 src/asp/Tools/point2dem.cc
  3. +74 −13 src/asp/Tools/point2dem.h
  4. +21 −31 src/asp/Tools/point2las.cc
  5. +6 −25 src/asp/Tools/point2mesh.cc
View
@@ -77,9 +77,9 @@ namespace vw {
template <class ViewT>
-BBox<float,3> point_image_bbox(ImageViewBase<ViewT> const& point_image) {
+BBox<3> pointcloud_bbox(ImageViewBase<ViewT> const& point_image) {
// Compute bounding box
- BBox<float,3> result;
+ BBox<3> result;
typename ViewT::pixel_accessor row_acc = point_image.impl().origin();
for( int32 row=0; row < point_image.impl().rows(); ++row ) {
typename ViewT::pixel_accessor col_acc = row_acc;
@@ -193,13 +193,13 @@ int main( int argc, char *argv[] ) {
ImageViewRef<Vector3> point_image = point_disk_image;
if (vm.count("center")) {
- BBox<float,3> bbox = point_image_bbox(point_disk_image);
+ BBox<float,3> bbox = pointcloud_bbox(point_disk_image);
std::cout << "\t--> Centering model around the origin.\n";
std::cout << "\t Initial point image bounding box: " << bbox << "\n";
Vector3 midpoint = (bbox.max() + bbox.min()) / 2.0;
std::cout << "\t Midpoint: " << midpoint << "\n";
point_image = point_image_offset(point_image, -midpoint);
- BBox<float,3> bbox2 = point_image_bbox(point_image);
+ BBox<float,3> bbox2 = pointcloud_bbox(point_image);
std::cout << "\t Re-centered point image bounding box: " << bbox2 << "\n";
}
View
@@ -73,7 +73,6 @@ struct Options : asp::BaseOptions {
BBox2 target_projwin;
BBox2i target_projwin_pixels;
uint32 fsaa;
- uint32 pix_size; // How many channels there are in PC
// Output
std::string out_prefix, output_file_type;
@@ -465,14 +464,16 @@ void do_software_rasterization( const ImageViewBase<ViewT>& proj_point_input,
// Write triangulation error image if requested
if ( opt.do_error ) {
- if (opt.pix_size == 4){
+ int num_channels = get_num_channles(opt.pointcloud_filename);
+
+ if (num_channels == 4){
// The error is a scalar.
DiskImageView<Vector4> point_disk_image(opt.pointcloud_filename);
ImageViewRef<double> error_channel = select_channel(point_disk_image,3);
rasterizer.set_texture( error_channel );
rasterizer_fsaa = generate_fsaa_raster( rasterizer, opt );
save_image(opt, rasterizer_fsaa, georef, "IntersectionErr");
- }else{
+ }else if (num_channels == 6){
// The error is a 3D vector. Convert it to NED coordinate system,
// and rasterize it.
DiskImageView<Vector6> point_disk_image(opt.pointcloud_filename);
@@ -492,6 +493,8 @@ void do_software_rasterization( const ImageViewBase<ViewT>& proj_point_input,
asp::combine_channels(opt.nodata_value,
rasterized[0], rasterized[1], rasterized[2]),
georef, "IntersectionErr");
+ }else{
+ vw_throw( ArgumentErr() << "Expecting the input point cloud to have points of size 4 or 6.");
}
}
@@ -530,23 +533,7 @@ int main( int argc, char *argv[] ) {
try {
handle_arguments( argc, argv, opt );
- {
- boost::scoped_ptr<SrcImageResource> src(DiskImageResource::open(opt.pointcloud_filename));
- int num_channels = src->channels();
- int num_planes = src->planes();
- opt.pix_size = num_channels*num_planes;
- VW_ASSERT( (opt.pix_size == 4 || opt.pix_size == 6),
- ArgumentErr() << "Expecting the input point cloud to have points of size 4 or 6.");
- }
-
- ImageViewRef<Vector3> point_image;
- if (opt.pix_size == 4){
- DiskImageView<Vector4> point_disk_image(opt.pointcloud_filename);
- point_image = select_points(point_disk_image);
- }else{
- DiskImageView<Vector6> point_disk_image(opt.pointcloud_filename);
- point_image = select_points(point_disk_image);
- }
+ ImageViewRef<Vector3> point_image = read_n_channels<3>(opt.pointcloud_filename);
// Apply an (optional) rotation to the 3D points before building the mesh.
if (opt.phi_rot != 0 || opt.omega_rot != 0 || opt.kappa_rot != 0) {
View
@@ -21,6 +21,9 @@
/// This header represents the overflow of small objects and image
/// transforms that point2dem specifically applies.
+// To do: Many of these utilities may need to go to a lower-level header file
+// dealing with point cloud manipulations.
+
#include <stdlib.h>
#include <vw/FileIO.h>
@@ -31,7 +34,7 @@
namespace vw {
// Erases a file suffix if one exists and returns the base string
- static std::string prefix_from_pointcloud_filename(std::string const& filename) {
+ std::string prefix_from_pointcloud_filename(std::string const& filename) {
std::string result = filename;
// First case: filenames that match <prefix>-PC.<suffix>
@@ -111,18 +114,76 @@ namespace vw {
PointTransFunc(t));
}
- // Imageview operator that extracts only the first 3 channels of the
- // point cloud. The remaining channels are the point cloud error
- // (scalar or vector).
- template <class VectorT>
- struct SelectPoints : public ReturnFixedType<Vector3> {
- Vector3 operator() (VectorT const& pt) const { return subvector(pt,0,3); }
+ // Imageview operator that extracts only the first n channels of an
+ // image with m channels.
+ template <int n, int m>
+ struct SelectPoints : public ReturnFixedType< Vector<double, n> > {
+ Vector<double, n> operator() (Vector<double, m> const& pt) const { return subvector(pt,0,n); }
};
-
- template <class VectorT>
- UnaryPerPixelView<DiskImageView<VectorT>, SelectPoints<VectorT> >
- inline select_points( ImageViewBase<DiskImageView<VectorT> > const& image ) {
- return UnaryPerPixelView<DiskImageView<VectorT>, SelectPoints<VectorT> >( image.impl(),
- SelectPoints<VectorT>() );
+
+ template <int n, int m>
+ UnaryPerPixelView<DiskImageView< Vector<double, m> >, SelectPoints<n, m> >
+ inline select_points( ImageViewBase<DiskImageView< Vector<double, m> > > const& image ) {
+ return UnaryPerPixelView<DiskImageView< Vector<double, m> >, SelectPoints<n, m> >( image.impl(),
+ SelectPoints<n, m>() );
+ }
+
+ int get_num_channles(std::string filename){
+ boost::scoped_ptr<SrcImageResource> src(DiskImageResource::open(filename));
+ int num_channels = src->channels();
+ int num_planes = src->planes();
+ return num_channels*num_planes;
}
+
+ // Given an image with each pixel a vector of size m, return the
+ // first n channels of that image. We must have 1 <= n <= m <= 6.
+ template<int n>
+ ImageViewRef< Vector<double, n> > read_n_channels(std::string filename){
+
+ int max_m = 6;
+ int m = get_num_channles(filename);
+
+ VW_ASSERT( 1 <= n,
+ ArgumentErr() << "Attempting to read " << n << " channels from an image.");
+ VW_ASSERT( n <= m,
+ ArgumentErr() << "Attempting to read " << n << " channels from an image with "
+ << m << " channels.");
+ VW_ASSERT( m <= max_m,
+ NoImplErr() << "Reading from images with more than "
+ << max_m << " channels is not implemented.");
+
+ ImageViewRef< Vector<double, n> > out_image;
+ if (m == 1) out_image = select_points<n, 1>(DiskImageView< Vector<double, 1> >(filename));
+ else if (m == 2) out_image = select_points<n, 2>(DiskImageView< Vector<double, 2> >(filename));
+ else if (m == 3) out_image = select_points<n, 3>(DiskImageView< Vector<double, 3> >(filename));
+ else if (m == 4) out_image = select_points<n, 4>(DiskImageView< Vector<double, 4> >(filename));
+ else if (m == 5) out_image = select_points<n, 5>(DiskImageView< Vector<double, 5> >(filename));
+ else if (m == 6) out_image = select_points<n, 6>(DiskImageView< Vector<double, 6> >(filename));
+
+ return out_image;
+ }
+
+ // Helper function that extracts the bounding box of a point cloud. It
+ // skips the point Vector3(), which, being at the center of the
+ // planet, is an invalid point.
+
+ template <class ViewT>
+ BBox3 pointcloud_bbox(ImageViewBase<ViewT> const& point_image) {
+ // Compute bounding box
+ BBox3 result;
+ typename ViewT::pixel_accessor row_acc = point_image.impl().origin();
+ for( int32 row=0; row < point_image.impl().rows(); ++row ) {
+ typename ViewT::pixel_accessor col_acc = row_acc;
+ for( int32 col=0; col < point_image.impl().cols(); ++col ) {
+ if (*col_acc != Vector3())
+ result.grow(*col_acc);
+ col_acc.next_col();
+ }
+ row_acc.next_row();
+ }
+ return result;
+ }
+
+
+
}
View
@@ -27,6 +27,7 @@
#include <asp/Core/Macros.h>
#include <asp/Core/Common.h>
+#include <asp/Tools/point2dem.h> // We share common functions with point2dem
#include <vw/FileIO.h>
#include <vw/Image.h>
@@ -37,8 +38,10 @@ namespace po = boost::program_options;
// Allows FileIO to correctly read/write these pixel types
namespace vw {
+ typedef Vector<float64,6> Vector6;
template<> struct PixelFormatID<Vector3> { static const PixelFormatEnum value = VW_PIXEL_GENERIC_3_CHANNEL; };
template<> struct PixelFormatID<Vector4> { static const PixelFormatEnum value = VW_PIXEL_GENERIC_4_CHANNEL; };
+ template<> struct PixelFormatID<Vector6> { static const PixelFormatEnum value = VW_PIXEL_GENERIC_6_CHANNEL; };
}
struct Options : asp::BaseOptions {
@@ -92,32 +95,19 @@ int main( int argc, char *argv[] ) {
try {
handle_arguments( argc, argv, opt );
- DiskImageView<Vector4> point_disk_image(opt.pointcloud_filename);
+ ImageViewRef<Vector3> point_image = read_n_channels<3>(opt.pointcloud_filename);
+ BBox3 cloud_bbox = pointcloud_bbox(point_image);
- // Find the bounding box of the points
- BBox3 point_boundaries;
- TerminalProgressCallback progress_bar("asp","Statistics: ");
- for ( int y = 0; y < point_disk_image.rows(); y++ ) {
- progress_bar.report_fractional_progress(y,point_disk_image.rows());
- for ( int x = 0; x < point_disk_image.cols(); x++ ) {
- Vector3 point = subvector(point_disk_image(x, y), 0, 3);
- if ( point == Vector3() ) continue; // skip no-data points
-
- point_boundaries.grow( point );
- }
- }
- progress_bar.report_finished();
-
- // The las format stores the values as 32 bit integers. So, if we
- // want to store a point x, we store round((x-offset)/scale), as
- // well as the offset and scale values. Here we decide the values
- // for offset and scale to lose minimum amount of precision. We
- // make the scale almost as large as it can be without causing
- // integer overflow.
- Vector3 offset = (point_boundaries.min() + point_boundaries.max())/2.0;
+ // The las format stores the values as 32 bit integers. So, for a
+ // given point, we store round((point-offset)/scale), as well as
+ // the offset and scale values. Here we decide the values for
+ // offset and scale to lose minimum amount of precision. We make
+ // the scale almost as large as it can be without causing integer
+ // overflow.
+ Vector3 offset = (cloud_bbox.min() + cloud_bbox.max())/2.0;
double maxInt = std::numeric_limits<int32>::max();
maxInt *= 0.95; // Just in case stay a bit away
- Vector3 scale = point_boundaries.size()/(2.0*maxInt);
+ Vector3 scale = cloud_bbox.size()/(2.0*maxInt);
liblas::Header header;
header.SetDataFormatId(liblas::ePointFormat1);
@@ -132,24 +122,24 @@ int main( int argc, char *argv[] ) {
lasFile = opt.out_prefix + ".las";
}
- vw_out() << "Writing LAS File: " << lasFile + "\n";
+ vw_out() << "Writing LAS file: " << lasFile + "\n";
std::ofstream ofs;
ofs.open(lasFile.c_str(), std::ios::out | std::ios::binary);
liblas::Writer writer(ofs, header);
- progress_bar.set_progress_text("Integerizing:");
- for (int row = 0; row < point_disk_image.rows(); row++){
- progress_bar.report_fractional_progress(row, point_disk_image.rows());
- for (int col = 0; col < point_disk_image.cols(); col++){
+ TerminalProgressCallback progress_bar("asp","LAS: ");
+ for (int row = 0; row < point_image.rows(); row++){
+ progress_bar.report_fractional_progress(row, point_image.rows());
+ for (int col = 0; col < point_image.cols(); col++){
- Vector3 point = subvector(point_disk_image(col, row), 0, 3);
+ Vector3 point = point_image(col, row);
if ( point == Vector3() ) continue; // skip no-data points
- point = elem_quot((point - offset),scale);
+ point = round( elem_quot((point - offset), scale) );
liblas::Point las_point;
- las_point.SetCoordinates(point[0], point[1], point[2]); // Takes doubles. Unclear if it scales the data.
+ las_point.SetCoordinates(point[0], point[1], point[2]);
writer.WritePoint(las_point);
}
@@ -28,7 +28,7 @@
#include <math.h>
//VisionWorkbench & ASP
-#include <asp/Tools/point2dem.h> // We share common functions with p2dem
+#include <asp/Tools/point2dem.h> // We share common functions with point2dem
#include <asp/Core/Macros.h>
#include <asp/Core/Common.h>
using namespace vw;
@@ -57,8 +57,10 @@ namespace po = boost::program_options;
// Allows FileIO to correctly read/write these pixel types
namespace vw {
+ typedef Vector<float64,6> Vector6;
template<> struct PixelFormatID<Vector3> { static const PixelFormatEnum value = VW_PIXEL_GENERIC_3_CHANNEL; };
template<> struct PixelFormatID<Vector4> { static const PixelFormatEnum value = VW_PIXEL_GENERIC_4_CHANNEL; };
+ template<> struct PixelFormatID<Vector6> { static const PixelFormatEnum value = VW_PIXEL_GENERIC_6_CHANNEL; };
}
struct Options : asp::BaseOptions {
@@ -413,26 +415,6 @@ osg::Node* build_mesh( vw::ImageViewBase<ViewT> const& point_image,
}
-// Helper function that extracts the point box of a point cloud.
-template <class ViewT>
-BBox<float,3> point_image_bbox(ImageViewBase<ViewT> const& point_image) {
- // Compute bounding box
- BBox<float,3> result;
- typename ViewT::pixel_accessor row_acc = point_image.impl().origin();
- for( int32 row=0; row < point_image.impl().rows(); ++row ) {
- typename ViewT::pixel_accessor col_acc = row_acc;
- for( int32 col=0; col < point_image.impl().cols(); ++col ) {
- if (*col_acc != Vector3())
- result.grow(*col_acc);
- col_acc.next_col();
- }
- row_acc.next_row();
- }
- return result;
-}
-
-
-
// MAIN
// ---------------------------------------------------------
@@ -489,13 +471,12 @@ int main( int argc, char *argv[] ){
try {
handle_arguments( argc, argv, opt );
- // Loading point cloud!
- DiskImageView<Vector4> point_disk_image(opt.pointcloud_filename);
- ImageViewRef<Vector3> point_image = select_points(point_disk_image);
+ // Loading point cloud
+ ImageViewRef<Vector3> point_image = read_n_channels<3>(opt.pointcloud_filename);
// Centering Option (helpful if you are experiencing round-off error...)
if (opt.center) {
- BBox<float,3> bbox = point_image_bbox(select_points(point_disk_image));
+ BBox3 bbox = pointcloud_bbox(point_image);
vw_out() << "\t--> Centering model around the origin.\n";
vw_out() << "\t Initial point image bounding box: " << bbox << "\n";
Vector3 midpoint = (bbox.max() + bbox.min()) / 2.0;

0 comments on commit 0ea574c

Please sign in to comment.