Building detection: e.g. coplanar points detection

[jacksonvoelkel]

First off - this is a fantastic package. I’m truly blown away by the power of it. I work mostly in urban forestry, and am wondering if it is feasible to implement any building / coplanar detection functions! I’ve seen studies that do this with a Douglas-Peucker on a flightline-by-flightine basis, though I imagine there are far more efficient was to do this.

Again, thank you for this!

[Jean-Romain]

I have these two papers in my bibliography database. I added them with the hope to implement building detection someday. But so far I don't have time to do that.

• Borkowski, A. (2016). ISPRS Journal of Photogrammetry and Remote Sensing 3D building reconstruction from ALS data using unambiguous decomposition into elementary structures, 118, 1–12. https://doi.org/10.1016/j.isprsjprs.2016.04.005
• Limberger, F. A., & Oliveira, M. M. (2015). Real-time detection of planar regions in unorganized point clouds. Pattern Recognition, 48(6), 2043–2053. https://doi.org/10.1016/j.patcog.2014.12.020

What you can do to contribute is providing me small datasets with buildings and/or link to academic papers related to building detection. But if I implement such feature it will be on free time and for fun. So it is not likely to happen in a close future.

[bi0m3trics]

@Jean-Romain What do you consider a "small dataset"? One that can be used on cran or one that contains a geographically small area? I can easily throw together urban and wildland samples (with or without powerlines). I tried to implement the Limberger & Oliveira work ( it's what pdal uses for its filters.approximatecoplanar) but it was way beyond my capabilities so I have to give up... I did however compile and run their code (available here, released under the GNU GPLv3) and play around with the parameters using my own datasets and I got somewhat nice (visual) results.

[Jean-Romain]

@bi0m3trics I mean datasets that contain buildings and are sharable by email (e.g up to 5 Mb). I don't have a single building in my own datasets so can't even make tests.

Thank for the link to the original source code. This may significantly reduce the time required to make an implementation of this algorithm. I could create a new package as I did for the CSF algorithm. Sadly the code seems associated with useless openGL code and needs a cleaning that is likely to be a not hard but annoying step. Also the source code comes with example datasets.

Anyway your link helps me so much. I now know that it is possible to add such feature in a not so far future.

[bi0m3trics]

Here are two "building/powerline" datasets for anyone to play/build code to work with... the classifications were done (quickly, and cleaned up a bit) with GlobalMapper v20 and they were part of a Coconino County acquisition in 2016 (EPSG: 26912). I figured it might help to have classified files to compare against, but be aware that there are clear misclassifications that I didn't bother to clean up...

• WilliamsAZ_Wildland.las
• WilliamsAZ_Urban.las

NAU's email limits attachments to 20mb (What?!??!?!) so if you need anything specific or if this doens't meet someone needs please let me know and I can find something else to provide.

[Jean-Romain]

Thank you @bi0m3trics this should be good enough to run tests and debug stuff someday.

The readme of Limberger & Oliveira states:

-> Compile the application using Release x64 configuration On Visual Studio, simply open the project file "3dkht.sln"

There is nothing for other compilers. I guess that if you compiled the source you have VS. I don't have VS and I don't have a computer with VS. Do you think you can try to convert the VS file project into a makefile for gcc. See this stackoverflow thread. It seems that make it so can do this job pretty easily. But it must be run on a Windows machine with VS.

[bi0m3trics]

Thanks for the links and sorry for the lapse in memory with respect to windows. I'll look into it and let you know how it goes...

[Jean-Romain]

I compiled from this repo https://github.com/Kleiber/3dkht . Very hard, with missing files and wrong hard coded paths. But anyway this is what I get. Is is how it looks like for you? Very hard to understand what it does. And no way to tune the parameters but modifying the source code.

Edit: Hoo I understood. Everything is hardcoded in settings.h. :sob:

Loading Point Cloud...          
File: /home/jr/Téléchargements/3dkht_datasets/pointclouds/Museum.txt"
Size: 179744
Loaded in: 0.21 seconds.

3-D KHT starting...
Subdividing...                  Done [0.033025s]
Voting...                       Done [0.010086s]
Detecting Peaks....             Done [0.032277s]
Hough Transform:                Total [0.149]   6.71141 fps
3-D KHT ended: 11 planes found

[bi0m3trics]

Yeah... it's awful but the version I have compiled classifies the points corresponding to the planes, so your top figure (the museum) looks like the following for me...

It looks like you have it running correctly, and here's what it looked like when I fed it a real lidar dataset.

It's been sometime since I looked at the code, but based on the paper, pdal, and the code I recall thinking that the heavy lifting was done in the void octree_t::subdivide( hough_settings &settings ) function of the octree_t.cpp file which is where the coplanarity test appears to be preformed for each point.

[Jean-Romain]

Your output is much better than mine. In the mean time the source code I used, even if it is not the one form the authors looks identical. We can continue this discussion by email. This is not a good place discuss about third party code exploration/understanding/...

[Jean-Romain]

A very rough and inefficient R implementation of PDAL approximate coplanar filter. I'd like to see PDAL output for the same dataset.

library(lidR)
las = readLAS("building_WilliamsAZ_Urban.las")
dtm = grid_terrain(las, 0.5, tin())
las = lasnormalize(las, dtm)
las = lasfilter(las, Classification != 2L)
plot(las)

k   = 10
XYZ = as.matrix(lidR:::coordinates3D(las))
knn = RANN::nn2(XYZ, XYZ, k = k)$nn.idx

l1 = l2 = l3 = numeric(nrow(XYZ))
for (i in 1:nrow(XYZ)) {
  nn = knn[i,]
  xyz = XYZ[nn,]
  eigen = rev(prcomp(xyz)[[1]])
  l1[i] = eigen[1]
  l2[i] = eigen[2]
  l3[i] = eigen[3]
}

th1 = 6
th2 = 6
las@data$Coplanar = l2 > th1*l1 & th2*l2 > l3
plot(las, color = Coplanar)

[bi0m3trics]

I’ll do that as soon as I’m back from a visit to the vet, for both datasets. Maybe back in an hour or so…

[bi0m3trics]

So that took way longer than I remembered, but the pdal pipeline structure is a bit of a pain and not at all intuitive.... So in pdal, if I set both thresholds to 6, I got horrible classifications. So the above image is with the defaults of thresh1 = 25 and thresh2 = 6, and with knn = 10. The resulting classified file is here: WilliamsAZ_Urban_pdal_Classified.las

[Jean-Romain]

@bi0m3trics It did not provide the same output with both threshold set to 6 because prcomp does not return the eigenvalues of the covariance matrix but the square roots of the eigenvalues. Below a faster implementation with the actual eigenvalues. The good new is that it provides almost the same output than PDAL. Not worst at least.

@jacksonvoelkel I can make an efficient and documented coplanar detection function by the end of the week.

@bi0m3trics why PDAL does not classify the ground as coplanar?

Two interesting points:

• The coplanar function works better with non-normalized data. Not really surprising btw.
• Once we have the eigenvalues of the covariance matrix we can generalize the conditional statement to detect powerline (see bellow)

Rcpp::sourceCpp(code = "
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

// [[Rcpp::export]]
SEXP eigen_value(arma::mat A) {  
arma::mat coeff;
arma::mat score;
arma::vec latent;
arma::princomp(coeff, score, latent, A);
return(Rcpp::wrap(latent));
}")

library(lidR)
las = readLAS("WilliamsAZ_Urban.las")
las = lasfilter(las, Classification != 2L)
plot(las)

k  = 8
XYZ = as.matrix(lidR:::coordinates3D(las))
knn = RANN::nn2(XYZ, XYZ, k = k)$nn.idx

l1 = l2 = l3 = numeric(nrow(XYZ))
for (i in 1:nrow(XYZ))
{
  nn = knn[i,]
  xyz = XYZ[nn,]
  eigen = eigen_value(xyz)
  l1[i] = eigen[3,1]
  l2[i] = eigen[2,1]
  l3[i] = eigen[1,1]
}

th1 = 30
th2 = 5
th3 = 10
las@data$Coplanar = l2 > th1*l1 & th2*l2 > l3
las@data$Colinear = th3*l2 < l3 & th3*l1 < l3
plot(las, color = Coplanar)
plot(las, color = Colinear)

[Erikpinter]

Hello, I am also very interested in seeing the implementation of building detection in lidR!

I use the lidR package maily with mobile LiDAR data as opposed to airborne and am mostly getting good results, although there are some limitations in lidR due to the nature of my data.

Would it be possible, to implement an algorithm in a way that it would also deliver good results in detecting vertical surfaces from mobile LiDAR data? I could provide test-datasets and would be happy to help in testing.

Regards, Erik

[Jean-Romain]

@Erikpinter vertical surfaces should be detected as well. However the issue is more related to computation time. All my algorithms rely on QuadTrees to find the neighborhood of a given point in a fast way. This is relevant for airborne lidar data because the point cloud is spread on XY mainly. The Z coordinate account for ~0% of the variation on a wide tile. But with T-lidar to be efficient the algorithm should rely on an Octree.

[bi0m3trics]

@Jean-Romain pdal DOES classify the ground as coplanar, but in my old pipeline that would change the classification from 2 to 6 for the ground points, which was undesirable, so I remove them prior to the call to filter.approximatecoplanar. I didn't bother tweaking that for the results I sent...

FYI - the code above throws a "Error in eigen(xyz) : non-square matrix in 'eigen'" in the call to eigen(xyz) in both linux and windows for me...

[Jean-Romain]

It is just a clash with object names. I updated the code above

[bi0m3trics]

Should've seen that. Maybe I haven't had enough coffee yet... ;)

[Jean-Romain]

I pushed (b2e229c) a C++ function lascoplanar in the lidRplugins package. It is documented but suggest to modifications. I plan to move it into lidR later. I must evaluate pros and cons about adding a dependency to RcppArmadillo or not. I could rewrite the eigenvalue computation to do not depend on RcppArmadillo. I will also implement #205 for lascoplanar.

I made my job. I let you guys make tests :wink:

library(lidR)

las = readLAS("WilliamsAZ_Urban.las")
las = lasfilter(las, Classification != 2L)
las = lidRplugins::lascoplanar(las)
las = lidRplugins::lascoplanar(las, k = 10, th1 = 7, th2 = 0)

plot(las, color = Coplanar)
plot(las, color = Colinear)

[bi0m3trics]

I had run this prior to your push... but it looks great so far! I wonder what hulls will look like once applied to the coplanar points after a neighborhood cleaning filter (to remove the "coplanar noise") has been applied... ;) Lucklily, you've created an environment where this can be done my anyone! Thanks for doing this... you may have officially removed my need for any commerical software with this last addition. $!

