Plotly sphere intensity coloring with intersection

Question

How can I color a sphere in plotly relative to its hover value, that value's position within a range of 0-500, and a heat-map palette?

Further, how can I color the intersection of two spheres based on the sum of each sphere's hover value?

E.g. if sphere 1 has hover value 150, and sphere 2 has 100, the intersection should have hover value 250 and a darker color.

library(plotly)

f <- function(x, y, z){
  x^2 + y^2 + z^2
}
R <- 2 # radius
x <- y <- z <- seq(-R, R, length.out = 100) 
g <- expand.grid(x = x, y = y, z = z)
voxel <- array(with(g, f(x, y, z)), dim = c(100, 100, 100))

library(misc3d)
cont <- computeContour3d(voxel, level = R^2, x = x, y = y, z = z)
idx <- matrix(0:(nrow(cont)-1), ncol=3, byrow=TRUE)

plot_ly() |> 
  add_trace(
    type = "mesh3d",
    x = cont[, 1], y = cont[, 2], z = cont[, 3],
    i = idx[, 1], j = idx[, 2], k = idx[, 3],
    opacity = 0.2,
    hoverinfo = "text",
    text = "150"
  ) |> 
  add_trace(
    type = "mesh3d",
    x = cont[, 1]+2, y = cont[, 2]+2, z = cont[, 3]+2,
    i = idx[, 1], j = idx[, 2], k = idx[, 3],
    opacity = 0.2,
    hoverinfo = "text",
    text = "43"
  ) 

Answer 1

First I would recommend the rmarchingcubes package instead of misc3d because it is more efficient.

library(plotly)
library(rmarchingcubes)

f <- function(x, y, z){
  x^2 + y^2 + z^2
}
R <- 2 # radius

n <- 100
x_ <- y_ <- z_ <- seq(-R, R, length.out = 100) 
Grid <- expand.grid(X = x_, Y = y_, Z = z_)
voxel <- with(Grid, array(f(X, Y, Z), dim = c(n, n, n)))
surf <- contour3d(voxel, level = R^2, x_, y_, z_)

vertices <- surf$vertices
indices  <- surf$triangles - 1

plot_ly(
  x = vertices[, 1], y = vertices[, 2], z = vertices[, 3],
  i = indices[, 1], j = indices[, 2], k = indices[, 3],
  type = "mesh3d"
) %>% layout(scene = list(aspectmode = "data"))

I don't understand what you mean by "the position of the sphere within a range of 0-500".

If you want to color the sphere according to some values attached to its points, you firstly need to define the function which maps the points to these values. For example, I will take the sine of the azimuthal angle, to get something periodic. Then I firstly compute these values:

# azimuthal angles
phi <- atan2(vertices[, 2], vertices[, 1])
# we color a vertex according to sin(phi)
values <- sin(phi)

These values range from -1 to 1. The colorRamp function is very nice to associate a color to a number; it requires numbers between 0 and 1, so I linearly map the values to this range by applying the function v -> (v + 1)/2:

# we map sin(phi) to (0, 1) in order to use the colorRamp function
u <- (values + 1) / 2
#
fPalette <- colorRamp(c("yellow", "green", "blue"))
vertexColors <- fPalette(u)

At this step, vertexColors is a matrix representing RGB values between 0 and 255, one has to convert them to hex coded colors:

vertexColors <- rgb(
  vertexColors[, 1], vertexColors[, 2], vertexColors[, 3], maxColorValue = 255
)

And then we can do the plot with these colors:

plot_ly(
  x = vertices[, 1], y = vertices[, 2], z = vertices[, 3],
  i = indices[, 1], j = indices[, 2], k = indices[, 3],
  vertexcolor = toRGB(vertexColors),
  type = "mesh3d"
) %>% layout(scene = list(aspectmode = "data"))

Regarding the intersection, this is not very clear to me. You can get an implicit equation g(x,y,z)=0 of the intersection of two isosurfaces f1(x,y,z)=0 and f2(x,y,z)=0 by taking g(x,y,z) = max(f1(x,y,z), f2(x,y,z)).

f1 <- function(x, y, z){
  x^2 + y^2 + z^2 - R^2
}
f2 <- function(x, y, z){
  (x-2)^2 + (y-2)^2 + (z-2)^2 - R^2
}
g <- function(x, y, z) {
  pmax(f1(x, y, z), f2(x, y, z))
}
voxel <- with(Grid, array(g(X, Y, Z), dim = c(n, n, n)))
surf <- contour3d(voxel, level = 0, x_, y_, z_)
vertices <- surf$vertices
indices  <- surf$triangles - 1
plot_ly(
  x = vertices[, 1], y = vertices[, 2], z = vertices[, 3],
  i = indices[, 1], j = indices[, 2], k = indices[, 3],
  type = "mesh3d"
) %>% layout(scene = list(aspectmode = "data"))

