ccolin/mod_geo-tp
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

merge curvature stuff into this

Browse Source
This commit is contained in:
papush!
2021-11-12 17:08:19 +01:00
parent d51945ecaf
commit d01ba47b65
19 changed files with 672 additions and 129 deletions
Download Patch File Download Diff File Expand all files Collapse all files

265
src/curvature.cpp Normal file
View File

@ -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));
}
}

163
src/curvature.h Normal file
View File

@ -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

7
src/main.cpp
View File

@ -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();

15
src/mesh_processor.cpp
View File

@ -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);
}
}

2
src/mesh_processor.h
View File

@ -12,7 +12,7 @@ class MeshProcessor : public QObject {
public:
MyMesh mesh;
std::vector<std::vector<HalfedgeHandle>> holes;
MyMesh patch;
MeshProcessor(const QString &path);

42
src/mesh_view.cpp
View File

@ -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 */

6
src/mesh_view.h
View File

@ -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);
};

8
src/mesh_viewer.cpp
View File

@ -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();

6
src/mesh_viewer.h
View File

@ -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);

18
src/my_mesh.h
View File

@ -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;

46
src/quad_patch.cpp Normal file
View File

@ -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;
}

11
src/quad_patch.h Normal file
View File

@ -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

4
src/shader.frag
View File

@ -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);
}

4
src/shader.vert
View File

@ -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);
}

128
src/smoothing.cpp
View File

@ -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
View File

@ -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