[Jean-Romain]

The next step is clustering the output. I think it is what is done in Limberger & Oliveira but it is a much harder step. Coplanar noise filter seems easier to me.

[Jean-Romain]

I put that here as a proof of concept for later usage but you can't use it yet it only works on my local version.

library(lidR)
las = readLAS("WilliamsAZ_Urban.las", select = "xyzc")
col = c("gray", "gray", "blue", "darkgreen", "darkgreen", "darkgreen", "red", "gray", "cyan", "darkgray", "gray", "pink", "pink", "purple", "pink")

# Erase old classification but ground points.
las@data[Classification != 2L, Classification := 3L]

# Perform approximate coplanar test on non ground points
las = lascoplanar(las, filter = ~Classification != 2L)

# Reclassify coplanar point as building
las@data[Coplanar == TRUE, Classification := 6L]

# Perform approximate colinear test on non ground point
las = lascoplanar(las, k = 10, th1 = 10, th2 = 0, filter = ~Classification != 2L)

# Reclassify colinear points as wire if not already classified as building
las@data[Colinear == TRUE & Classification != 6L, Classification := 14L]

# plot
plot(las, color = Classification, colorPalette = col[3:15])

Not as good as the original classification but not too bad. @bi0m3trics do you know how the original classification was made? Manually?

