152 lines
4.6 KiB
Text
152 lines
4.6 KiB
Text
## Only the edges
|
|
|
|
<div style="display:none">
|
|
|
|
> import Graphics.Rendering.OpenGL
|
|
> import Graphics.UI.GLUT
|
|
> import Data.IORef
|
|
> -- Use UNPACK data because it is faster
|
|
> -- The ! is for strict instead of lazy
|
|
> data Complex = C {-# UNPACK #-} !Float
|
|
> {-# UNPACK #-} !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
|
|
> -- set the background color (dark solarized theme)
|
|
> clearColor $= Color4 0 0.1686 0.2117 1
|
|
> clear [ColorBuffer] -- make the window black
|
|
> loadIdentity -- reset any transformation
|
|
> preservingMatrix drawMandelbrot
|
|
> swapBuffers -- refresh screen
|
|
>
|
|
> width = 320 :: GLfloat
|
|
> height = 320 :: GLfloat
|
|
|
|
|
|
</div>
|
|
|
|
This time, instead of drawing all points,
|
|
we will simply draw the edges of the Mandelbrot set.
|
|
The method I use is a rough approximation.
|
|
I consider the Mandelbrot set to be almost convex.
|
|
The result will be good enough for the purpose of this tutorial.
|
|
|
|
We change slightly the `drawMandelbrot` function.
|
|
We replace the `Points` by `LineLoop`
|
|
|
|
> 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
|
|
|
|
And now, we should change our list of points.
|
|
Instead of drawing every point of the visible surface,
|
|
we will choose only point on the surface.
|
|
|
|
> allPoints = positivePoints ++
|
|
> map (\(x,y,c) -> (x,-y,c)) (reverse positivePoints)
|
|
|
|
We only need to compute the positive point.
|
|
The Mandelbrot set is symmetric relatively to the abscisse axis.
|
|
|
|
> 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)
|
|
|
|
This function is interesting.
|
|
For those not used to the list monad here is a natural language version of this function:
|
|
|
|
<code class="no-highlight">
|
|
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))
|
|
</code>
|
|
|
|
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
|
|
|
|
No rocket science here. See the result now:
|
|
|
|
blogimage("HGLMandelEdges.png","The edges of the mandelbrot set")
|
|
|
|
<div style="display:none">
|
|
|
|
> 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)
|
|
|
|
</div>
|
|
|