Commit d68c5b28 authored by Sylvain Thery's avatar Sylvain Thery

Merge branch 'master' of cgogn:~vanhoey/CGoGN

parents aae2b76b 3b5c2a93
......@@ -62,7 +62,7 @@ void decimate(
std::vector<VertexAttribute<typename PFP::VEC3> *>& position,
unsigned int nbWantedVertices,
const FunctorSelect& selected = allDarts,
VertexAttribute<typename PFP::VEC3> *edgeErrors = NULL,
EdgeAttribute<typename PFP::REAL> *edgeErrors = NULL,
void (*callback_wrapper)(void*, const void*) = NULL, void *callback_object = NULL
) ;
......
......@@ -2010,9 +2010,6 @@ void EdgeSelector_Lightfield<PFP>::computeEdgeInfo(Dart d, EdgeInfo& einfo)
double alpha = alpha1 + alpha2 ;
if (isnan(alpha))
std::cerr << "Nan: " << m_frameN[d] << " ; " << m_frameN[dd] << " ; " << newFN << std::endl ;
assert(m_quadricHF.isValid() | !"EdgeSelector_Lightfield<PFP>::computeEdgeInfo: quadricHF is not valid") ;
Utils::QuadricHF<REAL> quadHF = m_quadricHF[d] ;
......@@ -2023,11 +2020,11 @@ void EdgeSelector_Lightfield<PFP>::computeEdgeInfo(Dart d, EdgeInfo& einfo)
const REAL& errLF = quadHF(newHF) ; // function coefficients
// Check if errated values appear
if (errG < -1e-10 || errAngle < -1e-10 || errLF < -1e-10)
if (errG < -1e-6 || errAngle < -1e-6 || errLF < -1e-6)
einfo.valid = false ;
else
{
einfo.it = edges.insert(std::make_pair(std::max(errG + errAngle + errLF, REAL(0)), d)) ;
einfo.it = edges.insert(std::make_pair(std::max(/*errG +*/ errAngle + errLF, REAL(0)), d)) ;
einfo.valid = true ;
}
}
......
......@@ -939,9 +939,6 @@ void HalfEdgeSelector_Lightfield<PFP>::computeHalfEdgeInfo(Dart d, HalfEdgeInfo&
double alpha = alpha1 + alpha2 ;
if (isnan(alpha))
std::cerr << "Nan: " << m_frameN[d] << " ; " << m_frameN[dd] << " ; " << newFN << std::endl ;
assert(m_quadricHF.isValid() | !"EdgeSelector_Lightfield<PFP>::computeEdgeInfo: quadricHF is not valid") ;
Utils::QuadricHF<REAL> quadHF = m_quadricHF[d] ;
......
......@@ -204,10 +204,10 @@ void Approximator_HemiFuncCoefs<PFP>::approximate(Dart d)
bool opt = m_quadricHF[d].findOptimizedCoefs(coefs) ;
// copy result
for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
this->m_approx[i][d] = opt ? coefs[i] : (m_coefs[i]->operator[](d) + m_coefs[i]->operator[](dd))/2 ; // if fail average coefs (TODO better)
this->m_approx[i][d] = coefs[i] ;
if (!opt)
std::cerr << "LightfieldApproximator::Inversion failed: not treated correctly" << std::endl ;
std::cerr << "LightfieldApproximator::Inversion failed: coefs are averaged" << std::endl ;
}
}
......@@ -337,8 +337,8 @@ void Approximator_HemiFuncCoefsHalfEdge<PFP>::approximate(Dart d)
const REAL gamma1 = ((B1 * T) > 0 ? 1 : -1) * acos( std::max(std::min(1.0f, T1 * T ), -1.0f)) ; // angle positif ssi
const REAL gamma2 = ((B2 * T) > 0 ? 1 : -1) * acos( std::max(std::min(1.0f, T2 * T ), -1.0f)) ; // -PI/2 < angle(i,j1) < PI/2 ssi i*j1 > 0
// Rotation dans le sens trigo autour de l'axe T (N1 --> N)
const REAL alpha1 = ((N * B1prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N1), -1.0f) ) ; // angle positif ssi
const REAL alpha2 = ((N * B2prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N2), -1.0f) ) ; // PI/2 < angle(j1',n) < -PI/2 ssi j1'*n < 0
const REAL alpha1 = (N == N1) ? 0 : ((N * B1prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N1), -1.0f) ) ; // angle positif ssi
const REAL alpha2 = (N == N2) ? 0 : ((N * B2prime) > 0 ? -1 : 1) * acos( std::max(std::min(1.0f, N * N2), -1.0f) ) ; // PI/2 < angle(j1',n) < -PI/2 ssi j1'*n < 0
assert (m_quadricHF.isValid() || !"LightfieldApproximator:approximate quadricHF is not valid") ;
......@@ -383,10 +383,13 @@ void Approximator_HemiFuncCoefsHalfEdge<PFP>::approximate(Dart d)
bool opt = m_quadricHF[d].findOptimizedCoefs(coefs) ;
// copy result
for (unsigned int i = 0 ; i < m_nbCoefs ; ++i)
this->m_approx[i][d] = opt ? coefs[i] : (m_coefs[i]->operator[](d) + m_coefs[i]->operator[](dd))/2 ; // if fail average coefs (TODO better)
this->m_approx[i][d] = coefs[i] ;
if (!opt)
std::cerr << "LightfieldApproximator::Inversion failed: not treated correctly" << std::endl ;
{
std::cerr << "LightfieldApproximatorHalfCollapse::Optimization failed (should never happen since no optim is done)" << std::endl ;
std::cout << alpha1 << std::endl ;
}
// Add second one for error evaluation
m_quadricHF[d] += Utils::QuadricHF<REAL>(coefs2, gamma2, alpha2) ;
......
......@@ -31,6 +31,9 @@
#include "OpenNL/sparse_matrix.h"
#include "OpenNL/full_vector.h"
#include <Eigen/Core>
#include <Eigen/Eigenvalues>
namespace CGoGN
{
......
......@@ -313,6 +313,7 @@ void computeCurvatureVertex_NormalCycles(
typedef typename PFP::VEC3 VEC3 ;
typedef typename PFP::REAL REAL ;
// collect the normal cycle tensor
Algo::Selection::Collector_WithinSphere<PFP> neigh(map, position, radius, thread) ;
neigh.collectAll(dart) ;
neigh.computeArea() ;
......@@ -321,14 +322,14 @@ void computeCurvatureVertex_NormalCycles(
typename PFP::MATRIX33 tensor(0) ;
// inside
// collect edges inside the neighborhood
const std::vector<Dart>& vd1 = neigh.getInsideEdges() ;
for (std::vector<Dart>::const_iterator it = vd1.begin(); it != vd1.end(); ++it)
{
const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(map, *it, position) ;
tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) ;
}
// border
// collect edges crossing the neighborhood's border
const std::vector<Dart>& vd2 = neigh.getBorder() ;
for (std::vector<Dart>::const_iterator it = vd2.begin(); it != vd2.end(); ++it)
{
......@@ -340,43 +341,40 @@ void computeCurvatureVertex_NormalCycles(
tensor /= neigh.getArea() ;
long int n = 3, lda = 3, info, lwork = 9 ;
char jobz='V', uplo = 'U' ;
float work[9] ;
float w[3] ;
float a[3*3] = {
tensor(0,0), 0.0f, 0.0f,
tensor(1,0), tensor(1,1), 0.0f,
tensor(2,0), tensor(2,1), tensor(2,2)
} ;
// Solve eigenproblem
ssyev_(&jobz, &uplo, (integer*)&n, a, (integer*)&lda, w, work, (integer*)&lwork, (integer*)&info) ;
// Check for convergence
if(info > 0)
std::cerr << "clapack ssyev_ failed to compute eigenvalues : exit with code " << info << std::endl ;
// sort eigen components : w[s[0]] has minimal absolute value ; kmin = w[s[1]] <= w[s[2]] = kmax
// solve eigen problem
Eigen::Matrix3f m3 ;
m3 << tensor(0,0) , tensor(0,1) , tensor(0,2) , tensor(1,0) , tensor(1,1) , tensor(1,2) , tensor(2,0) , tensor(2,1) , tensor(2,2) ;
Eigen::SelfAdjointEigenSolver<Eigen::Matrix3f> solver (m3);
Eigen::Vector3f ev = solver.eigenvalues();
Eigen::Matrix3f evec = solver.eigenvectors();
// sort eigen components : ev[s[0]] has minimal absolute value ; kmin = ev[s[1]] <= ev[s[2]] = kmax
int s[3] = {0, 1, 2} ;
int tmp ;
if (abs(w[s[2]]) < abs(w[s[1]])) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }
if (abs(w[s[1]]) < abs(w[s[0]])) { tmp = s[0] ; s[0] = s[1] ; s[1] = tmp ; }
if (w[s[2]] < w[s[1]]) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }
if (abs(ev[s[2]]) < abs(ev[s[1]])) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }
if (abs(ev[s[1]]) < abs(ev[s[0]])) { tmp = s[0] ; s[0] = s[1] ; s[1] = tmp ; }
if (ev[s[2]] < ev[s[1]]) { tmp = s[1] ; s[1] = s[2] ; s[2] = tmp ; }
kmin[dart] = w[s[1]] ;
kmax[dart] = w[s[2]] ;
// set curvatures from sorted eigen components
// warning : Kmin and Kmax are switched w.r.t. kmin and kmax
// normal direction : minimal absolute eigen value
VEC3& dirNormal = Knormal[dart] ;
dirNormal[0] = evec(0,s[0]);
dirNormal[1] = evec(1,s[0]);
dirNormal[2] = evec(2,s[0]);
if (dirNormal * normal[dart] < 0) dirNormal *= -1; // change orientation
// min curvature
kmin[dart] = ev[s[1]] ;
VEC3& dirMin = Kmin[dart] ;
dirMin[0] = a[3*s[2]];
dirMin[1] = a[3*s[2]+1];
dirMin[2] = a[3*s[2]+2]; // warning : Kmin and Kmax are switched
dirMin[0] = evec(0,s[2]);
dirMin[1] = evec(1,s[2]);
dirMin[2] = evec(2,s[2]);
// max curvature
kmax[dart] = ev[s[2]] ;
VEC3& dirMax = Kmax[dart] ;
dirMax[0] = a[3*s[1]];
dirMax[1] = a[3*s[1]+1];
dirMax[2] = a[3*s[1]+2]; // warning : Kmin and Kmax are switched
VEC3& dirNormal = Knormal[dart] ;
dirNormal[0] = a[3*s[0]];
dirNormal[1] = a[3*s[0]+1];
dirNormal[2] = a[3*s[0]+2];
if (dirNormal * normal[dart] < 0)
dirNormal *= -1; // change orientation
dirMax[0] = evec(0,s[1]);
dirMax[1] = evec(1,s[1]);
dirMax[2] = evec(2,s[1]);
}
......
......@@ -1010,12 +1010,21 @@ bool MeshTablesSurface<PFP>::importPlySLFgenericBin(const std::string& filename,
attrNames.push_back(positions.name()) ;
VertexAttribute<typename PFP::VEC3> *frame = new VertexAttribute<typename PFP::VEC3>[3] ;
if (tangent)
{
frame[0] = m_map.template addAttribute<typename PFP::VEC3, VERTEX>("frameT") ; // Tangent
attrNames.push_back(frame[0].name()) ;
}
if (binormal)
{
frame[1] = m_map.template addAttribute<typename PFP::VEC3, VERTEX>("frameB") ; // Bitangent
attrNames.push_back(frame[0].name()) ;
}
if (normal)
{
frame[2] = m_map.template addAttribute<typename PFP::VEC3, VERTEX>("frameN") ; // Normal
attrNames.push_back(frame[0].name()) ;
attrNames.push_back(frame[1].name()) ;
attrNames.push_back(frame[2].name()) ;
}
VertexAttribute<typename PFP::VEC3> *PBcoefs = NULL, *SHcoefs = NULL ;
if (PTM)
......@@ -1064,14 +1073,15 @@ bool MeshTablesSurface<PFP>::importPlySLFgenericBin(const std::string& filename,
fp.read((char*)properties,nbProps * propSize) ;
// positions
if (nbProps > 2)
positions[id] = VEC3(properties[0],properties[1],properties[2]) ; // position
if (tangent && binormal && normal) // == if (nbprops > 11)
for (unsigned int k = 0 ; k < 3 ; ++k) // frame
for (unsigned int l = 0 ; l < 3 ; ++l)
frame[k][id][l] = (typename PFP::REAL)(properties[3+(3*k+l)]) ;
/*
for (unsigned int k = 0 ; k < 3 ; ++k) // coefficients
for (unsigned int l = 0 ; l < nbCoefs ; ++l)
*/
for (unsigned int l = 0 ; l < nbCoefs ; ++l) // coefficients
for (unsigned int k = 0 ; k < 3 ; ++k)
{
......@@ -1080,6 +1090,7 @@ bool MeshTablesSurface<PFP>::importPlySLFgenericBin(const std::string& filename,
else /* if SH */
SHcoefs[l][id][k] = (typename PFP::REAL)(properties[12+(3*l+k)]) ;
}
unsigned int cur = 12+3*nbCoefs ;
for (unsigned int k = 0 ; k < nbRemainders ; ++k) // remaining data
remainders[k][id] = (typename PFP::REAL)(properties[cur + k]) ;
......
......@@ -65,7 +65,7 @@ private:
std::vector<VSplit<PFP>*> m_splits ;
unsigned int m_cur ;
Algo::Decimation::HalfEdgeApproximator<PFP, VEC3, EDGE>* m_positionApproximator ;
Algo::Decimation::Approximator<PFP, VEC3, EDGE>* m_positionApproximator ;
bool m_initOk ;
......
......@@ -108,7 +108,7 @@ ProgressiveMesh<PFP>::ProgressiveMesh(
if(! (*it)->init())
m_initOk = false ;
if((*it)->getApproximatedAttributeName() == "position")
m_positionApproximator = reinterpret_cast<Algo::Decimation::HalfEdgeApproximator<PFP, VEC3, EDGE>*>(*it) ;
m_positionApproximator = reinterpret_cast<Algo::Decimation::Approximator<PFP, VEC3, EDGE>*>(*it) ;
}
CGoGNout << "..done" << CGoGNendl ;
......
......@@ -583,6 +583,8 @@ private:
Eigen::MatrixXd m_A ; /*!< The first QuadricHF member matrix A */
Eigen::VectorXd m_b[3] ; /*!< The second QuadricHF member vector b */
double m_c[3] ; /*!< The third QuadricHF member scalar c */
std::vector<VEC3> m_coefs ; /*!< The coefficients in cas optim fails */
bool m_noAlphaRot ; /*!< If alpha = 0 then optim will fail */
/*!
* \brief method to evaluate the error for a given lightfield function
......@@ -593,18 +595,6 @@ private:
*/
REAL evaluate(const std::vector<VEC3>& coefs) const ;
/*!
* \brief method to deduce an optimal coefficients in space
* w.r.t. the current QuadricHF.
*
* \param coefs will contain the optimal result (if it can be computed)
*
* \return true if an optimal result was correctly computed
*/
bool optimize(std::vector<VEC3>& coefs) const ;
// Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) ;
/*!
* \brief method to build a rotate matrix (rotation in tangent plane)
* given angle gamma
......
......@@ -300,11 +300,13 @@ QuadricNd<REAL,N>::optimize(VECN& v) const
}
template <typename REAL>
QuadricHF<REAL>::QuadricHF()
QuadricHF<REAL>::QuadricHF():
m_noAlphaRot(false)
{}
template <typename REAL>
QuadricHF<REAL>::QuadricHF(int i)
QuadricHF<REAL>::QuadricHF(int i):
m_noAlphaRot(false)
{
m_A.resize(i,i) ;
for (unsigned int c = 0 ; c < 3 ; ++c)
......@@ -323,18 +325,17 @@ QuadricHF<REAL>::QuadricHF(const std::vector<VEC3>& v, const REAL& gamma, const
}
template <typename REAL>
QuadricHF<REAL>::QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REAL& alpha)
QuadricHF<REAL>::QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REAL& alpha):
m_noAlphaRot(fabs(alpha) < 1e-13)
{
const unsigned int nbcoefs = ((T[0].order() + 1) * (T[0].order() + 2)) / 2. ;
//m_A.resize(nbcoefs, nbcoefs) ;
// 2D rotation
const Geom::Matrix33d R = buildRotateMatrix(gamma) ;
Geom::Tensor3d* Trot = new Geom::Tensor3d[3] ;
for (unsigned int c = 0 ; c < 3 ; ++c)
Trot[c] = rotate(T[c],R) ;
std::vector<VEC3> coefsR = coefsFromTensors(Trot) ;
m_coefs = coefsFromTensors(Trot) ;
delete[] Trot ;
// parameterized integral on intersection
......@@ -348,7 +349,7 @@ QuadricHF<REAL>::QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REA
{
Eigen::VectorXd v ; v.resize(nbcoefs) ;
for (unsigned int i = 0 ; i < nbcoefs ; ++i) // copy into vector
v[i] = coefsR[i][c] ;
v[i] = m_coefs[i][c] ;
m_b[c] = integ_b * v ; // Vector b
m_c[c] = v.transpose() * (integ_c * v) ; // Constant c
......@@ -369,6 +370,8 @@ QuadricHF<REAL>::zero()
m_b[c].setZero() ;
m_c[c] = 0 ;
}
m_coefs.clear() ;
m_noAlphaRot = false ;
}
template <typename REAL>
......@@ -381,6 +384,8 @@ QuadricHF<REAL>::operator= (const QuadricHF<REAL>& q)
m_b[c] = q.m_b[c] ;
m_c[c] = q.m_c[c] ;
}
m_coefs = q.m_coefs ;
m_noAlphaRot = q.m_noAlphaRot ;
return *this ;
}
......@@ -399,6 +404,11 @@ QuadricHF<REAL>::operator+= (const QuadricHF<REAL>& q)
m_b[c] += q.m_b[c] ;
m_c[c] += q.m_c[c] ;
}
m_coefs.resize(m_coefs.size()) ;
for (unsigned int i = 0 ; i < m_coefs.size() ; ++i)
m_coefs[i] += q.m_coefs[i] ;
m_noAlphaRot &= q.m_noAlphaRot ;
return *this ;
}
......@@ -418,6 +428,12 @@ QuadricHF<REAL>::operator -= (const QuadricHF<REAL>& q)
m_c[c] -= q.m_c[c] ;
}
m_coefs.resize(m_coefs.size()) ;
for (unsigned int i = 0 ; i < m_coefs.size() ; ++i)
m_coefs[i] -= q.m_coefs[i] ;
m_noAlphaRot &= q.m_noAlphaRot ;
return *this ;
}
......@@ -425,6 +441,7 @@ template <typename REAL>
QuadricHF<REAL>&
QuadricHF<REAL>::operator *= (const REAL& v)
{
std::cout << "Warning: QuadricHF<REAL>::operator *= should not be used !" << std::endl ;
m_A *= v ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
......@@ -439,6 +456,7 @@ template <typename REAL>
QuadricHF<REAL>&
QuadricHF<REAL>::operator /= (const REAL& v)
{
std::cout << "Warning: QuadricHF<REAL>::operator /= should not be used !" << std::endl ;
const REAL& inv = 1. / v ;
(*this) *= inv ;
......@@ -457,7 +475,24 @@ template <typename REAL>
bool
QuadricHF<REAL>::findOptimizedCoefs(std::vector<VEC3>& coefs)
{
return optimize(coefs) ;
coefs.resize(m_b[0].size()) ;
if (fabs(m_A.determinant()) < 1e-10 )
{
coefs = m_coefs ;
return m_noAlphaRot ; // if true inversion failed (normal!) and m_coefs forms a valid solution
}
Eigen::MatrixXd Ainv = m_A.inverse() ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
Eigen::VectorXd tmp(m_b[0].size()) ;
tmp = Ainv * m_b[c] ;
for (unsigned int i = 0 ; i < m_b[c].size() ; ++i)
coefs[i][c] = tmp[i] ;
}
return true ;
}
template <typename REAL>
......@@ -480,28 +515,6 @@ QuadricHF<REAL>::evaluate(const std::vector<VEC3>& coefs) const
return (res[0] + res[1] + res[2]) / 3. ;
}
template <typename REAL>
bool
QuadricHF<REAL>::optimize(std::vector<VEC3>& coefs) const
{
coefs.resize(m_b[0].size()) ;
if (fabs(m_A.determinant()) < 1e-10 )
return false ;
Eigen::MatrixXd Ainv = m_A.inverse() ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
Eigen::VectorXd tmp(m_b[0].size()) ;
tmp = Ainv * m_b[c] ;
for (unsigned int i = 0 ; i < m_b[c].size() ; ++i)
coefs[i][c] = tmp[i] ;
}
return true ;
}
template <typename REAL>
Geom::Matrix33d
QuadricHF<REAL>::buildRotateMatrix(const REAL& gamma)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment