Commit 0dd92877 authored by Pierre Kraemer's avatar Pierre Kraemer

First version of half-edge radiance error metric driven decimation

parent 3eab3531
......@@ -41,13 +41,11 @@ private:
DartAttribute<HalfEdgeInfo, PFP2::MAP> halfEdgeInfo;
VertexAttribute<Utils::Quadric<PFP2::REAL>, PFP2::MAP> m_quadric;
VertexAttribute<PFP2::VEC3, PFP2::MAP> m_avgColor;
// VertexAttribute<PFP2::VEC3, PFP2::MAP> m_avgColor;
unsigned int m_nb_coefs;
SphericalFunctionIntegratorCartesian m_integrator;
double* m_n;
double* m_workspace;
std::multimap<float, Dart> halfEdges;
typename std::multimap<float, Dart>::iterator cur;
......@@ -56,31 +54,22 @@ private:
void updateHalfEdgeInfo(Dart d);
void computeHalfEdgeInfo(Dart d, HalfEdgeInfo& einfo);
void recomputeQuadric(const Dart d);
/*
PFP2::VEC3 computeGradientLFerror(const Dart& v0, const Dart& v1) const;
PFP2::VEC3 computeDirDerivativeLFerror(const Dart& v0, const Dart& v1);
PFP2::VEC3 computeSquaredLightfieldDifferenceNumerical(const Dart& d1, const Dart& d2, const PFP2::VEC3& N) const;
PFP2::VEC3 computeSquaredLightfieldDifferenceAnalytical(const Dart& d1, const Dart& d2) const;
PFP2::VEC3 computeGradient(const PFP2::VEC3& P0, const PFP2::VEC3& Pi, const PFP2::VEC3& Pj, const Dart& v0, const Dart& v1, const Dart& vi, const Dart& vj, unsigned int channel) const;
PFP2::VEC3 integrateDscalGrad(const PFP2::VEC3& d, const unsigned int& K, const PFP2::VEC3& N, const PFP2::VEC3& ei, const PFP2::VEC3& ej,
const PFP2::VEC3* coefs1, const PFP2::VEC3& T1, const PFP2::VEC3& B1, const PFP2::VEC3& N1, const PFP2::VEC3& avg1,
const PFP2::VEC3* coefsi, const PFP2::VEC3& Ti, const PFP2::VEC3& Bi, const PFP2::VEC3& Ni, const PFP2::VEC3& avgi,
const PFP2::VEC3* coefsj, const PFP2::VEC3& Tj, const PFP2::VEC3& Bj, const PFP2::VEC3& Nj, const PFP2::VEC3& avgj) const;
PFP2::VEC3 integrateDlf(const PFP2::VEC3& d, const unsigned int& K, const PFP2::VEC3& N, const PFP2::VEC3& ei, const PFP2::VEC3& ej,
const std::vector<PFP2::VEC3*> coefs0, const PFP2::VEC3& T0, const PFP2::VEC3& B0, const PFP2::VEC3& N0, const PFP2::VEC3& avg0,
const std::vector<PFP2::VEC3*> coefs1, const PFP2::VEC3& T1, const PFP2::VEC3& B1, const PFP2::VEC3& N1, const PFP2::VEC3& avg1,
const std::vector<PFP2::VEC3*> coefsi, const PFP2::VEC3& Ti, const PFP2::VEC3& Bi, const PFP2::VEC3& Ni, const PFP2::VEC3& avgi,
const std::vector<PFP2::VEC3*> coefsj, const PFP2::VEC3& Tj, const PFP2::VEC3& Bj, const PFP2::VEC3& Nj, const PFP2::VEC3& avgj) const;
PFP2::REAL computeIntegral(const double *avgi, const PFP2::VEC3& ti, const PFP2::VEC3& bi, const PFP2::VEC3& ni, unsigned int nbCoefs, const std::vector<double>& coefs) const;
static double dispScalGrad (double x, double y, double z, void* data);
static double dlf (double x, double y, double z, void* data);
static double evalF(double* N, double* avg, unsigned int nb, bool isSH, double* T, double* B, double* coefs, double& x, double& y, double& z);
static void cart2spherical(double u, double v, double& theta, double& phi);
static double SquaredDifferenceOfCartesianFunctions (double x, double y, double z, void* data);
static bool isInDomain(double x, double y, double z, void *data);
static double CartesianFunction (double x, double y, double z, void* data);
*/
typename PFP::REAL computeRadianceError(Dart d);
static bool isInHemisphere(double x, double y, double z, void* u)
{ // true iff [x,y,z] and u have the same direction
REAL* n = (REAL*)(u);
return x*n[0] + y*n[1] + z*n[2] >= 0.0;
}
static double SHEvalCartesian_Error(double x, double y, double z, void* u)
{
SH& e = *(SH*)(u);
VEC3 c = e.evaluate_at(x, y, z);
return c.norm2();
}
public:
HalfEdgeSelector_Radiance(
MAP& m,
......@@ -98,24 +87,19 @@ public:
m_positionApproximator(posApprox),
m_normalApproximator(normApprox),
m_radianceApproximator(radApprox),
m_nb_coefs(0),
m_n(NULL),
m_workspace(NULL)
m_nb_coefs(0)
{
halfEdgeInfo = m.template checkAttribute<HalfEdgeInfo, DART, PFP2::MAP>("halfEdgeInfo");
m_quadric = m.template checkAttribute<Utils::Quadric<PFP2::REAL>, VERTEX, PFP2::MAP>("QEMquadric");
m_avgColor = m.template checkAttribute<PFP2::VEC3, VERTEX, PFP2::MAP>("avgColor");
m_n = new double[3];
// m_avgColor = m.template checkAttribute<PFP2::VEC3, VERTEX, PFP2::MAP>("avgColor");
}
~HalfEdgeSelector_Radiance()
{
this->m_map.removeAttribute(halfEdgeInfo);
this->m_map.removeAttribute(m_quadric);
this->m_map.removeAttribute(m_avgColor);
// this->m_map.removeAttribute(m_avgColor);
m_integrator.Release();
delete[] m_n;
delete[] m_workspace;
}
Algo::Surface::Decimation::SelectorType getType() { return Algo::Surface::Decimation::S_OTHER; }
......@@ -149,6 +133,65 @@ public:
(*errors)[d] = -1;
}
}
class ErrorEvaluator
{
private:
double wA;
double wB;
double wC;
SH* sA;
SH* sB;
SH* sC;
SH* sRef;
VEC3 nA;
VEC3 nB;
VEC3 nC;
VEC3 nRef;
VEC3 cA;
VEC3 cB;
VEC3 cC;
VEC3 cRef;
public:
ErrorEvaluator(
double p_wA,
double p_wB,
double p_wC,
SH* p_sA,
SH* p_sB,
SH* p_sC,
SH* p_sRef,
VEC3 p_nA,
VEC3 p_nB,
VEC3 p_nC,
VEC3 p_nRef,
VEC3 p_cA,
VEC3 p_cB,
VEC3 p_cC,
VEC3 p_cRef
)
{
wA = p_wA; wB = p_wB; wC = p_wC;
sA = p_sA; sB = p_sB; sC = p_sC; sRef = p_sRef;
nA = p_nA; nB = p_nB; nC = p_nC; nRef = p_nRef;
cA = p_cA; cB = p_cB; cC = p_cC; cRef = p_cRef;
}
void set_eval_direction (double x, double y, double z) // fix the direction in which the error has to be evaluated
{
SH::set_eval_direction (x,y,z);
}
double evaluate () // evaluates at a fixed direction
{
VEC3 color_difference ( sRef->evaluate() - wA*sA->evaluate()- wB*sB->evaluate() - wC*sC->evaluate() );
return color_difference.squaredNorm();
}
};
};
} // namespace SCHNApps
......
......@@ -36,11 +36,7 @@ bool HalfEdgeSelector_Radiance<PFP>::init()
m_nb_coefs = SH::get_nb_coefs();
unsigned int ruleId = m_nb_coefs <= 6 ? 17 : (m_nb_coefs <= 10) ? 23 : 29 ;
m_integrator.Init(ruleId) ;
// allocate workspaces
m_workspace = new double[60+12*m_nb_coefs] ;
m_integrator.Init(29) ;
// init QEM quadrics
for (Vertex v : allVerticesOf(m))
......@@ -198,15 +194,13 @@ void HalfEdgeSelector_Radiance<PFP>::computeHalfEdgeInfo(Dart d, HalfEdgeInfo& h
// Compute errors
Utils::Quadric<REAL> quadGeom ;
quadGeom += m_quadric[d] ; // compute the sum of the
quadGeom += m_quadric[dd] ; // two vertices quadrics
quadGeom += m_quadric[v0] ; // compute the sum of the
quadGeom += m_quadric[v1] ; // two vertices quadrics
const REAL geom = quadGeom(newPos) ;
const REAL rad = 0;
// const REAL rad = computeDirDerivativeLFerror(v0,v1).norm()/sqrt(3) ;
// const REAL rad = computeGradientLFerror(v0,v1).norm()/sqrt(3) ;
const REAL rad = computeRadianceError(d);
const REAL t = 1.0 ;
const REAL t = 0.01 ;
const REAL err = t*geom + (1-t)*rad ;
// Check if errated values appear
......@@ -218,793 +212,57 @@ void HalfEdgeSelector_Radiance<PFP>::computeHalfEdgeInfo(Dart d, HalfEdgeInfo& h
heinfo.valid = true ;
}
}
/*
template <typename PFP>
VEC3 HalfEdgeSelector_Radiance<PFP>::computeGradientLFerror(const Dart& v0, const Dart& v1) const
{
MAP& m = this->m_map ;
const VEC3& P0 = this->m_approx[m_approxindex_pos]->getAttr(m_attrindex_pos)[v0] ;
const VEC3& P1 = this->m_approx[m_approxindex_pos]->getAttr(m_attrindex_pos)[v1] ;
const VEC3& N = this->m_approx[m_approxindex_FN]->getApprox(v1,m_attrindex_FN) ;
const VEC3& T1 = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[v1] ;
const VEC3& B1 = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[v1] ;
const VEC3& N1 = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[v1] ;
VEC3 coefs1[m_K] ;
for (unsigned int k = 0 ; k < m_K ; ++k)
{
coefs1[k] = this->m_approx[m_approxindex_HF[k]]->getAttr(m_attrindex_HF[k])[v1] ;
}
Traversor2VF<MAP> tf(m,v0) ; // all faces around vertex v0
VEC3 count ;
REAL areaSum = 0 ;
const VEC3 d = P1 - P0 ; // displacement vector
for (Dart fi = tf.begin() ; fi != tf.end() ; fi = tf.next()) // foreach "blue" face
{
// get the data
const Dart& vi = m.phi1(fi) ;
const Dart& vj = m.phi_1(fi) ;
const VEC3& Pi = this->m_position[vi] ;
const VEC3& Pj = this->m_position[vj] ;
const VEC3& Ti = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[vi] ;
const VEC3& Bi = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[vi] ;
const VEC3& Ni = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[vi] ;
const VEC3& Tj = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[vj] ;
const VEC3& Bj = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[vj] ;
const VEC3& Nj = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[vj] ;
VEC3 coefsi[m_K], coefsj[m_K] ;
for (unsigned int k = 0 ; k < m_K ; ++k)
{
coefsi[k] = this->m_approx[m_approxindex_HF[k]]->getAttr(m_attrindex_HF[k])[vi] ;
coefsj[k] = this->m_approx[m_approxindex_HF[k]]->getAttr(m_attrindex_HF[k])[vj] ;
}
// utils
const VEC3 ei = P0 - Pj ;
const VEC3 ej = Pi - P0 ;
//const VEC3 e0 = Pj - Pi ;
const REAL areaIJ0sq = (ei ^ ej).norm2() ;
const REAL areaIJ0 = sqrt(areaIJ0sq)/2. ;
areaSum += areaIJ0 ;
const VEC3 displacementE = integrateDscalGrad(d, m_K, N, ei, ej,
coefs1, T1, B1, N1, m_avgColor[v1],
coefsi, Ti, Bi, Ni, m_avgColor[vi],
coefsj, Tj, Bj, Nj, m_avgColor[vj]) ;
count += displacementE * areaIJ0 ;
}
const VEC3 lfDiff = computeSquaredLightfieldDifferenceNumerical(v0,v1,N) ;
return lfDiff*areaSum + count ;
}
template <typename PFP>
VEC3 HalfEdgeSelector_Radiance<PFP>::computeDirDerivativeLFerror(const Dart& v0, const Dart& v1)
typename PFP::REAL HalfEdgeSelector_Radiance<PFP>::computeRadianceError(Dart d)
{
MAP& m = this->m_map ;
const VEC3& P0 = this->m_approx[m_approxindex_pos]->getAttr(m_attrindex_pos)[v0] ;
const VEC3& P1 = this->m_approx[m_approxindex_pos]->getAttr(m_attrindex_pos)[v1] ;
const VEC3& N = this->m_approx[m_approxindex_FN]->getApprox(v1,m_attrindex_FN) ;
const VEC3& T1 = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[v1] ;
const VEC3& B1 = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[v1] ;
const VEC3& N1 = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[v1] ;
const VEC3& T0 = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[v0] ;
const VEC3& B0 = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[v0] ;
const VEC3& N0 = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[v0] ;
std::vector<VEC3*> coefs0, coefs1 ;
coefs0.resize(m_K) ;
coefs1.resize(m_K) ;
for (unsigned int k = 0 ; k < m_K ; ++k)
{
coefs0[k] = &(this->m_approx[m_approxindex_HF[k]]->getAttr(m_attrindex_HF[k])[v0]) ;
coefs1[k] = &(this->m_approx[m_approxindex_HF[k]]->getAttr(m_attrindex_HF[k])[v1]) ;
}
Traversor2VF<MAP> tf(m,v0) ; // all faces around vertex v0
VEC3 error ;
const VEC3 d = P1 - P0 ; // displacement vector
std::vector<VEC3*> coefsi, coefsj ;
coefsi.resize(m_K) ;
coefsj.resize(m_K) ;
for (Dart fi = tf.begin() ; fi != tf.end() ; fi = tf.next()) // foreach "blue" face
{
// get the data
const Dart& vi = m.phi1(fi) ;
const Dart& vj = m.phi_1(fi) ;
const VEC3& Pi = this->m_position[vi] ;
const VEC3& Pj = this->m_position[vj] ;
const VEC3& Ti = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[vi] ;
const VEC3& Bi = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[vi] ;
const VEC3& Ni = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[vi] ;
const VEC3& Tj = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[vj] ;
const VEC3& Bj = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[vj] ;
const VEC3& Nj = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[vj] ;
//VEC3 coefsi[m_K], coefsj[m_K] ;
for (unsigned int k = 0 ; k < m_K ; ++k)
{
coefsi[k] = &(this->m_approx[m_approxindex_HF[k]]->getAttr(m_attrindex_HF[k])[vi]) ;
coefsj[k] = &(this->m_approx[m_approxindex_HF[k]]->getAttr(m_attrindex_HF[k])[vj]) ;
}
// utils
const VEC3 ei = P1 - Pj ;
const VEC3 ej = Pi - P1 ;
// area x integration
error += (ei ^ ej).norm()/2 * integrateDlf(d, m_K, N, ei, ej,
coefs0, T0, B0, N0, m_avgColor[v0],
coefs1, T1, B1, N1, m_avgColor[v1],
coefsi, Ti, Bi, Ni, m_avgColor[vi],
coefsj, Tj, Bj, Nj, m_avgColor[vj]) ;
}
return error ;
}
template <typename PFP>
VEC3 HalfEdgeSelector_Radiance<PFP>::computeSquaredLightfieldDifferenceNumerical(const Dart& v0, const Dart& v1, const VEC3& N) const
{
// get two frames
const VEC3& T1 = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[v0] ;
const VEC3& T2 = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[v1] ;
const VEC3& B1 = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[v0] ;
const VEC3& B2 = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[v1] ;
const VEC3& N1 = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[v0] ;
const VEC3& N2 = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[v1] ;
// get coefs of v1 and v2
std::vector<VEC3> coefs1, coefs2 ;
coefs1.resize(m_K) ; coefs2.resize(m_K) ;
for (unsigned int i = 0 ; i < m_K ; ++i)
{
coefs1[i] = this->m_approx[m_approxindex_HF[i]]->getAttr(m_attrindex_HF[i])[v0] ;
coefs2[i] = this->m_approx[m_approxindex_HF[i]]->getAttr(m_attrindex_HF[i])[v1] ;
}
// fill workspace
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[i] = double(N1[i]) ; // n
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[4+i] = double(T1[i]) ;
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[7+i] = double(B1[i]) ;
m_workspace[10] = double(m_K) ;
// fill workspace 2
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[11+m_K+i] = double(N2[i]) ; // n
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[15+m_K+i] = double(T2[i]) ;
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[18+m_K+i] = double(B2[i]) ;
m_workspace[21+m_K] = double(m_K) ;
// Get new frame
// precondition : this->m_approx[m_approxindex_FN]->approximate should have been called !
//const VEC3& N = this->m_approx[m_approxindex_FN]->getApprox(v1,m_attrindex_FN) ; // get new N
//assert(N == N2 || !"HalfEdgeSelector_Radiance<PFP>::computeSquaredLightfieldDifferenceNumerical: only works with half-collapse") ;
m_n[0] = N[0] ;
m_n[1] = N[1] ;
m_n[2] = N[2] ;
// call integral function
double integ, ar ;
VEC3 integral ;
VEC3 area ;
for (unsigned int channel = 0 ; channel < 3 ; ++channel)
{
m_workspace[3] = m_avgColor[v0][channel];
m_workspace[14+m_K] = m_avgColor[v1][channel];
for (unsigned int i = 0 ; i < m_K ; ++i)
{
m_workspace[11+i] = coefs1[i][channel] ;
m_workspace[22+m_K+i] = coefs2[i][channel] ;
}
m_integrator.Compute(&integ, &ar, &SquaredDifferenceOfCartesianFunctions, m_workspace, &isInDomain, m_n) ;
integral[channel] = integ/(4*M_PI) ;
area[channel] = ar ;
}
return integral ;
}
template <typename PFP>
VEC3 HalfEdgeSelector_Radiance<PFP>::computeSquaredLightfieldDifferenceAnalytical(const Dart& v0, const Dart& v1) const
{
// get two frames
const VEC3& T1 = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[v0] ;
const VEC3& T2 = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[v1] ;
const VEC3& B1 = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[v0] ;
const VEC3& B2 = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[v1] ;
const VEC3& N1 = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[v0] ;
const VEC3& N2 = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[v1] ;
// Get new frame
// precondition : this->m_approx[m_approxindex_FN]->approximate should have been called !
const VEC3& N = this->m_approx[m_approxindex_FN]->getApprox(v1,m_attrindex_FN) ; // get new N
assert(N == N2 || !"HalfEdgeSelector_Radiance<PFP>::computeSquaredLightfieldDifferenceAnalytical: only works with half-collapse") ;
// compute new frame
VEC3 T ;
if (N2 != N1)
T = N2 ^ N1 ; // i is perpendicular to newNormal
else
T = N1 ^ VEC3(1,2,3) ; // second random vector
T.normalize() ;
// Compute D1' and D2'
VEC3 B1prime = N1 ^ T ;
B1prime.normalize() ;
VEC3 B2prime = N2 ^ T ;
B2prime.normalize() ;
// Rotation dans sens trigo dans le plan tangent autour de N (T1 --> T)
const REAL gamma1 = ((B1 * T) > 0 ? 1 : -1) * acos( std::max(std::min(REAL(1), T1 * T ), -REAL(1))) ; // angle positif ssi
const REAL gamma2 = ((B2 * T) > 0 ? 1 : -1) * acos( std::max(std::min(REAL(1), T2 * T ), -REAL(1))) ; // -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(REAL(1), N * N1), -REAL(1)) ) ; // angle positif ssi
const REAL alpha2 = ((N * B2prime) > 0 ? -1 : 1) * acos( std::max(std::min(REAL(1), N * N2), -REAL(1)) ) ; // PI/2 < angle(j1',n) < -PI/2 ssi j1'*n < 0
double alpha = fabs(alpha1 + alpha2) ;
// get coefs of v1 and v2
std::vector<VEC3> coefs1, coefs2, coefs ;
coefs1.resize(m_K) ; coefs2.resize(m_K) ;
for (unsigned int i = 0 ; i < m_K ; ++i)
{
coefs1[i] = this->m_approx[m_approxindex_HF[i]]->getAttr(m_attrindex_HF[i])[v0] ;
coefs2[i] = this->m_approx[m_approxindex_HF[i]]->getAttr(m_attrindex_HF[i])[v1] ;
}
Utils::QuadricHF<REAL> q(coefs2, gamma2, alpha2) ;
bool opt = q.findOptimizedCoefs(coefs) ; // coefs of d1's lightfield rotated around new local axis
q += Utils::QuadricHF<REAL>(coefs1, gamma1, alpha1) ;
if (!opt)
{
std::cerr << "HalfEdgeSelector_Radiance<PFP>::Optimization failed (should never happen since no optim is done)" << std::endl ;
std::cout << alpha1 << std::endl ;
}
VEC3 avgColDiff = m_avgColor[v0] - m_avgColor[v1] ;
for (unsigned int i = 0 ; i < 3 ; ++i)
avgColDiff[i] = pow(avgColDiff[i],2) ;
const VEC3 LFerr = q.evalR3(coefs) ;
VEC3 err = (alpha / M_PI) * avgColDiff + LFerr ;
for (unsigned int i = 0 ; i < 3 ; ++i)
err[i] = std::max(REAL(0),err[i]) ;
return err ;
}
template <typename PFP>
VEC3 HalfEdgeSelector_Radiance<PFP>::computeGradient(const VEC3& P0, const VEC3& Pi, const VEC3& Pj, const Dart& v0, const Dart& v1, const Dart& vi, const Dart& vj, unsigned int channel) const
{
const VEC3 ei = P0 - Pj ;
const VEC3 ej = Pi - P0 ;
// domaine sur lequel intégrer : hémisphère défini par n
const VEC3& ni = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[vi] ;
const VEC3& nj = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[vj] ;
const VEC3& n0 = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[v0] ;
VEC3 normal = (ni+nj+n0)/3. ; normal.normalize() ;
for (unsigned int i = 0 ; i < 3 ; ++i)
m_n[i] = normal[i] ;
std::vector<Dart> vertices ; vertices.resize(3) ;
vertices[0] = vi ;
vertices[1] = vj ;
vertices[2] = v1 ;
std::vector<double> coefs ; coefs.resize(m_K) ;
VEC3 res ;
for (unsigned int i = 0 ; i < 3 ; ++i)
{
const Dart& d = vertices[i] ;
const VEC3& t = this->m_approx[m_approxindex_FT]->getAttr(m_attrindex_FT)[d] ;
const VEC3& b = this->m_approx[m_approxindex_FB]->getAttr(m_attrindex_FB)[d] ;
const VEC3& n_d = this->m_approx[m_approxindex_FN]->getAttr(m_attrindex_FN)[d] ;
double avgCol = m_avgColor[d][channel] ;
for (unsigned int k = 0 ; k < m_K ; ++k)
{
coefs[k] = this->m_approx[m_approxindex_HF[k]]->getAttr(m_attrindex_HF[k])[d][channel] ;
}
res[i] = computeIntegral(&avgCol, t, b, n_d, m_K, coefs) ;
}
const REAL& ci = res[0] ;
const REAL& cj = res[1] ;
const REAL& c1 = res[2] ;
return (ei.norm2()*fabs(ci-c1) + (ei*ej)*fabs(cj-c1))*ej - (ej.norm2()*fabs(cj-c1) + (ei*ej)*fabs(ci-c1))*ei ;
}
template <typename PFP>
VEC3 HalfEdgeSelector_Radiance<PFP>::integrateDscalGrad(
const VEC3& d, const unsigned int& K, const VEC3& N, const VEC3& ei, const VEC3& ej,
const VEC3* coefs1, const VEC3& T1, const VEC3& B1, const VEC3& N1, const VEC3& avg1,
const VEC3* coefsi, const VEC3& Ti, const VEC3& Bi, const VEC3& Ni, const VEC3& avgi,
const VEC3* coefsj, const VEC3& Tj, const VEC3& Bj, const VEC3& Nj, const VEC3& avgj
) const
{
// set domain of integration
for (unsigned int i = 0 ; i < 3 ; ++i)
m_n[i] = N[i] ;
// prepare workspace
for (unsigned int i = 0 ; i < 3 ; ++i) // d
m_workspace[i] = d[i] ;
// nbCoefs
m_workspace[3] = K ;
for (unsigned int i = 0 ; i < 3 ; ++i) // ei
m_workspace[4+i] = ei[i] ;
for (unsigned int i = 0 ; i < 3 ; ++i) // ej
m_workspace[7+i] = ej[i] ;
// avg1, avgi, avgj
m_workspace[11] = avg1[0] ;
m_workspace[12] = avgi[0] ;
m_workspace[13] = avgj[0] ;
m_workspace[14] = avg1[1] ;
m_workspace[15] = avgi[1] ;
m_workspace[16] = avgj[1] ;
m_workspace[17] = avg1[2] ;
m_workspace[18] = avgi[2] ;
m_workspace[19] = avgj[2] ;
// N1, Ni, Nj
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[20+i] = N1[i] ;
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[23+i] = Ni[i] ;
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[26+i] = Nj[i] ;
// T1, Ti, Tj
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[29+i] = T1[i] ;
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[32+i] = Ti[i] ;
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[35+i] = Tj[i] ;
// B1, Bi, Bj
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[38+i] = B1[i] ;
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[41+i] = Bi[i] ;
for (unsigned int i = 0 ; i < 3 ; ++i)
m_workspace[44+i] = Bj[i] ;
// coefs1, coefsi, coefsj
for (unsigned int c = 0 ; c < 3 ; ++c)
for (unsigned int k = 0 ; k < K ; ++k)
{
m_workspace[47+c*K+k] = coefs1[k][c] ;
m_workspace[47+(3+c)*K+k] = coefsi[k][c] ;
m_workspace[47+(6+c)*K+k] = coefsj[k][c] ;
}
VEC3 res ;
for (unsigned int c = 0 ; c < 3 ; ++c)
{
m_workspace[10] = c ;
// exec integrator
double integral, area ;
m_integrator.Compute(&integral, &area, &dispScalGrad, m_workspace, &isInDomain, m_n) ;
res[c] = integral/(4*M_PI) ;
}
Dart d2 = m.phi2(d);
return res ;
}
const VEC3& p0 = m_position[d2];
VEC3& n0 = m_normal[d2];
const SH& r0 = m_radiance[d2];
template <typename PFP>
VEC3 HalfEdgeSelector_Radiance<PFP>::integrateDlf(
const VEC3& d, const unsigned int& K, const VEC3& N, const VEC3& ei, const VEC3& ej,
const std::vector<VEC3*> coefs0, const VEC3& T0, const VEC3& B0, const VEC3& N0, const VEC3& avg0,
const std::vector<VEC3*> coefs1, const VEC3& T1, const VEC3& B1, const VEC3& N1, const VEC3& avg1,
const std::vector<VEC3*> coefsi, const VEC3& Ti, const VEC3& Bi, const VEC3& Ni, const VEC3& avgi,
const std::vector<VEC3*> coefsj, const VEC3& Tj, const VEC3& Bj, const VEC3& Nj, const VEC3& avgj
) const
{
// set domain of integration
for (unsigned int i = 0 ; i < 3 ; ++i)
m_n[i] = N[i] ;
// prepare workspace
for (unsigned int i = 0 ; i < 3 ; ++i) // d
m_workspace[i] = d[i] ;
// nbCoefs
m_workspace[3] = K ;
for (unsigned int i = 0 ; i < 3 ; ++i) // ei
m_workspace[4+i] = ei[i] ;
for (unsigned int i = 0 ; i < 3 ; ++i) // ej
m_workspace[7+i] = ej[i] ;
// avg1, avgi, avgj
m_workspace[12] = avg0[0] ;
m_workspace[13] = avg1[0] ;
m_workspace[14] = avgi[0] ;
m_workspace[15] = avgj[0] ;
m_workspace[16] = avg0[1] ;
m_workspace[17] = avg1[1] ;