segment intersection search

Here's a simple method I learned from Gareth Rees for finding the intersection (or parallel / collinear status) of two line segments. Rees credits Ronald Goldman, though I'd imagine the technique goes further back.

	import Control.Arrow ((***))
	import Control.Monad (liftM2)
	import Data.Int (Int16)
	import Data.List (unfoldr)
	import Graphics.UI.SDL as SDL
	import Graphics.UI.SDL.Primitives (line, circle)
	import System.Random.Mersenne.Pure64 (pureMT, randomInt)
	(xres, yres) = (1600, 900)
	(red, white) = (0xFF0000FF, 0xFFFFFFFF)
	main = withInit [InitVideo] $ do
	w <- setVideoMode xres yres 32 []
	enableEvent SDLMouseMotion False
	setCaption "Line Intersection" "Line Intersection"
	run w . take 15 $ filter f segments
	where
	f pq = dist pq < 1000
	run w ss = do
	delay 256
	drawSegs w ss
	let ss' = map (vectorize . cst) ss
	drawPoints w $ liftM2 crossing ss' ss'
	SDL.flip w
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	_ -> run w ss
	drawSegs w = mapM_ f
	where
	f ((a,b), (c,d)) = line w a b c d $ Pixel red
	drawPoints w = mapM_ (f . (round *** round))
	where
	f (x,y) = circle w x y 5 $ Pixel white
	-------------------------------------------------------
	segments = h $ zip (f xres 13) $ f yres 23
	where
	f m = map (fromIntegral . (`mod` m)) . g . pureMT
	g = unfoldr $ fmap randomInt . Just
	h (p:ps) = (p, head ps) : h (tail ps)
	dist (p,q) = sqrt $ (a-c)^2 + (b-d)^2
	where
	[a,b,c,d] = map f [fst p, snd p, fst q, snd q]
	f = fromIntegral
	crossing (p,r) (q,s)
	| f = p .+. (t .^. r) -- exactly one intersection
	| True = (-5,-5) -- parallel, colinear, etc.
	where
	t = (q .-. p) .*. s / (r .*. s)
	u = (q .-. p) .*. r / (r .*. s)
	f = (r .*. s) /= 0 && h t && h u
	h n = 0 <= n && n <= 1
	(a,b) .+. (c,d) = (a+c, b+d)
	(a,b) .-. (c,d) = (a-c, b-d)
	n .^. (x,y) = (n*x, n*y) -- vector scaling
	(a,b) .*. (c,d) = a*d - b*c -- 2D cross product
	vectorize ((a,b), (c,d)) = ((a,b), (c-a, d-b))
	cst (p,q) = (f p, f q)
	where
	f = fromIntegral *** fromIntegral

view raw line_segment_intersections.hs hosted with ❤ by GitHub

today in oblique approaches

If your standard library offers complex numbers and you're not in any hurry, you can't ask for much simpler (or more obscure!) ways to compute Π than this.

	import Data.Complex (magnitude, Complex((:+)))
	pi = length $ takeWhile (< 2) $ map magnitude zs
	where
	zs = iterate (\z -> z^2 + ((-0.75) :+ 0.0000001)) 0

view raw mandelbrot_PI.hs hosted with ❤ by GitHub

ikeda map

Continuing the theme of strange attractors, here's the well-known one embedded in the Ikeda map. The thumbnail at left shows the central 'vortex' of the attractor, and links to a larger viewport.

