SOCR ≫ DSPA ≫ Topics ≫

This DSPA section Appendix.3.1 (Primitive Surfaces with and without Boundaries) is part of the DSPA Appendix on visualization of geometric and parametric surfaces. This DSPA Appendix (3) covers the following 3 topics:

1 (Section 3.1) Geometric Primitive Surfaces with and without Boundaries

1.1 Saddle Point Surface

A saddle point is a point on the surface represented by a function \(z=z(x,y)\), where the gradient slopes (partial derivatives) of orthogonal function components defining the surface are trivial (\(0\)), however the point does not represent a local extremum on both axes, as may be expected. We already showed a simple example of a saddle point earlier, which we rendered as a surface. Now we will illustrate the parametric definition of a surface with a saddle point and it’s triangulated representation as a mesh3d.

library(plotly)
library (geometry)

# Plotly layout 
axs <- list(
  backgroundcolor="rgb(200,200,200)", # gray
  gridcolor="rgb(255,255,255)",       # white
  showbackground=TRUE,
  zerolinecolor="rgb(255,255,255)"     # white
)

n <- 36
h <- 1/(n-1)
r = seq(h, 1, length.out=n)
theta = seq(0, 2*pi, length.out=100)
 
grid.df <- expand.grid(r=r, theta=theta)
 
x <- c(grid.df$r * cos(grid.df$theta), 0)
y <- c(grid.df$r * sin(grid.df$theta), 0)
z <- sin(x*y)
 
mat <- matrix(
  c(x,y,z), 
  ncol = 3,
  dimnames = list(NULL, c("x", "y", "z"))
)
 
triangulated <- delaunayn(mat[,1:2])
 
# now figure out the colormap
zmean <- apply(triangulated, MARGIN=1, function(row){mean(mat[row,3])})
 
library(scales)

facecolor = colour_ramp(
  colorRampPalette(c("pink", "purple"))(20)
)(rescale(x=zmean))
 
plot_ly(
  x=x, y=y, z=z,
  i=triangulated[,1]-1, j=triangulated[,2]-1, k=triangulated[,3]-1,
  facecolor=facecolor,
  type="mesh3d",
  opacity = 0.7,
  contour=list(show=TRUE, color="#000", width=15)
) %>%
  layout(
    title="Triangulated Saddle Point surface",
    scene=list(
      xaxis=axs,
      yaxis=axs,
      zaxis=axs,
      camera=list(
        eye=list(x=1.75,y=-0.7,z=0.75)
      )
    )
  )

1.2 Geometric Shapes: 3D Spherical Parameterization

Below, we show three complementary examples of rendering synthetic geometric shapes; convex shapes (cone, sphere) and a non-convex surface (complex). It’s worthwhile reviewing the fundamentals of the spherical coordinate system representation.

# library(plotly)

# sweep or define (u,v) spherical coordinate parameter ranges
phi <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(200 - 1)))
psi <- seq(from = 0, to = pi, by = ((pi - 0)/(200 - 1)))

#p <- plot_ly(x = ~x, y = ~y, z = ~z, type = 'surface', opacity=1,
#             contour=list(show=TRUE, color="#000", width=15, lwd=10))  %>%
#  layout(title = paste("Layout ", shape), 
#         scene = list(xaxis=x_label,yaxis=y_label, zaxis=z_label))
#p

# shape=="complex")
    # rendering (u,v) parametric surfaces requires x,y,z arguments to be 2D arrays
    # In out case, the three coordinates have to be 200*200 parameterized tensors/arrays
    r1 = 2 + sin(3 * phi + 5 * psi)    # r = 2 + sin(7phi+5psi)
    x1 = (r1 * cos(phi)) %o% sin(psi)    # x = r*cos(phi)*sin(psi)
    y1 = (r1 * sin(phi)) %o% sin(psi)     # y = r*sin(phi)*sin(psi)
    z1 = r1 %o% cos(psi)  # z = r*cos(psi)
    
#shape=="cone")
    h2= 10   # cone height
    r2 = seq(from = 0, to = h2, by = ((h2 - 0)/(200 - 1)))  # r = radius
    x2 = 3* ((h2 - r2)/h2 ) %o% rep(1, 200)             # x = 3*r
    y2 = 3* ((h2 - r2)/h2 ) %o% sin(phi)   # y = r*sin(phi)
    z2 = 3* ((h2 - r2)/h2 ) %o% cos(phi)   # z = r*cos(phi)

