Simple Ramer-Douglas-Peucker algorithm in clojure

rdp.clj

(ns rdp.core)

; Ramer-Douglas-Peucker algorithm based on http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm

(defn- dist [point start end]
  (let [num (- (* (- (first end) (first start))
                  (- (last start) (last point)))

               (* (- (first start) (first point))
                  (- (last end) (last start))))

        den (+ (Math/pow (- (first end) (first start)) 2)
               (Math/pow (- (last end) (last start)) 2))]
    (/ (Math/abs num) (Math/pow den 0.5))))

(defn exec
  "
    Usage:
      (exec [[1 2] [2 3] [3 4 ] [4 6] [4 2]] 0.5) => [[1 2] [4 6] [4 2]]
  "
  [points eps]
  (if (< (count points) 3)
    points
    (do
      (let [d-max (atom 0.0)
            index (atom -1)
            first-item (first points)
            last-item (last points)]
        (doall
          (map-indexed
            (fn [i item]
              (let [d (dist item first-item last-item)]
                (when (> d @d-max)
                  (reset! d-max d)
                  (reset! index (inc i)))))                 ; add 1 because we start at 0 but we're dropping 1 item
            (rest points)))

        (if (> @d-max eps)
          (concat
            (drop-last 1 (exec (take @index points) eps))
            (exec (drop @index points) eps))
          [first-item last-item])))))

TODO: remove atoms and use loop

Not sure if you've thought about this code snippet for awhile, but I came across this while doing some shapefile processing.

Anyway, I figured out a more concise way to do it without atoms:

(defn- point-dist-to-line
  "calculates perpendicular distance of point to line defined by two points"
  [point start end]
  (let [[[x0 y0] [x1 y1] [x2 y2]] [point start end]
         numerator (- (* (- x2 x1) (- y1 y0)) (* (- x1 x0) (- y2 y1)))
         denominator (Math/sqrt (+ (* (- x2 x1) (- x2 x1)) (* (- y2 y1) (- y2 y1))))]
    (Math/abs (/ numerator denominator))))

(defn- rdp-exec
  [points epsilon]
  (if (< (count points 3))
    points
    (let [[dmax index] (loop [dmax 0 index -1 i 1]
                                    (let [d (point-dist-to-line (get points i) (first points) (last points))]
                                      (if (< i (dec (count points)))
                                        (if (> d dmax)
                                          (recur d i (inc i))
                                          (recur dmax index (inc i)))
                                        (if (> d dmax)
                                          [d i]
                                          [dmax index]))))]
      (if (> dmax epsilon)
        (vec (mapcat #(rdp-exec (vec %) epsilon) (split-at index points)))
        ((juxt first last) points)))))

For some reason, mine keeps [3 4] in the output pair given your input though.

Anyway, happy Clojuring!