Friday, June 3, 2016

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

No comments:

Post a Comment