diff --git a/01_Introduction/hglmandel.lhs b/01_Introduction/hglmandel.lhs index a33961d..ba87e20 100644 --- a/01_Introduction/hglmandel.lhs +++ b/01_Introduction/hglmandel.lhs @@ -1,19 +1,26 @@ ## First version -Let us start with the general architecture. +I splitted even the first version in two parts. +The first being mostly some boilerplate. +Generally in Haskell you need to declare a lot of import lines. +This is something I find annoying. +In particular, it should be possible to create a special file, Import.hs +which make all the necessary import for you, as you generally need them all. +I understand why this is cleaner to force the programmer not to do so, but, each time I do a copy/paste, I feel something is wrong. +The second part, contain more interresting stuff. +Even in this part, there are some necessary boilerplate. +But it is due to the OpenGL library this time. -First, the import system of Haskell is full of boilerplate. -A bit like Java. Let's play the song of our people: + ### Let's play the song of our people > import Graphics.Rendering.OpenGL > import Graphics.UI.GLUT > import Data.IORef -Also, for efficiency reason, I won't use the default Haskell `Complex` data type. We declare our type: +For efficiency reason, I won't use the default Haskell `Complex` data type. > newtype Complex = C (Float,Float) deriving (Show,Eq) -and make it an element of the typeclass `Num`: > instance Num Complex where > fromInteger n = C (fromIntegral n,0.0) @@ -36,6 +43,9 @@ We declare some useful functions for manipulating complex numbers: > magnitude :: Complex -> Float > magnitude = real.abs + + ### Let us start + Well, up until here we didn't made something useful. Just a lot of boilerplate and default value. Sorry but it is not completely the end. @@ -130,20 +140,16 @@ Given two coordinates in pixels, it returns some integer value: > in > f (complex r i) 0 64 -It uses the main mandelbrot function for each complex c. +It uses the main mandelbrot function for each complex $c$. The mandelbrot set is the set of complex number c such that the following sequence does not escape to infinity. -Let us define: +Let us define $f_c: \mathbb{C} \to \mathbb{C}$ -
-$$ f_c : z \rightarrow z^2 + c $$ -
+$$ f_c(z) = z^2 + c $$ The serie is: -
$$ 0 \rightarrow f_c(0) \rightarrow f_c(f_c(0)) \rightarrow \cdots \rightarrow f^n_c(0) \rightarrow \cdots $$ -
Of course, instead of trying to test the real limit, we just make a test after a finite number of occurences. @@ -171,4 +177,6 @@ See what occurs if we make the window bigger: blogimage("hglmandel_v01_too_wide.png","The mandelbrot too wide, black lines and columns") Wep, we see some black lines. -Why? Simply because we drawed points and not surfaces. +Why? Simply because we drawed less point than there is on the surface. +We can repair this by drawing little squares instead of just points. +But, instead we will do something a bit different and unusual. diff --git a/03_Mandelbulb/Mandelbulb.lhs b/03_Mandelbulb/Mandelbulb.lhs new file mode 100644 index 0000000..c1298ce --- /dev/null +++ b/03_Mandelbulb/Mandelbulb.lhs @@ -0,0 +1,180 @@ + ## 3D Mandelbrot? + +Why only draw the edge? +It is clearly not as nice as drawing the complete surface. +Yeah, I know, but, as we use OpenGL, why not show something in 3D. + +But, complex number are only in 2D and there is no 3D equivalent to complex. +In fact, the only extension known are quaternions, 4D. +As I know almost nothing about quaternions, I will use some extended complex. +I am pretty sure this construction is not useful for numbers. +But it will be enough for us to create something nice. + +
+ +> import Graphics.Rendering.OpenGL +> import Graphics.UI.GLUT +> import Data.IORef + +> type ColoredPoint = (GLfloat,GLfloat,GLfloat,Color3 GLfloat) + +
+ +> newtype ExtComplex = C (Float,Float,Float) deriving (Show,Eq) +> instance Num ExtComplex where +> fromInteger n = C (fromIntegral n,0.0,0.0) +> C (x,y,u) * C (z,t,v) = C (z*x - y*t, y*z + x*t, 0.0) +> C (x,y,u) + C (z,t,v) = C (x+z, y+t, u+v) +> abs (C (x,y,z)) = C (sqrt (x*x + y*y + z*z),0.0,0.0) +> signum (C (x,y,z)) = C (signum x , 0.0, 0.0) +> extcomplex :: Float -> Float -> Float -> ExtComplex +> extcomplex x y z = C (x,y,z) +> +> real :: ExtComplex -> Float +> real (C (x,y,z)) = x +> +> im :: ExtComplex -> Float +> im (C (x,y,z)) = y +> +> strange :: ExtComplex -> Float +> strange (C (x,y,z)) = z +> +> magnitude :: ExtComplex -> 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 +> +> width = 32 :: GLfloat +> height = 32 :: GLfloat +> deep = 32 :: GLfloat + + + + +This time, instead of drawing all points, I'll simply want to draw the edges of the Mandelbrot set. +We change slightly the drawMandelbrot function. +We replace the `Points` by `LineLoop` + +> drawMandelbrot = +> -- We will print Points (not triangles for example) +> renderPrimitive TriangleStrip $ do +> mapM_ drawColoredPoint allPoints +> where +> drawColoredPoint (x,y,z,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 z + +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 :: [ColoredPoint] +> allPoints = order $ positivePoints ++ map (\(x,y,z,c) -> (x,-y,z,c)) (reverse positivePoints) + +> closestPoint :: ColoredPoint -> [ColoredPoint] -> ColoredPoint +> closestPoint x xs = clPoint x xs (width,height,deep,colorFromValue 0) +> where +> clPoint :: ColoredPoint -> [ColoredPoint] -> ColoredPoint -> ColoredPoint +> clPoint _ [] res = res +> clPoint p (p':xs) res = +> if dist p p' < dist p res +> then clPoint p xs p' +> else clPoint p xs res + +> dist (x,y,z,_) (x',y',z',_) = x*x' + y*y' + z*z' + +> order :: [ColoredPoint] -> [ColoredPoint] +> order [] = [] +> order (x:[]) = [x] +> order (x:xs) = +> let +> closest = closestPoint x xs +> newlist = filter (/=closest) xs +> in x:order (closest:newlist) + + +We only need to compute the positive point. +The mandelbrot set is symetric on the abscisse axis. + +> positivePoints :: [ColoredPoint] +> positivePoints = do +> x <- [-width..width] +> y <- [0..height] +> let z = findMaxOrdFor (mandel x y) 0 deep 6 -- log deep +> if z < 1 +> then [] +> else return (x/width,y/height,z/deep,colorFromValue $ mandel x y z) + + > else [(x/width,y/height,z/deep,colorFromValue $ mandel x y z), + > (x/width,y+1/height,z/deep,colorFromValue $ mandel x y z), + > (x+1/width,y+1/height,z/deep,colorFromValue $ mandel x y z), + > (x+1/width,y/height,z/deep,colorFromValue $ mandel x y z) ] + + +This function is interresting. +For those not used to the list monad here is a natural language version of this function: + +~~~ +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 create a simple dichotomic search: + +> 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 + +The new mandel function + +> 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 + + +
+ +> 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)) + +> f :: ExtComplex -> ExtComplex -> 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) + +
+ diff --git a/create_ymd.sh b/create_ymd.sh index 4855d17..96fdc3b 100755 --- a/create_ymd.sh +++ b/create_ymd.sh @@ -68,4 +68,4 @@ for fic in **/*.lhs(.N); do cat $fic ((contains_haskell)) && \ print -- "\n${fic:h}/${fic:t} \n" -done | perl -pe 'BEGIN{$/="";} s#((^>.*\n)+)#
\n\n$1\n
#mg' | perl -pe 's#^> ?##' | perl -pe 's/^ #/#/' +done | perl -pe 'BEGIN{$/="";} s#((^>.*\n)+)#
\n\n$1\n
\n#mg' | perl -pe 's#^> ?##' | perl -pe 's/^ #/#/'