Delaunay triangulation

delaunay.tq

let ! = Array.get

let x (x, _) = x
let y (_, y) = y
let add (x1,y1) (x2,y2) = x2+x1, y2+y1
let sub (x1,y1) (x2,y2) = x2-x1, y2-y1
let scale s (x,y) = s*x, s*y
let dot (x1,y1) (x2,y2) = x1*x2 + y1*y2
let length u = √(dot u u)
let cross (x1,y1) (x2,y2) = x1*y2 - x2*y1

(** This returns 0 if the four points are on the same circle. It returns >0 if going around the circle in the direction a -> b -> c then d is on the side of the circle pointed to by my left hand.  It returns >0 otherwise.  It's a fourth degree function in the given coordinates. *)
let inCircle t u v w = 
  let tw, uw, vw = sub t w, sub u w, sub v w in
  let twd = dot tw tw in
  let uwd = dot uw uw in
  let vwd = dot vw vw in
  x tw*(y uw*vwd - uwd*y vw) - y tw*(x uw*vwd - uwd*x vw) +
    twd*(x uw*y vw - y uw*x vw)

(** Radius of a circle given three points. *)
let radius a b c  = 
  let a = sub a c in
  let b = sub b c in
  length a * length b * length(sub a b) / (2*(abs(cross a b)))

(** Returns >0 if c on the left side of a->b, <0 if on the right, and 0 if on. *)
let leftside a b c = cross (sub a c) (sub b c)

let len2 u v =
  let d = sub u v in
  dot d d

let rec closest p (oi, oj as o) (bi, bj as b) dist = 
  if oj=Array.length p then b else
    let d = len2 (! p oi) (! p oj) in
    let next = if oi+1 < oj then oi+1, oj else 0, oj+1 in
    if bj = 0 || d<dist then closest p next o d else closest p next b dist

(* This returns an array of indices of this face in counterclockwise order. Think of i pointing up to j. All the other points of the face are to the left of this line.  We sort them by taking the dot product of the i->j vector with the normalized j->k vector, and sort it in decreasing order. *)
let orderFace p i j face nv =
  let v = Array.init (nv-2) [_ -> 0, 0] in
  let () =
    face
    @ Unsafe.Stack.iteri [k, h ->
      let a = sub (! p j) (! p i) in
      let b = sub (! p h) (! p j) in
      let dp = dot a b / length b in
      Array.Unsafe.set v k (dp, h)] in
  let () = Array.Unsafe.sortInPlaceWith [a, b -> compare(b, a)] v in
  Array.init nv [k -> if k=0 then i else if k=1 then j else y(! v (k-2))]

let rec findface p i j =
  let ac = Unsafe.Stack.empty() in
  let () =
    for 0 1 (Array.length p-1) () [(), k ->
      if k=i || k=j || leftside (! p i) (! p j) (! p k) <= 0 then () else
        Unsafe.Stack.peek ac
        @ [ None -> Unsafe.Stack.push k ac
          | Some l ->
              compare(inCircle (! p i) (! p j) (! p l) (! p k), 0)
              @ [ Less -> ()
                | Equal -> Unsafe.Stack.push k ac
                | Greater ->
                    let () = Unsafe.Stack.clear ac in
                    Unsafe.Stack.push k ac ] ] ] in
  ac

let rec processFace nv v k s t =
  if k = nv then () else
    let l = mod(k+1, nv) in
    let () = Unsafe.Set.add s (! v k, ! v l) in
    let () = Unsafe.Stack.push (! v l, ! v k) t in
    processFace nv v (k+1) s t

(** Takes as input an array of p of n points, and computes the Delaunay triangulation of them. *)
let delaunay p =
  let (bi, bj as b) = closest p (0,1) (0,0) 0 in
  let completed = Unsafe.Set.empty() in
  let todo = Unsafe.Stack.ofArray{b; (bj, bi)} in
  let ac = {} in
  let rec loop ac =
    Unsafe.Stack.pop todo
    @ [ None -> ac
      | Some(oi, oj as o) ->
          let ac =
            if Unsafe.Set.contains completed o then ac else
              let face = findface p oi oj in
              if Unsafe.Stack.count face = 0 then
                (* Convex hull edge *)
                let () = Unsafe.Set.add completed o in
                ac
              else
                let nv = 2 + Unsafe.Stack.count face in
                let v = orderFace p oi oj face nv in
                let () = processFace nv v 0 completed todo in
                Array.Unsafe.append ac v in
          loop ac ] in
  loop {}

let p = Array.init 100 [_ -> Random.next 800, Random.next 600]

let () =
  delaunay p
  @ Array.map [ns ->
    let ps = Array.map [n -> Array.get p n] ns in
    let tp = Array.fold [u, v -> add u v] (0, 0) ps @ scale (1 / # ps) in
    Array.map [p -> add (scale 0.9 p) (scale 0.1 tp)] ps
    @ polygon (Style(Fill Red, NoStroke))]
  @ Html.ofSvg (800, 600)
  @ yield