nevrome - 1 year ago 73
R Question

# Filling a 3D-Body with a systematic point raster

I have a set of 3D-Bodies. Each Body is defined by 8 points with three coordinates each. All of the bodies are cubical or approximately cubical. I would like to "fill" the cubes with a systematic point raster. The coordinates are stored in simple data.frames.

I developed the following code that does what I want for cubical bodies:

``````# libraries
library(rgl)

# define example cube with 8 points
excube <- data.frame(
x = c(1,1,1,1,5,5,5,5),
y = c(1,1,4,4,1,1,4,4),
z = c(4,8,4,8,4,8,4,8)
)

# cubeconst: fill cube (defined by 8 corner points) with a 3D-point-raster
cubeconst <- function(x, y, z, res) {
cube <- data.frame()
xvec = seq(min(x), max(x), res)
yvec = seq(min(y), max(y), res)
zvec = seq(min(z), max(z), res)
for (xpoint in 1:length(xvec)) {
for (ypoint in 1:length(yvec)) {
for (zpoint in 1:length(zvec)) {
cube <- rbind(cube, c(xvec[xpoint], yvec[ypoint], zvec[zpoint]))
}
}
}
colnames(cube) <- c("x", "y", "z")
return(cube)
}

# apply cubeconst to excube
fcube <- cubeconst(x = excube\$x, y = excube\$y, z = excube\$z, res = 0.5)

# plot result
plot3d(
fcube\$x,
fcube\$y,
fcube\$z,
type = "p",
xlab = "x",
ylab = "y",
zlab = "z"
)
``````

Now I'm searching for a solution to "fill" approximately cubical bodies like for example the following body:

``````# badcube
x = c(1,1,1,1,5,5,5,5),
y = c(1,1,4,4,1,1,4,4),
z = c(4,10,4,12,4,8,4,8)
)

plot3d(
col = "red",
size = 10,
type = "p",
xlab = "x",
ylab = "y",
zlab = "z"
)
``````

Maybe you can point me in the right direction.

You need to transform the hexahedron (wonky cube) to a unit cube. The following image shows what I mean, and gives us a numbering scheme for the vertices of the hexa. Vertex `2` is hidden behind the cube.

The transformation is from real space `x,y,z`, to a new coordinate system `u,v,w`, in which the hexa is a unit cube. The typical function used for hexa looks like this.

``````x = A + B*u + C*v + D*w + E*u*v + F*u*w + G*v*w + H*u*v*w
``````

Transformations for y and z coordinates are of the same form. You have 8 corners to your cube, so you can substitute these in to solve for coefficients `A,B,...`. The unit coordinates `u,v,w` are either `0` or `1` at every vertex, so this simplifies things a lot.

``````x0 = A                             // everything = zero
x1 = A + B                         // u = 1, others = zero
x2 = A + C                         // v = 1, ...
x4 = A + D                         // w = 1
x3 = A + B + C + E                 // u = v = 1
x5 = A + B + D + F                 // u = w = 1
x6 = A + C + D + G                 // v = w = 1
x7 = A + B + C + D + E + F + G + H // everything = 1
``````

You then have to solve for `A,B,...`. This is easy because you just forward substitute. `A` equals `x0`. `B` equals `x1 - A`, etc... You have to do this for `y` and `z` also, but if your language supports vector operations, this can probably be done in the same step as for `x`.

Once you have the coefficients, you can convert a point `u,v,w` to `x,y,z`. Now, if you have a point generation scheme which works on a 1x1x1 cube, you can transform the result to the origonal hex. You could retain the same triple-loop structure in your posted code, and vary `u,v,w` between `0` and `1` to create a grid of points within the hex.

I'm afraid I don't know `r`, so I can't give you any example code in that language. Here's a quick `python3` example, though, just to prove it works.

``````import matplotlib.pyplot as pp
import numpy as np
from mpl_toolkits.mplot3d import Axes3D

np.random.seed(0)

cube = np.array([
[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 1.0, 0.0],
[0.0, 0.0, 1.0], [1.0, 0.0, 1.0], [0.0, 1.0, 1.0], [1.0, 1.0, 1.0]])

hexa = cube + 0.5*np.random.random(cube.shape)

edges = np.array([
[0, 1], [0, 2], [1, 3], [2, 3],
[0, 4], [1, 5], [2, 6], [3, 7],
[4, 5], [4, 6], [5, 7], [6, 7]])

def cubeToHexa(hexa, u, v, w):
A = hexa[0]
B = hexa[1] - A
C = hexa[2] - A
D = hexa[4] - A
E = hexa[3] - A - B - C
F = hexa[5] - A - B - D
G = hexa[6] - A - C - D
H = hexa[7] - A - B - C - D - E - F - G
xyz = (
A +
B*u[...,np.newaxis] +
C*v[...,np.newaxis] +
D*w[...,np.newaxis] +
E*u[...,np.newaxis]*v[...,np.newaxis] +
F*u[...,np.newaxis]*w[...,np.newaxis] +
G*v[...,np.newaxis]*w[...,np.newaxis] +
H*u[...,np.newaxis]*v[...,np.newaxis]*w[...,np.newaxis])
return xyz[...,0], xyz[...,1], xyz[...,2]

fg = pp.figure()

temp = np.reshape(np.append(hexa[edges], np.nan*np.ones((12,1,3)), axis=1), (36,3))
ax.plot(temp[:,0], temp[:,1], temp[:,2], 'o-')

u, v, w = np.meshgrid(*[np.linspace(0, 1, 6)]*3)
x, y, z = cubeToHexa(hexa, u, v, w)
ax.plot(x.flatten(), y.flatten(), z.flatten(), 'o')

pp.show()
``````

I can't recall the exact justificiation for this form of the transformation. It's certainly easy to solve, and it has no squared terms, so lines in the directions of the `u,v,w` axes map to straight lines in `x,y,z`. This means your cube edges and faces are guaranteed to conform, as well as the corners. I lack the maths to prove it, though, and I couldn't find any googleable information either. My knowledge comes from a distant memory of textbooks on Finite Element Methods, where these sort of transformations are common. If you need more information, I suggest you start looking there.

Recommended from our users: Dynamic Network Monitoring from WhatsUp Gold from IPSwitch. Free Download