Centroid of a polygon in Clojure

Question

I am studying Clojure and functional programming. I wrote this function to compute the centroid of a polygon specified as a vector of points, e.g. [[10 20] [0 11] [50 30]]).

Here you can read how to compute the centroid of a polygon.

Can you give me hints and suggestions on how to improve this code from a functional coding style POV?

(defn centroid [p]
 "Calcola il centroide del poligono p"
 (let [six*area (* 6 (polygon-area p))
 n (count p)
 first-vertex (p 0)
 polygon (conj p first-vertex)
 x-terms (atom [])
 y-terms (atom [])]
 (dotimes [i n]
 (let [point_i (polygon i)
 point_i+1 (polygon (inc i))
 x_i (point-x point_i)
 y_i (point-y point_i)
 x_i+1 (point-x point_i+1)
 y_i+1 (point-y point_i+1)]
 (swap! x-terms conj (* (+ x_i x_i+1)
 (- (* x_i y_i+1)
 (* x_i+1 y_i))))
 (swap! y-terms conj (* (+ y_i y_i+1)
 (- (* x_i y_i+1)
 (* x_i+1 y_i))))))
 (make-point (/ (reduce + 0 @x-terms) six*area)
 (/ (reduce + 0 @y-terms) six*area))))

Answer 1

You could replace the dotimes/atoms with more functional loop/recur recursion.

(defn centroid [p]
 (let [six*area (* 6 (polygon-area p))
 n (count p)
 first-vertex (p 0)
 polygon (conj p first-vertex)]
 (loop [i 0
 x-terms []
 y-terms []]
 (if (< i n)
 (let [point (polygon i)
 point' (polygon (inc i))
 x (point-x point)
 y (point-y point)
 x' (point-x point')
 y' (point-y point')
 dxy (- (* x y')
 (* x' y))]
 (recur (inc i)
 (conj x-terms (* (+ x x') dxy))
 (conj y-terms (* (+ y y') dxy))))
 (make-point (/ (reduce + 0 x-terms) six*area)
 (/ (reduce + 0 y-terms) six*area))))))

You could also do it through sequences entirely:

(defn centroid [p]
 (let [six*area (* 6 (polygon-area p))
 polygon (map (juxt point-x point-y) p) ;; turns the points in [x y] pairs
 ;; `juxt` applies each function in order and outputs the results in a list.
 polygon' (drop 1 (cycle polygon)) ;; circular permutation of polygon 
 ;; `cycle` lazily repeats and concats the input sequence.
 term (fn [[[x y] [x' y']]] ;; destructuring a pair of
 ;; [x y] pairs
 (let [dxy (- (* x y')
 (* x' y))] 
 [(* dxy (+ x x')) (* dxy (+ y y'))]))
 [x-terms yterms] (->> (map list polygon polygon') ;; makes the [point point'] pairs
 ;; when `polygon` is exhausted, the generated sequence stops, so the fact that
 ;; `polygon'` is infinite has no effect.
 (map term)
 (apply (partial map list)))] ;; turns the sequence of pairs
 ;; into a pair of sequences
 ;; `apply` turns the list argument in an argument list that it feeds to the function
 ;; - here, a sequence of pairs. `partial` takes a function and some of its arguments
 ;; and returns it with its first arguments set - here, `map list ...`. `map list`
 ;; applied to a list of pairs will call `list` with the first elements of all pairs,
 ;; then with the seconds.
 (make-point (/ (reduce + 0 x-terms) six*area)
 (/ (reduce + 0 y-terms) six*area))))

EDIT: a(n unefficient) version with map galore:

(defn centroid [p]
 (let [six*area (* 6 (polygon-area p))
 polygon (map (juxt point-x point-y) p)
 polygon' (drop 1 (cycle polygon))
 x (map first polygon)
 x' (map first polygon')
 y (map second polygon)
 y' (map second polygon')
 x**x (map * x x')
 y**y (map * y y')
 x**y (map * x y')
 y**x (map * y x')
 dxy (map - x**y y**x)
 x-terms (map * x**x dxy)
 y-terms (map * y**y dxy)]
 (make-point (/ (reduce + 0 x-terms) six*area)
 (/ (reduce + 0 y-terms) six*area))))

A full loop/recur version:

(defn centroid [[h & t :as p]]
 (let [six*area (* 6 (polygon-area p))
 a (point-x h)
 b (point-y h)]
 (loop [x a
 y b
 t t
 x-term 0
 y-term 0]
 (if-let [[h & t] t]
 (let [x' (point-x h)
 y' (point-y h)
 dxy (- (* x y')
 (* x' y))
 x-term (+ x-term (* dxy (+ x x')))
 y-term (+ y-term (* dxy (+ y y')))]
 (recur x' y' t x-term y-term))
 (let [dxy (- (* x b)
 (* a y))
 x-term (+ x-term (* dxy (+ x a)))
 y-term (+ y-term (* dxy (+ y b)))]
 (make-point (/ x-term six*area)
 (/ y-term six*area)))))))

• Domenico De Felice: Wow, the sequences version is very compact though I will need more experience to understand it completely. Looking at the loop/recur version, I missed that those terms were equal, good catch. I could also accumulate the sum instead of the terms to sum at the end. It would consume less memory and CPU cycles but I don't know if it would be as clear as it is now. PS: in the first example the variable is defined as x**y and later referred to as x^y.
• omiel: Thanks, I made the correction. I included some comments in the sequences version, I hope that helps.
• Domenico De Felice: In this line from the 3rd example: let [ ... x' (map first polygon') ...], being polygon' an infinite list, wouldn't map loop forever?
• omiel: map returns (immediately) a lazy sequence whose terms will be evaluated when they are iterated through. We never actually iterate over the whole length of polygon' (neither on those of x' or y'), but on its "product" with polygon - or one of its derivated sequences (which are finite) - through map. This "product" has the same length as the shortest of its two parameters. reduce provokes a cascading evaluation of these lazy sequences, and only the terms that are needed are evaluated.
• Domenico De Felice: Thanks for the explanation, now I am sure I understand well how they works. However I was referring to y' (map second polygon') where the only (and shortest) list is polygon', an infinite list. I tried a similar code and it looped for ever..

Answer 2

A fairly concise functional version, incorporating the area calculation:

(defn centroid [p]
 (let [pairs (partition 2 1 (conj p (first p)))
 xp (mapv (fn [[[x1 y1] [x2 y2]]] (- (* x1 y2) (* x2 y1))) pairs)
 twice-area (apply + xp)
 dot-xp (fn [xs] (apply + (map * xs xp)))
 transpose (fn [vv] (apply mapv vector vv))
 pair-sums (mapv #(map (partial reduce +) (transpose %)) pairs)]
 (mapv (comp #(/ % (* 3 twice-area)) dot-xp) (transpose pair-sums))))

Then

(centroid [[10 20] [0 11] [50 30]])

... produces

[20 61/3]