3D voxel rendering

[Saadi4469]

Hello,

Is it possible to do 3D tree modeling using voxels like shown in Figure 1 (under section 2.3) of this paper using lidR? Currently, I tried the code below, and the results are good, but I was still wondering if it could be taken to the next level.

# Voxel metrics
TLS_LAS_vox_met = voxel_metrics(TLS_LAS_denoised, ~list(N = length(Z)), 0.01) # calculate voxel metrics
plot(TLS_LAS_vox_met , color="N", colorPalette = heat.colors(50), size = 1, bg = "white")

Thank you

Reference Hess, C., Härdtle, W., Kunz, M., Fichtner, A., & von Oheimb, G. (2018). A high‐resolution approach for the spatiotemporal analysis of forest canopy space using terrestrial laser scanning data. Ecology and evolution, 8(13), 6800-6811.

[Jean-Romain]

If you mean plotting cubes instead of points. No it is not possible with lidR because the rendering is excessively slow so I chose to display points. Now, if you put your hand into rgl, everything is doable

[Saadi4469]

You see in that figure the points connect and are rendered to show an actual 3D tree, so you mean it's possible to pass the voxel object to rgl and plot it the way its shown in that paper. If so then can you please share any good packages in R that have the ability to do so or should I just use rgl?

[Jean-Romain]

You must use "bare bone" rgl stuff.

• https://stackoverflow.com/questions/13982597/plotting-voxel-image-in-r
• https://stackoverflow.com/questions/35177603/r-visualize-many-3d-boxes-inside-a-big-box
• https://richbits.rbarnes.org/plotting-voxels-in-r.html

[Jean-Romain]

Finally, after many years of search, I discovered how to render many cubes more efficiently using shapelist3d(). For example the second example I gave you was:

library(rgl)
grd <- expand.grid(x=seq(0,10,2), y=seq(0,10,2), z=seq(0,10,2))
grd$dist <- sqrt(grd$x^2 + grd$y^2 + grd$z^2) # distance to coordinate 0,0,0
grd$col <- rainbow(ceiling(max(grd$dist+1)))[ceiling(grd$dist+1)]
grd$alpha <- rep(c(0.2, 1), each=nrow(grd)/2)
open3d()
for(i in seq(nrow(grd))){
  shade3d( translate3d( cube3d(col = grd$col[i]), grd$x[i], grd$y[i], grd$z[i]) , alpha=grd$alpha[i])
}

And can be written:

library(rgl)
grd <- expand.grid(x=seq(0,10,2), y=seq(0,10,2), z=seq(0,10,2))
grd$dist <- sqrt(grd$x^2 + grd$y^2 + grd$z^2) # distance to coordinate 0,0,0
grd$col <- rainbow(ceiling(max(grd$dist+1)))[ceiling(grd$dist+1)]
grd$alpha <- rep(c(0.2, 1), each=nrow(grd)/2)
n = nrow(grd)
u = vector("list", n)
for(i in seq(nrow(grd))){
  u[[i]] = translate3d( cube3d(col = grd$col[i], alpha = grd$alpha[i]), grd$x[i], grd$y[i], grd$z[i])
}

shapelist3d(u)

Which is 5 to 15 times faster. In this condition I guess I can try to add a voxel rendering in lidR.

[Jean-Romain]

I started doing something that looks good.

[Jean-Romain]

So I added a voxel rendering mode in v3.2.0 in the devel branch. You must install it that way

remotes::install_github("Jean-Romain/lidR", ref = "devel")

You can test it that way:

LASfile <- system.file("extdata", "MixedConifer.laz", package="lidR")
las <- readLAS(LASfile)
vm <- voxel_metrics(las, ~list(Intensity = mean(Intensity)), 2)
plot(vm, color = "Intensity", voxels = TRUE)

It is still experimental and your feedbacks are welcomed. Remember that it is much more computationally demanding and should be reserved to small scene with not too much voxels. My plan is to drastically improve the rendering by removing the voxels we can't see. There is no need to render voxels we can't see unless voxels are displayed transparent.

[Saadi4469]

Thank you for sharing this and keep up the good work.

Cheers

[Saadi4469]

After updating the package to 3.2.0 and using the find_trees function in manual mode, I can no longer drag the view with the left mouse button. Instead, it's again like I reported in another issue here, that I have to use the left mouse button to create a selection window, and even that window only removes a single tree top. Moreover, the middle mouse button stretches the view. And thus, the drag function is again missing.

remotes::install_github("Jean-Romain/lidR", ref = "devel")

[Jean-Romain]

My bad, I did not push the new code online in this branch. It should work now.

[Saadi4469]

Is it to possible to export the voxel object to ArcGIS Pro? The purpose would be to export the voxel object to ArcGIS Pro and render it there in order to overlay it on a base-map or raster as mentioned here. The results from 3D voxel rendering are really cool as seen here, but I would also like to overlay this 3D model over a base-map. Another option could be to embed this feature in lidR (with optional ability to add a north arrow, title, scale bar and etc to basically be able to save the 3D model in the form of a map as .PNG or other image formats).

Thank you

[Jean-Romain]

I've no idea if it is possible to export to ArcGIS. First we must know what kind of format are supported by ArcGIS. Putting a map in the rgl scene should be doable. See these resources:

• https://www.tylermw.com/adding-open-street-map-data-to-rayshader-maps-in-r/
• https://github.com/h-a-graham/rayvista

I already created publication ready figures with 3D DTM + point-cloud overlay like that:

Plotting voxels uses the same function than for point-cloud so it is not a big deal to display voxels on top a 3D map. It is only a matter of finding satellite images an applying it. Which is what the rayvista is doing.

[Saadi4469]

Thank you for sharing this.

[Jean-Romain]

I added an optimization to do not render hidden voxels. I'm closeing it is done for 3.2.0

[Saadi4469]

Thank you for adding this optimization. On a side note, can you please look into the comment I made here?

[Saadi4469]

Good day JRR, hope all is well, I found some more info on the 3D voxel format supported by ArcGIS Pro 3. When you have have the time can you kindly look at it? Maybe then you can suggest a way to bring the cool 3D stuff that lidR and bring it into ArcGIS Pro 3. This way the product can be shared via ESRI based web applications with other stakeholders and the general public.

What is a voxel layer?

Create a voxel layer (this is the method that ESRI uses maybe this can give some further insight too).

[bi0m3trics]

I apologize for enabling the side-tracking of this thread - Looks like there's already a package that reads/writes/uses netCDF file format (https://cran.r-project.org/web/packages/ncdf4/) in R using either version 3 or 4 of the netcdf library. Have you looked into this @Saadi4469? Not sure if it's the same one ESRI uses (would not be surprised if they tweaked it in some unconventional way) but you might start there to make sure your request isn't already implemented.

[Saadi4469]

@bi0m3trics thank you, yeah I know about netCDF4 but I don't know how to convert a 3D voxel object that lidR creates to a nc object then write that object to file.