#shape=="sphere") 
    r3 = 1                           # r = 1
    x3 = r3 * cos(phi) %o% sin(psi)   # x = r*cos(phi)*sin(psi)
    y3 = r3 * sin(phi) %o% sin(psi)   # y = r*sin(phi)*sin(psi)
    z3 = r3 * rep(1, 200) %o% cos(psi) # still need z to be 200*200 parameterized tensor/array

shape_names <- c("complex", "cone", "sphere")

# https://plot.ly/r/custom-buttons/

#p <- plot_ly(x = ~x, y = ~y, z = ~z, type = 'surface', opacity=1,
#             contour=list(show=TRUE, color="#000", width=15, lwd=10),
#             layout=layout_shapes)

# updatemenus component
updatemenus <- list(
  list(
    active = -1,
    type = 'buttons',
    buttons = list(
      list(
        label = shape_names[1],
        method = "update",
        args = list(list(visible = c(TRUE, FALSE, FALSE)),
                    list(title = shape_names[1]))),
      list(
        label = shape_names[2],
        method = "update",
        args = list(list(visible = c(FALSE, TRUE, FALSE)),
                    list(title = shape_names[2]))),
      list(
        label = shape_names[3],
        method = "update",
        args = list(list(visible = c(FALSE, FALSE, TRUE)),
                    list(title = shape_names[3])))
    )
  )
)

p <- plot_ly(hoverinfo="none", legendshow=FALSE, showscale = FALSE) %>%
  add_trace(x = ~x1, y = ~y1, z = ~z1, type = 'surface', opacity=1, visible=T,
             contour=list(show=TRUE, color="#000", width=15, lwd=10,
                          opacity=1.0, hoverinfo="none", legendshow=F)) %>%
  add_trace(x = ~x2, y = ~y2, z = ~z2, type='surface', opacity=1,visible=F,
             contour=list(show=TRUE, color="#000", width=15, lwd=10,
                          opacity=1.0, hoverinfo="none", legendshow=F)) %>%
  add_trace(x = ~x3, y = ~y3, z = ~z3, type = 'surface', opacity=0.7,visible=F,
             contour=list(show=TRUE, color="#000", width=15, lwd=10,
                          opacity=0.7, hoverinfo="none", legendshow=F)) %>%
  layout(title = "Choose a Shape", showlegend = FALSE,
         updatemenus = updatemenus)
p

The next example shows a bumpy surface parameterized in spherical coordinates. Complex-time indexed kimesurfaces, which extend time-series data in 5D spacekime, represent examples of such polar and spherical parameterizations of 2D manifolds. Mind the outer product \(\%o\%\) operator that takes two vectors and generates their outer product as a second order tensor (i.e., a matrix). This example displays the parametric surface modeled by this function: \[\rho=f(\phi,\theta) = 1 + \frac{1}{5}\left (\sin(5\phi)\ \%o\%\ \sin(7\theta)\right ) .\]

library(plotly)
phi <- seq(from = 0, to = 2*pi, length.out = 200)
theta <- seq(from = 0, to = pi/3, length.out = 200)

r1 = 1 + sin(5*phi) %o% (sin(7*theta)/5) # 2 + sin(3 * phi + 5 * theta)
x1 = r1 * (cos(phi) %o% sin(theta))    # x = r*cos(phi)*sin(theta)
y1 = r1 * (sin(phi) %o% sin(theta))     # y = r*sin(phi)*sin(theta)
z1 = r1 * (cos(theta) %o% rep(1,200))  # z = r*cos(theta)
surf_color1= phi %o% rep(1,200)
colorscale = cbind(seq(0, 1, by=1/(length(phi) - 1)), rainbow(length(phi)))

# Add Angular and Radial (polar-mesh) traces
# phi <- seq(from = 0, to = 2*pi, length.out = 10)
# theta <- seq(from = 0, to = pi/3, length.out = 8)
r = 1 + sin(5*phi) %o% (sin(7*theta)/5) # 2 + sin(3 * phi + 5 * theta)
x = r * (cos(phi) %o% sin(theta))    # x = r*cos(phi)*sin(theta)
y = r * (sin(phi) %o% sin(theta))     # y = r*sin(phi)*sin(theta)
z = r * (cos(theta) %o% rep(1,200))  # z = r*cos(theta)