But if you want to color each vertex of this intersection according to the two spheres values of this vertex, you need a "sphere value" associated to each point in the volume bounded by the sphere. Having some values only for the points at the surface is not enough, because this is the intersection of the volumes.

---

EDIT

The idea you advanced in a comment is a good one: separately plot the first sphere minus the intersection, the second sphere minus the intersection, and the intersection.

This is what the code below does. To compute a mesh of a sphere minus the intersection, I use the clipMesh3d function of the rgl package. I also use it to compute a mesh of the intersection.

I use the facecolor argument to set the color of a whole mesh: I assign the same color to all faces of the mesh. The documentation claims that one can use the color argument to set the color of a whole mesh but that does not work.

library(plotly)
library(rmarchingcubes)
library(rgl)

R <- 2
f1 <- function(x, y, z){
  x^2 + y^2 + z^2 - R^2
}
f2 <- function(x, y, z){
  (x-2)^2 + (y-2)^2 + (z-2)^2 - R^2
}

n <- 100
x_ <- y_ <- z_ <- seq(-R, R, length.out = n) 
Grid <- expand.grid(X = x_, Y = y_, Z = z_)
voxel <- with(Grid, array(f1(X, Y, Z), dim = c(n, n, n)))
surf <- contour3d(voxel, level = 0, x_, y_, z_)
mesh_sphere1 <- tmesh3d(
  vertices = t(surf$vertices), indices = t(surf$triangles), homogeneous =  FALSE
)
mesh1 <- clipMesh3d(mesh_sphere1, fn = f2, bound = 0)
mesh_halfintersection1 <- 
  clipMesh3d(mesh_sphere1, fn = f2, bound = 0, greater = FALSE)

x_ <- y_ <- z_ <- seq(-R, R, length.out = n) + 2
Grid <- expand.grid(X = x_, Y = y_, Z = z_)
voxel <- with(Grid, array(f2(X, Y, Z), dim = c(n, n, n)))
surf <- contour3d(voxel, level = 0, x_, y_, z_)
mesh_sphere2 <- tmesh3d(
  vertices = t(surf$vertices), indices = t(surf$triangles), homogeneous =  FALSE
)
mesh2 <- clipMesh3d(mesh_sphere2, fn = f1, bound = 0)
mesh_halfintersection2 <- 
  clipMesh3d(mesh_sphere2, fn = f1, bound = 0, greater = FALSE)

mesh3 <- merge(mesh_halfintersection1, mesh_halfintersection2)

vertices1 <- t(mesh1$vb[1:3,])
vertices2 <- t(mesh2$vb[1:3,])
vertices3 <- t(mesh3$vb[1:3,])
indices1 <- t(mesh1$it) - 1L
indices2 <- t(mesh2$it) - 1L
indices3 <- t(mesh3$it) - 1L

hover1 <- 350
hover2 <- 100
hover3 <- hover1 + hover2
fPalette <- colorRamp(c("yellow", "green", "blue"))
colors <- fPalette(c(hover1, hover2, hover3) / 500)
colors <- rgb(
  colors[, 1], colors[, 2], colors[, 3], maxColorValue = 255
)

plot_ly() |>
  add_trace(
    type = "mesh3d",
    x = vertices1[, 1], y = vertices1[, 2], z = vertices1[, 3],
    i = indices1[, 1], j = indices1[, 2], k = indices1[, 3],
    opacity = 0.2, facecolor = toRGB(rep(colors[1], nrow(indices1)))
  ) |>
  add_trace(
    type = "mesh3d",
    x = vertices2[, 1], y = vertices2[, 2], z = vertices2[, 3],
    i = indices2[, 1], j = indices2[, 2], k = indices2[, 3],
    opacity = 0.2, facecolor = toRGB(rep(colors[2], nrow(indices2)))
  ) |>
  add_trace(
    type = "mesh3d",
    x = vertices3[, 1], y = vertices3[, 2], z = vertices3[, 3],
    i = indices3[, 1], j = indices3[, 2], k = indices3[, 3],
    opacity = 0.2, facecolor = toRGB(rep(colors[3], nrow(indices3)))
  ) |>
  layout(scene = list(aspectmode = "data"))