NeoGeographyToolkit / StereoPipeline

The NASA Ames Stereo Pipeline is a suite of automated geodesy & stereogrammetry tools designed for processing planetary imagery captured from orbiting and landed robotic explorers on other planets.
Apache License 2.0
478 stars 168 forks source link

point2las output triangulation error as .las attribute #339

Closed bpurinton closed 2 years ago

bpurinton commented 2 years ago

I am interested in doing my own point cloud filtering prior to point2dem. For this I would like the raw x,y,z points without triangulation percentile filtering currently implemented in point2las (which I can turn off with "--remove-outliers-params 100 1"), BUT I would still like this triangulation error information to be included in the output *.las file. Then I could independently utilize the triangulation error for each point to create bespoke filtering algorithms for cleaning up the point clouds e.g., by weighting points based on triangulation errors.

As it stands, point2las only outputs the x,y,z information from the -PC.tif file. I propose a command line option to include a fourth attribute in the .las file: the fourth band of the *.tif, the triangulation error at each point. So you get out an x,y,z,e file; where e is triangulation error.

Would this be easy / possible to implement?

oleg-alexandrov commented 2 years ago

This is a good idea and something which was proposed before. I put it on my to-do list to implement this next week.

Maybe that is too late for you, but regretfully here we have funds only for a handful of projects, each with specific goals, and anything else which gets accomplished, including even bug fixing, is when there's time left.

You can also try to pull points straight from the PC.tif file. It can be done in Python with the GDAL package. The PC file is a four-band TIF, having x, y, z, and error. The x, y, z values need to have a shift added to them, which is in the PC file tif header. All this will be slow with Python but doable if there's no other option. And note that an (x, y, z) pixel with value (0, 0, 0) is declared invalid.

bpurinton commented 2 years ago

Glad to hear it isn‘t just me. I can wait for the feature, no problem. The slow Python+GDAL way is what I wanted to avoid. As always, thanks!

oleg-alexandrov commented 2 years ago

I added the ability to save the triangulation error to the las file. Things are somewhat constrained by the fact that LAS 1.2 (which is what we support so far) has no room for storing float values or 32-bit values apart from x, y, z which we need for the cloud proper. All it can do is store an intensity field using two bytes (https://www.asprs.org/a/society/committees/standards/asprs_las_format_v12.pdf).

The obvious thing to do is to multiply the non-negative triangulation error by some scale factor, round it to int, and cast it to uint16, with care to set all values larger than the largest uint16 to that value before casting, to avoid junk.

This is not ideal, but usable, as likely the triangulation errors are no more than 1-60 meters, and if the scale factor is, for example, 1000, that still allows for storing the triangulation error with submillimeter precision.

Hence I added to point2las the following option:

--triangulation-error-factor <float (default: 0)> If this factor is positive, save the point cloud triangulation error to the 2-byte LAS intensity field by storing min(round(factor*error), 65535). Resulting values that equal 65535 should be treated with caution.

This will be in the build at https://github.com/NeoGeographyToolkit/StereoPipeline/releases sometime over the weekend (the date will change appropriately).

I tested this, by using las2txt on the obtained las file with the option '--parse xyzi' and then creating a dem from that text file with point2dem and that last intensity value as height. I get an image very similar one would get otherwise with point2dem --errorimage from the original PC.tif file, save for the multiplying factor.

bpurinton commented 2 years ago

Great! Thanks for the update, this is useful

oleg-alexandrov commented 2 years ago

Sure. The daily build finally made it (it was held up mostly by unrelated things).