diff --git a/CgBase/CgEnums.h b/CgBase/CgEnums.h
index 6ef90ebb772c1e8e15e3e30428ef9dbd9b0b3287..f96432cebefeb81e0430948e42cfdf4fe75d406e 100644
--- a/CgBase/CgEnums.h
+++ b/CgBase/CgEnums.h
@@ -39,9 +39,8 @@ namespace Cg {
     } EventType;
 
 	typedef enum ClusteringType {
-		CgIncrementalNumber = 0x00,
-		CgIncrementalDiameter = 0x01,
-		CgHierarchical = 0x02,
+		CgIncremental = 0x00,
+		CgHierarchical = 0x01
 	} ClusteringType;
 
 
diff --git a/CgEvents/CgSimplifyEvent.cpp b/CgEvents/CgSimplifyEvent.cpp
index 461fbdada93d64b355f0e554cb942179e7f338fe..b1339b1912418c575da317597d104e61918014c6 100644
--- a/CgEvents/CgSimplifyEvent.cpp
+++ b/CgEvents/CgSimplifyEvent.cpp
@@ -7,7 +7,7 @@ CgBaseEvent* CgSimplifyEvent::clone()
         m_do_reset, 
         m_max_number_of_points, 
         m_max_diameter_of_cluster, 
-        m_max_variance_of_cluster, 
+        m_max_variance_of_cluster,
         m_clustering_type
     );
 }
@@ -17,6 +17,11 @@ Cg::EventType CgSimplifyEvent::getType()
     return m_type;
 }
 
+Cg::ClusteringType CgSimplifyEvent::getClusteringType() const { return m_clustering_type;}
+bool CgSimplifyEvent::doReset() const { return m_do_reset;}
+float CgSimplifyEvent::getMaxPoints() const {return m_max_number_of_points;}
+float CgSimplifyEvent::getMaxDiameter() const {return m_max_diameter_of_cluster;}
+float CgSimplifyEvent::getMaxVariance() const {return m_max_variance_of_cluster;}
 
 std::ostream& operator<<(std::ostream& os,const CgSimplifyEvent& e)
 {
@@ -38,9 +43,9 @@ CgSimplifyEvent::CgSimplifyEvent()
 CgSimplifyEvent::CgSimplifyEvent(
     Cg::EventType type,
     bool do_reset, 
-    unsigned int max_number_of_points, 
+    float max_number_of_points, 
     float max_diameter_of_cluster, 
-    float max_variance_of_cluster, 
+    float max_variance_of_cluster,
     Cg::ClusteringType clustering_type
 ) {
     m_type = type;
@@ -49,5 +54,4 @@ CgSimplifyEvent::CgSimplifyEvent(
     m_max_diameter_of_cluster = max_diameter_of_cluster;
     m_max_variance_of_cluster = max_variance_of_cluster;
     m_clustering_type = clustering_type;
-}
-
+}
\ No newline at end of file
diff --git a/CgEvents/CgSimplifyEvent.h b/CgEvents/CgSimplifyEvent.h
index f8d0d4f6e2e921e75151ac941b61e0b92908a16b..1fd9d4131454fb0f3ca54bd09c13050477659b48 100644
--- a/CgEvents/CgSimplifyEvent.h
+++ b/CgEvents/CgSimplifyEvent.h
@@ -13,13 +13,19 @@ public:
     CgSimplifyEvent(
         Cg::EventType type,
         bool do_reset, 
-        unsigned int max_number_of_points, 
+        float max_number_of_points, 
         float max_diameter_of_cluster, 
-        float max_variance_of_cluster, 
+        float max_variance_of_cluster,
         Cg::ClusteringType clustering_type
     );
     ~CgSimplifyEvent(){};
 
+    Cg::ClusteringType getClusteringType() const;
+    bool doReset() const;
+    float getMaxPoints() const;
+    float getMaxDiameter() const;
+    float getMaxVariance() const;
+
     //inherited
     Cg::EventType  getType();
     CgBaseEvent* clone();
@@ -29,9 +35,10 @@ public:
 private:
      Cg::EventType m_type;
      bool m_do_reset;
-     unsigned int m_max_number_of_points;
+     float m_max_number_of_points;
      float m_max_diameter_of_cluster;
      float m_max_variance_of_cluster;
+     bool m_use_number_of_points;
      Cg::ClusteringType m_clustering_type;
 };
 #endif // CG_SIMPLIFYEVENT
