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

Monday, March 28, 2016

today in oblique approaches

If your standard library offers complex numbers and you're not in any hurry, you can't ask for much simpler (or more obscure) ways to compute pi than this.
import Data.Complex (magnitude, Complex((:+)))
 
pi = length $ takeWhile (< 2) $ map magnitude zs
 where
  zs = iterate (\z -> z^2 + ((-0.75) :+ 0.0000001)) 0

Wednesday, February 17, 2016

ikeda map

Continuing the theme of strange attractors, here's the well-known one embedded in the Ikeda map. The thumbnail at right shows the central 'vortex' of the attractor, and links to a larger viewport.

In these images, I've plotted the real and imaginary components along the x and y axes respectively. But the more popular way to visualize this attractor adds an extra parameter to the system and is expressed in trigonometric functions. Such adaptation of the code below yields these results.

{-# LANGUAGE NoMonomorphismRestriction #-}
 
import Data.Complex
import Graphics.UI.SDL as SDL
  
(xRes, yRes) = (1366, 768)
[a, b, k, p] = [0.85, 0.9, 0.4, 7.7]
 
main = withInit [InitVideo] $ do
  w <- setVideoMode xRes yRes 32 [NoFrame]
  s <- createRGBSurface [] 1 1 32 0 0 0 0
  fillRect s Nothing $ Pixel 0xFFFFFF
  enableEvent SDLMouseMotion False
  setCaption "Ikeda" "Ikeda"
  render w s $ ikeda $ 0 :+ 0
  SDL.flip w
  run w
 
run w = do
  e <- pollEvent
  delay 64
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> run w
 
render w s = mapM_ draw
 where
  rect (x,y) = Just $ Rect (round x) (round y) 1 1
  g (x,y)    = (x + xRes / 24, y + yRes / 2.25)
  draw z     = blitSurface s Nothing w $ rect $ g p
   where
    p = (1050 * realPart z, 1050 * imagPart z)
 
ikeda = take 1500000 . filter g . iterate f
 where
  f z      = a + b * z * exp(i * (k-p) * (1 + abs z^2))
  i        = 0 :+ 1
  g (r:+i) = -35 < r && r < 35 && -18 < i && i < 18

Saturday, February 13, 2016

hénon attractor

French astronomer Michel Hénon reported on this chaotic dynamical system in 1976.
{-# LANGUAGE NoMonomorphismRestriction #-}

import Graphics.UI.SDL as SDL
 
(xRes, yRes) = (1366, 768)
(a, b)       = (1.4, 0.3)

main = withInit [InitVideo] $ do
  w <- setVideoMode xRes yRes 32 [NoFrame]
  s <- createRGBSurface [] 1 1 32 0 0 0 0
  fillRect s Nothing $ Pixel 0xFFFFFF
  enableEvent SDLMouseMotion False
  setCaption "Hénon" "Hénon"
  render w s henon
  SDL.flip w
  run w

run w = do
  e <- pollEvent
  delay 64
  case e of
   KeyUp (Keysym SDLK_ESCAPE _ _) -> return ()
   _                              -> run w

render w s = mapM_ draw
 where
  rect (x,y) = Just $ Rect (round x) (round y) 1 1
  g (x,y)    = (500*x + xRes / 2, 900*y + yRes / 2)
  draw (x,y) = blitSurface s Nothing w $ rect $ g (x,y)

henon = take 99999 $ filter g $ iterate f (0,0)
 where
  f (x,y) = (1 - a*x^2 + y, b*x)
  g (x,y) = -1.5 < x && x < 1.5 && -0.45 < y && y < 0.45

Tuesday, January 26, 2016

langton's ant

For a round of code golf, I wrote this spare implementation of Chris Langton's remarkably simple universal computer. If you want amenities like pause, random starting pattern or even quit, check out this more complete version.
import Control.Arrow ((***))
import Control.Monad (when, join)
import Data.Set (insert, delete, member, empty)
import Graphics.UI.SDL as SDL

(xres, yres, sq) = (1366, 768, 4)

main = withInit [InitVideo] $ do
  w <- setVideoMode xres yres 32 [NoFrame]
  enableEvent SDLMouseMotion False
  setCaption "Langton's Ant" "Langton's Ant"
  run w [1..] p (0,1) empty
 where
  p = (xres `div` 2 `div` sq, yres `div` 2 `div` sq)

run w (n:ns) p v ps = do
  when (n `mod` 7 == 0) render
  run w ns (move p $ g v) (g v) $ f p ps
 where
  b          = member p ps
  (f,g)      = if b then (delete, fl) else (insert, fr)
  fr (x,y)   = if x == 0 then (-y,x) else (y,x)
  fl (x,y)   = if x == 0 then (y,x)  else (y,-x)
  move (x,y) = (x+) *** (+y)
  render     = do
    fillRect w (Just $ Rect 0 0 xres yres) $ Pixel 0
    mapM_ (draw w . join (***) (* sq)) ps
    SDL.flip w
 
draw w p = f p =<< g [SWSurface] sq sq 32 0 0 0 0
 where
  rect x y  = Just $ Rect x y sq sq
  g         = createRGBSurface
  f (x,y) s = do fillRect s (rect 0 0) $ Pixel 0xFFFFFF
                 blitSurface s (rect 0 0) w $ rect x y