From fd9156223581bef8a613df38776de5fe60242b55 Mon Sep 17 00:00:00 2001 From: abunge Date: Sat, 11 Feb 2023 15:24:48 +0100 Subject: [PATCH] Added plane fitting --- src/Surface/diffgeo.cpp | 38 ++++++++++++++++++++++------ src/Surface/diffgeo.h | 5 ++++ src/SurfaceViewer.cpp | 52 +++++++++++++++++---------------------- src/run_surface_tests.cpp | 40 ++++++++++++++++-------------- 4 files changed, 79 insertions(+), 56 deletions(-) diff --git a/src/Surface/diffgeo.cpp b/src/Surface/diffgeo.cpp index d2807d3..0ddec9d 100644 --- a/src/Surface/diffgeo.cpp +++ b/src/Surface/diffgeo.cpp @@ -30,15 +30,15 @@ void setup_prolongation_matrix(pmp::SurfaceMesh &mesh, SparseMatrix &P) { Eigen::VectorXd w; std::vector tripletsA; - for (auto v : mesh.vertices()) { + for (auto v: mesh.vertices()) { tripletsA.emplace_back(v.idx(), v.idx(), 1.0); } unsigned int j = 0; - for (auto f : mesh.faces()) { + for (auto f: mesh.faces()) { w = area_weights[f]; unsigned int i = 0; - for (auto v : mesh.vertices(f)) { + for (auto v: mesh.vertices(f)) { tripletsA.emplace_back(nv + j, v.idx(), w(i)); i++; } @@ -80,17 +80,17 @@ void setup_face_point_properties(pmp::SurfaceMesh &mesh, unsigned int min_point) std::vector> trip; for (pmp::Face f: mesh.faces()) { - const int n = (int)mesh.valence(f); + const int n = (int) mesh.valence(f); poly.resize(n, 3); int i = 0; - for (pmp::Vertex v : mesh.vertices(f)) { + for (pmp::Vertex v: mesh.vertices(f)) { for (int h = 0; h < 3; h++) { poly.row(i)(h) = mesh.position(v)[h]; } i++; } if (min_point == Centroid) { - int val = (int)poly.rows(); + int val = (int) poly.rows(); w = Eigen::MatrixXd::Ones(val, 1); w /= (double) val; } else { @@ -108,9 +108,9 @@ void setup_face_point_properties(pmp::SurfaceMesh &mesh, unsigned int min_point) //-------------------------------------------------------------------------------- void find_area_minimizer_weights(const Eigen::MatrixXd &poly, - Eigen::VectorXd &weights) { + Eigen::VectorXd &weights) { - int val = (int)poly.rows(); + int val = (int) poly.rows(); Eigen::MatrixXd J(val, val); Eigen::VectorXd b(val); weights.resize(val); @@ -181,4 +181,26 @@ void find_area_minimizer_weights(const Eigen::MatrixXd &poly, } +//-------------------------------------------------------------------------------- + +void fit_plane_to_polygon(const Eigen::MatrixXd &poly, Eigen::Vector3d &normal, Eigen::Vector3d &origin) { + // Plane equation Ax+By+Cz+d =0; + int n = (int) poly.rows(); + Eigen::MatrixXd A(n, 3); + Eigen::VectorXd b(n); + + for (int i = 0; i < n; i++) { + A(i, 0) = poly(i, 0); + A(i, 1) = poly(i, 1); + A(i, 2) = 1.0; + b(i) = -poly(i, 2); + } + Eigen::MatrixXd AtA = A.transpose() * A; + Eigen::VectorXd Atb = A.transpose() * b; + Eigen::Vector3d x = AtA.llt().solve(Atb); + normal = Eigen::Vector3d(x(0),x(1),1.0); + normal.normalize(); + origin = Eigen::Vector3d(1.0,0.0,-x(0)-x(2)); +} + diff --git a/src/Surface/diffgeo.h b/src/Surface/diffgeo.h index b5b64b8..fc314f7 100644 --- a/src/Surface/diffgeo.h +++ b/src/Surface/diffgeo.h @@ -27,4 +27,9 @@ void setup_face_point_properties(pmp::SurfaceMesh &mesh, unsigned int min_point) void find_area_minimizer_weights(const Eigen::MatrixXd &poly, Eigen::VectorXd &weights); +//------------------------ Miscellaneous----------------------------------------------------------------- + +void fit_plane_to_polygon(const Eigen::MatrixXd &poly, Eigen::Vector3d &normal, Eigen::Vector3d &origin); + + //============================================================================= \ No newline at end of file diff --git a/src/SurfaceViewer.cpp b/src/SurfaceViewer.cpp index 57f2cef..0a7d35c 100644 --- a/src/SurfaceViewer.cpp +++ b/src/SurfaceViewer.cpp @@ -268,42 +268,34 @@ void Viewer::process_imgui() if (ImGui::Button("Check for non-planarity")) { auto vpoint = mesh_.get_vertex_property("v:point"); int ctr = 0; - double max_sum = 0.0; + double dist_sum = 0.0; for (auto f: mesh_.faces()) { - Point p0, p1, p2; - for(auto h :mesh_.halfedges(f)) { - p0 = vpoint[mesh_.from_vertex(h)]; - p1 = vpoint[mesh_.to_vertex(h)]; - p2 = vpoint[mesh_.to_vertex(mesh_.next_halfedge(h))]; - break; - } - Point n = normalize(cross(p0-p1,p2-p1)); - for (auto v: mesh_.vertices(f)) { + // fit plane to face + Eigen::MatrixXd poly(mesh_.valence(f),3); + int i =0; + for(auto v : mesh_.vertices(f)){ Point p = vpoint[v]; - if (abs(dot(n, p - p1)) > 0.0001) { - ctr += 1; - break; - } + poly(i,0) = p[0]; + poly(i,1) = p[1]; + poly(i,2) = p[2]; + i++; } - - double max = -1000000; - for(auto h :mesh_.halfedges(f)) { - p0 = vpoint[mesh_.from_vertex(h)]; - p1 = vpoint[mesh_.to_vertex(h)]; - p2 = vpoint[mesh_.to_vertex(mesh_.next_halfedge(h))]; - - n = normalize(cross(p0-p1,p2-p1)); - for (auto v: mesh_.vertices(f)) { - Point p = vpoint[v]; - if (abs(dot(n, p - p1)) > max) { - max = abs(dot(n, p - p1)); - } - } + Eigen::Vector3d n,o; + fit_plane_to_polygon(poly,n,o); + // compute mean distance to plane + double dist = 0.0; + for(auto v : mesh_.vertices(f)){ + Eigen::Vector3d p(vpoint[v][0],vpoint[v][1],vpoint[v][2]); + dist += abs(n.dot(p-o)); + } + dist /= (double)mesh_.valence(f); + if(dist > 0.0001){ + ctr +=1; } - max_sum += max; + dist_sum+=dist; } std::cout << "Nr. non planar faces: " << ctr << std::endl; - std::cout << "Mean maximum plane distance " << max_sum/mesh_.n_faces() << std::endl; + std::cout << "Mean plane distance " << dist_sum/(double)mesh_.n_faces() << std::endl; } } diff --git a/src/run_surface_tests.cpp b/src/run_surface_tests.cpp index 7f36b24..c318b2a 100644 --- a/src/run_surface_tests.cpp +++ b/src/run_surface_tests.cpp @@ -513,7 +513,7 @@ void write_convergence_data_csv(Function function, int lvl = 7, int start_lvl = } std::ofstream file_noise_quad(filename_); file_noise_quad - << "[AW11] l=2,[AW11] l=1,[AW11] l=0.5,[AW11] l=0.1,[dGBD20] l=2,[dGBD20] l=1,[dGBD20] l=0.5,[dGBD20] l=0.1,[BHKB20],[BBA21],[MKB08],Mean Max plane distance" + << "[AW11] l=2,[AW11] l=1,[AW11] l=0.5,[AW11] l=0.1,[dGBD20] l=2,[dGBD20] l=1,[dGBD20] l=0.5,[dGBD20] l=0.1,[BHKB20],[BBA21],[MKB08],Mean plane distance" << std::endl; for (int i = start_lvl; i < 9; i++) { std::string meshname = "../data/surface_meshes/sphere/hex2_noise_" + @@ -521,27 +521,31 @@ void write_convergence_data_csv(Function function, int lvl = 7, int start_lvl = mesh.read(meshname); std::cout << meshname << std::endl; auto vpoint = mesh.get_vertex_property("v:point"); - double max_sum = 0.0; + double dist_sum = 0.0; for (auto f: mesh.faces()) { - Point p0, p1, p2,n; - double max = -1000000; - for(auto h :mesh.halfedges(f)) { - p0 = vpoint[mesh.from_vertex(h)]; - p1 = vpoint[mesh.to_vertex(h)]; - p2 = vpoint[mesh.to_vertex(mesh.next_halfedge(h))]; - - n = normalize(cross(p0-p1,p2-p1)); - for (auto v: mesh.vertices(f)) { - Point p = vpoint[v]; - if (abs(dot(n, p - p1)) > max) { - max = abs(dot(n, p - p1)); - } - } + // fit plane to face + Eigen::MatrixXd poly(mesh.valence(f),3); + int j =0; + for(auto v : mesh.vertices(f)){ + Point p = vpoint[v]; + poly(j,0) = p[0]; + poly(j,1) = p[1]; + poly(j,2) = p[2]; + j++; + } + Eigen::Vector3d n,o; + fit_plane_to_polygon(poly,n,o); + // compute mean distance to plane + double dist = 0.0; + for(auto v : mesh.vertices(f)){ + Eigen::Vector3d p(vpoint[v][0],vpoint[v][1],vpoint[v][2]); + dist += abs(n.dot(p-o)); } - max_sum += max; + dist /= (double)mesh.valence(f); + dist_sum+=dist; } write_all_laplace_data(mesh, file_noise_quad, function, sphere_dist, l, m); - file_noise_quad << max_sum/(double)mesh.n_faces() << std::endl; + file_noise_quad << dist_sum/(double)mesh.n_faces() << std::endl; } file_noise_quad.close(); }