diff --git a/CgEvents/CgShowNormalsEvent.cpp b/CgEvents/CgShowNormalsEvent.cpp
index 84e5a142e01c5c75fa19d68954c26f81a5a4c097..6009d0b40187bf699a47b8782bef452bec1910a9 100644
--- a/CgEvents/CgShowNormalsEvent.cpp
+++ b/CgEvents/CgShowNormalsEvent.cpp
@@ -11,9 +11,9 @@ Cg::EventType CgShowNormalsEvent::getType()
 }
 
 
-std::ostream& operator<<(std::ostream& os,const CgShowNormalsEvent& e)
+std::ostream& operator<<(std::ostream& os, const CgShowNormalsEvent& e)
 {
-    os << "CgShowNormalsEvent" << " ";
+    os << "CgShowNormalsEvent " << &e <<" ";
     return os;
 }
 
diff --git a/CgEvents/CgTrackballEvent.cpp b/CgEvents/CgTrackballEvent.cpp
index ad09bc2d2f16793c8a2f2943298fc5f38b2fe948..6ae1fbd6522dda9fcda4e2e583260ea37046666b 100644
--- a/CgEvents/CgTrackballEvent.cpp
+++ b/CgEvents/CgTrackballEvent.cpp
@@ -19,7 +19,7 @@ Cg::EventType  CgTrackballEvent::getType()
 
 std::ostream& operator<<(std::ostream& os,const CgTrackballEvent& e)
 {
-    os << "CgTrackballEvent: "<< std::endl;
+    os << "CgTrackballEvent: "<< &e << std::endl;
     return os;
 }
 
