---
title: "Parametric surface in Haskell OpenGL, with surface normals"
author: "Stéphane Laurent"
date: "2018-10-23"
output:
  md_document:
    variant: markdown
    preserve_yaml: true
  html_document:
    keep_md: no
prettify: yes
linenums: yes
prettifycss: twitter-bootstrap
tags: haskell, graphics, opengl
highlighter: kate
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(fig.path = "./figures/OpenglParametricWithNormals-")
```

Similarly to [a previous post](https://laustep.github.io/stlahblog/posts/OpenglParametricPlot.html), 
I will show here how to draw a parametric surface with the Haskell OpenGL 
library, but this time we will include the surface normal at each vertex.

As the example of a surface, I take the stereographic projection of a Hopf torus. 
The parameterization is given by the function defined as follows in Haskell:

```haskell
type Point = (Double, Double, Double)

hopf :: Double -> Double -> Point
hopf v = (x1/(den-x4), x2/(den-x4), x3/(den-x4))
  where
    a = 0.44
    nlobes = 3
    a1 = pi/2 - (pi/2-a) * cos(u*nlobes)
    a2 = u + a*sin(2*u*nlobes)
    sina1 = sin a1
    p1 = cos a1
    p2 = sina1 * cos a2
    p3 = sina1 * sin a2
    cosphi = cos v
    sinphi = sin v
    x1 = cosphi*p3 + sinphi*p2
    x2 = cosphi*p2 - sinphi*p3
    x3 = sinphi * (1+p1)
    x4 = cosphi * (1+p1)
    den = sqrt(2*(1+p1))
```

for $0 \leq u < 2\pi$ and $0 \leq v < 2\pi$. 

We will evaluate this function at the vertices of a grid like the one shown below 
(we will see later why we show the six red triangles on this picture):

```{r grid, echo=FALSE, fig.width=8, fig.height=8}
par(mar=c(5,5.5,1,1))
plot(0,0, xlim=c(0,10), ylim=c(0,10), type="n", axes=FALSE, 
     xlab=NA, ylab=NA, asp=1)
mtext(expression(italic(u)~group("(", italic(n[u]) == 5, ")")), side=1, line=4)
mtext(expression(italic(v)~group("(", italic(n[v]) == 5, ")")), side=2, line=4)
grid()
axis(1, at=c(0,2,4,6,8,10), 
     labels = c(0, 
                expression(over(2~pi, italic(n[u]))),
                expression(over(2 %*% 2~pi, italic(n[u]))),
                expression(over(3 %*% 2~pi, italic(n[u]))),
                expression(over(4 %*% 2~pi, italic(n[u]))),
                expression(2~pi)), 
     padj=0.5)
axis(2, at=c(0,2,4,6,8,10), 
     labels = c(0, 
                expression(over(2~pi, italic(n[v]))),
                expression(over(2 %*% 2~pi, italic(n[v]))),
                expression(over(3 %*% 2~pi, italic(n[v]))),
                expression(over(4 %*% 2~pi, italic(n[v]))),
                expression(2~pi)), 
     padj=0.1)

text(6, 4, 
     expression(group("(",list(italic(u[i]), italic(v[j])),")")), 
     cex=0.9, adj=c(1,1))
text(8, 4, 
     expression(group("(",list(italic(u[i+1]), italic(v[j])),")")), 
     cex=0.9, adj=c(0,1))
text(6, 6, 
     expression(group("(",list(italic(u[i]), italic(v[j+1])),")")), 
     cex=0.9, adj=c(1,0))
text(8, 6, 
     expression(group("(",list(italic(u[i+1]), italic(v[j+1])),")")), 
     cex=0.9, adj=c(0,0))
