diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index deeb7fc2c6b4e62b6c6428c0f49ee80117f43ecf..ef6b3947183a92b92616f0fc2453e6f4d239c160 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -998,24 +998,118 @@ int GFace::genusGeom()
return nSeams - single_seams.size();
}
-void GFace::fillPointCloud(double maxDist, std::vector<SPoint3> *points,
+class IntSPoint2LessThan{
+ public:
+ bool operator()(const std::pair<int, SPoint2> &v1,
+ const std::pair<int, SPoint2> &v2) const
+ {
+ return v1.second < v2.second;
+ }
+};
+
+static int addSTLVertex(std::set<std::pair<int, SPoint2>, IntSPoint2LessThan> &all,
+ std::vector<SPoint2> &vertices, SPoint2 p)
+{
+ std::pair<int, SPoint2> tmp(0, p);
+ std::set<std::pair<int, SPoint2>, IntSPoint2LessThan>::iterator it = all.find(tmp);
+ if(it == all.end()){
+ vertices.push_back(p);
+ int index = vertices.size() - 1;
+ all.insert(std::pair<int, SPoint2>(index, p));
+ return index;
+ }
+ return it->first;
+}
+
+static void recurSTLTriangle(GFace *gf, double maxDist,
+ std::set<std::pair<int, SPoint2>, IntSPoint2LessThan> &all,
+ std::vector<SPoint2> &vertices,
+ std::vector<int> &triangles, int v0, int v1, int v2)
+{
+ SPoint2 &p0(vertices[v0]); GPoint gp0 = gf->point(p0);
+ SPoint2 &p1(vertices[v1]); GPoint gp1 = gf->point(p1);
+ SPoint2 &p2(vertices[v2]); GPoint gp2 = gf->point(p2);
+ int v01 = addSTLVertex(all, vertices, (p0 + p1) * 0.5);
+ int v12 = addSTLVertex(all, vertices, (p1 + p2) * 0.5);
+ int v20 = addSTLVertex(all, vertices, (p2 + p0) * 0.5);
+ SPoint2 &p01(vertices[v01]); GPoint gp01 = gf->point(p01);
+ SPoint2 &p12(vertices[v12]); GPoint gp12 = gf->point(p12);
+ SPoint2 &p20(vertices[v20]); GPoint gp20 = gf->point(p20);
+
+ if(gp0.distance(gp01) > maxDist ||
+ gp01.distance(gp20) > maxDist ||
+ gp20.distance(gp0) > maxDist){
+ recurSTLTriangle(gf, maxDist, all, vertices, triangles, v0, v01, v20);
+ }
+ else{
+ triangles.push_back(v0); triangles.push_back(v01); triangles.push_back(v20);
+ }
+
+ if(gp1.distance(gp12) > maxDist ||
+ gp12.distance(gp01) > maxDist ||
+ gp01.distance(gp1) > maxDist){
+ recurSTLTriangle(gf, maxDist, all, vertices, triangles, v1, v12, v01);
+ }
+ else{
+ triangles.push_back(v1); triangles.push_back(v12); triangles.push_back(v01);
+ }
+
+ if(gp2.distance(gp20) > maxDist ||
+ gp20.distance(gp12) > maxDist ||
+ gp12.distance(gp2) > maxDist){
+ recurSTLTriangle(gf, maxDist, all, vertices, triangles, v2, v20, v12);
+ }
+ else{
+ triangles.push_back(v2); triangles.push_back(v20); triangles.push_back(v12);
+ }
+
+ if(gp01.distance(gp12) > maxDist ||
+ gp12.distance(gp20) > maxDist ||
+ gp20.distance(gp01) > maxDist){
+ recurSTLTriangle(gf, maxDist, all, vertices, triangles, v01, v12, v20);
+ }
+ else{
+ triangles.push_back(v01); triangles.push_back(v12); triangles.push_back(v20);
+ }
+}
+
+bool GFace::fillPointCloud(double maxDist, std::vector<SPoint3> *points,
std::vector<SVector3> *normals)
{
if(!buildSTLTriangulation()){
Msg::Error("No STL triangulation available to fill point cloud");
- return;
+ return false;
}
+
+ if(!points) return false;
+
+ // FIXME: this is too complicated -- we should just split the
+ // triangles in one step (i.e., assume thay are mostly flat) and not
+ // worry about duplicate vertices on the triangle edges)
+
+ std::set<std::pair<int, SPoint2>, IntSPoint2LessThan> all_vertices;
+ for(unsigned int i = 0; i < stl_vertices.size(); i++)
+ all_vertices.insert(std::pair<int, SPoint2>(i, stl_vertices[i]));
+
+ std::vector<SPoint2> new_vertices;
+ new_vertices.insert(new_vertices.begin(), stl_vertices.begin(),
+ stl_vertices.end());
+
+ std::vector<int> new_triangles;
+ new_triangles.insert(new_triangles.begin(), stl_triangles.begin(),
+ stl_triangles.end());
- bool computePoints = points ? true : false;
- bool computeNormals = normals ? true : false;
- if(!computePoints && !computeNormals) return;
- /*
- int N = va_geom_triangles->getNumVertices();
- for (int i = 0; i < N; i += 3) {
- SPoint3 p((va_geom_triangles->getVertexArray(3 * i))[0],
- (va_geom_triangles->getVertexArray(3 * i))[1],
- (va_geom_triangles->getVertexArray(3 * i))[2]);
+ for(unsigned int i = 0; i < stl_triangles.size(); i += 3)
+ recurSTLTriangle(this, maxDist, all_vertices, new_vertices, new_triangles,
+ stl_triangles[i], stl_triangles[i + 1], stl_triangles[i + 2]);
+
+ for(unsigned int i = 0; i < new_vertices.size(); i++){
+ GPoint p(point(new_vertices[i]));
+ points->push_back(SPoint3(p.x(), p.y(), p.z()));
+ }
+
+ if(normals){
+ for(unsigned int i = 0; i < new_vertices.size(); i++)
+ normals->push_back(normal(new_vertices[i]));
}
- */
-
}
diff --git a/Geo/GFace.h b/Geo/GFace.h
index f6b033a93b7426e45c7dbc4e2ba14ded3ce6f16d..5bef52a9bf6bf2c20825e6c9a0eaed047122e915 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -232,8 +232,8 @@ class GFace : public GEntity
// add points (and optionalluy normals) in vectors so that two
// points are at most maxDist apart
- void fillPointCloud(double maxDist, std::vector<SPoint3> *points,
- std::vector<SVector3> *normals);
+ bool fillPointCloud(double maxDist, std::vector<SPoint3> *points,
+ std::vector<SVector3> *normals=0);
struct {
// do we recombine the triangles of the mesh?
@@ -267,7 +267,12 @@ class GFace : public GEntity
// of start/end points
std::vector<SPoint3> cross;
- // a vertex array containing an STL representation of the surface
+ // the STL mesh
+ std::vector<SPoint2> stl_vertices;
+ std::vector<int> stl_triangles;
+
+ // a vertex array containing a geometrical representation of the
+ // surface
VertexArray *va_geom_triangles;
// a array for accessing the transfinite vertices using a pair of
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index 96cd08ffd6523c2508ced27e932f258affcf08ca..34dda7daf68eab635e69cd06d9366e3d491e93d9 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -312,62 +312,65 @@ bool OCCFace::buildSTLTriangulation(bool force)
return true;
}
- TopLoc_Location loc;
- int p1, p2, p3;
+ stl_vertices.clear();
+ stl_triangles.clear();
+
Bnd_Box aBox;
- Standard_Boolean bWithShare = Standard_False;
- Standard_Real aDiscret, aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
- Standard_Real dX, dY, dZ, dMax, aCoeff;
BRepBndLib::Add(s, aBox);
+ Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
aBox.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
-
- dX = aXmax-aXmin;
- dY = aYmax-aYmin;
- dZ = aZmax-aZmin;
- dMax = dX;
- if(dY > dMax) {
- dMax = dY;
- }
- if(dZ > dMax) {
- dMax = dZ;
- }
-
- aCoeff = 0.01;
- aDiscret = aCoeff * dMax;
- BRepMesh_FastDiscret aMesher(aDiscret, 0.5, aBox, bWithShare, Standard_True,
+ Standard_Real dX = aXmax - aXmin;
+ Standard_Real dY = aYmax - aYmin;
+ Standard_Real dZ = aZmax - aZmin;
+ Standard_Real dMax = dX;
+ if(dY > dMax) dMax = dY;
+ if(dZ > dMax) dMax = dZ;
+ Standard_Real aCoeff = 0.01;
+ Standard_Real aDiscret = aCoeff * dMax;
+ BRepMesh_FastDiscret aMesher(aDiscret, 0.5, aBox, Standard_False, Standard_True,
Standard_False, Standard_True);
aMesher.Add(s);
+ TopLoc_Location loc;
Handle(Poly_Triangulation) triangulation = BRep_Tool::Triangulation(s, loc);
- if (triangulation.IsNull()){
+ if(triangulation.IsNull()){
Msg::Warning("OCC STL triangulation failed");
return false;
}
-
- int ntriangles = triangulation->NbTriangles();
- if(!triangulation->HasUVNodes()) return false;
+ if(!triangulation->HasUVNodes()){
+ Msg::Warning("OCC STL triangulation has no u,v coordinates");
+ return false;
+ }
- va_geom_triangles = new VertexArray(3, ntriangles);
-
+ for(int i = 1; i <= triangulation->NbNodes(); i++){
+ gp_Pnt2d p = (triangulation->UVNodes())(i);
+ stl_vertices.push_back(SPoint2(p.X(), p.Y()));
+ }
+
+ for(int i = 1; i <= triangulation->NbTriangles(); i++){
+ Poly_Triangle triangle = (triangulation->Triangles())(i);
+ int p1, p2, p3;
+ triangle.Get(p1, p2, p3);
+ stl_triangles.push_back(p1 - 1);
+ stl_triangles.push_back(p2 - 1);
+ stl_triangles.push_back(p3 - 1);
+ }
+
+ va_geom_triangles = new VertexArray(3, stl_triangles.size() / 3);
unsigned int c = CTX::instance()->color.geom.surface;
unsigned int col[4] = {c, c, c, c};
- for (int j = 1; j <= ntriangles; j++){
- Poly_Triangle triangle = (triangulation->Triangles())(j);
- triangle.Get(p1, p2, p3);
- gp_Pnt2d x1 = (triangulation->UVNodes())(p1);
- gp_Pnt2d x2 = (triangulation->UVNodes())(p2);
- gp_Pnt2d x3 = (triangulation->UVNodes())(p3);
- GPoint gp1 = point(x1.X(), x1.Y());
- GPoint gp2 = point(x2.X(), x2.Y());
- GPoint gp3 = point(x3.X(), x3.Y());
- SVector3 n[3];
- n[0] = normal(SPoint2(x1.X(), x1.Y()));
- n[1] = normal(SPoint2(x2.X(), x2.Y()));
- n[2] = normal(SPoint2(x3.X(), x3.Y()));
+ for (unsigned int i = 0; i < stl_triangles.size(); i += 3){
+ SPoint2 &p1(stl_vertices[stl_triangles[i]]);
+ SPoint2 &p2(stl_vertices[stl_triangles[i + 1]]);
+ SPoint2 &p3(stl_vertices[stl_triangles[i + 2]]);
+ GPoint gp1 = GFace::point(p1);
+ GPoint gp2 = GFace::point(p2);
+ GPoint gp3 = GFace::point(p3);
double x[3] = {gp1.x(), gp2.x(), gp3.x()};
double y[3] = {gp1.y(), gp2.y(), gp3.y()};
double z[3] = {gp1.z(), gp2.z(), gp3.z()};
+ SVector3 n[3] = {normal(p1), normal(p2), normal(p3)};
va_geom_triangles->add(x, y, z, n, col);
- }
+ }
va_geom_triangles->finalize();
return true;
}
diff --git a/Geo/SPoint2.h b/Geo/SPoint2.h
index fae78d7c7a43a892c921a337a5c7bb1449a17506..9ca2c07e70d4cfa1ae3ec774922ef8c9cf600e96 100644
--- a/Geo/SPoint2.h
+++ b/Geo/SPoint2.h
@@ -30,6 +30,13 @@ class SPoint2 {
void operator*=(double mult);
SPoint2 operator*(double mult);
operator double *() { return P; }
+ bool operator < (const SPoint2 &other) const
+ {
+ if(other.P[0] < P[0]) return true;
+ if(other.P[0] > P[0]) return false;
+ if(other.P[1] < P[1]) return true;
+ return false;
+ }
};
inline SPoint2 operator + (const SPoint2 &a, const SPoint2 &b)