diff --git a/.vscode/launch.json b/.vscode/launch.json
new file mode 100644
index 0000000000000000000000000000000000000000..224af64df30e1dcd38224cdaaf40e2f5321a6d71
--- /dev/null
+++ b/.vscode/launch.json
@@ -0,0 +1,16 @@
+{
+ // Use IntelliSense to learn about possible attributes.
+ // Hover to view descriptions of existing attributes.
+ // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387
+ "version": "0.2.0",
+ "configurations": [
+ {
+ "type": "lldb",
+ "request": "launch",
+ "name": "Debug",
+ "program": "${workspaceFolder}/<your program>",
+ "args": [],
+ "cwd": "${workspaceFolder}"
+ }
+ ]
+}
\ No newline at end of file
diff --git a/CgMath/CgDecomposition.h b/CgMath/CgDecomposition.h
index fd262a871d2671603333eb2d6e98b3ecbf74f870..7fdbfc3325068a53d3d5a84b6632babc766d5226 100644
--- a/CgMath/CgDecomposition.h
+++ b/CgMath/CgDecomposition.h
@@ -5,8 +5,6 @@
#include <iostream>
#include "glm/glm.hpp"
#include "Eigen/Dense"
-#include "CgMath/CgEigenDecomposition3x3.h"
-#include <CgMath/Eigen/SVD>
#include <CgMath/Eigen/Core>
using namespace Eigen;
diff --git a/CgMath/CgEigenDecomposition3x3.cpp b/CgMath/CgEigenDecomposition3x3.cpp
deleted file mode 100644
index 75735712e0dc549ca8b99fb1ed54eae161b4a184..0000000000000000000000000000000000000000
--- a/CgMath/CgEigenDecomposition3x3.cpp
+++ /dev/null
@@ -1,41 +0,0 @@
-#include "CgEigenDecomposition3x3.h"
-
-using namespace Eigen;
-
-CgEigenDecomposition3x3::CgEigenDecomposition3x3(glm::mat3 arg) {
- Eigen::Matrix3d A;
- A(0,0) = arg[0][0]; A(0,1) = arg[0][1]; A(0,2) = arg[0][2];
- A(1,0) = arg[1][0]; A(1,1) = arg[1][1]; A(1,2) = arg[1][2];
- A(2,0) = arg[2][0]; A(2,1) = arg[2][1]; A(2,2) = arg[2][2];
-
- // the covariance matrix is a symmetric real matrix and therefore
- // self-adjoint
- Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> es;
- es.compute(A);
-
- // The eigenvalues of a selfadjoint matrix are always real.
- Eigen::Matrix3d V = es.eigenvectors();
- Eigen::Vector3d res = es.eigenvalues();
-
- // eigenvectors are returned in the COLUMNS of
- // the Matrix
- // V(0,1) is the second element of the first ROW
- for(unsigned int row = 0; row < 3; row++) {
- for(unsigned int col = 0; col < 3; col++) {
- m_eigenvectors[col][row] = V(row, col);
- }
- }
-
- m_eigenvalues[0] = res(0);
- m_eigenvalues[1] = res(1);
- m_eigenvalues[2] = res(2);
-}
-
-
-glm::mat3 CgEigenDecomposition3x3::getEigenvectors() {
- return m_eigenvectors;
-}
-
-glm::vec3 CgEigenDecomposition3x3::getEigenvalues() {
- return m_eigenvalues;
-}
diff --git a/CgMath/CgEigenDecomposition3x3.h b/CgMath/CgEigenDecomposition3x3.h
deleted file mode 100644
index d9f481b09abff51530fac624fc509e140ea21365..0000000000000000000000000000000000000000
--- a/CgMath/CgEigenDecomposition3x3.h
+++ /dev/null
@@ -1,20 +0,0 @@
-#ifndef CGEIGENDECOMPOSITION3X3_H
-#define CGEIGENDECOMPOSITION3X3_H
-
-#include "glm/glm.hpp"
-#include "CgMath/Eigen/Dense"
-#include "CgMath/Eigen/Eigenvalues"
-#include "glm/glm.hpp"
-#include <iostream>
-
-class CgEigenDecomposition3x3 {
-public:
- CgEigenDecomposition3x3(glm::mat3 arg);
- glm::mat3 getEigenvectors();
- glm::vec3 getEigenvalues();
-private:
- glm::mat3 m_eigenvectors;
- glm::vec3 m_eigenvalues;
-};
-
-#endif // CGEIGENDECOMPOSITION3X3_H
diff --git a/CgMath/CgGMST.h b/CgMath/CgGMST.h
new file mode 100644
index 0000000000000000000000000000000000000000..9905ee46641f8cb1c7e6f5e65c411edb0036bf64
--- /dev/null
+++ b/CgMath/CgGMST.h
@@ -0,0 +1,16 @@
+
+#ifndef CGGMST_H
+#define CGGMST_H
+
+// benchmarking utility
+// takes a closure without params that returns void and runs it
+// amount times, printing the duration in milliseconds to stdout afterwards
+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();
+ 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;
+}
+
+#endif // CGGMST_H
\ No newline at end of file
diff --git a/CgSceneGraph/CgPointCloud.cpp b/CgSceneGraph/CgPointCloud.cpp
index 4a64ba93eaff5377cfc7c0f77e07b32ae91d3df9..d04c0693ed9de697ed2265f3eff62819e4a09e22 100644
--- a/CgSceneGraph/CgPointCloud.cpp
+++ b/CgSceneGraph/CgPointCloud.cpp
@@ -1,5 +1,4 @@
#include "CgPointCloud.h"
-#include "CgUtils/Timer.h"
// squared distance between a point and a ray
float sqrPointRayDist(glm::vec3 point, glm::vec3 origin, glm::vec3 direction) {
@@ -61,6 +60,8 @@ void CgPointCloud::init( std::string filename, bool cheat_normals) {
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);
@@ -78,9 +79,9 @@ void CgPointCloud::init( std::string filename, bool cheat_normals) {
//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.push_back(glm::vec3(0.0,1.0,0.0));
- }
+ //for(unsigned int i = 0; i < m_vertices.size(); i++) {
+ // m_vertex_colors[i] = glm::vec3(0.0,1.0,0.0);
+ //}
}
void CgPointCloud::buildKdTree(glm::vec3* begin, glm::vec3* end, int depth) {
@@ -106,43 +107,52 @@ void CgPointCloud::buildKdTree(glm::vec3* begin, glm::vec3* end, int depth) {
}
void CgPointCloud::estimateNormals() {
- unsigned int neighbourhood_size = 25;
- unsigned int thread_count = std::thread::hardware_concurrency();
+ unsigned int neighbourhood_size = 15;
// find point with highest y-value
- //std::vector<unsigned int> spanning_tree;
- //spanning_tree.resize(m_vertices.size());
- m_vertex_normals.clear();
- m_vertex_normals.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++) {
- // float curr = m_vertices[i].y;
- // if(max_y < curr) {
- // max_y = curr;
- // max_y_index = i;
- // }
- //}
-
- // the extremal point is selected as root (has no parent)
- //spanning_tree[max_y_index] = max_y_index;
-
- // build minimal spanning tree (for each index, record the index we came from)
- // traverse the tree, search neighborhoods
+ 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++) {
- std::vector<unsigned int> nh_indices = getNearestNeighborsFast(i, neighbourhood_size);
- for(unsigned int j = 0; j < neighbourhood_size; j++) nh[j] = m_vertices[nh_indices[j]];
- glm::vec3 normal = estimateNormalFromPCA(nh);
- m_vertex_normals[i] = glm::dot(normal, m_vertices[i] - m_centroid) > 0.0 ? normal : -normal;
+ 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;
@@ -150,7 +160,58 @@ void CgPointCloud::estimateNormals() {
threads.push_back(std::thread(worker, t * slice_size, slice_size + slice_rest));
for(t = 0; t < thread_count; t++) threads[t].join();
- // select normal pointing roughly in the same direction as parent's normal
+ 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;
+ m_vertex_colors[curr] = glm::vec3(
+ (normal.x + 1.0) * .5,
+ (normal.y + 1.0) * .5,
+ (normal.z + 1.0) * .5
+ );
+
+ // 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) {
@@ -235,7 +296,7 @@ std::vector<unsigned int> CgPointCloud::getNearestNeighborsFast(unsigned int cur
glm::vec3 curr_point = *(begin + half_size);
// insert current point into set
- candidates.push_back({glm::distance2(query, curr_point), begin + half_size - first});
+ 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);
@@ -249,7 +310,7 @@ std::vector<unsigned int> CgPointCloud::getNearestNeighborsFast(unsigned int cur
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.begin()).first;
+ double neighbourhood_radius = candidates[0].first;
if(
begin != begin + half_size
&& lower_aabb.position[axis] + lower_aabb.extent[axis] >= query[axis] - neighbourhood_radius
diff --git a/CgSceneGraph/CgPointCloud.h b/CgSceneGraph/CgPointCloud.h
index df391353c74345b9382b9d62b0f08f14fab329aa..fe5ef5e52b7aa0346150088dda30283e4a4d4b60 100644
--- a/CgSceneGraph/CgPointCloud.h
+++ b/CgSceneGraph/CgPointCloud.h
@@ -20,6 +20,7 @@
#include <glm/glm.hpp>
#include <glm/gtx/norm.hpp>
#include <glm/gtc/matrix_transform.hpp>
+#include "CgUtils/Timer.h"
class CgPointCloud : public CgBasePointCloud {
public:
diff --git a/CgUtils/Timer.h b/CgUtils/Timer.h
index 1aee94726c44d5e12498ddeefeafe3eb2955c5c3..799c6f80d5806e2d74beaa29ebcb92d0b8011547 100644
--- a/CgUtils/Timer.h
+++ b/CgUtils/Timer.h
@@ -9,7 +9,7 @@
// benchmarking utility
// takes a closure without params that returns void and runs it
// amount times, printing the duration in milliseconds to stdout afterwards
-void run_timed(std::string tag, unsigned int amount, std::function<void ()> fn) {
+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();
auto t2 = std::chrono::high_resolution_clock::now();
diff --git a/ExerciseVC.pro b/ExerciseVC.pro
index c623731de3ce1f8fb608023e43b299c1a359356a..e2a40a362350d07db812fe9be302d11013fd2850 100644
--- a/ExerciseVC.pro
+++ b/ExerciseVC.pro
@@ -10,7 +10,6 @@ SOURCES += main.cpp \
CgEvents/CgPickRayEvent.cpp \
CgEvents/CgKdTreeEvent.cpp \
CgEvents/CgSplatEvent.cpp \
- CgMath/CgEigenDecomposition3x3.cpp \
CgQtViewer/CGQtGLRenderWidget.cpp \
CgQtViewer/CgQtGui.cpp \
CgBase/CgObservable.cpp \
@@ -36,7 +35,6 @@ HEADERS += \
CgEvents/CgPickRayEvent.h \
CgEvents/CgKdTreeEvent.h \
CgEvents/CgSplatEvent.h \
- CgMath/CgEigenDecomposition3x3.h \
CgMath/CgDecomposition.h \
CgMath/Eigen/Core \
CgMath/Eigen/Eigen \