{-# 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

## Friday, February 24, 2017

### bézier curves

De Casteljau's algorithm is a fast, numerically stable way to rasterize Bézier curves (I was disappointed to see Wikipedia already gives succinct Haskell for it!). Clicking appends nodes to the curve's defining polygon; you can also drag any node to alter the curve.

## Wednesday, February 22, 2017

### 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 this algorithm, noted by A.M Andrew, sorts the set lexicographically and finds the upper and lower hull chains separately. This 'monotone chain' technique is often preferred, since it's easier to do robustly.

I've also written a couple JavaScript examples, the Graham scan and a slower (but constant-space and adaptable to higher dimensions) method called the Jarvis march.

{-# 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)

## Wednesday, February 15, 2017

### more Π obscurantism

Here's an illustration of an approach to pi 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, pi = n^{2}/ 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..]

Subscribe to:
Posts (Atom)