{
_face = Face;
_attractorShape = Attractor;
- //_mapGrid = 200;
init();
}
bool BLSURFPlugin_Attractor::init(){
- MESSAGE("construction d'un attracteur");
Standard_Real u0,v0;
- int i,j ;
- //_known.clear();
+ int i,j,i0,j0 ;
+ _known.clear();
_trial.clear();
Handle(Geom_Surface) aSurf = BRep_Tool::Surface(_face);
- Trial_Pnt TPnt(3,0); //TPnt(0,0,0);
-
- // Discretization of the parameters
- MESSAGE("Calcul des bornes de la surface : ");
- //aSurf->Bounds(_u1, _u2, _v1, _v2); // unusable in the generic case because the surface may be infinite (ok for prototype on a Sphere)
+ Trial_Pnt TPnt(3,0);
- // Calculation of the bounds for the generic case
+ // Calculation of the bounds of the face
ShapeAnalysis::GetFaceUVBounds(_face,_u1,_u2,_v1,_v2);
MESSAGE("u1 = "<<_u1<<" ,u2 = "<<_u2);
MESSAGE("v1 = "<<_v1<<" ,v2 = "<<_v2);
_gridU = floor (_u2 - _u1) / _step;
_gridV = floor (_v2 - _v1) / _step;
-// double Ustep = (_u2 - _u1) / _mapGrid;
-// double Vstep = (_v2 - _v1) / _mapGrid;
for (i=0; i<=_gridU; i++){
_vectU.push_back(_u1+i*(_u2-_u1)/_gridU) ;
}
for (j=0; j<=_gridV; j++){
_vectV.push_back(_v1+j*(_v2-_v1)/_gridV) ;
}
-
- // Determination of the starting point
- gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(_attractorShape));
- GeomAPI_ProjectPointOnSurf projector( P, aSurf );
- projector.LowerDistanceParameters(u0,v0);
- MESSAGE("u0 = "<<u0<<" ,v0 = "<<v0);
- //gp_Pnt2d P = BRep_Tool::Parameters(TopoDS::Vertex(_attractorShape),_face);
- //u0 = P.X();
- //v0 = P.Y();
- int i0 = floor ( (u0 - _u1) * _gridU / (_u2 - _u1) + 0.5 );
- int j0 = floor ( (v0 - _v1) * _gridV / (_v2 - _v1) + 0.5 );
- MESSAGE("i0 = "<<i0<<" ,j0 = "<<j0);
- TPnt[0]=0.;//.dist=0.; // Set distance to 0.
- TPnt[1]=i0;//.Ti=i0;
- TPnt[2]=j0;//.Tj=j0;
- _trial.insert(TPnt); // Move starting point to _trial
- // Initialization of _DMap
- //std::vector<double> temp(_mapGrid+1,std::numeric_limits<double>::infinity()); // Set distance of all "far" points to Infinity
- std::vector<double> temp(_gridV+1,1.e+6); // Provisoirement infini = 10000
+ // Initialization of _DMap and _known
+ std::vector<double> temp(_gridV+1,std::numeric_limits<double>::infinity()); // Set distance of all "far" points to Infinity
for (i=0; i<=_gridU; i++){
_DMap.push_back(temp);
}
for (i=0; i<=_gridU; i++){
_known.push_back(temp2);
}
- MESSAGE("DMap(0,0) : "<< _DMap[0][0]);
- _DMap[i0][j0] = 0.; // Set distance of starting point to 0.
-
- // _buildMap(); // Computes the distance for all points of the discrete surface
-}
+ // Determination of the starting points
+ if (_attractorShape.ShapeType() == TopAbs_VERTEX){
+ gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(_attractorShape));
+ GeomAPI_ProjectPointOnSurf projector( P, aSurf );
+ projector.LowerDistanceParameters(u0,v0);
+ MESSAGE("u0 = "<<u0<<" ,v0 = "<<v0);
+ i0 = floor ( (u0 - _u1) * _gridU / (_u2 - _u1) + 0.5 );
+ j0 = floor ( (v0 - _v1) * _gridV / (_v2 - _v1) + 0.5 );
+ TPnt[0]=0.; // Set the distance of the starting point to 0.
+ TPnt[1]=i0;
+ TPnt[2]=j0;
+ _DMap[i0][j0] = 0.;
+ _trial.insert(TPnt); // Move starting point to _trial
+ }
+ else if (_attractorShape.ShapeType() == TopAbs_EDGE){
+ gp_Pnt2d P2;
+ double first;
+ double last;
+ Handle(Geom2d_Curve) Curve2d = BRep_Tool::CurveOnSurface(TopoDS::Edge(_attractorShape), _face, first, last);
+ int N = 2 * (last - first) / _step;
+ for (i=0; i<=N; i++){
+ P2 = Curve2d->Value(first + i * (last-first) / N);
+ i0 = floor( (P2.X() - _u1) * _gridU / (_u2 - _u1) + 0.5 );
+ j0 = floor( (P2.Y() - _v1) * _gridV / (_v2 - _v1) + 0.5 );
+ MESSAGE("i0 = "<<i0<<" , j0 = "<<j0)
+ TPnt[0] = 0.;
+ TPnt[1] = i0;
+ TPnt[2] = j0;
+ _DMap[i0][j0] = 0.;
+ _trial.insert(TPnt);
+ }
+ }
+}
void BLSURFPlugin_Attractor::SetParameters(double Start_Size, double End_Size, double Action_Radius, double Constant_Radius){
_startSize = Start_Size;
_endSize = End_Size;
_actionRadius = Action_Radius;
_constantRadius = Constant_Radius;
- MESSAGE("_startSize = "<<_startSize<<"; _endSize = "<<_endSize<<" ; _actionRadius = "<<_actionRadius<<" ; _constantRadius = "<<_constantRadius)
}
double BLSURFPlugin_Attractor::_distance(double u, double v){
// BLSURF seems to perform a linear interpolation so it's sufficient to give it a non-continuous distance map
int i = floor ( (u - _u1) * _gridU / (_u2 - _u1) + 0.5 );
int j = floor ( (v - _v1) * _gridV / (_v2 - _v1) + 0.5 );
-
-// int i = floor ( (u - _u1) * _mapGrid / (_u2 - _u1) );
-// int j = floor ( (v - _v1) * _mapGrid / (_v2 - _v1) );
return _DMap[i][j];
-
}
double BLSURFPlugin_Attractor::GetSize(double u, double v){
- double myDist = 0.5 * (_distance(u,v) - _constantRadius + abs(_distance(u,v) - _constantRadius));
+ double myDist = 0.5 * (_distance(u,v) - _constantRadius + fabs(_distance(u,v) - _constantRadius));
+ if (myDist<0){
+ MESSAGE("Warning myDist<0 : myDist= "<<myDist)
+ }
switch(_type)
{
case TYPE_EXP:
MESSAGE("building the map");
int i, j, k, n;
int count = 0;
- //TPointSet::iterator pIt;
int ip, jp, kp, np;
int i0, j0;
gp_Pnt P;
double Dist = 0.;
bool Dist_changed;
IJ_Pnt Current_Pnt(2,0);
- Trial_Pnt TPnt(3,0);//TPnt(0,0,0);
+ Trial_Pnt TPnt(3,0);
TTrialSet::iterator min;
TTrialSet::iterator found;
- //TPointSet::iterator kn_found;
Handle(Geom_Surface) aSurf = BRep_Tool::Surface(_face);
- while (_trial.size() > 0 && count < 1000000){
- count++;
- MESSAGE("boucle while : "<<count);
- min = _trial.begin(); // Get trial point with min distance from start
- i0 = (*min)[1];//min->Ti;//static_cast<int>((*min)[1]);
- j0 = (*min)[2];//min->Tj;//static_cast<int>((*min)[2]);
- MESSAGE("i0 = "<<i0<<" ,j0 = "<<j0<<" *min[0] = "<<(*min)[0]);
- //MESSAGE("d(i0,j0) = "<<min->dist);//(*min)[0]);
-// Current_Pnt[0] = i0; // Define it as the current point
-// Current_Pnt[1] = j0;
- //if (i0 == 23 && j0 == 49){
- //MESSAGE(" point 23 49 insere dans known");
- //break;
- //}
-// pIt =_known.find( Current_Pnt );
- /*if (*pIt == test){
- MESSAGE(" point 23 49 insere dans known");
- break;
- }*/
-// if (pIt != _known.end()){
-// MESSAGE("point "<<(*pIt)[0]<<" , "<<(*pIt)[1]<<" deja dans _known -> probleme");
-// break;
-// }
-// _known.insert(Current_Pnt); // Move it to _known and remove it from _trial
- _known[i0][j0] = true; // Move point to "Known"
- //MESSAGE(Current_Pnt[0]<<" , "<<Current_Pnt[1]<<" insere dans known")
- //MESSAGE("known contient : "<<_known.size()<<"elements");
- //MESSAGE(" d("<<i0<<" , "<<j0<<") ="<<_DMap[i0][j0]);
- //_DMap[i0][j0]=(*min)[0];//min->dist;//(*min)[0];
- //MESSAGE("avt erase trial contient : "<<_trial.size()<<"elts");
-// MESSAGE((*min)[0]<<", "<<(*min)[1]<<", "<<(*min)[2]<<" supprime de Trial");
- _trial.erase(min);
- //MESSAGE("trial contient mtnt : "<<_trial.size()<<"elts");
- // ------------------------ // Loop on neighbours of the trial min
+ // While there are points in "Trial" (representing a kind of advancing front), loop on them ------------------
+ while (_trial.size() > 0 ){
+ min = _trial.begin(); // Get trial point with min distance from start
+ i0 = (*min)[1];
+ j0 = (*min)[2];
+ _known[i0][j0] = true; // Move it to "Known"
+ _trial.erase(min); // Remove it from "Trial"
+ // Loop on neighbours of the trial min ----------------------------------------------------------------------
for (i=i0 - 1 ; i <= i0 + 1 ; i++){
- //MESSAGE("boucle for sur i");
- if (!aSurf->IsUPeriodic()){ // Periodic conditions in U
- //MESSAGE("surface is not U periodique");
+ if (!aSurf->IsUPeriodic()){ // Periodic conditions in U
if (i > _gridU ){
break; }
else if (i < 0){
i++; }
}
- ip = (i + _gridU + 1) % (_gridU+1); // We get a periodic index ip=modulo(i,N+2) so that i=-1->ip=N; i=0 -> ip=0 ; ... ; i=N+1 -> ip=0;
+ ip = (i + _gridU + 1) % (_gridU+1); // We get a periodic index ip=modulo(i,N+2) so that i=-1->ip=N; i=0 -> ip=0 ; ... ; i=N+1 -> ip=0;
for (j=j0 - 1 ; j <= j0 + 1 ; j++){
- //MESSAGE("boucle for sur j");
- if (!aSurf->IsVPeriodic()){ // Periodic conditions in V .
- //MESSAGE("surface is not V periodique");
+ if (!aSurf->IsVPeriodic()){ // Periodic conditions in V .
if (j > _gridV ){
break; }
else if (j < 0){
j++;
- if (j<0){
- MESSAGE("j = "<<j);
- //break;
- }
}
}
jp = (j + _gridV + 1) % (_gridV+1);
- //jp=j;
-// Current_Pnt[0] = ip; // Define the neighbour as current point
-// Current_Pnt[1] = jp;
- //MESSAGE("i = "<<i<<" ; j = "<<j);
- //MESSAGE("ip = "<<ip<<" ; jp = "<<jp);
-// pIt = _known.find( Current_Pnt );
- //MESSAGE("found a point in Known with i = "<<(*pIt)[0]<<"j = "<<(*pIt)[1]);
- if (!_known[ip][jp]){ //.find( Current_Pnt ) == _known.end()){ // If the distance is not known
- //MESSAGE("Le voisin "<<ip<<","<<jp<<" de min(_trial) n'est pas dans _known");
- aSurf->D1(_vectU[ip],_vectV[jp],P,D1U,D1V); // Calculate the metric at (i,j)
- Guu = D1U.X()*D1U.X() + D1U.Y()*D1U.Y() + D1U.Z()*D1U.Z(); // Guu = ||dS/du||**2 G(i,j)= | ||dS/du||**2 * |
- Gvv = D1V.X()*D1V.X() + D1V.Y()*D1V.Y() + D1V.Z()*D1V.Z(); // Gvv = ||dS/dv||**2 | <dS/du,dS/dv> ||dS/dv||**2 |
+
+ if (!_known[ip][jp]){ // If the distance is not known yet
+ aSurf->D1(_vectU[ip],_vectV[jp],P,D1U,D1V); // Calculate the metric at (i,j)
+ // G(i,j) = | ||dS/du||**2 * |
+ // | <dS/du,dS/dv> ||dS/dv||**2 |
+ Guu = D1U.X()*D1U.X() + D1U.Y()*D1U.Y() + D1U.Z()*D1U.Z(); // Guu = ||dS/du||**2
+ Gvv = D1V.X()*D1V.X() + D1V.Y()*D1V.Y() + D1V.Z()*D1V.Z(); // Gvv = ||dS/dv||**2
Guv = D1U.X()*D1V.X() + D1U.Y()*D1V.Y() + D1U.Z()*D1V.Z(); // Guv = Gvu = < dS/du,dS/dv >
- //Guv = 0; // test
- if (Guv > 1e-15){
- MESSAGE("WARNING : Guv != 0 ; Guv = "<<Guv);
- }
- //MESSAGE("u ("<<ip<<") = "<< _vectU[ip]<<"v ("<<jp<<") = "<< _vectV[jp]);
- if (count == 2500) {
- MESSAGE("Guu = " << Guu<<" , Gvv = "<<Gvv);
- }
- D_Ref = _DMap[ip][jp]; // Set a ref. distance to the value in DMap (may be infinite or uncertain)
- //MESSAGE("D_Ref = "<<D_Ref);
- TPnt[0] = D_Ref;//.dist=D_Ref; // Store the point as a trial point
- TPnt[1] = ip;//.Ti=ip;//[1] = static_cast<double>(ip);
- TPnt[2] = jp;//.Tj=jp;//[2] = static_cast<double>(jp);
+ D_Ref = _DMap[ip][jp]; // Set a ref. distance of the point to its value in _DMap (may be infinite or uncertain)
+ TPnt[0] = D_Ref;
+ TPnt[1] = ip;
+ TPnt[2] = jp;
Dist_changed = false;
- // -------------------------------------------- // Calculate the min distance among the neighbours
+ // Loop on neighbours to calculate the min distance from them ---------------------------------------------------------------------------------
for (k=i - 1 ; k <= i + 1 ; k++){
- //MESSAGE("boucle for sur k");
- if (!aSurf->IsUPeriodic()){ // Periodic conditions in U
+ if (!aSurf->IsUPeriodic()){ // Periodic conditions in U
if(k > _gridU ){
- MESSAGE("Break");
break;
}
else if (k < 0){
k++; }
}
- kp = (k + _gridU + 1) % (_gridU+1); // periodic index
- //MESSAGE(" k = "<<k<<" ; kp = "<<kp);
+ kp = (k + _gridU + 1) % (_gridU+1); // periodic index
for (n=j - 1 ; n <= j + 1 ; n++){
- //MESSAGE("boucle for sur n");
- if (!aSurf->IsVPeriodic()){ // Periodic conditions in V
- if(n > _gridV){
- MESSAGE("Break");
+ if (!aSurf->IsVPeriodic()){ // Periodic conditions in V
+ if(n > _gridV){
break;
}
else if (n < 0){
n++; }
}
- np = (n + _gridV + 1) % (_gridV+1);
- //MESSAGE(" n = "<<n<<" ; np = "<<np);
-// Current_Pnt[0] = kp;
-// Current_Pnt[1] = np;
- if (_known[kp][np]){//.find( Current_Pnt ) != _known.end()){ // If distance of the neighbour is known
- // Calculate the distance from (k,n)
- //MESSAGE("La distance du point ( "<<kp<<","<<np<<" ) de min(_trial) est connue");
+ np = (n + _gridV + 1) % (_gridV+1);
+ if (_known[kp][np]){ // If the distance of the neighbour is known
+ // Calculate the distance from (k,n)
du = (k-i) * (_u2 - _u1) / _gridU;
dv = (n-j) * (_v2 - _v1) / _gridV;
- Dist = _DMap[kp][np] + sqrt( Guu * du*du + Gvv * dv*dv + 2*Guv * du*dv); // ds**2 = du'Gdu + dv'Gdv + 2*du'Gdv (G is symetrical)
- //MESSAGE("D_Ref = "<<D_Ref);
- //MESSAGE("Dist = "<<Dist);
- //bool comp = Dist < D_Ref;
- //MESSAGE("Dist < D_Ref ="<<comp);
- if (Dist < D_Ref) { // If smaller than ref. distance -> update ref. distance
- //MESSAGE("update de la distance");
+ Dist = _DMap[kp][np] + sqrt( Guu * du*du + 2*Guv * du*dv + Gvv * dv*dv ); // ds**2 = du'Gdu + 2*du'Gdv + dv'Gdv (G is always symetrical)
+ if (Dist < D_Ref) { // If smaller than ref. distance -> update ref. distance
D_Ref = Dist;
Dist_changed = true;
- if (count == 2500) {
- MESSAGE("du = "<<du<<" , dv = "<<dv<<" , D_Ref = "<<D_Ref);
- MESSAGE("_DMap["<<kp<<"]["<<np<<"] = "<<_DMap[kp][np]);
- MESSAGE("_DMap[kp][np] + sqrt( Guu * du*du + Gvv * dv*dv ) ="<<_DMap[kp][np]<<" + sqrt( "<<Guu<<" * "<<du<<"**2 + "<<Gvv<<" * "<<dv<<"**2" );
- //MESSAGE("k - i = "<<k-i<<" , n - j = "<<n-j);
- //MESSAGE("(_v2 - _v1) = "<<_v2 - _v1 );
- }
}
}
- //else
- //MESSAGE("La distance du point ( "<<kp<<","<<np<<" ) de min(_trial) est INconnue");
}
- }
- if (Dist_changed) { // If distance has been updated, update _trial
- //MESSAGE("trial contient : "<<_trial.size()<<"elements");
- //MESSAGE("Dist updated");
+ } // End of the loop on neighbours --------------------------------------------------------------------------------------------------------------
+ if (Dist_changed) { // If distance has been updated, update _trial
found=_trial.find(TPnt);
if (found != _trial.end()){
- //MESSAGE("pnt : "<<(*found)[0]<<" , "<<(*found)[1]<<" , "<<(*found)[2]<<" efface de trial (avt update)");
- //MESSAGE("trial contient : "<<_trial.size()<<"elements");
- _trial.erase(found);
+ _trial.erase(found); // Erase the point if it was already in _trial
}
- TPnt[0] = D_Ref;//.dist = D_Ref;
- TPnt[1] = ip;//.Ti = ip;//[1] = static_cast<double>(ip);
- TPnt[2] = jp;//.Tj = jp;//[2] = static_cast<double>(jp);
- _DMap[ip][jp] = D_Ref;
- //MESSAGE("TPnt[0] = "<<TPnt[0]);
- //MESSAGE("updated distance d("<<ip<<","<<jp<<") = "<<D_Ref);
- //MESSAGE("Guu = "<<Guu<<","<<"Gvv = "<<Gvv);
- //found=_trial.find(TPnt);
- //MESSAGE("Recherche de : "<<TPnt[0]<<" , "<<TPnt[1]<<" , "<<TPnt[2]);
- //if (found != _trial.end()){
- //MESSAGE("pnt : "<<(*found)[0]<<" , "<<(*found)[1]<<" , "<<(*found)[2]<<" trouve");
- //MESSAGE("trial contient : "<<_trial.size()<<"elements");
- //}
- std::pair< TTrialSet::iterator, bool> result = _trial.insert(TPnt);
-// if (result.second) {
-// MESSAGE("Pnt : "<<TPnt[0]<<" , "<<TPnt[1]<<" , "<<TPnt[2]<<" insere ds trial");
-// }
-// else {
-// MESSAGE(" !!!!!!! Pnt : "<<TPnt[0]<<" , "<<TPnt[1]<<" , "<<TPnt[2]<<" pas insere ds trial");
-// }
-
- //MESSAGE("trial contient : "<<_trial.size()<<"elements");
+ TPnt[0] = D_Ref;
+ TPnt[1] = ip;
+ TPnt[2] = jp;
+ _DMap[ip][jp] = D_Ref; // Set it distance to the minimum distance found during the loop above
+ _trial.insert(TPnt); // Insert it (or reinsert it) in _trial
}
- //MESSAGE("trial contient : "<<_trial.size()<<"elements");
} // if
} // for
} // for
_known.clear();
_trial.clear();
_isMapBuilt = true;
-
+ MESSAGE("_gridU = "<<_gridU<<" , _gridV = "<<_gridV)
} // end of BuildMap()
