nevrome - 4 months ago 22

R Question

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`

badcube <- 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,10,4,12,4,8,4,8)

)

# plot badcube

plot3d(

badcube$x,

badcube$y,

badcube$z,

col = "red",

size = 10,

type = "p",

xlab = "x",

ylab = "y",

zlab = "z"

)

Maybe you can point me in the right direction.

Answer

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()
ax = fg.add_subplot(111, projection='3d')
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.