# 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"