Compare commits
2 Commits
37fbe47aa8
...
d01ba47b65
Author | SHA1 | Date |
---|---|---|
papush! | d01ba47b65 | |
papush! | d51945ecaf |
|
@ -22,7 +22,11 @@ target_sources(${PROJECT_NAME} PRIVATE
|
|||
src/smoothing.cpp
|
||||
src/smoothing.h
|
||||
src/util.cpp
|
||||
src/util.h)
|
||||
src/util.h
|
||||
src/curvature.cpp
|
||||
src/curvature.h
|
||||
src/quad_patch.cpp
|
||||
src/quad_patch.h)
|
||||
target_link_libraries(${PROJECT_NAME} PRIVATE
|
||||
Qt5::Core Qt5::Gui Qt5::Widgets
|
||||
OpenMeshCore
|
||||
|
|
File diff suppressed because it is too large
Load Diff
|
@ -8128,6 +8128,7 @@ vn 0.8166 -0.5571 -0.1512
|
|||
vn 0.4879 -0.8093 -0.3271
|
||||
vn 0.4835 -0.8315 -0.2736
|
||||
vn 0.7111 -0.7029 -0.0163
|
||||
usemtl None
|
||||
s 1
|
||||
f 1//1 2//1 3//1
|
||||
f 4//2 5//2 6//2
|
||||
|
|
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
@ -1,30 +1,37 @@
|
|||
# Blender v2.90.1 OBJ File: ''
|
||||
# www.blender.org
|
||||
o Cube
|
||||
v -1.000000 -1.000000 1.000000
|
||||
v -1.000000 1.000000 1.000000
|
||||
v -1.000000 -1.000000 -1.000000
|
||||
v -1.000000 1.000000 -1.000000
|
||||
v 1.000000 -1.000000 1.000000
|
||||
v 1.000000 1.000000 1.000000
|
||||
v 1.000000 -1.000000 -1.000000
|
||||
v 1.000000 1.000000 -1.000000
|
||||
vn -1.0000 0.0000 0.0000
|
||||
vn 0.0000 0.0000 -1.0000
|
||||
vn 1.0000 0.0000 0.0000
|
||||
vn 0.0000 0.0000 1.0000
|
||||
vn 0.0000 -1.0000 0.0000
|
||||
v 1.000000 -1.000000 -1.000000
|
||||
v 1.000000 1.000000 1.000000
|
||||
v 1.000000 -1.000000 1.000000
|
||||
v -1.000000 1.000000 -1.000000
|
||||
v -1.000000 -1.000000 -1.000000
|
||||
v -1.000000 1.000000 1.000000
|
||||
v -1.000000 -1.000000 1.000000
|
||||
vt 0.625000 0.500000
|
||||
vt 0.875000 0.500000
|
||||
vt 0.875000 0.750000
|
||||
vt 0.625000 0.750000
|
||||
vt 0.375000 0.750000
|
||||
vt 0.625000 1.000000
|
||||
vt 0.375000 1.000000
|
||||
vt 0.375000 0.000000
|
||||
vt 0.625000 0.000000
|
||||
vt 0.625000 0.250000
|
||||
vt 0.375000 0.250000
|
||||
vt 0.125000 0.500000
|
||||
vt 0.375000 0.500000
|
||||
vt 0.125000 0.750000
|
||||
vn 0.0000 1.0000 0.0000
|
||||
vn 0.0000 0.0000 1.0000
|
||||
vn -1.0000 0.0000 0.0000
|
||||
vn 0.0000 -1.0000 0.0000
|
||||
vn 1.0000 0.0000 0.0000
|
||||
vn 0.0000 0.0000 -1.0000
|
||||
usemtl Material
|
||||
s off
|
||||
f 2//1 3//1 1//1
|
||||
f 4//2 7//2 3//2
|
||||
f 8//3 5//3 7//3
|
||||
f 6//4 1//4 5//4
|
||||
f 7//5 1//5 3//5
|
||||
f 4//6 6//6 8//6
|
||||
f 2//1 4//1 3//1
|
||||
f 4//2 8//2 7//2
|
||||
f 8//3 6//3 5//3
|
||||
f 6//4 2//4 1//4
|
||||
f 7//5 5//5 1//5
|
||||
f 4//6 2//6 6//6
|
||||
f 1/1/1 5/2/1 7/3/1 3/4/1
|
||||
f 4/5/2 3/4/2 7/6/2 8/7/2
|
||||
f 8/8/3 7/9/3 5/10/3 6/11/3
|
||||
f 6/12/4 2/13/4 4/5/4 8/14/4
|
||||
f 2/13/5 1/1/5 3/4/5 4/5/5
|
||||
f 6/11/6 5/10/6 1/1/6 2/13/6
|
||||
|
|
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
File diff suppressed because it is too large
Load Diff
|
@ -0,0 +1,265 @@
|
|||
#include "curvature.h"
|
||||
#include "util.h"
|
||||
|
||||
|
||||
void Courbures::normales_locales() {
|
||||
_mesh.request_vertex_normals();
|
||||
size_t i;
|
||||
for (VertexHandle vh : _mesh.vertices()) {
|
||||
MyMesh::Normal normal(0,0,0);
|
||||
i = 0;
|
||||
for (MyMesh::FaceHandle vf : _mesh.vf_range(vh)) {
|
||||
i++ ;
|
||||
normal += _mesh.calc_face_normal(vf);
|
||||
}
|
||||
if (i != 0) {
|
||||
normal /= i;
|
||||
normal /= norm(normal);
|
||||
}
|
||||
_mesh.set_normal(vh, normal);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
std::vector<MyMesh::VertexHandle> Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh)
|
||||
{
|
||||
OpenMesh::VPropHandleT<bool> vprop_flag;
|
||||
_mesh.add_property(vprop_flag, "vprop_flag");
|
||||
|
||||
// Initialisation
|
||||
for (MyMesh::VertexIter v_it = _mesh.vertices_begin(); v_it != _mesh.vertices_end(); ++v_it)
|
||||
{
|
||||
_mesh.property(vprop_flag, *v_it) = false ;
|
||||
}
|
||||
// Circulateur sur le premier cercle
|
||||
std::vector<MyMesh::VertexHandle> neigh, neigh2 ;
|
||||
|
||||
_mesh.property(vprop_flag, vh) = true ;
|
||||
for(MyMesh::VertexVertexIter vv_it = _mesh.vv_iter(vh);vv_it.is_valid();++vv_it)
|
||||
{
|
||||
neigh.push_back(*vv_it) ; // ajout du point à la liste
|
||||
_mesh.property(vprop_flag, *vv_it) = true ;
|
||||
}
|
||||
|
||||
// Parcours du premier cercle et ajout du second cercle par circulateurs
|
||||
for (int i=0; i<neigh.size(); i++)
|
||||
{
|
||||
MyMesh::VertexHandle vh = neigh.at(i) ;
|
||||
for(MyMesh::VertexVertexIter vv_it = _mesh.vv_iter(vh);vv_it.is_valid();++vv_it)
|
||||
{
|
||||
if (!_mesh.property(vprop_flag, *vv_it)) // sommet non encore rencontré
|
||||
neigh2.push_back(*vv_it) ; // ajout du point à la liste
|
||||
_mesh.property(vprop_flag, *vv_it) = true ;
|
||||
}
|
||||
}
|
||||
|
||||
// Concaténation des deux cercles
|
||||
neigh.insert(neigh.end(), neigh2.begin(), neigh2.end()) ;
|
||||
|
||||
_mesh.remove_property(vprop_flag);
|
||||
|
||||
return neigh ;
|
||||
}
|
||||
// std::vector<MyMesh::VertexHandle> Courbures::get_two_neighborhood(const MyMesh::VertexHandle vh) {
|
||||
// OpenMesh::VPropHandleT<bool> vprop_flag;
|
||||
// _mesh.add_property(vprop_flag, "vprop_flag");
|
||||
|
||||
// // Initialisation
|
||||
// for (VertexHandle vh : _mesh.vertices()) {
|
||||
// _mesh.property(vprop_flag, vh) = false;
|
||||
// }
|
||||
// // Circulateur sur le premier cercle
|
||||
// std::vector<MyMesh::VertexHandle> neigh, neigh2;
|
||||
|
||||
// _mesh.property(vprop_flag, vh) = true;
|
||||
// for(VertexHandle vv : _mesh.vv_range(vh)) {
|
||||
// neigh.push_back(vv); // ajout du point à la liste
|
||||
// _mesh.property(vprop_flag, vv) = true;
|
||||
// }
|
||||
|
||||
// // Parcours du premier cercle et ajout du second cercle par circulateurs
|
||||
// for (size_t i = 0; i < neigh.size(); i++) {
|
||||
// MyMesh::VertexHandle vh = neigh.at(i) ;
|
||||
// for (VertexHandle vv : _mesh.vv_range(vh)) {
|
||||
// if (!_mesh.property(vprop_flag, vv)) {
|
||||
// neigh2.push_back(vv);
|
||||
// }
|
||||
// _mesh.property(vprop_flag, vv) = true;
|
||||
// }
|
||||
// }
|
||||
|
||||
// // Concaténation des deux cercles
|
||||
// neigh.insert(neigh.end(), neigh2.begin(), neigh2.end());
|
||||
// _mesh.remove_property(vprop_flag);
|
||||
// return neigh;
|
||||
// }
|
||||
|
||||
|
||||
MyQuad Courbures::fit_quad(MyMesh::VertexHandle vh)
|
||||
{
|
||||
MyQuad q ;
|
||||
|
||||
std::vector<MyMesh::VertexHandle> neigh = get_two_neighborhood(vh) ;
|
||||
|
||||
// std::cout << "-- Neighborhood" << std::endl ;
|
||||
for (int i = 0; i< neigh.size(); i++)
|
||||
{
|
||||
MyMesh::Point p = _mesh.point(neigh.at(i)) ;
|
||||
// std::cout << p[0] << ", " << p[1] << ", " << p[2] << " ; " << std::endl ;
|
||||
}
|
||||
|
||||
// Calcul de la matrice de changement de base
|
||||
MyMesh::Normal n = _mesh.normal(vh) ;
|
||||
MyMesh::Point p = _mesh.point(vh);
|
||||
// std::cout << "-- normale" << std::endl ;
|
||||
// std::cout << n[0] << ", " << n[1] << ", " << n[2] << std::endl ;
|
||||
// std::cout << "-- point" << std::endl ;
|
||||
// std::cout << p[0] << ", " << p[1] << ", " << p[2] << std::endl ;
|
||||
|
||||
Eigen::Vector3d ne(n[0], n[1], n[2]), Oz(0,0,1), axis;
|
||||
Eigen::Vector3d p_e(p[0], p[1], p[2]), pi_e ;
|
||||
|
||||
axis = ne.cross(Oz) ;
|
||||
double sina = axis.norm(), cosa = ne.dot(Oz), angle ;
|
||||
if (sina >= 0)
|
||||
angle = acos(cosa) ;
|
||||
else
|
||||
angle = -acos(cosa) ;
|
||||
axis = axis.normalized() ;
|
||||
|
||||
Eigen::AngleAxisd r(angle, axis) ;
|
||||
// std::cout << "-- rotation" << std:: endl ;
|
||||
// std::cout << r.matrix()(0,0) << ", " << r.matrix()(0, 1) << ", " << r.matrix()(0,2) << "; " ;
|
||||
// std::cout << r.matrix()(1,0) << ", " << r.matrix()(1, 1) << ", " << r.matrix()(1,2) << "; " ;
|
||||
// std::cout << r.matrix()(2,0) << ", " << r.matrix()(2, 1) << ", " << r.matrix()(2,2) << std::endl ;
|
||||
Eigen::Translation3d t(-p_e[0], -p_e[1], -p_e[2]) ;
|
||||
Eigen::Transform<double, 3, Eigen::Affine> ch_base = r * t ;
|
||||
|
||||
|
||||
// Calcul de la matrice / vecteur de moindres carrés linéaires (Eigen)
|
||||
|
||||
if (neigh.size() >= 5)
|
||||
{
|
||||
int n(neigh.size()) ;
|
||||
Eigen::MatrixXd A(n,5);
|
||||
Eigen::VectorXd B(n);
|
||||
|
||||
// std::cout << "-- Après changement de base" << std::endl ;
|
||||
for(int i=0; i<neigh.size(); i++)
|
||||
{
|
||||
MyMesh::Point p = _mesh.point(neigh.at(i)) ;
|
||||
pi_e << p[0], p[1], p[2] ;
|
||||
pi_e = ch_base * pi_e ; // Application du changement de base
|
||||
// std::cout << pi_e[0] << ", " << pi_e[1] << ", " << pi_e[2] << std::endl ;
|
||||
|
||||
B[i] = -pi_e[2]; // -zi
|
||||
A(i, 0) = pi_e[0] * pi_e[0]; // xi^2
|
||||
A(i, 1) = pi_e[0] * pi_e[1]; // xi*yi
|
||||
A(i, 2) = pi_e[1] * pi_e[1]; // yi^2
|
||||
A(i, 3) = pi_e[0]; // xi
|
||||
A(i, 4) = pi_e[1]; // yi
|
||||
}
|
||||
|
||||
// Résolution aux moindres carrés par SVD
|
||||
Eigen::VectorXd coef(5) ;
|
||||
coef = A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B) ;
|
||||
MyQuad q(coef[0], coef[1], coef[2], coef[3], coef[4], r, t) ;
|
||||
// std::cout << "-- quad : " << std::endl ;
|
||||
// std::cout << q[0] << ", " << q[1] << ", " << q[2] << ", " << q[3] << ", " << q[4] << ", " << q[5] << ", " << std::endl ;
|
||||
return q ;
|
||||
}
|
||||
else
|
||||
{
|
||||
std::cout << "Quad fitting : not enough neighbors" ;
|
||||
throw "Quad fitting : not enough neighbors" ;
|
||||
}
|
||||
}
|
||||
// MyQuad Courbures::fit_quad(MyMesh::VertexHandle vh) {
|
||||
// std::vector<MyMesh::VertexHandle> neigh = get_two_neighborhood(vh);
|
||||
// if (neigh.size() < 5) throw "Quad fitting: not enough neighbors";
|
||||
|
||||
// // Calcul de la matrice de changement de base
|
||||
// Eigen::Vector3d Oz(0,0,1);
|
||||
// Eigen::Vector3d axis = _mesh.normal(vh).cross(Oz);
|
||||
// double sina = axis.norm();
|
||||
// double cosa = _mesh.normal(vh).dot(Oz);
|
||||
// double angle;
|
||||
// if (sina >= 0) {
|
||||
// angle = acos(cosa);
|
||||
// } else {
|
||||
// angle = -acos(cosa);
|
||||
// }
|
||||
// axis = axis.normalized();
|
||||
// Eigen::AngleAxisd rot(angle, axis);
|
||||
// Eigen::Translation3d trans(-_mesh.point(vh));
|
||||
// Eigen::Transform<double, 3, Eigen::Affine> ch_base = rot * trans;
|
||||
|
||||
// // Calcul de la matrice / vecteur de moindres carrés linéaires (Eigen)
|
||||
// Eigen::MatrixXd A(neigh.size(), 5);
|
||||
// Eigen::VectorXd B(neigh.size());
|
||||
// for(size_t i = 0; i < neigh.size(); i++) {
|
||||
// MyMesh::Point point = _mesh.point(neigh[i]);
|
||||
// point = ch_base * point; // Application du changement de base
|
||||
// B[i] = -point[2]; // -zi
|
||||
// A(i, 0) = point[0] * point[0]; // xi^2
|
||||
// A(i, 1) = point[0] * point[1]; // xi*yi
|
||||
// A(i, 2) = point[1] * point[1]; // yi^2
|
||||
// A(i, 3) = point[0]; // xi
|
||||
// A(i, 4) = point[1]; // yi
|
||||
// }
|
||||
|
||||
// // Résolution aux moindres carrés par SVD
|
||||
// Eigen::VectorXd coef(5);
|
||||
// coef = A.bdcSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(B);
|
||||
// return MyQuad(coef[0], coef[1], coef[2], coef[3], coef[4], rot, trans);
|
||||
// }
|
||||
|
||||
|
||||
void Courbures::compute_KH() {
|
||||
_mesh.add_property(vprop_K, "vprop_K");
|
||||
_mesh.add_property(vprop_H, "vprop_H");
|
||||
_mesh.add_property(vprop_quad, "vprop_quad");
|
||||
|
||||
for (VertexHandle vh : _mesh.vertices()) {
|
||||
MyQuad q = fit_quad(vh);
|
||||
_mesh.property(vprop_quad, vh) = q;
|
||||
|
||||
// K = det (K_P) = 4 a_0 a_2 - a_1 ^ 2
|
||||
_mesh.property(vprop_K, vh) = 4 * q[0] * q[2] - q[1]*q[1];
|
||||
|
||||
// H = Tr (K_P) = a_0 + a_2
|
||||
_mesh.property(vprop_H, vh) = q[0] + q[2];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Courbures::set_fixed_colors() {
|
||||
size_t i = 0;
|
||||
for (VertexHandle vh : _mesh.vertices()) {
|
||||
if (i % 2)
|
||||
_mesh.set_color(vh, MyMesh::Color(0, 1, 0));
|
||||
else
|
||||
_mesh.set_color(vh, MyMesh::Color(0, 0, 1));
|
||||
i++;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Courbures::set_K_colors() {
|
||||
MyStats<double> dist;
|
||||
for (VertexHandle vh : _mesh.vertices()) {
|
||||
dist.push_back(_mesh.property(vprop_K, vh));
|
||||
}
|
||||
double m = dist.mean();
|
||||
double sigma = dist.stdev(m);
|
||||
std::cout << "K min : " << dist.min()
|
||||
<< " - max : " << dist.max() << std::endl ;
|
||||
std::cout << "K mean : " << m << " - sigma : " << sigma << std::endl;
|
||||
|
||||
for (VertexHandle vh : _mesh.vertices()) {
|
||||
double tmp = _mesh.property(vprop_K, vh) - m;
|
||||
tmp = clamp<double>(tmp, -sigma, +sigma);
|
||||
tmp = (tmp + sigma) / (2 * sigma);
|
||||
_mesh.set_color(vh, MyMesh::Color(tmp, .23, .5));
|
||||
}
|
||||
}
|
|
@ -0,0 +1,163 @@
|
|||
#ifndef CURVATURE_H
|
||||
#define CURVATURE_H
|
||||
|
||||
#include "my_mesh.h"
|
||||
|
||||
|
||||
class MyQuad {
|
||||
double _coefs[5] ; // a_0 x^2 + a1 xy + a2 y^2 + a3 x + a4 y + a5
|
||||
public:
|
||||
Eigen::AngleAxisd _r ;
|
||||
Eigen::Translation3d _t ;
|
||||
|
||||
MyQuad(double *data,
|
||||
const Eigen::AngleAxisd &r = Eigen::AngleAxisd(0, Eigen::Vector3d(0,0,1)),
|
||||
const Eigen::Translation3d &t = Eigen::Translation3d(0,0,0))
|
||||
: _r(r), _t(t)
|
||||
{
|
||||
for (size_t i=0; i<5; i++)
|
||||
_coefs[i] = data[i] ;
|
||||
}
|
||||
MyQuad(const Eigen::VectorXd &v,
|
||||
const Eigen::AngleAxisd &r = Eigen::AngleAxisd(0, Eigen::Vector3d(0,0,1)),
|
||||
const Eigen::Translation3d &t = Eigen::Translation3d(0,0,0))
|
||||
: _r(r), _t(t)
|
||||
{
|
||||
for (size_t i=0; i<5; i++)
|
||||
_coefs[i] = v[i] ;
|
||||
}
|
||||
MyQuad(double a0=0, double a1=0, double a2=0, double a3=0, double a4=0,
|
||||
const Eigen::AngleAxisd &r = Eigen::AngleAxisd(0, Eigen::Vector3d(0,0,1)),
|
||||
const Eigen::Translation3d &t = Eigen::Translation3d(0,0,0))
|
||||
: _r(r), _t(t)
|
||||
{
|
||||
_coefs[0] = a0 ;
|
||||
_coefs[1] = a1 ;
|
||||
_coefs[2] = a2 ;
|
||||
_coefs[3] = a3 ;
|
||||
_coefs[4] = a4 ;
|
||||
}
|
||||
MyQuad(const MyQuad &q) : _r(q._r), _t(q._t)
|
||||
{
|
||||
for (size_t i=0; i<5; i++)
|
||||
_coefs[i] = q._coefs[i] ;
|
||||
}
|
||||
|
||||
MyQuad &operator=(const MyQuad &q) {
|
||||
_r = q._r;
|
||||
_t = q._t;
|
||||
for (size_t i=0; i<5; i++)
|
||||
_coefs[i] = q._coefs[i] ;
|
||||
return *this;
|
||||
}
|
||||
|
||||
double & operator[] (size_t i) {return _coefs[i] ; }
|
||||
|
||||
double quad_fun (double x, double y)
|
||||
{
|
||||
return _coefs[0] * x*x + _coefs[1] * x*y + _coefs[2] * y*y + _coefs[3] * x + _coefs[4] * y ;
|
||||
}
|
||||
|
||||
double quad_fun (const OpenMesh::Vec2d &v)
|
||||
{
|
||||
|
||||
return quad_fun(v[0], v[1]) ;
|
||||
}
|
||||
|
||||
};
|
||||
|
||||
|
||||
template <typename T>
|
||||
class MyStats {
|
||||
private:
|
||||
std::vector<T> _distrib ;
|
||||
public:
|
||||
MyStats () {} ;
|
||||
|
||||
void push_back (T data)
|
||||
{
|
||||
_distrib.push_back(data) ;
|
||||
}
|
||||
|
||||
T min ()
|
||||
{
|
||||
T tmp = _distrib.at(0) ;
|
||||
for (size_t i=1 ; i<_distrib.size(); i++)
|
||||
{
|
||||
if (_distrib.at(i) < tmp)
|
||||
tmp = _distrib.at(i) ;
|
||||
}
|
||||
return tmp ;
|
||||
}
|
||||
|
||||
T max ()
|
||||
{
|
||||
T tmp = _distrib.at(0) ;
|
||||
for (size_t i=1 ; i<_distrib.size(); i++)
|
||||
{
|
||||
if (_distrib.at(i) > tmp)
|
||||
tmp = _distrib.at(i) ;
|
||||
}
|
||||
return tmp ;
|
||||
}
|
||||
|
||||
T mean ()
|
||||
{
|
||||
T acc(0) ;
|
||||
std::cout << "acc : " << acc << std::endl;
|
||||
if (_distrib.size() > 0)
|
||||
{
|
||||
for(size_t i=0; i<_distrib.size(); i++)
|
||||
{
|
||||
acc += _distrib.at(i) ;
|
||||
}
|
||||
return acc/(_distrib.size()) ;
|
||||
}
|
||||
else
|
||||
return acc ;
|
||||
}
|
||||
|
||||
T stdev ()
|
||||
{
|
||||
return stdev(mean());
|
||||
}
|
||||
|
||||
T stdev (T m)
|
||||
{
|
||||
T acc(0), tmp ;
|
||||
if (_distrib.size() > 0)
|
||||
{
|
||||
for(size_t i=0; i<_distrib.size(); i++)
|
||||
{
|
||||
tmp = _distrib.at(i) - m ;
|
||||
acc += tmp * tmp ;
|
||||
}
|
||||
return sqrt(acc/(_distrib.size())) ;
|
||||
}
|
||||
else
|
||||
return acc ;
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
class Courbures {
|
||||
private:
|
||||
MyMesh &_mesh ;
|
||||
|
||||
public:
|
||||
OpenMesh::VPropHandleT<double> vprop_K;
|
||||
OpenMesh::VPropHandleT<double> vprop_H;
|
||||
OpenMesh::VPropHandleT<MyQuad> vprop_quad;
|
||||
|
||||
Courbures(MyMesh &mesh) : _mesh(mesh) {}
|
||||
|
||||
void set_fixed_colors() ;
|
||||
void normales_locales() ;
|
||||
std::vector<MyMesh::VertexHandle> get_two_neighborhood(MyMesh::VertexHandle vh);
|
||||
MyQuad fit_quad(MyMesh::VertexHandle vh) ;
|
||||
void compute_KH() ;
|
||||
void set_K_colors() ;
|
||||
};
|
||||
|
||||
|
||||
#endif
|
|
@ -21,7 +21,8 @@ int main(int argc, char *argv[]) {
|
|||
QObject::connect(mesh_viewer, &MeshViewer::initialized,
|
||||
[&]() {
|
||||
if (mesh_processor) {
|
||||
mesh_viewer->addMesh(*mesh_processor);
|
||||
mesh_viewer->addMesh(mesh_processor->mesh);
|
||||
mesh_viewer->addMesh(mesh_processor->patch);
|
||||
}
|
||||
});
|
||||
|
||||
|
@ -34,11 +35,11 @@ int main(int argc, char *argv[]) {
|
|||
QObject::connect(&main_window, &MainWindow::open,
|
||||
[&](const QString &path) {
|
||||
if (mesh_processor) {
|
||||
mesh_viewer->removeMesh(*mesh_processor);
|
||||
mesh_viewer->removeMesh(mesh_processor->mesh);
|
||||
delete mesh_processor;
|
||||
}
|
||||
mesh_processor = new MeshProcessor(path);
|
||||
mesh_viewer->addMesh(*mesh_processor);
|
||||
mesh_viewer->addMesh(mesh_processor->mesh);
|
||||
});
|
||||
main_window.show();
|
||||
return app.exec();
|
||||
|
|
|
@ -1,6 +1,8 @@
|
|||
#include "quad_patch.h"
|
||||
#include "mesh_processor.h"
|
||||
#include "util.h"
|
||||
#include "smoothing.h"
|
||||
#include "curvature.h"
|
||||
#include <QGenericMatrix>
|
||||
#include <unordered_set>
|
||||
|
||||
|
@ -31,9 +33,16 @@ MeshProcessor::MeshProcessor(const QString &path) {
|
|||
qWarning() << "Failed to read" << path;
|
||||
return;
|
||||
}
|
||||
// holes = findHoles(mesh);
|
||||
mesh.update_normals();
|
||||
Courbures courb(mesh);
|
||||
courb.compute_KH();
|
||||
courb.set_K_colors();
|
||||
VertexHandle vh = mesh.vertex_handle(16);
|
||||
MyQuad q = mesh.property(courb.vprop_quad, vh);
|
||||
patch = quad_patch(q, mesh, vh);
|
||||
// mesh.holes = findHoles(mesh);
|
||||
// fillHoles();
|
||||
smooth(mesh);
|
||||
// smooth(mesh);
|
||||
}
|
||||
|
||||
|
||||
|
@ -60,7 +69,7 @@ void fillHoleDumb(MyMesh &mesh, std::vector<HalfedgeHandle> &hole) {
|
|||
}
|
||||
|
||||
void MeshProcessor::fillHoles() {
|
||||
for (auto hole : holes) {
|
||||
for (auto hole : mesh.holes) {
|
||||
fillHoleDumb(mesh, hole);
|
||||
}
|
||||
}
|
|
@ -12,7 +12,7 @@ class MeshProcessor : public QObject {
|
|||
|
||||
public:
|
||||
MyMesh mesh;
|
||||
std::vector<std::vector<HalfedgeHandle>> holes;
|
||||
MyMesh patch;
|
||||
|
||||
MeshProcessor(const QString &path);
|
||||
|
||||
|
|
|
@ -4,21 +4,22 @@
|
|||
#include <QOpenGLFunctions_2_1>
|
||||
|
||||
|
||||
MeshView::MeshView(const MeshProcessor &mesh_processor, QOpenGLShaderProgram &program)
|
||||
MeshView::MeshView(const MyMesh &mesh, QOpenGLShaderProgram &program)
|
||||
:vertex_buffer(QOpenGLBuffer::VertexBuffer),
|
||||
index_buffer(QOpenGLBuffer::IndexBuffer),
|
||||
holes_index_buffer(QOpenGLBuffer::IndexBuffer),
|
||||
n_faces(mesh_processor.mesh.n_faces()),
|
||||
n_vertices(mesh_processor.mesh.n_vertices()),
|
||||
mesh_processor(mesh_processor) {
|
||||
const MyMesh &mesh = mesh_processor.mesh;
|
||||
|
||||
QVector<GLfloat> vertices(n_vertices * 3);
|
||||
n_faces(mesh.n_faces()),
|
||||
n_vertices(mesh.n_vertices()),
|
||||
mesh(mesh) {
|
||||
QVector<GLfloat> vertices(n_vertices * 6);
|
||||
size_t i = 0;
|
||||
for (const VertexHandle it : mesh.vertices()) {
|
||||
vertices[3*i + 0] = mesh.point(it)[0];
|
||||
vertices[3*i + 1] = mesh.point(it)[1];
|
||||
vertices[3*i + 2] = mesh.point(it)[2];
|
||||
vertices[6*i + 0] = mesh.point(it)[0];
|
||||
vertices[6*i + 1] = mesh.point(it)[1];
|
||||
vertices[6*i + 2] = mesh.point(it)[2];
|
||||
vertices[6*i + 3] = mesh.color(it)[0];
|
||||
vertices[6*i + 4] = mesh.color(it)[1];
|
||||
vertices[6*i + 5] = mesh.color(it)[2];
|
||||
i++;
|
||||
}
|
||||
|
||||
|
@ -36,7 +37,7 @@ MeshView::MeshView(const MeshProcessor &mesh_processor, QOpenGLShaderProgram &pr
|
|||
|
||||
QVector<GLuint> holes_indices;
|
||||
n_boundary_halfedges = 0;
|
||||
for (std::vector<HalfedgeHandle> hole : mesh_processor.holes) {
|
||||
for (std::vector<HalfedgeHandle> hole : mesh.holes) {
|
||||
for (HalfedgeHandle it : hole) {
|
||||
holes_indices.push_back(mesh.from_vertex_handle(it).idx());
|
||||
holes_indices.push_back(mesh.to_vertex_handle(it).idx());
|
||||
|
@ -48,18 +49,24 @@ MeshView::MeshView(const MeshProcessor &mesh_processor, QOpenGLShaderProgram &pr
|
|||
{QOpenGLVertexArrayObject::Binder binder(&vao);
|
||||
vertex_buffer.create();
|
||||
vertex_buffer.bind();
|
||||
vertex_buffer.setUsagePattern(QOpenGLBuffer::StreamDraw);
|
||||
vertex_buffer.allocate(vertices.constData(), n_vertices * 3 * sizeof (GLfloat));
|
||||
vertex_buffer.setUsagePattern(QOpenGLBuffer::DynamicDraw);
|
||||
vertex_buffer.allocate(vertices.constData(), n_vertices * 6 * sizeof (GLfloat));
|
||||
index_buffer.create();
|
||||
index_buffer.bind();
|
||||
index_buffer.setUsagePattern(QOpenGLBuffer::StreamDraw);
|
||||
index_buffer.setUsagePattern(QOpenGLBuffer::DynamicDraw);
|
||||
index_buffer.allocate(indices.constData(), n_faces * 3 * sizeof (GLuint));
|
||||
holes_index_buffer.create();
|
||||
holes_index_buffer.bind();
|
||||
holes_index_buffer.setUsagePattern(QOpenGLBuffer::StreamDraw);
|
||||
holes_index_buffer.setUsagePattern(QOpenGLBuffer::DynamicDraw);
|
||||
holes_index_buffer.allocate(holes_indices.constData(), n_boundary_halfedges * 2 * sizeof (GLuint));
|
||||
program.setAttributeBuffer("pos", GL_FLOAT, 0, 3);
|
||||
|
||||
vertex_buffer.bind();
|
||||
program.setAttributeBuffer("pos", GL_FLOAT, 0,
|
||||
3, 6 * sizeof (GLfloat));
|
||||
program.enableAttributeArray("pos");
|
||||
program.setAttributeBuffer("col", GL_FLOAT, 3 * sizeof (GLfloat),
|
||||
3, 6 * sizeof (GLfloat));
|
||||
program.enableAttributeArray("col");
|
||||
}
|
||||
}
|
||||
|
||||
|
@ -72,7 +79,6 @@ MeshView::~MeshView() {
|
|||
|
||||
|
||||
void MeshView::paint(QOpenGLShaderProgram &program) {
|
||||
const MyMesh &mesh = mesh_processor.mesh;
|
||||
QOpenGLContext *ctx = QOpenGLContext::currentContext();
|
||||
QOpenGLExtraFunctions *glf = ctx->extraFunctions();
|
||||
|
||||
|
@ -82,7 +88,7 @@ void MeshView::paint(QOpenGLShaderProgram &program) {
|
|||
qFatal("Failed to get OpenGL 2.1 functions");
|
||||
|
||||
program.setUniformValue("model", mesh.transform);
|
||||
program.setUniformValue("col", mesh.color.redF(), mesh.color.greenF(), mesh.color.blueF());
|
||||
// program.setUniformValue("col", mesh.color.redF(), mesh.color.greenF(), mesh.color.blueF());
|
||||
|
||||
{QOpenGLVertexArrayObject::Binder binder(&vao);
|
||||
/* Mesh */
|
||||
|
|
|
@ -1,7 +1,7 @@
|
|||
#ifndef MESH_VIEW_H
|
||||
#define MESH_VIEW_H
|
||||
|
||||
#include "mesh_processor.h"
|
||||
#include "my_mesh.h"
|
||||
#include <QOpenGLExtraFunctions>
|
||||
#include <QOpenGLShaderProgram>
|
||||
#include <QOpenGLVertexArrayObject>
|
||||
|
@ -18,9 +18,9 @@ class MeshView {
|
|||
unsigned n_boundary_halfedges;
|
||||
|
||||
public:
|
||||
const MeshProcessor &mesh_processor;
|
||||
const MyMesh &mesh;
|
||||
|
||||
MeshView(const MeshProcessor &mesh_processor, QOpenGLShaderProgram &program);
|
||||
MeshView(const MyMesh &mesh_processor, QOpenGLShaderProgram &program);
|
||||
~MeshView();
|
||||
void paint(QOpenGLShaderProgram &program);
|
||||
};
|
||||
|
|
|
@ -83,19 +83,19 @@ void MeshViewer::paintGL() {
|
|||
}
|
||||
|
||||
|
||||
void MeshViewer::addMesh(const MeshProcessor &mesh_processor) {
|
||||
void MeshViewer::addMesh(const MyMesh &mesh) {
|
||||
Q_ASSERT(isValid());
|
||||
makeCurrent();
|
||||
meshes.emplace_back(mesh_processor, program);
|
||||
meshes.emplace_back(mesh, program);
|
||||
doneCurrent();
|
||||
update();
|
||||
}
|
||||
|
||||
|
||||
void MeshViewer::removeMesh(const MeshProcessor &mesh_processor) {
|
||||
void MeshViewer::removeMesh(const MyMesh &mesh) {
|
||||
makeCurrent();
|
||||
meshes.remove_if([&](const MeshView &mv) {
|
||||
return &mv.mesh_processor == &mesh_processor;
|
||||
return &mv.mesh == &mesh;
|
||||
});
|
||||
doneCurrent();
|
||||
update();
|
||||
|
|
|
@ -2,7 +2,7 @@
|
|||
#define MESH_VIEWER_H
|
||||
|
||||
#define GL_GLEXT_PROTOTYPES
|
||||
#include "mesh_processor.h"
|
||||
#include "my_mesh.h"
|
||||
#include "mesh_view.h"
|
||||
#include <list>
|
||||
#include <QMatrix4x4>
|
||||
|
@ -39,8 +39,8 @@ public:
|
|||
void paintGL() override;
|
||||
|
||||
public slots:
|
||||
void addMesh(const MeshProcessor &mesh);
|
||||
void removeMesh(const MeshProcessor &mesh);
|
||||
void addMesh(const MyMesh &mesh);
|
||||
void removeMesh(const MyMesh &mesh);
|
||||
|
||||
protected:
|
||||
virtual void mousePressEvent(QMouseEvent *e);
|
||||
|
|
|
@ -10,11 +10,18 @@
|
|||
#include <OpenMesh/Core/Geometry/EigenVectorT.hh>
|
||||
|
||||
|
||||
// struct MyTraits : public OpenMesh::DefaultTraits {
|
||||
// using Point = Eigen::Vector3<qreal>;
|
||||
// using Normal = Eigen::Vector3<qreal>;
|
||||
// using Color = Eigen::Vector3<qreal>;
|
||||
// VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Color);
|
||||
// HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge);
|
||||
// FaceAttributes(OpenMesh::Attributes::Normal);
|
||||
// EdgeAttributes(OpenMesh::Attributes::Color);
|
||||
// };
|
||||
struct MyTraits : public OpenMesh::DefaultTraits {
|
||||
using Point = Eigen::Vector3<qreal>;
|
||||
using Normal = Eigen::Vector3<qreal>;
|
||||
using Color = Eigen::Vector3<qreal>;
|
||||
VertexAttributes(OpenMesh::Attributes::Normal);
|
||||
typedef OpenMesh::Vec3f Color;
|
||||
VertexAttributes(OpenMesh::Attributes::Normal | OpenMesh::Attributes::Color);
|
||||
HalfedgeAttributes(OpenMesh::Attributes::PrevHalfedge);
|
||||
FaceAttributes(OpenMesh::Attributes::Normal);
|
||||
EdgeAttributes(OpenMesh::Attributes::Color);
|
||||
|
@ -23,7 +30,8 @@ struct MyTraits : public OpenMesh::DefaultTraits {
|
|||
class MyMesh : public OpenMesh::TriMesh_ArrayKernelT<MyTraits> {
|
||||
public:
|
||||
QMatrix4x4 transform;
|
||||
QColor color = {127, 127, 127};
|
||||
std::vector<std::vector<HalfedgeHandle>> holes;
|
||||
// QColor color = {127, 127, 127};
|
||||
};
|
||||
|
||||
typedef MyMesh::FaceHandle FaceHandle;
|
||||
|
|
|
@ -0,0 +1,46 @@
|
|||
#include "quad_patch.h"
|
||||
|
||||
|
||||
MyMesh quad_patch(MyQuad q, MyMesh mesh, VertexHandle vh) {
|
||||
MyMesh patch;
|
||||
Eigen::Vector3d point
|
||||
(mesh.point(vh)[0], mesh.point(vh)[1], mesh.point(vh)[2]);
|
||||
Eigen::Transform<double, 3, Eigen::Affine> ch_base = q._r * q._t;
|
||||
Eigen::Transform<double, 3, Eigen::Affine> inv_ch_base = ch_base.inverse();
|
||||
point = ch_base * point;
|
||||
double xmin = point[0], xmax = point[0];
|
||||
double ymin = point[1], ymax = point[1];
|
||||
for (VertexHandle vi : mesh.vv_range(vh)) {
|
||||
Eigen::Vector3d point
|
||||
(mesh.point(vi)[0], mesh.point(vi)[1], mesh.point(vi)[2]);
|
||||
point = ch_base * point;
|
||||
xmin = std::min(xmin, point[0]);
|
||||
xmax = std::max(xmax, point[0]);
|
||||
ymin = std::min(ymin, point[1]);
|
||||
ymax = std::max(ymax, point[1]);
|
||||
}
|
||||
double xstep = (xmax-xmin)/64.;
|
||||
double ystep = (ymax-ymin)/64.;
|
||||
for (size_t y = 0; y < 64; y++) {
|
||||
for (size_t x = 0; x < 64; x++) {
|
||||
double dx = xmin + x * xstep;
|
||||
double dy = ymin + y * ystep;
|
||||
Eigen::Vector3d point(dx, dy, -q.quad_fun(dx, dy));
|
||||
point = inv_ch_base * point;
|
||||
patch.new_vertex(Point(point[0], point[1], point[2]));
|
||||
}
|
||||
}
|
||||
for (VertexHandle vhi : patch.vertices()) {
|
||||
patch.set_color(vhi, MyMesh::Color(0, 1, .2));
|
||||
size_t i = vhi.idx();
|
||||
size_t x = i % 64;
|
||||
size_t y = i / 64;
|
||||
if (x == 63 || y == 63) continue;
|
||||
patch.add_face(vhi, patch.vertex_handle(i + 64),
|
||||
patch.vertex_handle(i + 1));
|
||||
patch.add_face(patch.vertex_handle(i + 1),
|
||||
patch.vertex_handle(i + 64),
|
||||
patch.vertex_handle(i + 65));
|
||||
}
|
||||
return patch;
|
||||
}
|
|
@ -0,0 +1,11 @@
|
|||
#ifndef QUAD_PATCH_H
|
||||
#define QUAD_PATCH_H
|
||||
|
||||
#include "curvature.h"
|
||||
#include "my_mesh.h"
|
||||
|
||||
|
||||
MyMesh quad_patch(MyQuad q, MyMesh mesh, VertexHandle vh);
|
||||
|
||||
|
||||
#endif
|
|
@ -1,3 +1,5 @@
|
|||
varying vec3 frag_col;
|
||||
|
||||
uniform vec3 col;
|
||||
uniform vec3 wf_col;
|
||||
uniform bool wireframe;
|
||||
|
@ -7,5 +9,5 @@ void main() {
|
|||
if (wireframe)
|
||||
gl_FragColor = vec4(wf_col, alpha);
|
||||
else
|
||||
gl_FragColor = vec4(col, alpha);
|
||||
gl_FragColor = vec4(frag_col, alpha);
|
||||
}
|
||||
|
|
|
@ -1,9 +1,13 @@
|
|||
attribute vec3 pos;
|
||||
attribute vec3 col;
|
||||
|
||||
varying vec3 frag_col;
|
||||
|
||||
uniform mat4 proj;
|
||||
uniform mat4 view;
|
||||
uniform mat4 model;
|
||||
|
||||
void main() {
|
||||
frag_col = col;
|
||||
gl_Position = proj * view * model * vec4(pos, 1.0);
|
||||
}
|
||||
|
|
|
@ -28,6 +28,7 @@ Point laplace_beltrami(const MyMesh &mesh, VertexHandle vert) {
|
|||
Point sum {0, 0, 0};
|
||||
qreal area = 0;
|
||||
Point p = mesh.point(vert);
|
||||
qreal count = 0;
|
||||
for (HalfedgeHandle v_vi : mesh.voh_range(vert)) {
|
||||
Point pi = mesh.point(mesh.to_vertex_handle(v_vi));
|
||||
HalfedgeHandle vi_v = mesh.opposite_halfedge_handle(v_vi);
|
||||
|
@ -37,74 +38,81 @@ Point laplace_beltrami(const MyMesh &mesh, VertexHandle vert) {
|
|||
Point pb = mesh.point(mesh.to_vertex_handle(vi_vb));
|
||||
qreal a = angle_between(pi, pa, p);
|
||||
qreal b = angle_between(pi, pb, p);
|
||||
sum += (cotan(a) + cotan(b)) * (pi - p);
|
||||
sum += (cotan(a) + cotan(b)) * (p - pi);
|
||||
area += triangle_area(p, pi, pb);
|
||||
count++;
|
||||
}
|
||||
area /= 3.;
|
||||
return sum / (2.*area);
|
||||
return sum / (2.*count);
|
||||
}
|
||||
|
||||
|
||||
Eigen::SparseMatrix<qreal> laplacian_matrix(const MyMesh &mesh) {
|
||||
size_t n_verts = mesh.n_vertices();
|
||||
Eigen::SparseMatrix<qreal> mass(n_verts, n_verts);
|
||||
Eigen::SparseMatrix<qreal> cotangent(n_verts, n_verts);
|
||||
for (VertexHandle vert : mesh.vertices()) {
|
||||
if (mesh.is_boundary(vert)) continue;
|
||||
qreal sum = 0;
|
||||
qreal area = 0;
|
||||
qreal count = 0;
|
||||
Point p = mesh.point(vert);
|
||||
for (HalfedgeHandle v_vi : mesh.voh_range(vert)) {
|
||||
VertexHandle vi = mesh.to_vertex_handle(v_vi);
|
||||
Point pi = mesh.point(vi);
|
||||
HalfedgeHandle vi_v = mesh.opposite_halfedge_handle(v_vi);
|
||||
HalfedgeHandle v_va = mesh.next_halfedge_handle(vi_v);
|
||||
Point pa = mesh.point(mesh.to_vertex_handle(v_va));
|
||||
HalfedgeHandle vi_vb = mesh.next_halfedge_handle(v_vi);
|
||||
Point pb = mesh.point(mesh.to_vertex_handle(vi_vb));
|
||||
qreal a = angle_between(pi, pa, p);
|
||||
qreal b = angle_between(pi, pb, p);
|
||||
qreal w = -(cotan(a) + cotan(b)) / 2.;
|
||||
sum += w;
|
||||
cotangent.insert(vert.idx(), vi.idx()) = w;
|
||||
area += triangle_area(p, pi, pb);
|
||||
count++;
|
||||
}
|
||||
area /= 3.;
|
||||
cotangent.insert(vert.idx(), vert.idx()) = -sum;
|
||||
mass.insert(vert.idx(), vert.idx()) = 1. / area;
|
||||
}
|
||||
return mass * cotangent;
|
||||
}
|
||||
|
||||
|
||||
void smooth(MyMesh &mesh) {
|
||||
auto new_pos = OpenMesh::makeTemporaryProperty<VertexHandle, Point>(mesh);
|
||||
for (VertexHandle v : mesh.vertices()) {
|
||||
if (!mesh.is_boundary(v)) {
|
||||
new_pos[v] = mesh.point(v) + laplace_beltrami(mesh, v);
|
||||
}
|
||||
}
|
||||
for (VertexHandle v : mesh.vertices()) {
|
||||
if (!mesh.is_boundary(v)) {
|
||||
mesh.set_point(v, new_pos[v]);
|
||||
}
|
||||
}
|
||||
|
||||
// // Approche matricielle
|
||||
// size_t n_verts = mesh.n_vertices();
|
||||
// Eigen::SparseMatrix<qreal> mass(n_verts, n_verts);
|
||||
// Eigen::SparseMatrix<qreal> cotangent(n_verts, n_verts);
|
||||
// for (VertexHandle vert : mesh.vertices()) {
|
||||
// if (mesh.is_boundary(vert)) continue;
|
||||
// qreal sum = 0;
|
||||
// qreal area = 0;
|
||||
// qreal count = 0;
|
||||
// Point v = mesh.point(vert);
|
||||
// for (HalfedgeHandle v_vi : mesh.voh_range(vert)) {
|
||||
// VertexHandle vi_h = mesh.to_vertex_handle(v_vi);
|
||||
// Point vi = mesh.point(vi_h);
|
||||
|
||||
// HalfedgeHandle vi_v = mesh.opposite_halfedge_handle(v_vi);
|
||||
// HalfedgeHandle v_va = mesh.next_halfedge_handle(vi_v);
|
||||
// Point va = mesh.point(mesh.to_vertex_handle(v_va));
|
||||
|
||||
// HalfedgeHandle vi_vb = mesh.next_halfedge_handle(v_vi);
|
||||
// Point vb = mesh.point(mesh.to_vertex_handle(vi_vb));
|
||||
// qreal a = angle_between(vi, va, v);
|
||||
// qreal b = angle_between(v, vb, vi);
|
||||
// qreal w = -(1./qTan(a) + 1./qTan(b)) / 2.;
|
||||
// sum += w;
|
||||
// cotangent.insert(vert.idx(), vi_h.idx()) = w;
|
||||
// area += triangle_area(v, vi, vb);
|
||||
// count++;
|
||||
// auto new_pos = OpenMesh::makeTemporaryProperty<VertexHandle, Point>(mesh);
|
||||
// for (VertexHandle v : mesh.vertices()) {
|
||||
// if (!mesh.is_boundary(v)) {
|
||||
// new_pos[v] = mesh.point(v) - laplace_beltrami(mesh, v);
|
||||
// }
|
||||
// cotangent.insert(vert.idx(), vert.idx()) = -sum;
|
||||
// mass.insert(vert.idx(), vert.idx()) = 1. / area;
|
||||
// }
|
||||
// Eigen::SparseMatrix<qreal> laplacian = mass * cotangent;
|
||||
// Eigen::VectorX<qreal> X(n_verts), Y(n_verts), Z(n_verts);
|
||||
// for (VertexHandle vert : mesh.vertices()) {
|
||||
// if (mesh.is_boundary(vert)) continue;
|
||||
// size_t id = vert.idx();
|
||||
// Point p = mesh.point(vert);
|
||||
// X[id] = p[0];
|
||||
// Y[id] = p[1];
|
||||
// Z[id] = p[2];
|
||||
// }
|
||||
// X = laplacian * X;
|
||||
// Y = laplacian * Y;
|
||||
// Z = laplacian * Z;
|
||||
// for (VertexHandle vert : mesh.vertices()) {
|
||||
// if (mesh.is_boundary(vert)) continue;
|
||||
// size_t id = vert.idx();
|
||||
// mesh.set_point(vert, {X[id], Y[id], Z[id]});
|
||||
// for (VertexHandle v : mesh.vertices()) {
|
||||
// if (!mesh.is_boundary(v)) {
|
||||
// mesh.set_point(v, new_pos[v]);
|
||||
// }
|
||||
// }
|
||||
|
||||
// Approche matricielle
|
||||
Eigen::SparseMatrix<qreal> laplacian = laplacian_matrix(mesh);
|
||||
size_t n_verts = mesh.n_vertices();
|
||||
Eigen::VectorX<qreal> X(n_verts), Y(n_verts), Z(n_verts);
|
||||
for (VertexHandle vert : mesh.vertices()) {
|
||||
if (mesh.is_boundary(vert)) continue;
|
||||
size_t id = vert.idx();
|
||||
Point p = mesh.point(vert);
|
||||
X(id) = p[0];
|
||||
Y(id) = p[1];
|
||||
Z(id) = p[2];
|
||||
}
|
||||
X = laplacian * X;
|
||||
Y = laplacian * Y;
|
||||
Z = laplacian * Z;
|
||||
for (VertexHandle vert : mesh.vertices()) {
|
||||
if (mesh.is_boundary(vert)) continue;
|
||||
size_t id = vert.idx();
|
||||
Point offset {X(id), Y(id), Z(id)};
|
||||
mesh.set_point(vert, mesh.point(vert) - offset);
|
||||
}
|
||||
}
|
||||
|
|
10
src/util.h
10
src/util.h
|
@ -65,4 +65,12 @@ qreal triangle_area(Point p1, Point p2, Point p3) {
|
|||
qreal cotan(const qreal x);
|
||||
|
||||
|
||||
#endif
|
||||
template <typename T>
|
||||
constexpr const T& clamp(const T& v, const T& lo, const T& hi) {
|
||||
if (v < lo) return lo;
|
||||
if (v > hi) return hi;
|
||||
return v;
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
|
|
Loading…
Reference in New Issue