# Radial (Time) mesh curves
count <- 50
phi <- seq(from = 0, to = 2*pi, length.out = count)
theta <- seq(from = 0, to = pi/3, length.out = count)
x <- array(1, dim=c(count*count))
y <- array(1, dim=c(count*count))
z <- array(1, dim=c(count*count))
r <- array(1, dim=c(count*count))
opaque <-  array(1, dim=c(count*count))

for (i in 1:count) {  # phi
  for (j in 1:count) { # theta
    r[(i-1)*count+j] <- 1 + sin(5*phi[i]) * (sin(7*theta[j])/5)
    x[(i-1)*count+j] <- r[(i-1)*count+j] * (cos(phi[i]) * sin(theta[j]))
    y[(i-1)*count+j] <- r[(i-1)*count+j] * (sin(phi[i]) * sin(theta[j]))
    z[(i-1)*count+j] <- r[(i-1)*count+j] * cos(theta[j])
    opaque[(i-1)*count+j] <- ifelse(j==10, 0, 1)
  }
}
length(x); length(y); length(r); dim(z); dim(opaque)
## [1] 2500
## [1] 2500
## [1] 2500
## [1] 2500
## [1] 2500
# Angular (Phase) mesh curves
count <- 50 # theta
phi <- seq(from = 0, to = 2*pi, length.out = count)
theta <- seq(from = 0, to = pi/3, length.out = count)
x2 <- array(1, dim=c(count*count))
y2 <- array(1, dim=c(count*count))
z2 <- array(1, dim=c(count*count))
r2 <- array(1, dim=c(count*count))
for (i in 1:count) {  # phi
  for (j in 1:count) { # theta
    r2[(i-1)*count+j] <- 1 + sin(5*phi[i]) * (sin(7*theta[j])/5)
    x2[(i-1)*count+j] <- r2[(i-1)*count+j] * (cos(phi[i]) * sin(theta[j]))
    y2[(i-1)*count+j] <- r2[(i-1)*count+j] * (sin(phi[i]) * sin(theta[j]))
    z2[(i-1)*count+j] <- r2[(i-1)*count+j] * cos(theta[j])
  }
}
length(x2); length(y2); length(r2); dim(z2)
## [1] 2500
## [1] 2500
## [1] 2500
## [1] 2500
phase_seq5 <- seq(from=25, to=count^2, by=count)  # phase mesh indexing patterns
phase_seq15 <- seq(from=40, to=count^2, by=count)
epsilon=0.03   # slight offset (up shift) of contour meshes

p <- plot_ly(hoverinfo="none", showscale = FALSE) %>%
    add_trace(x = ~x1, y = ~y1, z = ~r1,
              surfacecolor=~surf_color1, colorscale=colorscale, #Phase-based color
              type = 'surface', opacity=1, visible=T,
              contour=list(x = list(highlight = FALSE),
                           y = list(highlight = FALSE),
                           z = list( highlight = TRUE, highlightcolor = "blue"),
                           color="#000", width=15, lwd=10,
                           opacity=1.0, hoverinfo="none")) %>%
    # Add a few Radial Traces
    add_trace(x = x[(count+1):(2*count)], y = y[(count+1):(2*count)], 
              z = r[(count+1):(2*count)]+epsilon, 
            type = 'scatter3d', mode = 'lines', opacity = 1, # opaque, 
            line = list(color="#000", width=7, lwd=7), name = 'Time') %>%
    add_trace(x = x[(10*count+1):(11*count)], y = y[(10*count+1):(11*count)], 
              z = r[(10*count+1):(11*count)]+epsilon, 
            type = 'scatter3d', mode = 'lines', opacity = 1, # opaque, 
            line = list(color="#000", width=7, lwd=7), name = 'Time') %>%
    add_trace(x = x[(32*count+1):(33*count)], y = y[(32*count+1):(33*count)], 
              z = r[(32*count+1):(33*count)]+epsilon, 
            type = 'scatter3d', mode = 'lines', opacity = 1, # opaque, 
            line = list(color="#000", width=7, lwd=7), name = 'Time') %>%
    
    # Add a few Phase (Angular) Traces
    add_trace(x = x2[phase_seq5], y = y2[phase_seq5], z = r2[phase_seq5]+epsilon, 
            type = 'scatter3d', mode = 'lines', opacity = 1, # opaque, 
            line = list(color="darkgray", width=7, lwd=7), name = 'Phase') %>%
    add_trace(x = x2[phase_seq15], y = y2[phase_seq15], 
              z = r2[phase_seq15]+epsilon, 
            type = 'scatter3d', mode = 'lines', opacity = 1, # opaque, 
            line = list(color="darkgray", width=7, lwd=7), name = 'Phase') %>%
  
    # Configure the scene layout
    layout(title = "Spherical Surface Parameterization", showlegend = FALSE,
           scene = list(aspectmode = "manual", 
                        aspectratio = list(x=1, y=1, z=0.3))) # 1:1:1 aspect ratio
