src/BLSURFPlugin/BLSURFPlugin_Attractor.cxx

   1 //  Copyright (C) 2007-2010  CEA/DEN, EDF R&D
   2 //
   3 //  This library is free software; you can redistribute it and/or
   4 //  modify it under the terms of the GNU Lesser General Public
   5 //  License as published by the Free Software Foundation; either
   6 //  version 2.1 of the License.
   7 //
   8 //  This library is distributed in the hope that it will be useful,
   9 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 //  Lesser General Public License for more details.
  12 //
  13 //  You should have received a copy of the GNU Lesser General Public
  14 //  License along with this library; if not, write to the Free Software
  15 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 //  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 // ---
  21 // File    : BLSURFPlugin_Attractor.cxx
  22 // Authors : Renaud Nédélec (OCC)
  23 // ---
  24 //
  25 // The idea of the algorithm used to calculate the distance on a
  26 // non-euclidian parametric surface has been found in the ref. below:
  27 //
  28 // Ref:"Accurate Anisotropic Fast Marching for Diffusion-Based Geodesic Tractography"
  29 // S. Jbabdi, P. Bellec, R. Toro, Daunizeau, M. Pélégrini-Issac, and H. Benali1
  30 //
  31
  32 #include "BLSURFPlugin_Attractor.hxx"
  33 #include <utilities.h>
  34 #include <algorithm>
  35 #include <cmath>
  36
  37 // cascade include
  38 #include "ShapeAnalysis.hxx"
  39 #include "ShapeConstruct_ProjectCurveOnSurface.hxx"
  40 #include <Precision.hxx>
  41
  42 BLSURFPlugin_Attractor::BLSURFPlugin_Attractor ()
  43   : _face(),
  44   _attractorShape(),
  45   _attEntry(),
  46   _gridU(0),
  47   _gridV(0),
  48   _vectU(),
  49   _vectV(),
  50   _DMap(),
  51   _known(),
  52   _trial(),
  53   _u1 (0.),
  54   _u2 (0.),
  55   _v1 (0.),
  56   _v2 (0.),
  57   _startSize(-1),
  58   _endSize(-1),
  59   _actionRadius(-1),
  60   _constantRadius(-1),
  61   _type(-1),
  62   _isMapBuilt(false),
  63   _isEmpty(true){ MESSAGE("construction of a void attractor"); }
  64
  65 BLSURFPlugin_Attractor::BLSURFPlugin_Attractor (const TopoDS_Face& Face, const TopoDS_Shape& Attractor, const std::string& attEntry)
  66   : _face(),
  67   _attractorShape(),
  68   _attEntry(attEntry),
  69   _gridU(),
  70   _gridV(),
  71   _vectU(),
  72   _vectV(),
  73   _DMap(),
  74   _known(),
  75   _trial(),
  76   _u1 (0.),
  77   _u2 (0.),
  78   _v1 (0.),
  79   _v2 (0.),
  80   _startSize(-1),
  81   _endSize(-1),
  82   _actionRadius(-1),
  83   _constantRadius(-1),
  84   _type(0),
  85   _isMapBuilt(false),
  86   _isEmpty(false)
  87 {
  88   _face = Face;
  89   _attractorShape = Attractor;
  90
  91   init();
  92 }
  93
  94 bool BLSURFPlugin_Attractor::init(){
  95   Standard_Real u0,v0;
  96   int i,j,i0,j0 ;
  97   _known.clear();
  98   _trial.clear();
  99   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(_face);
 100
 101   // Calculation of the bounds of the face
 102   ShapeAnalysis::GetFaceUVBounds(_face,_u1,_u2,_v1,_v2);
 103   MESSAGE("u1 = "<<_u1<<" ,u2  = "<<_u2);
 104   MESSAGE("v1 = "<<_v1<<" ,v2  = "<<_v2);
 105 //   _gridU = floor (_u2 - _u1) / _step;
 106 //   _gridV = floor (_v2 - _v1) / _step;
 107
 108   _gridU = 300;
 109   _gridV = 300;
 110
 111   for (i=0; i<=_gridU; i++){
 112     _vectU.push_back(_u1+i*(_u2-_u1)/_gridU) ;
 113   }
 114   for (j=0; j<=_gridV; j++){
 115     _vectV.push_back(_v1+j*(_v2-_v1)/_gridV) ;
 116   }
 117
 118   // Initialization of _DMap and _known
 119   std::vector<double> temp(_gridV+1,std::numeric_limits<double>::infinity());  // Set distance of all "far" points to Infinity
 120   for (i=0; i<=_gridU; i++){
 121     _DMap.push_back(temp);
 122   }
 123   std::vector<bool> temp2(_gridV+1,false);
 124   for (i=0; i<=_gridU; i++){
 125     _known.push_back(temp2);
 126   }
 127
 128
 129   // Determination of the starting points
 130   TopExp_Explorer anEdgeExp(_attractorShape, TopAbs_EDGE, TopAbs_FACE);
 131   TopExp_Explorer aVertExp(_attractorShape, TopAbs_VERTEX, TopAbs_EDGE);
 132
 133   for(; anEdgeExp.More(); anEdgeExp.Next()){
 134     const TopoDS_Edge& anEdge = TopoDS::Edge(anEdgeExp.Current());
 135     edgeInit(aSurf, anEdge);
 136   }
 137
 138   for(; aVertExp.More(); aVertExp.Next()){
 139     const TopoDS_Vertex& aVertex = TopoDS::Vertex(aVertExp.Current());
 140     Trial_Pnt TPnt(3,0);
 141     gp_Pnt P = BRep_Tool::Pnt(aVertex);
 142     GeomAPI_ProjectPointOnSurf projector( P, aSurf );
 143     projector.LowerDistanceParameters(u0,v0);
 144     MESSAGE("u0 = "<<u0<<" ,v0  = "<<v0);
 145     i0 = floor ( (u0 - _u1) * _gridU / (_u2 - _u1) + 0.5 );
 146     j0 = floor ( (v0 - _v1) * _gridV / (_v2 - _v1) + 0.5 );
 147     TPnt[0]=0.;                                                                // Set the distance of the starting point to 0.
 148     TPnt[1]=i0;
 149     TPnt[2]=j0;
 150     _DMap[i0][j0] = 0.;
 151     _trial.insert(TPnt);         // Move starting point to _trial
 152   }
 153
 154 }
 155
 156 void BLSURFPlugin_Attractor::edgeInit(Handle(Geom_Surface) theSurf, const TopoDS_Edge& anEdge){
 157   gp_Pnt2d P2;
 158   double first;
 159   double last;
 160   int i,j,i0,j0;
 161   Trial_Pnt TPnt(3,0);
 162   Handle(Geom2d_Curve) aCurve2d;
 163   Handle(Geom_Curve) aCurve3d = BRep_Tool::Curve (anEdge, first, last);
 164   ShapeConstruct_ProjectCurveOnSurface curveProjector;
 165   curveProjector.Init(theSurf, Precision::Confusion());
 166   curveProjector.PerformAdvanced (aCurve3d, first, last, aCurve2d);
 167   //int N = 20 * (last - first) / _step;  // If the edge is a circle : 4>Pi so the number of points on the edge should be good -> 5 for ellipses
 168   int N = 1200;
 169   MESSAGE("Initialization of the starting points")
 170   for (i=0; i<=N; i++){
 171     P2 = aCurve2d->Value(first + i * (last-first) / N);
 172     i0 = floor( (P2.X() - _u1) * _gridU / (_u2 - _u1) + 0.5 );
 173     j0 = floor( (P2.Y() - _v1) * _gridV / (_v2 - _v1) + 0.5 );
 174 //     MESSAGE("i0 = "<<i0<<" , j0 = "<<j0)
 175     TPnt[0] = 0.;
 176     TPnt[1] = i0;
 177     TPnt[2] = j0;
 178     _DMap[i0][j0] = 0.;
 179     _trial.insert(TPnt);
 180   }
 181   MESSAGE("_trial.size() = "<<_trial.size())
 182 }
 183
 184
 185 void BLSURFPlugin_Attractor::SetParameters(double Start_Size, double End_Size, double Action_Radius, double Constant_Radius){
 186   MESSAGE("BLSURFPlugin_Attractor::SetParameters")
 187   _startSize = Start_Size;
 188   _endSize = End_Size;
 189   _actionRadius = Action_Radius;
 190   _constantRadius = Constant_Radius;
 191 }
 192
 193 double BLSURFPlugin_Attractor::_distance(double u, double v){
 194
 195 //   // calcul des coins du carre
 196 //   double stepU = (_u2 - _u1) / _gridU ;
 197 //   double stepV = (_v2 - _v1) / _gridV ;
 198 //
 199 //   int i_sup =  floor ( fabs(u - _u1) / stepU ) + 1;
 200 //   int j_sup =  floor ( fabs(v - _v1) / stepV ) + 1;
 201 //   i_sup = std::min( i_sup, _gridU);
 202 //   j_sup = std::min( j_sup, _gridV);
 203 //
 204 //
 205 // //   int i_inf = floor ( (u - _u1) * _mapGrid / (_u2 - _u1) );
 206 // //   int j_inf = floor ( (v - _v1) * _mapGrid / (_v2 - _v1) );
 207 //   int i_inf = i_sup - 1;
 208 //   int j_inf = j_sup - 1;
 209 //
 210 //   double u_sup = _vectU[i_sup];
 211 //   double v_sup = _vectV[j_sup];
 212 //
 213 //   double u_inf = _vectU[i_inf];
 214 //   double v_inf = _vectV[j_inf];
 215 // //
 216 // // //   MESSAGE("i_inf , i_sup, j_inf, j_sup = "<<i_inf<<" , "<<i_sup<<" , "<<j_inf<<" , "<<j_sup)
 217 // // //   MESSAGE("u = "<<u<<", _u1 ="<<_u1)
 218 // // //   MESSAGE("(u - _u1) / stepU = "<< (u - _u1) / stepU)
 219 //   double d1 = _DMap[i_sup][j_sup];
 220 //   double d2 = _DMap[i_sup][j_inf];
 221 //   double d3 = _DMap[i_inf][j_sup];
 222 //   double d4 = _DMap[i_inf][j_inf];
 223 // //
 224 // //   double d = 0.25 * (d1 + d2 + d3 + d4);
 225 // //   //double d = d1;
 226 // //  //
 227 // // //   if (fabs(v_inf-v_sup) < 1e-14 || fabs(u_inf-u_sup) < 1e-14 ){
 228 // // //     MESSAGE("division par zero v_inf-v_sup = "<< fabs(v_inf-v_sup)<<" , u_inf-u_sup"<<fabs(u_inf-u_sup))
 229 // // //     MESSAGE("v_inf = "<< v_inf<<" , v_sup"<<v_sup)
 230 // // //     MESSAGE("u_inf = "<< u_inf<<" , u_sup"<<u_sup)
 231 // // //   }
 232 //   double P_inf =   d4 * ( 1 + (u_inf - u) / stepU + (v_inf - v) / stepV )
 233 //                  + d3 * ( (v - v_inf) / stepV )
 234 //                  + d2 * ( (u - u_inf) / stepU ) ;
 235 //
 236 //   double P_sup =   d1 * ( 1 + (u - u_sup) / stepU + (v - v_sup) / stepV )
 237 //                  + d3 * ( (u_sup - u) / stepU )
 238 //                  + d2 * ( (v_sup - v) / stepV ) ;
 239 //
 240 //   // calcul de la distance (interpolation lineaire)
 241 //   bool P_switch = (u+v > u_sup+v_sup);
 242 //   double d;
 243 //   if (P_switch){
 244 //     d = P_inf;
 245 //   }
 246 //   else {
 247 //     d = P_sup;
 248 //   }
 249 //
 250 //   return d;
 251
 252   //   BLSURF seems to perform a linear interpolation so it's sufficient to give it a non-continuous distance map
 253   int i = floor ( (u - _u1) * _gridU / (_u2 - _u1) + 0.5 );
 254   int j = floor ( (v - _v1) * _gridV / (_v2 - _v1) + 0.5 );
 255
 256   return _DMap[i][j];
 257 }
 258
 259
 260 double BLSURFPlugin_Attractor::GetSize(double u, double v){
 261   double myDist = 0.5 * (_distance(u,v) - _constantRadius + fabs(_distance(u,v) - _constantRadius));
 262 //   if (myDist<0){
 263 //     MESSAGE("Warning myDist<0 : myDist= "<<myDist)
 264 //   }
 265   switch(_type)
 266   {
 267     case TYPE_EXP:
 268       if (fabs(_actionRadius) <= std::numeric_limits<double>::epsilon()){
 269         if (myDist <= std::numeric_limits<double>::epsilon()){
 270           return _startSize;
 271         }
 272         else {
 273           return _endSize;
 274         }
 275       }
 276       else{
 277         return _endSize - (_endSize - _startSize) * exp(- myDist * myDist / (_actionRadius * _actionRadius) );
 278       }
 279       break;
 280     case TYPE_LIN:
 281         return _startSize + ( 0.5 * (_distance(u,v) - _constantRadius + abs(_distance(u,v) - _constantRadius)) ) ;
 282       break;
 283   }
 284 }
 285
 286
 287 void BLSURFPlugin_Attractor::BuildMap(){
 288
 289   MESSAGE("building the map");
 290   int i, j, k, n;
 291   int count = 0;
 292   int ip, jp, kp, np;
 293   int i0, j0;
 294   gp_Pnt P;
 295   gp_Vec D1U,D1V;
 296   double Guu, Gvv, Guv;         // Components of the local metric tensor
 297   double du, dv;
 298   double D_Ref = 0.;
 299   double Dist = 0.;
 300   bool Dist_changed;
 301   IJ_Pnt Current_Pnt(2,0);
 302   Trial_Pnt TPnt(3,0);
 303   TTrialSet::iterator min;
 304   TTrialSet::iterator found;
 305   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(_face);
 306
 307   // While there are points in "Trial" (representing a kind of advancing front), loop on them -----------------------------------------------------------
 308   while (_trial.size() > 0 ){
 309     min = _trial.begin();                        // Get trial point with min distance from start
 310     i0 = (*min)[1];
 311     j0 = (*min)[2];
 312     _known[i0][j0] = true;                       // Move it to "Known"
 313 //     MESSAGE("_DMap["<<i0<<"]["<<j0<<"] = "<<_DMap[i0][j0])
 314     _trial.erase(min);                           // Remove it from "Trial"
 315     // Loop on neighbours of the trial min --------------------------------------------------------------------------------------------------------------
 316     for (i=i0 - 1 ; i <= i0 + 1 ; i++){
 317       if (!aSurf->IsUPeriodic()){                          // Periodic conditions in U
 318         if (i > _gridU ){
 319           break; }
 320         else if (i < 0){
 321           i++; }
 322       }
 323       ip = (i + _gridU + 1) % (_gridU+1);                  // We get a periodic index :
 324       for (j=j0 - 1 ; j <= j0 + 1 ; j++){                  //    ip=modulo(i,N+2) so that  i=-1->ip=N; i=0 -> ip=0 ; ... ; i=N+1 -> ip=0;
 325         if (!aSurf->IsVPeriodic()){                        // Periodic conditions in V .
 326           if (j > _gridV ){
 327             break; }
 328           else if (j < 0){
 329             j++;
 330           }
 331         }
 332         jp = (j + _gridV + 1) % (_gridV+1);
 333
 334         if (!_known[ip][jp]){                              // If the distance is not known yet
 335           aSurf->D1(_vectU[ip],_vectV[jp],P,D1U,D1V);      // Calculate the metric at (i,j)
 336           // G(i,j)  =  | ||dS/du||**2          *     |
 337           //            | <dS/du,dS/dv>  ||dS/dv||**2 |
 338           Guu = D1U.X()*D1U.X() +  D1U.Y()*D1U.Y() + D1U.Z()*D1U.Z();    // Guu = ||dS/du||**2
 339           Gvv = D1V.X()*D1V.X() +  D1V.Y()*D1V.Y() + D1V.Z()*D1V.Z();    // Gvv = ||dS/dv||**2
 340           Guv = D1U.X()*D1V.X() +  D1U.Y()*D1V.Y() + D1U.Z()*D1V.Z();    // Guv = Gvu = < dS/du,dS/dv >
 341           D_Ref = _DMap[ip][jp];                           // Set a ref. distance of the point to its value in _DMap
 342           TPnt[0] = D_Ref;                                 // (may be infinite or uncertain)
 343           TPnt[1] = ip;
 344           TPnt[2] = jp;
 345           Dist_changed = false;
 346           // Loop on neighbours to calculate the min distance from them ---------------------------------------------------------------------------------
 347           for (k=i - 1 ; k <= i + 1 ; k++){
 348             if (!aSurf->IsUPeriodic()){                              // Periodic conditions in U
 349               if(k > _gridU ){
 350                 break;
 351               }
 352               else if (k < 0){
 353                 k++; }
 354             }
 355             kp = (k + _gridU + 1) % (_gridU+1);                      // periodic index
 356             for (n=j - 1 ; n <= j + 1 ; n++){
 357               if (!aSurf->IsVPeriodic()){                            // Periodic conditions in V
 358                 if(n > _gridV){
 359                   break;
 360                 }
 361                 else if (n < 0){
 362                   n++; }
 363               }
 364               np = (n + _gridV + 1) % (_gridV+1);
 365               if (_known[kp][np]){                                   // If the distance of the neighbour is known
 366                                                                      // Calculate the distance from (k,n)
 367                 du = (k-i) * (_u2 - _u1) / _gridU;
 368                 dv = (n-j) * (_v2 - _v1) / _gridV;
 369                 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)
 370                 if (Dist < D_Ref) {                                  // If smaller than ref. distance  ->  update ref. distance
 371                   D_Ref = Dist;
 372                   Dist_changed = true;
 373                 }
 374               }
 375             }
 376           } // End of the loop on neighbours --------------------------------------------------------------------------------------------------------------
 377           if (Dist_changed) {                              // If distance has been updated, update _trial
 378             found=_trial.find(TPnt);
 379             if (found != _trial.end()){
 380               _trial.erase(found);                         // Erase the point if it was already in _trial
 381             }
 382             TPnt[0] = D_Ref;
 383             TPnt[1] = ip;
 384             TPnt[2] = jp;
 385             _DMap[ip][jp] = D_Ref;                         // Set it distance to the minimum distance found during the loop above
 386             _trial.insert(TPnt);                           // Insert it (or reinsert it) in _trial
 387           }
 388         } // end if (!_known[ip][jp])
 389       } // for
 390     } // for
 391   } // while (_trial)
 392   _known.clear();
 393   _trial.clear();
 394   _isMapBuilt = true;
 395   MESSAGE("_gridU = "<<_gridU<<" , _gridV = "<<_gridV)
 396 } // end of BuildMap()
 397
 398
 399