----- isHidden: false menupriority: 1 kind: article created_at: 2012-06-15T19:16:43+02:00 title: Un example progressif avec Haskell subtitle: Une extension de l'ensemble de Mandelbrot en 3D et en OpenGL author_name: Yann Esposito author_uri: yannesposito.com tags: - Haskell - programming - functional - tutorial - fractal ----- blogimage("BenoitBMandelbrot.jpg","The B in Benoît B. Mandelbrot stand for Benoît B. Mandelbrot") begindiv(intro) %tlal Un exemple progressif d'utilisation d'Haskell. Vous pourrez voir un ensemble de Mandelbrot étendu à la troisième dimension. De plus le code sera très propre. Les détails de rendu sont séparés dans un module externe. Le code descriptif intéressant est concentré dans un environnement pur et fonctionnel. Vous pouvez vous inspirer de ce code utilisant le paradigme fonctional dans tous les languages. >
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 = maxZeroIndex (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)
positivePoints =
for all x in the range [-width..width]
let y be smallest number s.t. mandel x y > 0
if y is on 0 then don't return a point
else return the value corresonding to (x,y,color for (x+iy))
In fact using the list monad you write like if you consider only one element at a time and the computation is done non deterministically.
To find the smallest number such that `mandel x y > 0` we use a simple dichotomy:
-- given f min max nbtest,
-- considering
-- - f is an increasing function
-- - f(min)=0
-- - f(max)≠0
-- then maxZeroIndex f min max nbtest returns x such that
-- f(x - ε)=0 and f(x + ε)≠0
-- where ε=(max-min)/2^(nbtest+1)
maxZeroIndex func minval maxval 0 = (minval+maxval)/2
maxZeroIndex func minval maxval n =
if (func medpoint) /= 0
then maxZeroIndex func minval medpoint (n-1)
else maxZeroIndex 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' = maxZeroIndex (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) = maxZeroIndex (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) = maxZeroIndex (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]
-- given f min max nbtest,
-- considering
-- - f is an increasing function
-- - f(min)=0
-- - f(max)≠0
-- then maxZeroIndex f min max nbtest returns x such that
-- f(x - ε)=0 and f(x + ε)≠0
-- where ε=(max-min)/2^(nbtest+1)
maxZeroIndex func minval maxval 0 = (minval+maxval)/2
maxZeroIndex func minval maxval n =
if (func medpoint) /= 0
then maxZeroIndex func minval medpoint (n-1)
else maxZeroIndex 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 interactions and functions 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 for world w
-- is the list of one unique element
-- The element is an YObject
-- more precisely the XYFunc Function3D Box3D
-- where the Function3D is the type
-- Point -> Point -> Maybe (Point,Color)
-- and its value here is ((shape w) res)
-- and the Box3D value is defbox
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 = maxZeroIndex (ymandel x y) 0 1 20
in
if and [ maxZeroIndex (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))
-- given f min max nbtest,
-- considering
-- - f is an increasing function
-- - f(min)=0
-- - f(max)≠0
-- then maxZeroIndex f min max nbtest returns x such that
-- f(x - ε)=0 and f(x + ε)≠0
-- where ε=(max-min)/2^(nbtest+1)
maxZeroIndex :: (Fractional a,Num a,Num b,Eq b) =>
(a -> b) -> a -> a -> Int -> a
maxZeroIndex _ minval maxval 0 = (minval+maxval)/2
maxZeroIndex func minval maxval n =
if (func medpoint) /= 0
then maxZeroIndex func minval medpoint (n-1)
else maxZeroIndex 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
Main loop -> World -> Function -> List of Objects -> Atoms -> Display
The World state could change.
The compiler can no more optimize the computation for us.
We have to manually explain when to redraw the shape.
To optimize we must do some things in a lower level.
Mostly the program remains the same,
but it 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 (0-eps, 0-eps, 0-eps)
, maxPoint = makePoint3D (0+eps, 0+eps, 0+eps)
, resolution = 0.02 }
, told = 0
-- We declare cache directly this time
, cache = objectFunctionFromWorld initialWorld
}
where eps=2
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) }}