diff --git a/CgMath/CgCovariance.h b/CgMath/CgCovariance.h
deleted file mode 100644
index 8157154776cb36500cf438a2335bead675b4fcc8..0000000000000000000000000000000000000000
--- a/CgMath/CgCovariance.h
+++ /dev/null
@@ -1,33 +0,0 @@
-#ifndef CGCOVARIANCE_H
-#define CGCOVARIANCE_H
-
-#include <vector>
-#include "glm/glm.hpp"
-#include "Eigen/Dense"
-
-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 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;
-}
-
-#endif // CGCOVARIANCE_H
\ No newline at end of file
diff --git a/CgMath/CgDecomposition.h b/CgMath/CgDecomposition.h
new file mode 100644
index 0000000000000000000000000000000000000000..fd262a871d2671603333eb2d6e98b3ecbf74f870
--- /dev/null
+++ b/CgMath/CgDecomposition.h
@@ -0,0 +1,66 @@
+#ifndef CGCOVARIANCE_H
+#define CGCOVARIANCE_H
+
+#include <vector>
+#include <iostream>
+#include "glm/glm.hpp"
+#include "Eigen/Dense"
+#include "CgMath/CgEigenDecomposition3x3.h"
+#include <CgMath/Eigen/SVD>
+#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/CgSceneGraph/CgPointCloud.cpp b/CgSceneGraph/CgPointCloud.cpp
index 0b774e119fb0384baa1338cea85de682928c2e93..4a64ba93eaff5377cfc7c0f77e07b32ae91d3df9 100644
--- a/CgSceneGraph/CgPointCloud.cpp
+++ b/CgSceneGraph/CgPointCloud.cpp
@@ -61,13 +61,14 @@ void CgPointCloud::init( std::string filename, bool cheat_normals) {
ObjLoader loader;
loader.load(filename);
loader.getPositionData(m_vertices);
+ std::cout << "loaded " << m_vertices.size() << " points" << std::endl;
+ std::cout << "using " << std::thread::hardware_concurrency() << " threads" << std::endl;
m_centroid = getCentroid(m_vertices);
- buildKdTree(&m_vertices.front(), &(*(m_vertices.end() - 1)), 0);
+ run_timed("build kd", 1, [&](){buildKdTree(&m_vertices.front(), &(*(m_vertices.end() - 1)), 0);});
calculateAABB();
- // run_timed(1, [&]() {estimateNormals();});
- estimateNormals();
- std::cout << "got " << m_vertices.size() << " points" << std::endl;
+ run_timed("est. normals", 1, [&](){estimateNormals();});
+ //estimateNormals();
// calculate local coordinate system for splats
// (arbitrary orientation of ellipse in plane)
@@ -80,12 +81,6 @@ void CgPointCloud::init( std::string filename, bool cheat_normals) {
for(unsigned int i = 0;i < m_vertices.size(); i++) {
m_vertex_colors.push_back(glm::vec3(0.0,1.0,0.0));
}
-
- unsigned int k = 50;
- // "fast" if k is relatively small, if k is on the order
- // of m_vertices.size(), this is slower than brute force.
- std::vector<unsigned int> neighbors = getNearestNeighborsFast(0, k);
- setColors(neighbors, glm::vec3(0.0, 0.0, 1.0));
}
void CgPointCloud::buildKdTree(glm::vec3* begin, glm::vec3* end, int depth) {
@@ -112,6 +107,7 @@ 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();
// find point with highest y-value
//std::vector<unsigned int> spanning_tree;
@@ -133,22 +129,26 @@ void CgPointCloud::estimateNormals() {
// build minimal spanning tree (for each index, record the index we came from)
// traverse the tree, search neighborhoods
- for(unsigned int i = 0; i < m_vertices.size(); i++) {
- //std::cout << "alive, now " << i << std::endl;
+ std::vector<std::thread> threads;
+ threads.reserve(thread_count);
+ auto worker = [&](
+ const unsigned int start_index,
+ const unsigned int count
+ ) {
+ 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);
- std::vector<glm::vec3> nh;
- nh.resize(neighbourhood_size);
- for(unsigned int j = 0; j < neighbourhood_size; j++) {
- nh[j] = m_vertices[nh_indices[j]];
- }
- glm::mat3 cov_mat = getCovarianceMatrix(nh);
- CgEigenDecomposition3x3 eigen(cov_mat);
- // The eigenvalues/vectors are sorted in increasing order, and we want the smallest one.
- glm::vec3 normal = eigen.getEigenvectors()[0];
- m_vertex_normals[i] = glm::dot(normal, m_vertices[i] - m_centroid) > 0.0
- ? normal
- : -normal;
- }
+ 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;
+ }
+ };
+ 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();
// select normal pointing roughly in the same direction as parent's normal
}
diff --git a/CgSceneGraph/CgPointCloud.h b/CgSceneGraph/CgPointCloud.h
index 0103dc03d17dfc28b1c311148892edb073f6692c..df391353c74345b9382b9d62b0f08f14fab329aa 100644
--- a/CgSceneGraph/CgPointCloud.h
+++ b/CgSceneGraph/CgPointCloud.h
@@ -11,13 +11,12 @@
#include <utility>
#include <functional>
#include <tuple>
-#include <set>
+#include <thread>
#include "CgBase/CgBasePointCloud.h"
#include "CgSceneGraph/CgAABB.h"
#include "CgBase/CgEnums.h"
#include "CgUtils/ObjLoader.h"
-#include "CgMath/CgCovariance.h"
-#include "CgMath/CgEigenDecomposition3x3.h"
+#include "CgMath/CgDecomposition.h"
#include <glm/glm.hpp>
#include <glm/gtx/norm.hpp>
#include <glm/gtc/matrix_transform.hpp>
diff --git a/CgSceneGraph/CgSceneControl.cpp b/CgSceneGraph/CgSceneControl.cpp
index 870c95c9b50bf832d574658a3b59f5669c09e41a..7c9a79d7bde71f0038a9c8d659b97e1367b897f8 100644
--- a/CgSceneGraph/CgSceneControl.cpp
+++ b/CgSceneGraph/CgSceneControl.cpp
@@ -21,11 +21,7 @@
#include <glm/gtc/matrix_transform.hpp>
#include <glm/ext.hpp>
#include "CgUtils/ObjLoader.h"
-#include "CgMath/CgEigenDecomposition3x3.h"
#include <string>
-#include <CgMath/Eigen/SVD>
-#include <CgMath/Eigen/Core>
-using namespace Eigen;
CgSceneControl::CgSceneControl()
{
@@ -119,34 +115,6 @@ void CgSceneControl::calculatePickRay(float x, float y) {
m_renderer->init(m_select_ray);
}
-
-void CgSceneControl::calculateEigenDecomposition3x3() {
-
- std::vector<glm::vec3> testvecs{
- glm::vec3(1.0, 1.0, 1.0),
- glm::vec3(1.0, 1.0, -1.0),
- glm::vec3(1.0, -1.0, 1.0),
- glm::vec3(-1.0, 4.3, 1.0),
- glm::vec3(1.0, -2.15, -1.0),
- glm::vec3(-.45, 1.0, -1.0),
- glm::vec3(-1.0, -1.0, 1.0),
- glm::vec3(-1.0, -2.0, -1.0),
- };
-
- glm::mat3 testmat = getCovarianceMatrix(testvecs);
-
- std::cout << "CovMat is\n" << glm::to_string(testmat) << std::endl;
-
- CgEigenDecomposition3x3 eigen(testmat);
- glm::mat3 eigenvectors= eigen.getEigenvectors();
- glm::vec3 eigenvalues = eigen.getEigenvalues();
-
- std::cout << "evec0: " << glm::to_string(eigenvectors[0]) << std::endl;
- std::cout << "evec1: " << glm::to_string(eigenvectors[1]) << std::endl;
- std::cout << "evec2: " << glm::to_string(eigenvectors[2]) << std::endl;
- std::cout << "eval:" << glm::to_string(eigenvalues) << std::endl;
-}
-
void CgSceneControl::calculateSingularValueDecomposition() {
//init some arbitrary test matrix
MatrixXd C(2,3);
diff --git a/CgSceneGraph/CgSceneControl.h b/CgSceneGraph/CgSceneControl.h
index 741aefb55b858e3db902c78fa0f0d2bc64ad7c1a..4fd2c38cb1e329780f9356b6e6de433c7dda4c86 100644
--- a/CgSceneGraph/CgSceneControl.h
+++ b/CgSceneGraph/CgSceneControl.h
@@ -39,7 +39,6 @@ private:
CgPolyLine* m_select_ray;
std::vector<CgQuad*> m_kd_tree_mesh;
- void calculateEigenDecomposition3x3();
void calculateSingularValueDecomposition();
void createKdTreeViz(int max_depth);
void handlePickEvent(float x, float y);
diff --git a/CgUtils/Timer.h b/CgUtils/Timer.h
index 9f1a4cced61f02fda6c3156047f6cbedc2bdc0b3..1aee94726c44d5e12498ddeefeafe3eb2955c5c3 100644
--- a/CgUtils/Timer.h
+++ b/CgUtils/Timer.h
@@ -9,12 +9,12 @@
// 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(unsigned int amount, std::function<void ()> fn) {
+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: " << ms_double.count() / double(amount) << " ms" << std::endl;
+ std::cout << "average for " << tag << ": "<< ms_double.count() / double(amount) << " ms" << std::endl;
}
#endif // TIMER_H
\ No newline at end of file
diff --git a/ExerciseVC.pro b/ExerciseVC.pro
index c7f87359d4602a70143360916e94010ef623de4d..c623731de3ce1f8fb608023e43b299c1a359356a 100644
--- a/ExerciseVC.pro
+++ b/ExerciseVC.pro
@@ -1,7 +1,7 @@
QT += core gui opengl widgets
TEMPLATE = app
TARGET = PointViewer
-QMAKE_CXXFLAGS += -g -std=c++20
+QMAKE_CXXFLAGS += -std=c++20 -pthread
CONFIG += c++20
SOURCES += main.cpp \
@@ -37,7 +37,7 @@ HEADERS += \
CgEvents/CgKdTreeEvent.h \
CgEvents/CgSplatEvent.h \
CgMath/CgEigenDecomposition3x3.h \
- CgMath/CgCovariance.h \
+ CgMath/CgDecomposition.h \
CgMath/Eigen/Core \
CgMath/Eigen/Eigen \
CgMath/Eigen/SVD \