diff --git a/CgMath/CgPointMath.h b/CgMath/CgPointMath.h
index edcebed1811baf142509b8f4b9f365a1a67834f1..dfc8e1b67c4d51799a871041f0d753a5b4ef2ef5 100644
--- a/CgMath/CgPointMath.h
+++ b/CgMath/CgPointMath.h
@@ -15,13 +15,14 @@ using namespace Eigen;
 #include "CgAABB.h"
 
 struct PCA {
+    double radius;
     glm::vec3 centroid;
     glm::vec3 evec0;
     glm::vec3 evec1;
     glm::vec3 evec2;
-    float eval0;
-    float eval1;
-    float eval2;
+    double eval0;
+    double eval1;
+    double eval2;
 };
 
 /*
@@ -35,7 +36,9 @@ class CgPointMath {
         static CgAABB calculateAABB(const std::vector<glm::vec3>& vertices);
         static glm::mat3 getCovarianceMatrix(const std::vector<glm::vec3> &points);
         static PCA calculatePCA(const std::vector<glm::vec3> &points);
+        static double varianceFromPCA(PCA pca);
         static glm::vec3 getPerpendicularVector(glm::vec3 arg);
+        static glm::vec3 normalToColor(glm::vec3 normal);
     private:
         CgPointMath(){}
 };
@@ -103,7 +106,7 @@ inline float CgPointMath::sqrBoxRayDist(CgAABB box, glm::vec3 origin, glm::vec3
 
 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];
+    for(auto p : points) center += p;
     center /= double(points.size());
     return center;
 }
@@ -156,10 +159,12 @@ inline PCA CgPointMath::calculatePCA(const std::vector<glm::vec3> &points) {
     unsigned int size = points.size();
     Eigen::MatrixXd mat(size, 3);
     glm::vec3 centroid = calculateCentroid(points);
+    double radius = 0.0;
     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];
+        radius = std::max(radius, double(glm::distance(centroid, points[i])));
     }
     Eigen::Matrix3d res = (1 / double(size)) * (mat.transpose() * mat);
 
@@ -174,6 +179,7 @@ inline PCA CgPointMath::calculatePCA(const std::vector<glm::vec3> &points) {
     Eigen::Vector3d evec1 = es.eigenvectors().col(1);
     Eigen::Vector3d evec2 = es.eigenvectors().col(2);
     return {
+        radius,
         centroid,
         glm::vec3(evec0(0), evec0(1), evec0(2)),
         glm::vec3(evec1(0), evec1(1), evec1(2)),
@@ -184,6 +190,10 @@ inline PCA CgPointMath::calculatePCA(const std::vector<glm::vec3> &points) {
     };
 }
 
+inline double CgPointMath::varianceFromPCA(PCA pca) {
+    return pca.eval0 / (pca.eval0 + pca.eval1 + pca.eval2);
+}
+
 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.);
@@ -192,4 +202,12 @@ inline glm::vec3 CgPointMath::getPerpendicularVector(glm::vec3 arg) {
   return glm::normalize(glm::vec3(-arg.y, arg.x, 0.0));
 }
 
+inline glm::vec3 CgPointMath::normalToColor(glm::vec3 normal) {
+    return glm::vec3(
+        (normal.x + 1.0) * 0.5, 
+        (normal.y + 1.0) * 0.5,
+        (normal.z + 1.0) * 0.5
+    );
+}
+
 #endif // CGPOINTMATH_H
\ No newline at end of file
diff --git a/CgMath/CgPointWrangler.h b/CgMath/CgPointWrangler.h
index 88c631587149f6a4d18a130765901847d1cab59d..4e5460e094138760c2af0534ad4f46428f0e1bd7 100644
--- a/CgMath/CgPointWrangler.h
+++ b/CgMath/CgPointWrangler.h
@@ -31,10 +31,8 @@ class CgPointWrangler {
     public:
     
         // rearrange the vec3 between begin and end so they form a kd-tree
-        static void buildKdTree(
-            glm::vec3* begin, glm::vec3* end, 
-            int depth
-        );
+        static void buildKdTree(glm::vec3* begin, glm::vec3* end, unsigned int depth);
+        static void buildKdTree(PCA* begin, PCA* end, unsigned int depth);
 
         // traverse a kd-tree between begin and end in BFO and call a function on each node
         static void traverseBFO(
@@ -51,8 +49,28 @@ class CgPointWrangler {
 
         // naive estimation of normals
         static std::vector<std::vector<unsigned int>> getKNeighbourhoods(unsigned int neighbourhood_size, const std::vector<glm::vec3>& vertices);
-        static std::vector<glm::vec3> estimateNormals(const std::vector<glm::vec3>& vertices, const std::vector<std::vector<unsigned int>>& neighbourhoods);
-        static void alignNormals();
+        static std::vector<std::vector<glm::vec3>> kNeighbourhoodsToClusters(const std::vector<glm::vec3>& vertices, const std::vector<std::vector<unsigned int>>& neighbourhoods);
+        static std::vector<std::pair<unsigned int, unsigned int>> makeNaiveSpanningTree(const std::vector<std::vector<unsigned int>>& neighbourhoods, const std::vector<PCA>& pca);
+        static std::vector<PCA> performPcaOnClusters(std::vector<std::vector<glm::vec3>> clusters); 
+        static std::vector<glm::vec3> alignNormals(
+            glm::vec3 ref_normal, // reference the first normal will be aligned to 
+            const std::vector<std::pair<unsigned int, unsigned int>> spanning_tree, // list of edges (from, to)
+            const std::vector<PCA>& pca // contains the unaligned normals
+        );
+
+        static std::vector<PCA> performHierarchicalSimplification(
+            std::vector<glm::vec3> original_points,
+            float max_cluster_size,
+            float max_cluster_variance,
+            float max_cluster_radius
+        );
+
+        static std::vector<PCA> performIncrementalSimplification(
+            std::vector<glm::vec3> original_points,
+            float max_cluster_size,
+            float max_cluster_variance,
+            float max_cluster_radius
+        );
 
         static std::vector<unsigned int> getNearestNeighborsSlow(
             unsigned int current_point, 
@@ -67,21 +85,30 @@ class CgPointWrangler {
             CgAABB aabb
         );
 
+        static std::vector<unsigned int> getNearestNeighborsFast(
+            glm::vec3 query, 
+            unsigned int k, 
+            const std::vector<glm::vec3>& vertices,
+            CgAABB aabb
+        );
+
     private:
         // 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);
+        template<typename T>
+        static unsigned int halfSize(const T* begin, const T* end);
         
         CgPointWrangler(){}
 };
 
-inline unsigned int CgPointWrangler::halfSize(const glm::vec3* begin, const glm::vec3* end) {
+template<typename T>
+inline unsigned int CgPointWrangler::halfSize(const T* begin, const T* 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) {
+inline void CgPointWrangler::buildKdTree(glm::vec3* begin, glm::vec3* end, unsigned int depth) {
     unsigned int half_size = halfSize(begin, end);
 
     // none or one element
@@ -90,7 +117,7 @@ inline void CgPointWrangler::buildKdTree(glm::vec3* begin, glm::vec3* end, int d
     // 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) {
+    std::nth_element(begin, begin + half_size, end, [=](const glm::vec3& a, const glm::vec3& b){
         return a[depth % 3] < b[depth % 3];
     });
 
@@ -99,6 +126,17 @@ inline void CgPointWrangler::buildKdTree(glm::vec3* begin, glm::vec3* end, int d
     buildKdTree(begin + half_size + 1, end, depth + 1);
 }
 
+// see the version for glm::vec3
+inline void CgPointWrangler::buildKdTree(PCA* begin, PCA* end, unsigned int depth) {
+    unsigned int half_size = halfSize(begin, end);
+    if(half_size == 0) return;
+    std::nth_element(begin, begin + half_size, end, [=](const PCA& a, const PCA& b){
+        return a.centroid[depth % 3] < b.centroid[depth % 3];
+    });
+    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, 
@@ -202,28 +240,72 @@ inline std::vector<std::vector<unsigned int>> CgPointWrangler::getKNeighbourhood
         const unsigned int count
     ) {
         for(unsigned int i = start_index; i < start_index + count; i++) {
-            neighbourhoods[i] = getNearestNeighborsFast(i, neighbourhood_size, vertices, aabb);
+            neighbourhoods[i] = getNearestNeighborsFast(i, neighbourhood_size, vertices, aabb);    
         }
     });
 
     return neighbourhoods;
 }
 
-inline std::vector<glm::vec3> CgPointWrangler::estimateNormals(
+inline std::vector<std::vector<glm::vec3>> CgPointWrangler::kNeighbourhoodsToClusters(
     const std::vector<glm::vec3>& vertices,
     const std::vector<std::vector<unsigned int>>& neighbourhoods
 ) {
-    std::vector<glm::vec3> normals;
-    normals.resize(vertices.size());
+    std::vector<std::vector<glm::vec3>> clusters;
+    clusters.resize(neighbourhoods.size());
+
+    run_threaded(neighbourhoods.size(), [&](
+        const unsigned int start_index, 
+        const unsigned int count
+    ) {
+        for(unsigned int i = start_index; i < start_index + count; i++) {
+            std::vector<unsigned int> nh_indices = neighbourhoods[i];
+            clusters[i].resize(nh_indices.size());
+            for(unsigned int j = 0; j < nh_indices.size(); j++) {
+                unsigned int index = nh_indices[j];
+                glm::vec3 p = vertices[index];
+                clusters[i][j] = p;
+            };
+        }
+    });
+
+    return clusters;
+}
+
+inline std::vector<PCA> CgPointWrangler::performPcaOnClusters(std::vector<std::vector<glm::vec3>> clusters) {
+    std::vector<PCA> pca;
+    pca.resize(clusters.size());
+
+    // worker for doing the PCAs
+    run_threaded(pca.size(), [&](
+        const unsigned int start_index, 
+        const unsigned int count
+    ) {
+        for(unsigned int i = start_index; i < start_index + count; i++) {
+            std::vector<glm::vec3> nh = clusters[i];
+            pca[i] = CgPointMath::calculatePCA(nh);
+        }
+    });
+
+    return pca;
+}
+
+inline std::vector<std::pair<unsigned int, unsigned int>> CgPointWrangler::makeNaiveSpanningTree(
+    const std::vector<std::vector<unsigned int>>& neighbourhoods,
+    const std::vector<PCA>& pca
+) {
+    std::cout << pca.size() << " pca elements for spanning tree" << std::endl;
+    std::vector<std::pair<unsigned int, unsigned int>> spanning_tree;
+    spanning_tree.clear();
+    spanning_tree.reserve(pca.size() + 1);
 
     // find point with highest y-value
     std::vector<bool> visited;
-    visited.resize(vertices.size());
+    visited.resize(pca.size(), false);
     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;
+    for(unsigned int i = 0; i < pca.size(); i++) {
+        float curr = pca[i].centroid.y;
         if(max_y < curr) {
             max_y = curr;
             max_y_index = i;
@@ -231,84 +313,225 @@ inline std::vector<glm::vec3> CgPointWrangler::estimateNormals(
     }
     visited[max_y_index] = true;
 
-    // used as a stack for traversal order
-    std::vector<unsigned int> spanning_tree;
+    std::vector<unsigned int> stack;
     // the "best" parent found for each vertex
     std::vector<std::pair<unsigned int, double>> scores;
-    scores.resize(vertices.size());
-
-    // worker for doing the PCAs
-    run_threaded(vertices.size(), [&](
-        const unsigned int start_index, 
-        const unsigned int count
-    ) {
-        // thread-local buffer for points of one neighbourhood
-        std::vector<glm::vec3> nh;
-        for(unsigned int i = start_index; i < start_index + count; i++) {
-            std::vector<unsigned int> nh_indices = neighbourhoods[i];
-            nh.clear();
-            nh.resize(nh_indices.size());
-            for(int j = 0; j < nh_indices.size(); j++) nh[j] = vertices[nh_indices[j]];
-            normals[i] = CgPointMath::calculatePCA(nh).evec0;
-        }
-    });
+    scores.resize(pca.size());
 
-    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;
+    auto visit = [&](unsigned int curr) {
+        glm::vec3 normal = pca[curr].evec0;
 
         // 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 < curr_nh_ind.size(); i++) {
+        for(unsigned int i = curr_nh_ind.size() - 1; i > 0; i--) {
             unsigned int index = curr_nh_ind[i];
-            // score by how far away it is from the local tangent plane
+            // score by how far away it is from the current point
             // 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 alignment = 1.0 - std::abs(glm::dot(pca[index].evec0, normal));
+            double distance = glm::distance(pca[index].centroid, pca[curr].centroid);
             double score = distance * alignment;
             if(!visited[index]) {
                 scores[index].first = curr;
                 scores[index].second = score;
-                spanning_tree.push_back(index);
+                stack.push_back(index);
                 visited[index] = true;
-            } 
-            if(scores[index].second > score) {
-                // update current neighbors parent to current node
+            } else if(scores[index].second > score) {
+                // update current neighbors parent to current node, since it's closer
                 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]);
+    // get started by aligning the normal of the highest point
+    spanning_tree.push_back({max_y_index, max_y_index});
+    visit(max_y_index); 
 
     // start processing
-    while(!spanning_tree.empty()) {
-        unsigned int curr = spanning_tree.back();
+    while(!stack.empty()) {
+        unsigned int curr = stack.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);
+        stack.pop_back();
+        spanning_tree.push_back({parent, curr});
+        visit(curr);
     }
 
-    return normals;
+    return spanning_tree;
 }
 
-inline void CgPointWrangler::alignNormals(
-   // glm::vec3 ref_normal, 
-   // const std::vector<unsigned int>& edges, 
-   // std::vector<glm::vec3>& normals
+inline std::vector<glm::vec3> CgPointWrangler::alignNormals(
+    glm::vec3 ref_normal, // reference the first normal will be aligned to
+    const std::vector<std::pair<unsigned int, unsigned int>> spanning_tree, // list of edges (from, to)
+    const std::vector<PCA>& pca // contains the unaligned normals
 ) {
+    std::cout << spanning_tree.size() << " spanning tree edges" << std::endl;
+    std::vector<glm::vec3> normals;
+    normals.resize(pca.size());
+    for(unsigned int i = 0; i < pca.size(); i++) normals[i] = pca[i].evec0;
+    // align first normal with ref
+    unsigned int start_index = spanning_tree[0].first;
+    normals[start_index] = glm::dot(normals[start_index], ref_normal) > 0.0 
+        ? normals[start_index] 
+        : -(normals[start_index]);
+
+    // align the rest 
+    // (spanning tree should ensure that we don't align to a normal that's not yet aligned itself)
+    for(unsigned int i = 1; i < spanning_tree.size(); i++) {
+        unsigned int parent_index = spanning_tree[i].first;
+        unsigned int own_index = spanning_tree[i].second;
+        normals[own_index] = glm::dot(normals[own_index], normals[parent_index]) > 0.0 
+            ? normals[own_index] 
+            : -(normals[own_index]);
+    }
 
+    return normals;
+}
+
+inline std::vector<PCA> CgPointWrangler::performHierarchicalSimplification(
+    std::vector<glm::vec3> original_points,
+    float max_cluster_size,
+    float max_cluster_variance,
+    float max_cluster_radius
+){
+    unsigned int MAX_POINTS = std::min(100.0f, max_cluster_size);
+    unsigned int MIN_POINTS = 6;
+
+    std::vector<unsigned int> hist;
+    hist.resize(MAX_POINTS + 1, 0);
+
+    // used by the threads to define the clusters
+    std::vector<unsigned int> indices;
+    indices.resize(original_points.size());
+    for(unsigned int i = 0; i < original_points.size(); i++) indices[i] = i;
+
+    // used to define the work portions
+    std::vector<std::pair<unsigned int, unsigned int>> jobs; // list of ranges of indices
+    jobs.push_back({0, indices.size()}); // initially, use entire point cloud
+    std::vector<PCA> pca; // finished clusters
+
+    while(!jobs.empty()) {
+        std::pair<unsigned int, unsigned int> job = jobs.back();
+        jobs.pop_back();
+
+        // get cluster and check if it fits constraints
+        std::vector<glm::vec3> cluster;
+        for(unsigned int i = job.first; i < job.second; i++) cluster.push_back(original_points[indices[i]]);
+        PCA current_pca = CgPointMath::calculatePCA(cluster);
+        double variance = CgPointMath::varianceFromPCA(current_pca);
+        if(cluster.size() <= MIN_POINTS  || (variance < max_cluster_variance && current_pca.radius < max_cluster_radius && cluster.size() <= MAX_POINTS)) {
+        // push cluster as-is into cluster vector
+        hist[cluster.size()] += 1;
+        pca.push_back(current_pca);
+        } else {
+        // split it & recurse
+        // indices vector is not behind mutex because 
+        // the threads work on mutually exclusive ranges
+        unsigned int lower = job.first;
+        unsigned int upper = job.second;
+        while(lower != upper) {
+            // separate indices depending on which side of the centroid-normal plane
+            // their points are
+            if(glm::dot(original_points[indices[lower]] - current_pca.centroid, current_pca.evec2) > 0) {
+            lower += 1;
+            } else {
+            upper -= 1;
+            unsigned int t = indices[upper];
+            indices[upper] = indices[lower];
+            indices[lower] = t;
+            }
+        }
+        
+        jobs.push_back({job.first, lower});
+        jobs.push_back({lower, job.second});
+        }
+    }
+
+    // clustering is done 
+    std::cout << pca.size() << " clusters generated, histogram (cluster_size:count):" << std::endl;
+    for(unsigned int i = 0; i < hist.size(); i++) std::cout << " |" << i + 1 << ":" << hist[i];
+    std::cout << std::endl;
+    return pca;
+}
+
+inline std::vector<PCA> CgPointWrangler::performIncrementalSimplification(
+    std::vector<glm::vec3> original_points,
+    float max_cluster_size,
+    float max_cluster_variance,
+    float max_cluster_radius
+){
+    unsigned int MAX_POINTS = std::min(100.0f, max_cluster_size);
+    std::vector<glm::vec3> remaining_points;
+    std::vector<unsigned int> current_nh;
+    std::vector<bool> used;
+    used.resize(original_points.size(), false);
+    std::vector<PCA> pca;
+    glm::vec3 next_seed_point = original_points[0];
+    std::vector<unsigned int> hist;
+    hist.resize(MAX_POINTS + 1, 0);
+    while(original_points.size() > 0) { // while there are still points to consider...
+        unsigned int k;
+        CgAABB aabb = CgPointMath::calculateAABB(original_points);
+        PCA current_pca;
+        current_pca.centroid = next_seed_point;
+        // walk up to the largest cluster that fits the constraints
+        // for each size:
+        for(k = 1; k <= MAX_POINTS + 1; k++) {
+        // get neighbourhood indices
+        std::vector<unsigned int> nh = CgPointWrangler::getNearestNeighborsFast(
+            current_pca.centroid, k,
+            original_points, aabb
+        );
+
+        // get neighbourhood points from indices
+        std::vector<glm::vec3> cluster(nh.size());
+        for(unsigned int i = 0; i < nh.size(); i++) cluster[i] = original_points[nh[i]];
+
+        // do pca and check if cluster fits constraint
+        PCA cluster_pca = CgPointMath::calculatePCA(cluster);
+        float variance = CgPointMath::varianceFromPCA(cluster_pca);
+        if(variance > max_cluster_variance || cluster_pca.radius > max_cluster_radius) {
+            // the point that's last in the cluster is furthest away from the seed
+            // use it as next seed
+            next_seed_point = cluster[cluster.size()- 1];
+            hist[nh.size()] += 1;
+            break;
+        }
+        current_nh = nh;
+        current_pca = cluster_pca;
+        // if there are not enough points left to break the constraint
+        if(k > nh.size()) break;
+        }
+
+        // push new point information
+        pca.push_back(current_pca);
+
+        // if we have all the rest of the points in the current cluster, we're done
+        if(current_nh.size() == original_points.size()) break;
+
+        // copy unused points over
+        used.clear();
+        used.resize(original_points.size(), false);
+        remaining_points.clear();
+        remaining_points.reserve(original_points.size());
+        for(unsigned int i = 0; i < current_nh.size(); i++) used[current_nh[i]] = true;
+        for(unsigned int i = 0; i < original_points.size(); i++) {
+        if(!used[i]) remaining_points.push_back(original_points[i]);
+        used[i] = false;
+        }
+        // build new kd-tree in remaining points.
+        original_points = remaining_points;
+        CgPointWrangler::buildKdTree(original_points.data(), original_points.data() + original_points.size(), 0);
+    }
+
+    // clustering is done 
+    std::cout << pca.size() << " clusters generated, histogram (cluster_size:count):" << std::endl;
+    for(unsigned int i = 0; i < hist.size(); i++) std::cout << " |" << i + 1 << ":" << hist[i];
+    std::cout << std::endl;
+
+    return pca;
 }
 
 inline std::vector<unsigned int> CgPointWrangler::getNearestNeighborsSlow(
@@ -324,14 +547,10 @@ inline std::vector<unsigned int> CgPointWrangler::getNearestNeighborsSlow(
         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;
 
-        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);
-    }
+    for(unsigned int i = 0; i < k; i++) erg.push_back(distances[i].second);
 
     return erg;
 }
@@ -341,15 +560,28 @@ inline std::vector<unsigned int> CgPointWrangler::getNearestNeighborsFast(
     unsigned int k,
     const std::vector<glm::vec3>& vertices,
     CgAABB aabb
+){
+    glm::vec3 query = vertices[current_point];
+    return getNearestNeighborsFast(
+        query, k, vertices, aabb
+    );
+}
+
+inline std::vector<unsigned int> CgPointWrangler::getNearestNeighborsFast(
+    glm::vec3 query, 
+    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;
+    if(k == 0) 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();
+    // comp function for max heap
     auto cmp = [&](const Pair &left, const Pair &right) { return (left.first < right.first);};
 
     // using an ordered set for simplicitly, this may be slow
@@ -377,6 +609,7 @@ inline std::vector<unsigned int> CgPointWrangler::getNearestNeighborsFast(
         candidates.push_back({glm::distance(query, curr_point), begin + half_size - first});
         std::push_heap(candidates.begin(), candidates.end(), cmp);
         if(candidates.size() > k) {
+            // pop the candidate with the largest distance
             std::pop_heap(candidates.begin(), candidates.end(), cmp);
             candidates.pop_back();
         }
@@ -387,7 +620,8 @@ inline std::vector<unsigned int> CgPointWrangler::getNearestNeighborsFast(
         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
+        // neighbourhood bounding sphere 
+        // (candidates is a max-heap with the first element being furthest away from query)
         double neighbourhood_radius = candidates[0].first;
         if(
             begin < begin + half_size
@@ -396,7 +630,7 @@ inline std::vector<unsigned int> CgPointWrangler::getNearestNeighborsFast(
             nodes.push_back({begin, begin + half_size, lower_aabb, depth + 1});
         }
         if(
-            begin + half_size < end
+            begin + half_size + 1 < 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});
@@ -405,8 +639,7 @@ inline std::vector<unsigned int> CgPointWrangler::getNearestNeighborsFast(
 
     // 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;
-
+    for(unsigned int i = 0; i < candidates.size(); i++)  erg[i] = candidates[i].second;
     return erg;
 }
 
diff --git a/CgQtViewer/CgQtGui.cpp b/CgQtViewer/CgQtGui.cpp
index 6fe4600521f5c52322fd3b8d78879a72cf942886..93e54519b0f43800b2204213df6101f413226d8b 100644
--- a/CgQtViewer/CgQtGui.cpp
+++ b/CgQtViewer/CgQtGui.cpp
@@ -193,7 +193,7 @@ void CgQtGui::createSimplifyOptionPanel(QWidget* parent) {
     tab2_control->addWidget(max_number_of_points_spinbox);
     max_number_of_points_spinbox->setMinimum(1);
     max_number_of_points_spinbox->setMaximum(100);
-    max_number_of_points_spinbox->setValue(15);
+    max_number_of_points_spinbox->setValue(30);
     tab2_control->addWidget(max_number_of_points_spinbox);
     max_number_of_points_spinbox->setAlignment(Qt::AlignTop);
     max_number_of_points_spinbox->setSizePolicy(QSizePolicy::Ignored, QSizePolicy::Fixed);
@@ -206,10 +206,11 @@ void CgQtGui::createSimplifyOptionPanel(QWidget* parent) {
 
     max_diameter_of_cluster_spinbox = new QDoubleSpinBox();
     tab2_control->addWidget(max_diameter_of_cluster_spinbox);
-    max_diameter_of_cluster_spinbox->setMinimum(0.01);
+    max_diameter_of_cluster_spinbox->setDecimals(3);
+    max_diameter_of_cluster_spinbox->setMinimum(0.001);
     max_diameter_of_cluster_spinbox->setMaximum(10.0);
-    max_diameter_of_cluster_spinbox->setValue(0.05);
-    max_diameter_of_cluster_spinbox->setSingleStep(0.05);
+    max_diameter_of_cluster_spinbox->setValue(0.1);
+    max_diameter_of_cluster_spinbox->setSingleStep(0.001);
     max_diameter_of_cluster_spinbox->setAlignment(Qt::AlignTop);
     max_diameter_of_cluster_spinbox->setSizePolicy(QSizePolicy::Ignored, QSizePolicy::Fixed);
 
@@ -221,26 +222,27 @@ void CgQtGui::createSimplifyOptionPanel(QWidget* parent) {
 
     max_variance_of_cluster_spinbox = new QDoubleSpinBox();
     tab2_control->addWidget(max_variance_of_cluster_spinbox);
-    max_variance_of_cluster_spinbox->setMinimum(0.01);
+    max_variance_of_cluster_spinbox->setDecimals(3);
+    max_variance_of_cluster_spinbox->setMinimum(0.001);
     max_variance_of_cluster_spinbox->setMaximum(5.0);
-    max_variance_of_cluster_spinbox->setValue(0.05);
-    max_variance_of_cluster_spinbox->setSingleStep(0.05);
+    max_variance_of_cluster_spinbox->setValue(0.001);
+    max_variance_of_cluster_spinbox->setSingleStep(0.001);
     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);
+    number_of_points_checkbox->setChecked(true);
     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);
+    diameter_of_cluster_checkbox->setChecked(true);
     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);
+    variance_of_cluster_checkbox->setChecked(true);
     tab2_control->addWidget(variance_of_cluster_checkbox);
 
     /* select clustering type button group */
