diff --git a/00_Introduction.lhs b/00_Introduction.lhs index 1075c68..4c130ac 100644 --- a/00_Introduction.lhs +++ b/00_Introduction.lhs @@ -39,7 +39,9 @@ and something nice to see in 3D. Here is what you'll end with: -blogimage("GoldenMandelbulb.png","A golden mandelbulb") +blogfigure("GoldenMandelbulb.png","The entire Mandelbulb") +blogfigure("3DMandelbulbDetail.png","A Mandelbulb detail") +blogfigure("3DMandelbulbDetail2.png","Another detail of the Mandelbulb") And here are the intermediate steps: diff --git a/01_Introduction/hglmandel.lhs b/01_Introduction/hglmandel.lhs index 41b6531..f8f5d69 100644 --- a/01_Introduction/hglmandel.lhs +++ b/01_Introduction/hglmandel.lhs @@ -140,7 +140,7 @@ Given two coordinates in pixels, it returns some integer value: > f (complex r i) 0 64 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. +The Mandelbrot set is the set of complex number \\(c\\) such that the following sequence does not escape to infinity. Let us define \\(f_c: \mathbb{C} \to \mathbb{C}\\) diff --git a/02_Edges/HGLMandelEdge.lhs b/02_Edges/HGLMandelEdge.lhs index b56a1f9..f292b99 100644 --- a/02_Edges/HGLMandelEdge.lhs +++ b/02_Edges/HGLMandelEdge.lhs @@ -55,7 +55,7 @@ The method I use is a rough approximation. I consider the Mandelbrot set to be almost convex. The result will be good enough. -We change slightly the drawMandelbrot function. +We change slightly the `drawMandelbrot` function. We replace the `Points` by `LineLoop` > drawMandelbrot = @@ -77,12 +77,12 @@ we will choose only point on the surface. > map (\(x,y,c) -> (x,-y,c)) (reverse positivePoints) We only need to compute the positive point. -The Mandelbrot set is symmetric on the abscisse axis. +The Mandelbrot set is symmetric relatively to the abscisse axis. > positivePoints :: [(GLfloat,GLfloat,Color3 GLfloat)] > positivePoints = do > x <- [-width..width] -> let y = findMaxOrdFor (mandel x) 0 height (log2 height) +> 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) @@ -92,22 +92,30 @@ The Mandelbrot set is symmetric on the abscisse axis. This function is interesting. 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 use a simple dichotomy: -> findMaxOrdFor func minval maxval 0 = (minval+maxval)/2 -> findMaxOrdFor func minval maxval n = +> -- 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 findMaxOrdFor func minval medpoint (n-1) -> else findMaxOrdFor func medpoint maxval (n-1) +> 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: diff --git a/03_Mandelbulb/Mandelbulb.lhs b/03_Mandelbulb/Mandelbulb.lhs index 59af425..ad90261 100644 --- a/03_Mandelbulb/Mandelbulb.lhs +++ b/03_Mandelbulb/Mandelbulb.lhs @@ -10,7 +10,7 @@ But it will be enough for us to create something that look nice. This section is quite long, but don't be afraid, most of the code is some OpenGL boilerplate. -For those you want to skim, +If you just want to skim this section, here is a high level representation: > - OpenGL Boilerplate @@ -263,7 +263,7 @@ depthPoints = do x <- [-width..width] y <- [-height..height] let - depthOf x' y' = findMaxOrdFor (mandel x' y') 0 deep logdeep + 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 @@ -292,7 +292,7 @@ Here is a harder to read but shorter and more generic rewritten function: > 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 +> 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 @@ -324,11 +324,21 @@ The rest of the program is very close to the preceding one.
-> findMaxOrdFor func minval maxval 0 = (minval+maxval)/2 -> findMaxOrdFor func minval maxval n = +> -- 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 func minval maxval 0 = (minval+maxval)/2 +> maxZeroIndex func minval maxval n = > if (func medpoint) /= 0 -> then findMaxOrdFor func minval medpoint (n-1) -> else findMaxOrdFor func medpoint maxval (n-1) +> then maxZeroIndex func minval medpoint (n-1) +> else maxZeroIndex func medpoint maxval (n-1) > where medpoint = (minval+maxval)/2 I made the color slightly brighter diff --git a/04_Mandelbulb/Mandelbulb.lhs b/04_Mandelbulb/Mandelbulb.lhs index 270eab6..9bad8ec 100644 --- a/04_Mandelbulb/Mandelbulb.lhs +++ b/04_Mandelbulb/Mandelbulb.lhs @@ -42,7 +42,7 @@ This is similar to the preceding section. > 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 +> 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 @@ -56,11 +56,20 @@ This is similar to the preceding section. > -- 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 = +> +> -- 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 findMaxOrdFor func minval medpoint (n-1) -> else findMaxOrdFor func medpoint maxval (n-1) +> then maxZeroIndex func minval medpoint (n-1) +> else maxZeroIndex func medpoint maxval (n-1) > where medpoint = (minval+maxval)/2 > > colorFromValue n = diff --git a/05_Mandelbulb/Mandelbulb.lhs b/05_Mandelbulb/Mandelbulb.lhs index 08fe85c..a480b41 100644 --- a/05_Mandelbulb/Mandelbulb.lhs +++ b/05_Mandelbulb/Mandelbulb.lhs @@ -2,21 +2,23 @@ Some points: -1. OpenGL and GLUT are linked to the C library. +1. OpenGL and GLUT is done in C. In particular the `mainLoop` function is a direct link to the C library (FFI). - This function if so far from the pure spirit of functional languages. + This function is clearly far from the functional paradigm. Could we make this better? We will have two choices, or create our own `mainLoop` function to make it more functional. Or deal with the imperative nature of the GLUT `mainLoop` function. - As a goal of this article is to understand how to deal with existing library and particularly the one coming from imperative language we will continue to use the `mainLoop` function. -2. Or main problem come from user interaction. - If you ask the Internet, about how to deal with user interaction with a functional paradigm, the main answer is to use _functional reactive programming_ (FRP). - I read very few about FRP, and I might be completely wrong when I say that it is about creating a DSL where atoms are time functions. - While I'm writing these lines, I don't know if I'll do something looking close to that. - For now I'll simply try to resolve the first problem. + As one of the goal of this article is to understand how to deal with existing libraries and particularly the one coming from imperative languages, we will continue to use the `mainLoop` function. +2. Our main problem come from user interaction. + If you ask "the Internet", + about how to deal with user interaction with a functional paradigm, + the main answer is to use _functional reactive programming_ (FRP). + I won't use FRP in this article. + Instead, I'll use a simpler while less effective way to deal with user interaction. + But The method I'll use will be as pure and functional as possible. -Then here is how I imagine things should go. -First, what the main loop should look like: +Here is how I imagine things should go. +First, what the main loop should look like if we could make our own: functionalMainLoop = @@ -89,6 +91,14 @@ We simply have to provide the definition of some functions. > 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 @@ -110,8 +120,8 @@ These function are used to update the world state. > zdir :: Point3D > zdir = makePoint3D (0,0,1) -Note `(-*<)` is scalar product. -Also note we could add Point3D as numbers. +Note `(-*<)` is the scalar product (`α -*< (x,y,z) = (αx,αy,αz)`). +Also note we could add two Point3D. > rotate :: Point3D -> Scalar -> World -> World > rotate dir angleValue world = @@ -188,9 +198,9 @@ Because we consider partial functions > shapeFunc :: Scalar -> Function3D > shapeFunc res x y = > let -> z = findMaxOrdFor (ymandel x y) 0 1 20 +> z = maxZeroIndex (ymandel x y) 0 1 20 > in -> if and [ findMaxOrdFor (ymandel (x+xeps) (y+yeps)) 0 1 20 < 0.000001 | +> 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)) @@ -207,13 +217,21 @@ With the color function. The rest is similar to the preceding sections. -> findMaxOrdFor :: (Fractional a,Num a,Num b,Eq b) => +> -- 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 -> 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) +> 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 > > ymandel :: Point -> Point -> Point -> Point diff --git a/06_Mandelbulb/Mandelbulb.lhs b/06_Mandelbulb/Mandelbulb.lhs index d8d7c97..1b05472 100644 --- a/06_Mandelbulb/Mandelbulb.lhs +++ b/06_Mandelbulb/Mandelbulb.lhs @@ -1,7 +1,8 @@ ## Optimization From the architecture stand point all is clear. -If you read the code, you'll see I didn't made everything perfect, for example, I didn't coded nicely the lights. +If you read the code of `YGL.hs`, you'll see I didn't made everything perfect. +For example, I didn't finished the code of the lights. But I believe it is a good first step and it will be easy to go further. The separation between rendering and world behavior is clear. Unfortunately the program of the preceding session is extremely slow. @@ -124,13 +125,14 @@ Our initial world state is slightly changed: > , 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 } +> , 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 We use the `YGL.getObject3DFromShapeFunction` function directly. This way instead of providing `XYFunc`, we provide directly a list of Atoms. @@ -183,9 +185,9 @@ All the rest is exactly the same. > shapeFunc :: Scalar -> Function3D > shapeFunc res x y = > let -> z = findMaxOrdFor (ymandel x y) 0 1 20 +> z = maxZeroIndex (ymandel x y) 0 1 20 > in -> if and [ findMaxOrdFor (ymandel (x+xeps) (y+yeps)) 0 1 20 < 0.000001 | +> 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 0) @@ -198,13 +200,21 @@ All the rest is exactly the same. > in > makeColor (t n) (t (n+5)) (t (n+10)) > -> findMaxOrdFor :: (Fractional a,Num a,Num b,Eq b) => +> -- 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 -> findMaxOrdFor _ minval maxval 0 = (minval+maxval)/2 -> findMaxOrdFor func minval maxval n = +> maxZeroIndex _ minval maxval 0 = (minval+maxval)/2 +> maxZeroIndex func minval maxval n = > if func medpoint /= 0 -> then findMaxOrdFor func minval medpoint (n-1) -> else findMaxOrdFor func medpoint maxval (n-1) +> then maxZeroIndex func minval medpoint (n-1) +> else maxZeroIndex func medpoint maxval (n-1) > where medpoint = (minval+maxval)/2 > > ymandel :: Point -> Point -> Point -> Point @@ -212,6 +222,8 @@ All the rest is exactly the same.
+And you can also consider small changes in other source files. + - [`YGL.hs`](code/06_Mandelbulb/YGL.hs), the 3D rendering framework - [`Mandel`](code/06_Mandelbulb/Mandel.hs), the mandel function - [`ExtComplex`](code/06_Mandelbulb/ExtComplex.hs), the extended complexes diff --git a/07_Conclusion/Conclusion.lhs b/07_Conclusion/Conclusion.lhs new file mode 100644 index 0000000..96797c8 --- /dev/null +++ b/07_Conclusion/Conclusion.lhs @@ -0,0 +1,19 @@ + ## Conclusion + +In conclusion, if you add the courage to read this long article, +you should see how to organize your code in a functional way. + +As we can use imperative style in a functional language, +know you can use functional style in imperative languages. + +The usage of Haskell made it very simple to use our own data type. +The gigantic advantage is to code while being perfectly pure. + +In the two last sections, the code is completely pure and functional. +Furthermore I don't use `GLfloat`, `Color3` or any other OpenGL type. +Therefore, if tomorrow I want to use another library I could do this simply by updating the `YGL` module. + +The `YGL` module is a "wrapper". +It is a clean separator between the imperative paradigm and functional paradigm. + +Furthermore I demonstrated you can make user interaction and display 3D objects in real time using only a pure functional paradigm.