polygon(c(6,6,8,8), c(4,6,6,4))
segments(6, 4, 6, 6, col="red", lty="dashed", lwd=2)
segments(6, 4, 8, 4, col="red", lty="dashed", lwd=2)
segments(8, 4, 6, 6, col="red", lty="dashed", lwd=2)
segments(6, 4, 8, 2, col="red", lty="dashed", lwd=2)
segments(6, 4, 6, 2, col="red", lty="dashed", lwd=2)
segments(8, 2, 6, 2, col="red", lty="dashed", lwd=2)
segments(6, 4, 4, 4, col="red", lty="dashed", lwd=2)
segments(6, 4, 4, 6, col="red", lty="dashed", lwd=2)
segments(4, 4, 4, 6, col="red", lty="dashed", lwd=2)
segments(8, 4, 8, 2, col="red", lty="dashed", lwd=2)
segments(6, 2, 4, 4, col="red", lty="dashed", lwd=2)
segments(4, 6, 6, 6, col="red", lty="dashed", lwd=2)
```


We write a function that evaluates the values of a parametrization at the point 
of this grid and put them in an array:

```haskell
import           Data.Array      (Array, (!), array)
import qualified Data.Array as A

frac :: Int -> Int -> Double
frac p q = realToFrac p / realToFrac q

allVertices :: (Double -> Double -> Point) -> (Int,Int) -> Array (Int,Int) Point
allVertices f (n_u, n_v) = array ((0,0), (n_u-1,n_v-1)) associations
  where
  u_ = [2*pi * frac i n_u | i <- [0 .. n_u-1]]
  v_ = [2*pi * frac i n_v | i <- [0 .. n_v-1]]
  indices = [(i,j) | i <- [0 .. n_u-1], j <- [0 .. n_v-1]]
  g (i,j) = ((i,j), f (u_ !! i) (v_ !! j))
  associations = map g indices
```

These values are the surface vertices. 
Now, we write a function that approximates the surface normal at vertex $(i,j)$. 
This normal approximately is the average of the normals of the six triangles 
incident to the vertex. 

```haskell
type Vector = (Double, Double, Double)

triangleNormal :: (Point, Point, Point) -> Vector
triangleNormal ((x1,x2,x3), (y1,y2,y3), (z1,z2,z3)) = (a/norm, b/norm, c/norm)
  where
    (a, b, c) = crossProd (z1-x1, z2-x2, z3-x3) (y1-x1, y2-x2, y3-x3) 
    crossProd (a1,a2,a3) (b1,b2,b3) = (a2*b3-a3*b2, a3*b1-a1*b3, a1*b2-a2*b1)
    norm = sqrt(a*a + b*b + c*c)

averageNormals :: Vector -> Vector -> Vector -> Vector -> Vector -> Vector -> Vector
averageNormals (x1,y1,z1) (x2,y2,z2) (x3,y3,z3) (x4,y4,z4) (x5,y5,z5) (x6,y6,z6) = 
  ((x1+x2+x3+x4+x5+x6)/6, (y1+y2+y3+y4+y5+y6)/6, (z1+z2+z3+z4+z5+z6)/6)

normal_ij :: Array (Int,Int) Point -> (Int, Int) -> Vector
normal_ij vertices (i,j) = averageNormals n1 n2 n3 n4 n5 n6
  where
  ((_,_), (n_u',n_v')) = A.bounds vertices
  im1 = if i==0 then n_u' else i-1
  ip1 = if i==n_u' then 0 else i+1
  jm1 = if j==0 then n_v' else j-1
  jp1 = if j==n_v' then 0 else j+1
  n1 = triangleNormal (vertices ! (i,j), vertices ! (i,jp1), vertices ! (ip1,j))
  n2 = triangleNormal (vertices ! (i,j), vertices ! (ip1,jm1), vertices ! (i,jm1))
  n3 = triangleNormal (vertices ! (i,j), vertices ! (im1,j), vertices ! (im1,jp1))
  n4 = triangleNormal (vertices ! (i,j), vertices ! (ip1,j), vertices ! (ip1,jm1))
  n5 = triangleNormal (vertices ! (i,j), vertices ! (i,jm1), vertices ! (im1,j))
  n6 = triangleNormal (vertices ! (i,j), vertices ! (im1,jp1), vertices ! (i,jp1))