@@ -250,10 +252,10 @@ void CgQtGui::createSimplifyOptionPanel(QWidget* parent) {
     clustering_type_buttongroup = new QButtonGroup(subBox);
     clustering_type_buttongroup->setExclusive(true);
 
-    QRadioButton* radiobutton1 = new QRadioButton( "Incremental");
+    QRadioButton* radiobutton1 = new QRadioButton( "Incremental (slow!)");
     QRadioButton* radiobutton2 = new QRadioButton( "Hierarchical");
 
-    radiobutton1->setChecked(true);
+    radiobutton2->setChecked(true);
 
     clustering_type_buttongroup->addButton(radiobutton1, int(Cg::CgIncremental));
     clustering_type_buttongroup->addButton(radiobutton2, int(Cg::CgHierarchical));
diff --git a/CgSceneGraph/CgPointCloud.cpp b/CgSceneGraph/CgPointCloud.cpp
index 332f93626d2e97b5eab1751c27e2b5930c559742..d3dc6dee791142a5e230b1e1ab6d5369c3ee7853 100644
--- a/CgSceneGraph/CgPointCloud.cpp
+++ b/CgSceneGraph/CgPointCloud.cpp
@@ -5,15 +5,13 @@ CgPointCloud::CgPointCloud(
   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
+  std::vector<glm::vec2> splat_scaling
 ): 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);