\ No newline at end of file
diff --git a/CgMath/CgAABB.cpp b/CgMath/CgAABB.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ac5bd1d5f2af878ace6a73cd5784bcb9439b3f7f
--- /dev/null
+++ b/CgMath/CgAABB.cpp
@@ -0,0 +1,29 @@
+#include "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) {
+    double delta = value - position[axis];
+
+    glm::vec3 new_position1 = position;
+    glm::vec3 new_extent1 = extent;
+    new_position1[axis] -= (extent[axis] - delta) * 0.5;
+    new_extent1[axis] += delta;
+    new_extent1[axis] *= .5;
+
+    glm::vec3 new_position2 = position;
+    glm::vec3 new_extent2 = extent;
+    new_position2[axis] += (extent[axis] + delta) * 0.5; 
+    new_extent2[axis] -= delta;
+    new_extent2[axis] *= .5;
+    
+    return {CgAABB(new_position1, new_extent1), CgAABB(new_position2, new_extent2)};
+}
\ No newline at end of file
diff --git a/CgMath/CgAABB.h b/CgMath/CgAABB.h
new file mode 100644
index 0000000000000000000000000000000000000000..83fa063751066a5dd0f7ac87a05fa6b89d10cff9
--- /dev/null
+++ b/CgMath/CgAABB.h
@@ -0,0 +1,19 @@
+#ifndef CGAABB_H
+#define CGAABB_H
+
+#include <utility>
+#include <glm/glm.hpp>
+
+// a class representing an axis-aligned bounding box
+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 disjoint 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/CgMath/CgDecomposition.h b/CgMath/CgDecomposition.h
deleted file mode 100644
index 7fdbfc3325068a53d3d5a84b6632babc766d5226..0000000000000000000000000000000000000000
--- a/CgMath/CgDecomposition.h
+++ /dev/null
@@ -1,64 +0,0 @@
-#ifndef CGCOVARIANCE_H
-#define CGCOVARIANCE_H
-
-#include <vector>
-#include <iostream>
-#include "glm/glm.hpp"
-#include "Eigen/Dense"
-#include <CgMath/Eigen/Core>
-using namespace Eigen;
-
-inline glm::vec3 getCentroid(const std::vector<glm::vec3> &points) {
-  glm::vec3 center(0.);
-  for(unsigned int i = 0; i < points.size(); i++) center += points[i];
-  center /= double(points.size());
-  return center;
-}
-
-inline void calculateEigenDecomposition3x3() {
-
-}
-
-inline glm::mat3 getCovarianceMatrix(const std::vector<glm::vec3> &points) {
-    unsigned int size = points.size();
-    Eigen::MatrixXd mat(size, 3);
-    glm::vec3 centroid = getCentroid(points);
-    for(unsigned int i = 0; i < size; i++) {
-        mat(i, 0) = points[i][0] - centroid[0];
-        mat(i, 1) = points[i][1] - centroid[1];
-        mat(i, 2) = points[i][2] - centroid[2];
-    }
-    Eigen::Matrix3d res = (1 / double(size)) * (mat.transpose() * mat);
-
-    glm::mat3 ret(
-      res(0, 0), res(0, 1), res(0, 2), 
-      res(1, 0), res(1, 1), res(1, 2), 
-      res(2, 0), res(2, 1), res(2, 2)
-    );
-    return ret;
-}
-
-inline glm::vec3 estimateNormalFromPCA(const std::vector<glm::vec3> &points) {
-  // calculate covariance matrix
-  unsigned int size = points.size();
-  Eigen::MatrixXd mat(size, 3);
-  glm::vec3 centroid = getCentroid(points);
-  for(unsigned int i = 0; i < size; i++) {
-      mat(i, 0) = points[i][0] - centroid[0];
-      mat(i, 1) = points[i][1] - centroid[1];
-      mat(i, 2) = points[i][2] - centroid[2];
-  }
-  Eigen::Matrix3d res = (1 / double(size)) * (mat.transpose() * mat);
-
-  // get eigendecomposition
-  // the covariance matrix is a symmetric real matrix and therefore 
-  // self-adjoint
-  Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> es(res);
-
-  // The eigenvalues of a selfadjoint matrix are always real.
-  // the eigenvalues are sorted in ascending order
-  Eigen::Vector3d normal = es.eigenvectors().col(0);
-  return glm::vec3(normal(0), normal(1), normal(2));
-}
-
-#endif // CGCOVARIANCE_H
\ No newline at end of file
diff --git a/CgMath/CgPointMath.h b/CgMath/CgPointMath.h
new file mode 100644
index 0000000000000000000000000000000000000000..b550ebb5f7c9f4e5260ba04c30bd5045ed404d33
--- /dev/null
+++ b/CgMath/CgPointMath.h
@@ -0,0 +1,176 @@
+#ifndef CGPOINTMATH_H
+#define CGPOINTMATH_H
+
+#include <vector>
+#include <algorithm>
+#include <limits>
+
+#include <glm/glm.hpp>
+#include <glm/gtx/norm.hpp>
+#include <glm/gtc/matrix_transform.hpp>
+#include "Eigen/Dense"
+#include <CgMath/Eigen/Core>
+using namespace Eigen;
+
+#include "CgAABB.h"
+
+/*
+* simple geometric calculations used in multiple places
+*/
+class CgPointMath {
+    public:
+        static float sqrPointRayDist(glm::vec3 point, glm::vec3 origin, glm::vec3 direction);
+        static float sqrBoxRayDist(CgAABB box, glm::vec3 origin, glm::vec3 direction, glm::vec3 inv_direction);
+        static glm::vec3 calculateCentroid(const std::vector<glm::vec3> &points);
+        static CgAABB calculateAABB(const std::vector<glm::vec3>& vertices);
+        static glm::mat3 getCovarianceMatrix(const std::vector<glm::vec3> &points);
+        static glm::vec3 estimateNormalFromPCA(const std::vector<glm::vec3> &points);
+        static glm::vec3 getPerpendicularVector(glm::vec3 arg);
+    private:
+        CgPointMath(){}
+};
+
+inline float CgPointMath::sqrPointRayDist(glm::vec3 point, glm::vec3 origin, glm::vec3 direction)  {
+    glm::vec3 delta = point - origin;
+    float t = std::max(0.0f, glm::dot(direction, delta) / glm::length2(direction));
+    return glm::length2(direction * t - delta);
+}
+
+inline float CgPointMath::sqrBoxRayDist(CgAABB box, glm::vec3 origin, glm::vec3 direction, glm::vec3 inv_direction) {
+    glm::vec3 bmin = box.position - box.extent;
+    glm::vec3 bmax = box.position + box.extent;
+
+    // calculate distance along the ray for intersection with each
+    // of the boxes planes
+    float tx1 = (bmin.x - origin.x) * inv_direction.x;
+    float tx2 = (bmax.x - origin.x) * inv_direction.x;
+    float ty1 = (bmin.y - origin.y) * inv_direction.y;
+    float ty2 = (bmax.y - origin.y) * inv_direction.y;
+    float tz1 = (bmin.z - origin.z) * inv_direction.z;
+    float tz2 = (bmax.z - origin.z) * inv_direction.z;
+
+    // figure out which of the planes is closest
+    float minX = std::min(tx1, tx2);
+    float maxX = std::max(tx1, tx2);
+    float minY = std::min(ty1, ty2);
+    float maxY = std::max(ty1, ty2);
+    float minZ = std::min(tz1, tz2);
+    float maxZ = std::max(tz1, tz2);
+
+    float p1 = std::max(minX, std::max(minY, minZ));
+    float p2 = std::min(maxX, std::min(maxY, maxZ));
+
+    // clamp in ray direction
+    p1 = std::max(0.0f, p1);
+    p2 = std::max(0.0f, p2);
+
+    // closest point on box
+    float boxX = (origin.x + direction.x * p1 + origin.x + direction.x * p2) / 2;
+    float boxY = (origin.y + direction.y * p1 + origin.y + direction.y * p2) / 2;
+    float boxZ = (origin.z + direction.z * p1 + origin.z + direction.z * p2) / 2;
+
+    boxX = std::max(std::min(boxX, bmax.x), bmin.x);
+    boxY = std::max(std::min(boxY, bmax.y), bmin.y);
+    boxZ = std::max(std::min(boxZ, bmax.z), bmin.z);
+
+    // closest point on ray
+    float t = (boxX - origin.x) * direction.x + (boxY - origin.y) * direction.y + (boxZ - origin.z) * direction.z;
+    t =  t / (direction.x * direction.x + direction.y * direction.y + direction.z * direction.z);
+    t = std::max(0.0f, t);
+    
+    float rayX = origin.x + direction.x * t;
+    float rayY = origin.y + direction.y * t;
+    float rayZ = origin.z + direction.z * t;
+
+    // deltas between ray point and box point
+    float dx = rayX - boxX;
+    float dy = rayY - boxY;
+    float dz = rayZ - boxZ;
+
+    // squared dist
+    return dx * dx + dy * dy + dz * dz;
+}
+
+inline glm::vec3 CgPointMath::calculateCentroid(const std::vector<glm::vec3> &points) {
+    glm::vec3 center(0.);
+    for(unsigned int i = 0; i < points.size(); i++) center += points[i];
+    center /= double(points.size());
+    return center;
+}
+
+inline CgAABB CgPointMath::calculateAABB(const std::vector<glm::vec3>& vertices) {
+    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();
+
+    for(unsigned int i=0; i < vertices.size(); i++) {
+        glm::vec3 curr = vertices[i];
+        max_x = max_x < curr[0] ? curr[0] : max_x;
+        min_x = min_x > curr[0] ? curr[0] : min_x;
+        max_y = max_y < curr[1] ? curr[1] : max_y;
+        min_y = min_y > curr[1] ? curr[1] : min_y;
+        max_z = max_z < curr[2] ? curr[2] : max_z;
+        min_z = min_z > curr[2] ? curr[2] : min_z;
+    }
+
+    return CgAABB(
+        glm::vec3((max_x + min_x) * 0.5, (max_y + min_y) * 0.5, (max_z + min_z) * 0.5),
+        glm::vec3((max_x - min_x) * 0.5, (max_y - min_y) * 0.5, (max_z - min_z) * 0.5)
+    );
+}
+
+inline glm::mat3 CgPointMath::getCovarianceMatrix(const std::vector<glm::vec3> &points) {
+    unsigned int size = points.size();
+    Eigen::MatrixXd mat(size, 3);
+    glm::vec3 centroid = calculateCentroid(points);
+    for(unsigned int i = 0; i < size; i++) {
+        mat(i, 0) = points[i][0] - centroid[0];
+        mat(i, 1) = points[i][1] - centroid[1];
+        mat(i, 2) = points[i][2] - centroid[2];
+    }
+    Eigen::Matrix3d res = (1 / double(size)) * (mat.transpose() * mat);
+
+    glm::mat3 ret(
+    res(0, 0), res(0, 1), res(0, 2), 
+    res(1, 0), res(1, 1), res(1, 2), 
+    res(2, 0), res(2, 1), res(2, 2)
+    );
+    return ret;
+}
+
+inline glm::vec3 CgPointMath::estimateNormalFromPCA(const std::vector<glm::vec3> &points) {
+    // calculate covariance matrix
+    unsigned int size = points.size();
+    Eigen::MatrixXd mat(size, 3);
+    glm::vec3 centroid = calculateCentroid(points);
+    for(unsigned int i = 0; i < size; i++) {
+        mat(i, 0) = points[i][0] - centroid[0];
+        mat(i, 1) = points[i][1] - centroid[1];
+        mat(i, 2) = points[i][2] - centroid[2];
+    }
+    Eigen::Matrix3d res = (1 / double(size)) * (mat.transpose() * mat);
+
+    // get eigendecomposition
+    // the covariance matrix is a symmetric real matrix and therefore 
+    // self-adjoint
+    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> es(res);
+
+    // The eigenvalues of a selfadjoint matrix are always real.
+    // the eigenvalues are sorted in ascending order
+    Eigen::Vector3d normal = es.eigenvectors().col(0);
+    return glm::vec3(normal(0), normal(1), normal(2));
+}
+
+// calculates an arbitrary verctor perpendicular to the given one
+inline glm::vec3 CgPointMath::getPerpendicularVector(glm::vec3 arg) {
+  if((arg.x == 0.0) && (arg.y == 0.0)) {
+    if(arg.z == 0.0) return glm::vec3(0.);
+    return glm::vec3(0.0, 1.0, 0.0);
+  }
+  return glm::normalize(glm::vec3(-arg.y, arg.x, 0.0));
+}
+
+#endif // CGPOINTMATH_H
\ No newline at end of file
diff --git a/CgMath/CgPointWrangler.h b/CgMath/CgPointWrangler.h
new file mode 100644
index 0000000000000000000000000000000000000000..d06f6fa79116ea5f9855dca232ad016bf07475e3
--- /dev/null
+++ b/CgMath/CgPointWrangler.h
@@ -0,0 +1,397 @@
+#ifndef CGPOINTWRANGLER_H
+#define CGPOINTWRANGLER_H
+
+// stl
+#include <vector>
+#include <queue>
+#include <limits>
+#include <string>
+#include <algorithm>
+#include <iostream>
+#include <numeric>
+#include <utility>
+#include <functional>
+#include <tuple>
+#include <thread>
+
+// glm & eigen
+#include <glm/glm.hpp>
+#include <glm/gtx/norm.hpp>
+#include <glm/gtc/matrix_transform.hpp>
+
+// this
+#include "CgAABB.h"
+#include "CgPointMath.h"
+
+/*
+* more complex point cloud algorithms & data structures
+*/
+class CgPointWrangler {
+    public:
+
+        // calculate offset between begin and half-way pointer between begin and end
+        // requires end >= begin
+        static unsigned int halfSize(const glm::vec3* begin, const glm::vec3* end);
+    
+        // rearrange the vec3 between begin and end so they form a kd-tree
+        static void buildKdTree(
+            glm::vec3* begin, glm::vec3* end, 
+            int depth
+        );
+
+        // traverse a kd-tree between begin and end in BFO and call a function on each node
+        static void traverseBFO(
+            glm::vec3* begin, glm::vec3* end, 
+            std::function<void(glm::vec3*, const unsigned int)> fn
+        );
+
+        // get the point that's closest to the ray in the kd-tree between begin and end
+        static int getClosestPointToRay(
+            const glm::vec3* begin, const glm::vec3* end,
+            const CgAABB aabb,
+            glm::vec3 origin, glm::vec3 direction
+        );
+
+        // naive estimation of normals
+        static std::vector<glm::vec3> estimateNormals(unsigned int neighbourhood_size, const std::vector<glm::vec3>& vertices);
+
+        static std::vector<unsigned int> getNearestNeighborsSlow(
+            unsigned int current_point, 
+            unsigned int k,
+            const std::vector<glm::vec3>& vertices
+        );
+
+        static std::vector<unsigned int> getNearestNeighborsFast(
+            unsigned int current_point, 
+            unsigned int k, 
+            const std::vector<glm::vec3>& vertices,
+            CgAABB aabb
+        );
+
+    private:
+        CgPointWrangler(){}
+};
+
+inline unsigned int CgPointWrangler::halfSize(const glm::vec3* begin, const glm::vec3* end) {
+    return (end - begin) / 2;
+}
+
+// rearrange the vec3 between begin and end (exclusive) so they form a kd-tree
+// requires begin < end
+inline void CgPointWrangler::buildKdTree(glm::vec3* begin, glm::vec3* end, int depth) {
+    unsigned int half_size = halfSize(begin, end);
+
+    // none or one element
+    if(half_size == 0) return;
+
+    // partially sort and ensure median element 
+    // is in the middle (runs in O(n))
+    // about 4x faster on tyra.obj than std::sort
+    std::nth_element(begin, begin + half_size, end, [=](glm::vec3 a, glm::vec3 b) {
+        return a[depth % 3] < b[depth % 3];
+    });
+
+    // split array in two (excluding median) and recurse
+    buildKdTree(begin, begin + half_size, depth + 1);
+    buildKdTree(begin + half_size + 1, end, depth + 1);
+}
+
+inline void CgPointWrangler::traverseBFO(
+    glm::vec3* begin, 
+    glm::vec3* end, 
+    std::function<void(glm::vec3*, const unsigned int)> fn
+) {
+    // queue of pairs of vec3 pointers
+    // (range of points in the current node and its children)
+    std::queue<std::pair<glm::vec3*, glm::vec3*>> q;
+    if(begin == end) return;
+    unsigned int depth = 0;
+    q.push({begin, end});
+    while(!q.empty()) {
+        int level_count = q.size();
+        for(int i = 0; i < level_count; i++) {
+            auto curr = q.front();
+            q.pop();
+            glm::vec3* begin = curr.first;
+            glm::vec3* end = curr.second;
+            unsigned int half_size = halfSize(begin, end);
+            glm::vec3* curr_pointer = begin + half_size;
+            fn(curr_pointer, depth);
+            if(end - begin <= 1) continue;
+            // push left & right half of array to reconstruct splits
+            if(begin <= begin + half_size) q.push({begin, begin + half_size});
+            if(begin + half_size < end) q.push({begin + half_size + 1, end});
+        }
+        depth += 1;
+    }
+}
+
+inline int CgPointWrangler::getClosestPointToRay(
+    const glm::vec3* begin, 
+    const glm::vec3* end,
+    CgAABB aabb,
+    glm::vec3 origin, 
+    glm::vec3 direction
+) {
+  if(begin == end) return -1;
+
+  glm::vec3 inv_direction = glm::vec3(1.0 / direction.x, 1.0 / direction.y, 1.0 / direction.z);
+  // sensible starting point
+  float erg_sqr_dist = std::numeric_limits<float>::infinity();
+  int erg = -1;
+
+  // stack of traversed nodes (start index, end index, aabb, depth)
+  std::vector<std::tuple<const glm::vec3*, const glm::vec3*, CgAABB, unsigned int>> nodes;
+  nodes.push_back({begin, end, aabb, 0});
+  
+  // traverse the tree as long as there's something to do.
+  while(!nodes.empty()) {
+    // unpack current node & pop it
+    const glm::vec3* cur_begin;
+    const glm::vec3* cur_end;
+    CgAABB cur_aabb;
+    unsigned int cur_depth;
+    std::tie(cur_begin, cur_end, cur_aabb, cur_depth) = nodes.back();
+    auto half_size = halfSize(cur_begin, cur_end);
+    nodes.pop_back();
+
+    glm::vec3 cur_point = *(cur_begin + half_size);
+    // insert current point into set if it's closer to ray than last one
+    float cand_dist = CgPointMath::sqrPointRayDist(cur_point, origin, direction);
+    if(cand_dist < erg_sqr_dist) {
+      erg_sqr_dist = cand_dist;
+      erg = cur_begin + half_size - begin;
+    }
+
+    // the aabb of the sub-nodes
+    unsigned int axis = cur_depth % 3;
+    CgAABB lower_aabb;
+    CgAABB higher_aabb;
+    std::tie(lower_aabb, higher_aabb) = cur_aabb.split(axis, cur_point[axis]);
+    // only push children if their aabb is closer to ray than 
+    // the closest point we found until now and there are
+    // points in them
+    if (
+        half_size != 0
+        && CgPointMath::sqrBoxRayDist(lower_aabb, origin, direction, inv_direction) < erg_sqr_dist
+    ) {
+        nodes.push_back({cur_begin, cur_begin + half_size, lower_aabb, cur_depth + 1});
+    }
+    if (
+        cur_begin + half_size < cur_end
+        && CgPointMath::sqrBoxRayDist(higher_aabb, origin, direction, inv_direction) < erg_sqr_dist
+    ) {
+        nodes.push_back({cur_begin + half_size + 1, cur_end, higher_aabb, cur_depth + 1});
+    }
+  }
+
+  return erg;
+}
+
+// naive estimation of normals
+inline std::vector<glm::vec3> CgPointWrangler::estimateNormals(unsigned int neighbourhood_size, const std::vector<glm::vec3>& vertices) {
+    std::vector<glm::vec3> normals;
+    normals.resize(vertices.size());
+    CgAABB aabb = CgPointMath::calculateAABB(vertices);
+
+    // find point with highest y-value
+    std::vector<bool> visited;
+    visited.resize(vertices.size());
+    float max_y = -std::numeric_limits<float>::infinity();
+    unsigned int max_y_index = 0;
+    for(int i = 0; i < vertices.size(); i++) {
+        visited[i] = false;
+        float curr = vertices[i].y;
+        if(max_y < curr) {
+            max_y = curr;
+            max_y_index = i;
+        }
+    }
+    visited[max_y_index] = true;
+
+    // used as a stack for traversal order
+    std::vector<unsigned int> spanning_tree;
+    // the "best" parent found for each vertex
+    std::vector<std::pair<unsigned int, double>> scores;
+    scores.resize(vertices.size());
+    // all k-neighbourhoods
+    std::vector<std::vector<unsigned int>> neighbourhoods;
+    neighbourhoods.resize(vertices.size());
+
+    // record all neighbourhoods and calculate (possibly flipped) normals
+    unsigned int thread_count = std::thread::hardware_concurrency();
+    std::vector<std::thread> threads;
+    threads.reserve(thread_count);
+
+    auto worker = [&](
+        const unsigned int start_index, 
+        const unsigned int count
+    ) {
+        // thread-local buffer for points of one neighbourhood
+        std::vector<glm::vec3> nh(neighbourhood_size);
+        for(unsigned int i = start_index; i < start_index + count; i++) {
+            neighbourhoods[i] = getNearestNeighborsFast(i, neighbourhood_size, vertices, aabb);
+            for(int j = 0; j < neighbourhood_size; j++) {
+                nh[j] = vertices[neighbourhoods[i][j]];
+            }
+            normals[i] = CgPointMath::estimateNormalFromPCA(nh);
+        }
+    };
+
+    unsigned int slice_size = vertices.size() / thread_count;
+    unsigned int slice_rest = vertices.size() % thread_count;
+    unsigned int t = 0;
+    for(t = 0; t < thread_count - 1; t++) threads.push_back(std::thread(worker, t * slice_size, slice_size));
+    threads.push_back(std::thread(worker, t * slice_size, slice_size + slice_rest));
+    for(t = 0; t < thread_count; t++) threads[t].join();
+
+    auto visit = [&](unsigned int curr, glm::vec3 parent_normal, glm::vec3 parent_pos) {
+        // flip current normal to conform to parent
+        glm::vec3 normal = normals[curr];
+        normal = glm::dot(normal, parent_normal) > 0.0 ? normal : -normal;
+        normals[curr] = normal;
+
+        // retrieve neighbourhood indices
+        auto curr_nh_ind = neighbourhoods[curr];
+
+        // find out which of the neighbours we didn't process yet
+        // and score them to find out in which order to process
+        for(unsigned int i = 0; i < neighbourhood_size; i++) {
+            unsigned int index = curr_nh_ind[i];
+            // score by how far away it is from the local tangent plane
+            // and how aligned the normals are.
+            double alignment = 1.0 - std::abs(glm::dot(normals[index], normal));
+            double distance = std::abs(glm::dot(vertices[index] - vertices[curr], normal));
+            double score = distance * alignment;
+            if(!visited[index]) {
+                scores[index].first = curr;
+                scores[index].second = score;
+                spanning_tree.push_back(index);
+                visited[index] = true;
+            } 
+            if(scores[index].second > score) {
+                // update current neighbors parent to current node
+                scores[index].first = curr;
+                scores[index].second = score;
+            }
+        }
+    };
+
+    // get started by aligning the normal of 
+    // the highest point away from the centroid
+    glm::vec3 centroid = CgPointMath::calculateCentroid(vertices);
+    visit(max_y_index, vertices[max_y_index] - centroid, vertices[max_y_index]);
+
+    // start processing
+    while(!spanning_tree.empty()) {
+        unsigned int curr = spanning_tree.back();
+        unsigned int parent = scores[curr].first;
+        spanning_tree.pop_back();
+        glm::vec3 parent_normal = normals[parent];
+        glm::vec3 parent_pos = vertices[parent];
+        visit(curr, parent_normal, parent_pos);
+    }
+
+    return normals;
+}
+
+inline std::vector<unsigned int> CgPointWrangler::getNearestNeighborsSlow(
+    unsigned int current_point, 
+    unsigned int k, 
+    const std::vector<glm::vec3>& vertices
+) {
+    glm::vec3 q = vertices[current_point];
+
+    std::vector<std::pair<double,int>> distances;
+
+    for(unsigned int i = 0; i < vertices.size(); i++) {
+        double dist = glm::distance(vertices[i],q);
+        distances.push_back(std::make_pair(dist,i));
+    }
+
+        std::sort(distances.begin(), distances.end());
+
+        std::vector<unsigned int> erg;
+
+    for(unsigned int i = 0; i < k; i++) {
+        erg.push_back(distances[i].second);
+    }
+
+    return erg;
+}
+
+inline std::vector<unsigned int> CgPointWrangler::getNearestNeighborsFast(
+    unsigned int current_point, 
+    unsigned int k,
+    const std::vector<glm::vec3>& vertices,
+    CgAABB aabb
+) {
+    std::vector<unsigned int> erg(k);
+    if(k == 0 || vertices.size() < k) return erg;
+    // pair of distance + index into backing vec3 array 
+    using Pair = std::pair<double, unsigned int>;
+
+    glm::vec3 query = vertices[current_point];
+    const glm::vec3* first = vertices.data();
+    const glm::vec3* last = vertices.data() + vertices.size();
+    auto cmp = [&](const Pair &left, const Pair &right) { return (left.first < right.first);};
+
+    // using an ordered set for simplicitly, this may be slow
+    // because we frequently erase elements, causing the rb-tree
+    // to be rebuilt.
+    std::vector<Pair> candidates;
+    candidates.reserve(k + 1);
+
+    // stack of traversed nodes (start index, end index, aabb, depth)
+    std::vector<std::tuple<const glm::vec3*, const glm::vec3*, CgAABB, unsigned int>> nodes;
+    nodes.push_back({first, last, aabb, 0});
+    
+    // traverse the tree as long as there's something to do.
+    while(!nodes.empty()) {
+        // unpack current node & pop it
+        const glm::vec3* begin;
+        const glm::vec3* end;
+        CgAABB aabb;
+        unsigned int depth;
+        std::tie(begin, end, aabb, depth) = nodes.back();
+        nodes.pop_back();
+        unsigned int half_size = halfSize(begin, end);
+        glm::vec3 curr_point = *(begin + half_size);
+        // insert current point distance & index into heap
+        candidates.push_back({glm::distance(query, curr_point), begin + half_size - first});
+        std::push_heap(candidates.begin(), candidates.end(), cmp);
+        if(candidates.size() > k) {
+            std::pop_heap(candidates.begin(), candidates.end(), cmp);
+            candidates.pop_back();
+        }
+
+        // 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 children if its aabb intersects with the 
+        // neighbourhood bounding sphere
+        double neighbourhood_radius = candidates[0].first;
+        if(
+            begin < begin + half_size
+            && (lower_aabb.position[axis] + lower_aabb.extent[axis] >= query[axis] - neighbourhood_radius || k > candidates.size())
+        ) {
+            nodes.push_back({begin, begin + half_size, lower_aabb, depth + 1});
+        }
+        if(
+            begin + half_size < end
+            && (higher_aabb.position[axis] - higher_aabb.extent[axis] <= query[axis] + neighbourhood_radius || k > candidates.size())
+        ) {
+            nodes.push_back({begin + half_size + 1, end, higher_aabb, depth + 1});
+        }
+    }
+
+    // copy all nearest neighbour indices we found into
+    // a vector to return
+    for(unsigned int i = 0; i < k; i++) erg[i] = candidates[i].second;
+
+    return erg;
+}
+
+#endif // CGPOINTWRANGLER_H
\ No newline at end of file
diff --git a/CgQtViewer/CgQtGui.cpp b/CgQtViewer/CgQtGui.cpp
index edea8d18fd523afee01558ac69f5f28d0f8eca37..6fe4600521f5c52322fd3b8d78879a72cf942886 100644
--- a/CgQtViewer/CgQtGui.cpp
+++ b/CgQtViewer/CgQtGui.cpp
@@ -160,7 +160,7 @@ void CgQtGui::createKdTreeOptionPanel(QWidget* parent) {
     tab1_control->addWidget(max_kd_depth_spinbox);
     max_kd_depth_spinbox->setMinimum(1);
     max_kd_depth_spinbox->setMaximum(50);
-    max_kd_depth_spinbox->setValue(7);
+    max_kd_depth_spinbox->setValue(3);
     connect(max_kd_depth_spinbox, SIGNAL(valueChanged(int) ), this, SLOT(slotKdDepthSpinboxChanged()) );
     tab1_control->addWidget(max_kd_depth_spinbox);
     max_kd_depth_spinbox->setAlignment(Qt::AlignTop);
@@ -210,7 +210,6 @@ void CgQtGui::createSimplifyOptionPanel(QWidget* parent) {
     max_diameter_of_cluster_spinbox->setMaximum(10.0);
     max_diameter_of_cluster_spinbox->setValue(0.05);
     max_diameter_of_cluster_spinbox->setSingleStep(0.05);
-    tab2_control->addWidget(max_diameter_of_cluster_spinbox);
     max_diameter_of_cluster_spinbox->setAlignment(Qt::AlignTop);
     max_diameter_of_cluster_spinbox->setSizePolicy(QSizePolicy::Ignored, QSizePolicy::Fixed);
 
@@ -226,31 +225,42 @@ void CgQtGui::createSimplifyOptionPanel(QWidget* parent) {
     max_variance_of_cluster_spinbox->setMaximum(5.0);
     max_variance_of_cluster_spinbox->setValue(0.05);
     max_variance_of_cluster_spinbox->setSingleStep(0.05);
-    tab2_control->addWidget(max_variance_of_cluster_spinbox);
     max_variance_of_cluster_spinbox->setAlignment(Qt::AlignTop);
     max_variance_of_cluster_spinbox->setSizePolicy(QSizePolicy::Ignored, QSizePolicy::Fixed);
 
+    number_of_points_checkbox = new QCheckBox("use number of points");
+    number_of_points_checkbox->setCheckable(true);
+    number_of_points_checkbox->setChecked(false);
+    tab2_control->addWidget(number_of_points_checkbox);
+
+    diameter_of_cluster_checkbox = new QCheckBox("use cluster diameter");
+    diameter_of_cluster_checkbox->setCheckable(true);
+    diameter_of_cluster_checkbox->setChecked(false);
+    tab2_control->addWidget(diameter_of_cluster_checkbox);
+
+    variance_of_cluster_checkbox = new QCheckBox("use variance");
+    variance_of_cluster_checkbox->setCheckable(true);
+    variance_of_cluster_checkbox->setChecked(false);
+    tab2_control->addWidget(variance_of_cluster_checkbox);
+
     /* select clustering type button group */
 
-    QGroupBox* myGroupBox = new QGroupBox("Clustering Type: ");
+    QGroupBox* myGroupBox = new QGroupBox("Clustering type: ");
 
     clustering_type_buttongroup = new QButtonGroup(subBox);
     clustering_type_buttongroup->setExclusive(true);
 
-    QRadioButton* radiobutton1 = new QRadioButton( "Incremental w/ n");
-    QRadioButton* radiobutton2 = new QRadioButton( "Incremental w/ d");
-    QRadioButton* radiobutton3 = new QRadioButton( "Hierarchical");
+    QRadioButton* radiobutton1 = new QRadioButton( "Incremental");
+    QRadioButton* radiobutton2 = new QRadioButton( "Hierarchical");
 
     radiobutton1->setChecked(true);
 
-    clustering_type_buttongroup->addButton(radiobutton1, int(Cg::CgIncrementalNumber));
-    clustering_type_buttongroup->addButton(radiobutton2, int(Cg::CgIncrementalDiameter));
-    clustering_type_buttongroup->addButton(radiobutton3, int(Cg::CgHierarchical));
+    clustering_type_buttongroup->addButton(radiobutton1, int(Cg::CgIncremental));
+    clustering_type_buttongroup->addButton(radiobutton2, int(Cg::CgHierarchical));
 
     QVBoxLayout *vbox = new QVBoxLayout;
     vbox->addWidget(radiobutton1);
     vbox->addWidget(radiobutton2);
-    vbox->addWidget(radiobutton3);
     vbox->addStretch(1);
     myGroupBox->setLayout(vbox);
     subBox->addWidget(myGroupBox);
@@ -270,12 +280,24 @@ void CgQtGui::createSimplifyOptionPanel(QWidget* parent) {
 }
 
 void CgQtGui::sendSimplifyEvent(bool reset) {
+    float number_of_points = number_of_points_checkbox->isChecked()
+        ? max_number_of_points_spinbox->value()
+        : std::numeric_limits<float>::infinity();
+
+    double diameter = diameter_of_cluster_checkbox->isChecked()
+        ? max_diameter_of_cluster_spinbox->value()
+        : std::numeric_limits<float>::infinity();
+
+    double variance = variance_of_cluster_checkbox->isChecked()
+        ? max_variance_of_cluster_spinbox->value()
+        : std::numeric_limits<float>::infinity();
+
     CgBaseEvent* e = new CgSimplifyEvent(
         Cg::CgSimplifyEvent,
         reset,
-        max_number_of_points_spinbox->value(),
-        max_diameter_of_cluster_spinbox->value(),
-        max_variance_of_cluster_spinbox->value(),
+        number_of_points,
+        diameter,
+        variance,
         (Cg::ClusteringType) clustering_type_buttongroup->checkedId()
     );
     notifyObserver(e);
diff --git a/CgQtViewer/CgQtGui.h b/CgQtViewer/CgQtGui.h
index c3a159dbf76b71c46ed40b70e7b2f7294d873c05..f39f0225816182bf0fe28e5f7bb5ae25dc745d03 100644
--- a/CgQtViewer/CgQtGui.h
+++ b/CgQtViewer/CgQtGui.h
@@ -42,6 +42,7 @@
 #define WINDOW_H
 
 #include <QWidget>
+#include <limits>
 #include "../CgBase/CgObservable.h"
 
 QT_BEGIN_NAMESPACE
@@ -108,6 +109,10 @@ private:
     QDoubleSpinBox* max_variance_of_cluster_spinbox;
     QButtonGroup* clustering_type_buttongroup;
 
+    QCheckBox* number_of_points_checkbox;
+    QCheckBox* diameter_of_cluster_checkbox;
+    QCheckBox* variance_of_cluster_checkbox;
+
     bool m_use_spats;
     bool m_show_pickray;
 
diff --git a/CgSceneGraph/CgAABB.cpp b/CgSceneGraph/CgAABB.cpp
deleted file mode 100644
index b2e3ad0fda65274414be0df7e8da9c48dd40c7e7..0000000000000000000000000000000000000000
--- a/CgSceneGraph/CgAABB.cpp
+++ /dev/null
@@ -1,114 +0,0 @@
-#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)};
-}
-
-// idea stolen from https://tavianator.com/2011/ray_box.html
-bool CgAABB::intersectsRay(glm::vec3 origin, glm::vec3 inv_direction) {
-    glm::vec3 bmin = position - extent;
-    glm::vec3 bmax = position + extent;
-
-    // t-val for intersection with each plane
-    glm::vec3 tlo = (bmin - origin) * inv_direction;
-    glm::vec3 thi = (bmax - origin) * inv_direction;
-
-    // t-vals that intersect first for each dimension
-    glm::vec3 tmin = glm::vec3(
-        std::min(tlo.x, thi.x),
-        std::min(tlo.y, thi.y),
-        std::min(tlo.z, thi.z)
-    );
-    // t-vals that intersect last for each dimension
-    glm::vec3 tmax = glm::vec3(
-        std::max(tlo.x, thi.x),
-        std::max(tlo.y, thi.y),
-        std::max(tlo.z, thi.z)
-    );
-
-    return glm::compMax(tmin) >= glm::compMin(tmax);
-}
-
-float CgAABB::sqrDistanceToRay(glm::vec3 origin, glm::vec3 direction, glm::vec3 inv_direction) {
-    glm::vec3 bmin = position - extent;
-    glm::vec3 bmax = position + extent;
-
-    // calculate distance along the ray for intersection with each
-    // of the boxes planes
-    float tx1 = (bmin.x - origin.y) * inv_direction.x;
-    float tx2 = (bmax.x - origin.y) * inv_direction.x;
-    float ty1 = (bmin.y - origin.y) * inv_direction.y;
-    float ty2 = (bmax.y - origin.y) * inv_direction.y;
-    float tz1 = (bmin.z - origin.z) * inv_direction.z;
-    float tz2 = (bmax.z - origin.z) * inv_direction.z;
-
-    // figure out which of the planes is closest
-    float minX = std::min(tx1, tx2);
-    float maxX = std::max(tx1, tx2);
-    float minY = std::min(ty1, ty2);
-    float maxY = std::max(ty1, ty2);
-    float minZ = std::min(tz1, tz2);
-    float maxZ = std::max(tz1, tz2);
-
-    float p1 = std::max(minX, std::max(minY, minZ));
-    float p2 = std::min(maxX, std::min(maxY, maxZ));
-
-    // clamp in ray direction
-    p1 = std::max(0.0f, p1);
-    p2 = std::max(0.0f, p2);
-
-    // closest point on box
-    float boxX = (origin.x + direction.x * p1 + origin.x + direction.x * p2) / 2;
-    float boxY = (origin.y + direction.y * p1 + origin.y + direction.y * p2) / 2;
-    float boxZ = (origin.z + direction.z * p1 + origin.z + direction.z * p2) / 2;
-
-    boxX = std::max(std::min(boxX, bmax.x), bmin.x);
-    boxY = std::max(std::min(boxY, bmax.y), bmin.y);
-    boxZ = std::max(std::min(boxZ, bmax.z), bmin.z);
-
-    // closest point on ray
-    float t = (boxX - origin.x) * direction.x + (boxY - origin.y) * direction.y + (boxZ - origin.z) * direction.z;
-    t =  t / (direction.x * direction.x + direction.y * direction.y + direction.z * direction.z);
-    t = std::max(0.0f, t);
-    
-    float rayX = origin.x + direction.x * t;
-    float rayY = origin.y + direction.y * t;
-    float rayZ = origin.z + direction.z * t;
-
-    // deltas between ray point and box point
-    float dx = rayX - boxX;
-    float dy = rayY - boxY;
-    float dz = rayZ - boxZ;
-
-    // squared dist
-    return dx * dx + dy * dy + dz * dz;
-}
-
-bool CgAABB::intersectsPoint(glm::vec3 point) {
-    return true 
-    && (point.x <= position.x + extent.x && point.x >= position.x + extent.x)
-    && (point.y <= position.y + extent.y && point.y >= position.y + extent.y)
-    && (point.z <= position.z + extent.z && point.z >= position.z + extent.z);
-}
\ No newline at end of file
diff --git a/CgSceneGraph/CgAABB.h b/CgSceneGraph/CgAABB.h
deleted file mode 100644
index 7737180be8f4d1bae190ea2c84a47c945dda6ebd..0000000000000000000000000000000000000000
--- a/CgSceneGraph/CgAABB.h
+++ /dev/null
@@ -1,27 +0,0 @@
-#ifndef CGAABB_H
-#define CGAABB_H
-
-#include <utility>
-#include <algorithm>
-#include <glm/glm.hpp>
-#include <glm/gtx/component_wise.hpp>
-
-// a class representing an axis-aligned bounding box
-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);
-    // ray-box collision
-    bool intersectsRay(glm::vec3 origin, glm::vec3 inv_direction);
-    // point-box collision
-    bool intersectsPoint(glm::vec3 point);
-    // return the squared distance between the box and a ray
-    float sqrDistanceToRay(glm::vec3 origin, glm::vec3 direction, glm::vec3 inv_direction);
-};
-
-
-#endif // CGAABB_H
\ No newline at end of file
diff --git a/CgSceneGraph/CgPointCloud.cpp b/CgSceneGraph/CgPointCloud.cpp
index 629b4caa6f165390f9791f3b17823f5f3c5a5cbe..a7cec1c16e21e43c2de0bf4ceaba3b78c8c48da0 100644
--- a/CgSceneGraph/CgPointCloud.cpp
+++ b/CgSceneGraph/CgPointCloud.cpp
@@ -1,23 +1,32 @@
 #include "CgPointCloud.h"
 
-// squared distance between a point and a ray
-float sqrPointRayDist(glm::vec3 point, glm::vec3 origin, glm::vec3 direction) {
-  glm::vec3 delta = point - origin;
-  float t = std::max(0.0f, glm::dot(direction, delta) / glm::length2(direction));
-  return glm::length2(direction * t - delta);
-}
-
-CgPointCloud::CgPointCloud(): CgPointCloud::CgPointCloud(51) {}
-
-CgPointCloud::CgPointCloud(int id):
-m_type(Cg::PointCloud),
-m_id(id) {
-    m_vertices.push_back(glm::vec3(0.0,0.0,0.0));
-    m_vertex_normals.push_back(glm::vec3(0.0,0.0,1.0));
-    m_vertex_colors.push_back(glm::vec3(0.0,0.0,1.0));
-    calculateSplatOrientations();
-}
-
+CgPointCloud::CgPointCloud(
+  std::vector<glm::vec3> vertices,
+  std::vector<glm::vec3> vertex_normals,
+  std::vector<glm::vec3> vertex_colors,
+  std::vector<glm::mat4> splat_orientations,
+  std::vector<glm::vec2> splat_scaling,
+  std::vector<unsigned int> splat_indices
+): CgPointCloud::CgPointCloud(51) {
+    m_vertices = vertices;
+    m_vertex_normals = vertex_normals;
+    m_vertex_colors = vertex_colors;
+    m_splat_orientations = splat_orientations;
+    m_splat_scaling = splat_scaling;
+    m_splat_indices = splat_indices;
+
+    m_centroid = CgPointMath::calculateCentroid(m_vertices);
+    m_aabb = CgPointMath::calculateAABB(m_vertices);
+}
+
+CgPointCloud::CgPointCloud(const CgPointCloud& other) : CgPointCloud::CgPointCloud(
+  other.m_vertices,
+  other.m_vertex_normals,
+  other.m_vertex_colors,
+  other.m_splat_orientations,
+  other.m_splat_scaling,
+  other.m_splat_indices
+){}
 
 CgPointCloud::~CgPointCloud() {
     m_vertices.clear();
@@ -26,7 +35,6 @@ CgPointCloud::~CgPointCloud() {
     m_splat_indices.clear();
 }
 
-
 void CgPointCloud::calculateSplatOrientations() {
   // calculate local coordinate system for splats (arbitrary orientation of ellipse in plane)
   // replace this if you have the real coordinate system, use up vector = y-Axis of your local coordinate system instead of getPerpendicularVector(...)
@@ -35,357 +43,21 @@ void CgPointCloud::calculateSplatOrientations() {
   m_splat_scaling.clear();
   m_splat_indices.clear();
 
-  for(unsigned int i=0;i<m_vertices.size();i++)
-    {
-      glm::mat4 lookAt_matrix(glm::lookAt(glm::vec3(m_vertices[i]),glm::vec3(m_vertices[i]-m_vertex_normals[i]),getPerpendicularVector(m_vertex_normals[i])));
-      m_splat_orientations.push_back(lookAt_matrix);
-      m_splat_scaling.push_back(glm::vec2(0.02,0.005));
-
-      // use all points for splatting by default
-      m_splat_indices.push_back(i);
-    }
-
-
-}
-
-void CgPointCloud::init( std::string filename, bool cheat_normals) {
-    m_vertices.clear();
-    m_vertex_normals.clear();
-    m_vertex_colors.clear();
-    m_splat_orientations.clear();
-    m_splat_scaling.clear();
-    m_splat_indices.clear();
-
-    // load obj File
-    ObjLoader loader;
-    loader.load(filename);
-    loader.getPositionData(m_vertices);
-    m_vertex_normals.resize(m_vertices.size());
-    m_vertex_colors.resize(m_vertices.size());
-    std::cout << "loaded " << m_vertices.size() << " points" << std::endl;
-    std::cout << "using " << std::thread::hardware_concurrency() << " threads" << std::endl;
-    m_centroid = getCentroid(m_vertices);
-
-    run_timed("build kd", 1, [&](){buildKdTree(&m_vertices.front(), &(*(m_vertices.end() - 1)), 0);});
-    calculateAABB();
-    run_timed("est. normals", 1, [&](){estimateNormals();});
-    //estimateNormals();
-
-    // calculate local coordinate system for splats 
-    // (arbitrary orientation of ellipse in plane)
-    // replace this if you have the real coordinate 
-    // system, use up vector = y-Axis of your local 
-    // coordinate system instead of getPerpendicularVector(...)
-    //calculateSplatOrientations();
-
-    //add a standard color for each point if lighting turned off
-    //for(unsigned int i = 0; i < m_vertices.size(); i++) {
-    //  m_vertex_colors[i] = glm::vec3(0.0,1.0,0.0);
-    //}
-    resetColors(glm::vec3(0.0, 0.75, 0.0));
-}
-
-void CgPointCloud::buildKdTree(glm::vec3* begin, glm::vec3* end, int depth) {
-    auto half_size = (end - begin) / 2;
-
-    if(begin == end || half_size == 0) {
-        // zero or one element, nothing to do 
-        return;
-    } else { 
-        // more than one element!
-
-        // partially sort and ensure median element 
-        // is in the middle (runs in O(n))
-        // about 4x faster on tyra.obj than std::sort
-        std::nth_element(begin, begin + half_size, end, [=](glm::vec3 a, glm::vec3 b) {
-            return a[depth % 3] < b[depth % 3];
-        });
-
-        // split array in two (excluding median) and recurse
-        buildKdTree(begin, begin + half_size, depth + 1);
-        buildKdTree(begin + half_size + 1, end, depth + 1);
-    }
-}
-
-void CgPointCloud::estimateNormals() {
-  unsigned int neighbourhood_size = 15;
-
-  // find point with highest y-value
-  std::vector<bool> visited;
-  visited.resize(m_vertices.size());
-  float max_y = -std::numeric_limits<float>::infinity();
-  unsigned int max_y_index = 0;
-  for(int i = 0; i < m_vertices.size(); i++) {
-    visited[i] = false;
-    float curr = m_vertices[i].y;
-    if(max_y < curr) {
-      max_y = curr;
-      max_y_index = i;
-    }
-  }
-  visited[max_y_index] = true;
-
-  // used as a stack for traversal order
-  std::vector<unsigned int> spanning_tree;
-  // the "best" parent found for each vertex
-  std::vector<std::pair<unsigned int, double>> scores;
-  scores.resize(m_vertices.size());
-  // all k-neighbourhoods
-  std::vector<std::vector<unsigned int>> neighbourhoods;
-  neighbourhoods.resize(m_vertices.size());
-
-  // record all neighbourhoods and calculate (possibly flipped) normals
-  unsigned int thread_count = std::thread::hardware_concurrency();
-  std::vector<std::thread> threads;
-  threads.reserve(thread_count);
-
-  auto worker = [&](
-    const unsigned int start_index, 
-    const unsigned int count
-  ) {
-    // buffer for points of one neighbourhood
-    std::vector<glm::vec3> nh(neighbourhood_size);
-    for(unsigned int i = start_index; i < start_index + count; i++) {
-      neighbourhoods[i] = getNearestNeighborsFast(i, neighbourhood_size);
-      for(int j = 0; j < neighbourhood_size; j++) {
-        nh[j] = m_vertices[neighbourhoods[i][j]];
-      }
-      m_vertex_normals[i] = estimateNormalFromPCA(nh);
-    }
-  };
-
-  unsigned int slice_size = m_vertices.size() / thread_count;
-  unsigned int slice_rest = m_vertices.size() % thread_count;
-  unsigned int t = 0;
-  for(t = 0; t < thread_count - 1; t++) threads.push_back(std::thread(worker, t * slice_size, slice_size));
-  threads.push_back(std::thread(worker, t * slice_size, slice_size + slice_rest));
-  for(t = 0; t < thread_count; t++) threads[t].join();
-
-  auto visit = [&](unsigned int curr, glm::vec3 parent_normal, glm::vec3 parent_pos) {
-    // flip current normal to conform to parent
-    glm::vec3 normal = m_vertex_normals[curr];
-    normal = glm::dot(normal, parent_normal) > 0.0 ? normal : -normal;
-    m_vertex_normals[curr] = normal;
-
-    // retrieve neighbourhood indices
-    auto curr_nh_ind = neighbourhoods[curr];
-
-    // find out which of the neighbours we didn't process yet
-    // and score them to find out in which order to process
-    for(unsigned int i = 0; i < neighbourhood_size; i++) {
-      unsigned int index = curr_nh_ind[i];
-      // score by how far away it is from the local tangent plane
-      // and how aligned the normals are.
-      double alignment = 1.0 - std::abs(glm::dot(m_vertex_normals[index], normal));
-      double distance = std::abs(glm::dot(m_vertices[index] - m_vertices[curr], normal));
-      double score = distance * alignment;
-      if(!visited[index]) {
-        scores[index].first = curr;
-        scores[index].second = score;
-        spanning_tree.push_back(index);
-        visited[index] = true;
-      } 
-      if(scores[index].second > score) {
-        // update current neighbors parent to current node
-        scores[index].first = curr;
-        scores[index].second = score;
-      }
-    }
-  };
-
-  // get started by aligning the normal of 
-  // the highest point away from the centroid
-  visit(max_y_index, m_vertices[max_y_index] - m_centroid, m_vertices[max_y_index]);
-
-  // start processing
-  while(!spanning_tree.empty()) {
-    unsigned int curr = spanning_tree.back();
-    unsigned int parent = scores[curr].first;
-    spanning_tree.pop_back();
-    glm::vec3 parent_normal = m_vertex_normals[parent];
-    glm::vec3 parent_pos = m_vertices[parent];
-    visit(curr, parent_normal, parent_pos);
-  }
-
-  m_vertex_colors[max_y_index] = glm::vec3(1.0,0,0);
-}
-
-void CgPointCloud::traverseBFO(std::function<void(glm::vec3*, const unsigned int)> fn) {
-    // queue of pairs of vec3 pointers
-    // (range of points in the current node and its children)
-    std::queue<std::pair<glm::vec3*, glm::vec3*>> q;
-    if(m_vertices.empty()) return;
-    unsigned int depth = 0;
-    q.push({&m_vertices.front(), &(*(m_vertices.end() - 1))});
-    while(!q.empty()) {
-        int level_count = q.size();
-        for(int i = 0; i < level_count; i++) {
-            auto curr = q.front();
-            glm::vec3* begin = curr.first;
-            glm::vec3* end = curr.second;
-            unsigned int half_size = (end - begin) / 2;
-            glm::vec3* curr_pointer = begin + half_size;
-            fn(curr_pointer, depth);
-            // push left & right half of array to reconstruct splits
-            if(half_size != 0) q.push({begin, begin + half_size});
-            if(begin + half_size != end) q.push({begin + half_size + 1, end});
-            q.pop();
-        }
-        depth += 1;
-    }
-}
-
-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;
-
   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());
-
-    std::vector<unsigned int> erg;
-
-  for(unsigned int i = 0; i < k; i++) {
-       erg.push_back(distances[i].second);
-  }
-
-  return erg;
-}
-
-std::vector<unsigned int> CgPointCloud::getNearestNeighborsFast(unsigned int current_point, unsigned int k) {
-
-  std::vector<unsigned int> erg(k);
-  if(k == 0 || m_vertices.size() < k) 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 = [&](const P &left, const P &right) { return (left.first < right.first);};
-
-  // using an ordered set for simplicitly, this may be slow
-  // because we frequently erase elements, causing the rb-tree
-  // to be rebuilt.
-  std::vector<P> candidates;
-  candidates.reserve(k + 1);
-
-  // 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
-    candidates.push_back({glm::distance(query, curr_point), begin + half_size - first});
-    std::push_heap(candidates.begin(), candidates.end(), cmp);
-    if(candidates.size() > k) {
-      std::pop_heap(candidates.begin(), candidates.end(), cmp);
-      candidates.pop_back();
-    }
-
-    // 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 = candidates[0].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});
-    }
-  }
-
-  // copy all nearest neighbour indices we found into
-  // a vector to return
-  for(unsigned int i = 0; i < k; i++) erg[i] = candidates[i].second;
-
-  return erg;
-}
-
-int CgPointCloud::getClosestPointToRay(glm::vec3 origin, glm::vec3 direction) {
-  if(m_vertices.size() == 0) return -1;
-
-  glm::vec3 inv_direction = glm::vec3(1.0 / direction.x, 1.0 / direction.y, 1.0 / direction.z);
-  glm::vec3* first = &m_vertices.front();
-  glm::vec3* last = &(*(m_vertices.end() - 1));
-  // sensible starting point
-  float erg_sqr_dist = std::numeric_limits<float>::infinity();
-  int erg = -1;
-
-  // 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 if it's closer to ray than last one
-    float cand_dist = sqrPointRayDist(curr_point, origin, direction);
-    if(cand_dist < erg_sqr_dist) {
-      erg_sqr_dist = cand_dist;
-      erg = begin + half_size - first;
-    }
-
-    // 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 children if their aabb is closer to ray than 
-    // the closest point we found until now and there are
-    // points in them
-    if (
-      begin != begin + half_size
-      && lower_aabb.sqrDistanceToRay(origin, direction, inv_direction) < erg_sqr_dist
-    ) {
-      nodes.push_back({begin, begin + half_size, lower_aabb, depth + 1});
-    }
-    if (
-      end != begin + half_size
-      && higher_aabb.sqrDistanceToRay(origin, direction, inv_direction) < erg_sqr_dist
-    ) {
-      nodes.push_back({begin + half_size + 1, end, higher_aabb, depth + 1});
-    }
+    // eye, center, up
+    // look from vertex position to the point the normal points at
+    glm::mat4 lookAt_matrix(glm::lookAt(
+      glm::vec3(m_vertices[i]),
+      glm::vec3(m_vertices[i] - m_vertex_normals[i]),
+      CgPointMath::getPerpendicularVector(m_vertex_normals[i])
+    ));
+    m_splat_orientations.push_back(lookAt_matrix);
+    m_splat_scaling.push_back(glm::vec2(0.02,0.02));
+
+    // use all points for splatting by default
+    m_splat_indices.push_back(i);
   }
 
