Convex Hull with Aki-Toussaint heuristics
;; Calculates the convex hull of a set of points using the Graham scan
;; algorithm.
(ns i27.geom.convex-hull
(:use [i27.geom.misc :only (angle cross-product)]))
(defn presort-points [pts]
"Presorts the cartesian points in descending order by angle.
The point with lowest y value is last in the vector."
(let [pts (distinct pts)] ;; must be distinct or else errors ...
(letfn [(find-start [[x1 y1 :as pt1] [x2 y2 :as pt2]]
(cond (< y1 y2) pt1
(= y1 y2) (if (< x1 x2) pt1 pt2)
true pt2))]
(let [pt (reduce find-start pts)
npts (remove (fn [[x y :as cpt]]
(and (= (first pt) x) (= (last pt) y))) pts)]
(conj (apply vector
(sort (fn [pt1 pt2]
(> (angle pt pt1) (angle pt pt2))) npts))
pt)))))
(defn convex-hull
"Calculates the convex hull given a set of points in the cartesian plane.
Uses the Graham scan algorithm. Performance is O(n lg n) time."
[pts]
(when (>= (count pts) 3)
(let [start (apply vector
(reverse (subvec pts (- (count pts) 3) (count pts))))
other (reverse (subvec pts 0 (- (count pts) 3)))]
(letfn [(popstack [stk pt]
(if (not (> (cross-product
(last (pop stk)) (last stk) pt) 0))
stk (recur (pop stk) pt)))
(scan [res pt]
(conj (popstack res pt) pt))]
(reduce scan start other)))))
(defn elim-points
"Prepare data with Aki-Toussaint heuristics before passing to convex hull function.
Performance is O(n) time regardless of the convex hull algorithm used."
[pts]
(when (>= (count pts) 4)
(letfn [(find-quad
[[[x1 y1 :as pt1] [x3 y3 :as pt3]
[x2 y2 :as pt2] [x4 y4 :as pt4] :as quad] [xc yc :as cpt]]
[(if (> x1 xc) cpt pt1) (if (> y3 yc) cpt pt3)
(if (< x2 xc) cpt pt2) (if (< y4 yc) cpt pt4)])]
(let [poly (distinct (reduce find-quad (take 4 pts) pts))
;; generate the line segments that form the convex polygon
lines (conj (reduce (fn [lines pt] (conj lines [(last (last lines)) pt]))
[[(first poly) (second poly)]] (rest (rest poly)))
[(last poly) (first poly)])
;; filter out all points that fall inside the convex polygon
pts (filter
(fn [tpt]
(not (reduce (fn [b v] (and b (neg? v))) true
(map (fn [pt]
(cross-product tpt (first pt) (last pt)))
lines))))
pts)]
pts))))
Timing without heuristics :
java version "1.6.0_18" OpenJDK Runtime Environment (IcedTea6 1.8) (6b18-1.8-0ubuntu1) OpenJDK 64-Bit Server VM (build 14.0-b16, mixed mode) generating 1000000 random points ... "Elapsed time: 896.39526 msecs" sorting 1000000 points ... "Elapsed time: 52142.589419 msecs" calculating convex hull for 1000000 points ... "Elapsed time: 5370.579308 msecs" [[246336 0] [516315 0] [332671 0] [920929 0] [992271 1]] - 35 points "Elapsed time: 57547.947404 msecs"
Timing with heuristics :
java version "1.6.0_18" OpenJDK Runtime Environment (IcedTea6 1.8) (6b18-1.8-0ubuntu1) OpenJDK 64-Bit Server VM (build 14.0-b16, mixed mode) generating 1000000 random points ... "Elapsed time: 897.92648 msecs" heuristic elimination of points ... "Elapsed time: 264.116189 msecs" sorting 510685 points ... "Elapsed time: 25988.023392 msecs" calculating convex hull for 510685 points ... "Elapsed time: 3159.032778 msecs" [[246336 0] [516315 0] [332671 0] [920929 0] [992271 1]] - 35 points "Elapsed time: 32544.642604 msecs"