@@ -24,15 +22,13 @@ CgPointCloud::CgPointCloud(const CgPointCloud& other) : CgPointCloud::CgPointClo
   other.m_vertex_normals,
   other.m_vertex_colors,
   other.m_splat_orientations,
-  other.m_splat_scaling,
-  other.m_splat_indices
+  other.m_splat_scaling
 ){}
 
 CgPointCloud::~CgPointCloud() {
     m_vertices.clear();
     m_vertex_normals.clear();
     m_vertex_colors.clear();
-    m_splat_indices.clear();
 }
 
 void CgPointCloud::setColors(std::vector<unsigned int> indices, glm::vec3 color) {
@@ -44,7 +40,7 @@ void CgPointCloud::resetColors(glm::vec3 color) {
 }
 
 void CgPointCloud::showNormalColors() {
-  for(int i = 0; i < m_vertices.size(); i++) {
+  for(unsigned int i = 0; i < m_vertices.size(); i++) {
     glm::vec3 normal = m_vertex_normals[i];
     m_vertex_colors[i] = glm::vec3(
       (normal.x + 1.0) * .5, 
diff --git a/CgSceneGraph/CgPointCloud.h b/CgSceneGraph/CgPointCloud.h
index 2d7cb9f92a24345ae9fcc4e3c304e2ed85a25fdb..c47c2ec57b324ac64ec29b61c95dcc0dd9c8c76a 100644
--- a/CgSceneGraph/CgPointCloud.h
+++ b/CgSceneGraph/CgPointCloud.h
@@ -28,8 +28,7 @@ public:
     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
+    std::vector<glm::vec2> splat_scaling
   );
   // copy ctor
   CgPointCloud(const CgPointCloud& other);
@@ -43,7 +42,6 @@ public:
   const std::vector<glm::vec3>& getVertexColors() const;
   const std::vector<glm::mat4>& getSplatOrientations() const;
   const std::vector<glm::vec2>& getSplatScalings() const;
-  const std::vector<unsigned int>& getSplatIndices() const;
   const CgAABB getAABB() const;
   const glm::vec3 getCentroid() const;
 
@@ -75,7 +73,6 @@ inline const std::vector<glm::vec3>& CgPointCloud::getVertices() const {return m
 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; };
diff --git a/CgSceneGraph/CgSceneControl.cpp b/CgSceneGraph/CgSceneControl.cpp
index 3373dc9db748c91bc47dd6233d48b830103ea5de..8717eebc7fb81bea2689f6dfa2bdd8e6c8187a1d 100644
--- a/CgSceneGraph/CgSceneControl.cpp
+++ b/CgSceneGraph/CgSceneControl.cpp
@@ -57,7 +57,7 @@ CgSceneControl::~CgSceneControl() {
         delete m_disc;
     if(m_select_ray!=NULL)
         delete m_select_ray;
-    for(int i = 0; i < m_kd_tree_mesh.size(); i++) delete m_kd_tree_mesh[i];
+    for(unsigned int i = 0; i < m_kd_tree_mesh.size(); i++) delete m_kd_tree_mesh[i];
     m_kd_tree_mesh.clear();
 }
 
@@ -106,158 +106,78 @@ void CgSceneControl::calculatePickRay(float x, float y) {
 // down to max_depth
 void CgSceneControl::createKdTreeViz(int max_depth) {
   // clear old quads
-  for(int i = 0; i < m_kd_tree_mesh.size(); i++) delete m_kd_tree_mesh[i];
+  for(unsigned int i = 0; i < m_kd_tree_mesh.size(); i++) delete m_kd_tree_mesh[i];
   m_kd_tree_mesh.clear();
-
+  if(m_pointcloud == nullptr || max_depth == -1) return;
+  unsigned int max_depth_u = max_depth;
   // recalculate mesh
-  if(m_pointcloud != nullptr && max_depth != -1) {
-    std::vector<double> splits;
-    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;
-    unsigned int depth = 0;
-    while(index < splits.size() && depth < max_depth) {
-      CgAABB current_aabb = aabb_queue.front();
-      double current_split = splits[index];
-      unsigned int current_axis = depth % 3;
-      unsigned int next_axis = (depth + 1) % 3;
-      aabb_queue.pop();
-
-      // create quad
-      glm::vec3 span1 = current_aabb.extent;
-      glm::vec3 span2 = current_aabb.extent;
-      glm::vec3 quad_position = current_aabb.position;
-      quad_position[current_axis] = current_split;
-      span1[current_axis] = 0.0;
-      span2[current_axis] = 0.0;
-      span2[next_axis] = -span2[next_axis];
-      CgQuad* new_quad = new CgQuad(1000 + index, quad_position, span1, span2);
-      m_renderer->init(new_quad);
-      m_kd_tree_mesh.push_back(new_quad);
-
-      // push children of current aabb
-      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;
-    }
+  std::vector<double> splits;
+  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;
+  unsigned int depth = 0;
+  while(index < splits.size() && depth < max_depth_u) {
+    CgAABB current_aabb = aabb_queue.front();
+    double current_split = splits[index];
+    unsigned int current_axis = depth % 3;
+    unsigned int next_axis = (depth + 1) % 3;
+    aabb_queue.pop();
+
+    // create quad
+    glm::vec3 span1 = current_aabb.extent;
+    glm::vec3 span2 = current_aabb.extent;
+    glm::vec3 quad_position = current_aabb.position;
+    quad_position[current_axis] = current_split;
+    span1[current_axis] = 0.0;
+    span2[current_axis] = 0.0;
+    span2[next_axis] = -span2[next_axis];
+    CgQuad* new_quad = new CgQuad(1000 + index, quad_position, span1, span2);
+    m_renderer->init(new_quad);
+    m_kd_tree_mesh.push_back(new_quad);
+
+    // push children of current aabb
+    std::pair<CgAABB, CgAABB> children = current_aabb.split(current_axis, current_split);
+    aabb_queue.push(children.first);
+    aabb_queue.push(children.second);
+
+    index += 1;
+    unsigned int v = (0x01 << (depth + 1)) - 1;
+    if(index >= v) depth += 1;
   }
 }
 
 void CgSceneControl::handlePickEvent(float x, float y) {
   calculatePickRay(x, y);
-  if(m_pointcloud != nullptr) {
-    std::vector<glm::vec3> points = m_select_ray->getVertices();
-    glm::vec3 start = points[0];
-    glm::vec3 direction = points[1] - start;
-    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));
-      m_renderer->init(m_pointcloud);
-      m_renderer->redraw();
-    }
-  }
-}
-
-void CgSceneControl::handleIncrementalSimplifyEvent(const CgSimplifyEvent& ev) {
-  unsigned int MAX_POINTS = 50;
-  std::cout << ev << std::endl;
-  // get current points
-  std::vector<glm::vec3> original_points = m_pointcloud->getVertices();
-  std::vector<glm::vec3> remaining_points;
-  std::vector<glm::vec3> current_cluster;
-  std::vector<bool> used;
-  used.resize(original_points.size(), false);
-  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_scalings;
-  std::vector<unsigned int> new_splat_indices;
-  unsigned int next_seed = 0;
-
-  while(true) { // while there are still points to consider...
-    std::vector<unsigned int> nh = CgPointWrangler::getNearestNeighborsFast(
-      next_seed, MAX_POINTS + 1,
-      original_points,
-      CgPointMath::calculateAABB(original_points)
-    );
-    unsigned int i;
-    PCA cluster_pca;
-    current_cluster.push_back(original_points[next_seed]);
-    used[next_seed] = true;
-    for(i = 1; i < MAX_POINTS; i++) {
-      current_cluster.push_back(original_points[nh[i]]);
-      used[nh[i]] = true;
-      PCA current_pca = CgPointMath::calculatePCA(current_cluster);
-      float variance = current_pca.eval0 / (current_pca.eval0 + current_pca.eval1 + current_pca.eval2);
-      if(variance > ev.getMaxVariance()) break;
-      cluster_pca = current_pca;
-    }
-    next_seed = i + 1;
-    current_cluster.clear();
-    // push new point information
-    new_points.push_back(cluster_pca.centroid);
-    new_normals.push_back(cluster_pca.evec0);
-    new_colors.push_back(glm::vec3(0.5));
-    new_splat_orientations.push_back(glm::lookAt(
-      cluster_pca.centroid,
-      cluster_pca.centroid - cluster_pca.evec0,
-      CgPointMath::getPerpendicularVector(cluster_pca.evec0)
-    ));
-    new_splat_scalings.push_back(glm::vec2(0.02));
-    new_splat_indices.push_back(new_points.size() - 1);
-
-    // copy unused points over
-    for(unsigned int i = 0; i < used.size(); i++) {
-      if(!used[i]) remaining_points.push_back(original_points[i]);
-      if(i == next_seed) next_seed = remaining_points.size() - 1;
-      used[i] = false;
-    }
-    original_points = remaining_points;
-
-    if(original_points.size() == 0) break;
-    // build new kd-tree in remaining points.
-    CgPointWrangler::buildKdTree(original_points.data(), original_points.data() + original_points.size(), 0);
-    remaining_points.clear();
-    used.resize(original_points.size(), false);
-  }
-  std::cout << new_points.size() << " clusters" << std::endl;
-  // 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_scalings,
-    new_splat_indices
+  if(m_pointcloud == nullptr) return;
+  std::vector<glm::vec3> points = m_select_ray->getVertices();
+  glm::vec3 start = points[0];
+  glm::vec3 direction = points[1] - start;
+  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
   );
-}
-
-void CgSceneControl::handleHierarchicalSimplifyEvent(const CgSimplifyEvent& ev) {
-  std::cout << "hierarchical" << std::endl;
+  if(picked_index == -1) return;
+  unsigned int picked_index_u = picked_index;
+  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_u}, glm::vec3(1.0, 1.0, 0.0));
+  m_renderer->init(m_pointcloud);
+  m_renderer->redraw();
 }
 
 void CgSceneControl::handleEvent(CgBaseEvent* e) {
@@ -369,15 +289,58 @@ void CgSceneControl::handleEvent(CgBaseEvent* e) {
 
     if(e->getType() & Cg::CgSimplifyEvent && m_pointcloud != nullptr) {
       CgSimplifyEvent* ev = (CgSimplifyEvent*) e;
+      std::cout << *ev << std::endl;
       if(ev->doReset()) {
         CgPointCloud* copy = new CgPointCloud(*m_original_pointcloud);
         delete m_pointcloud;
         m_pointcloud = copy;
-      } else if(ev->getClusteringType() == Cg::CgIncremental) {
-        run_timed("incr. simplify", 1, [&](){handleIncrementalSimplifyEvent(*ev);});
       } else {
-        handleHierarchicalSimplifyEvent(*ev);
+        std::vector<glm::vec3> original_points = m_pointcloud->getVertices();
+        std::vector<PCA> pca = ev->getClusteringType() == Cg::CgIncremental
+          ? CgPointWrangler::performIncrementalSimplification(
+            original_points, 
+            ev->getMaxPoints(), 
+            ev->getMaxVariance(),
+            ev->getMaxDiameter()
+          )
+          : CgPointWrangler::performHierarchicalSimplification(
+            original_points,
+            ev->getMaxPoints(), 
+            ev->getMaxVariance(),
+            ev->getMaxDiameter()
+          );
+
+        // use new information to create a new pointcloud
+        std::vector<glm::vec3> new_points(pca.size());
+        std::vector<glm::mat4> new_splat_orientations(pca.size());
+        std::vector<glm::vec2> new_splat_scalings(pca.size());
+        std::vector<glm::vec3> new_colors(pca.size());
+        CgPointWrangler::buildKdTree(pca.data(), pca.data() + pca.size(), 0);
+        for(unsigned int i = 0; i < pca.size(); i++) new_points[i] = pca[i].centroid;
+        std::vector<std::vector<unsigned int>> neighbourhoods = CgPointWrangler::getKNeighbourhoods(5, new_points);
+        std::vector<std::pair<unsigned int, unsigned int>> spanning_tree = CgPointWrangler::makeNaiveSpanningTree(neighbourhoods, pca);
+        std::vector<glm::vec3> new_normals = CgPointWrangler::alignNormals(glm::vec3(0.0,1.0,0.0), spanning_tree, pca);
+
+        for(unsigned int i = 0; i < new_points.size(); i++) {
+          new_colors[i] = CgPointMath::normalToColor(new_normals[i]);
+          new_splat_orientations[i] = glm::lookAt(
+            new_points[i],
+            new_points[i] - new_normals[i],
+            CgPointMath::getPerpendicularVector(new_normals[i])
+          );
+          new_splat_scalings[i] = glm::vec2(std::max(pca[i].radius, 0.02));
+        }
+
+        delete m_pointcloud;
+        m_pointcloud = new CgPointCloud(
+          new_points,
+          new_normals,
+          new_colors,
+          new_splat_orientations,
+          new_splat_scalings
+        );
       }
+
       m_renderer->init(m_pointcloud);
       m_renderer->redraw();
     }
@@ -408,11 +371,13 @@ void CgSceneControl::handleEvent(CgBaseEvent* e) {
 
           std::vector<glm::vec3> positions;
           std::vector<std::vector<unsigned int>> neighbourhoods;
+          std::vector<std::vector<glm::vec3>> clusters;
+          std::vector<PCA> pca;
+          std::vector<std::pair<unsigned int, unsigned int>> spanning_tree;
           std::vector<glm::vec3> normals;
           std::vector<glm::vec3> colors;
           std::vector<glm::mat4> splat_orientations;
           std::vector<glm::vec2> splat_scalings;
-          std::vector<unsigned int> splat_indices;
           loader->getPositionData(positions);
 
           std::cout << "loaded " << positions.size() << " points" << std::endl;
@@ -423,22 +388,32 @@ void CgSceneControl::handleEvent(CgBaseEvent* e) {
           run_timed("get neighbourhoods", 1, 
             [&](){neighbourhoods = CgPointWrangler::getKNeighbourhoods(15, positions);}
           );
-          run_timed("est. normals", 1, 
-            [&](){normals = CgPointWrangler::estimateNormals(positions, neighbourhoods);}
+          run_timed("get clusters", 1, 
+            [&](){clusters = CgPointWrangler::kNeighbourhoodsToClusters(positions, neighbourhoods);}
+          );
+          run_timed("do PCA", 1,
+            [&](){pca = CgPointWrangler::performPcaOnClusters(clusters);}
+          );
+          // we don't want to use centroids, but the actual points after loading.
+          for(unsigned int i = 0; i < pca.size(); i++) pca[i].centroid = positions[i];
+
+          run_timed("get spanning tree", 1,
+            [&](){spanning_tree = CgPointWrangler::makeNaiveSpanningTree(neighbourhoods, pca);}
+          );
+          run_timed("align normals", 1, 
+            [&](){normals = CgPointWrangler::alignNormals(glm::vec3(0.0, 1.0, 0.0), spanning_tree, pca);}
           );
           colors.resize(positions.size());
           splat_orientations.resize(positions.size());
           splat_scalings.resize(positions.size());
-          splat_indices.resize(positions.size());
           for(unsigned int i = 0; i < positions.size(); i++) {
-            colors[i] = (normals[i] + glm::vec3(1.0)) * .5;
+            colors[i] = CgPointMath::normalToColor(normals[i]);
             splat_orientations[i] = glm::lookAt(
               positions[i],
               positions[i] - normals[i],
               CgPointMath::getPerpendicularVector(normals[i])
             );
             splat_scalings[i] = glm::vec2(0.01, 0.01);
-            splat_indices[i] = i;
           }
 
           m_pointcloud = new CgPointCloud(
@@ -446,9 +421,14 @@ void CgSceneControl::handleEvent(CgBaseEvent* e) {
             normals,
             colors,
             splat_orientations,
-            splat_scalings,
-            splat_indices
+            splat_scalings
           );
+
+          std::cout 
+          << "dx=" << 2*m_pointcloud->getAABB().extent.x
+          << " dy=" << 2*m_pointcloud->getAABB().extent.y
+          << " dz=" << 2*m_pointcloud->getAABB().extent.z
+          << std::endl;
           
           // save away a copy
           m_original_pointcloud = new CgPointCloud(*m_pointcloud);
@@ -543,7 +523,7 @@ void CgSceneControl::renderObjects() {
     }
 
     if(m_kd_tree_mesh.size() != 0) {
-      for(int i = 0; i < m_kd_tree_mesh.size(); i++) {
+      for(unsigned int i = 0; i < m_kd_tree_mesh.size(); i++) {
         m_current_transformation = glm::translate(m_current_transformation, -m_center);
         glm::mat4 mv_matrix = m_lookAt_matrix * m_current_transformation;
         m_renderer->render(m_kd_tree_mesh[i], mv_matrix, m_proj_matrix);
@@ -559,16 +539,14 @@ void CgSceneControl::renderObjects() {
 
       //const std::vector<glm::vec3>& vertex_colors = m_pointcloud->getVertexColors();
       const std::vector<glm::vec3>& vertices = m_pointcloud->getVertices();
-      const std::vector<unsigned int>& splat_indices= m_pointcloud->getSplatIndices();
       const std::vector<glm::mat4>& orientations = m_pointcloud->getSplatOrientations();
       const std::vector<glm::vec2>& scalings = m_pointcloud->getSplatScalings();
 
       glm::mat4 mv_matrix;
 
-      for(unsigned int i = 0; i < splat_indices.size(); i++) {
-        unsigned int splat_index = splat_indices[i];
-        glm::mat4 splat_orientation = orientations[splat_index];
-        glm::vec2 splat_scaling = scalings[splat_index];
+      for(unsigned int i = 0; i < vertices.size(); i++) {
+        glm::mat4 splat_orientation = orientations[i];
+        glm::vec2 splat_scaling = scalings[i];
 
         m_current_transformation = glm::translate(m_current_transformation, -m_center);
         m_current_transformation = m_current_transformation * glm::inverse(splat_orientation);
diff --git a/CgSceneGraph/CgSceneControl.h b/CgSceneGraph/CgSceneControl.h
index c0560364764fcceaf54c4485790b67bc93970d2f..12a6b9f4428f6eecae786b6cb102acb04ab894e2 100644
--- a/CgSceneGraph/CgSceneControl.h
+++ b/CgSceneGraph/CgSceneControl.h
@@ -56,8 +56,6 @@ private:
     void calculateSingularValueDecomposition();
     void createKdTreeViz(int max_depth);
     void handlePickEvent(float x, float y);
-    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;
diff --git a/CgUtils/Timer.h b/CgUtils/Timer.h
index 3e13464a880b4b04abe70668627e09e0c8db1862..d294d3562d12feddf7f754603bdac9bfbbd42ebf 100644
--- a/CgUtils/Timer.h
+++ b/CgUtils/Timer.h
@@ -12,7 +12,7 @@
 // amount times, printing the duration in milliseconds to stdout afterwards
 inline void run_timed(std::string tag, unsigned int amount, std::function<void ()> fn) {
     auto t1 = std::chrono::high_resolution_clock::now();
-    for(int i = 0; i < amount; i++) fn();
+    for(unsigned 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 << "average for " << tag << ": "<< ms_double.count() / double(amount) << " ms" << std::endl;