CODEWITHSUNDEEP
rgeometrycoordinate-transformationlidar
Radde1683
Radde1683

Reputation: 37

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

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.

enter image description here

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).

enter image description here

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?

Upvotes: 1

Views: 1722

Answers (1)

Allan Cameron
Allan Cameron

Reputation: 173793

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")

enter image description here

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

enter image description here

View from side

enter image description here

Upvotes: 2

Related Questions