Commit 32f4ff4e authored by Pierre Kraemer's avatar Pierre Kraemer

Merge branch 'develop' into 'develop'

Develop

See merge request !43
parents 4dc1bdba 79056382
......@@ -70,8 +70,8 @@ int main(int argc, char **argv)
std::cout << " NB Faces "<< Algo::Topo::getNbOrbits<FACE>(myMap) << std::endl;
std::cout << " NB Vertices "<< nbVertices << std::endl;
std::vector<VertexAttribute<typename PFP::VEC3, MAP> *> attr;
attr.push_back(&position);
std::vector<VertexAttribute<typename PFP::VEC3, MAP> > attr;
attr.push_back(position);
Algo::Surface::Decimation::decimate<PFP>(myMap, Algo::Surface::Decimation::S_QEM, Algo::Surface::Decimation::A_QEM, attr, nbVertices * 0.05) ;
VertexAttribute<PFP::VEC3, MAP> normal = myMap.addAttribute<PFP::VEC3,VERTEX,MAP>( "normal") ;
......
......@@ -52,8 +52,8 @@ int main(int argc, char **argv)
Algo::Surface::Modelisation::LoopSubdivision<PFP>(myMap, position) ;
unsigned int nbVertices = Algo::Topo::getNbOrbits<VERTEX>(myMap) ;
std::vector<VertexAttribute<typename PFP::VEC3, MAP> *> attr;
attr.push_back(&position);
std::vector<VertexAttribute<typename PFP::VEC3, MAP> > attr;
attr.push_back(position);
Algo::Surface::Decimation::decimate<PFP>(myMap, Algo::Surface::Decimation::S_QEM, Algo::Surface::Decimation::A_QEM, attr, nbVertices * 0.1) ;
Algo::Surface::Modelisation::LoopSubdivision<PFP>(myMap, position) ;
......
......@@ -21,8 +21,35 @@
</property>
</widget>
</item>
<item row="0" column="1">
<widget class="QComboBox" name="combo_positionVBO">
<item row="5" column="1">
<spacer name="verticalSpacer">
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>20</width>
<height>40</height>
</size>
</property>
</spacer>
</item>
<item row="4" column="0" colspan="2">
<widget class="QPushButton" name="button_decimate">
<property name="text">
<string>Decimate</string>
</property>
</widget>
</item>
<item row="1" column="0">
<widget class="QLabel" name="label_3">
<property name="text">
<string>Normal :</string>
</property>
</widget>
</item>
<item row="1" column="1">
<widget class="QComboBox" name="combo_normalVBO">
<property name="sizePolicy">
<sizepolicy hsizetype="Expanding" vsizetype="Fixed">
<horstretch>0</horstretch>
......@@ -36,15 +63,21 @@
</item>
</widget>
</item>
<item row="1" column="0">
<widget class="QLabel" name="label_3">
<property name="text">
<string>Normal :</string>
<item row="2" column="0" colspan="2">
<widget class="QSlider" name="slider_decimationGoal">
<property name="maximum">
<number>100</number>
</property>
<property name="value">
<number>50</number>
</property>
<property name="orientation">
<enum>Qt::Horizontal</enum>
</property>
</widget>
</item>
<item row="1" column="1">
<widget class="QComboBox" name="combo_normalVBO">
<item row="0" column="1">
<widget class="QComboBox" name="combo_positionVBO">
<property name="sizePolicy">
<sizepolicy hsizetype="Expanding" vsizetype="Fixed">
<horstretch>0</horstretch>
......@@ -58,18 +91,12 @@
</item>
</widget>
</item>
<item row="2" column="1">
<spacer name="verticalSpacer">
<property name="orientation">
<enum>Qt::Vertical</enum>
</property>
<property name="sizeHint" stdset="0">
<size>
<width>20</width>
<height>40</height>
</size>
<item row="3" column="0" colspan="2">
<widget class="QCheckBox" name="checkbox_halfCollapse">
<property name="text">
<string>Half Collapse</string>
</property>
</spacer>
</widget>
</item>
</layout>
</widget>
......
#ifndef SPHERICALFUNCTIONINTEGRATORCARTESIAN_H
#define SPHERICALFUNCTIONINTEGRATORCARTESIAN_H
#include "sphere_lebedev_rule.h"
typedef double (*CartesianFunction)(double x, double y, double z, void* userData); // Prototype for the function to be evaluated
typedef bool (*CartesianDomain)(double x, double y, double z, void* userData); // Prototype for the domain definition (true: inside, false: outside)
class SphericalFunctionIntegratorCartesian
{
public:
SphericalFunctionIntegratorCartesian();
~SphericalFunctionIntegratorCartesian();
static const unsigned int maxRuleId = 65; // Span of rule id (inclusive)
static inline bool RuleAvailable(unsigned int ruleId); // States that the quadrature rule is available
static inline unsigned int RuleOrder(unsigned int ruleId); // Number of points used in the quadrature
static inline unsigned int RulePrecision(unsigned int ruleId); // Max degree of exactly integrated polynomial
void Init(unsigned int ruleId); // Rule to use for subsequent integration - allocates quadrature samples and weights
void Release(); // Release cached informations
/*
* Integrates a function over a user-specified domain of the full 2D-sphere
*
* outIntegral: resulting value
* outArea: area of the integration domain, as specified by the dom function
* f: function to integrate
* userDataFunction: user callback value passed to f during evaluation
* dom: Implicit domain definition (true when inside, false outside)
* userDataDomain: user callback value passed to dom during evaluation
*
*/
void Compute(double* outIntegral, double* outArea, CartesianFunction f, void* userDataFunction, CartesianDomain dom, void* userDataDomain) const;
protected:
unsigned int rId;
unsigned int rOrder;
double* quadValues; // [x_i] then [y_i], then [z_i], then [w_i] values
};
bool SphericalFunctionIntegratorCartesian::RuleAvailable(unsigned int ruleId)
{
return available_table(ruleId) == 1;
}
unsigned int SphericalFunctionIntegratorCartesian::RuleOrder(unsigned int ruleId)
{
return order_table(ruleId);
}
unsigned int SphericalFunctionIntegratorCartesian::RulePrecision(unsigned int ruleId)
{
return precision_table(ruleId);
}
#endif
#ifndef _HALFEDGESELECTOR_RADIANCE_H_
#define _HALFEDGESELECTOR_RADIANCE_H_
#include "Algo/Decimation/selector.h"
#include "Algo/Decimation/approximator.h"
#include "Utils/qem.h"
#include "Utils/sphericalHarmonics.h"
#include "SphericalFunctionIntegratorCartesian.h"
namespace CGoGN
{
namespace SCHNApps
{
template <typename PFP>
class HalfEdgeSelector_Radiance : public Algo::Surface::Decimation::Selector<PFP>
{
public:
typedef typename PFP::MAP MAP ;
typedef typename PFP::REAL REAL ;
typedef typename PFP::VEC3 VEC3 ;
typedef typename Utils::SphericalHarmonics<REAL, VEC3> SH;
private:
VertexAttribute<VEC3, MAP>& m_position;
VertexAttribute<VEC3, MAP>& m_normal;
VertexAttribute<SH, MAP>& m_radiance;
Algo::Surface::Decimation::Approximator<PFP, VEC3, DART>& m_positionApproximator;
Algo::Surface::Decimation::Approximator<PFP, VEC3, DART>& m_normalApproximator;
Algo::Surface::Decimation::Approximator<PFP, SH, DART>& m_radianceApproximator;
typedef struct
{
typename std::multimap<float, Dart>::iterator it;
bool valid;
static std::string CGoGNnameOfType() { return "LightfieldGradHalfEdgeInfo"; }
} LightfieldHalfEdgeInfo;
typedef NoTypeNameAttribute<LightfieldHalfEdgeInfo> HalfEdgeInfo;
DartAttribute<HalfEdgeInfo, PFP2::MAP> halfEdgeInfo;
VertexAttribute<Utils::Quadric<PFP2::REAL>, PFP2::MAP> m_quadric;
// VertexAttribute<PFP2::VEC3, PFP2::MAP> m_avgColor;
unsigned int m_nb_coefs;
SphericalFunctionIntegratorCartesian m_integrator;
std::multimap<float, Dart> halfEdges;
typename std::multimap<float, Dart>::iterator cur;
void initHalfEdgeInfo(Dart d);
void updateHalfEdgeInfo(Dart d, bool recompute);
void computeHalfEdgeInfo(Dart d, HalfEdgeInfo& einfo);
void recomputeQuadric(const Dart d);
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,
VertexAttribute<VEC3, MAP>& pos,
VertexAttribute<VEC3, MAP>& norm,
VertexAttribute<SH, MAP>& rad,
Algo::Surface::Decimation::Approximator<PFP, VEC3, DART>& posApprox,
Algo::Surface::Decimation::Approximator<PFP, VEC3, DART>& normApprox,
Algo::Surface::Decimation::Approximator<PFP, SH, DART>& radApprox
) :
Algo::Surface::Decimation::Selector<PFP2>(m),
m_position(pos),
m_normal(norm),
m_radiance(rad),
m_positionApproximator(posApprox),
m_normalApproximator(normApprox),
m_radianceApproximator(radApprox),
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");
}
~HalfEdgeSelector_Radiance()
{
this->m_map.removeAttribute(halfEdgeInfo);
this->m_map.removeAttribute(m_quadric);
// this->m_map.removeAttribute(m_avgColor);
m_integrator.Release();
}
Algo::Surface::Decimation::SelectorType getType() { return Algo::Surface::Decimation::S_OTHER; }
bool init();
bool nextEdge(Dart& d) const;
void updateBeforeCollapse(Dart d);
void updateAfterCollapse(Dart d2, Dart dd2);
void updateWithoutCollapse() { }
void getEdgeErrors(EdgeAttribute<PFP2::REAL, PFP2::MAP> *errors) const
{
assert(errors != NULL || !"HalfEdgeSelector_Radiance::getEdgeErrors requires non null vertexattribute argument") ;
if (!errors->isValid())
std::cerr << "HalfEdgeSelector_Radiance::getEdgeErrors requires valid edgeattribute argument" << std::endl ;
assert(halfEdgeInfo.isValid());
TraversorE<PFP2::MAP> travE(this->m_map);
for(Dart d = travE.begin() ; d != travE.end() ; d = travE.next())
{
Dart dd = this->m_map.phi2(d);
if (halfEdgeInfo[d].valid)
{
(*errors)[d] = halfEdgeInfo[d].it->first;
}
if (halfEdgeInfo[dd].valid && halfEdgeInfo[dd].it->first < (*errors)[d])
{
(*errors)[d] = halfEdgeInfo[dd].it->first;
}
if (!(halfEdgeInfo[d].valid || halfEdgeInfo[dd].valid))
(*errors)[d] = -1;
}
}
};
} // namespace SCHNApps
} // namespace CGoGN
#include "halfEdgeSelectorRadiance.hpp"
#endif // _HALFEDGESELECTOR_RADIANCE_H_
namespace CGoGN
{
namespace SCHNApps
{
template <typename PFP>
bool HalfEdgeSelector_Radiance<PFP>::init()
{
MAP& m = this->m_map ;
assert(m_positionApproximator.getType() == Algo::Surface::Decimation::A_hQEM
|| m_positionApproximator.getType() == Algo::Surface::Decimation::A_hHalfCollapse
|| !"Approximator for selector (HalfEdgeSelector_Radiance) must be of a half-edge approximator") ;
assert(m_normalApproximator.getType() == Algo::Surface::Decimation::A_hQEM
|| m_normalApproximator.getType() == Algo::Surface::Decimation::A_hHalfCollapse
|| !"Approximator for selector (HalfEdgeSelector_Radiance) must be of a half-edge approximator") ;
assert(m_radianceApproximator.getType() == Algo::Surface::Decimation::A_hQEM
|| m_radianceApproximator.getType() == Algo::Surface::Decimation::A_hHalfCollapse
|| !"Approximator for selector (HalfEdgeSelector_Radiance) must be of a half-edge approximator") ;
if(m_positionApproximator.getApproximatedAttributeName() != m_position.name())
{
return false;
}
if(m_normalApproximator.getApproximatedAttributeName() != m_normal.name())
{
return false;
}
if(m_radianceApproximator.getApproximatedAttributeName() != m_radiance.name())
{
return false;
}
m_nb_coefs = SH::get_nb_coefs();
m_integrator.Init(29) ;
// init QEM quadrics
for (Vertex v : allVerticesOf(m))
{
Utils::Quadric<REAL> q ; // create one quadric
m_quadric[v] = q ; // per vertex
}
for (Face f : allFacesOf(m))
{
Dart d = f.dart;
Dart d1 = m.phi1(d) ; // for each triangle,
Dart d_1 = m.phi_1(d) ; // initialize the quadric of the triangle
Utils::Quadric<REAL> q(m_position[d], m_position[d1], m_position[d_1]) ;
m_quadric[d] += q ; // and add the contribution of
m_quadric[d1] += q ; // this quadric to the ones
m_quadric[d_1] += q ; // of the 3 incident vertices
}
// init multimap for each Half-edge
halfEdges.clear() ;
for(Dart d = m.begin() ; d != m.end() ; m.next(d))
{
initHalfEdgeInfo(d) ; // init the edges with their optimal position
} // and insert them in the multimap according to their error
cur = halfEdges.begin() ; // init the current edge to the first one
return true ;
}
template <typename PFP>
bool HalfEdgeSelector_Radiance<PFP>::nextEdge(Dart& d) const
{
if(cur == halfEdges.end() || halfEdges.empty())
return false ;
d = (*cur).second ;
return true ;
}
template <typename PFP>
void HalfEdgeSelector_Radiance<PFP>::updateBeforeCollapse(Dart d)
{
MAP& m = this->m_map ;
Dart dd = m.phi2(d);
HalfEdgeInfo& edgeE = halfEdgeInfo[d];
if(edgeE.valid)
{
edgeE.valid = false ;
halfEdges.erase(edgeE.it) ;
}
edgeE = halfEdgeInfo[dd];
if(edgeE.valid)
{
edgeE.valid = false ;
halfEdges.erase(edgeE.it) ;
}
edgeE = halfEdgeInfo[m.phi1(d)];
if(edgeE.valid)
{
edgeE.valid = false ;
halfEdges.erase(edgeE.it) ;
}
edgeE = halfEdgeInfo[m.phi_1(d)];
if(edgeE.valid)
{
edgeE.valid = false ;
halfEdges.erase(edgeE.it) ;
}
edgeE = halfEdgeInfo[m.phi1(dd)];
if(edgeE.valid)
{
edgeE.valid = false ;
halfEdges.erase(edgeE.it) ;
}
edgeE = halfEdgeInfo[m.phi_1(dd)];
if(edgeE.valid)
{
edgeE.valid = false ;
halfEdges.erase(edgeE.it) ;
}
}
template <typename PFP>
void HalfEdgeSelector_Radiance<PFP>::recomputeQuadric(const Dart d)
{
MAP& m = this->m_map;
Utils::Quadric<REAL>& q = m_quadric[d];
q.zero() ;
Traversor2VF<MAP> tf(m, d);
for (Dart f = tf.begin() ; f != tf.end() ; f = tf.next())
{
q += Utils::Quadric<REAL>(m_position[f], m_position[m.phi1(f)], m_position[m.phi_1(f)]) ;
}
}
template <typename PFP>
void HalfEdgeSelector_Radiance<PFP>::updateAfterCollapse(Dart d2, Dart dd2)
{
MAP& m = this->m_map ;
Dart stop = m.phi2(m.phi_1(d2));
recomputeQuadric(d2);
Dart it = dd2;
do
{
recomputeQuadric(m.phi1(it));
it = m.phi2(m.phi_1(it));
} while (it != stop);
DartMarkerStore<MAP> dm(m);
Traversor2VVaE<MAP> tv(m, d2);
for (Dart v = tv.begin() ; v != tv.end() ; v = tv.next())
{
dm.mark(v);
updateHalfEdgeInfo(v, true);
}
it = dd2;
do
{
Traversor2VVaE<MAP> tv2(m, m.phi1(it));
for (Dart v = tv2.begin() ; v != tv2.end() ; v = tv2.next())
{
dm.mark(v);
updateHalfEdgeInfo(v, true);
}
it = m.phi2(m.phi_1(it));
} while (it != stop);
Traversor2VE<MAP> te(m, d2);
for (Dart v = te.begin() ; v != te.end() ; v = te.next())
{
if (!dm.isMarked(v))
updateHalfEdgeInfo(v, false);
}
it = dd2;
do
{
Traversor2VE<MAP> te2(m, m.phi1(it));
for (Dart v = te2.begin() ; v != te2.end() ; v = te2.next())
{
if (!dm.isMarked(v))
updateHalfEdgeInfo(v, false);
}
it = m.phi2(m.phi_1(it));
} while (it != stop);
cur = halfEdges.begin() ; // set the current edge to the first one
}
template <typename PFP>
void HalfEdgeSelector_Radiance<PFP>::initHalfEdgeInfo(Dart d)
{
MAP& m = this->m_map ;
HalfEdgeInfo heinfo ;
if(m.edgeCanCollapse(d))
computeHalfEdgeInfo(d, heinfo) ;
else
heinfo.valid = false ;
halfEdgeInfo[d] = heinfo ;
}
template <typename PFP>
void HalfEdgeSelector_Radiance<PFP>::updateHalfEdgeInfo(Dart d, bool recompute)
{
MAP& m = this->m_map ;
HalfEdgeInfo& heinfo = halfEdgeInfo[d] ;
if(recompute)
{
if (heinfo.valid)
halfEdges.erase(heinfo.it);
if (m.edgeCanCollapse(d))
computeHalfEdgeInfo(d, heinfo);
else
heinfo.valid = false;
}
else
{
if (m.edgeCanCollapse(d))
{ // if the edge can be collapsed now
if (!heinfo.valid) // but it was not before
computeHalfEdgeInfo(d, heinfo) ;
}
else
{ // if the edge cannot be collapsed now
if (heinfo.valid) // and it was before
{
halfEdges.erase(heinfo.it) ;
heinfo.valid = false ;
}
}
}
}
template <typename PFP>
void HalfEdgeSelector_Radiance<PFP>::computeHalfEdgeInfo(Dart d, HalfEdgeInfo& heinfo)
{
MAP& m = this->m_map ;
Dart dd = m.phi1(d) ;
// compute all approximated attributes
m_positionApproximator.approximate(d) ;
m_normalApproximator.approximate(d) ;
m_radianceApproximator.approximate(d) ;
// Get all approximated attributes
// New position
const VEC3& newPos = m_positionApproximator.getApprox(d) ;
const VEC3& newN = m_normalApproximator.getApprox(d) ;
const SH& newRad = m_radianceApproximator.getApprox(d) ;
// WARNING : in the following, we consider the half-edge contraction v0 --> v1
const Dart& v0 = dd ;
const Dart& v1 = d ;
assert(newPos == m_position[v1]) ;
// Compute errors
Utils::Quadric<REAL> quadGeom ;
quadGeom += m_quadric[v0] ; // compute the sum of the
quadGeom += m_quadric[v1] ; // two vertices quadrics
const REAL geom = quadGeom(newPos) ;
const REAL rad = computeRadianceError(d);
const REAL t = 0.01 ;
const REAL err = t*geom + (1-t)*rad ;
// Check if errated values appear
if (err < -1e-6)
heinfo.valid = false ;
else
{
heinfo.it = this->halfEdges.insert(std::make_pair(std::max(err, REAL(0)), d)) ;
heinfo.valid = true ;
}
}
template <typename PFP>
typename PFP::REAL HalfEdgeSelector_Radiance<PFP>::computeRadianceError(Dart d)
{
MAP& m = this->m_map ;
Dart d2 = m.phi2(d);
const VEC3& p0 = m_position[d2];
VEC3& n0 = m_normal[d2];
const SH& r0 = m_radiance[d2];
const VEC3& p1 = m_position[d];
const SH& r1 = m_radiance[d];
REAL error(0);
Dart it1 = m.phi2(m.phi_1(d2));
Dart it2 = m.phi2(m.phi_1(it1));
while (it2 != d2)
{
const VEC3& p2 = m_position[m.phi1(it1)];
const SH& r2 = m_radiance[m.phi1(it1)];
const VEC3& p3 = m_position[m.phi1(it2)];
const SH& r3 = m_radiance[m.phi1(it2)];
VEC3 e = p0 - p1;
VEC3 e2 = p2 - p1;
VEC3 e3 = p3 - p1;
REAL tArea = Geom::triangleArea(p1, p2, p3);
REAL alpha2 = ((e3 * e3) * (e2 * e) - (e2 * e3) * (e3 * e)) / (4. * tArea * tArea) ;
REAL alpha3 = ((e2 * e2) * (e3 * e) - (e2 * e3) * (e2 * e)) / (4. * tArea * tArea) ;
REAL alpha1 = 1. - alpha2 - alpha3;
SH diffRad(r0);
diffRad -= (r1 * alpha1 + r2 * alpha2 + r3 * alpha3);
double integral;
double area;
m_integrator.Compute(&integral, &area, HalfEdgeSelector_Radiance<PFP>::SHEvalCartesian_Error, &diffRad, HalfEdgeSelector_Radiance<PFP>::isInHemisphere, n0.data());
error += tArea * integral / area;
it1 = it2;
it2 = m.phi2(m.phi_1(it1));
}
return error;
}
} // namespace SCHNApps
} // namespace CGoGN
#ifndef _MESHTABLESURFACE_RADIANCE_H_
#define _MESHTABLESURFACE_RADIANCE_H_
#include "Algo/Import/import.h"
#include "Algo/Import/import2tables.h"
namespace CGoGN
{
namespace SCHNApps
{