p

1.3 Cartesian Parameterization

Let’s demonstrate how to visualize as a surface in 3D a 2-parameter (Cartesian parametrization) function:

\[z=f(x,y) = (x-3)^2 +(y+4)^2 + x\ y .\] First, we will introduce the Cartesian parametrization of the \((x,y)\in D=\{-10\leq x,y\leq 10\} \subset\Re^2\) domain. This will be done by steps of size \(1\). Next we will express the symbolic definition of the function \(z=f(x,y)\). To plot the function using plotly, we will need to define the surface density (height, \(z\)) as an outer product defined on the \((x,y)\) grid, see outer() function. In the previous section, we showed how to use polar coordinate parametric grids. Mind that the \(z\) array needs to be transformed, as R stores matrices in column-wise order (not row-wise). Hence the function we actually plot is \(t(z)=z^t\).

library(plotly)
radius <- 10
x <- seq(from=-radius, to=radius, by=1)
y <- seq(from=-radius, to=radius, by=1) 
z_fun <- function(x,y) { (x-3)^2 +(y+4)^2 + x*y }
z    <- outer(x, y, FUN="z_fun")

# check surface function
dim(z); z[11+7,11-7]
## [1] 21 21
## [1] -24
plot_ly() %>%   add_trace(x=x, y=y, z=t(z), type="surface", opacity=0.7) %>%
  # add vertical axis at the minimum
  add_trace(x=7, y=-8, z=c(-24:350), type="scatter3d", mode="lines", 
              line = list(width = 6, color="navy blue"), name="Z",
              hoverinfo="none") %>%
  # add a Ball centered at arg_min
  add_trace(x=7, y=-8, z=-24, type="scatter3d", mode="markers")

1.4 Non-Convex Surfaces

1.4.1 Torus

A horizontal torus revolves around the z axis and is specified by two radii \(a\) and \(r\). A simple toral surface parameterization is

\[x = (a + r \cos(\phi)) \cos(\psi),\] \[y = (a + r \cos(\phi)) \sin(\psi),\] \[z = r \sin(\phi),\] where \(\phi,\psi \in [0, 2\pi)\) represent the angles about the z axis and about it’s core.

When \(a > r\), the torus does not intersect itself and the points on the torus satisfy \[\left ( a - \sqrt{x^2 + y^2}\right )^2 + z^2 = r^2.\] The 2D-torus is a 2-manifold (a surface) generated by revolving a circle in 3D about a given axis co-planar with the circle. Torus is the shape of an inner-tube; it is different from a solid torus, such as a doughnut or a bagel. The Torus is the boundary of a solid torus. It is homeomorphic to the Cartesian product of two circles \(S^1 \times S^1\) and represents a compact 2-manifold of genus 1.

# TORUS
a <- 6   # Torus radius (from torus gravitational center)
r <- 2   # Tube radius (from torus circular core)

phi <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(200 - 1)))
psi <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(200 - 1)))

xOuter <- outer(phi, psi, function(phi, psi) { cos(phi)*(a+r*cos(psi)) })
yOuter <- outer(phi, psi, function(phi, psi) { sin(phi)*(a+r*cos(psi)) })
zOuter <- outer(phi, psi, function(phi, psi) { r*sin(psi) })

