Sunday, June 28, 2015

triangle geometry

Illustrates the centroid, incircle, and circumcircle of a random triangle; the viewport centers on the circumcenter. All functions are defined in terms of vertex triples. The triangle interior is filled using barycentric coordinates.
import Control.Arrow ((***))
import Control.Monad (ap, liftM2, void, when)
import Data.List (unfoldr)
import Graphics.UI.SDL as SDL
import Graphics.UI.SDL.Primitives (circle, filledCircle, pixel)
import System.Random.Mersenne.Pure64 (newPureMT, randomInt)
 
(xres, yres) = (1600, 900)
 
main = withInit [InitVideo] $ do
  w   <- setVideoMode xres yres 32 [NoFrame]
  rng <- newPureMT
  enableEvent SDLMouseMotion False
  setCaption "Triangle" "Triangle"
  render w rng $ liftM2 (,) [0.. xres] [0.. yres]
  SDL.flip w
  pause w
 
pause w = do
  delay 128
  e <- pollEvent
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> pause w

render w rng ps = do
  let pxl (x,y) = void . pixel w x y $ Pixel 0x0000FFFF
  let (x,y)     = centroid tri
  let (p,q)     = circumcenter tri
  let (a,b,r)   = incircle tri
  let rad       = dist (p,q) $ head tri
  let f         = round . (+ (fromIntegral xres / 2 - p))
  let g         = round . (+ (fromIntegral yres / 2 - q))
  let h         = (f . fromIntegral) *** (g . fromIntegral)
  mapM_ (ap (when . within tri) (pxl . h)) ps
  filledCircle w (f x) (g y) 2 $ Pixel 0xFFFF00FF
  filledCircle w (f a) (g b) 2 $ Pixel 0xFF00FF
  circle w (f p) (g q) rad $ Pixel 0xFF0000FF
  circle w (f a) (g b) (round r) $ Pixel 0xFF00FF
 where
  tri = head $ filter p $ triples rng
  p t = let n = area t in n > 15000 && n < 40000

------------------------------------------------------------

within [(ax,ay), (bx,by), (cx,cy)] (x', y') = f [a,b,g]
 where
  (x,y) = (fromIntegral x', fromIntegral y')
  a     = ((by-cy) * (x-cx)  + (cx-bx) * (y-cy)) / d
  b     = ((cy-ay) * (x-cx)  + (ax-cx) * (y-cy)) / d
  d     =  (by-cy) * (ax-cx) + (cx-bx) * (ay-cy)
  g     = 1 - a - b;
  f     = all (> 0)

circumcenter [(ax,ay), (bx,by), (cx,cy)] = (x/d, y/d)
 where
  x = ((ax-cx) * (ax+cx) + (ay-cy) * (ay+cy)) / 2 * (by-cy)
    - ((bx-cx) * (bx+cx) + (by-cy) * (by+cy)) / 2 * (ay-cy)
  y = ((bx-cx) * (bx+cx) + (by-cy) * (by+cy)) / 2 * (ax-cx)
    - ((ax-cx) * (ax+cx) + (ay-cy) * (ay+cy)) / 2 * (bx-cx)
  d =  (ax-cx) * (by-cy) - (bx-cx) * (ay-cy)

incircle [(x0,y0), (x1,y1), (x2,y2)] = (x/d, y/d, r)
 where
  x = a*x0 + b*x1 + c*x2
  y = a*y0 + b*y1 + c*y2
  a = fromIntegral $ dist (x1,y1) (x2,y2)
  b = fromIntegral $ dist (x2,y2) (x0,y0)
  c = fromIntegral $ dist (x0,y0) (x1,y1)
  d = a + b + c
  r = let s = d/2 in sqrt ((s-a)*(s-b)*(s-c) / s)

centroid [(ax,ay), (bx,by), (cx,cy)] = (x,y)
 where
  (x,y) = ((ax+bx+cx) / 3, (ay+by+cy) / 3)

area [(ax,ay), (bx,by), (cx,cy)] = abs $ (a+b+c) / 2
 where
  a = ax * (by - cy)
  b = bx * (cy - ay)
  c = cx * (ay - by)

dist (p,p') (q,q') = round . sqrt $ (p-q)^2 + (p'-q')^2

triples rng = let f xs = take 3 xs : f (drop 3 xs) in f ps
 where
  ps          = go $ map ((+ 10) . (`mod` 580)) rands
  rands       = unfoldr (Just . randomInt) rng
  go (x:y:ns) = (fromIntegral x, fromIntegral y) : go ns