diff --git a/CgSceneGraph/CgAABB.cpp b/CgSceneGraph/CgAABB.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4d8354ef2929a9e0e32881753930680d6ca8848f
--- /dev/null
+++ b/CgSceneGraph/CgAABB.cpp
@@ -0,0 +1,27 @@
+#include "CgSceneGraph/CgAABB.h"
+
+CgAABB::CgAABB(glm::vec3 pos, glm::vec3 ext){
+    position = pos;
+    extent = ext;
+}
+
+CgAABB::CgAABB() {
+    position = glm::vec3(0.0);
+    extent = glm::vec3(0.0);
+}
+
+std::pair<CgAABB, CgAABB> CgAABB::split(int axis, double value) {
+    glm::vec3 new_position1 = position;
+    glm::vec3 new_position2 = position;
+    double delta = value - position[axis];
+    new_position1[axis] -= (extent[axis] - delta) * 0.5;
+    new_position2[axis] += (extent[axis] + delta) * 0.5; 
+
+    glm::vec3 new_extent1 = extent;
+    glm::vec3 new_extent2 = extent;
+    new_extent1[axis] += delta;
+    new_extent2[axis] -= delta;
+    new_extent1[axis] *= .5;
+    new_extent2[axis] *= .5;
+    return {CgAABB(new_position1, new_extent1), CgAABB(new_position2, new_extent2)};
+}
\ No newline at end of file
diff --git a/CgSceneGraph/CgAABB.h b/CgSceneGraph/CgAABB.h
index 48f4c32d90728b4f8d6645e6f8f6068414066f45..d2026bfce251fc08c8f1abbe3d1b53d76ad9f42e 100644
--- a/CgSceneGraph/CgAABB.h
+++ b/CgSceneGraph/CgAABB.h
@@ -1,11 +1,18 @@
 #ifndef CGAABB_H
 #define CGAABB_H
 
+#include <utility>
 #include <glm/glm.hpp>
 