plot_ly(hoverinfo="none", showscale = FALSE) %>%
    add_trace(x=~xOuter, y=~yOuter, z=~zOuter, type='surface', opacity=1, visible=T) %>%
    layout(title = "Torus", showlegend = FALSE,
          scene = list(aspectmode = "manual", aspectratio = list(x=2, y=2, z=1)))

We can introduce a checkerboard pattern on a coarser (parametric) surface representation of the torus.

# Checkerboard TORUS
a <- 6   # Torus radius (from torus gravitational center)
r <- 2   # Tube radius (from torus circular core)

phi <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(20 - 1)))
psi <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(20 - 1)))

xOuter <- outer(phi, psi, function(phi, psi) { cos(phi)*(a+r*cos(psi)) })
yOuter <- outer(phi, psi, function(phi, psi) { sin(phi)*(a+r*cos(psi)) })
zOuter <- outer(phi, psi, function(phi, psi) { r*sin(psi) })

# checkerboard color pattern
colorPair <- c('white','black')
surf_colors <- array("", dim = c(length(phi), length(phi)))
surf_colors_bin <- array(0, dim = c(length(phi), length(phi)))
# Check checkerboard pattern
for (j in 1:length(phi)) {
    for (i in 1:length(phi)) {
        surf_colors[j,i] <- colorPair[1+ (i+j) %% length(colorPair)]
        # surf_colors[j,i] <- colorPair[1+ ((-1+i+(j)/2) %% 10)]
        surf_colors_bin[j,i] <- ifelse(surf_colors[j,i]=="white", 1, 0)
    }
}

# Check checkerboard pattern
image(matrix(unclass(as.factor(surf_colors)), nrow=length(phi)))

plot_ly(hoverinfo="none", showscale = FALSE) %>%
    add_trace(x=~xOuter, y=~yOuter, z=~zOuter, type='surface', opacity=1, visible=T,
              surfacecolor = surf_colors_bin, opacity=0.9) %>%
    layout(title = "Torus", showlegend = FALSE,
          scene = list(aspectmode = "manual", aspectratio = list(x=2, y=2, z=1)))

Let’s go a step further and visualize a parametric curve on the torus.

a <- 6   # Torus radius (from torus gravitational center)
r <- 2   # Tube radius (from torus circular core)

# phi <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(20 - 1)))
# psi <- seq(from = 0, to = 2*pi, by = ((2*pi - 0)/(20 - 1)))

xOuter <- outer(phi, psi, function(phi, psi) { cos(phi)*(a+r*cos(psi)) })
yOuter <- outer(phi, psi, function(phi, psi) { sin(phi)*(a+r*cos(psi)) })
zOuter <- outer(phi, psi, function(phi, psi) { r*sin(psi) })

# checkerboard color pattern
colorPair <- c('white','black')
surf_colors <- array("", dim = c(length(phi), length(phi)))
# Check checkerboard pattern
for (j in 1:length(phi)) {
    for (i in 1:length(phi)) {
        surf_colors[j,i] <- colorPair[1+ (i+j) %% length(colorPair)]
        # surf_colors[j,i] <- colorPair[1+ ((-1+i+(j)/2) %% 10)]
    }
}
# Check checkerboard pattern
image(matrix(unclass(as.factor(surf_colors)), nrow=length(phi)))

### Kime tube on torus
#### Parametric curve on torus
# k1 <- phi -pi
# k2 <- ((phi-1)^3 + 5*(phi-1)^2 + 2*(phi-5) - 8)/10
xCurveOuter <- outer(phi, psi, function(phi, psi) { 
  cos(phi -pi)*(a+r*cos(((phi -pi-1)^3 + 5*(phi -pi-1)^2 + 2*(phi -pi-5) - 8)/10)) })
yCurveOuter <- outer(phi, psi, function(phi, psi) { 
  sin(phi -pi)*(a+r*cos(((phi -pi-1)^3 + 5*(phi -pi-1)^2 + 2*(phi -pi-5) - 8)/10)) })
zCurveOuter <- outer(phi, psi, function(phi, psi) { 
  r*sin(((phi -pi-1)^3 + 5*(phi -pi-1)^2 + 2*(phi -pi-5) - 8)/10) })

