

 This page uses content from Wikipedia. The original article was at K-d tree. The list of authors can be seen in the page history. As with Rosetta Code, the text of Wikipedia is available under the GNU FDL. (See links for details on variance)

A k-d tree (short for k-dimensional tree) is a space-partitioning data structure for organizing points in a k-dimensional space. k-d trees are a useful data structure for several applications, such as searches involving a multidimensional search key (e.g. range searches and nearest neighbor searches).
k-d trees are a special case of binary space partitioning trees.

k-d trees are not suitable, however, for efficiently finding the nearest neighbor in high dimensional spaces. As a general rule, if the dimensionality is k, the number of points in the data, N, should be N ≫ 2k.
Otherwise, when k-d trees are used with high-dimensional data, most of the points in the tree will be evaluated and the efficiency is no better than exhaustive search, and other methods such as approximate nearest-neighbor are used instead.

Task: Construct a k-d tree and perform a nearest neighbor search for two example data sets:

 The Wikipedia example data of [(2,3), (5,4), (9,6), (4,7), (8,1), (7,2)].
 1000 3-d points uniformly distributed in a 3-d cube.
For the Wikipedia example, find the nearest neighbor to point (9, 2)
For the random data, pick a random location and find the nearest neighbor.

In addition, instrument your code to count the number of nodes visited in the nearest neighbor search. Count a node as visited if any field of it is accessed.

Output should show the point searched for, the point found,
the distance to the point, and the number of nodes visited.

There are variant algorithms for constructing the tree.
You can use a simple median strategy or implement something more efficient.
Variants of the nearest neighbor search include nearest N neighbors, approximate nearest neighbor, and range searches.
You do not have to implement these.
The requirement for this task is specifically the nearest single neighbor.
Also there are algorithms for inserting, deleting, and balancing k-d trees.
These are also not required for the task.

 
(in-package cl-user)
 
(defvar *random-target*
  (list (float (/ (random 1000) 100))
        (float (/ (random 1000) 100))
        (float (/ (random 1000) 100))))
 
(defun distance (n target &key (semi nil))
  "distance node to target:
returns squared euclidean distance, or squared semi distance if option set"
  (if semi
      (expt (- (nth (first n) (second n))
               (nth (first n) target))
            2)
      (reduce #'+
              (mapcar (lambda (x y) (* (- x y) (- x y)))
                      (second n)
                      target))))
 
(defun target< (n target)
  "returns true if target is to its left in axis dim"
  (< (nth (first n) target)
     (nth (first n) (second n))))
 
(defun next-node (n target &key (opposite nil))
  "return the next child when nn searching, return opposing child if option
oppose set"
  (if (or (and (target< n target) (not opposite))
          (and (not (target< n target)) opposite))
      (third n)
      (fourth n)))
 
(defun make-kdtree (axis data dims)
  "a kdtree is a binary tree where nodes are:
terminal: (axis data-point),  or
branch:   (axis split-point (left-kdtree) (right-kdtree))"
  (cond ((null data) nil)   
        ((eql (length data) 1)          ; singleton?
         (list axis (first data)))      ; terminal node
        ;; else branch node:
        ;; #pts=odd  splits list into 2 even parts with sp in middle
        ;; #pts=even splits list into 2 uneven parts with shorter length first (but never nil)
        (t
         (let ((sd (sort (copy-list data)
                         #'<
                         :key (lambda (x) (nth axis x)))) ; sort the axis ordinates 
               (sp (truncate (length data) 2)) ; get mid pt
               (nxta (mod (1+ axis) dims)))
           (list axis
                 (nth sp sd)
                 (make-kdtree nxta
                              (subseq sd 0 sp)
                              dims)
                 (make-kdtree nxta
                              (subseq sd (1+ sp))
                              dims))))))
 
(defun visit-kdtree (kdt &key (node-function nil))
  "depth first visit all nodes in kdtree and optionally apply a function to
each node visited"
  (when kdt
    (when node-function (funcall node-function kdt))
    (visit-kdtree (third kdt) :node-function node-function)
    (visit-kdtree (fourth kdt) :node-function node-function)))
 
(defun count-nodes (kdt)
  "count of the terminal nodes"
  (cond ((null kdt) 0)
        ((= (length kdt) 2) 1)
        (t (+ 1
              (count-nodes (third kdt))
              (count-nodes (fourth kdt))))))
 
(defvar *hits*)
 
(defun nn-kdtree (kdt node-stack target)
  "nearest neighbour search"
  (when kdt
    ;; stage 1 - find the 'closest' terminal node using insertion logic
    (let ((best (do ((node kdt (next-node node target)))
                    ((not (next-node node target))
                     (incf *hits*)
                     node)              ; return first best est.
                  (push node node-stack)
                  (incf *hits*))))
      ;; stage 2 - unwind the path, at each node if node is closer then make it best
      (do ((node (pop node-stack) (pop node-stack)))
          ((null node) best)         ; return nearest pt                    
        ;; iteration: update best if node is closer
        (when (< (distance node target) (distance best target))
          (setf best node))
        ;; venture down opposing side if split point is inside HS
        (let ((opposing-best
                (if (< (distance node target :semi t) (distance best target))
                    (nn-kdtree (next-node node target :opposite t)
                               (list)
                               target)
                    nil)))              ; otherwise ignore this subtree
          (when (and opposing-best
                     (< (distance opposing-best target) (distance best target)))
            (setf best opposing-best)))))))
 
(defun process (data target
                &key (render nil)
                &aux (dims (length target)))
  "process one set of data & optionally display tree"
  (let* ((*hits* 0)
         (kdt (make-kdtree 0 data dims))
         (nn (nn-kdtree kdt (list) target)))
    (when render
      (visit-kdtree kdt
                    :node-function
                    (lambda (n)
                      (format t
                              "~A node: axis:~A point: ~A target:~A semi-distance-sqd:~A euclidean-distance-sqd:~A~%"
                              (if (not (next-node n target))
                                  "TERMINAL"
                                  "BRANCH")
                              (first n)
                              (second n) 
                              target
                              (distance n target :semi t)
                              (distance n target)))))        
    (format t "~%NN to ~A is ~A, distance ~A [tree has ~A nodes, ~A were visited.]~%"
            target
            (second nn)
            (sqrt (distance nn target))
            (count-nodes kdt)
            *hits*)))
 
(defun main ()
  ;; TASK 1 - nn search small set of 2D points
  (process '((2 3) (5 4) (9 6) (4 7) (8 1) (7 2))
           '(9 2)
           :render t)
 
  ;; TASK 2 - nn search 1000 coordinate points in 3D space
  (process (let ((ll (list)))
             (dotimes (i 10)
               (dotimes (j 10)
                 (dotimes (k 10)
                   (push (list i j k) ll))))
             ll)
           *random-target*))
 
 