-struct CgAABB {
+class CgAABB {
+public:
     glm::vec3 position;
     glm::vec3 extent;
+    CgAABB();
+    CgAABB(glm::vec3 pos, glm::vec3 ext);
+    // split a aabb on an axis into two aabb whose union occupies the same space
+    std::pair<CgAABB, CgAABB> split(int axis, double value);
 };
 
+
 #endif // CGAABB_H
\ No newline at end of file
diff --git a/CgSceneGraph/CgPointCloud.cpp b/CgSceneGraph/CgPointCloud.cpp
index c419e3f0b6816e4bcf0071d5ab44518e914c319b..6e5f60dc5ac06567fe2db3d940eaf21e34a2b89d 100644
--- a/CgSceneGraph/CgPointCloud.cpp
+++ b/CgSceneGraph/CgPointCloud.cpp
@@ -1,6 +1,5 @@
 #include "CgPointCloud.h"
 
-
 CgPointCloud::CgPointCloud():
 CgPointCloud::CgPointCloud(51)
 {
@@ -72,14 +71,17 @@ void CgPointCloud::init( std::string filename, bool cheat_normals)
     // system, use up vector = y-Axis of your local 
     // coordinate system instead of getPerpendicularVector(...)
     calculateSplatOrientations();
+    calculateAABB();
 
     //add a standard color for each point if lighting turned off
-    for(unsigned int i=0;i<m_vertices.size();i++) {
+    for(unsigned int i = 0;i < m_vertices.size(); i++) {
       m_vertex_colors.push_back(glm::vec3(0.0,1.0,0.0));
     }
 
-    unsigned int k=50;
-    std::vector<int> neighbors = getNearestNeighbors(0,k);
+    unsigned int k = 50;
+    // "fast" if k is relatively small, if k is on the order 
+    // of m_vertices.size(), this is slower than brute force.
+    std::vector<unsigned int> neighbors = getNearestNeighborsFast(0, k);
 
     for(unsigned int i = 0;i < k; i++) {
       m_vertex_colors[neighbors[i]] = glm::vec3(0.0,0.0,1.0);
@@ -133,39 +135,110 @@ void CgPointCloud::traverseBFO(std::function<void(glm::vec3*, const unsigned int
     }
 }
 
-std::vector<int> CgPointCloud::getNearestNeighbors(int current_point,unsigned int k) {
+std::vector<unsigned int> CgPointCloud::getNearestNeighborsSlow(unsigned int current_point,unsigned int k) {
 
   glm::vec3 q = m_vertices[current_point];
 
   std::vector<std::pair<double,int>> distances;
 
-  // very inefficient, just to show that it works for rendering colored neighborhood
-  // use min heap for real purposes
+  for(unsigned int i = 0; i < m_vertices.size(); i++) {
+    double dist = glm::distance(m_vertices[i],q);
+    distances.push_back(std::make_pair(dist,i));
+  }
 
+    std::sort(distances.begin(), distances.end());
 
-  for(unsigned int i=0;i<m_vertices.size();i++)
-    {
-      double dist=glm::distance(m_vertices[i],q);
+    std::vector<unsigned int> erg;
 
-      distances.push_back(std::make_pair(dist,i));
-    }
+  for(unsigned int i = 0; i < k; i++) {
+       erg.push_back(distances[i].second);
+  }
 
-    std::sort(distances.begin(),distances.end());
+  return erg;
+}
 
-    std::vector<int> erg;
+std::vector<unsigned int> CgPointCloud::getNearestNeighborsFast(unsigned int current_point, unsigned int k) {
+
+  std::vector<unsigned int> erg;
+  if(k == 0 || m_vertices.size() == 0) return erg;
+  // pair of distance + index into backing vec3 array 
+  using P = std::pair<double, unsigned int>;
+
+  glm::vec3 query = m_vertices[current_point];
+  glm::vec3* first = &m_vertices.front();
+  glm::vec3* last = &(*(m_vertices.end() - 1));
+  auto cmp = [&](P left, P right) {
+    // std::set sees elements as equivalent if neither 
+    // compares less than the other
+    // -> needs to take index into consideration since 
+    // different points can have same distance.
+    return left.first != right.first
+      ? left.first < right.first
+      : left.second < right.second;
+  };
+
+  // using an ordered set for simplicitly, this may be slow
+  // because we frequently erase elements, causing the rb-tree
+  // to be rebuilt.
+  std::set<P, decltype(cmp)> set(cmp);
+
+  // stack of traversed nodes (start index, end index, aabb, depth)
+  std::vector<std::tuple<glm::vec3*, glm::vec3*, CgAABB, unsigned int>> nodes;
+  nodes.push_back({first, last, m_aabb, 0});
+  
+  // traverse the tree as long as there's something to do.
+  while(!nodes.empty()) {
+    // unpack current node & pop it
+    glm::vec3* begin;
+    glm::vec3* end;
+    CgAABB aabb;
+    unsigned int depth;
+    std::tie(begin, end, aabb, depth) = nodes.back();
+    auto half_size = (end - begin) / 2;
+    nodes.pop_back();
+
+    glm::vec3 curr_point = *(begin + half_size);
+    // insert current point into set
+    set.insert({glm::distance(query, curr_point), begin + half_size - first});
+    // make sure our set only contains k elements at most
+    while(set.size() > k) set.erase(std::prev(set.end()));
+
+    // the aabb of the sub-nodes
+    unsigned int axis = depth % 3;
+    CgAABB lower_aabb;
+    CgAABB higher_aabb;
+    std::tie(lower_aabb, higher_aabb) = aabb.split(axis, curr_point[axis]);
+    // only push lower child if its aabb intersects with the 
+    // neighbourhood we found 
+    double neighbourhood_radius = (*std::prev(set.end())).first;
+    if(
+      begin != begin + half_size
+      && lower_aabb.position[axis] + lower_aabb.extent[axis] > query[axis] - neighbourhood_radius
+    ) {
+      nodes.push_back({begin, begin + half_size, lower_aabb, depth + 1});
+    }
+    // only push higher if it intersects with
+    // the neighbourhood we found.
+    if(
+      begin + half_size != end
+      && higher_aabb.position[axis] - higher_aabb.extent[axis] < query[axis] + neighbourhood_radius
+    ) {
+      nodes.push_back({begin + half_size + 1, end, higher_aabb, depth + 1});
+    }
+  }
 
-   for(unsigned int i=0;i<k;i++)
-    {
-       erg.push_back(distances[i].second);
-     }
+  // copy all nearest neighbour indices we found into
+  // a vector to return
+  std::for_each(set.cbegin(), set.cend(), [&](P entry) {
+    erg.push_back(entry.second);
+  });
 
-    return erg;
-  }
+  return erg;
+}
 
 
 // calculates an arbitrary verctor perpendicular to the given one
-glm::vec3 CgPointCloud::getPerpendicularVector(glm::vec3 arg)
-{
+glm::vec3 CgPointCloud::getPerpendicularVector(glm::vec3 arg) {
   if((arg[0]==0.0)&&(arg[1]==0.0))
     {
     if(arg[2]==0.0)
@@ -176,8 +249,7 @@ glm::vec3 CgPointCloud::getPerpendicularVector(glm::vec3 arg)
   return glm::normalize(glm::vec3(-arg[1],arg[0],0.0));
 }
 
-const glm::vec3 CgPointCloud::getCenter() const
-{
+const glm::vec3 CgPointCloud::getCenter() const {
   glm::vec3 center(0.);
   for(unsigned int i=0;i<m_vertices.size();i++) {
       center+=m_vertices[i];
@@ -186,15 +258,13 @@ const glm::vec3 CgPointCloud::getCenter() const
   return center;
 }
 
-// returns the point clouds AABB
-const CgAABB CgPointCloud::getAABB() const {
+void CgPointCloud::calculateAABB() {
   double max_x = -std::numeric_limits<double>::infinity();
   double max_y = -std::numeric_limits<double>::infinity();
   double max_z = -std::numeric_limits<double>::infinity();
   double min_x = std::numeric_limits<double>::infinity();
   double min_y = std::numeric_limits<double>::infinity();
   double min_z = std::numeric_limits<double>::infinity();
-  CgAABB aabb;
 
   for(unsigned int i=0; i < m_vertices.size(); i++) {
     glm::vec3 curr = m_vertices[i];
@@ -205,9 +275,14 @@ const CgAABB CgPointCloud::getAABB() const {
     max_z = max_z < curr[2] ? curr[2] : max_z;
     min_z = min_z > curr[2] ? curr[2] : min_z;
   }
-  aabb.position = glm::vec3((max_x + min_x) * 0.5, (max_y + min_y) * 0.5, (max_z + min_z) * 0.5);
-  aabb.extent = glm::vec3((max_x - min_x) * 0.5, (max_y - min_y) * 0.5, (max_z - min_z) * 0.5);
-  return aabb;
+
+  m_aabb.position =  glm::vec3((max_x + min_x) * 0.5, (max_y + min_y) * 0.5, (max_z + min_z) * 0.5);
+  m_aabb.extent = glm::vec3((max_x - min_x) * 0.5, (max_y - min_y) * 0.5, (max_z - min_z) * 0.5);
+}
+
+// returns the point clouds AABB
+const CgAABB CgPointCloud::getAABB() const {
+  return m_aabb;
 }
 
 const std::vector<glm::vec2>& CgPointCloud::getSplatScalings() const
diff --git a/CgSceneGraph/CgPointCloud.h b/CgSceneGraph/CgPointCloud.h
index ddbfba204b3082b2074a2fcbfeff26287e7730f5..251ec4fcfabf04b9b2ac4a67f191d738c50c34c9 100644
--- a/CgSceneGraph/CgPointCloud.h
+++ b/CgSceneGraph/CgPointCloud.h
@@ -4,21 +4,22 @@
 #include <vector>
 #include <queue>
 #include <limits>
-#include <glm/glm.hpp>
 #include <string>
 #include <algorithm>
 #include <iostream>
 #include <numeric>
 #include <utility>
 #include <functional>
+#include <tuple>
+#include <set>
 #include "CgBase/CgBasePointCloud.h"
 #include "CgSceneGraph/CgAABB.h"
 #include "CgBase/CgEnums.h"
 #include "CgUtils/ObjLoader.h"
+#include <glm/glm.hpp>
 #include <glm/gtc/matrix_transform.hpp>
 
-class CgPointCloud : public CgBasePointCloud
-{
+class CgPointCloud : public CgBasePointCloud {
 public:
   CgPointCloud();
   CgPointCloud(int id);
@@ -29,7 +30,6 @@ public:
   Cg::ObjectType getType() const;
   unsigned int getID() const;
 
-
   //inherited from CgBasePointCloud
 
   // vertex positions in local coordinates
@@ -68,12 +68,14 @@ private:
 
     //for demonstration: find local coordinate system (normal plus arbitrary tangent spanning directions)
     void calculateSplatOrientations();
+    void calculateAABB();
 
     // for demonstration: for a given normal direction find an arbitrary vector to span the tangent plane
     glm::vec3 getPerpendicularVector(glm::vec3 arg);
 
     // for demonstration purposes, very inefficient
-    std::vector<int> getNearestNeighbors(int current_point, unsigned int k);
+    std::vector<unsigned int> getNearestNeighborsSlow(unsigned int current_point, unsigned int k);
+    std::vector<unsigned int> getNearestNeighborsFast(unsigned int current_point, unsigned int k);
 
     // rearrange the vec3 between begin and end so they form a kd-tree
     void buildKdTree(glm::vec3* begin, glm::vec3* end, int depth);
@@ -81,6 +83,7 @@ private:
     std::vector<glm::vec3> m_vertices;
     std::vector<glm::vec3> m_vertex_normals;
     std::vector<glm::vec3> m_vertex_colors;
+    CgAABB m_aabb;
 
     // indices of vertices for which a splat will be rendered
     std::vector<unsigned int> m_splat_indices;
diff --git a/CgSceneGraph/CgSceneControl.cpp b/CgSceneGraph/CgSceneControl.cpp
index 997e33991b5d00b64a0864a2ea22003203400c0d..a27161fba06be2f4b19c414d4c9c06297d0f41cb 100644
--- a/CgSceneGraph/CgSceneControl.cpp
+++ b/CgSceneGraph/CgSceneControl.cpp
@@ -63,6 +63,7 @@ CgSceneControl::~CgSceneControl()
         delete m_disc;
     if(m_select_ray!=NULL)
         delete m_select_ray;
+    m_kd_tree_mesh.clear();
 }
 
 
@@ -204,20 +205,9 @@ void CgSceneControl::createKdTreeViz(int max_depth) {
       m_kd_tree_mesh.push_back(new_quad);
 
       // push children of current aabb
-      glm::vec3 new_position1 = current_aabb.position;
-      glm::vec3 new_position2 = current_aabb.position;
-      double delta = current_split - current_aabb.position[current_axis];
-      new_position1[current_axis] -= (current_aabb.extent[current_axis] - delta) * 0.5;
-      new_position2[current_axis] += (current_aabb.extent[current_axis] + delta) * 0.5; 
-
-      glm::vec3 new_extent1 = current_aabb.extent;
-      glm::vec3 new_extent2 = current_aabb.extent;
-      new_extent1[current_axis] += delta;
-      new_extent2[current_axis] -= delta;
-      new_extent1[current_axis] *= .5;
-      new_extent2[current_axis] *= .5;
-      aabb_queue.push({new_position1, new_extent1});
-      aabb_queue.push({new_position2, new_extent2});
+      std::pair<CgAABB, CgAABB> children = current_aabb.split(current_axis, current_split);
+      aabb_queue.push(children.first);
+      aabb_queue.push(children.second);
 
       index += 1;
       if(index >= ((0x01 << (depth + 1)) - 1)) depth += 1;
diff --git a/CgUtils/Timer.h b/CgUtils/Timer.h
index fa35f43c85ecdbf40cd021cc91c28213e1843fb0..9f1a4cced61f02fda6c3156047f6cbedc2bdc0b3 100644
--- a/CgUtils/Timer.h
+++ b/CgUtils/Timer.h
@@ -14,7 +14,7 @@ void run_timed(unsigned int amount, std::function<void ()> fn) {
     for(int i = 0; i < amount; i++) fn();
     auto t2 = std::chrono::high_resolution_clock::now();
     std::chrono::duration<double, std::milli> ms_double = t2 - t1;
-    std::cout << ms_double.count() << "ms" << std::endl;
+    std::cout << "average: " << ms_double.count() / double(amount) << " ms" << std::endl;
 }
 
 #endif // TIMER_H
\ No newline at end of file
diff --git a/ExerciseVC.pro b/ExerciseVC.pro
index f7cb3432e8ffd2006512908262d5de3db9ae01dd..7bf2a922d822eb1fa44e0feef561d756fa12e9ca 100644
--- a/ExerciseVC.pro
+++ b/ExerciseVC.pro
@@ -26,6 +26,7 @@ SOURCES += main.cpp \
     CgQtViewer/CgTrackball.cpp \
     CgEvents/CgWindowResizeEvent.cpp \
     CgSceneGraph/CgTriangleMesh.cpp \
+    CgSceneGraph/CgAABB.cpp \
     CgUtils/ObjLoader.cpp \
     CgEvents/CgTrackballEvent.cpp