Commit 2ce8c117 authored by Pierre Kraemer's avatar Pierre Kraemer

update neighborhood area computation in normal cycle curvature algorithm

parent ca12b904
......@@ -74,9 +74,15 @@ typename PFP::REAL vertexBarycentricArea(typename PFP::MAP& map, Vertex v, const
template <typename PFP>
typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Vertex v, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position) ;
template <typename PFP>
typename PFP::REAL edgeArea(typename PFP::MAP& map, Edge e, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position);
template <typename PFP>
void computeAreaFaces(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, FaceAttribute<typename PFP::REAL, typename PFP::MAP>& face_area) ;
template <typename PFP>
void computeAreaEdges(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edge_area);
template <typename PFP>
void computeOneRingAreaVertices(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, VertexAttribute<typename PFP::REAL, typename PFP::MAP>& vertex_area) ;
......@@ -88,9 +94,16 @@ void computeVoronoiAreaVertices(typename PFP::MAP& map, const VertexAttribute<ty
namespace Parallel
{
template <typename PFP>
typename PFP::REAL totalArea(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position) ;
template <typename PFP>
void computeAreaFaces(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, FaceAttribute<typename PFP::REAL, typename PFP::MAP>& area);
template <typename PFP>
void computeAreaEdges(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& area);
template <typename PFP>
void computeOneRingAreaVertices(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, VertexAttribute<typename PFP::REAL, typename PFP::MAP>& area);
......@@ -99,8 +112,10 @@ void computeBarycentricAreaVertices(typename PFP::MAP& map, const VertexAttribut
template <typename PFP>
void computeVoronoiAreaVertices(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, VertexAttribute<typename PFP::REAL, typename PFP::MAP>& area);
} // namespace Parallel
} // namespace Geometry
} // namespace Surface
......
......@@ -74,15 +74,19 @@ typename PFP::REAL convexFaceArea(typename PFP::MAP& map, Face d, const VertexAt
template <typename PFP>
typename PFP::REAL totalArea(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position)
{
typename PFP::REAL area(0) ;
if (CGoGN::Parallel::NumberOfThreads > 1)
{
return Parallel::totalArea<PFP>(map, position);
}
typename PFP::REAL area(0);
foreach_cell<FACE>(map, [&] (Face f)
{
area += convexFaceArea<PFP>(map, f, position);
}
,AUTO);
return area ;
return area;
}
template <typename PFP>
......@@ -137,12 +141,23 @@ typename PFP::REAL vertexVoronoiArea(typename PFP::MAP& map, Vertex v, const Ver
}
template <typename PFP>
void computeAreaFaces(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, FaceAttribute<typename PFP::REAL, typename PFP::MAP>& face_area)
typename PFP::REAL edgeArea(typename PFP::MAP& map, Edge e, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position)
{
typename PFP::REAL area(0) ;
foreach_incident2<FACE>(map, e, [&] (Face f)
{
area += convexFaceArea<PFP>(map, f, position) / map.faceDegree(f) ;
});
return area ;
}
template <typename PFP>
void computeAreaFaces(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, FaceAttribute<typename PFP::REAL, typename PFP::MAP>& face_area)
{
if (CGoGN::Parallel::NumberOfThreads > 1)
{
Parallel::computeAreaFaces<PFP>(map,position,face_area);
Parallel::computeAreaFaces<PFP>(map, position, face_area);
return;
}
......@@ -151,7 +166,22 @@ void computeAreaFaces(typename PFP::MAP& map, const VertexAttribute<typename PFP
face_area[f] = convexFaceArea<PFP>(map, f, position) ;
}
,AUTO);
}
template <typename PFP>
void computeAreaEdges(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edge_area)
{
if (CGoGN::Parallel::NumberOfThreads > 1)
{
Parallel::computeAreaEdges<PFP>(map, position, edge_area);
return;
}
foreach_cell<EDGE>(map, [&] (Edge e)
{
edge_area[e] = edgeArea<PFP>(map, e, position) ;
}
,AUTO);
}
template <typename PFP>
......@@ -215,6 +245,24 @@ void computeVoronoiAreaVertices(typename PFP::MAP& map, const VertexAttribute<ty
namespace Parallel
{
template <typename PFP>
typename PFP::REAL totalArea(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position)
{
std::vector<typename PFP::REAL> areas;
areas.resize(CGoGN::Parallel::NumberOfThreads);
CGoGN::Parallel::foreach_cell<FACE>(map, [&] (Face f, unsigned int thr)
{
areas[thr-1] += convexFaceArea<PFP>(map, f, position);
}
,AUTO);
typename PFP::REAL area(0);
for (typename PFP::REAL a : areas)
area += a;
return area;
}
template <typename PFP>
void computeAreaFaces(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, FaceAttribute<typename PFP::REAL, typename PFP::MAP>& area)
{
......@@ -239,6 +287,14 @@ void computeAreaFaces(typename PFP::MAP& map, const VertexAttribute<typename PFP
});
}
template <typename PFP>
void computeAreaEdges(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& area)
{
CGoGN::Parallel::foreach_cell<EDGE>(map, [&] (Edge e, unsigned int /*thr*/)
{
area[e] = edgeArea<PFP>(map, e, position) ;
});
}
template <typename PFP>
void computeOneRingAreaVertices(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, VertexAttribute<typename PFP::REAL, typename PFP::MAP>& area)
......
......@@ -57,6 +57,58 @@ inline typename PFP::REAL edgeLength(typename PFP::MAP& map, Dart d, const Verte
return v.norm() ;
}
namespace Parallel
{
template <typename PFP>
typename PFP::REAL meanEdgeLength(typename PFP::MAP& map, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position)
{
std::vector<typename PFP::REAL> lengths;
std::vector<unsigned int> nbedges;
lengths.resize(CGoGN::Parallel::NumberOfThreads);
nbedges.resize(CGoGN::Parallel::NumberOfThreads);
CGoGN::Parallel::foreach_cell<EDGE>(map, [&] (Edge e, unsigned int thr)
{
lengths[thr-1] += edgeLength<PFP>(map, e, position);
++nbedges[thr-1];
}
,AUTO);
typename PFP::REAL length(0);
unsigned int nbe = 0;
for (unsigned int i = 0; i < CGoGN::Parallel::NumberOfThreads; ++i)
{
length += lengths[i];
nbe += nbedges[i];
}
return length / nbe;
}
} // namespace Parallel
template <typename PFP>
inline typename PFP::REAL meanEdgeLength(typename PFP::MAP& map,const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position)
{
if (CGoGN::Parallel::NumberOfThreads > 1)
{
return Parallel::meanEdgeLength<PFP>(map, position);
}
typename PFP::REAL length(0);
unsigned int nbe = 0;
foreach_cell<EDGE>(map, [&] (Edge e)
{
length += edgeLength<PFP>(map, e, position);
++nbe;
}
,AUTO);
return length / nbe;
}
template <typename PFP>
inline float angle(typename PFP::MAP& map, Dart d1, Dart d2, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position)
{
......
......@@ -113,6 +113,7 @@ void computeCurvatureVertices_NormalCycles(
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
......@@ -127,6 +128,7 @@ void computeCurvatureVertex_NormalCycles(
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
......@@ -160,6 +162,7 @@ void computeCurvatureVertices_NormalCycles_Projected(
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
......@@ -174,6 +177,7 @@ void computeCurvatureVertex_NormalCycles_Projected(
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
......@@ -189,6 +193,7 @@ void computeCurvatureVertices_NormalCycles(
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
......@@ -202,6 +207,7 @@ void computeCurvatureVertex_NormalCycles(
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
......@@ -215,6 +221,7 @@ void computeCurvatureVertices_NormalCycles_Projected(
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
......@@ -228,6 +235,7 @@ void computeCurvatureVertex_NormalCycles_Projected(
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
......@@ -245,6 +253,7 @@ void computeCurvatureVertices_NormalCycles(
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
......@@ -258,6 +267,7 @@ void computeCurvatureVertices_NormalCycles_Projected(
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position,
const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& normal,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgeangle,
const EdgeAttribute<typename PFP::REAL, typename PFP::MAP>& edgearea,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmax,
VertexAttribute<typename PFP::REAL, typename PFP::MAP>& kmin,
VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& Kmax,
......
This diff is collapsed.
......@@ -83,7 +83,7 @@ bool areTrianglesInIntersection(typename PFP::MAP& map, Face t1, Face t2, const
* et v1 et v2 les deux sommets de l'arête
*/
template <typename PFP>
bool intersectionSphereEdge(typename PFP::MAP& map, typename PFP::VEC3& center, typename PFP::REAL radius, Edge e, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, typename PFP::REAL& alpha) ;
bool intersectionSphereEdge(typename PFP::MAP& map, const typename PFP::VEC3& center, typename PFP::REAL radius, Edge e, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, typename PFP::REAL& alpha) ;
} // namespace Geometry
......
......@@ -238,7 +238,7 @@ bool areTrianglesInIntersection(typename PFP::MAP& map, Face t1, Face t2, const
}
template <typename PFP>
bool intersectionSphereEdge(typename PFP::MAP& map, typename PFP::VEC3& center, typename PFP::REAL radius, Edge e, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, typename PFP::REAL& alpha)
bool intersectionSphereEdge(typename PFP::MAP& map, const typename PFP::VEC3& center, typename PFP::REAL radius, Edge e, const VertexAttribute<typename PFP::VEC3, typename PFP::MAP>& position, typename PFP::REAL& alpha)
{
typedef typename PFP::VEC3 VEC3 ;
typedef typename PFP::REAL REAL ;
......
......@@ -122,16 +122,12 @@ public:
template <typename PPFP>
friend std::ostream& operator<<(std::ostream &out, const Collector<PPFP>& c);
virtual REAL computeArea (const VertexAttribute<VEC3, MAP>& /*pos*/)
virtual REAL computeArea (const VertexAttribute<VEC3, MAP>& /*pos*/, const EdgeAttribute<REAL, MAP>& /*edgearea*/)
{
assert(!"Warning: Collector<PFP>::computeArea() should be overloaded in non-virtual derived classes");
return 0.0;
}
virtual void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& /*pos*/, const EdgeAttribute<REAL, MAP>& /*edgeangle*/, typename PFP::MATRIX33&)
{
assert(!"Warning: Collector<PFP>::computeNormalCyclesTensor() should be overloaded in non-virtual derived classes");
}
virtual void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& /*pos*/, typename PFP::MATRIX33&)
virtual void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& /*pos*/, const EdgeAttribute<REAL, MAP>& /*edgeangle*/, const EdgeAttribute<REAL, MAP>& /*edgearea*/, typename PFP::MATRIX33&)
{
assert(!"Warning: Collector<PFP>::computeNormalCyclesTensor() should be overloaded in non-virtual derived classes");
}
......@@ -162,9 +158,8 @@ public:
void collectAll(Dart d);
void collectBorder(Dart d);
REAL computeArea(const VertexAttribute<VEC3, MAP>& pos);
void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& pos, const EdgeAttribute<REAL, MAP>& edgeangle, typename PFP::MATRIX33&);
void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& pos, typename PFP::MATRIX33&);
REAL computeArea(const VertexAttribute<VEC3, MAP>& pos, const EdgeAttribute<REAL, MAP>& edgearea);
void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& pos, const EdgeAttribute<REAL, MAP>& edgeangle, const EdgeAttribute<REAL, MAP>& edgearea, typename PFP::MATRIX33&);
};
/*********************************************************
......@@ -194,10 +189,8 @@ public:
void collectAll(Dart d);
void collectBorder(Dart d);
REAL computeArea(const VertexAttribute<VEC3, MAP>& pos);
void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& pos, const EdgeAttribute<REAL, MAP>& edgeangle, typename PFP::MATRIX33&);
void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& pos, typename PFP::MATRIX33&);
REAL computeArea(const VertexAttribute<VEC3, MAP>& pos, const EdgeAttribute<REAL, MAP>& edgearea);
void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& pos, const EdgeAttribute<REAL, MAP>& edgeangle, const EdgeAttribute<REAL, MAP>& edgearea, typename PFP::MATRIX33&);
};
/*********************************************************
......@@ -234,9 +227,8 @@ public:
void collectAll(Dart d);
void collectBorder(Dart d);
REAL computeArea(const VertexAttribute<VEC3, MAP>& pos);
void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& pos, const EdgeAttribute<REAL, MAP>& edgeangle, typename PFP::MATRIX33&);
void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& pos, typename PFP::MATRIX33&);
REAL computeArea(const VertexAttribute<VEC3, MAP>& pos, const EdgeAttribute<REAL, MAP>& edgearea);
void computeNormalCyclesTensor (const VertexAttribute<VEC3, MAP>& pos, const EdgeAttribute<REAL, MAP>& edgeangle, const EdgeAttribute<REAL, MAP>& edgearea, typename PFP::MATRIX33&);
};
/*********************************************************
......
This diff is collapsed.
......@@ -272,9 +272,19 @@ void Surface_DifferentialProperties_Plugin::computeCurvature(
if(!edgeAngle.isValid())
edgeAngle = mh->addAttribute<PFP2::REAL, EDGE>("edgeAngle");
EdgeAttribute<PFP2::REAL, PFP2::MAP> edgeArea = mh->getAttribute<PFP2::REAL, EDGE>("edgeArea");
if(!edgeArea.isValid())
edgeArea = mh->addAttribute<PFP2::REAL, EDGE>("edgeArea");
PFP2::MAP* map = mh->getMap();
Algo::Surface::Geometry::computeAnglesBetweenNormalsOnEdges<PFP2>(*map, position, edgeAngle);
Algo::Surface::Geometry::computeCurvatureVertices_NormalCycles_Projected<PFP2>(*map, 0.01f * mh->getBBdiagSize(), position, normal, edgeAngle, kmax, kmin, Kmax, Kmin, Knormal);
Algo::Surface::Geometry::computeAreaEdges<PFP2>(*map, position, edgeArea);
PFP2::REAL meanEdgeLength = Algo::Surface::Geometry::meanEdgeLength<PFP2>(*map, position);
float radius = 2.0f * meanEdgeLength;
Algo::Surface::Geometry::computeCurvatureVertices_NormalCycles_Projected<PFP2>(*map, radius, position, normal, edgeAngle, edgeArea, kmax, kmin, Kmax, Kmin, Knormal);
computeCurvatureLastParameters[mapName] =
ComputeCurvatureParameters(
......
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