diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 03899fca412344a61247b25b431c05218d538a6e..ef4cf58442d69f439fa110cbb3fa234f21274635 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -185,10 +185,11 @@ void Centerline::importFile(std::string fileName){
mod = new GModel();
mod->readVTK(fileName.c_str());
mod->removeDuplicateMeshVertices(1.e-8);
-
std::vector<GEntity*> entities ;
- mod->getEntities(entities) ;
-
+ mod->getEntities(entities) ;
+
+ current->setAsCurrent();
+
int maxN = 0.0;
for(unsigned int i = 0; i < entities.size(); i++){
if( entities[i]->dim() == 1){
@@ -209,14 +210,13 @@ void Centerline::importFile(std::string fileName){
}
}
- //splitEdges
- splitEdges(maxN);
-
- current->setAsCurrent();
+ createBranches(maxN);
+
}
-void Centerline::splitEdges(int maxN){
+
+void Centerline::createBranches(int maxN){
//sort colored lines and create edges
std::vector<std::vector<MLine*> > color_edges;
@@ -283,15 +283,13 @@ void Centerline::splitEdges(int maxN){
else itl++;
}
if (vB == vE) { Msg::Error("Begin and end points branch are the same \n");break;}
- orderMLines(myLines); //, junctions, vBB, vEE);
+ orderMLines(myLines);
std::vector<Branch> children;
- Branch newBranch ={ tag++, myLines, computeLength(myLines), vB, vE, children};
+ Branch newBranch ={ tag++, myLines, computeLength(myLines), vB, vE, children, 1.e6, 0.0};
//printf("VB = %d %d \n", vB->getNum(), vE->getNum());
edges.push_back(newBranch) ;
}
}
-
- //create new edges
printf("edges =%d new =%d \n", color_edges.size(), edges.size());
//create children
@@ -305,6 +303,10 @@ void Centerline::splitEdges(int maxN){
edges[i].children = myChildren;
}
+ //compute radius
+ distanceToLines();
+ computeRadii();
+
//print info
// for(unsigned int i = 0; i < edges.size(); ++i) {
// printf("EDGE =%d tag=%d length = %g childs = %d \n", i, edges[i].tag, edges[i].length, edges[i].children.size());
@@ -321,6 +323,69 @@ void Centerline::splitEdges(int maxN){
}
+void Centerline::distanceToLines(){
+
+ std::vector<SPoint3> pts;
+ std::set<SPoint3> pts_set;
+ std::vector<double> distances;
+ std::vector<double> distancesE;
+ std::vector<SPoint3> closePts;
+
+ GModel *current = GModel::current();
+ std::vector<GEntity*> _entities ;
+ current->getEntities(_entities) ;
+ for(unsigned int i = 0; i < _entities.size(); i++){
+ if( _entities[i]->dim() != 2) continue;
+ for(unsigned int j = 0; j < _entities[i]->getNumMeshElements(); j++){
+ MElement *e = _entities[i]->getMeshElement(j);
+ for (int k = 0; k < e->getNumVertices(); k++){
+ MVertex *v = e->getVertex(k);
+ pts_set.insert(SPoint3(v->x(), v->y(),v->z()));
+ }
+ }
+ }
+ std::set<SPoint3>::iterator its = pts_set.begin();
+ for (; its != pts_set.end(); its++){
+ pts.push_back(*its);
+ }
+ if (pts.size() == 0) {Msg::Error("Centerline needs a GModel with a surface \n"); exit(1);}
+
+ for(unsigned int i = 0; i < lines.size(); i++){
+ std::vector<double> iDistances;
+ std::vector<SPoint3> iClosePts;
+ std::vector<double> iDistancesE;
+ MLine *l = lines[i];
+ MVertex *v1 = l->getVertex(0);
+ MVertex *v2 = l->getVertex(1);
+ SPoint3 p1(v1->x(), v1->y(), v1->z());
+ SPoint3 p2(v2->x(), v2->y(), v2->z());
+ signedDistancesPointsLine(iDistances, iClosePts, pts, p1,p2);
+
+ double minRad = std::abs(iDistances[0]);
+ for (unsigned int kk = 1; kk< pts.size(); kk++) {
+ if (std::abs(iDistances[kk]) < minRad)
+ minRad = std::abs(iDistances[kk]);
+ }
+ radiusl.insert(std::make_pair(lines[i], minRad));
+ }
+
+}
+void Centerline::computeRadii(){
+
+ for(unsigned int i = 0; i < edges.size(); ++i) {
+ std::vector<MLine*> lines = edges[i].lines;
+ for(unsigned int j = 0; j < lines.size(); ++j) {
+ MLine *l = lines[j];
+ std::map<MLine*,double>::iterator itr = radiusl.find(l);
+ if (itr != radiusl.end()){
+ edges[i].minRad = std::min(itr->second, edges[i].minRad);
+ edges[i].maxRad = std::min(itr->second, edges[i].maxRad);
+ }
+ else printf("ARGG line not found \n");
+ }
+ }
+
+}
void Centerline::buildKdTree(){
FILE * f = fopen("myPOINTS.pos","w");
@@ -388,6 +453,7 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge){
}
void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge){
+ Msg::Error("This anisotropic operator of CenterlineField is not implemnted yet ");
return;
}
@@ -395,28 +461,6 @@ void Centerline::printSplit() const{
FILE * f = fopen("mySPLIT.pos","w");
fprintf(f, "View \"\"{\n");
-
- // for(unsigned int i = 0; i < lines.size(); ++i){
- // MLine *l = lines[i];
- // std::map<MLine*,int>::const_iterator itc = colorl.find(l);
- // std::map<MVertex*,int>::const_iterator it0 =
- // colorp.find(l->getVertex(0));
- // std::map<MVertex*,int>::const_iterator it1 =
- // colorp.find(l->getVertex(1));
- // fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
- // l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(),
- // l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(),
- // (double)itc->second,(double)itc->second);
- // }
- // std::map<MVertex*,int>::const_iterator itp = colorp.begin();
- // while (itp != colorp.end()){
- // MVertex *v = itp->first;
- // fprintf(f, "SP(%g,%g,%g){%g};\n",
- // v->x(), v->y(), v->z(),
- // 0.0);
- // itp++;
- // }
-
for(unsigned int i = 0; i < edges.size(); ++i){
std::vector<MLine*> lines = edges[i].lines;
for(unsigned int k = 0; k < lines.size(); ++k){
@@ -427,7 +471,6 @@ void Centerline::printSplit() const{
(double)i, (double)i);
}
}
-
fprintf(f,"};\n");
fclose(f);
@@ -456,7 +499,19 @@ void Centerline::printSplit() const{
}
fprintf(f3,"};\n");
fclose(f3);
-
+
+ FILE * f4 = fopen("myRADII.pos","w");
+ fprintf(f4, "View \"\"{\n");
+ for(unsigned int i = 0; i < lines.size(); ++i){
+ MLine *l = lines[i];
+ std::map<MLine*,double>::const_iterator itc = radiusl.find(l);
+ fprintf(f, "SL(%g,%g,%g,%g,%g,%g){%g,%g};\n",
+ l->getVertex(0)->x(), l->getVertex(0)->y(), l->getVertex(0)->z(),
+ l->getVertex(1)->x(), l->getVertex(1)->y(), l->getVertex(1)->z(),
+ itc->second,itc->second);
+ }
+ fprintf(f4,"};\n");
+ fclose(f4);
}
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index 79f60fa1d96394274f0fae92391a674feafee864..b0e5480cebdaaaa4b52be977366b7f6f2fb7b31b 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -23,6 +23,7 @@ class GEntity;
#include <ANN/ANN.h>
class ANNkd_tree;
+// A branch of a 1D tree
struct Branch{
int tag;
std::vector<MLine*> lines;
@@ -30,10 +31,20 @@ struct Branch{
MVertex *vB;
MVertex *vE;
std::vector<Branch> children;
+ double minRad;
+ double maxRad;
};
+// This class takes as input A 1D mesh which is the centerline
+// of a tubular 2D surface mesh
+// It computes a mesh size field function of the distance to the centerlines
+// and a cross field that is the direction of the centerline
+// It splits the tubular structure in many mesh partitions
+// along planes perpendicuar to the centerlines
+
class Centerline : public Field{
+ protected:
ANNkd_tree *kdtree;
ANNpointArray nodes;
ANNidxArray index;
@@ -41,14 +52,23 @@ class Centerline : public Field{
GModel *mod;
std::string fileName;
+ //all (unique) lines of centerlines
+ std::vector<MLine*> lines;
+ //the stuctured tree of the centerlines
+ std::vector<Branch> edges;
+ //the radius of the surface mesh at a given line
+ std::map<MLine*,double> radiusl;
+ //the junctions of the tree
+ std::set<MVertex*> junctions;
+ //some colors (int) for all points and lines
+ std::map<MVertex*,int> colorp;
+ std::map<MLine*,int> colorl;
+
public:
Centerline(std::string fileName);
Centerline();
~Centerline();
- void importFile(std::string fileName);
- void splitEdges(int maxN);
-
virtual bool isotropic () const {return false;}
virtual const char *getName()
{
@@ -62,18 +82,35 @@ class Centerline : public Field{
" using the following script:\n\n"
"vmtk vmtkcenterlines -seedselector openprofiles -ifile mysurface.stl -ofile centerlines.vtp --pipe vmtksurfacewriter -ifile centerlines.vtp -ofile centerlines.vtk\n";
}
+
+ //isotropic operator for mesh size field function of distance to centerline
double operator() (double x, double y, double z, GEntity *ge=0);
+ //anisotropic operator
void operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge=0);
- void printSplit() const;
+ //import the 1D mesh of the centerlines (in vtk format)
+ //and fill the vector of lines
+ void importFile(std::string fileName);
+
+ //refine the 1D mesh to have many points on the 1D centerlines
+ //for annKDTree distance computations
void buildKdTree();
- std::vector<MLine*> lines;
- std::vector<Branch> edges;
- std::map<MVertex*,int> colorp;
- std::map<MLine*,int> colorl;
+ //Creates the branch structure (topology, connectivity) from the
+ //vector of lines
+ void createBranches(int maxN);
+
+ //Computes for the Branches the min and maxRadius of the tubular structure
+ //this function needs teh current GModel
+ void computeRadii();
+
+ //Computes for each MLine the minRadius
+ void distanceToLines();
+
+ //Print for debugging
+ void printSplit() const;
+
- std::set<MVertex*> junctions;
};
#endif