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"

The code.

Leave a Reply

Your email address will not be published. Required fields are marked *