```

Now we write a function that takes the array of surface vertices as input and 
returns an array containing the surface normals:

```haskell
allNormals :: Array (Int,Int) Point -> Array (Int,Int) Vector
allNormals vertices = array bounds associations
  where
  bounds = A.bounds vertices
  indices = A.indices vertices  
  associations = map (\(i,j) -> ((i,j), normal_ij vertices (i,j))) indices
```

Let's say that a surface triangle whose each vertex is attached to the  
corresponding surface normal is a n-triangle. 
To each vertex $(i,j)$, we associate two n-triangles: 
the lower n-triangle for vertices $(i,j)$-$(i+1,j)$-$(i,j+1)$ and 
the upper n-triangle for vertices $(i+1,j+1)$-$(i,j+1)$-$(i+1,j)$. 
We write a function that takes as input the two arrays (vertices and normals), 
an index $(i,j)$, and that returns the two n-triangles:

```haskell
import Graphics.Rendering.OpenGL.GL (Normal3 (..), Vertex3 (..))

type NPoint = (Vertex3 Double, Normal3 Double)
type NTriangle = (NPoint, NPoint, NPoint)

pointToVertex3 :: Point -> Vertex3 Double
pointToVertex3 (x,y,z) = Vertex3 x y z

vectorToNormal3 :: Vector -> Normal3 Double
vectorToNormal3 (x,y,z) = Normal3 x y z

triangles_ij :: Array (Int,Int) Point -> Array (Int,Int) Vector 
            -> (Int, Int) -> (Int, Int)
            -> (NTriangle, NTriangle)
triangles_ij vertices normals (n_u,n_v) (i,j) = 
  (((a,na), (b,nb), (c,nc)), ((c,nc), (b,nb), (d,nd)))
  where
  ip1 = if i==n_u-1 then 0 else i+1
  jp1 = if j==n_v-1 then 0 else j+1
  a = pointToVertex3 $ vertices ! (i,j)
  na = vectorToNormal3 $ normals ! (i,j)
  c = pointToVertex3 $ vertices ! (i,jp1)
  nc = vectorToNormal3 $ normals ! (i,jp1)
  d = pointToVertex3 $ vertices ! (ip1,jp1)
  nd = vectorToNormal3 $ normals ! (ip1,jp1)
  b = pointToVertex3 $ vertices ! (ip1,j)
  nb = vectorToNormal3 $ normals ! (ip1,j)
```

Finally, we write a function returning the list of all pairs of n-triangles:

```haskell
allTriangles :: (Int,Int) -> [(NTriangle,NTriangle)]
allTriangles (n_u,n_v) =
  map (triangles_ij vertices normals (n_u,n_v)) indices
  where
  vertices = allVertices hopf (n_u,n_v)
  normals = allNormals vertices
  indices = [(i,j) | i <- [0 .. n_u-1], j <- [0 .. n_v-1]]
```

Done. It remains to write the OpenGL side:

```haskell
import Data.IORef
import Graphics.Rendering.OpenGL.GL
import Graphics.UI.GLUT

hopfTorus :: [(NTriangle,NTriangle)]
hopfTorus = allTriangles (400,400)

data Context = Context
    {
      contextRot1      :: IORef GLfloat
    , contextRot2      :: IORef GLfloat
    , contextRot3      :: IORef GLfloat
    , contextTriangles :: IORef [(NTriangle,NTriangle)]
    }

white,black,pink :: Color4 GLfloat
white      = Color4    1   1   1    1
black      = Color4    0   0   0    1
pink       = Color4    1   0   0.5  1

