Transform 3D Lidar-Points from cartesian coordinate-system to spherical coordinate-system and a stereographic projection

Question

I'm working on my master thesis to simulate hemispherical photographs from Lidar-Data. So my main goal is to project 3D Points (X,Y,Z) which are in a cartesian coordinate system to a stereographic projection (See picture 1, from: here). The coordinate-system of my Pointcloud is transformed so that the center point is located a (0,0,0) and all z-values are positive.

I'm coding in RStudio and I first tried to achieve a spherical projection of the Pointcloud by using the formula for cartesian to spherical coordinates listed on wikipedia.

r     <- sqrt(x^2 + y^2  + z^2)
theta <- acos(z/r)
phi   <- atan2(y,x)

I also tried to do it with the cart2sph-function from the RPackage pracma which should do the same. But unfortunately the results with both methods don't look nearly how I wanted them to be. The points don't fit on a plane half-sphere but just seem to be oriented along the z-axis (See picture 2).

Has anyone suggestions how to achieve the stereographic projection of the Lidar Pointcloud in RStudio? Is there maybe even a package with some functions to do a coordinate transformation? Unfortunately I don't have many coding experiences, I would be very thankful for every help.

Edit

I plot the cartesian points on a hemisphere using Allans script. The result looks like this using the software Cloudcompare:

Image: Hemisphere

I also did a stereographic projection with the transformed cartesian coordinates:

x2 <- (x / (1 + z))
y2 <- (y / (1 + z))
project_2d <- data.frame(x2, y2)

Image: Stereographic projection

Still I wonder why the outer edge of the plots don't fit a perfect circle since all tree trunks have the same min(z) value (I used a clipping box). So all tree trunks should line up on the "equator" of the hemisphere. Do you have any clue?

Hi Radde. It looks as though your points are all on the y-z plane (though it's difficult to tell with this particular image). Are you sure the transformed variables aren't correct on this plane? — Allan Cameron
hey, the points are oriented along the x-z plane and the x zero point splits the point cloud in two halves. But still I can't see any spherical-like transformation in my result. — Radde1683
If I understand, you project the scene onto a hemisphere, then the latter onto a plane by stereography. What are you showing us on the picture 2 (alias 4) ? — Yves Daoust
And why do you mention a spherical projection and show the formulas of conversion to spherical coordinates ? — Yves Daoust
Hey Yves, thanks for your answer. I want to project my local points from a carthesian coordinate system on a spherical half plane first and later on a stereographic 2D-Plane (like shown in picture 1). Don't I have to a conversion to spherical coordinates to do so? — Radde1683

Allan Cameron Allan Cameron · Accepted Answer · 2020-09-24T13:44:00

It's difficult to show this result in 3D, but I'll try my best.

Suppose we have a collection of points in the shape of a cross suspended over the camera, paralell to the ground:

x <- c(rep(seq(-10, 10, 1), 2), rep(-1, 21), rep(1, 21))
y <- c(rep(-1, 21), rep(1, 21), rep(seq(-10, 10, 1), 2))
z <- rep(3, 84)

We can try to show this in rgl, adding a red dot to indicate camera position:

library(rgl)

plot3d(x, y, z, zlim = c(0, 10))
points3d(0, 0, 0, col = "red")

Using the formula you found for transformation to spherical co-ordinates we can convert these 3D points from x, y, z to theta, phi and r. However, to project them onto a hemisphere in cartesian co-ordinates, we need have r at a fixed distance from the camera. We then need to calculate the x, y, z co-ordinates of the points on that hemisphere from the theta, phi, and fixed r. This function should do that:

projection <- function(x, y, z) {
  r     <- sqrt(x^2 + y^2  + z^2)
  theta <- acos(z / r)
  phi   <- atan2(y, x)
  y <- theta * cos(phi)
  x <- theta * sin(phi)
  z <- sqrt((pi / 2)^2 - x^2 - y^2)
  data.frame(x, y, z)
}

So now when we do:

df <- projection(x, y, z)

We can plot again to show the points projected onto a hemisphere:

plot3d(df$x, df$y, df$z, zlim = c(0, pi))

View from below

View from side

Transform 3D Lidar-Points from cartesian coordinate-system to spherical coordinate-system and a stereographic projection

1 Answers