-  return erg;
 }
 
 void CgPointCloud::setColors(std::vector<unsigned int> indices, glm::vec3 color) {
@@ -405,48 +77,4 @@ void CgPointCloud::showNormalColors() {
       (normal.z + 1.0) * .5
     );
   }
-}
-
-// calculates an arbitrary verctor perpendicular to the given one
-glm::vec3 CgPointCloud::getPerpendicularVector(glm::vec3 arg) {
-  if((arg.x == 0.0) && (arg.y == 0.0)) {
-    if(arg.z == 0.0) return glm::vec3(0.);
-    return glm::vec3(0.0, 1.0, 0.0);
-  }
-  return glm::normalize(glm::vec3(-arg.y, arg.x, 0.0));
-}
-
-const glm::vec3 CgPointCloud::getCenter() const {
-  return m_centroid;
-}
-
-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();
-
-  for(unsigned int i=0; i < m_vertices.size(); i++) {
-    glm::vec3 curr = m_vertices[i];
-    max_x = max_x < curr[0] ? curr[0] : max_x;
-    min_x = min_x > curr[0] ? curr[0] : min_x;
-    max_y = max_y < curr[1] ? curr[1] : max_y;
-    min_y = min_y > curr[1] ? curr[1] : min_y;
-    max_z = max_z < curr[2] ? curr[2] : max_z;
-    min_z = min_z > curr[2] ? curr[2] : min_z;
-  }
-
-  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 {
-  return m_splat_scaling;
 }
