---
title: "Clipping an isosurface to a ball, and more"
author: "Stéphane Laurent"
date: '2022-01-14'
tags: geometry, R, rgl, misc3d, graphics, python, pyvista
rbloggers: yes
output:
  md_document:
    variant: markdown
    preserve_yaml: true
  html_document:
    highlight: kate
    keep_md: no
highlighter: pandoc-solarized
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

We will firstly show how to clip an isosurface to a ball with R, and then, more 
generally, how to clip a surface to an arbitrary region. In the last part we 
show how to achieve the same with Python.

# Using R

## The Togliatti isosurface, clipped to a box

The Togliatti surface is the isosurface $f(x, y, z) = 0$, where the function 
$f$ is defined as follows in R:

```{r fTogliatti}
f <- function(x,y,z){
  64*(x-1)*
    (x^4-4*x^3-10*x^2*y^2-4*x^2+16*x-20*x*y^2+5*y^4+16-20*y^2) - 
    5*sqrt(5-sqrt(5))*(2*z-sqrt(5-sqrt(5)))*(4*(x^2+y^2-z^2)+(1+3*sqrt(5)))^2
}
```

To plot an isosurface in R, there is the **misc3d** package. Below we plot the 
Togliatti surface clipped to the box $[-5,5] \times [-5,5] \times [-4,4]$.

```{r Togliatti_box, eval=FALSE}
library(misc3d)
# make grid
nx <- 220L; ny <- 220L; nz <- 220L
x <- seq(-5, 5, length.out = nx) 
y <- seq(-5, 5, length.out = ny)
z <- seq(-4, 4, length.out = nz) 
G <- expand.grid(x = x, y = y, z = z)
# calculate voxel
voxel <- array(with(G, f(x, y, z)), dim = c(nx, ny, nz))
# compute isosurface
surf1 <- 
  computeContour3d(voxel, maxvol = max(voxel), level = 0, x = x, y = y, z = z)
# make triangles
triangles <- makeTriangles(surf1, smooth = TRUE)
```

```{r drawTogliatti_box, eval=FALSE}
drawScene.rgl(triangles, color = "yellowgreen")
```

![](figures/Togliatti_box.gif)

It is fun, but you will see later that it is prettier when clipped to a ball. 
And it is desirable to get rid of the isolated components at the top.

The simplest way to clip the surface to a ball consists in using the `mask` 
argument of the function `computeContour3d`. We use $4.8$ as the radius.

```{r drawTogliatti_mask, eval=FALSE}
mask <- array(with(G, x^2+y^2+z^2 < 4.8^2), dim = c(nx, ny, nz))
surf2 <- computeContour3d(
  voxel, maxvol = max(voxel), level = 0, mask = mask, x = x, y = y, z = z
)
triangles <- makeTriangles(surf2, smooth = TRUE)
drawScene.rgl(triangles, color = "orange")
```

![](figures/Togliatti_mask.gif)

That's not bad, but you see the problem: the borders are not smooth.


## Resorting to spherical coordinates

Using spherical coordinates will allow us to get the surface clipped to the 
ball with smooth borders. Here is the method:

```{r drawTogliatti_spherical, eval=FALSE}
# Togliatti surface equation with spherical coordinates
h <- function(ρ, θ, ϕ){
  x <- ρ * cos(θ) * sin(ϕ)
  y <- ρ * sin(θ) * sin(ϕ)
  z <- ρ * cos(ϕ)
  f(x, y, z)
}
# make grid
nρ <- 300L; nθ <- 300L; nϕ <- 300L
ρ <- seq(0, 4.8, length.out = nρ) # ρ runs from 0 to the desired radius
θ <- seq(0, 2*pi, length.out = nθ)
ϕ <- seq(0, pi, length.out = nϕ) 
G <- expand.grid(ρ=ρ, θ=θ, ϕ=ϕ)
# calculate voxel
voxel <- array(with(G, h(ρ, θ, ϕ)), dim = c(nρ, nθ, nϕ))
# calculate isosurface
surf3 <- 
  computeContour3d(voxel, maxvol = max(voxel), level = 0, x = ρ, y = θ, z = ϕ)
# transform to Cartesian coordinates
surf3 <- t(apply(surf3, 1L, function(ρθϕ){
  ρ <- ρθϕ[1L]; θ <- ρθϕ[2L]; ϕ <- ρθϕ[3L] 
  c(
    ρ * cos(θ) * sin(ϕ),
    ρ * sin(θ) * sin(ϕ),
    ρ * cos(ϕ)
  )
}))
# draw isosurface
drawScene.rgl(makeTriangles(surf3, smooth=TRUE), color = "violetred")
```

