add implicit hole filling base code
This commit is contained in:
parent
1659341370
commit
8177b52927
|
|
@ -54,5 +54,248 @@ void fillHolesDumb(MyMesh &mesh) {
|
|||
}
|
||||
|
||||
|
||||
using namespace std ;
|
||||
using namespace MeshReconstruction;
|
||||
|
||||
/* ************** Implicit RBF ************** */
|
||||
|
||||
Implicit_RBF::Implicit_RBF(vector<float> alpha, vector<float> beta, vector<MyMesh::Point> center) : _alpha(alpha), _beta(beta), _center(center)
|
||||
{
|
||||
_n = _alpha.size() ;
|
||||
_d = _beta.size() ;
|
||||
if (_center.size() != _n)
|
||||
throw runtime_error("Inconsistent size of alpha and centers in Implicit_RBF constructor.");
|
||||
}
|
||||
|
||||
// Computation of the value of the implicit surface at point X
|
||||
float Implicit_RBF::val(MyMesh::Point X) const
|
||||
{
|
||||
float res = 0 ;
|
||||
// Computation of the sum of RBF at centers
|
||||
for(int i=0; i<_n; i++)
|
||||
{
|
||||
res += _alpha.at(i) * myphi((X-_center.at(i)).norm()) ;
|
||||
}
|
||||
// Computation of the polynomial part
|
||||
for(int j=0; j<_d; j++)
|
||||
{
|
||||
res += _beta.at(j) * myp(j, X) ;
|
||||
}
|
||||
return res ;
|
||||
}
|
||||
|
||||
/* ************** Hole Filling ************** */
|
||||
|
||||
Hole_Filling::Hole_Filling(MyMesh &mesh) : _mesh(mesh)
|
||||
{
|
||||
_mesh.add_property(_vprop, "vprop_flag");
|
||||
cout << "Starting hole filling ..." << endl ;
|
||||
}
|
||||
|
||||
// ***** Computation of boundary and its neighborhood
|
||||
|
||||
MyMesh::HalfedgeHandle Hole_Filling::find_boundary_edge()
|
||||
{
|
||||
MyMesh::HalfedgeIter he_it = _mesh.halfedges_begin();
|
||||
while ( (he_it != _mesh.halfedges_end()) && (!_mesh.is_boundary(*he_it)))
|
||||
{
|
||||
++he_it ;
|
||||
}
|
||||
if (he_it != _mesh.halfedges_end())
|
||||
return *he_it ;
|
||||
else
|
||||
throw std::runtime_error("Boundary HE does not exist") ;
|
||||
}
|
||||
|
||||
vector<MyMesh::VertexHandle> Hole_Filling::find_boundary(MyMesh::HalfedgeHandle heh)
|
||||
{
|
||||
MyMesh::HalfedgeHandle heh_ini = heh ;
|
||||
vector<MyMesh::VertexHandle> boundary ;
|
||||
|
||||
// Follow (and memorize) boundary edges
|
||||
do
|
||||
{
|
||||
boundary.push_back(_mesh.to_vertex_handle(heh));
|
||||
heh = _mesh.next_halfedge_handle(heh);
|
||||
} while (heh != heh_ini) ;
|
||||
return boundary ;
|
||||
}
|
||||
|
||||
void Hole_Filling::init_mark_boundary(const vector<MyMesh::VertexHandle> & bnd)
|
||||
{
|
||||
for (MyMesh::VertexIter v_it = _mesh.vertices_begin() ; v_it != _mesh.vertices_end(); ++v_it)
|
||||
_mesh.property(_vprop, *v_it) = false ;
|
||||
for (int i=0; i<bnd.size(); ++i)
|
||||
_mesh.property(_vprop, bnd.at(i)) = true ;
|
||||
}
|
||||
|
||||
vector<MyMesh::VertexHandle> Hole_Filling::next_neighbors(const vector<MyMesh::VertexHandle> & bnd)
|
||||
{
|
||||
// Visit bnd vertices to find and mark next circle
|
||||
vector<MyMesh::VertexHandle> next_bnd ;
|
||||
for (int i=0; i<bnd.size(); i++)
|
||||
{
|
||||
for (MyMesh::VertexVertexIter vv_it = _mesh.vv_iter(bnd.at(i));vv_it.is_valid();++vv_it)
|
||||
{
|
||||
if (_mesh.property(_vprop, *vv_it) == false) // new vertex
|
||||
{
|
||||
_mesh.property(_vprop, *vv_it) = true ;
|
||||
next_bnd.push_back(*vv_it);
|
||||
}
|
||||
}
|
||||
}
|
||||
return next_bnd ;
|
||||
}
|
||||
|
||||
// ***** Computation of RBF
|
||||
|
||||
pair<pair<Eigen::MatrixXd &, Eigen::VectorXd &>, vector<MyMesh::Point> &> Hole_Filling::compute_approx_mat(vector<MyMesh::VertexHandle> vlist)
|
||||
{
|
||||
const int n(vlist.size()), d(10) ;
|
||||
Eigen::MatrixXd & A = *(new Eigen::MatrixXd(3*n+d,3*n+d)) ;
|
||||
Eigen::MatrixXd Phi(3*n,3*n);
|
||||
Eigen::MatrixXd P(3*n,d);
|
||||
Eigen::VectorXd & B = *(new Eigen::VectorXd(3*n+d));
|
||||
|
||||
vector<MyMesh::Point> & pts_list = *(new vector<MyMesh::Point>);
|
||||
//Append vertices to pts_list
|
||||
for (int i=0; i<n; i++)
|
||||
{
|
||||
pts_list.push_back(_mesh.point(vlist.at(i))) ;
|
||||
}
|
||||
//Append vertices+normals to pts_list
|
||||
for (int i=0; i<n; i++)
|
||||
{
|
||||
pts_list.push_back(_mesh.point(vlist.at(i)) + _mesh.normal(vlist.at(i))) ;
|
||||
}
|
||||
//Append vertices-normals to pts_list
|
||||
for (int i=0; i<n; i++)
|
||||
{
|
||||
pts_list.push_back(_mesh.point(vlist.at(i)) - _mesh.normal(vlist.at(i))) ;
|
||||
}
|
||||
|
||||
int nn = pts_list.size() ;
|
||||
// Compute corresponding B vector (0 / 1 / -1 )
|
||||
B << Eigen::VectorXd::Zero(n), Eigen::VectorXd::Ones(n), -Eigen::VectorXd::Ones(n), Eigen::VectorXd::Zero(d) ;
|
||||
|
||||
// Fill Phi matrix
|
||||
// TODO
|
||||
|
||||
// Fill P matrix
|
||||
//TODO
|
||||
|
||||
// Set final A matrix
|
||||
/* A = Phi | P
|
||||
* P' | 0 */
|
||||
|
||||
// TODO
|
||||
|
||||
cout << "size of pts_list : " << nn << endl ;
|
||||
cout << "End computation of matrices" << endl ;
|
||||
|
||||
return {{A,B},pts_list} ;
|
||||
}
|
||||
|
||||
pair<vector<float>&, vector<float>&> Hole_Filling::solve_approx(const pair<Eigen::MatrixXd &, Eigen::VectorXd &> &p, int n, int d)
|
||||
{
|
||||
Eigen::MatrixXd & A = p.first ;
|
||||
Eigen::VectorXd & B = p.second ;
|
||||
|
||||
Eigen::VectorXd res = A.householderQr().solve(B) ;
|
||||
|
||||
cout << "End of solver" << endl ;
|
||||
cout << "res : " << res.head(10) << endl ;
|
||||
vector<float> & alpha = *(new vector<float>);
|
||||
vector<float> & beta = *(new vector<float>);
|
||||
|
||||
if (res.size() != (n+d))
|
||||
{
|
||||
cout << "taille du res : " << res.size() << endl ;
|
||||
throw std::runtime_error("Error in solve_approx") ;
|
||||
}
|
||||
for (int i=0; i<n; i++)
|
||||
{
|
||||
alpha.push_back(res(i)) ;
|
||||
}
|
||||
|
||||
for (int j=0; j<d; j++)
|
||||
{
|
||||
beta.push_back(res(n+j)) ;
|
||||
}
|
||||
|
||||
return {alpha, beta} ;
|
||||
}
|
||||
|
||||
// ***** Hole filling
|
||||
|
||||
void Hole_Filling::fill_hole(string out)
|
||||
{
|
||||
// TODO !!!
|
||||
cout << "Au travail !" << endl ;
|
||||
}
|
||||
|
||||
// ***** IO
|
||||
|
||||
void Hole_Filling::colorize_prop()
|
||||
{
|
||||
for (MyMesh::VertexIter v_it = _mesh.vertices_begin(); v_it != _mesh.vertices_end() ; ++v_it)
|
||||
{
|
||||
if(_mesh.property(_vprop, *v_it) == true)
|
||||
_mesh.set_color(*v_it, MyMesh::Color(255, 0, 0)) ;
|
||||
else
|
||||
_mesh.set_color(*v_it, MyMesh::Color(200, 200, 200)) ;
|
||||
}
|
||||
}
|
||||
|
||||
void Hole_Filling::colorize_verts(const vector<MyMesh::VertexHandle> &vlist)
|
||||
{
|
||||
for (MyMesh::VertexIter v_it = _mesh.vertices_begin(); v_it != _mesh.vertices_end() ; ++v_it)
|
||||
{
|
||||
_mesh.set_color(*v_it, MyMesh::Color(200, 200, 200)) ;
|
||||
}
|
||||
|
||||
for(int i=0; i<vlist.size(); i++)
|
||||
{
|
||||
_mesh.set_color(vlist.at(i), MyMesh::Color(255, 0, 0)) ;
|
||||
}
|
||||
}
|
||||
|
||||
Rect3 Hole_Filling::estimate_BB(const vector<MyMesh::VertexHandle> &vlist)
|
||||
{
|
||||
MyMesh::Point minp = _mesh.point(vlist.at(0)), maxp=minp ;
|
||||
for (int i=0; i<vlist.size(); i++)
|
||||
{
|
||||
minp = min(minp, _mesh.point(vlist.at(i))) ;
|
||||
maxp = max(maxp, _mesh.point(vlist.at(i))) ;
|
||||
}
|
||||
MyMesh::Point sizep(maxp-minp), centerp = (minp+maxp)/2 ;
|
||||
minp = centerp - _scale*sizep/2 ;
|
||||
sizep *= _scale ;
|
||||
|
||||
Rect3 domain ;
|
||||
domain.min = {minp[0], minp[1], minp[2]} ;
|
||||
domain.size = {sizep[0], sizep[1], sizep[2]} ;
|
||||
return domain ;
|
||||
}
|
||||
|
||||
void Hole_Filling::poly_n_out(const Implicit_RBF &implicit, Rect3 domain, string filename)
|
||||
{
|
||||
auto implicitEq = [&implicit](Vec3 const& pos)
|
||||
{
|
||||
MyMesh::Point X(pos.x, pos.y, pos.z) ;
|
||||
return implicit.val(X) ;
|
||||
};
|
||||
|
||||
|
||||
Vec3 cubeSize(domain.size*(_discr)) ;
|
||||
cout << "mind " << domain.min << endl ;
|
||||
cout << "sized" << domain.size << endl ;
|
||||
cout << "cubesize "<< cubeSize << endl ;
|
||||
|
||||
auto mesh = MarchCube(implicitEq, domain, cubeSize);
|
||||
WriteObjFile(mesh, filename);
|
||||
}
|
||||
|
||||
|
||||
void fillHolesImplicit(MyMesh &mesh) {
|
||||
}
|
||||
|
|
|
|||
|
|
@ -3,11 +3,125 @@
|
|||
|
||||
#include "my_mesh.h"
|
||||
#include <vector>
|
||||
#include <iostream>
|
||||
#include <cmath>
|
||||
#include <Eigen/Core>
|
||||
#include <Eigen/Geometry>
|
||||
#include <Eigen/Dense>
|
||||
#include <MeshReconstruction.h>
|
||||
#include <IO.h>
|
||||
|
||||
|
||||
void fillHoleDumb(MyMesh &mesh, std::vector<HalfedgeHandle> &hole);
|
||||
void fillHolesDumb(MyMesh &mesh);
|
||||
void fillHolesImplicit(MyMesh &mesh);
|
||||
|
||||
using namespace std ;
|
||||
using namespace MeshReconstruction;
|
||||
|
||||
// ******************************
|
||||
// Function computing polynomial P_i (degree 2 polynomials in 3 variables)
|
||||
inline float myp(int i, MyMesh::Point X) {
|
||||
float x = X[0] ;
|
||||
float y = X[1] ;
|
||||
float z = X[2] ;
|
||||
switch (i) {
|
||||
case 0:
|
||||
return x*x ;
|
||||
break;
|
||||
case 1:
|
||||
return y*y ;
|
||||
break;
|
||||
case 2:
|
||||
return z*z ;
|
||||
break;
|
||||
case 3:
|
||||
return x*y ;
|
||||
break;
|
||||
case 4:
|
||||
return y*z ;
|
||||
break;
|
||||
case 5:
|
||||
return x*z ;
|
||||
break;
|
||||
case 6:
|
||||
return x ;
|
||||
break;
|
||||
case 7:
|
||||
return y ;
|
||||
break;
|
||||
case 8:
|
||||
return z;
|
||||
break;
|
||||
case 9:
|
||||
return 1;
|
||||
break;
|
||||
default:
|
||||
throw std::runtime_error("Error on indice i : unknown polynomial P.");
|
||||
}
|
||||
}
|
||||
|
||||
// ******************************
|
||||
// Function computing thin spline RBF
|
||||
inline float myphi(float r) { if (r == 0) return 0 ; else return r*r*log(r) ; }
|
||||
|
||||
// ******************************
|
||||
// Class encoding an implicit function built from the basis of RBF
|
||||
// f(X) = Sum_i=0^n-1 alpha_i * phi(|| X-c_i ||) + Sum_j=0^9 beta_j * P_j(X)
|
||||
class Implicit_RBF {
|
||||
private:
|
||||
int _n ;
|
||||
int _d ;
|
||||
vector<float> _alpha, _beta ;
|
||||
vector<MyMesh::Point> _center ;
|
||||
|
||||
public:
|
||||
Implicit_RBF (vector<float> alpha, vector<float> beta, vector<MyMesh::Point> center) ;
|
||||
float val(MyMesh::Point) const ;
|
||||
};
|
||||
|
||||
// ******************************
|
||||
// Class encoding all the functions required for implicit hole filling
|
||||
class Hole_Filling
|
||||
{
|
||||
private:
|
||||
MyMesh &_mesh ;
|
||||
OpenMesh::VPropHandleT<bool> _vprop ;
|
||||
|
||||
// Constantes servant au rendu
|
||||
// TODO : à mettre dans des boutons ...
|
||||
const float _scale = 4 ; // Boîte englobante de rendu = _scale * BB des centres
|
||||
const float _discr = 1./40. ; // BB discrétisée en 1/_discr voxels
|
||||
public:
|
||||
Hole_Filling(MyMesh &mesh);
|
||||
inline ~Hole_Filling() { std::cout << "Ending hole filling" << std::endl ; }
|
||||
|
||||
// Computation of boundary and its neighborhood
|
||||
MyMesh::HalfedgeHandle find_boundary_edge() ;
|
||||
vector<MyMesh::VertexHandle> find_boundary(MyMesh::HalfedgeHandle heh) ;
|
||||
void init_mark_boundary(const vector<MyMesh::VertexHandle> & bnd) ;
|
||||
vector<MyMesh::VertexHandle> next_neighbors(const vector<MyMesh::VertexHandle> & bnd) ;
|
||||
|
||||
// Computation of RBF
|
||||
pair<pair<Eigen::MatrixXd &,Eigen::VectorXd &>,vector<MyMesh::Point> &> compute_approx_mat(vector<MyMesh::VertexHandle> vlist) ;
|
||||
pair<vector<float>&, vector<float>&> solve_approx(const pair<Eigen::MatrixXd &, Eigen::VectorXd &> &p, int n, int d) ;
|
||||
|
||||
// Hole filling
|
||||
void fill_hole(string out) ;
|
||||
|
||||
// IO
|
||||
void colorize_prop() ;
|
||||
void colorize_verts(const vector<MyMesh::VertexHandle> &vlist) ;
|
||||
Rect3 estimate_BB(const vector<MyMesh::VertexHandle> &vlist) ;
|
||||
void poly_n_out (const Implicit_RBF & implicit, Rect3 domain, std::string filename) ;
|
||||
};
|
||||
|
||||
// Other ....
|
||||
|
||||
inline MyMesh::Point min (MyMesh::Point X, MyMesh::Point Y) { MyMesh::Point P ; for (int i=0; i<3; i++) P[i] = min(X[i], Y[i]) ; return P ;}
|
||||
inline MyMesh::Point max (MyMesh::Point X, MyMesh::Point Y) { MyMesh::Point P ; for (int i=0; i<3; i++) P[i] = max(X[i], Y[i]) ; return P ;}
|
||||
|
||||
inline ostream & operator<< (ostream & out, Vec3 v) { out << v.x << ", " << v.y << ", " << v.z ; return out ;}
|
||||
|
||||
|
||||
#endif
|
||||
Loading…
Reference in New Issue