To compare with original classification

library(lidR)
las = readLAS("WilliamsAZ_Urban.las", select = "xyzc")
col = c("gray", "gray", "blue", "darkgreen", "darkgreen", "darkgreen", "red", "gray", "cyan", "darkgray", "gray", "pink", "pink", "purple", "pink")
plot(las, color = Classification, colorPalette = col[1:15])

[bi0m3trics]

I did it in GlobalMapper v20 using their classification tools (~95%) in tge lidar module and then I cleaned it up manually where there were misclassification errors within the large planes or linear features.

[Jean-Romain]

Branch devel have two functions lascoplanar and lascolinear (parallelized at c++ level). I'm not sure if I will keep them "as is". @komazsofi suggested in #217 other eigen values related metrics that are, in my opinion simpler. For example lascoplanar based on Limberger, F. A., & Oliveira expects two threshold values (th) to assess coplanarity:

λ2 > (th1×λ1) && (th2×λ2) > λ3

But @komazsofi metrics are simpler:

(λ2 - λ3) / λ1 < th

Any opinion/advice on the question?

[komazsofi]

I think eigenvalue metrics can be used more widely it is not restricted to only extract the buildings. I have used it because it has shown that using these metrics it is possible to classify, trees, buildings etc. from ALS data using machine learning (e.g. his paper). But I am not familiar with the other approach which was suggested.

