#include "BLSURFPlugin_Attractor.hxx"
+#include <utilities.h>
BLSURFPlugin_Attractor::BLSURFPlugin_Attractor ()
: _face(),
_u1 (0.),
_u2 (0.),
_v1 (0.),
- _v2 (0.){}
+ _v2 (0.){MESSAGE("construction d'un attracteur vide");}
BLSURFPlugin_Attractor::BLSURFPlugin_Attractor (TopoDS_Face Face, TopoDS_Shape Attractor)
: _face(),
{
_face = Face;
_attractorShape = Attractor;
- _mapGrid = 50;
+ _mapGrid = 200;
init();
}
bool BLSURFPlugin_Attractor::init(){
-
+ MESSAGE("construction d'un attracteur");
Standard_Real u0,v0;
int i,j ;
_known.clear();
_trial.clear();
Handle(Geom_Surface) aSurf = BRep_Tool::Surface(_face);
- Trial_Pnt TPnt(0,0,0);
+ 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)
+ MESSAGE("u1 = "<<_u1<<" ,u2 = "<<_u2);
+ MESSAGE("v1 = "<<_v1<<" ,v2 = "<<_v2);
double Ustep = (_u2 - _u1) / _mapGrid;
double Vstep = (_v2 - _v1) / _mapGrid;
for (i=0; i<=_mapGrid; i++){
_vectU.push_back(_u1+i*(_u2-_u1)/_mapGrid) ;
- _vectV.push_back(_u1+i*(_u2-_u1)/_mapGrid) ;
+ _vectV.push_back(_v1+i*(_v2-_v1)/_mapGrid) ;
}
// Determination of the starting point
- GeomAPI_ProjectPointOnSurf projector( point, surface );
+ 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) * _mapGrid / (_u2 - _u1) + 0.5 );
int j0 = floor ( (v0 - _v1) * _mapGrid / (_v2 - _v1) + 0.5 );
- TPnt.dist=0.; // Set distance to 0.
- TPnt.Ti=i0;
- TPnt.Tj=j0;
+ 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(_mapGrid+1,std::numeric_limits<double>::infinity()); // Set distance of all "far" points to Infinity
+ std::vector<double> temp(_mapGrid+1,1.e+6); // Provisoirement infini = 10000
for (i=0; i<=_mapGrid; i++){
_DMap.push_back(temp);
}
+ 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
double BLSURFPlugin_Attractor::GetSize(double u, double v){
int i = floor ( (u - _u1) * _mapGrid / (_u2 - _u1) + 0.5 );
int j = floor ( (v - _v1) * _mapGrid / (_v2 - _v1) + 0.5 );
- return _DMap[i][j]; // more generally: function of _DMap
+ return 0.1 * ( 1. + ( 0.5 * (_DMap[i][j] - 20. + abs(_DMap[i][j] - 20.)) ) ); // more generally: function of _DMap
}
-bool BLSURFPlugin_Attractor::_buildMap(){
- int i, j, k, n;
- int ip, jp, kp, np;
+bool BLSURFPlugin_Attractor::_buildMap(){
+
+ 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;
gp_Vec D1U,D1V;
double Guu, Gvv; //Diagonal components of the local metric tensor
double du, dv;
double D_Ref = 0.;
double Dist = 0.;
- bool Dist_changed = false;
- IJ_Pnt Current_Pnt(0,0);
- Trial_Pnt TPnt(0,0,0);
+ bool Dist_changed;
+ IJ_Pnt Current_Pnt(2,0);
+ Trial_Pnt TPnt(3,0);//TPnt(0,0,0);
TTrialSet::iterator min;
+ TTrialSet::iterator found;
+ //TPointSet::iterator kn_found;
Handle(Geom_Surface) aSurf = BRep_Tool::Surface(_face);
- while (_trial.size() > 0){
-
+ while (_trial.size() > 0 && count < 40000){
+ count++;
+ MESSAGE("boucle while : "<<count);
min = _trial.begin(); // Get trial point with min distance from start
- int i0 = min->Ti;
- int j0 = min->Tj;
- Current_Pnt.Pi = i0; // Define it as the current point
- Current_Pnt.Pj = j0;
+ 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
+ //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
- for (i=i0 - 1 ; i < i0 + 1 ; i++){
+ for (i=i0 - 1 ; i <= i0 + 1 ; i++){
+ //MESSAGE("boucle for sur i");
if (!aSurf->IsUPeriodic()){ // Periodic conditions in U
- if (i > _mapGrid + 1){
+ //MESSAGE("surface is not U periodique");
+ if (i > _mapGrid ){
break; }
else if (i < 0){
i++; }
}
- ip = i % _mapGrid+2; // We get a periodic index ip=modulo(i,N+2) so that i=-1->ip=N+1; i=0 -> ip=0 ; ... ; i=N+2 -> ip=0;
- for (j=j0 - 1 ; j < j0 + 1 ; j++){
- if (!aSurf->IsVPeriodic()){ // Periodic conditions in V .
- if (j > _mapGrid + 1){
+ ip = (i + _mapGrid + 1) % (_mapGrid+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 (j > _mapGrid ){
break; }
else if (j < 0){
- j++; }
+ j++;
+ if (j<0){
+ MESSAGE("j = "<<j);
+ break;
+ }
+ }
}
- jp = j % _mapGrid+2;
- Current_Pnt.Pi = ip; // Define the neighbour as current point
- Current_Pnt.Pj = jp;
+ //jp = (j + _mapGrid + 1) % (_mapGrid+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.find( Current_Pnt ) == _known.end()){ // If the distance is not known
- aSurf->D1(_vectU[i],_vectV[j],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 0 | assuming that the isolines are orthogonal in 3D space
- Gvv = D1V.X()*D1V.X() + D1V.Y()*D1V.Y() + D1V.Z()*D1V.Z(); // Gvv = ||dS/dv||**2 | 0 ||dS/dv||**2 |
+ //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 0 | assuming that the isolines are orthogonal in 3D space
+ Gvv = D1V.X()*D1V.X() + D1V.Y()*D1V.Y() + D1V.Z()*D1V.Z(); // Gvv = ||dS/dv||**2 | 0 ||dS/dv||**2 |
+ //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)
- TPnt.dist = D_Ref; // Store the point as a trial point
- TPnt.Ti = ip;
- TPnt.Tj = jp;
+ //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);
+ Dist_changed = false;
// -------------------------------------------- // Calculate the min distance among the neighbours
- for (k=i - 1 ; k < i + 1 ; k++){
+ for (k=i - 1 ; k <= i + 1 ; k++){
+ //MESSAGE("boucle for sur k");
if (!aSurf->IsUPeriodic()){ // Periodic conditions in U
- if(k > _mapGrid + 1){
- break; }
+ if(k > _mapGrid ){
+ MESSAGE("Break");
+ break;
+ }
else if (k < 0){
k++; }
}
- kp = k % _mapGrid+2; // periodic index
- for (n=j - 1 ; n < j + 1 ; n++){
- if (!aSurf->IsUPeriodic()){ // Periodic conditions in V
- if(n > _mapGrid + 1){
- break; }
+ kp = (k + _mapGrid + 1) % (_mapGrid+1); // periodic index
+ //MESSAGE(" k = "<<k<<" ; kp = "<<kp);
+ for (n=j - 1 ; n <= j + 1 ; n++){
+ //MESSAGE("boucle for sur n");
+ if (!aSurf->IsVPeriodic()){ // Periodic conditions in V
+ if(n > _mapGrid){
+ MESSAGE("Break");
+ break;
+ }
else if (n < 0){
n++; }
}
- np = k % _mapGrid+2;
- Current_Pnt.Pi = kp;
- Current_Pnt.Pj = np;
+ np = (n + _mapGrid + 1) % (_mapGrid+1);
+ //MESSAGE(" n = "<<n<<" ; np = "<<np);
+ Current_Pnt[0] = kp;
+ Current_Pnt[1] = np;
if (_known.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");
du = (k-i) * (_u2 - _u1) / _mapGrid;
dv = (n-j) * (_v2 - _v1) / _mapGrid;
- Dist = _DMap[kp][np] + Guu * du*du + Gvv * dv*dv; // ds**2 = du'Gdu + dv'Gdv (we assume Guv=Gvu=0)
+ Dist = _DMap[kp][np] + sqrt( Guu * du*du + Gvv * dv*dv ); // ds**2 = du'Gdu + dv'Gdv (we assume Guv=Gvu=0)
+ //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");
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
- _trial.erase(TPnt);
- TPnt.dist = D_Ref;
- TPnt.Ti = i;
- TPnt.Tj = j;
- _trial.insert(TPnt);
+ //MESSAGE("trial contient : "<<_trial.size()<<"elements");
+ //MESSAGE("Dist updated");
+ 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);
+ }
+ 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");
}
+ //MESSAGE("trial contient : "<<_trial.size()<<"elements");
} // if
} // for
} // for
