Source code in the racket programming language

#lang racket

(require plot)

;; Performs clustering of points in a grid
;; using the nearest neigbour approach and shows
;; clusters in different colors
(define (plot-Voronoi-diagram point-list)
  (define pts
    (for*/list ([x (in-range 0 1 0.005)]
                [y (in-range 0 1 0.005)])
      (vector x y)))
  
  (define clusters (clusterize pts point-list))
  
  (plot
   (append
    (for/list ([r (in-list clusters)] [i (in-naturals)])
      (points (rest r) #:color i #:sym 'fullcircle1))
    (list (points point-list #:sym 'fullcircle5 #:fill-color 'white)))))

;; Divides the set of points into clusters
;; using given centroids
(define (clusterize data centroids)
  (for*/fold ([res (map list centroids)]) ([x (in-list data)])
    (define c (argmin (curryr (metric) x) centroids))
    (dict-set res c (cons x (dict-ref res c)))))


(define (euclidean-distance a b)
  (for/sum ([x (in-vector a)] [y (in-vector b)]) 
    (sqr (- x y))))

(define (manhattan-distance a b)
  (for/sum ([x (in-vector a)] [y (in-vector b)]) 
    (abs (- x y))))

(define metric (make-parameter euclidean-distance))


;; Plots the Voronoi diagram as a contour plot of
;; the classification function built for a set of points
(define (plot-Voronoi-diagram2 point-list)
  (define n (length point-list))
  (define F (classification-function point-list))
  (plot 
   (list
    (contour-intervals (compose F vector) 0 1 0 1
                       #:samples 300
                       #:levels n
                       #:colors (range n)
                       #:contour-styles '(solid)
                       #:alphas '(1))
    (points point-list #:sym 'fullcircle3))))

;; For a set of centroids returns a function 
;; which finds the index of the centroid nearest 
;; to a given point
(define (classification-function centroids)
  (define tbl 
    (for/hash ([p (in-list centroids)] [i (in-naturals)])
      (values p i)))
  (λ (x)
    (hash-ref tbl (argmin (curry (metric) x) centroids))))


(define pts
  (for/list ([i 50]) (vector (random) (random))))

(display (plot-Voronoi-diagram pts))

(display (plot-Voronoi-diagram2 pts))

(parameterize ([metric manhattan-distance])
  (display (plot-Voronoi-diagram2 pts)))

;; Using the classification function it is possible to plot Voronoi diagram in 3D.
(define pts3d (for/list ([i 7]) (vector (random) (random) (random))))
(plot3d (list
         (isosurfaces3d (compose (classification-function pts3d) vector) 
                        0 1 0 1 0 1
                        #:line-styles '(transparent)
                        #:samples 100
                        #:colors (range 7)
                        #:alphas '(1))
         (points3d pts3d #:sym 'fullcircle3)))