--- title: "Parametric surface with OpenGL in Haskell" author: "Stéphane Laurent" date: "2018-03-27" 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/OpenglParametric-")  In this blog post, I will show how to draw a parametric surface with the Haskell OpenGL library. We will take as example the following function with two variables$u$and$v$and depending on a parameter$n$: $$f_n(u,v) = \begin{pmatrix} \sin\Bigl(n v + \dfrac{\pi}{2}\Bigr) \left(\Bigl(1-\dfrac{v}{2\pi}\Bigr)(1 - \cos u) -0.2 \right) \\ \dfrac{4v}{\pi} \Bigl(1-\dfrac{v}{2\pi}\Bigr) \sin u \\ \cos\Bigl(n v + \dfrac{\pi}{2}\Bigr) \left(\Bigl(1-\dfrac{v}{2\pi}\Bigr)(1 - \cos u) +0.2 \right) \end{pmatrix}$$ for$0 \leq u \leq 2\pi$and$0 \leq v \leq 2\pi$. It draws a conical spiral (see the images below). Firslty, we define this function in Haskell: haskell function :: Double -- parameter n -> Double -> Double -- variables u and v -> [Double] -- result vector function n u v = [x, y, z] where x = sin(n*v+0.5*pi) * ((1-0.5*v/pi) * (1-cos u) - 0.2) y = 4*v/pi + (1-0.5*v/pi) * sin u z = cos(n*v+0.5*pi) * ((1-0.5*v/pi) * (1-cos u) + 0.2)  We will evaluate the function$f_n$at the corners of each square of a grid like the one show below: {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))  This is done by the Haskell function quad below, which returns a pair made of a four-tuple of Vertex3 and the normal. To get the normal, we use the helper function triangleNormal which computes the cross product. haskell import Graphics.Rendering.OpenGL.GL (Normal3 (..), Vertex3 (..)) triangleNormal :: Floating a => (Vertex3 a, Vertex3 a, Vertex3 a) -> Normal3 a triangleNormal (Vertex3 x1 x2 x3, Vertex3 y1 y2 y3, Vertex3 z1 z2 z3) = Normal3 (a/norm) (b/norm) (c/norm) where (a, b, c) = crossProd (y1-x1, y2-x2, y3-x3) (z1-x1, z2-x2, z3-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) quad :: (Double -> Double -> Double -> [Double]) -- the function -> Double -- parameter n -> [Double] -> [Double] -- sequences of u and v -> Int -> Int -- indices -> ((Vertex3 Double, Vertex3 Double, Vertex3 Double, Vertex3 Double), Normal3 Double) quad f n u_ v_ i j = ((a, b, c, d), norm) where (a,b,c,d) = ( toVx3$ f n (u_!!i) (v_!!j) , toVx3 $f n (u_!!i) (v_!!(j+1)) , toVx3$ f n (u_!!(i+1)) (v_!!(j+1)) , toVx3 $f n (u_!!(i+1)) (v_!!j) ) norm = triangleNormal (a, b, c) toVx3 x = Vertex3 (x!!0) (x!!1) (x!!2)  Now we compute all the quad's for two given integers defining the grid: haskell allQuads :: Int -> Int -- numbers of subdivisions for u and v -> Double -- parameter n -> [((Vertex3 Double, Vertex3 Double, Vertex3 Double, Vertex3 Double), Normal3 Double)] allQuads n_u n_v n = map (uncurry (quad function n seq_u seq_v)) [(i,j) | i <- [0 .. n_u-1], j <- [0 .. n_v-1]] where seq_u,seq_v :: [Double] seq_u = [2*pi * frac i n_u | i <- [0 .. n_u]] seq_v = [2*pi * frac i n_v | i <- [0 .. n_v]] frac :: Int -> Int -> Double frac p q = realToFrac p / realToFrac q  Now we make a standard OpenGL program. We allow the user to use the keyboard for some controls: rotating the scene and increasing/decreasing the parameter$n$. haskell module ConicalSpiral.ConicalSpiral2 where import ConicalSpiral.Data2 import Data.IORef import Graphics.Rendering.OpenGL.GL import Graphics.UI.GLUT white,black,red :: Color4 GLfloat white = Color4 1 1 1 1 black = Color4 0 0 0 1 red = Color4 1 0 0 1 display :: IORef GLfloat -> IORef GLfloat -> IORef GLfloat -- rotations -> IORef GLdouble -- parameter n -> IORef GLdouble -- zoom -> DisplayCallback display rot1 rot2 rot3 n zoom = do clear [ColorBuffer, DepthBuffer] n' <- get n r1 <- get rot1 r2 <- get rot2 r3 <- get rot3 z <- get zoom (_, size) <- get viewport let surface = allQuads 200 200 n' loadIdentity resize z size rotate r1$ Vector3 1 0 0 rotate r2 $Vector3 0 1 0 rotate r3$ Vector3 0 0 1 renderPrimitive Quads $do materialDiffuse FrontAndBack$= red mapM_ drawQuad surface swapBuffers where drawQuad ((v1,v2,v3,v4),norm) = do normal norm vertex v1 vertex v2 vertex v3 vertex v4 resize :: Double -> Size -> IO () resize zoom s@(Size w h) = do viewport $= (Position 0 0, s) matrixMode$= Projection loadIdentity perspective 45.0 (w'/h') 1.0 100.0 lookAt (Vertex3 0 0 (-15+zoom)) (Vertex3 0 2 0) (Vector3 0 1 0) matrixMode $= Modelview 0 where w' = realToFrac w h' = realToFrac h keyboard :: IORef GLfloat -> IORef GLfloat -> IORef GLfloat -- rotations -> IORef GLdouble -- parameter n -> IORef GLdouble -- zoom -> KeyboardCallback keyboard rot1 rot2 rot3 n zoom c _ = case c of 'b' -> n$~! subtract 0.1 'n' -> n $~! (+ 0.1) 'e' -> rot1$~! subtract 1 'r' -> rot1 $~! (+1) 't' -> rot2$~! subtract 1 'y' -> rot2 $~! (+1) 'u' -> rot3$~! subtract 1 'i' -> rot3 $~! (+1) 'm' -> zoom$~! (+1) 'l' -> zoom $~! subtract 1 'q' -> leaveMainLoop _ -> return () idle :: IdleCallback idle = postRedisplay Nothing main :: IO () main = do _ <- getArgsAndInitialize _ <- createWindow "Conical Surface" windowSize$= Size 500 500 initialDisplayMode $= [RGBAMode, DoubleBuffered, WithDepthBuffer] clearColor$= white materialAmbient FrontAndBack $= black materialShininess FrontAndBack$= 95 materialSpecular FrontAndBack $= white lighting$= Enabled lightModelTwoSide $= Enabled light (Light 0)$= Enabled position (Light 0) $= Vertex4 0 0 (-100) 1 ambient (Light 0)$= black diffuse (Light 0) $= white specular (Light 0)$= white depthFunc $= Just Less shadeModel$= Smooth n <- newIORef 5.0 rot1 <- newIORef 0.0 rot2 <- newIORef 0.0 rot3 <- newIORef 0.0 zoom <- newIORef 0.0 displayCallback $= display rot1 rot2 rot3 n zoom reshapeCallback$= Just (resize 0) keyboardCallback $= Just (keyboard rot1 rot2 rot3 n zoom) idleCallback$= Just idle putStrLn "*** Haskell OpenGL Conical Spiral ***\n\ \ To quit, press q.\n\ \ Scene rotation:\n\ \ e, r, t, y, u, i\n\ \ Zoom: l, m\n\ \ Increase/decrease parameter: n, b \n\ \" mainLoop  And this is the result: ![](figures/OpenGLredCone.png) We have chosen to render the conical spiral in red. We can include colors. A possible way to do so is to compute the average value of the quad, normalize it and consider the result as a RGB code: haskell import Graphics.Rendering.OpenGL.GL (Color4 (..), GLfloat, Vertex3 (..)) type Quad = (Vertex3 Double, Vertex3 Double, Vertex3 Double, Vertex3 Double) colMeans :: Quad -> (Double, Double, Double) colMeans quad = (x/m , y/m, z/m) where (a, b, c, d) = quad Vertex3 a1 a2 a3 = a Vertex3 b1 b2 b3 = b Vertex3 c1 c2 c3 = c Vertex3 d1 d2 d3 = d x = abs (a1 + b1 + c1 + d1) y = abs (a2 + b2 + c2 + d2) z = abs (a3 + b3 + c3 + d3) m = x + y + z quadColor :: Quad -> Color4 GLfloat quadColor quad = Color4 (realToFrac x) (realToFrac y) (realToFrac z) 1 where (x,y,z) = colMeans quad  Then, in our display function: haskell ...... renderPrimitive Quads $mapM_ drawQuad surface swapBuffers where drawQuad ((v1,v2,v3,v4),norm) = do materialDiffuse FrontAndBack$= quadColor (v1,v2,v3,v4) normal norm vertex v1 vertex v2 vertex v3 vertex v4  And this is the result: ![](figures/OpenGLcoloredCone.png) The full code is available in [this Github repo](https://github.com/stla/openglParametricPlot).