Commit f61a4b22 authored by Kenneth Vanhoey's avatar Kenneth Vanhoey

lightfield decimation correction

parent 87b0e875
...@@ -69,7 +69,7 @@ void Approximator_Frame<PFP>::approximate(Dart d) ...@@ -69,7 +69,7 @@ void Approximator_Frame<PFP>::approximate(Dart d)
segmentNew -= this->m_position[d] ; segmentNew -= this->m_position[d] ;
// Orthogonal projection of v0-v onto v0-v1 : get coefficient t // Orthogonal projection of v0-v onto v0-v1 : get coefficient t
REAL t = (segment * segmentNew) / (segment.norm()*segment.norm()) ; REAL t = (segment * segmentNew) / segment.norm() ;
t = std::max (std::min (t , REAL(1)) , REAL(0) ) ; // clamp it to [0,1] t = std::max (std::min (t , REAL(1)) , REAL(0) ) ; // clamp it to [0,1]
VEC3 n1, n2 ; VEC3 n1, n2 ;
......
...@@ -320,13 +320,13 @@ template <unsigned int DIM, typename T> ...@@ -320,13 +320,13 @@ template <unsigned int DIM, typename T>
inline Vector<DIM,T> slerp(const Vector<DIM,T> &v1, const Vector<DIM,T> &v2, const T &t) { inline Vector<DIM,T> slerp(const Vector<DIM,T> &v1, const Vector<DIM,T> &v2, const T &t) {
Vector<DIM,T> res ; Vector<DIM,T> res ;
T omega = acos(v1 * v2) ; T omega = acos(v1 * v2) ; // assume v1,v2 normalized
T den = sin(omega) ; T den = sin(omega) ;
if (-1e-8 < den && den < 1e-8) if (-1e-8 < den && den < 1e-8)
return t < 0.5 ? v1 : v2 ; return t < 0.5 ? v1 : v2 ;
T f1 = sin(( T(1)-t ) * omega ) / den ; T f1 = sin(( T(1)-t ) * omega ) / den ; // assume 0 <= t <= 1
T f2 = sin( t * omega ) / den ; T f2 = sin( t * omega ) / den ;
res += f1 * v1 ; res += f1 * v1 ;
......
...@@ -41,9 +41,9 @@ public: ...@@ -41,9 +41,9 @@ public:
QuadricRGBfunctions(); QuadricRGBfunctions();
QuadricRGBfunctions(int i); QuadricRGBfunctions(int i);
QuadricRGBfunctions(const QuadricRGBfunctions&); QuadricRGBfunctions(const QuadricRGBfunctions&);
QuadricRGBfunctions(const RGBFUNCTIONS&, REAL gamma = REAL(0), REAL alpha = REAL(0)); QuadricRGBfunctions(const RGBFUNCTIONS&, const REAL gamma = REAL(0), const REAL alpha = REAL(0)) ;
virtual ~QuadricRGBfunctions(); virtual ~QuadricRGBfunctions() {} ;
REAL operator() (const RGBFUNCTIONS&) const; REAL operator() (const RGBFUNCTIONS&) const;
...@@ -51,8 +51,8 @@ public: ...@@ -51,8 +51,8 @@ public:
void operator += (const QuadricRGBfunctions&) ; void operator += (const QuadricRGBfunctions&) ;
void operator -= (const QuadricRGBfunctions&) ; void operator -= (const QuadricRGBfunctions&) ;
void operator *= (REAL v) ; void operator *= (const REAL v) ;
void operator /= (REAL v) ; void operator /= (const REAL v) ;
void zero () ; void zero () ;
...@@ -60,10 +60,10 @@ public: ...@@ -60,10 +60,10 @@ public:
friend std::istream& operator>> (std::istream &in, const QuadricRGBfunctions&) {return in;}; friend std::istream& operator>> (std::istream &in, const QuadricRGBfunctions&) {return in;};
private : private :
void buildIntegralMatrix_A(MATRIX66 &, REAL alpha) const; void buildIntegralMatrix_A(MATRIX66 &, const REAL alpha) const;
void buildIntegralMatrix_b(MATRIX66 &, REAL alpha) const; void buildIntegralMatrix_b(MATRIX66 &, const REAL alpha) const;
void buildIntegralMatrix_c(MATRIX66 &, REAL alpha) const; void buildIntegralMatrix_c(MATRIX66 &, const REAL alpha) const;
void buildRotateMatrix(MATRIX66 &N, REAL gamma) const; void buildRotateMatrix(MATRIX66 &N, const REAL gamma) const;
}; };
......
...@@ -16,10 +16,12 @@ std::string QuadricRGBfunctions<REAL>::CGoGNnameOfType() { ...@@ -16,10 +16,12 @@ std::string QuadricRGBfunctions<REAL>::CGoGNnameOfType() {
template <typename REAL> template <typename REAL>
QuadricRGBfunctions<REAL>::QuadricRGBfunctions() { QuadricRGBfunctions<REAL>::QuadricRGBfunctions() {
for (unsigned col = RED; col < BLUE+1 ; ++col) { for (unsigned col = RED; col < BLUE+1 ; ++col) {
for (unsigned i = 0; i < 6; ++i) { for (unsigned i = 0; i < 6; ++i) {
b[col][i] = (REAL)0;
for (unsigned j = 0; j < 6; ++j) for (unsigned j = 0; j < 6; ++j)
A[col](i,j) = REAL(0); A[col](i,j) = REAL(0) ;
b[col][i] = REAL(0) ;
} }
c[col] = REAL(0); c[col] = REAL(0);
...@@ -34,76 +36,55 @@ QuadricRGBfunctions<REAL>::QuadricRGBfunctions(int i) { ...@@ -34,76 +36,55 @@ QuadricRGBfunctions<REAL>::QuadricRGBfunctions(int i) {
template <typename REAL> template <typename REAL>
QuadricRGBfunctions<REAL>::QuadricRGBfunctions(const QuadricRGBfunctions& q) { QuadricRGBfunctions<REAL>::QuadricRGBfunctions(const QuadricRGBfunctions& q) {
for (unsigned col = RED; col < BLUE+1 ; ++col) { for (unsigned col = RED; col < BLUE+1 ; ++col) {
for (unsigned i = 0; i < 6; ++i) { for (unsigned i = 0; i < 6; ++i) {
b[col][i] = q.b[col][i];
for (unsigned j = 0; j < 6; ++j) for (unsigned j = 0; j < 6; ++j)
A[col](i,j) = q.A[col](i,j); A[col](i,j) = q.A[col](i,j) ;
b[col][i] = q.b[col][i] ;
} }
c[col] = q.c[col]; c[col] = q.c[col] ;
} }
} }
template <typename REAL> template <typename REAL>
QuadricRGBfunctions<REAL>::QuadricRGBfunctions(const RGBFUNCTIONS& cf, REAL gamma, REAL alpha) { QuadricRGBfunctions<REAL>::QuadricRGBfunctions(const RGBFUNCTIONS& cf, const REAL gamma, const REAL alpha) {
MATRIX66 R1,R2_A,R2_b,R2_c; MATRIX66 R1,R2_A,R2_b,R2_c;
// Matrice de rotation 1
buildRotateMatrix(R1,gamma);
// std::cout <<"R1 : " << R1 << std::endl ;
// std::cout <<"alpha : " << alpha << std::endl ;
// Matrice de rotation 2 + intégrale
buildIntegralMatrix_A(R2_A,alpha);
buildIntegralMatrix_b(R2_b,alpha);
buildIntegralMatrix_b(R2_c,alpha);
// std::cout << "A : " << R2_A << std::endl ; buildRotateMatrix(R1,gamma); // Rotation 1
// std::cout << "b : " << R2_b << std::endl ;
// std::cout << "c : " << R2_c << std::endl ;
buildIntegralMatrix_A(R2_A,alpha); // Parameterized integral matrix A
buildIntegralMatrix_b(R2_b,alpha); // Parameterized integral matrix b
buildIntegralMatrix_c(R2_c,alpha); // Parameterized integral matrix c
// Quadrique (A,b,c) tel que L*A*Lt - 2*b*Lt + c = ERROR // Quadric (A,b,c) => L*A*Lt - 2*b*Lt + c = ERROR
for (unsigned col = RED; col < BLUE+1; ++col) { for (unsigned col = RED; col < BLUE+1; ++col) {
// Rotation 1 Geom::Vector<6,REAL> function; // get function coefficients
Geom::Vector<6,REAL> function;
if (!cf.getSubVectorH(col,0,function)) if (!cf.getSubVectorH(col,0,function))
assert(!"QuadricRGBfunctions::constructor") ; assert(!"QuadricRGBfunctions::constructor") ;
VEC6 coefs = R1 * function ; // Multiply coefs VEC6 coefs = R1 * function ; // Rotation 1
// Calcul de A (intégrale)
A[col] = R2_A ;
// Calcul de b (intégrale+rotation sur 1 vec)
b[col] = R2_b * function ;
// Calcul de c (intégrale+rotation sur 2 vec) A[col] = R2_A ; // Matrix A : integral
c[col] = function * (R2_c * function) ; b[col] = R2_b * function ; // Vector b : integral + rotation on 1 vector
c[col] = function * (R2_c * function) ; // Scalar c : integral + rotation on 2 vectors
} }
} }
template <typename REAL>
QuadricRGBfunctions<REAL>::~QuadricRGBfunctions() {
}
template <typename REAL> template <typename REAL>
REAL QuadricRGBfunctions<REAL>::operator() (const RGBFUNCTIONS& cf) const { REAL QuadricRGBfunctions<REAL>::operator() (const RGBFUNCTIONS& cf) const {
REAL res = REAL(0); REAL res = REAL(0);
for (unsigned col = RED; col < BLUE+1; ++col) { for (unsigned col = RED; col < BLUE+1; ++col) {
Geom::Vector<6,REAL> function ; Geom::Vector<6,REAL> function ; // Get function coefficients
if (!cf.getSubVectorH(col,0,function)) if (!cf.getSubVectorH(col,0,function))
assert (!"QuadricRGBfunctions::getSubVectorH") ; assert (!"QuadricRGBfunctions::getSubVectorH") ;
VEC6 Al = A[col] * function;
REAL res_local = REAL(0) ; REAL res_local = REAL(0) ;
res_local += (function * Al); // l*A*lt res_local += function * (A[col] * function) ; // l*A*lt
res_local -= 2 * (function * b[col]); res_local -= 2 * (function * b[col]) ; // -2*l*b
res_local += c[col]; res_local += c[col] ; // c
// res = l*A*lT - 2*l*b + c
res += res_local; res += res_local;
} }
...@@ -114,19 +95,15 @@ template <typename REAL> ...@@ -114,19 +95,15 @@ template <typename REAL>
bool QuadricRGBfunctions<REAL>::findOptimizedRGBfunctions(RGBFUNCTIONS& cf) const { bool QuadricRGBfunctions<REAL>::findOptimizedRGBfunctions(RGBFUNCTIONS& cf) const {
for (unsigned col = RED; col < BLUE+1 ; ++col) { for (unsigned col = RED; col < BLUE+1 ; ++col) {
MATRIX66 Ainv,M; MATRIX66 Ainv ;
M = A[col];
REAL det = M.invert(Ainv);
REAL det = A[col].invert(Ainv) ; // Invert matrix
if(det > -1e-8 && det < 1e-8) if(det > -1e-8 && det < 1e-8)
return false ; return false ; // invert failed
VEC6 coefs;
coefs = Ainv * b[col]; VEC6 coefs = Ainv * b[col]; // function = A^(-1) * b
if (!cf.setSubVectorH(col,0,coefs)) if (!cf.setSubVectorH(col,0,coefs)) // save in argument cf
assert (!"QuadricRGBfunctions::findOptimizedRGBfunctions(cf) setSubVector failed") ; assert (!"QuadricRGBfunctions::findOptimizedRGBfunctions(cf) setSubVector failed") ;
} }
...@@ -134,18 +111,99 @@ bool QuadricRGBfunctions<REAL>::findOptimizedRGBfunctions(RGBFUNCTIONS& cf) cons ...@@ -134,18 +111,99 @@ bool QuadricRGBfunctions<REAL>::findOptimizedRGBfunctions(RGBFUNCTIONS& cf) cons
} }
template <typename REAL> template <typename REAL>
void QuadricRGBfunctions<REAL>::buildIntegralMatrix_A(MATRIX66 &M, REAL alpha) const { void QuadricRGBfunctions<REAL>::buildIntegralMatrix_A(MATRIX66 &M, const REAL alpha) const {
// Int(phi=0..pi)(theta=0..pi-alpha) variables^2 dTheta dPhi if alpha > 0
// Int(phi=0..pi)(theta= -alpha..pi) variables^2 dTheta dPhi if alpha < 0
REAL alphaAbs = alpha > 0 ? alpha : -alpha; REAL alphaAbs = alpha > 0 ? alpha : -alpha;
const REAL pi = 3.141592 ; const REAL pi = 3.141592 ;
const REAL cosAbs = cos(alphaAbs); // = - cos(alpha) si alpha < 0
const REAL cosinus = cos(alpha); const REAL cosinus = cos(alpha);
const REAL sinus = sin(alpha); const REAL sinus = sin(alpha);
const REAL sinAbs = sin(alphaAbs); // = - sin(alpha) si alpha < 0
// Line 1 // Line 1
M(0,0) = 2.0*(pi-alphaAbs)/5.0f; M(0,0) = 2.0*(pi-alphaAbs)/5.0f;
M(0,1) = 2.0 * (pi - alphaAbs - cosinus*sinAbs) / 15.0 ;
M(0,2) = 0;
M(0,3) = 0;
M(0,4) = sinAbs*pi / 8.0;
M(0,5) = 2.0*(pi-alphaAbs)/3.0;
// Line 2
M(1,0) = M(0,1);
M(1,1) = (-4.0 * sinAbs*cosinus*cosinus*cosinus + 6.0 * (pi - cosinus*sinAbs - alphaAbs)) / 15.0 ;
M(1,2) = REAL(0);
M(1,3) = REAL(0);
M(1,4) = (sinus*cosinus*cosinus*pi + 2*sinus*pi) / 8.0;
M(1,5) = 2.0*(pi-cosinus*sinAbs-alphaAbs)/3.0;
// LINE 3
M(2,0) = REAL(0);
M(2,1) = REAL(0);
M(2,2) = 2.0*(pi - alphaAbs - cosinus*sinAbs)/15.0;
M(2,3) = pi*sinus/8.0f;
M(2,4) = REAL(0) ;
M(2,5) = REAL(0) ;
// Line 4
M(3,0) = REAL(0) ;
M(3,1) = REAL(0) ;
M(3,2) = M(2,3);
M(3,3) = 2.0*(pi-alphaAbs)/3.0 ;
M(3,4) = REAL(0) ;
M(3,5) = REAL(0) ;
// Line 5
M(4,0) = M(0,4) ;
M(4,1) = M(1,4) ;
M(4,2) = REAL(0) ;
M(4,3) = REAL(0) ;
M(4,4) = 2.0 * (pi - cosinus*sinAbs - alphaAbs ) / 3.0 ;
M(4,5) = pi*sinus / 2.0 ;
// Line 6
M(5,0) = M(0,5) ;
M(5,1) = M(1,5) ;
M(5,2) = REAL(0) ;
M(5,3) = REAL(0) ;
M(5,4) = M(4,5) ;
M(5,5) = 2.0*(pi-alphaAbs);
/* // Line 1
M(0,0) = 2.0*(pi-alphaAbs)/5.0f;
M(0,1) = 2.0 * (pi - alphaAbs - cosAbs*sinus) / 15.0 ; M(0,1) = 2.0 * (pi - alphaAbs - cosAbs*sinus) / 15.0 ;
M(0,2) = 0; M(0,2) = 0;
...@@ -220,53 +278,60 @@ void QuadricRGBfunctions<REAL>::buildIntegralMatrix_A(MATRIX66 &M, REAL alpha) c ...@@ -220,53 +278,60 @@ void QuadricRGBfunctions<REAL>::buildIntegralMatrix_A(MATRIX66 &M, REAL alpha) c
M(5,4) = M(4,5) ; M(5,4) = M(4,5) ;
M(5,5) = 2.0*(pi-alphaAbs); M(5,5) = 2.0*(pi-alphaAbs);
*/
} }
template <typename REAL> template <typename REAL>
void QuadricRGBfunctions<REAL>::buildIntegralMatrix_b(MATRIX66 &M, REAL alpha) const { void QuadricRGBfunctions<REAL>::buildIntegralMatrix_b(MATRIX66 &M, const REAL alpha) const {
// Int(phi=0..pi)(theta=0..pi-alpha) variables*variablesRotated dTheta dPhi * coefs if alpha > 0
// Int(phi=0..pi)(theta= -alpha..pi) variables*variablesRotated dTheta dPhi * coefs if alpha < 0
REAL alphaAbs = alpha > 0 ? alpha : -alpha; REAL alphaAbs = alpha > 0 ? alpha : -alpha;
const REAL pi = 3.141592 ; const REAL pi = 3.141592 ;
const REAL cosAbs = cos(alphaAbs); // = - cos(alpha) si alpha < 0 const REAL cosinus = cos(alpha) ;
const REAL cosinus = cos(alpha); const REAL cos2 = cosinus*cosinus ;
const REAL sinus = sin(alpha); const REAL cos3 = cos2*cosinus ;
const REAL cos4 = cos3*cosinus ;
const REAL cos5 = cos4*cosinus ;
const REAL sinus = sin(alpha) ;
const REAL sinAbs = sin(alphaAbs) ; // = - sin(alpha) si alpha < 0
// Line 1 // Line 1
M(0,0) = 2.0*(pi-alphaAbs)/5.0f; M(0,0) = 2.0*(pi-alphaAbs)/5.0f;
M(0,1) = 2.0 * (pi - alphaAbs - cosAbs*sinus) / 15.0 ; M(0,1) = ( 6*cosinus*sinAbs - 8*sinAbs*cos3 - 2*alphaAbs + 2*pi ) / 15.0 ;
M(0,2) = 0; M(0,2) = 0;
M(0,3) = 0; M(0,3) = 0;
M(0,4) = sinus*pi / 8.0; M(0,4) = (sinus*pi + 2*pi*cosinus*sinus) / 8.0 ;
M(0,5) = 2.0*(pi-alphaAbs)/3.0; M(0,5) = 2.0*(pi-alphaAbs)/3.0;
// Line 2 // Line 2
M(1,0) = M(0,1); M(1,0) = 2 * (pi - cosinus*sinAbs - alphaAbs ) / 15.0 ;
M(1,1) = 2 * (pi - alphaAbs + 2*(cosAbs*cosAbs*pi-alphaAbs*cosAbs*cosAbs))/15.0 - 2.0*cosAbs*sinus/3.0 ; M(1,1) = ( 6*cosinus*sinAbs - 2*alphaAbs + 2*pi - 16*sinAbs*cos5+4*cos2*pi - 4*alphaAbs * cos2 ) / 15.0 ;
M(1,2) = 0; M(1,2) = 0;
M(1,3) = 0; M(1,3) = 0;
M(1,4) = (sinus*pi - 2*sinus*cosinus*pi) / 8.0; M(1,4) = (sinus*pi + 2* (pi*sinus*cos3+pi*cosinus*sinus)) / 8.0 ;
M(1,5) = 2.0*(pi-cosAbs*sinus-alphaAbs)/3.0; M(1,5) = 2.0*(pi-cosinus*sinAbs-alphaAbs) / 3.0;
// LINE 3 // LINE 3
M(2,0) = REAL(0); M(2,0) = REAL(0);
M(2,1) = REAL(0); M(2,1) = REAL(0);
M(2,2) = 2.0*(cosinus*pi - cosinus*alphaAbs + (alpha < 0 ? sinus : -sinus))/15.0; M(2,2) = 2 * (sinAbs - cosinus*alphaAbs + cosinus*pi - 2*cos2*sinAbs) / 15.0 ;
M(2,3) = -pi*sinus/8.0f; M(2,3) = pi*sinus / 8.0f;
M(2,4) = REAL(0) ; M(2,4) = REAL(0) ;
...@@ -277,7 +342,7 @@ void QuadricRGBfunctions<REAL>::buildIntegralMatrix_b(MATRIX66 &M, REAL alpha) c ...@@ -277,7 +342,7 @@ void QuadricRGBfunctions<REAL>::buildIntegralMatrix_b(MATRIX66 &M, REAL alpha) c
M(3,1) = REAL(0) ; M(3,1) = REAL(0) ;
M(3,2) = - M(2,3); M(3,2) = ( sinus*pi + 2*pi*cosinus*sinus ) / 8.0 ;
M(3,3) = 2.0*(pi-alphaAbs)/3.0 ; M(3,3) = 2.0*(pi-alphaAbs)/3.0 ;
...@@ -286,76 +351,83 @@ void QuadricRGBfunctions<REAL>::buildIntegralMatrix_b(MATRIX66 &M, REAL alpha) c ...@@ -286,76 +351,83 @@ void QuadricRGBfunctions<REAL>::buildIntegralMatrix_b(MATRIX66 &M, REAL alpha) c
M(3,5) = REAL(0) ; M(3,5) = REAL(0) ;
// Line 5 // Line 5
M(4,0) = - M(0,4) ; M(4,0) = pi*sinus / 8.0 ;
M(4,1) = (2*pi*cosinus*sinus - sinus*pi ) / 8.0 ; M(4,1) = (sinus*pi + 4 * sinus * cos4 + 2*pi*cosinus*sinus) / 8.0 ;
M(4,2) = REAL(0) ; M(4,2) = REAL(0) ;
M(4,3) = REAL(0) ; M(4,3) = REAL(0) ;
M(4,4) = 2.0 * (cosinus*pi - cosinus*alphaAbs + (alpha < 0 ? sinus : -sinus) ) / 3.0 ; M(4,4) = 2*(sinAbs - cosinus*alphaAbs + cosinus*pi - 2*cos2*sinAbs) / 3.0 ;
M(4,5) = pi*sinus / 2.0 ; M(4,5) = pi*sinus / 2.0 ;
// Line 6 // Line 6
M(5,0) = M(0,5) ; M(5,0) = M(0,5) ;
M(5,1) = M(1,5) ; M(5,1) = 2*cosinus*sinAbs + 2*(pi-4*sinAbs*cos3-alphaAbs) / 3.0 ;
M(5,2) = REAL(0) ; M(5,2) = REAL(0) ;
M(5,3) = REAL(0) ; M(5,3) = REAL(0) ;
M(5,4) = M(4,5) ; M(5,4) = sinus*pi / 2.0 + pi*cosinus*sinus ;
M(5,5) = 2.0*(pi-alphaAbs); M(5,5) = 2.0*(pi-alphaAbs);
} }
template <typename REAL> template <typename REAL>
void QuadricRGBfunctions<REAL>::buildIntegralMatrix_c(MATRIX66 &M, REAL alpha) const { void QuadricRGBfunctions<REAL>::buildIntegralMatrix_c(MATRIX66 &M, const REAL alpha) const {
// coefs * Int(phi=0..pi)(theta=0..pi-alpha) variablesRotated^2 dTheta dPhi * coefs if alpha > 0
// coefs * Int(phi=0..pi)(theta= -alpha..pi) variablesRotated^2 dTheta dPhi * coefs if alpha < 0
REAL alphaAbs = alpha > 0 ? alpha : -alpha; REAL alphaAbs = alpha > 0 ? alpha : -alpha;
const REAL pi = 3.141592 ; const REAL pi = 3.141592 ;
const REAL cosAbs = cos(alphaAbs); // = - cos(alpha) si alpha < 0 const REAL cosinus = cos(alpha);
const REAL cos = cos(alpha); const REAL cos2 = cosinus*cosinus ;
const REAL sin = sin(alpha); const REAL cos3 = cos2*cosinus ;
const REAL cos5 = cos2*cos3 ;
const REAL cos7 = cos2*cos5 ;
const REAL sinus = sin(alpha);
const REAL sinAbs = sin(alphaAbs); // = - sin(alpha) si alpha < 0
// Line 1 // Line 1
M(0,0) = 2.0*(pi-alphaAbs)/5.0f; M(0,0) = 2.0*(pi-alphaAbs)/5.0f;
M(0,1) = 2.0 * (pi - alphaAbs - cosAbs*sin) / 15.0 ;