plot_ly(x=~as.vector(xCurveOuter), y=~as.vector(yCurveOuter), 
        z=~as.vector(zCurveOuter), type="scatter3d", mode="markers+lines",
        line=list(color='red', width=15))
plot_ly(hoverinfo="none") %>%
  add_trace(x=~xOuter, y=~yOuter, z=~zOuter, type='surface', 
            opacity=0.1, visible=TRUE, showlegend= FALSE, showscale = FALSE) %>%
  add_trace(x=~as.vector(xCurveOuter), y=~as.vector(yCurveOuter), z=~as.vector(zCurveOuter),
            name="Tube", type='scatter3d', mode='markers+lines', showscale = FALSE,
            hoverinfo = 'Tube', line=list(color='red', width=15), showlegend = FALSE) %>%
  # add projection on plane x=8.5
  add_trace(x = 8.5, y = ~as.vector(yCurveOuter), z=~as.vector(zCurveOuter), 
            type="scatter3d", mode="lines", line = list(width = 2, dash="solid", color = ~c, 
            colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>%
  # add projection on plane y=-8.5
  add_trace(x = ~as.vector(xCurveOuter), y = -8.5, z = ~as.vector(xCurveOuter), 
            type="scatter3d", mode="lines", line = list(width = 2, dash="solid", color = ~c,
            colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>%
  # add projection on plane z=-2.5
  add_trace(x = ~as.vector(xCurveOuter), y = ~as.vector(yCurveOuter), z = ~-2.5, 
            type="scatter3d", mode="lines", line = list(width = 2, dash="solid", color = ~c,
            colorscale = list(c(0,'gray'), c(1,'black'))), name="Z", hoverinfo="none") %>%
  layout(title = "Torus", showlegend = FALSE,
          scene = list(aspectmode = "manual", aspectratio = list(x=2, y=2, z=1)))

1.4.2 Objective function optimization

Non-convex optimization is NP-hard problem that is very difficult to solve in general. The example below illustrates such non-convex optimization using complex surfaces. The plot_ly graph below shows an interactive surface plot of the objective function along with the trajectories of several optimization-solutions. Mind the different optimization trajectories that may either smoothly approach the local minima (surface sulci or valleys) or rapidly switch between neighboring surface slopes or even jump across different crests.

# 1. Initialize the gradient descent optimization process
# install.packages("numDeriv")
library(numDeriv)         # to access the gradient calculation function: grad()

gradDescent = function(f, x0, max.iter=500, stepSize=0.01, stoppingCriterion=0.001) {
  n = length(x0)
  xmat = matrix(0, nrow=n, ncol=max.iter)
  xmat[,1] = x0
  
  for (k in 2:max.iter) {
    # Calculate the current gradient
    currentGrad = grad(f, xmat[,k-1]) 
    
    # Check Stopping criterion
    if (all(abs(currentGrad) < stoppingCriterion)) {
      k = k-1; break
    }
    
    # Move in the opposite direction of the gradient
    xmat[,k] = xmat[,k-1] - stepSize * currentGrad
  }

  xmat = xmat[,1:k]     # remove the initial guess (x0)
  # return the final solution, the solution path, min.value, and the number of iterations
  return(list(x=xmat[,k], xmat=xmat, min.val=f(xmat[,k]), k=k))
}

# 2. Define the objective function and use a manual optimization gradient-descent scheme to get solutions
f1 = function(x) {
  return((1/2*x[1]^2-1/4*x[2]^2+3)*cos(2*x[1]+1-exp(x[2])))
}

x0 = c(0.5,0.5)
optim.f0 = gradDescent(f1,x0,stepSize=0.01)
x1 = c(-0.1,-1.3)
optim.f1 = gradDescent(f1,x1,stepSize=0.01,max.iter=400)
x2 = c(-0.5,-1.3)
optim.f2 = gradDescent(f1,x2,stepSize=0.01,max.iter=400)
x3 = c(-0.2,1.4)
optim.f3 = gradDescent(f1,x3,stepSize=0.01,max.iter=400)
x4 = c(-0.5,-0.5)
optim.f4 = gradDescent(f1,x4,stepSize=0.01,max.iter=400)
x5 = c(-1.7, 1.45)
optim.f5 = gradDescent(f1,x5,stepSize=0.01,max.iter=400)

# 3. PLot the 3D scene
library(magrittr)
library(plotly)

radius <- 2
x <- seq(from=-2*radius, to=2*radius+1, by=0.05)
y <- seq(from=-radius, to=radius+1, by=0.05) 
z_fun <- function(x,y) { return((1/2*x^2-1/4*y^2+3)*cos(2*x+1-exp(y))) }
z <- outer(x, y, FUN="z_fun")

plot_ly() %>%   
  add_trace(x=x, y=y, z=t(z), type="surface", opacity=0.7, showlegend = FALSE) %>%
  add_trace(x = optim.f0$xmat[1,], y = optim.f0$xmat[2,], z = apply(optim.f0$xmat,2,f1), 
            type = 'scatter3d', mode = 'lines', opacity = 1, 
            line = list(width = 4, color = "red"), name = 'Solution 1') %>%
  add_trace(x = optim.f1$xmat[1,], y = optim.f1$xmat[2,], z = apply(optim.f1$xmat,2,f1), 
            type = 'scatter3d', mode = 'lines', opacity = 1, 
            line = list(width = 4, color = "green"), name = 'Solution 2') %>%
  add_trace(x = optim.f2$xmat[1,], y = optim.f2$xmat[2,], z = apply(optim.f2$xmat,2,f1), 
            type = 'scatter3d', mode = 'lines', opacity = 1, 
            line = list(width = 4, color = "blue"), name = 'Solution 3') %>%
  add_trace(x = optim.f3$xmat[1,], y = optim.f3$xmat[2,], z = apply(optim.f3$xmat,2,f1), 
            type = 'scatter3d', mode = 'lines', opacity = 1, 
            line = list(width = 4, color = "purple"), name = 'Solution 4') %>%
  add_trace(x = optim.f4$xmat[1,], y = optim.f4$xmat[2,], z = apply(optim.f4$xmat,2,f1), 
            type = 'scatter3d', mode = 'lines', opacity = 1, 
            line = list(width = 4, color = "orange"), name = 'Solution 5') %>%
  add_trace(x = optim.f5$xmat[1,], y = optim.f5$xmat[2,], z = apply(optim.f5$xmat,2,f1), 
            type = 'scatter3d', mode = 'lines', opacity = 1, 
            line = list(width = 4, color = "gray"), name = 'Solution 6')

1.4.3 Human face mask surface

The (x,y,z) coordinates of a 3D surface representing a face mask are available here (CSV). We can display this (non-convex) surface in 3D using plot_ly().

# 1. Load the vertex-based surface data (x,y,z) locations on the surface are obtained by parameterizing 
#    real 3D face scans
library(plotly)
vertexMaskDataXYZ <- read.csv("https://umich.instructure.com/files/13417342/download?download_frd=1", header = T)
dim(vertexMaskDataXYZ) # [1] 7159    3
## [1] 7159    3
# 2. Plot the 3D scene
library(magrittr)
library(plotly)

y <- 200*sin(seq(0, 2*pi, length.out = 200))*cos(seq(0, 2*pi, length.out = 200))

plot_ly() %>%   
  add_trace(x=vertexMaskDataXYZ[ , 1], y=vertexMaskDataXYZ[ , 2], 
            z=vertexMaskDataXYZ[ , 3], type="mesh3d", # colors = colorRamp(rainbow(8)), 
            opacity=0.7, name="Face Mask", hoverinfo="none",
            contour=list(show=TRUE, color="#000", width=15)) %>%
  # trace the boundary of the saddle point surface
  add_trace(x=~vertexMaskDataXYZ[ , 1], y=~vertexMaskDataXYZ[ , 2], 
            z=15, type="scatter3d", mode="lines", 
              line = list(width = 10, color="red"), name="Surface (X,Y) Domain" # , hoverinfo="none"
            ) %>%
  # trace the main Z-axis
  add_trace(x=-8, y=~y, z= 20, type="scatter3d", mode="lines", 
              line = list(width = 20, color="navy blue"), name="Z",
              name="Mask Pole") %>%
  layout(title = "Human Face Surface", showlegend = FALSE,
         scene = list(
          xaxis = list(title = "X"),
          yaxis = list(title = "Y"),
          zaxis = list(title = "Z")
        ))

2 References

SOCR Resource Visitor number Web Analytics SOCR Email