Skip to content

Commit

Permalink
Added plane fitting
Browse files Browse the repository at this point in the history
  • Loading branch information
@Chonnik
Chonnik committed Feb 11, 2023
1 parent ccfd308 commit fd91562
Show file tree
Hide file tree
Showing 4 changed files with 79 additions and 56 deletions.
38 changes: 30 additions & 8 deletions src/Surface/diffgeo.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,15 +30,15 @@ void setup_prolongation_matrix(pmp::SurfaceMesh &mesh, SparseMatrix &P) {
Eigen::VectorXd w;

std::vector<Triplet> 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++;
}
Expand Down Expand Up @@ -80,17 +80,17 @@ void setup_face_point_properties(pmp::SurfaceMesh &mesh, unsigned int min_point)
std::vector<Eigen::Triplet<double>> 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 {
Expand All @@ -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);
Expand Down Expand Up @@ -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));
}


5 changes: 5 additions & 0 deletions src/Surface/diffgeo.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);


//=============================================================================
52 changes: 22 additions & 30 deletions src/SurfaceViewer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -268,42 +268,34 @@ void Viewer::process_imgui()
if (ImGui::Button("Check for non-planarity")) {
auto vpoint = mesh_.get_vertex_property<Point>("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;
}
}

Expand Down
40 changes: 22 additions & 18 deletions src/run_surface_tests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -513,35 +513,39 @@ 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_" +
std::to_string(i) + ".off";
mesh.read(meshname);
std::cout << meshname << std::endl;
auto vpoint = mesh.get_vertex_property<Point>("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();
}
Expand Down

0 comments on commit fd91562

Please sign in to comment.