bézier curves

De Casteljau's algorithm is a fast, numerically stable way to rasterize Bézier curves. This code implements an interactive demo: click to append nodes to the curve's defining polygon, and drag any node to alter the curve.

	{-# LANGUAGE NoMonomorphismRestriction #-}
	import Control.Monad (when)
	import Control.Arrow ((***))
	import Data.List (findIndex)
	import Graphics.UI.SDL as SDL hiding (init)
	import Graphics.UI.SDL.Primitives (circle, line)
	dim = 400
	main = withInit [InitVideo] $ do
	w <- setVideoMode dim dim 32 []
	enableEvent SDLMouseMotion False
	setCaption "Bézier Curves" "Bézier Curves"
	loop w []
	plot w ps = do
	fillRect w (Just $ Rect 0 0 dim dim) $ Pixel 0xFF222255
	mapM_ (f 2 0xFFFFFFFF . zz) [head ps, last ps]
	when b $ mapM_ (f 3 0x888888FF . zz) controls
	when b $ mapM_ (f 2 0xBBBBBBFF . zz) controls
	where
	f r c (x,y) = circle w x y r $ Pixel c
	(b, controls) = (length ps > 2, tail $ init ps)
	limn w [_] = SDL.flip w
	limn w ((a,b):(x,y):ps) = do
	line w a b x y $ Pixel 0xFFFFFFFF
	limn w $ (x,y) : ps
	loop w ps = do
	delay 128
	event <- pollEvent
	case event of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	MouseButtonDown x y _ -> click x y
	_ -> loop w ps
	where
	click x y = let p = rr (x,y) in
	case findIndex ((10 >) . dist p) ps of
	Just i -> drag w i ps
	Nothing -> do
	let ps' = p : ps
	plot w ps' >> SDL.flip w
	when (length ps' > 2) $ render w ps'
	loop w ps'
	drag w i ps = do
	delay 16
	(x,y,_) <- getMouseState
	event <- pollEvent
	let ps' = swap i ps $ rr (x,y)
	plot w ps'
	when (length ps' > 2) $ render w ps'
	case event of
	MouseButtonUp x y _ -> loop w ps'
	_ -> drag w i ps'
	render w ps = limn w $ map zz curve
	where
	curve = map (casteljau ps) [0, 0.001.. 1]
	casteljau [p] t = p
	casteljau ps t = casteljau ps' t
	where
	ps' = zipWith (g t) ps $ tail ps
	g t (a,b) (c,d) = (f t a c, f t b d)
	f t a b = (1 - t) * a + t * b
	swap 0 ps p = p : tail ps
	swap i ps p = take i ps ++ p : drop (i+1) ps
	dist (a,b) (c,d) = sqrt $ (a-c)^2 + (b-d)^2
	rr = fromIntegral *** fromIntegral
	zz = round *** round

view raw interactive_bezier.hs hosted with ❤ by GitHub

convex hulls

This code demonstrates the Graham scan, an O(n log n) method for finding the convex hull (smallest enclosing polygon) of a planar point set. It's a great example of exploiting order: it works by sorting the set by angle about a known extreme point, allowing the hull points to be found in linear time.

A variation on the algorithm, noted by A.M. Andrew, sorts the set lexicographically and finds the upper and lower hull chains separately. This 'monotone chain' version is often preferred, since it's easier to do robustly.

	{-# LANGUAGE NoMonomorphismRestriction #-}
	import Control.Arrow ((***))
	import Data.List (maximumBy, delete, sort, sortBy, unfoldr)
	import Data.Ord (comparing)
	import Graphics.UI.SDL as SDL
	import Graphics.UI.SDL.Primitives (filledCircle, line)
	import System.Random.Mersenne.Pure64 (newPureMT, randomDouble)
	res = 250
	main = withInit [InitVideo] $ do
	w <- setVideoMode res res 32 []
	ps <- randPoints
	enableEvent SDLMouseMotion False
	setCaption "Graham Scan" "Graham Scan"
	fillRect w (Just $ Rect 0 0 res res) $ Pixel 0
	limn w $ map (round *** round) $ hull ps
	plot w ps
	pause
	plot w ps = do
	mapM_ (f . (round *** round)) ps
	SDL.flip w
	where
	f (x,y) = filledCircle w x y 1 $ Pixel 0xFFFFFFFF
	limn w ps = f $ ps ++ [head ps]
	where
	f [_] = return ()
	f ((a,b):(x,y):ps) = do
	line w a b x y $ Pixel 0xFF0000FF
	f $ (x,y) : ps
	pause = do
	delay 128
	e <- pollEvent
	case e of
	KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
	_ -> pause
	hull qs = go (drop 2 ps) $ reverse $ take 2 ps
	where
	o = bottomRightP qs
	ps = o : sortBy (ccw o) (delete o qs)
	go [] qs = qs
	go (p:ps) s@(a:b:qs)
	| ccw a b p /= GT = go (p:ps) $ b:qs
	| otherwise = go ps $ p:s
	ccw (ax, ay) (bx, by) (cx, cy)
	| d < 0 = LT
	| d > 0 = GT
	| True = EQ
	where
	d = (bx - ax) * (cy - ay) - (by - ay) * (cx - ax)
	bottomRightP = maximumBy (comparing snd) . sort
	randPoints = fmap f newPureMT
	where
	f = uncurry zip . splitAt 20 . g
	g = map (* res) . unfoldr (Just . randomDouble)

view raw graham-scan.hs hosted with ❤ by GitHub

more Π obscurantism

Here's an illustration of an approach to Π pointed out by Kevin Brown. Define f(n) as the nearest greater or equal multiple of n-1, then of n-2, etc (yielding OEIS sequence 2491). Then, inverting a result found by Duane Broline and Daniel Loeb, Π = n² / f(n).

But as you can see from the comment, the series converges very slowly!

	main = print $ e**2 / f e e 1
	where
	e = 900000 -- yields 3.1416003...
	f n 1 _ = n
	f n k l = f n' (k - 1) $ n' / k
	where
	n' = head $ dropWhile (< n) $ map (* k) [l..]

view raw kevin_brown_PI.hs hosted with ❤ by GitHub