Algorithm or command line tool to decimate point cloud of terrain points?

Question

I need to take a larger (more dense) than needed list of lidar survey points (longitude, latitude, and elevation) for terrain definition and decimate it based on a 2 dimensional grid. The idea would be to end up with points based on a NxN (i.e. 1 meter x 1 meter) dimension grid using the longitude, latitude (x,y) values, therefore eliminating the points that are more than are needed. The goal is to determine what the elevation is at each point in the grid after the decimation, not use elevation as part of the decimation rule itself.

An actual or precisely structured grid is not necessary or the goal here, I only use the grid terminology to best approximate what I envision as the remainder of the cloud of points after reducing it in a manner that we have always have a point within a certain radius (i.e. 1 meter). It is possible there is a better term to use than grid.

I would like to either code/script this myself in a scripting or programming language if I can start with a decimation algorithm or use a command line tool from a project that may already exist that can do this that can run on Ubuntu and called from our application as system call. The approach should not require using a GUI based type of software or tool to solve this. It needs to be part of an automated set of steps.

The data currently exists in a tab separated values file but I could load the data into a sqlite database file if using an database/sql query driven algorithm would be better/faster. The ideal scripting language would be ruby or python but can be any really and if there exists C/C++/C# libraries for this already then we could wrap those for our needs.

Ideas?

Update Clarifying the use of the result of this decimated list: Given a user's location (known by latitude and longitude), what is the closest point in the list and in turn its elevation? We can do this now of course, but we have more data than is necessary so we just want to relax the density of the data so that if we can find the closest point within a tolerance distance (i.e. 1 meter) if able to use a decimated list vs the full list. The latitude, longitude values in the list are in decimal GPS (i.e. 38.68616190027656, -121.11013105991036)

Answer 1

PART 1: decimated version

• Load data

Load the data from the tabular file (change sep according to the separator you are using):

# installed as dependency
import pandas as pd  
# https://github.com/daavoo/pyntcloud
from pyntcloud import PyntCloud  

dense = PyntCloud(pd.read_csv("example.tsv",
                     sep='\t',
                     names=["x","y","z"]))

This is how it looks the example I created:

• Build VoxelGrid

Asuming that the latitude and longitude in your file are in meters you can generate a grid as follows:

grid_id = dense.add_structure("voxelgrid",
                           sizes=[1, 1,None],
                           bb_cuboid=False)
voxelgrid = dense.voxelgrids[grid_id]

This voxelgrid has a size of 1 along the x (latitude) and y (longitude) dimensions.

• Build decimated version

  decimated = dense.get_sample("voxelgrid_centroids", voxelgrid=grid_id)

decimated is a numpy (N,3) array. You can store it for later use in a SQL database, etc.

---

PART 2: Query

Option A: query voxelgrid

• Get mean altitudes for each grid cell

You can know get a vector with the mean z (altitude) value for each cell in the grid:

z_mean = voxelgrid.get_feature_vector(mode="z_mean")

• Query the grid with users's location:

  users_location = np.random.rand(100000, 2)

Add a column of zeros because query requires 3D (This doesn't affect the results):

users_location = np.c_[ users_location, np.zeros(users_location.shape[0]) ]

Get in wich cell each user is:

users_cell = voxelgrid.query(users_location)

And finally, get the altitude corresponding to each user:

users_altitude = z_mean[users_cell]

Option B: Use decimated version for query

• Build a KDTree of decimated:

  from scipy.spatial import cKDTree kdt = cKDTree(decimated)

• Query the KDTree with user locations:

  users_location = np.random.rand(100000, 2) users_location = np.c_[ users_location, np.zeros(users_location.shape[0])

  distances, indices = kdt.query(user_locations, k=1, n_jobs=-1)

---

Extra, you can save and laod the voxelgrid with pickle:

pickle.dump(voxelgrid, open("voxelgrid.pkl", "wb"))

voxelgrid = pickle.load(open("voxelgrid.pkl", "rb"))

Answer 2

If you have a point cloud as a text file (.xyz) a simple and fast solution is to take a random sample from the file using shuf.

10 million points in a xyz-file equals 10 million lines of text. You can run:

shuf input.xyz -l 5000000 -o out.xyz

You have decimated the file to half the original size.