Introduction

In recent years, there was great progress in the development of LIDAR detectors that resulted in huge amounts of raw data. LIDARs are now more accurate and provide a much higher resolution than even 10 years ago. Aerial-based LIDARs have become an efficient method for receiving information about the environment. Nevertheless, the data you get is in fact only a collection of sparse points that may provide some nice visualization but require extensive processing when it comes to more practical purposes. Unfortunately, as of now, the progress of computer vision mostly concerns structured two-dimensional data (photo, video). The current methods do not generalize to multidimensional sparse data like Point Cloud that we receive after a basic LIDAR data pre-processing. Extensive research is being done in the field. We

should particularly mention PCL – a graf of libraries developed by a great international community to provide instruments for both 2D and 3D data for a wide variety of applications. Unfortunately, at the moment it is not an easy task to apply the library to some solution of interest. It often means that you have to dive into the library too deep. It makes sense for production-grade products that need high scalability. But it may be too costly for a PoC development. Therefore, I decided to try what can be done with point cloud data using a simple approach and pretty standard Python libraries (PCL can be used from Python but only so far, since only small subsets can be integrated seamlessly).

Data for the Experiment

As an example let’s use the data generated by an aerial-based LIDAR for the detection of power lines. Power lines are often clearly visible in point cloud visualization. However, mapping the points that belong to a power line requires a lot of manual efforts. On the other hand, simple geometrical consideration may provide us with means to greatly simplify or even automate such processing.

The power line on the picture is actually a collection of points which form certain geometrical patterns. To simplify further categorization, I decided to check how we can form clusters out of these points.

Software Used for the Experiment

I will use NumPy, Sklearn, Laspy and SciPy libraries for the purpose of cluster formation and matplotlib for visualization.

import laspy
import scipy
import numpy as np
import matplotlib.pyplot as plt
from sklearn.cluster import DBSCAN
from sklearn import metrics
from sklearn import preprocessing
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import path

Laspy is great for handling point cloud data in Python. We read point cloud data from a las file and check the shape of the actual dataset.

# Open a file in read mode:
inFile = laspy.file.File(“./LAS/simple.las”)
# Grab a numpy dataset of our clustering dimensions:
dataset = np.vstack([inFile.x, inFile.y, inFile.z]).transpose()
dataset.shape

(5942479, 3) — our point cloud consists of 5942479 points. It is not enough if you want to get to small details. But the number is too big if you try to convert this DataFrame into a three-dimensional NumPy array, as in this case, we will get a huge 5942479³ = 2.09*10²⁰ array. It will use an enormous amount of RAM for the storage of really sparse data. The obvious thing to do is to use a NumPy sparse array. But as a matter of fact, a sparse array works great for 2D but fails with 3D data. Matrix manipulation functions are not fully compatible with sparse 3D matrices. We have to stick to working with a DataFrame instead of a NumPy array due to memory requirements.

Eliminating Out-of-Scope Points

We need to find a method to eliminate those points that are not power lines. Power lines are placed high above the ground for safety reasons and to ensure their optimal performance. But for the rugged terrain, we have to take into account that some ground points can be higher than power lines in different parts of the image due to the ground incline. To avoid this let’s divide our point cloud to small vertical parts.

%%time
def frange(start, stop, step):
  i = start
  while i < stop:
    yield i
    i += step
#ground points grid filter
n = 100 #grid step
dataset_Z_filtered = dataset[[0]]
zfiltered = (dataset[:, 2].max() — dataset[:, 2].min())/10 #setting height filtered from ground
print(‘zfiltered =’, zfiltered)
xstep = (dataset[:, 0].max() — dataset[:, 0].min())/n
ystep = (dataset[:, 1].max() — dataset[:, 1].min())/n
for x in frange (dataset[:, 0].min(), dataset[:, 0].max(), xstep):
    for y in frange (dataset[:, 1].min(), dataset[:, 1].max(), ystep):
     datasetfiltered = dataset[(dataset[:,0] > x)
                             &(dataset[:, 0] < x+xstep)
                             &(dataset[:, 1] > y)
                             &(dataset[:, 1] < y+ystep)]
      if datasetfiltered.shape[0] > 0:
      datasetfiltered = datasetfiltered[datasetfiltered[:, 2]
                        >(datasetfiltered[:, 2].min()+ zfiltered)]
        if datasetfiltered.shape[0] > 0:
        dataset_Z_filtered = np.concatenate((dataset_Z_filtered,
                                               datasetfiltered))
