# Centroid of a polygon in Clojure

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]]).

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))))


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)))))))

• 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. – Domenico De Felice Jan 11 '14 at 12:39
• Thanks, I made the correction. I included some comments in the sequences version, I hope that helps. – omiel Jan 11 '14 at 13:14
• In this line from the 3rd example: let [ ... x' (map first polygon') ...], being polygon' an infinite list, wouldn't map loop forever? – Domenico De Felice Jan 11 '14 at 21:38
• 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. – omiel Jan 11 '14 at 22:53
• 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.. – Domenico De Felice Jan 12 '14 at 8:29

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]