In these images, I've plotted the real and imaginary components along the x and y axes respectively. But the more popular way to visualize this attractor adds an extra parameter to the system and is expressed in trigonometric functions. Such adaptation of the code below yields these results.

	{-# LANGUAGE NoMonomorphismRestriction #-}
	import Data.Complex
	import Graphics.UI.SDL as SDL
	(xRes, yRes) = (1366, 768)
	[a, b, k, p] = [0.85, 0.9, 0.4, 7.7]
	main = withInit [InitVideo] $ do
	w <- setVideoMode xRes yRes 32 [NoFrame]
	s <- createRGBSurface [] 1 1 32 0 0 0 0
	fillRect s Nothing $ Pixel 0xFFFFFF
	enableEvent SDLMouseMotion False
	setCaption "Ikeda" "Ikeda"
	render w s $ ikeda $ 0 :+ 0
	SDL.flip w
	run w
	run w = do
	e <- pollEvent
	delay 64
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	_ -> run w
	render w s = mapM_ draw
	where
	rect (x,y) = Just $ Rect (round x) (round y) 1 1
	g (x,y) = (x + xRes / 24, y + yRes / 2.25)
	draw z = blitSurface s Nothing w $ rect $ g p
	where
	p = (1050 * realPart z, 1050 * imagPart z)
	ikeda = take 1500000 . filter g . iterate f
	where
	f z = a + b * z * exp(i * (k-p) * (1 + abs z^2))
	i = 0 :+ 1
	g (r:+i) = -35 < r && r < 35 && -18 < i && i < 18

view raw ikeda-map.hs hosted with ❤ by GitHub

hénon attractor

French astronomer Michel Hénon reported on this strange, fractal attractor in 1976. Since then, it has been among the most studied examples of chaotic dynamical systems.

	{-# LANGUAGE NoMonomorphismRestriction #-}
	import Graphics.UI.SDL as SDL
	(xRes, yRes) = (1366, 768)
	(a, b) = (1.4, 0.3)
	main = withInit [InitVideo] $ do
	w <- setVideoMode xRes yRes 32 [NoFrame]
	s <- createRGBSurface [] 1 1 32 0 0 0 0
	fillRect s Nothing $ Pixel 0xFFFFFF
	enableEvent SDLMouseMotion False
	setCaption "Hénon" "Hénon"
	render w s henon
	SDL.flip w
	run w
	run w = do
	e <- pollEvent
	delay 64
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	_ -> run w
	render w s = mapM_ draw
	where
	rect (x,y) = Just $ Rect (round x) (round y) 1 1
	g (x,y) = (500*x + xRes / 2, 900*y + yRes / 2)
	draw (x,y) = blitSurface s Nothing w $ rect $ g (x,y)
	henon = take 99999 $ filter g $ iterate f (0,0)
	where
	f (x,y) = (1 - a*x^2 + y, b*x)
	g (x,y) = -1.5 < x && x < 1.5 && -0.45 < y && y < 0.45

view raw michel_henon.hs hosted with ❤ by GitHub

langton's ant

For a round of code golf, I wrote this spare implementation of Chris Langton's remarkably simple universal computer. If you want amenities like pause, random starting pattern or even quit, check out this more complete version.

	import Control.Arrow ((***))
	import Control.Monad (when, join)
	import Data.Set (insert, delete, member, empty)
	import Graphics.UI.SDL as SDL
	(xres, yres, sq) = (1366, 768, 4)
	main = withInit [InitVideo] $ do
	w <- setVideoMode xres yres 32 [NoFrame]
	enableEvent SDLMouseMotion False
	setCaption "Langton's Ant" "Langton's Ant"
	run w [1..] p (0,1) empty
	where
	p = (xres `div` 2 `div` sq, yres `div` 2 `div` sq)
	run w (n:ns) p v ps = do
	when (n `mod` 7 == 0) render
	run w ns (move p $ g v) (g v) $ f p ps
	where
	b = member p ps
	(f,g) = if b then (delete, fl) else (insert, fr)
	fr (x,y) = if x == 0 then (-y,x) else (y,x)
	fl (x,y) = if x == 0 then (y,x) else (y,-x)
	move (x,y) = (x+) *** (+y)
	render = do
	fillRect w (Just $ Rect 0 0 xres yres) $ Pixel 0
	mapM_ (draw w . join (***) (* sq)) ps
	SDL.flip w
	draw w p = f p =<< g [SWSurface] sq sq 32 0 0 0 0
	where
	rect x y = Just $ Rect x y sq sq
	g = createRGBSurface
	f (x,y) s = do fillRect s (rect 0 0) $ Pixel 0xFFFFFF
	blitSurface s (rect 0 0) w $ rect x y

view raw langton_ant.hs hosted with ❤ by GitHub