[Jean-Romain]

I finally opt for a scalable function that could accept more geometries or more methods to estimate local shapes.

In addition stdshapemetrics proposed by @komazsofi in #217 is an option do make a coarse estimation at the pixel level.

las = readLAS("~/Documents/ALS data/Fichiers las divers/building_WilliamsAZ_Urban_normalized.las")

# Find all plane
las = lasdetectshape(las, "plane", th = c(25, 6), k = 8)
plot(las, color = "Coplanar", colorPalette = c("blue", "red"))

# Find all plane but do no process ground points (faster)
las = lasdetectshape(las, "plane", th = c(25, 6), k = 8, filter = ~Classification != 2L)
plot(las, color = "Coplanar", colorPalette = c("blue", "red"))

# Find all line but do not process ground points 
las = lasdetectshape(las, "line", th = 10, k = 8, filter = ~Classification != 2L)
plot(las, color = "Colinear", colorPalette = c("blue", "red"))

[WalidGharianiEAGLE]

@Jean-Romain Ive been trying to apply the same recommended workflow for potential Trees/Buildings Segmentation using the shape_detection function. I got the following output where the fig in the left is from the normalized data without processing ground points (Here much of the buildings point clouds are filtered although I did beforehand a pmf, Im not sure why) and on the right is the output from the segmentation using shp_plane. Clearly something went wrong in this workflow to accurate filtering and segmentation. Here I provided the script and the data. https://github.com/WalidGharianiEAGLE/MB3_lidar/blob/lidR/lidR_processingR.R .I would appreciate any suggestion, Thanks for the great package!

[Jean-Romain]

Tested with a 10th of the point filter = "-keep_random_fraction 0.1"

1. The roof are classified as ground. You can't expect anything after that. Try csf(rigidness = 3, cloth_resolution = 1) it is much faster much better in you study case (still imperfect).
2. You normalized the point cloud. The consequence of such step is that original geometries are deformed. Not a lot here because the terrain is relatively flat but maybe enough to not be enable to detect planes or lines
3. Ask questions on https://gis.stackexchange.com/questions/tagged/lidr with the tag lidr

[WalidGharianiEAGLE]

@Jean-Romain Thank you for the great feedback! I did experimented with 10th of the point clouds. As expected It was much efficient for computation. Using the Cloth Simulation Filter for segmentation (fig in the left) worked much faster and even much accurate than the pmf (fig in the right), Interesting cause as far as I know the pmf is more robust and here the opposite is true!