Search code examples
rmatlab3d

Plotting a 3D object in R


I'm trying to create a 3D shape plot, as is done in MATLAB using the patch command, like this one:

enter image description here

I've created my coordinate data using expand.grid:

library(pracma)
n1=60
n2=60
n3=60
data = expand.grid(x=linspace(-3,3,n1),y=linspace(-3,3,n2),z=linspace(-3,3,n3))

I've used this data to fill a 3D matrix to obtain a 3D matrix V, which defines the shape:

someData <- rep(0, 60*60*60)  
V <- array(someData, c(60, 60, 60))
count = 1
lambda = c(2,1,1)
for (ii in 1:n1){
    for (jj in 1:n2){
        for (kk in 1:n3){
            V[ii, jj, kk] = t(as.matrix(lambda))%*%as.matrix(c(data$x[count],data$y[count],data$z[count]))
            count = count+1
        }
    }
}

However, I'm struggling to plot this in R.

Reproducible MATLAB example:

n1 = 60;
n2 = 60;
n3 = 60;
[y,x,z] = ndgrid(linspace(-3,3,n1),linspace(-3,3,n2),linspace(-3,3,n3));
V = zeros(n1, n2, n3);
lambda = [2;1;1];
for ii = 1: n1
    for jj = 1:n2
        for kk =1:n3
            V(ii, jj, kk) =  lambda'*[abs(y(ii, jj, kk)); abs(x(ii, jj, kk)); abs(z(ii, jj, kk))];
        end
    end
end
p = patch(isosurface(x,y,z,V,1));
p.FaceColor = 'cyan';
p.EdgeColor = 'none';
view(3);
camlight
axis equal
lighting gouraud
box on

My attempt in R:

n1=60
n2=60
n3=60
someData <- rep(0, 60*60*60)  
V <- array(someData, c(60, 60, 60))
data = expand.grid(x=linspace(-3,3,n1),y=linspace(-3,3,n2),z=linspace(-3,3,n3))
lambda = c(2,1,1)
count = 1
for (ii in 1:n1){
    for (jj in 1:n2){
        for (kk in 1:n3){
            V[ii, jj, kk] =  t(as.matrix(lambda))%*%as.matrix(c(data$x[count],data$y[count],data$z[count]))
                 count = count+1
        }
    }
}
M = melt(V)
plot_ly(type = "mesh3d",
        x = data$x,
        y = data$y,
        z = data$z,
        i = M$Var1,
        j = M$Var2,
        k = M$Var3,
        facecolor = rep(toRGB(colorRampPalette(c("navy", "blue"))(6)), each = 2)
)

However, I just get an empty coordinate axis, with no shape


Solution

  • You can use the rmarchingcubes package to compute the isosurface and the rgl package to plot in 3D. The shape I get with the code below looks like a diamond. I didn't try to improve your code (there are possible improvements).

    library(pracma)
    library(rmarchingcubes) # to compute isosurface
    library(rgl) # to plot in 3D
    
    # your code with a couple of modifications
    n1 <- 60
    n2 <- 60
    n3 <- 60
    x <- linspace(-3, 3, n1)
    y <- linspace(-3, 3, n2)
    z <- linspace(-3, 3, n3)
    data <- expand.grid(x = x, y = y, z = z)
    someData <- rep(0, 60 * 60 * 60)
    V <- array(someData, c(60, 60, 60))
    count <- 1
    lambda <- c(2, 1, 1)
    for(ii in 1:n1) {
      for(jj in 1:n2) {
        for(kk in 1:n3) {
          V[ii, jj, kk] <- t(as.matrix(lambda)) %*% # you forgot the absolute value:
            abs(as.matrix(c(data$x[count], data$y[count], data$z[count])))
          count <- count + 1
        }
      }
    }
    
    # compute isosurface
    isosurface <- contour3d(
      griddata = V, level = 1, x = x, y = y, z = z
    )
    
    # make rgl mesh
    mesh <- tmesh3d(
      vertices = t(isosurface[["vertices"]]),
      indices = t(isosurface[["triangles"]]),
      normals = isosurface[["normals"]],
      homogeneous = FALSE
    )
    
    # plot
    shade3d(mesh, color = "turquoise")
    

    EDIT

    A big improvement to calculate V:

    Mat <- abs(t(as.matrix(data)))
    V <- array(t(lambda) %*% Mat, dim = c(n1, n2, n3))
    

    EDIT

    If you really want to use plotly:

    library(plotly)
    vs = isosurface[["vertices"]]
    idx = isosurface[["triangles"]] - 1
    plot_ly(
      x = vs[,1], y = vs[,2], z = vs[,3],
      i = idx[,1], j = idx[,2], k = idx[,3],
      type = "mesh3d"
    ) %>% layout(scene = list(aspectmode = "data"))