print(‘dataset_Z_filtered shape’, dataset_Z_filtered.shape)

With the help of this simple approach, we can greatly reduce the number of points in the cloud in no time even using moderate computing power. In our case, it was the reduction of the number of points by an order of magnitude in 3 minutes — not bad for a few lines of code, given that we made no real efforts for any optimization.

dataset_Z_filtered shape (169862, 3)
CPU times: user 3min 16s, sys: 7.14 ms, total: 3min 16s
Wall time: 3min 16s

Now we’ll be using a much smaller filtered dataset.

Exploring Our Data

Let’s explore our data:

print(“Examining Point Format: “)
pointformat = inFile.point_format
for spec in inFile.point_format:
print(spec.name)
Examining Point Format:
X
Y
Z
intensity
flag_byte
raw_classification
scan_angle_rank
user_data
pt_src_id
gps_time

During my experiments, I try to use the 4D representation of data (X, Y, Z, and intensity) but the results do not improve over 3D (X, Y, Z) so let’s stick to the latter subset of data.

print(‘Z range =’, dataset[:, 2].max() — dataset[:, 2].min())
print(‘Z max =’, dataset[:, 2].max(), ‘Z min =’, dataset[:, 2].min())
print(‘Y range =’, dataset[:, 1].max() — dataset[:, 1].min())
print(‘Y max =’, dataset[:, 1].max(), ‘Y min =’, dataset[:, 1].min())
print(‘X range =’, dataset[:, 0].max() — dataset[:, 0].min())
print(‘X max =’, dataset[:, 0].max(), ‘X min =’, dataset[:, 0].min())
Z range = 149.81
Z max = 181.78999908447264 Z min = 31.979999084472652
Y range = 622.9700000002049
Y max = 2576396.509974365 Y min = 2575773.539974365
X range = 556.4400000000605
X max = 711882.7199987792 X min = 711326.2799987792

Data Normalization

As you can see, the values are in different ranges. For better results, we should normalize the dataset.

dataset = preprocessing.normalize(dataset)

Clustering the Points That Are Geometrically Close

Now we are ready to process our data. Our power lines are actually spatial clusters of points so it is natural to try a clustering algorithm. After a short investigation, I found that DBSCAN provided by the Sklearn library works best out-of-the-box.

clustering = DBSCAN(eps=2, min_samples=5, leaf_size=30).fit(dataset)

Now let’s visualize our results.

core_samples_mask = np.zeros_like(clustering.labels_, dtype=bool)
core_samples_mask[clustering.core_sample_indices_] = True
labels = clustering.labels_
# Number of clusters in labels, ignoring noise if present.
n_clusters_ = len(set(labels)) — (1 if -1 in labels else 0)
n_noise_ = list(labels).count(-1)
print(‘Estimated number of clusters: %d’ % n_clusters_)
print(‘Estimated number of noise points: %d’ % n_noise_)
Estimated number of clusters: 501
Estimated number of noise points: 1065

Visualization

Most of our points were grouped into clusters. Let’s see what it looks like in practice:

# Black removed and is used for noise instead.
fig = plt.figure(figsize=[100, 50])
ax = fig.add_subplot(111, projection=’3d’)
unique_labels = set(labels)
colors = [plt.cm.Spectral(each)
  for each in np.linspace(0, 1, len(unique_labels))]
  for k, col in zip(unique_labels, colors):
   if k == -1:
# Black used for noise.
col = [0, 0, 0, 1]
  class_member_mask = (labels == k)
  xyz = dataset[class_member_mask & core_samples_mask]
  ax.scatter(xyz[:, 0], xyz[:, 1], xyz[:, 2], c=col, marker=”.”)
plt.title(‘Estimated number of cluster: %d’ % n_clusters_)
plt.show()

It is clear now that simple geometrical consideration and a pretty standard clustering method help us to simplify the point categorization using moderate computing resources.

Each cluster of points can be categorized separately if needed.

Results

Our experiment showed that a combination of geometrical consideration and standard Python libraries can result in a significant reduction of efforts needed to categorize raw point cloud data for further usage.

Author:
Alex Simkiv, Data Scientist|Co-founder at MindCraft
Information Technology & Data Science