Commit b90b69f1 authored by Sauvage's avatar Sauvage

collector one ring : ajout de l'aire et du tenseur Normal cycles

parent c8b7f315
......@@ -318,10 +318,11 @@ void computeCurvatureVertex_NormalCycles(
// collect the normal cycle tensor
Algo::Selection::Collector_WithinSphere<PFP> neigh(map, position, radius, thread) ;
// Algo::Selection::Collector_OneRing<PFP> neigh(map, thread) ;
neigh.collectAll(dart) ;
MATRIX tensor(0) ;
neigh.computeNormalCyclesTensor(edgeangle,tensor);
neigh.computeNormalCyclesTensor(position, edgeangle,tensor);
// solve eigen problem
Eigen::SelfAdjointEigenSolver<E_MATRIX> solver (Utils::convertRef<E_MATRIX>(tensor));
......
......@@ -113,11 +113,11 @@ public:
template <typename PPFP>
friend std::ostream& operator<<(std::ostream &out, const Collector<PPFP>& c);
REAL computeArea () {
REAL computeArea (const VertexAttribute<VEC3>& pos) {
assert(!"Warning: Collector<PFP>::computeArea() should be overloaded in non-virtual derived classes");
return 0.0;
}
void computeNormalCyclesTensor (const EdgeAttribute<REAL>&, typename PFP::MATRIX33&) {assert(!"Warning: Collector<PFP>::computeNormalCyclesTensor() should be overloaded in non-virtual derived classes"); }
void computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, const EdgeAttribute<REAL>& edgeangle, typename PFP::MATRIX33&) {assert(!"Warning: Collector<PFP>::computeNormalCyclesTensor() should be overloaded in non-virtual derived classes"); }
};
......@@ -135,11 +135,18 @@ public:
template <typename PFP>
class Collector_OneRing : public Collector<PFP>
{
protected:
typedef typename PFP::VEC3 VEC3 ;
typedef typename PFP::REAL REAL ;
public:
Collector_OneRing(typename PFP::MAP& m, unsigned int thread=0):
Collector<PFP>(m, thread) {}
void collectAll(Dart d);
void collectBorder(Dart d);
REAL computeArea(const VertexAttribute<VEC3>& pos);
void computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, const EdgeAttribute<REAL>&edgeangle, typename PFP::MATRIX33&);
};
/*********************************************************
......@@ -175,8 +182,8 @@ public:
void collectAll(Dart d);
void collectBorder(Dart d);
REAL computeArea();
void computeNormalCyclesTensor (const EdgeAttribute<REAL>&, typename PFP::MATRIX33&);
REAL computeArea(const VertexAttribute<VEC3>& pos);
void computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, const EdgeAttribute<REAL>& edgeangle, typename PFP::MATRIX33&);
};
/*********************************************************
......
......@@ -139,6 +139,35 @@ void Collector_OneRing<PFP>::collectBorder(Dart d)
this->border.push_back(this->map.phi1(it));
}
template <typename PFP>
typename PFP::REAL Collector_OneRing<PFP>::computeArea(const VertexAttribute<VEC3>& pos)
{
assert(this->isInsideCollected || !"computeArea: inside cells have not been collected.") ;
REAL area = 0;
for (std::vector<Dart>::const_iterator it = this->insideFaces.begin(); it != this->insideFaces.end(); ++it)
area += Algo::Geometry::triangleArea<PFP>(this->map, *it, pos);
return area;
}
template <typename PFP>
void Collector_OneRing<PFP>::computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, const EdgeAttribute<REAL>& edgeangle, typename PFP::MATRIX33& tensor){
assert(this->isInsideCollected || !"computeNormalCyclesTensor: inside cells have not been collected.") ;
tensor.zero() ;
// collect edges inside the neighborhood
for (std::vector<Dart>::const_iterator it = this->insideEdges.begin(); it != this->insideEdges.end(); ++it)
{
const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(this->map, *it, pos) ;
tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) ;
}
tensor /= computeArea(pos) ;
}
/*********************************************************
* Collector Within Sphere
*********************************************************/
......@@ -258,61 +287,61 @@ void Collector_WithinSphere<PFP>::collectBorder(Dart d)
template <typename PFP>
typename PFP::REAL Collector_WithinSphere<PFP>::computeArea()
typename PFP::REAL Collector_WithinSphere<PFP>::computeArea(const VertexAttribute<VEC3>& pos)
{
assert(this->isInsideCollected || !"computeArea: inside cells have not been collected.") ;
VEC3 centerPosition = position[this->centerDart];
VEC3 centerPosition = pos[this->centerDart];
REAL area = 0;
for (std::vector<Dart>::const_iterator it = this->insideFaces.begin(); it != this->insideFaces.end(); ++it)
area += Algo::Geometry::triangleArea<PFP>(this->map, *it, this->position);
area += Algo::Geometry::triangleArea<PFP>(this->map, *it, pos);
for (std::vector<Dart>::const_iterator it = this->border.begin(); it != this->border.end(); ++it)
{
const Dart f = this->map.phi1(*it); // we know that f is outside
const Dart g = this->map.phi1(f);
if (Geom::isPointInSphere(this->position[g], centerPosition, this->radius))
if (Geom::isPointInSphere(pos[g], centerPosition, this->radius))
{ // only f is outside
typename PFP::REAL alpha, beta;
Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, this->radius, *it, this->position, alpha);
Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, this->radius, this->map.phi2(f), this->position, beta);
area += (alpha+beta - alpha*beta) * Algo::Geometry::triangleArea<PFP>(this->map, *it, this->position);
Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, this->radius, *it, pos, alpha);
Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, this->radius, this->map.phi2(f), pos, beta);
area += (alpha+beta - alpha*beta) * Algo::Geometry::triangleArea<PFP>(this->map, *it, pos);
}
else
{ // f and g are outside
typename PFP::REAL alpha, beta;
Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, this->radius, *it, this->position, alpha);
Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, this->radius, this->map.phi2(g), this->position, beta);
area += alpha * beta * Algo::Geometry::triangleArea<PFP>(this->map, *it, this->position);
Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, this->radius, *it, pos, alpha);
Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, this->radius, this->map.phi2(g), pos, beta);
area += alpha * beta * Algo::Geometry::triangleArea<PFP>(this->map, *it, pos);
}
}
return area;
}
template <typename PFP>
void Collector_WithinSphere<PFP>::computeNormalCyclesTensor (const EdgeAttribute<REAL>& edgeangle, typename PFP::MATRIX33& tensor){
void Collector_WithinSphere<PFP>::computeNormalCyclesTensor (const VertexAttribute<VEC3>& pos, const EdgeAttribute<REAL>& edgeangle, typename PFP::MATRIX33& tensor){
assert(this->isInsideCollected || !"computeNormalCyclesTensor: inside cells have not been collected.") ;
VEC3 centerPosition = position[this->centerDart];
VEC3 centerPosition = pos[this->centerDart];
tensor.zero() ;
// collect edges inside the neighborhood
for (std::vector<Dart>::const_iterator it = this->insideEdges.begin(); it != this->insideEdges.end(); ++it)
{
const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(this->map, *it, position) ;
const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(this->map, *it, pos) ;
tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) ;
}
// collect edges crossing the neighborhood's border
for (std::vector<Dart>::const_iterator it = this->border.begin(); it != this->border.end(); ++it)
{
const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(this->map, *it, position) ;
const VEC3 e = Algo::Geometry::vectorOutOfDart<PFP>(this->map, *it, pos) ;
REAL alpha ;
Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, radius, *it, position, alpha) ;
Algo::Geometry::intersectionSphereEdge<PFP>(this->map, centerPosition, radius, *it, pos, alpha) ;
tensor += Geom::transposed_vectors_mult(e,e) * edgeangle[*it] * (1 / e.norm()) * alpha ;
}
tensor /= computeArea() ;
tensor /= computeArea(pos) ;
}
......
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