![](figures/Togliatti_spherical.gif)

Now the surface is pretty nice. The borders are smooth.


## Another way: using the rgl package 

Nowadays, in the **rgl** package, there is the function `clipMesh3d` which 
allows to *clip a mesh to a general region*. In order to use it, we need a 
`mesh3d` object. There is an unexported function in **misc3d** which allows to 
easily get a `mesh3d` object; it is called `t2ve`.

We start with the isosurface clipped to the box:

```{r rglmesh, eval=FALSE}
triangles <- makeTriangles(surf1)
# convert to rgl `mesh3d`
mesh <- misc3d:::t2ve(triangles)
rglmesh <- tmesh3d(mesh$vb, mesh$ib)
```

To define the region to which we want to clip the mesh, one has to introduce a 
function and to give a bound for this function. Here the region is 
$x^2 + y^2 + z^2 < 4.8^2$, so we proceed as follows:

```{r rglclippedmesh, eval=FALSE}
fn <- function(x, y, z) x^2 + y^2 + z^2 # or function(xyz) rowSums(xyz^2)
clippedmesh <- addNormals(clipMesh3d(
  rglmesh, fn, bound = 4.8^2, greater = FALSE
))
```

It is a bit slow. Probably the algorithm is not very easy. But we get our 
desired result:

```{r drawClippedMesh, eval=FALSE}
shade3d(clippedmesh, color = "magenta")
```

![](figures/Togliatti_clipped.gif)


# Using Python 

With the wonderful Python library **PyVista**, we proceed as follows to create 
the mesh of the Togliatti isosurface:

```python
from math import sqrt
import numpy as np
import pyvista as pv

def f(x, y, z):
    return (
        64
        * (x - 1)
        * (
            x ** 4
            - 4 * x ** 3
            - 10 * x ** 2 * y ** 2
            - 4 * x ** 2
            + 16 * x
            - 20 * x * y ** 2
            + 5 * y ** 4
            + 16
            - 20 * y ** 2
        )
        - 5
        * sqrt(5 - sqrt(5))
        * (2 * z - sqrt(5 - sqrt(5)))
        * (4 * (x ** 2 + y ** 2 - z ** 2) + (1 + 3 * sqrt(5))) ** 2
    )

# generate data grid for computing the values
X, Y, Z = np.mgrid[(-5):5:250j, (-5):5:250j, (-4):4:250j]
# create a structured grid
grid = pv.StructuredGrid(X, Y, Z)
# compute and assign the values
values = f(X, Y, Z)
grid.point_data["values"] = values.ravel(order="F")
# compute the isosurface f(x, y, z) = 0
isosurf = grid.contour(isosurfaces=[0])
mesh = isosurf.extract_geometry()
```

To plot the mesh clipped to the box:

```python
# surface clipped to the box:
mesh.plot(smooth_shading=True, color="yellowgreen", specular=15)
```

The `remove_points` method is similar to the `mask` approach with **misc3d** 
(it produces non-smooth borders):

```python
# surface clipped to the ball of radius 4.8, with the help of `remove_points`:
lengths = np.linalg.norm(mesh.points, axis=1)
toremove = lengths >= 4.8
masked_mesh, idx = mesh.remove_points(toremove)
masked_mesh.plot(smooth_shading=True, color="orange", specular=15)
```

To get smooth borders, you have to use the `clip_scalar` method:

```python
# surface clipped to the ball of radius 4.8, with the help of `clip_scalar`:
mesh["dist"] = lengths
clipped_mesh = mesh.clip_scalar("dist", value=4.8)
clipped_mesh.plot(
    smooth_shading=True, cmap="inferno", window_size=[512, 512],
    show_scalar_bar=False, specular=15, show_axes=False, zoom=1.2,
    background="#363940"
)
```

![](figures/Togliatti_python.png)