----- isHidden: false menupriority: 1 kind: article created_at: 2012-04-30T19:17:53+02:00 title: Haskell Working Program subtitle: An OpenGL 3D extension of the Mandelbrot set author_name: Yann Esposito author_uri: yannesposito.com tags: - Haskell - programming - functional - tutorial ----- blogimage("BenoitBMandelbrot.jpg","The B in Benoît B. Mandelbrot stand for Benoît B. Mandelbrot") begindiv(intro) %tldr You will see how to go from theory to a real application using Haskell. >
import Graphics.Rendering.OpenGL
import Graphics.UI.GLUT
import Data.IORef
newtype Complex = C (Float,Float) deriving (Show,Eq)
instance Num Complex where
fromInteger n = C (fromIntegral n,0.0)
C (x,y) * C (z,t) = C (z*x - y*t, y*z + x*t)
C (x,y) + C (z,t) = C (x+z, y+t)
abs (C (x,y)) = C (sqrt (x*x + y*y),0.0)
signum (C (x,y)) = C (signum x , 0.0)
complex :: Float -> Float -> Complex
complex x y = C (x,y)
real :: Complex -> Float
real (C (x,y)) = x
im :: Complex -> Float
im (C (x,y)) = y
magnitude :: Complex -> Float
magnitude = real.abs
main :: IO ()
main = do
-- GLUT need to be initialized
(progname,_) <- getArgsAndInitialize
-- We will use the double buffered mode (GL constraint)
initialDisplayMode $= [DoubleBuffered]
-- We create a window with some title
createWindow "Mandelbrot Set with Haskell and OpenGL"
-- Each time we will need to update the display
-- we will call the function 'display'
displayCallback $= display
-- We enter the main loop
mainLoop
display = do
clear [ColorBuffer] -- make the window black
loadIdentity -- reset any transformation
preservingMatrix drawMandelbrot
swapBuffers -- refresh screen
drawMandelbrot =
-- We will print Points (not triangles for example)
renderPrimitive Points $ do
mapM_ drawColoredPoint allPoints
where
drawColoredPoint (x,y,c) = do
color c -- set the current color to c
-- then draw the point at position (x,y,0)
-- remember we're in 3D
vertex $ Vertex3 x y 0
width = 320 :: GLfloat
height = 320 :: GLfloat
allPoints :: [(GLfloat,GLfloat,Color3 GLfloat)]
allPoints = [ (x/width,y/height,colorFromValue $ mandel x y) |
x <- [-width..width],
y <- [-height..height]]
colorFromValue n =
let
t :: Int -> GLfloat
t i = 0.5 + 0.5*cos( fromIntegral i / 10 )
in
Color3 (t n) (t (n+5)) (t (n+10))
mandel x y =
let r = 2.0 * x / width
i = 2.0 * y / height
in
f (complex r i) 0 64
f :: Complex -> Complex -> Int -> Int
f c z 0 = 0
f c z n = if (magnitude z > 2 )
then n
else f c ((z*z)+c) (n-1)
drawMandelbrot =
-- We will print Points (not triangles for example)
renderPrimitive LineLoop $ do
mapM_ drawColoredPoint allPoints
where
drawColoredPoint (x,y,c) = do
color c -- set the current color to c
-- then draw the point at position (x,y,0)
-- remember we're in 3D
vertex $ Vertex3 x y 0
allPoints = positivePoints ++
map (\(x,y,c) -> (x,-y,c)) (reverse positivePoints)
positivePoints :: [(GLfloat,GLfloat,Color3 GLfloat)]
positivePoints = do
x <- [-width..width]
let y = findMaxOrdFor (mandel x) 0 height (log2 height)
if y < 1 -- We don't draw point in the absciss
then []
else return (x/width,y/height,colorFromValue $ mandel x y)
where
log2 n = floor ((log n) / log 2)
findMaxOrdFor func minval maxval 0 = (minval+maxval)/2
findMaxOrdFor func minval maxval n =
if (func medpoint) /= 0
then findMaxOrdFor func minval medpoint (n-1)
else findMaxOrdFor func medpoint maxval (n-1)
where medpoint = (minval+maxval)/2
data ExtComplex = C (GLfloat,GLfloat,GLfloat)
deriving (Show,Eq)
instance Num ExtComplex where
-- The shape of the 3D mandelbrot
-- will depend on this formula
C (x,y,z) * C (x',y',z') = C (x*x' - y*y' - z*z',
x*y' + y*x' + z*z',
x*z' + z*x' )
-- The rest is straightforward
fromInteger n = C (fromIntegral n, 0, 0)
C (x,y,z) + C (x',y',z') = C (x+x', y+y', z+z')
abs (C (x,y,z)) = C (sqrt (x*x + y*y + z*z), 0, 0)
signum (C (x,y,z)) = C (signum x, 0, 0)
main :: IO ()
main = do
-- GLUT need to be initialized
(progname,_) <- getArgsAndInitialize
-- We will use the double buffered mode (GL constraint)
-- We also Add the DepthBuffer (for 3D)
initialDisplayMode $=
[WithDepthBuffer,DoubleBuffered,RGBMode]
-- We create a window with some title
createWindow "3D HOpengGL Mandelbrot"
-- We add some directives
depthFunc $= Just Less
windowSize $= Size 500 500
-- Some state variables (I know it feels BAD)
angle <- newIORef ((35,0)::(GLfloat,GLfloat))
zoom <- newIORef (2::GLfloat)
campos <- newIORef ((0.7,0)::(GLfloat,GLfloat))
-- Function to call each frame
idleCallback $= Just idle
-- Function to call when keyboard or mouse is used
keyboardMouseCallback $=
Just (keyboardMouse angle zoom campos)
-- Each time we will need to update the display
-- we will call the function 'display'
-- But this time, we add some parameters
displayCallback $= display angle zoom campos
-- We enter the main loop
mainLoop
idle = postRedisplay Nothing
modVar v f = do
v' <- get v
v $= (f v')
mapFst f (x,y) = (f x, y)
mapSnd f (x,y) = ( x,f y)
keyboardMouse angle zoom campos key state modifiers position =
-- We won't use modifiers nor position
kact angle zoom campos key state
where
-- reset view when hitting space
kact a z p (Char ' ') Down = do
a $= (0,0) -- angle
z $= 1 -- zoom
p $= (0,0) -- camera position
-- use of hjkl to rotate
kact a _ _ (Char 'h') Down = modVar a (mapFst (+0.5))
kact a _ _ (Char 'l') Down = modVar a (mapFst (+(-0.5)))
kact a _ _ (Char 'j') Down = modVar a (mapSnd (+0.5))
kact a _ _ (Char 'k') Down = modVar a (mapSnd (+(-0.5)))
-- use o and i to zoom
kact _ z _ (Char 'o') Down = modVar z (*1.1)
kact _ z _ (Char 'i') Down = modVar z (*0.9)
-- use sdfe to move the camera
kact _ _ p (Char 's') Down = modVar p (mapFst (+0.1))
kact _ _ p (Char 'f') Down = modVar p (mapFst (+(-0.1)))
kact _ _ p (Char 'd') Down = modVar p (mapSnd (+0.1))
kact _ _ p (Char 'e') Down = modVar p (mapSnd (+(-0.1)))
-- any other keys does nothing
kact _ _ _ _ _ = return ()
display angle zoom position = do
-- set the background color (dark solarized theme)
clearColor $= Color4 0 0.1686 0.2117 1
clear [ColorBuffer,DepthBuffer]
-- Transformation to change the view
loadIdentity -- reset any transformation
-- tranlate
(x,y) <- get position
translate $ Vector3 x y 0
-- zoom
z <- get zoom
scale z z z
-- rotate
(xangle,yangle) <- get angle
rotate xangle $ Vector3 1.0 0.0 (0.0::GLfloat)
rotate yangle $ Vector3 0.0 1.0 (0.0::GLfloat)
-- Now that all transformation were made
-- We create the object(s)
preservingMatrix drawMandelbrot
swapBuffers -- refresh screen
nbDetails = 200 :: GLfloat
width = nbDetails
height = nbDetails
deep = nbDetails
drawMandelbrot = do
-- We will print Points (not triangles for example)
renderPrimitive Triangles $ do
mapM_ drawColoredPoint allPoints
where
drawColoredPoint (x,y,z,c) = do
color c
vertex $ Vertex3 x y z
depthPoints :: [ColoredPoint]
depthPoints = do
x <- [-width..width]
y <- [-height..height]
let
depthOf x' y' = findMaxOrdFor (mandel x' y') 0 deep logdeep
logdeep = floor ((log deep) / log 2)
z1 = depthOf x y
z2 = depthOf (x+1) y
z3 = depthOf (x+1) (y+1)
z4 = depthOf x (y+1)
c1 = mandel x y (z1+1)
c2 = mandel (x+1) y (z2+1)
c3 = mandel (x+1) (y+1) (z3+1)
c4 = mandel x (y+1) (z4+1)
p1 = ( x /width, y /height, z1/deep, colorFromValue c1)
p2 = ((x+1)/width, y /height, z2/deep, colorFromValue c2)
p3 = ((x+1)/width,(y+1)/height, z3/deep, colorFromValue c3)
p4 = ( x /width,(y+1)/height, z4/deep, colorFromValue c4)
if (and $ map (>=57) [c1,c2,c3,c4])
then []
else [p1,p2,p3,p1,p3,p4]
If you look at the function above, you see a lot of common patterns.
Haskell is very efficient to make this better.
Here is a harder to read but shorter and more generic rewritten function:
depthPoints :: [ColoredPoint]
depthPoints = do
x <- [-width..width]
y <- [-height..height]
let
neighbors = [(x,y),(x+1,y),(x+1,y+1),(x,y+1)]
depthOf (u,v) = findMaxOrdFor (mandel u v) 0 deep logdeep
logdeep = floor ((log deep) / log 2)
-- zs are 3D points with found depth
zs = map (\(u,v) -> (u,v,depthOf (u,v))) neighbors
-- ts are 3D pixels + mandel value
ts = map (\(u,v,w) -> (u,v,w,mandel u v (w+1))) zs
-- ps are 3D opengl points + color value
ps = map (\(u,v,w,c') ->
(u/width,v/height,w/deep,colorFromValue c')) ts
-- If the point diverged too fast, don't display it
if (and $ map (\(_,_,_,c) -> c>=57) ts)
then []
-- Draw two triangles
else [ps!!0,ps!!1,ps!!2,ps!!0,ps!!2,ps!!3]
allPoints :: [ColoredPoint]
allPoints = planPoints ++ map inverseDepth planPoints
where
planPoints = depthPoints
inverseDepth (x,y,z,c) = (x,y,-z+1/deep,c)
mandel x y z =
let r = 2.0 * x / width
i = 2.0 * y / height
s = 2.0 * z / deep
in
f (extcomplex r i s) 0 64
import YBoiler -- Most the OpenGL Boilerplate
import Mandel -- The 3D Mandelbrot maths
main :: IO ()
main = yMainLoop "3D Mandelbrot" (\_ -> allPoints)
nbDetails = 200 :: GLfloat
width = nbDetails
height = nbDetails
deep = nbDetails
allPoints :: [ColoredPoint]
allPoints = planPoints ++ map inverseDepth planPoints
where
planPoints = depthPoints ++ map inverseHeight depthPoints
inverseHeight (x,y,z,c) = (x,-y,z,c)
inverseDepth (x,y,z,c) = (x,y,-z+1/deep,c)
depthPoints :: [ColoredPoint]
depthPoints = do
x <- [-width..width]
y <- [0..height]
let
neighbors = [(x,y),(x+1,y),(x+1,y+1),(x,y+1)]
depthOf (u,v) = findMaxOrdFor (ymandel u v) 0 deep 7
-- zs are 3D points with found depth
zs = map (\(u,v) -> (u,v,depthOf (u,v))) neighbors
-- ts are 3D pixels + mandel value
ts = map (\(u,v,w) -> (u,v,w,ymandel u v (w+1))) zs
-- ps are 3D opengl points + color value
ps = map (\(u,v,w,c') ->
(u/width,v/height,w/deep,colorFromValue c')) ts
-- If the point diverged too fast, don't display it
if (and $ map (\(_,_,_,c) -> c>=57) ts)
then []
-- Draw two triangles
else [ps!!0,ps!!1,ps!!2,ps!!0,ps!!2,ps!!3]
findMaxOrdFor func minval maxval 0 = (minval+maxval)/2
findMaxOrdFor func minval maxval n =
if (func medpoint) /= 0
then findMaxOrdFor func minval medpoint (n-1)
else findMaxOrdFor func medpoint maxval (n-1)
where medpoint = (minval+maxval)/2
colorFromValue n =
let
t :: Int -> GLfloat
t i = 0.7 + 0.3*cos( fromIntegral i / 10 )
in
((t n),(t (n+5)),(t (n+10)))
ymandel x y z = mandel (2*x/width) (2*y/height) (2*z/deep) 64
functionalMainLoop =
Read user inputs and provide a list of actions
Apply all actions to the World
Display one frame
repetere aeternum
Clearly, ideally we should provide only three parameters to this main loop function:
- an initial World state
- a mapping between the user interaction and function which modify the world
- a function taking two parameters: time and world state and render a new world without user interaction.
Here is a real working code, I've hidden most display functions.
The YGL, is a kind of framework to display 3D functions.
But it can easily be extended to many kind of representation.
import YGL -- Most the OpenGL Boilerplate
import Mandel -- The 3D Mandelbrot maths
-- Centralize all user input interaction
inputActionMap :: InputMap World
inputActionMap = inputMapFromList [
(Press 'k' , rotate xdir 5)
,(Press 'i' , rotate xdir (-5))
,(Press 'j' , rotate ydir 5)
,(Press 'l' , rotate ydir (-5))
,(Press 'o' , rotate zdir 5)
,(Press 'u' , rotate zdir (-5))
,(Press 'f' , translate xdir 0.1)
,(Press 's' , translate xdir (-0.1))
,(Press 'e' , translate ydir 0.1)
,(Press 'd' , translate ydir (-0.1))
,(Press 'z' , translate zdir 0.1)
,(Press 'r' , translate zdir (-0.1))
,(Press '+' , zoom 1.1)
,(Press '-' , zoom (1/1.1))
,(Press 'h' , resize 1.2)
,(Press 'g' , resize (1/1.2))
]
-- I prefer to set my own name for these types
data World = World {
angle :: Point3D
, scale :: Scalar
, position :: Point3D
, shape :: Scalar -> Function3D
, box :: Box3D
, told :: Time -- last frame time
}
instance DisplayableWorld World where
winTitle _ = "The YGL Mandelbulb"
camera w = Camera {
camPos = position w,
camDir = angle w,
camZoom = scale w }
objects w = [XYFunc ((shape w) res) defbox]
where
res = resolution $ box w
defbox = box w
xdir :: Point3D
xdir = makePoint3D (1,0,0)
ydir :: Point3D
ydir = makePoint3D (0,1,0)
zdir :: Point3D
zdir = makePoint3D (0,0,1)
rotate :: Point3D -> Scalar -> World -> World
rotate dir angleValue world =
world {
angle = (angle world) + (angleValue -*< dir) }
translate :: Point3D -> Scalar -> World -> World
translate dir len world =
world {
position = (position world) + (len -*< dir) }
zoom :: Scalar -> World -> World
zoom z world = world {
scale = z * scale world }
resize :: Scalar -> World -> World
resize r world = world {
box = (box world) {
resolution = sqrt ((resolution (box world))**2 * r) }}
main :: IO ()
main = yMainLoop inputActionMap idleAction initialWorld
-- We initialize the world state
-- then angle, position and zoom of the camera
-- And the shape function
initialWorld :: World
initialWorld = World {
angle = makePoint3D (-30,-30,0)
, position = makePoint3D (0,0,0)
, scale = 0.8
, shape = shapeFunc
, box = Box3D { minPoint = makePoint3D (-2,-2,-2)
, maxPoint = makePoint3D (2,2,2)
, resolution = 0.16 }
, told = 0
}
idleAction :: Time -> World -> World
idleAction tnew world = world {
angle = (angle world) + (delta -*< zdir)
, told = tnew
}
where
anglePerSec = 5.0
delta = anglePerSec * elapsed / 1000.0
elapsed = fromIntegral (tnew - (told world))
shapeFunc :: Scalar -> Function3D
shapeFunc res x y =
let
z = findMaxOrdFor (ymandel x y) 0 1 20
in
if and [ findMaxOrdFor (ymandel (x+xeps) (y+yeps)) 0 1 20 < 0.000001 |
val <- [res], xeps <- [-val,val], yeps<-[-val,val]]
then Nothing
else Just (z,colorFromValue ((ymandel x y z) * 64))
colorFromValue :: Point -> Color
colorFromValue n =
let
t :: Point -> Scalar
t i = 0.7 + 0.3*cos( i / 10 )
in
makeColor (t n) (t (n+5)) (t (n+10))
findMaxOrdFor :: (Fractional a,Num a,Num b,Eq b) =>
(a -> b) -> a -> a -> Int -> a
findMaxOrdFor _ minval maxval 0 = (minval+maxval)/2
findMaxOrdFor func minval maxval n =
if func medpoint /= 0
then findMaxOrdFor func minval medpoint (n-1)
else findMaxOrdFor func medpoint maxval (n-1)
where medpoint = (minval+maxval)/2
ymandel :: Point -> Point -> Point -> Point
ymandel x y z = fromIntegral (mandel x y z 64) / 64
Constant Function -> Constant List of Triangles -> Display
Now we have
World -> Function -> List of Objects -> Atoms -> Display
And the World state could change.
Then it is no more straightforward for the compiler to understand
when not to recompute the entire list of atoms.
Then to optimize we will have to make things a little less separate.
We must control the flow of atom generation.
Mostly the program is the same as before, but instead of providing a
function, we will provide the list of atoms directly.
data World = World {
angle :: Point3D
, anglePerSec :: Scalar
, scale :: Scalar
, position :: Point3D
, box :: Box3D
, told :: Time
-- We replace shape by cache
, cache :: [YObject]
}
instance DisplayableWorld World where
winTitle _ = "The YGL Mandelbulb"
camera w = Camera {
camPos = position w,
camDir = angle w,
camZoom = scale w }
-- We update our objects instanciation
objects = cache
-- We initialize the world state
-- then angle, position and zoom of the camera
-- And the shape function
initialWorld :: World
initialWorld = World {
angle = makePoint3D (30,30,0)
, anglePerSec = 5.0
, position = makePoint3D (0,0,0)
, scale = 1.0
, box = Box3D { minPoint = makePoint3D (-2,-2,-2)
, maxPoint = makePoint3D (2,2,2)
, resolution = 0.03 }
, told = 0
-- We declare cache directly this time
, cache = objectFunctionFromWorld initialWorld
}
objectFunctionFromWorld :: World -> [YObject]
objectFunctionFromWorld w = [Atoms atomList]
where atomListPositive =
getObject3DFromShapeFunction
(shapeFunc (resolution (box w))) (box w)
atomList = atomListPositive ++
map negativeTriangle atomListPositive
negativeTriangle (ColoredTriangle (p1,p2,p3,c)) =
ColoredTriangle (negz p1,negz p3,negz p2,c)
where negz (P (x,y,z)) = P (x,y,-z)
resize :: Scalar -> World -> World
resize r world =
tmpWorld { cache = objectFunctionFromWorld tmpWorld }
where
tmpWorld = world { box = (box world) {
resolution = sqrt ((resolution (box world))**2 * r) }}