--- title: "Drawing slices of a hypersurface with R" author: "Stéphane Laurent" date: '2020-05-01' tags: R, graphics, rgl, geometry, maths 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, collapse = TRUE) ``` Let $s \colon I \times J \times K \to \mathbb{R}^4$ be a parameterization of a hypersurface $\mathcal{S}$, where $I,J,K \subset \mathbb{R}$ are some intervals. I'm going to show how to draw the cross-section of $\mathcal{S}$ by a hyperplane with R. For the illustration, we consider the [*tiger*](http://hi.gher.space/wiki/Tiger): ```{r, attr.source='.numberLines'} R1 = 2; R2 = 2; r = 0.5 s <- function(u, v, w){ rbind( cos(u) * (R1 + r*cos(w)), sin(u) * (R1 + r*cos(w)), cos(v) * (R2 + r*sin(w)), sin(v) * (R2 + r*sin(w)) ) } ``` Take a hyperplane: $$ \mathcal{P}\colon \quad \langle \mathbf{a}, \mathbf{x} \rangle = b, $$ let $\vec{\mathbf{n}} = \frac{\mathbf{a}}{\Vert\mathbf{a}\Vert}$ be a unit normal vector of $\mathcal{P}$, and $\mathbf{x}_0$ be an arbitrary point in $\mathcal{P}$. ```{r, attr.source='.numberLines'} a = c(1, 1, 1, 1); b = 2 # plane x+y+z+w = 2 x0 = c(b, b, b, b)/4 # a point in this plane nrml <- a/sqrt(c(crossprod(a))) # unit normal ``` Compute a mesh $\mathcal{M}_0$ of the isosurface $$ \bigl(s(u,v,w) - \mathbf{x}_0\bigr) \cdot \vec{\mathbf{n}} = 0. $$ ```{r, attr.source='.numberLines'} library(misc3d) f <- function(u, v, w){ c(crossprod(s(u, v, w), nrml)) } u_ <- v_ <- w_ <- seq(0, 2*pi, length.out = 100L) g <- expand.grid(u = u_, v = v_, w = w_) voxel <- array(with(g, f(u,v,w)), dim = c(100L,100L,100L)) surf <- computeContour3d(voxel, level = sum(x0*nrml), x = u_, y = v_, z = w_) trgls <- makeTriangles(surf) mesh0 <- misc3d:::t2ve(trgls) ``` Denote by $\mathcal{V}\mathcal{S}_0 \subset I \times J \times K$ the set of vertices of $\mathcal{M}_0$, and set $\mathcal{V}\mathcal{S} = s(\mathcal{V}\mathcal{S}_0) \subset \mathbb{R}^4$. ```{r, attr.source='.numberLines'} VS0 <- mesh0$vb VS <- s(VS0[1L,], VS0[2L,], VS0[3L,]) ``` Let $R$ be a rotation in $\mathbb{R}^4$ which sends $\vec{\mathbf{n}} =: \vec{\mathbf{v}}_1$ to the vector $(0,0,0,1) =: \vec{\mathbf{v}}_2$. One can take $R$ corresponding to the matrix $$ \frac{2}{{(\vec{\mathbf{v}}_1+\vec{\mathbf{v}}_2)}' (\vec{\mathbf{v}}_1+\vec{\mathbf{v}}_2)} (\vec{\mathbf{v}}_1+\vec{\mathbf{v}}_2) {(\vec{\mathbf{v}}_1+\vec{\mathbf{v}}_2)}' - I_4. $$ ```{r, attr.source='.numberLines'} rotationMatrix4D <- function(v1, v2){ v1 <- v1 / sqrt(c(crossprod(v1))) v2 <- v2 / sqrt(c(crossprod(v2))) 2*tcrossprod(v1+v2)/c(crossprod(v1+v2)) - diag(4L) } Rot <- rotationMatrix4D(nrml, c(0,0,0,1)) ``` Now define $\mathcal{V}\mathcal{S}' = R(\mathcal{V}\mathcal{S}) \subset \mathbb{R}^4$. Then all points in $\mathcal{V}\mathcal{S}'$ are equal on their fourth coordinate (up to numerical errors in R): ```{r, attr.source='.numberLines'} VSprime <- Rot %*% VS head(t(VSprime)) ``` Finally, define $\mathcal{V}\mathcal{S}'' \subset \mathbb{R}^3$ as the set obtained by removing the fourth coordinates of the elements of $\mathcal{V}\mathcal{S}'$, and define the mesh $\mathcal{M}$ whose set of vertices is $\mathcal{V}\mathcal{S}''$ and with the same edges as $\mathcal{M}_0$: ```{r, attr.source='.numberLines', eval=FALSE} library(rgl) mesh <- tmesh3d( vertices = VSprime[-4L,], indices = mesh0$ib, homogeneous = FALSE, normals = ? ) ``` What about the normals? If you have an implicit equation defining $\mathcal{S}$, that is, $\mathcal{S} = \iota^{-1}(0)$ with $\iota\colon \mathbb{R}^4 \to \mathbb{R}$, then a normal to $\mathcal{S}$ at a point $\mathbf{x} \in \mathbb{R}^4$ is given by the gradient of $\iota$ at $\mathbf{x}$. For the tiger, we know an implicit equation, and it is not difficult to get the gradient: ```{r, attr.source='.numberLines'} sNormal <- function(XYZT){ x <- XYZT[1L,]; y <- XYZT[2L,]; z <- XYZT[3L,]; t <- XYZT[4L,] rbind( x * (1 - R1/sqrt(x^2+y^2)), y * (1 - R1/sqrt(x^2+y^2)), z * (1 - R2/sqrt(z^2+t^2)), t * (1 - R2/sqrt(z^2+t^2)) ) } Normals <- sNormal(VS) ``` Once you get the normals: - project them to the hyperplane $\mathcal{P}$; - apply the rotation $R$ to the projected normals; - remove the fourth coordinates (all equal); - if necessary, negate the normals. The projection of $\mathbf{x} \in \mathbb{R}^4$ to the hyperplane $\mathcal{P}$ is given by $$ \mathbf{x} - \frac{\langle \mathbf{a}, \mathbf{x} \rangle - b}{\Vert \mathbf{a} \Vert^2} \mathbf{a}. $$ ```{r, echo=FALSE} library(rgl) ``` ```{r, attr.source='.numberLines'} # projection onto hyperplane = b projection <- function(a, b, X){ X - tcrossprod(a/c(crossprod(a)), colSums(a*X)-b) } mesh <- tmesh3d( vertices = VSprime[-4L,], indices = mesh0$ib, homogeneous = FALSE, normals = -t((Rot %*% projection(a, b, Normals))[-4L,]) ) ``` This works: ```{r, eval=FALSE} shade3d(mesh, color = "darkmagenta") ``` ![](figures/tiger1.png) Here is another way to get the normals. The normal at the point $s(u,v,w)$ is $$ \frac{\partial s}{\partial u}(u,v,w) \times \frac{\partial s}{\partial v}(u,v,w) \times \frac{\partial s}{\partial w}(u,v,w) $$ where $\cdot \times \cdot \times \cdot$ is the *ternary cross-product* in $\mathbb{R}^4$, defined by $$ \vec v_1 \times \vec v_2 \times \vec v_3 = \left\vert\begin{matrix} \vec i & \vec j & \vec k & \vec l \\ v_{1x} & v_{1y} & v_{1z} & v_{1t} \\ v_{2x} & v_{2y} & v_{2z} & v_{2t} \\ v_{3x} & v_{3y} & v_{3z} & v_{3t} \end{matrix}\right\vert. $$ ```{r, attr.source='.numberLines'} crossProd4D <- function(v1, v2, v3){ M <- rbind(v1, v2, v3) c(det(M[,-1L]), -det(M[,-2L]), det(M[,-3L]), -det(M[,-4L])) } sNormal <- function(uvw){ u <- uvw[1L]; v <- uvw[2L]; w <- uvw[3L] Du <- c((R1 + r*cos(w))*c(-sin(u),cos(u)), 0, 0) Dv <- c(0, 0, (R2 + r*sin(w))*c(-sin(v),cos(v))) Dw <- r * c(-sin(w)*c(cos(u),sin(u)), cos(w)*c(cos(v),sin(v))) crossProd4D(Du, Dv, Dw) } Normals <- apply(VS0, 2L, sNormal) ``` Then you can calculate the normals in this way and proceed as before: ```{r, attr.source='.numberLines'} mesh <- tmesh3d( vertices = VSprime[-4L,], indices = mesh0$ib, homogeneous = FALSE, normals = t((Rot %*% projection(a, b, Normals))[-4L,]) ) ``` Here is how to do an animation: ```{r, attr.source='.numberLines', eval=FALSE} b_ <- seq(-11.5, 11.5, length.out = 60L) open3d(windowRect = c(100, 100, 612, 612), zoom = 0.8) bg3d(rgb(54, 57, 64, maxColorValue = 255)) view3d(45, 40) for(i in 1L:length(b_)){ x0 <- rep(b_[i]/4, 4L) surf <- computeContour3d(voxel, level = sum(x0*nrml), x = u_, y = v_, z = w_) trgls <- makeTriangles(surf) mesh0 <- misc3d:::t2ve(trgls) VS0 <- mesh0$vb VS <- s(VS0[1L,], VS0[2L,], VS0[3L,]) Normals <- sNormal(VS) mesh <- tmesh3d( vertices = (Rot %*% VS)[-4L,], indices = mesh0$ib, homogeneous = FALSE, normals = -t((Rot %*% projection(a, b_[i], Normals))[-4L,]) ) shade3d(mesh, color = "firebrick3") snapshot3d(sprintf("pic%03d.png", i)) clear3d() } for(i in 1L:59L){ file.copy(sprintf("pic%03d.png", 60-i), sprintf("pic%03d.png", 60+i)) } # run gifski command <- "gifski --fps 12 pic*.png -o slicedTiger.gif" system(command) # cleaning pngfiles <- list.files(pattern = "^pic.*png$") file.remove(pngfiles) ``` ![](figures/tiger2.gif) ## Toroidal hyperboloid Let's give another example, a *toroidal hyperboloid*. This is a quadric with implicit equation $$ \frac{x^2}{a^2} + \frac{y^2}{b^2} - \frac{z^2}{c^2} - \frac{t^2}{d^2} = 1, $$ and a parameterization of this quadric is $$ \begin{array}{ccc} s \colon & (0,2\pi) \times (0,2\pi) \times (0, +\infty[ & \longrightarrow & \mathbb{R}^4 \\ & (u,v,w) & \longmapsto & \begin{pmatrix} a \cos u \cosh w \\ b \sin u \cosh w \\ c \cos v \sinh w \\ d \sin v \sinh w \end{pmatrix} \end{array}. $$ ![](figures/toroidalHyperboloid.gif)