\ No newline at end of file
diff --git a/CgSceneGraph/CgPointCloud.h b/CgSceneGraph/CgPointCloud.h
index b530e4d546cd14ef594aafe62633454cc1a283c9..d7d9decf1e83c26d63e22ac637c2ae3099b72e74 100644
--- a/CgSceneGraph/CgPointCloud.h
+++ b/CgSceneGraph/CgPointCloud.h
@@ -13,84 +13,51 @@
 #include <tuple>
 #include <thread>
 #include "CgBase/CgBasePointCloud.h"
-#include "CgSceneGraph/CgAABB.h"
+#include "CgMath/CgAABB.h"
+#include "CgMath/CgPointWrangler.h"
+#include "CgMath/CgPointMath.h"
 #include "CgBase/CgEnums.h"
-#include "CgUtils/ObjLoader.h"
-#include "CgMath/CgDecomposition.h"
 #include <glm/glm.hpp>
 #include <glm/gtx/norm.hpp>
 #include <glm/gtc/matrix_transform.hpp>
-#include "CgUtils/Timer.h"
 
 class CgPointCloud : public CgBasePointCloud {
 public:
-  CgPointCloud();
-  CgPointCloud(int id);
-
+  CgPointCloud(
+    std::vector<glm::vec3> vertices,
+    std::vector<glm::vec3> vertex_normals,
+    std::vector<glm::vec3> vertex_colors,
+    std::vector<glm::mat4> splat_orientations,
+    std::vector<glm::vec2> splat_scaling,
+    std::vector<unsigned int> splat_indices
+  );
+  // copy ctor
+  CgPointCloud(const CgPointCloud& other);
   ~CgPointCloud();
 
   //inherited from CgBaseRenderableObject
   Cg::ObjectType getType() const;
   unsigned int getID() const;
-
-  //inherited from CgBasePointCloud
-
-  // vertex positions in local coordinates
   const std::vector<glm::vec3>& getVertices() const;
-
-  // normal directions in local coordinates (normalized)
   const std::vector<glm::vec3>& getVertexNormals() const;
-
-  // vertex colors
   const std::vector<glm::vec3>& getVertexColors() const;
-
-  // a local "look at matrix" for normal direction and the spanning directions of the tangent plane
   const std::vector<glm::mat4>& getSplatOrientations() const;
-
-  // for rendering ellipes use a different scaling in the two spanning directions of the tangent plane
   const std::vector<glm::vec2>& getSplatScalings() const;
-
-  // get the indices of the splats which should be rendered
-   const std::vector<unsigned int>& getSplatIndices() const;
-
-  // read a dataset from file, can cheat the normals, i.e read the mormals from file
-  void init( std::string filename, bool cheat_normals=false);
-
-  // nearest-neighbour search in the point cloud
-  std::vector<unsigned int> getNearestNeighborsSlow(unsigned int current_point, unsigned int k);
-  std::vector<unsigned int> getNearestNeighborsFast(unsigned int current_point, unsigned int k);
-
-  // get the point that's closest to the ray
-  int getClosestPointToRay(glm::vec3 origin, glm::vec3 direction);
+  const std::vector<unsigned int>& getSplatIndices() const;
+  const CgAABB getAABB() const;
+  const glm::vec3 getCentroid() const;
 
   // change colors of a list of points
   void setColors(std::vector<unsigned int> indices, glm::vec3 color);
   void showNormalColors();
   void resetColors(glm::vec3 color);
 
-  // the center of gravity of the object, for rendering
-  const glm::vec3 getCenter() const;
-
-  // returns the point clouds AABB
-  const CgAABB getAABB() const;
-
-  // traverse the kd-tree in BFO and call a function on each node
-  void traverseBFO(std::function<void(glm::vec3*, const unsigned int)> fn);
-
 private:
 
-    // the following demonstration methods have to be replaced by your own calculations
+    CgPointCloud(int id): m_type(Cg::PointCloud), m_id(id) {};
 
     //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);
-
-    // rearrange the vec3 between begin and end so they form a kd-tree
-    void buildKdTree(glm::vec3* begin, glm::vec3* end, int depth);
-    void estimateNormals();
 
     std::vector<glm::vec3> m_vertices;
     std::vector<glm::vec3> m_vertex_normals;
@@ -111,11 +78,13 @@ private:
 
 inline Cg::ObjectType  CgPointCloud::getType() const {return m_type;}
 inline unsigned int CgPointCloud::getID() const {return m_id;}
-inline const std::vector<glm::vec3>& CgPointCloud::getVertices() const{return m_vertices;}
-inline const std::vector<glm::vec3>& CgPointCloud::getVertexNormals() const{return m_vertex_normals;}
-inline const std::vector<glm::vec3>& CgPointCloud::getVertexColors() const{return m_vertex_colors;}
-inline const std::vector<glm::mat4>& CgPointCloud::getSplatOrientations() const{return m_splat_orientations;}
-inline const std::vector<unsigned int>& CgPointCloud::getSplatIndices() const{return m_splat_indices;}
-
+inline const std::vector<glm::vec3>& CgPointCloud::getVertices() const {return m_vertices;}
+inline const std::vector<glm::vec3>& CgPointCloud::getVertexNormals() const {return m_vertex_normals;}
+inline const std::vector<glm::vec3>& CgPointCloud::getVertexColors() const {return m_vertex_colors;}
+inline const std::vector<glm::mat4>& CgPointCloud::getSplatOrientations() const {return m_splat_orientations;}
+inline const std::vector<unsigned int>& CgPointCloud::getSplatIndices() const {return m_splat_indices;}
+inline const std::vector<glm::vec2>& CgPointCloud::getSplatScalings() const { return m_splat_scaling; }
+inline const CgAABB CgPointCloud::getAABB() const { return m_aabb; };
+inline const glm::vec3 CgPointCloud::getCentroid() const { return m_centroid; };
 
 #endif // CGPOINTCLOUD_H
\ No newline at end of file
diff --git a/CgSceneGraph/CgSceneControl.cpp b/CgSceneGraph/CgSceneControl.cpp
index 73d1d79d2a6dc35b8aee780026deb837cabc1baf..c92078abdb6203b1ceb243aa714852b227145805 100644
--- a/CgSceneGraph/CgSceneControl.cpp
+++ b/CgSceneGraph/CgSceneControl.cpp
@@ -13,8 +13,7 @@
 #include "CgUtils/ObjLoader.h"
 #include <string>
 
-CgSceneControl::CgSceneControl()
-{
+CgSceneControl::CgSceneControl() {
     m_pointcloud = nullptr;
     resetTransform();
 
@@ -49,9 +48,7 @@ void CgSceneControl::resetTransform() {
     m_current_transformation = glm::mat4(1.0);
 }
 
-
-CgSceneControl::~CgSceneControl()
-{
+CgSceneControl::~CgSceneControl() {
     if(m_pointcloud!=NULL)
         delete m_pointcloud;
     if(m_triangle_mesh!=NULL)
@@ -105,39 +102,6 @@ void CgSceneControl::calculatePickRay(float x, float y) {
   m_renderer->init(m_select_ray);
 }
 
-void CgSceneControl::calculateSingularValueDecomposition() {
-      //init some arbitrary test matrix
-      MatrixXd C(2,3);
-      C(0,0)=3.0;C(0,1)=2.0;C(0,2)=2.0;
-      C(1,0)=2.0;C(1,1)=3.0;C(1,2)=-2.0;
-
-      std::cout << C << std::endl;
-
-      //compute decomposition
-      JacobiSVD<Eigen::MatrixXd> svd(C,ComputeThinU|ComputeThinV);
-      Eigen::MatrixXd U=svd.matrixU();
-      Eigen::MatrixXd V=svd.matrixV();
-      Eigen::VectorXd SV=svd.singularValues();
-
-      //build a diagonal matrix out of the singular values
-      Eigen::MatrixXd S(2,2);
-      S.setZero();
-      S(0,0)=SV(0);S(1,1)=SV(1);
-
-      std::cout << U << std::endl;
-      std::cout << S << std::endl;
-      std::cout << V.transpose() << std::endl;
-
-
-    // compute Moore-Penrose inverse now
-
-     S(0,0)=1.0/S(0,0);S(1,1)=1.0/S(1,1);
-     Eigen::MatrixXd C_plus = V*S*U.transpose();
-
-     std::cout << std::endl;
-     std::cout << C_plus << std::endl;
-}
-
 // deletes the old kd-tree-mesh and creates a new one 
 // down to max_depth
 void CgSceneControl::createKdTreeViz(int max_depth) {
@@ -148,9 +112,12 @@ void CgSceneControl::createKdTreeViz(int max_depth) {
   // recalculate mesh
   if(m_pointcloud != nullptr && max_depth != -1) {
     std::vector<double> splits;
-    m_pointcloud->traverseBFO([&](glm::vec3* point, unsigned int depth) {
-        splits.push_back((*point)[depth % 3]);
-    });
+    std::vector<glm::vec3> vertices = m_pointcloud->getVertices();
+    CgPointWrangler::traverseBFO(
+      &(vertices.front()),
+      &(*(vertices.end() - 1)),
+      [&](glm::vec3* point, unsigned int depth) {splits.push_back((*point)[depth % 3]);}
+    );
     std::queue<CgAABB> aabb_queue;
     aabb_queue.push(m_pointcloud->getAABB());
     unsigned int index = 0;
@@ -191,9 +158,20 @@ void CgSceneControl::handlePickEvent(float x, float y) {
     std::vector<glm::vec3> points = m_select_ray->getVertices();
     glm::vec3 start = points[0];
     glm::vec3 direction = points[1] - start;
-    int picked_index = m_pointcloud->getClosestPointToRay(start, direction);
-    if(picked_index != -1) { // -1 => point found
-      std::vector<unsigned int> neighbors = m_pointcloud->getNearestNeighborsFast(picked_index, 50);
+    const glm::vec3* begin = &(m_pointcloud->getVertices().front());
+    const glm::vec3* end = &(*(m_pointcloud->getVertices().end() - 1));
+    int picked_index = CgPointWrangler::getClosestPointToRay(
+      begin, end, 
+      m_pointcloud->getAABB(), 
+      start, direction
+    );
+    if(picked_index != -1) { // -1 => no point found
+      std::vector<unsigned int> neighbors = CgPointWrangler::getNearestNeighborsFast(
+        picked_index, 
+        50, 
+        m_pointcloud->getVertices(), 
+        m_pointcloud->getAABB()
+      );
       m_pointcloud->resetColors(glm::vec3(0.0, 0.75, 0.0));
       m_pointcloud->setColors(neighbors, glm::vec3(0.75, 0.0, 0.0));
       m_pointcloud->setColors({picked_index}, glm::vec3(1.0, 1.0, 0.0));
@@ -203,8 +181,46 @@ void CgSceneControl::handlePickEvent(float x, float y) {
   }
 }
 
-void CgSceneControl::handleSimplifyEvent(CgSimplifyEvent ev) {
+void CgSceneControl::handleIncrementalSimplifyEvent(const CgSimplifyEvent& ev) {
   std::cout << ev << std::endl;
+  // get current points
+  std::vector<glm::vec3> points = m_pointcloud->getVertices();
+  std::vector<glm::vec3> normals = m_pointcloud->getVertexNormals();
+  std::vector<glm::vec3> colors = m_pointcloud->getVertexColors();
+  std::vector<glm::mat4> splat_orientations = m_pointcloud->getSplatOrientations();
+  std::vector<glm::vec2> splat_scaling = m_pointcloud->getSplatScalings();
+  std::vector<unsigned int> splat_indices = m_pointcloud->getSplatIndices();
+
+  delete m_pointcloud;
+  // simplify by filtering the points
+  std::vector<glm::vec3> new_points;
+  std::vector<glm::vec3> new_normals;
+  std::vector<glm::vec3> new_colors;
+  std::vector<glm::mat4> new_splat_orientations;
+  std::vector<glm::vec2> new_splat_scaling;
+  std::vector<unsigned int> new_splat_indices;
+  for(unsigned int i = 0; i < points.size(); i += 2) {
+    new_points.push_back(points[i]);
+    new_normals.push_back(normals[i]);
+    new_colors.push_back(colors[i]);
+    new_splat_orientations.push_back(splat_orientations[i]);
+    new_splat_scaling.push_back(splat_scaling[i]);
+    new_splat_indices.push_back(splat_indices[i]);
+  }
+
+  // re-init the new point cloud with the new point list
+  m_pointcloud = new CgPointCloud(
+    new_points,
+    new_normals,
+    new_colors,
+    new_splat_orientations,
+    new_splat_scaling,
+    new_splat_indices
+  );
+}
+
+void CgSceneControl::handleHierarchicalSimplifyEvent(const CgSimplifyEvent& ev) {
+  std::cout << "hierarchical" << std::endl;
 }
 
 void CgSceneControl::handleEvent(CgBaseEvent* e) {
@@ -263,16 +279,6 @@ void CgSceneControl::handleEvent(CgBaseEvent* e) {
     if(e->getType() & Cg::CgKeyEvent) {
         CgKeyEvent* ev = (CgKeyEvent*)e;
 
-        if(ev->text() == "e") {
-          // example usage of eigen library to compute the eigen-decomposition of a matrix
-          calculateEigenDecomposition3x3();
-        }
-
-        if(ev->text() == "s") {
-          // example usage of eigen library to compute the eigen-decomposition of a matrix
-          calculateSingularValueDecomposition();
-        }
-
         if(ev->text() == "c") {
           resetTransform();
           m_renderer->redraw();
@@ -326,7 +332,15 @@ void CgSceneControl::handleEvent(CgBaseEvent* e) {
 
     if(e->getType() & Cg::CgSimplifyEvent && m_pointcloud != nullptr) {
       CgSimplifyEvent* ev = (CgSimplifyEvent*) e;
-      handleSimplifyEvent(*ev);
+      if(ev->doReset()) {
+        CgPointCloud* copy = new CgPointCloud(*m_original_pointcloud);
+        delete m_pointcloud;
+        m_pointcloud = copy;
+      } else if(ev->getClusteringType() == Cg::CgIncremental) {
+        handleIncrementalSimplifyEvent(*ev);
+      } else {
+        handleHierarchicalSimplifyEvent(*ev);
+      }
       m_renderer->init(m_pointcloud);
       m_renderer->redraw();
     }
@@ -338,35 +352,54 @@ void CgSceneControl::handleEvent(CgBaseEvent* e) {
     }
 
     if(e->getType() & Cg::CgLoadPointCloudEvent) {
-
         CgLoadPointCloudEvent* ev = (CgLoadPointCloudEvent*)e;
-        ObjLoader* loader= new ObjLoader();
+        ObjLoader* loader = new ObjLoader();
 
-        std::string filename =ev->FileName();
-        if(filename !="")
-          {
-            loader->load(filename);
+        std::string filename = ev->FileName();
+        if(filename != "") {
+          loader->load(filename);
 
-            if(m_triangle_mesh!=NULL)
-              {
-                delete m_triangle_mesh;
-                m_triangle_mesh=NULL;
-              }
-            if(m_pointcloud!=NULL)
-              {
-                delete m_pointcloud;
-                m_pointcloud=NULL;
-              }
+          if(m_triangle_mesh != nullptr) {
+              delete m_triangle_mesh;
+              m_triangle_mesh = nullptr;
+          }
 
-            m_pointcloud= new CgPointCloud();
-            m_pointcloud->init(ev->FileName(),true);
+          if(m_pointcloud != nullptr) {
+              delete m_pointcloud;
+              m_pointcloud = nullptr;
+          }
 
-            m_center=m_pointcloud->getCenter();
-            m_renderer->init(m_pointcloud);
-            m_renderer->redraw();
+          std::vector<glm::vec3> positions;
+          std::vector<glm::vec3> normals;
+          std::vector<glm::vec3> colors;
+          loader->getPositionData(positions);
+
+          std::cout << "loaded " << positions.size() << " points" << std::endl;
+          std::cout << "using " << std::thread::hardware_concurrency() << " threads" << std::endl;
+          run_timed("build kd", 1, [&](){CgPointWrangler::buildKdTree(positions.data(), positions.data() + positions.size(), 0);});
+          run_timed("est. normals", 1, [&](){normals = CgPointWrangler::estimateNormals(15, positions);});
+          colors.resize(positions.size());
+          for(unsigned int i = 0; i < positions.size(); i++) {
+            colors[i] = (normals[i] + glm::vec3(1.0)) * .5;
           }
-    }
 
+          m_pointcloud = new CgPointCloud(
+            positions,
+            normals,
+            colors,
+            std::vector<glm::mat4>(positions.size()),
+            std::vector<glm::vec2>(positions.size()),
+            std::vector<unsigned int>(positions.size())
+          );
+          
+          // save away a copy
+          m_original_pointcloud = new CgPointCloud(*m_pointcloud);
+
+          m_center = m_pointcloud->getCentroid();
+          m_renderer->init(m_pointcloud);
+          m_renderer->redraw();
+      }
+    }
 
     if(e->getType() & Cg::CgLoadMeshEvent) {
 
diff --git a/CgSceneGraph/CgSceneControl.h b/CgSceneGraph/CgSceneControl.h
index f8d6bfd82838f24d524fa30627dc6e21e08cbbf3..20689b6dd265f5d672f6ace89db0d7cad02f9a4a 100644
--- a/CgSceneGraph/CgSceneControl.h
+++ b/CgSceneGraph/CgSceneControl.h
@@ -3,7 +3,7 @@
 
 #include "CgBase/CgObserver.h"
 #include "CgBase/CgBaseSceneControl.h"
-#include "CgAABB.h"
+#include "CgMath/CgAABB.h"
 #include <glm/glm.hpp>
 #include <vector>
 #include <queue>
@@ -21,6 +21,7 @@
 #include "CgEvents/CgPickRayEvent.h"
 #include "CgEvents/CgKdTreeEvent.h"
 #include "CgEvents/CgSimplifyEvent.h"
+#include "CgUtils/Timer.h"
 
 class CgBaseEvent;
 class CgBaseRenderer;
@@ -29,7 +30,6 @@ class CgTriangleMesh;
 class CgPolyLine;
 class CgQuad;
 
-
 class CgSceneControl : public CgObserver, public CgBaseSceneControl
 {
 public:
@@ -54,7 +54,8 @@ private:
     void calculateSingularValueDecomposition();
     void createKdTreeViz(int max_depth);
     void handlePickEvent(float x, float y);
-    void handleSimplifyEvent(CgSimplifyEvent ev);
+    void handleIncrementalSimplifyEvent(const CgSimplifyEvent& ev);
+    void handleHierarchicalSimplifyEvent(const CgSimplifyEvent& ev);
 
     // one single Disc used as Object to render the spats (in own local coordinates)
     CgTriangleMesh* m_disc;
@@ -62,6 +63,7 @@ private:
 
     // there is no scenegraph, either a point cloud or a mesh is used
     CgPointCloud* m_pointcloud;
+    CgPointCloud* m_original_pointcloud;
     CgTriangleMesh* m_triangle_mesh;
 
     // all the things that you would normally find in a scenegraph
diff --git a/CgUtils/ObjLoader.cpp b/CgUtils/ObjLoader.cpp
index 6d6b279f97d17fd67c259017fda152ac89fc6d40..00a4078a84900156ded62e67e48afee014a31d2f 100644
--- a/CgUtils/ObjLoader.cpp
+++ b/CgUtils/ObjLoader.cpp
@@ -324,7 +324,7 @@ void ObjLoader::getFaceIndexData(std::vector<unsigned int>& arg_indices)
 void ObjLoader::getPositionData(std::vector<glm::vec3>& arg_verts)
 {
     arg_verts.clear();
-    arg_verts=Positions;
+    arg_verts = Positions;
 }
 
 void ObjLoader::getNormalData(std::vector<glm::vec3>& arg_normals)
diff --git a/ExerciseVC.pro b/ExerciseVC.pro
index 0fdf709a4e0ca3a0e0c80bab96c7862a6073cf42..b2feb7d7818535433017aec333965234b3dfb25f 100644
--- a/ExerciseVC.pro
+++ b/ExerciseVC.pro
@@ -1,7 +1,7 @@
 QT       += core gui opengl widgets
 TEMPLATE = app
 TARGET = PointViewer
-QMAKE_CXXFLAGS += -std=c++20 -pthread
+QMAKE_CXXFLAGS += -g -std=c++20 -pthread
 CONFIG += c++20
 
 SOURCES += main.cpp \
@@ -22,12 +22,12 @@ SOURCES += main.cpp \
     CgSceneGraph/CgPolyLine.cpp \
     CgSceneGraph/CgQuad.cpp \
     CgSceneGraph/CgSceneControl.cpp \
+    CgSceneGraph/CgTriangleMesh.cpp \
+    CgMath/CgAABB.cpp \
     CgEvents/CgKeyEvent.cpp \
     CgQtViewer/CgQtGlBufferObject.cpp \
     CgQtViewer/CgTrackball.cpp \
     CgEvents/CgWindowResizeEvent.cpp \
-    CgSceneGraph/CgTriangleMesh.cpp \
-    CgSceneGraph/CgAABB.cpp \
     CgUtils/ObjLoader.cpp \
     CgEvents/CgTrackballEvent.cpp
 
@@ -39,10 +39,12 @@ HEADERS += \
     CgEvents/CgShowNormalsEvent.h \
     CgEvents/CgSimplifyEvent.h \
     CgEvents/CgSplatEvent.h \
-    CgMath/CgDecomposition.h \
+    CgMath/CgPointWrangler.h \
+    CgMath/CgPointMath.h \
     CgMath/Eigen/Core \
     CgMath/Eigen/Eigen \
     CgMath/Eigen/SVD \
+    CgMath/CgAABB.h \
     CgQtViewer/CgQtGLRenderWidget.h \
     CgQtViewer/CgQtGui.h \
     CgBase/CgObserver.h \
@@ -56,6 +58,7 @@ HEADERS += \
     CgSceneGraph/CgPolyLine.h \
     CgSceneGraph/CgQuad.h \
     CgSceneGraph/CgSceneControl.h \
+    CgSceneGraph/CgTriangleMesh.h \
     CgEvents/CgKeyEvent.h \
     CgBase/CgBaseRenderer.h \
     CgBase/CgBaseRenderableObject.h \
@@ -67,8 +70,6 @@ HEADERS += \
     CgQtViewer/CgQtGlBufferObject.h \
     CgQtViewer/CgTrackball.h \
     CgEvents/CgWindowResizeEvent.h \
-    CgSceneGraph/CgTriangleMesh.h \
-    CgSceneGraph/CgAABB.h \
     CgUtils/ObjLoader.h \
     CgUtils/Timer.h \
     CgBase/CgBaseImage.h \