display :: Context -> IORef GLdouble -> DisplayCallback
display context zoom alpha = do
  clear [ColorBuffer, DepthBuffer]
  r1 <- get (contextRot1 context)
  r2 <- get (contextRot2 context)
  r3 <- get (contextRot3 context)
  ntriangles <- get (contextTriangles context)
  let ntriangles' = unzip ntriangles
      lowerTriangles = fst ntriangles'
      upperTriangles = snd ntriangles'
  z <- get zoom
  loadIdentity
  (_, size) <- get viewport
  resize z size
  rotate r1 $ Vector3 1 0 0
  rotate r2 $ Vector3 0 1 0
  rotate r3 $ Vector3 0 0 1
  renderPrimitive Triangles $ mapM_ drawTriangle lowerTriangles
  renderPrimitive Triangles $ mapM_ drawTriangle lowerTriangles
  swapBuffers
  where
    drawTriangle ((v1,n1),(v2,n2),(v3,n3)) = do
      materialDiffuse Front $= pink
      normal n1
      vertex v1
      normal n2
      vertex v2
      normal n3
      vertex v3

resize :: GLdouble -> Size -> IO ()
resize zoom s@(Size w h) = do
  viewport $= (Position 0 0, s)
  matrixMode $= Projection
  loadIdentity
  perspective 45.0 (realToFrac w / realToFrac h) 1.0 100.0
  lookAt (Vertex3 0 0 (24+zoom)) (Vertex3 0 0 0) (Vector3 0 1 0)
  matrixMode $= Modelview 0

keyboard :: IORef GLfloat -> IORef GLfloat -> IORef GLfloat -- rotations
         -> IORef GLdouble -- zoom
         -> IORef [(NTriangle,NTriangle)]
         -> KeyboardCallback
keyboard rot1 rot2 rot3 zoom triangles c _ = do
  case c of
    'e' -> rot1 $~! subtract 2
    'r' -> rot1 $~! (+2)
    't' -> rot2 $~! subtract 2
    'y' -> rot2 $~! (+2)
    'u' -> rot3 $~! subtract 2
    'i' -> rot3 $~! (+2)
    'm' -> zoom $~! (+0.1)
    'l' -> zoom $~! subtract 0.1
    'q' -> leaveMainLoop
    _   -> return ()
  postRedisplay Nothing

main :: IO ()
main = do
  _ <- getArgsAndInitialize
  _ <- createWindow "Hopf torus"
  windowSize $= Size 500 500
  initialDisplayMode $= [RGBAMode, DoubleBuffered, WithDepthBuffer]
  clearColor $= black
  materialAmbient Front $= black
  lighting $= Enabled
  light (Light 0) $= Enabled
  position (Light 0) $= Vertex4 0 0 (-1000) 1
  ambient (Light 0) $= white
  diffuse (Light 0) $= white
  specular (Light 0) $= white
  depthFunc $= Just Less
  shadeModel $= Smooth
  rot1 <- newIORef 0.0
  rot2 <- newIORef 0.0
  rot3 <- newIORef 0.0
  zoom <- newIORef 0.0
  nlobes' <- newIORef nlobes
  hopfTorus' <- newIORef hopfTorus
  displayCallback $= display Context {contextRot1 = rot1,
                                      contextRot2 = rot2,
                                      contextRot3 = rot3,
                                      contextTriangles = hopfTorus'}
                             zoom 
  reshapeCallback $= Just (resize 0)
  keyboardCallback $= Just (keyboard rot1 rot2 rot3 zoom hopfTorus')
  idleCallback $= Nothing
  putStrLn "*** Hopf torus ***\n\
        \    To quit, press q.\n\
        \    Scene rotation: e, r, t, y, u, i\n\
        \    Zoom: l, m\n\
        \"
  mainLoop
```

And this is the result:

![](./figures/OpenglParametricWithNormals-Result.png)


The full code is available in [this Github repo](https://github.com/stla/opengl-HopfTorus). 