Usage of libpointmatcher at NASA Ames and notes on large memory datasets

[oleg-alexandrov]

Hi Francois and Stephane,

I've mentioned earlier that in my group at NASA Ames we use libpointmatcher for aligning 3D terrains on Earth, Moon, and in the future Mars. Since this is not an application you probably had in mind when you wrote the software, I thought I'd describe a bit how we use it and also mention how it scales for large data, and how I tweaked it to scale better.

First, a bit of background. The larger project I work on is Ames Stereo Pipeline (http://ti.arc.nasa.gov/tech/asr/intelligent-robotics/ngt/stereo/, https://github.com/NeoGeographyToolkit/StereoPipeline; the alignment tool using libpointmatcher is at https://github.com/NeoGeographyToolkit/StereoPipeline/blob/master/src/asp/Tools/pc_align.cc).

This stereo software is used to do 3D terrain reconstruction from two views obtained by a satellite. Historically it was used for Mars and Moon, but recently it is being used also for Earth (although in full honesty it could work better on Earth, it handles well dead stuff, like rock and ice, but not so well vegetation, urban environments, clouds, and water).

Such 3D terrains are usually pretty accurate, but not always perfectly aligned, for example due to errors in camera position and orientation in space. Luckily one can use ground control points, which are a sparse accurate dataset of measurements, to align the 3D terrain obtained from stereo to ground truth. This even works on the Moon, where a spacecraft named LRO (https://en.wikipedia.org/wiki/Lunar_Reconnaissance_Orbiter) has an instrument named LOLA shooting laser pulses at the moon and getting accurate but sparse 3D measurements. libpointmatcher is being used for such dense from stereo to sparse and accurate 3D alignment.

libpointmatcher works awesomely well, with both translations and rotations (point-to-plane beats point-to-point ICP for rotations, but for small translations they work about the same).

Two applications we already use libpointmatcher for is studying the flow of ice in Greenland (the work is done at the University of Washington), and alignment of 3D terrain on Moon to LRO data mentioned above.

I made some tweaks to libpointmatcher to handle better large datasets. The alignment tool now uses 8 GB of RAM for a terrain with 100 million points (we'd wish we could handle terrains say 25 times as large). I'd like to describe the (rather simple optimizations I did.

• Resampling is happening on the fly as data is being loaded form disk, a given point is either kept or tossed.
• The libpointmatcher input.features datastructure is pre-allocated in advance, rather than using push-back into vectors and then copy.
• All intermediate copies of data are removed.
• Shift to the centroid of the reference data points is done in-place.
• Normals are computed only for point-to-plane and the identity filter is used for point-to-point (I don't use the yaml file and the user does not get to choose the filter type, the software makes the choices).
• Filtering is done in-place. For example, for SamplingSurfaceDataPointsFilter, the normals are just attached to existing datastructure rather than creating a new point cloud with normals. I implemented this only when resampling ratio is 1 (so no resampling) but it is easy to do even resampling in-place if necessary. Here's the diff of this change: https://github.com/oleg-alexandrov/libpointmatcher/commit/747f8da9f9744a361bf99a65b1c5b988f5babe07

More work would need to be done to bring this code in-line with your style conventions and other filters would have to be modified likewise for consistency, but I hope this diff is enough to give an idea of what could be done if you'd ever decide to support larger datasets.

There's only one more optization which could be done, but it would be tricky and I won't do it. libpointmatcher could be templated by two data types, say PointMatcher<T, S>. T would be the working variables type, say double, and S would be the input data type, say float. This because having the input data in float format is enough most of the time (by doing a careful initial shift), while float is not enough for doing accurate calculations.

This was a long post. Thanks for reading. And thanks again for the great library.

[oleg-alexandrov]

In the second to last paragraph the html markup stole some text, I meant PointMatcher<T,S > instead of existing PointMatcher<T>

[stephanemagnenat]

Dear Oleg,

Thank you very much for keeping us informed of your use of libpointmatcher. It is warming that our software is useful to you, and I am very happy to be able to bring a tiny contribution to planetary science.

I have to discuss with François, but in principle we are very willing to reintegrate your changes into our master branch. Could you provide a pull request against the latest git master with your changes? This would allow code review and ease the integration process.

kind regards

[stephanemagnenat]

A note about very large point clouds: to save memory, when representing kd-trees in memory, libnabo stores the cut dimension and the index of the cutting point in a single 32-bit word, which in 3D allows for about 1 billion points in the cloud. We could have a huge version, but it would double the memory footprint of the kd-tree. Come back to me if you need that, we would discuss how to make it happen.

[pomerlef]

That is very good to know! I'm curious to know how did you came across the library?

We can link to your project on our README page (you can do the same on your side if you wish).

As for the feature request, we indeed didn't have use cases for such large point clouds but it is very nice to know that we fulfil (mostly) your needs. I will put in place a wish list for the next release but I need to say that I would like to have more documentation/tutorials first. We will try to find more manpower internally but it seems to be harder than expected.

[oleg-alexandrov]

(I commented here earlier but somehow my comment went missing. Odd.)

I will spend some time next week trying to make a robust and easy to merge pull request for in-place filtering. Other changes I made would be harder to integrate, for example, shifting the reference point cloud by the center of mass in-place would require changes to the library interface which may be undesirable.

I found libpointmatcher from Google search. BTW, I added a link to libpointmatcher and a citation of your paper to the Wikipedia ICP article (http://en.wikipedia.org/wiki/Iterative_closest_point). (That article could by way benefit from some attention from experts.)

Thank you for linking to our project. I added a citation to your 2013 paper and a link to github's project page in our documentation where we describe the alignment tool (http://byss.arc.nasa.gov/stereopipeline/daily_build/asp_book.pdf). Soon that will be part of the official release docs.

[oleg-alexandrov]

Sorry for the delay. I made a pull request for my changes making all filters in DataPointsFiltersImpl.cpp work in-place, with one less temporary point cloud. The pull request is at

https://github.com/ethz-asl/libpointmatcher/pull/17

The name of the new function is inPlaceFilter(...). I kept the original filter() function for backward compatibility but it just does a call to the new inPlaceFilter() as it can be seen in the code.

I tested carefully all changes I made. Particularly,SamplingSurfaceNormalDataPointsFilter required three test-cases, as that code has several possible paths. My testsuite, containing test cases for all filters, is in libpointmatcher/examples/tmptest/test_all.sh. It can be run before and after my changes. In retrospect, I should not have put that testsuite in the same pull request. It may perhaps be useful for you but it would not belong in the final

One important note. The new implementation of SamplingSurfaceNormalDataPointsFilter returns the same results as before, but perhaps in different order (otherwise it is hard to copy the data in-place). The ICP algorithms don't care about order of data in the clouds,and give identical results as before.

However, the order of data returned by SamplingSurfaceNormalDataPointsFilter matters if another filter is applied on top of its output, such as ShadowDataPointsFilter, so the result of this second filter is slightly different if is applied on top of the new SamplingSurfaceNormalDataPointsFilter. I would not worry about this as this effect is quite expected. In my tests above, I avoid applying another filter on top of SamplingSurfaceNormalDataPointsFilter, and I test all filters which need a precomputed normal by using instead SurfaceNormalDataPointsFilter().

Let me know what you think so far. Thanks!

[oleg-alexandrov]

Closing this as the discussion moved to #17