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

No comments:

Post a Comment