Commit 960f62f8 authored by Kenneth Vanhoey's avatar Kenneth Vanhoey
Browse files

qem for lightfield + cleanup of QEM file

parent 8e313b85
......@@ -22,6 +22,13 @@
* *
*******************************************************************************/
/**
* @file qem.h
* Header file for Quadric Error Metric classes.
*/
#ifndef __QEM__
#define __QEM__
......@@ -36,11 +43,11 @@
// Eigen includes
#include <Eigen/Dense>
#define CONST_VAL -5212368.54127
#define CONST_VAL -5212368.54127 // random value
namespace CGoGN {
//namespace Utils {
namespace Utils {
template <typename REAL>
class Quadric
......@@ -52,123 +59,42 @@ public:
typedef Geom::Vector<4,REAL> VEC4 ;
typedef Geom::Matrix<4,4,double> MATRIX44 ; // double is crucial here !
Quadric()
{
A.zero() ;
}
Quadric() ;
Quadric(int i) ;
Quadric(int i)
{
A.zero() ;
}
Quadric(VEC3& p1, VEC3& p2, VEC3& p3) ;
Quadric(VEC3& p1, VEC3& p2, VEC3& p3)
{
// compute the equation of the plane of the 3 points
Geom::Plane3D<REAL> plane(p1, p2, p3) ;
const VEC3& n = plane.normal() ;
void zero() ;
Geom::Vector<4,double> p = Geom::Vector<4,double>(n[0],n[1],n[2],plane.d()) ;
A = Geom::transposed_vectors_mult(p,p) ;
}
void operator= (const Quadric<REAL>& q) ;
Quadric& operator+= (const Quadric<REAL>& q) ;
Quadric& operator -= (const Quadric<REAL>& q) ;
Quadric& operator *= (REAL v) ;
Quadric& operator /= (REAL v) ;
void zero()
{
A.zero() ;
}
void operator= (const Quadric<REAL>& q)
{
A = q.A ;
}
Quadric& operator+= (const Quadric<REAL>& q)
{
A += q.A ;
return *this ;
}
Quadric& operator -= (const Quadric<REAL>& q)
{
A -= q.A ;
return *this ;
}
Quadric& operator *= (REAL v)
{
A *= v ;
return *this ;
}
Quadric& operator /= (REAL v)
{
A /= v ;
return *this ;
}
REAL operator() (const VEC4& v) const
{
return evaluate(v) ;
}
REAL operator() (const VEC3& v) const
{
VEC4 hv(v[0], v[1], v[2], 1.0f) ;
return evaluate(hv) ;
}
REAL operator() (const VEC4& v) const ;
REAL operator() (const VEC3& v) const ;
friend std::ostream& operator<<(std::ostream& out, const Quadric<REAL>& q)
{
out << q.A ;
return out ;
}
} ;
friend std::istream& operator>>(std::istream& in, Quadric<REAL>& q)
{
in >> q.A ;
return in ;
}
} ;
bool findOptimizedPos(VEC3& v)
{
VEC4 hv ;
bool b = optimize(hv) ;
if(b)
{
v[0] = hv[0] ;
v[1] = hv[1] ;
v[2] = hv[2] ;
}
return b ;
}
bool findOptimizedPos(VEC3& v) ;
private:
MATRIX44 A ;
REAL evaluate(const VEC4& v) const
{
// Double computation is crucial for stability
Geom::Vector<4, double> Av = A * v ;
return v * Av ;
}
bool optimize(VEC4& v) const
{
if (std::isnan(A(0,0)))
return false ;
MATRIX44 A2(A) ;
for(int i = 0; i < 3; ++i)
A2(3,i) = 0.0f ;
A2(3,3) = 1.0f ;
MATRIX44 Ainv ;
REAL det = A2.invert(Ainv) ;
if(det > -1e-6 && det < 1e-6)
return false ;
REAL evaluate(const VEC4& v) const ;
VEC4 right(0,0,0,1) ;
v = Ainv * right ;
return true;
}
bool optimize(VEC4& v) const ;
} ;
template <typename REAL, unsigned int N>
......@@ -180,100 +106,33 @@ public:
typedef Geom::Vector<N,REAL> VECN ;
typedef Geom::Vector<N+1,REAL> VECNp ;
QuadricNd()
{
A.zero() ;
b.zero() ;
c = 0 ;
}
QuadricNd() ;
QuadricNd(int i)
{
A.zero() ;
b.zero() ;
c = 0 ;
}
QuadricNd(int i) ;
QuadricNd(const VECN& p1_r, const VECN& p2_r, const VECN& p3_r)
{
const Geom::Vector<N,double>& p1 = p1_r ;
const Geom::Vector<N,double>& p2 = p2_r ;
const Geom::Vector<N,double>& p3 = p3_r ;
QuadricNd(const VECN& p1_r, const VECN& p2_r, const VECN& p3_r) ;
Geom::Vector<N,double> e1 = p2 - p1 ; e1.normalize() ;
Geom::Vector<N,double> e2 = (p3 - p1) - (e1*(p3-p1))*e1 ; e2.normalize() ;
void zero() ;
A.identity() ;
A -= Geom::transposed_vectors_mult(e1,e1) + Geom::transposed_vectors_mult(e2,e2) ;
void operator= (const QuadricNd<REAL,N>& q) ;
b = (p1*e1)*e1 + (p1*e2)*e2 - p1 ;
QuadricNd& operator+= (const QuadricNd<REAL,N>& q) ;
c = p1*p1 - pow((p1*e1),2) - pow((p1*e2),2) ;
}
QuadricNd& operator -= (const QuadricNd<REAL,N>& q) ;
void zero()
{
A.zero() ;
b.zero() ;
c = 0 ;
}
QuadricNd& operator *= (REAL v) ;
void operator= (const QuadricNd<REAL,N>& q)
{
A = q.A ;
b = q.b ;
c = q.c ;
}
QuadricNd& operator+= (const QuadricNd<REAL,N>& q)
{
A += q.A ;
b += q.b ;
c += q.c ;
return *this ;
}
QuadricNd& operator -= (const QuadricNd<REAL,N>& q)
{
A -= q.A ;
b -= q.b ;
c -= q.c ;
return *this ;
}
QuadricNd& operator *= (REAL v)
{
A *= v ;
b *= v ;
c *= v ;
return *this ;
}
QuadricNd& operator /= (REAL v)
{
A /= v ;
b /= v ;
c /= v ;
return *this ;
}
REAL operator() (const VECNp& v) const
{
VECN hv ;
for (unsigned int i = 0 ; i < N ; ++i)
hv[i] = v[i] ;
QuadricNd& operator /= (REAL v) ;
return evaluate(v) ;
}
REAL operator() (const VECNp& v) const ;
REAL operator() (const VECN& v) const
{
return evaluate(v) ;
}
REAL operator() (const VECN& v) const ;
friend std::ostream& operator<<(std::ostream& out, const QuadricNd<REAL,N>& q)
{
out << "(" << q.A << ", " << q.b << ", " << q.c << ")" ;
return out ;
}
} ;
friend std::istream& operator>>(std::istream& in, QuadricNd<REAL,N>& q)
{
......@@ -281,12 +140,9 @@ public:
in >> q.b ;
in >> q.c ;
return in ;
}
} ;
bool findOptimizedVec(VECN& v)
{
return optimize(v) ;
}
bool findOptimizedVec(VECN& v) ;
private:
// Double computation is crucial for stability
......@@ -294,28 +150,9 @@ private:
Geom::Vector<N,double> b ;
double c ;
REAL evaluate(const VECN& v) const
{
Geom::Vector<N, double> v_d = v ;
return v_d*A*v_d + 2.*(b*v_d) + c ;
}
bool optimize(VECN& v) const
{
if (std::isnan(A(0,0)))
return false ;
Geom::Matrix<N,N,double> Ainv ;
double det = A.invert(Ainv) ;
if(det > -1e-6 && det < 1e-6)
return false ;
v.zero() ;
v -= Ainv * b ;
REAL evaluate(const VECN& v) const ;
return true ;
}
bool optimize(VECN& v) const ;
} ;
template <typename REAL>
......@@ -326,168 +163,36 @@ public:
typedef Geom::Vector<3,REAL> VEC3 ;
QuadricHF()
{}
QuadricHF() ;
QuadricHF(int i)
{
m_A.resize(i,i) ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
m_b[c].resize(i) ;
m_c[c] = 0 ;
}
}
QuadricHF(const std::vector<VEC3>& v, const REAL& gamma, const REAL& alpha)
{
*this = QuadricHF(tensorsFromCoefs(v), gamma, alpha) ;
}
QuadricHF(int i) ;
QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REAL& alpha)
{
const unsigned int nbcoefs = ((T[0].order() + 1) * (T[0].order() + 2)) / 2. ;
//m_A.resize(nbcoefs, nbcoefs) ;
// 2D rotation
const Geom::Matrix33d R = rotateMatrix(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) ;
// parameterized integral on intersection
// build A, b and c
m_A = buildIntegralMatrix_A(alpha,nbcoefs) ; // Parameterized integral matrix A
Eigen::MatrixXd integ_b = buildIntegralMatrix_B(alpha,nbcoefs) ; // Parameterized integral matrix b
Eigen::MatrixXd integ_c = buildIntegralMatrix_C(alpha,nbcoefs) ; // Parameterized integral matrix c
// for (unsigned int i = 0 ; i < nbcoefs ; ++i)
// for (unsigned int j = 0 ; j < nbcoefs ; ++j)
// std::cout << "(" << i << "," << j << ")=" << m_A(i,j) << std::endl ;
//
// for (unsigned int i = 0 ; i < nbcoefs ; ++i)
// for (unsigned int j = 0 ; j < nbcoefs ; ++j)
// std::cout << "(" << i << "," << j << ")=" << integ_b(i,j) << std::endl ;
//
// for (unsigned int i = 0 ; i < nbcoefs ; ++i)
// for (unsigned int j = 0 ; j < nbcoefs ; ++j)
// std::cout << "(" << i << "," << j << ")=" << integ_c(i,j) << std::endl ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
Eigen::VectorXd v ; v.resize(nbcoefs) ;
for (unsigned int i = 0 ; i < nbcoefs ; ++i) // copy into vector
v[i] = coefsR[i][c] ;
m_b[c] = integ_b * v ; // Vector b
m_c[c] = v.transpose() * (integ_c * v) ; // Constant c
}
}
~QuadricHF()
{
//delete m_A ;
//for (unsigned int c = 0 ; c < 3 ; ++c)
// delete m_b[c] ;
}
QuadricHF(const std::vector<VEC3>& v, const REAL& gamma, const REAL& alpha) ;
void zero()
{
m_A.setZero() ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
m_b[c].setZero() ;
m_c[c] = 0 ;
}
}
QuadricHF& operator= (const QuadricHF<REAL>& q)
{
m_A = q.m_A ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
m_b[c] = q.m_b[c] ;
m_c[c] = q.m_c[c] ;
}
QuadricHF(const Geom::Tensor3d* T, const REAL& gamma, const REAL& alpha) ;
return *this ;
}
~QuadricHF() ;
QuadricHF& operator+= (const QuadricHF<REAL>& q)
{
assert(((m_A.cols() == q.m_A.cols()) && (m_A.rows() == q.m_A.rows()) && m_b[0].size() == q.m_b[0].size()) || !"Incompatible add of matrices") ;
if (!(m_A.cols() == q.m_A.cols()) && (m_A.rows() == q.m_A.rows()) && (m_b[0].size() == q.m_b[0].size()))
return *this ;
m_A += q.m_A ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
m_b[c] += q.m_b[c] ;
m_c[c] += q.m_c[c] ;
}
return *this ;
}
QuadricHF& operator -= (const QuadricHF<REAL>& q)
{
assert(((m_A.cols() == q.m_A.cols()) && (m_A.rows() == q.m_A.rows()) && m_b[0].size() == q.m_b[0].size()) || !"Incompatible substraction of matrices") ;
if ((m_A.cols() == q.m_A.cols()) && (m_A.rows() == q.m_A.rows()) && (m_b[0].size() == q.m_b[0].size()))
return *this ;
m_A -= q.m_A ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
m_b[c] -= q.m_b[c] ;
m_c[c] -= q.m_c[c] ;
}
return *this ;
}
QuadricHF& operator *= (const REAL& v)
{
m_A *= v ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
m_b[c] *= v ;
m_c[c] *= v ;
}
return *this ;
}
QuadricHF& operator /= (const REAL& v)
{
const REAL& inv = 1. / v ;
(*this) *= inv ;
return *this ;
}
// REAL operator() (const VECNp& v) const
// {
// VECN hv ;
// for (unsigned int i = 0 ; i < 3 ; ++i)
// hv[i] = v[i] ;
//
// return evaluate(v) ;
// }
REAL operator() (const std::vector<VEC3>& coefs) const
{
return evaluate(coefs) ;
}
void zero() ;
QuadricHF& operator= (const QuadricHF<REAL>& q) ;
QuadricHF& operator+= (const QuadricHF<REAL>& q) ;
QuadricHF& operator -= (const QuadricHF<REAL>& q) ;
QuadricHF& operator *= (const REAL& v) ;
QuadricHF& operator /= (const REAL& v) ;
// REAL operator() (const VECNp& v) const ;
REAL operator() (const std::vector<VEC3>& coefs) const ;
friend std::ostream& operator<<(std::ostream& out, const QuadricHF<REAL>& q)
{
// TODO out << "(" << q.m_A << ", " << q.m_b << ", " << q.m_c << ")" ;
return out ;
}
} ;
friend std::istream& operator>>(std::istream& in, QuadricHF<REAL>& q)
{
......@@ -496,22 +201,9 @@ public:
// in >> q.b ;
// in >> q.c ;
return in ;
}
} ;
bool findOptimizedCoefs(std::vector<VEC3>& coefs)
{
// unsigned int nbcoefs = 6 ;
// for (unsigned int i = 0 ; i < nbcoefs ; ++i)
// for (unsigned int j = 0 ; j < nbcoefs ; ++j)
// std::cout << "A(" << i << "," << j << ")=" << m_A(i,j) << std::endl ;
//
// for (unsigned int i = 0 ; i < nbcoefs ; ++i)
// std::cout << "b(" << i << ")=" << m_b[0][i] << std::endl ;
//
// std::cout << "c=" << m_c[0] << std::endl ;
return optimize(coefs) ;
}
bool findOptimizedCoefs(std::vector<VEC3>& coefs) ;
private:
// Double computation is crucial for stability
......@@ -519,767 +211,30 @@ private:
Eigen::VectorXd m_b[3] ;
double m_c[3] ;
REAL evaluate(const std::vector<VEC3>& coefs) const
{
VEC3 res ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
Eigen::VectorXd tmp(coefs.size()) ;
for (unsigned int i = 0 ; i < coefs.size() ; ++i)
tmp[i] = coefs[i][c] ;
res[c] = tmp.transpose() * m_A * tmp ; // A
res[c] -= 2. * (m_b[c]).transpose() * tmp ; // - 2b
res[c] += m_c[c] ; // + c
}
res /= 2*M_PI ; // max integral value over hemisphere
return (res[0] + res[1] + res[2]) / 3. ;
}
bool optimize(std::vector<VEC3>& coefs) const
{
coefs.resize(m_b[0].size()) ;
REAL evaluate(const std::vector<VEC3>& coefs) const ;
if (fabs(m_A.determinant()) < 1e-6 )
return false ;
bool optimize(std::vector<VEC3>& coefs) const ;
Eigen::MatrixXd Ainv = m_A.inverse() ;
Geom::Matrix33d rotateMatrix(const REAL& gamma) ;
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] ;
}
Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R) ;
return true ;
}
Eigen::MatrixXd buildIntegralMatrix_A(const REAL& alpha, unsigned int size) ;
Eigen::MatrixXd buildIntegralMatrix_B(const REAL& alpha, unsigned int size) ;
Geom::Matrix33d rotateMatrix(const REAL& gamma)
{
Geom::Matrix33d R ;
R.identity() ;
R(0,0) = cos(gamma) ;
R(0,1) = sin(gamma) ;
R(1,0) = -sin(gamma) ;
R(1,1) = cos(gamma) ;
Eigen::MatrixXd buildIntegralMatrix_C(const REAL& alpha, unsigned int size) ;
return R ;
}
void fillTensor(Geom::Tensor3d& T) ;
Geom::Tensor3d rotate(const Geom::Tensor3d& T, const Geom::Matrix33d& R)
{
Geom::Tensor3d Tp(T.order()) ;
std::vector<unsigned int> p ; p.resize(T.order(), 0) ;
for (unsigned int i = 0 ; i < T.nbElem() ; ++i)
{
REAL S = 0 ;
std::vector<unsigned int> q ; q.resize(T.order(), 0) ;
for (unsigned int j = 0 ; j < T.nbElem() ; ++j)
{
REAL P = T[j] ;
for (unsigned int k = 0 ; k < T.order() ; ++k)
P *= R(q[k],p[k]) ;