2016-03-20 19 views
2

我有一套3D物體。每個物體由8個點定義,每個點有三個座標。所有的身體都是立方體或近似立方體的。我想用系統化的點柵格「填充」立方體。座標存儲在簡單的data.frames中。用系統化的點光柵填充3D物體

我開發了以下代碼,做什麼,我想爲立方形機構:

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

現在我正在尋找一個解決方案,以「補」大約立方形機構例如像了以下機身:

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

也許你可以指點我正確的方向。

+0

請您詳細說明您的「系統點柵格」是什麼意思?我把你的問題看成是問如何將一個扭曲的立方體分成N×N×N個較小的立方體。它是否正確? – Bill

+0

@Bill是的,我想你明白了。算法的結果應該是扭曲立方體限制範圍內的點列表。這些點可以是隨機分佈的,但我更喜歡它們是等距離的。將NxNxN分成更小的立方體可能是實現這一目標的一種方法。但是如何? – nevrome

回答

2

您需要將六面體(wonky cube)轉換爲單位立方體。下圖顯示了我的意思,並給出了hexa頂點的編號方案。頂點2隱藏在立方體後面。

enter image description here

轉變是從實際空間x,y,z,到一個新的座標系u,v,w,其中,所述六是單位立方體。用於hexa的典型功能如下所示。

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

y和z座標變換的形式是相同的。你有8個角落到你的立方體,所以你可以用這些替代來解決係數A,B,...。單位座標u,v,w在每個頂點都是01,所以這可以簡化很多事情。

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 

然後您必須解決A,B,...。這很簡單,因爲你只需轉發替代品。 A等於x0B等於x1 - A等等...您也必須爲yz這樣做,但是如果您的語言支持矢量操作,則可能與x的步驟相同。

一旦你有係數,你可以將一個點u,v,w轉換爲x,y,z。現在,如果您有一個適用於1x1x1立方體的點生成方案,則可以將結果轉換爲原始十六進制。您可以在發佈的代碼中保留相同的三重循環結構,並在01之間改變u,v,w,以在十六進制內創建一個網格點。

恐怕我不知道r,所以我不能給你任何使用該語言的示例代碼。儘管這裏有一個快速的python3例子,但它只是爲了證明它的工作原理。

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

我不記得這種形式的轉換的確切理由。這當然很容易解決,並且它沒有平方項,所以u,v,w軸的方向上的線映射爲x,y,z中的直線。這意味着您的立方體邊和麪保證符合以及角落。儘管我缺乏數學證明,但我也找不到任何可用於Google的信息。我的知識來自對有限元方法的教科書的遙遠記憶,其中這些變換是常見的。如果您需要更多信息,我建議您開始尋找。

+0

太棒了。我花了兩個小時在R中重建這個非常棒的解決方案 - 但現在它可以工作。非常感謝!明天我會清理R函數並將其發佈到這裏。 – nevrome

1

由於票據的解釋和例子我能拿出R中的以下解決方案:

# libraries 
library(rgl) 

# create heavily distorted cube - hexahedron 
hexatest <- data.frame(
    x = c(0,1,0,4,5,5,5,5), 
    y = c(1,1,4,4,1,1,4,4), 
    z = c(4,8,4,9,4,8,4,6) 
) 

# cubetohexa: Fills hexahedrons with a systematic point raster 
cubetohexa <- function(hexa, res){ 

    # create new coordinate system (u,v,w) 
    resvec <- seq(0, 1, res) 
    lres <- length(resvec) 

    u <- c() 
    for (p1 in 1:lres) { 
    u2 <- c() 
    for (p2 in 1:lres) { 
     u2 <- c(u2, rep(resvec[p2], lres)) 
    } 
    u <- c(u,u2) 
    } 

    v <- c() 
    for (p1 in 1:lres) { 
    v <- c(v, rep(resvec[p1], lres^2)) 
    } 

    w <- rep(resvec, lres^2) 

    # transformation 
    A <- as.numeric(hexa[1,]) 
    B <- as.numeric(hexa[2,]) - A 
    C <- as.numeric(hexa[3,]) - A 
    D <- as.numeric(hexa[5,]) - A 
    E <- as.numeric(hexa[4,]) - A - B - C 
    F <- as.numeric(hexa[6,]) - A - B - D 
    G <- as.numeric(hexa[7,]) - A - C - D 
    H <- as.numeric(hexa[8,]) - A - B - C - D - E - F - G 

    A <- matrix(A, ncol = 3, nrow = lres^3, byrow = TRUE) 
    B <- matrix(B, ncol = 3, nrow = lres^3, byrow = TRUE) 
    C <- matrix(C, ncol = 3, nrow = lres^3, byrow = TRUE) 
    D <- matrix(D, ncol = 3, nrow = lres^3, byrow = TRUE) 
    E <- matrix(E, ncol = 3, nrow = lres^3, byrow = TRUE) 
    F <- matrix(F, ncol = 3, nrow = lres^3, byrow = TRUE) 
    G <- matrix(G, ncol = 3, nrow = lres^3, byrow = TRUE) 
    H <- matrix(H, ncol = 3, nrow = lres^3, byrow = TRUE) 

    for (i in 1:(lres^3)) { 
    B[i,] <- B[i,] * u[i] 
    C[i,] <- C[i,] * v[i] 
    D[i,] <- D[i,] * w[i] 
    E[i,] <- E[i,] * u[i] * v[i] 
    F[i,] <- F[i,] * u[i] * w[i] 
    G[i,] <- G[i,] * v[i] * w[i] 
    H[i,] <- H[i,] * u[i] * v[i] * w[i] 
    } 

    m <- data.frame(A+B+C+D+E+F+G+H) 
    colnames(m) <- c("x", "y", "z") 

    # output 
    return(m) 
} 

# apply cubetohexa to hexatest 
cx <- cubetohexa(hexatest, 0.1) 

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

編輯:

此功能與RCPP現在實現我的R包recexcavAAR

+0

很高興:) – Bill