src/BLSURFPlugin/BLSURFPlugin_Attractor.cxx

   1 // Copyright (C) 2007-2021  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, or (at your option) any later version.
   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 #include <GeomLib_IsPlanarSurface.hxx>
  42
  43 BLSURFPlugin_Attractor::BLSURFPlugin_Attractor ()
  44   : _face(),
  45   _attractorShape(),
  46   _attEntry(),
  47   _vectU(),
  48   _vectV(),
  49   _DMap(),
  50   _known(),
  51   _trial(),
  52   _type(-1),
  53   _gridU(0),
  54   _gridV(0),
  55   _u1 (0.),
  56   _u2 (0.),
  57   _v1 (0.),
  58   _v2 (0.),
  59   _startSize(-1),
  60   _endSize(-1),
  61   _actionRadius(-1),
  62   _constantRadius(-1),
  63   _isMapBuilt(false),
  64   _isEmpty(true){ MESSAGE("construction of a void attractor"); }
  65
  66 BLSURFPlugin_Attractor::BLSURFPlugin_Attractor (const TopoDS_Face& Face, const TopoDS_Shape& Attractor, const std::string& attEntry)
  67   : _face(),
  68     _attractorShape(),
  69     _attEntry(attEntry),
  70     _vectU(),
  71     _vectV(),
  72     _DMap(),
  73     _known(),
  74     _trial(),
  75     _type(0),
  76     _gridU(),
  77     _gridV(),
  78     _u1 (0.),
  79     _u2 (0.),
  80     _v1 (0.),
  81     _v2 (0.),
  82     _startSize(-1),
  83     _endSize(-1),
  84     _actionRadius(-1),
  85     _constantRadius(-1),
  86     _isMapBuilt(false),
  87     _isEmpty(false)
  88 {
  89   _face = Face;
  90   _attractorShape = Attractor;
  91
  92   init();
  93 }
  94
  95 bool BLSURFPlugin_Attractor::init()
  96 {
  97   Standard_Real u0,v0;
  98   int i,j,i0,j0 ;
  99   _known.clear();
 100   _trial.clear();
 101   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(_face);
 102
 103   _distance = &BLSURFPlugin_Attractor::_distanceFromMap;
 104
 105   if ( GeomLib_IsPlanarSurface( aSurf ).IsPlanar() &&
 106        _attractorShape.ShapeType() == TopAbs_VERTEX )
 107   {
 108     // a specific case, the map is not needed
 109     gp_Pnt P = BRep_Tool::Pnt( TopoDS::Vertex( _attractorShape ));
 110     GeomAPI_ProjectPointOnSurf projector( P, aSurf );
 111     if ( projector.IsDone() && projector.NbPoints() == 1 )
 112     {
 113       projector.LowerDistanceParameters(u0,v0);
 114
 115       _attractorPnt = aSurf->Value( u0,v0 );
 116       _plane        = aSurf;
 117       _distance     = &BLSURFPlugin_Attractor::_distanceFromPoint;
 118       _isMapBuilt   = true;
 119       _isEmpty      = false;
 120       return true;
 121     }
 122   }
 123
 124   // Calculation of the bounds of the face
 125   ShapeAnalysis::GetFaceUVBounds(_face,_u1,_u2,_v1,_v2);
 126
 127   _gridU = 300;
 128   _gridV = 300;
 129
 130   for (i=0; i<=_gridU; i++){
 131     _vectU.push_back(_u1+i*(_u2-_u1)/_gridU) ;
 132   }
 133   for (j=0; j<=_gridV; j++){
 134     _vectV.push_back(_v1+j*(_v2-_v1)/_gridV) ;
 135   }
 136
 137   // Initialization of _DMap and _known
 138   std::vector<double> temp(_gridV+1,std::numeric_limits<double>::infinity());  // Set distance of all "far" points to Infinity
 139   for (i=0; i<=_gridU; i++){
 140     _DMap.push_back(temp);
 141   }
 142   std::vector<bool> temp2(_gridV+1,false);
 143   for (i=0; i<=_gridU; i++){
 144     _known.push_back(temp2);
 145   }
 146
 147
 148   // Determination of the starting points
 149   TopExp_Explorer anEdgeExp(_attractorShape, TopAbs_EDGE, TopAbs_FACE);
 150   TopExp_Explorer aVertExp(_attractorShape, TopAbs_VERTEX, TopAbs_EDGE);
 151
 152   for(; anEdgeExp.More(); anEdgeExp.Next()){
 153     const TopoDS_Edge& anEdge = TopoDS::Edge(anEdgeExp.Current());
 154     edgeInit(aSurf, anEdge);
 155   }
 156
 157   for(; aVertExp.More(); aVertExp.Next()){
 158     const TopoDS_Vertex& aVertex = TopoDS::Vertex(aVertExp.Current());
 159     Trial_Pnt TPnt(3,0);
 160     gp_Pnt P = BRep_Tool::Pnt(aVertex);
 161     GeomAPI_ProjectPointOnSurf projector( P, aSurf );
 162     projector.LowerDistanceParameters(u0,v0);
 163     i0 = int( floor ( (u0 - _u1) * _gridU / (_u2 - _u1) + 0.5 ));
 164     j0 = int( floor ( (v0 - _v1) * _gridV / (_v2 - _v1) + 0.5 ));
 165     TPnt[0]=0.;                                                                // Set the distance of the starting point to 0.
 166     TPnt[1]=double( i0 );
 167     TPnt[2]=double( j0 );
 168     _DMap[i0][j0] = 0.;
 169     _trial.insert(TPnt);                                                       // Move starting point to _trial
 170   }
 171
 172   return true;
 173 }
 174
 175 // check that i and j are inside the bounds of the grid to avoid out of bounds errors
 176 // in affectation of the grid's vectors
 177 void BLSURFPlugin_Attractor::avoidOutOfBounds(int& i, int& j)
 178 {
 179   if (i > _gridU)
 180     i = _gridU;
 181   if (i < 0)
 182     i = 0;
 183   if (j > _gridV)
 184     j = _gridV;
 185   if (j < 0)
 186     j = 0;
 187 }
 188
 189 void BLSURFPlugin_Attractor::edgeInit(Handle(Geom_Surface) theSurf, const TopoDS_Edge& anEdge)
 190 {
 191   gp_Pnt2d P2;
 192   double first;
 193   double last;
 194   int i,i0,j0;
 195   Trial_Pnt TPnt(3,0);
 196   Handle(Geom2d_Curve) aCurve2d;
 197   Handle(Geom_Curve) aCurve3d = BRep_Tool::Curve (anEdge, first, last);
 198   ShapeConstruct_ProjectCurveOnSurface curveProjector;
 199   curveProjector.Init(theSurf, Precision::Confusion());
 200   curveProjector.Perform (aCurve3d, first, last, aCurve2d);
 201
 202   int N = 1200;
 203   for (i=0; i<=N; i++){
 204     P2 = aCurve2d->Value(first + i * (last-first) / N);
 205     i0 = int(floor( (P2.X() - _u1) * _gridU / (_u2 - _u1) + 0.5 ));
 206     j0 = int(floor( (P2.Y() - _v1) * _gridV / (_v2 - _v1) + 0.5 ));
 207
 208     // Avoid out of bounds errors when the ends of the edge are outside the face
 209     avoidOutOfBounds(i0, j0);
 210
 211     TPnt[0] = 0.;
 212     TPnt[1] = i0;
 213     TPnt[2] = j0;
 214     _DMap[i0][j0] = 0.;
 215     _trial.insert(TPnt);
 216   }
 217 }
 218
 219
 220 void BLSURFPlugin_Attractor::SetParameters(double Start_Size, double End_Size, double Action_Radius, double Constant_Radius)
 221 {
 222   _startSize = Start_Size;
 223   _endSize = End_Size;
 224   _actionRadius = Action_Radius;
 225   _constantRadius = Constant_Radius;
 226 }
 227
 228 double BLSURFPlugin_Attractor::_distanceFromPoint(double u, double v)
 229 {
 230   return _attractorPnt.Distance( _plane->Value( u, v ));
 231 }
 232
 233 double BLSURFPlugin_Attractor::_distanceFromMap(double u, double v)
 234 {
 235   //   MG-CADSurf seems to perform a linear interpolation so it's sufficient to give it a non-continuous distance map
 236   int i = int(floor ( (u - _u1) * _gridU / (_u2 - _u1) + 0.5 ));
 237   int j = int(floor ( (v - _v1) * _gridV / (_v2 - _v1) + 0.5 ));
 238
 239   // Avoid out of bounds errors in _DMap
 240   avoidOutOfBounds(i, j);
 241
 242   return _DMap[i][j];
 243 }
 244
 245 double BLSURFPlugin_Attractor::GetSize(double u, double v)
 246 {
 247   const double attrDist = (this->*_distance)(u,v);
 248   const double myDist = 0.5 * (attrDist - _constantRadius + fabs(attrDist - _constantRadius));
 249   switch(_type)
 250   {
 251     case TYPE_EXP:
 252       if (fabs(_actionRadius) <= std::numeric_limits<double>::epsilon()){
 253         if (myDist <= std::numeric_limits<double>::epsilon()){
 254           return _startSize;
 255         }
 256         else {
 257           return _endSize;
 258         }
 259       }
 260       else{
 261         return _endSize - (_endSize - _startSize) * exp(- myDist * myDist / (_actionRadius * _actionRadius) );
 262       }
 263       break;
 264     case TYPE_LIN:
 265         return _startSize + ( 0.5 * (attrDist - _constantRadius + fabs(attrDist - _constantRadius)) ) ;
 266       break;
 267   }
 268   return -1;
 269 }
 270
 271
 272 void BLSURFPlugin_Attractor::BuildMap()
 273 {
 274   MESSAGE("building the map");
 275   int i, j, k, n;
 276   //int count = 0;
 277   int ip, jp, kp, np;
 278   int i0, j0;
 279   gp_Pnt P;
 280   gp_Vec D1U,D1V;
 281   double Guu, Gvv, Guv;         // Components of the local metric tensor
 282   double du, dv;
 283   double D_Ref = 0.;
 284   double Dist = 0.;
 285   bool Dist_changed;
 286   IJ_Pnt Current_Pnt(2,0);
 287   Trial_Pnt TPnt(3,0);
 288   TTrialSet::iterator min;
 289   TTrialSet::iterator found;
 290   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(_face);
 291
 292   // While there are points in "Trial" (representing a kind of advancing front), loop on them -----------------------------------------------------------
 293   while (_trial.size() > 0 ) {
 294     min = _trial.begin();                        // Get trial point with min distance from start
 295     i0 = int( (*min)[1] );
 296     j0 = int( (*min)[2] );
 297     // Avoid out of bounds errors in _known affectations
 298     avoidOutOfBounds(i0, j0);
 299     _known[i0][j0] = true;                       // Move it to "Known"
 300     _trial.erase(min);                           // Remove it from "Trial"
 301
 302     // Loop on neighbours of the trial min --------------------------------------------------------------------------------------------------------------
 303     for (i=i0 - 1 ; i <= i0 + 1 ; i++){
 304       if (!aSurf->IsUPeriodic()){                          // Periodic conditions in U
 305         if (i > _gridU ){
 306           break; }
 307         else if (i < 0){
 308           i++; }
 309       }
 310       ip = (i + _gridU + 1) % (_gridU+1);                  // We get a periodic index :
 311       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;
 312         if (!aSurf->IsVPeriodic()){                        // Periodic conditions in V .
 313           if (j > _gridV ){
 314             break; }
 315           else if (j < 0){
 316             j++;
 317           }
 318         }
 319         jp = (j + _gridV + 1) % (_gridV+1);
 320
 321         if (!_known[ip][jp]){                              // If the distance is not known yet
 322           aSurf->D1(_vectU[ip],_vectV[jp],P,D1U,D1V);      // Calculate the metric tensor at (i,j)
 323           // G(i,j)  =  | ||dS/du||**2          *     |
 324           //            | <dS/du,dS/dv>  ||dS/dv||**2 |
 325           Guu = D1U.X()*D1U.X() +  D1U.Y()*D1U.Y() + D1U.Z()*D1U.Z();    // Guu = ||dS/du||**2
 326           Gvv = D1V.X()*D1V.X() +  D1V.Y()*D1V.Y() + D1V.Z()*D1V.Z();    // Gvv = ||dS/dv||**2
 327           Guv = D1U.X()*D1V.X() +  D1U.Y()*D1V.Y() + D1U.Z()*D1V.Z();    // Guv = Gvu = < dS/du,dS/dv >
 328           D_Ref = _DMap[ip][jp];                           // Set a ref. distance of the point to its value in _DMap
 329           TPnt[0] = D_Ref;                                 // (may be infinite or uncertain)
 330           TPnt[1] = ip;
 331           TPnt[2] = jp;
 332           Dist_changed = false;
 333
 334           // Loop on neighbours to calculate the min distance from them ---------------------------------------------------------------------------------
 335           for (k=i - 1 ; k <= i + 1 ; k++){
 336             if (!aSurf->IsUPeriodic()){                              // Periodic conditions in U
 337               if(k > _gridU ){
 338                 break;
 339               }
 340               else if (k < 0){
 341                 k++; }
 342             }
 343             kp = (k + _gridU + 1) % (_gridU+1);                      // periodic index
 344             for (n=j - 1 ; n <= j + 1 ; n++){
 345               if (!aSurf->IsVPeriodic()){                            // Periodic conditions in V
 346                 if(n > _gridV){
 347                   break;
 348                 }
 349                 else if (n < 0){
 350                   n++; }
 351               }
 352               np = (n + _gridV + 1) % (_gridV+1);
 353               if (_known[kp][np]){                                   // If the distance of the neighbour is known
 354                                                                      // Calculate the distance from (k,n)
 355                 du = (k-i) * (_u2 - _u1) / _gridU;
 356                 dv = (n-j) * (_v2 - _v1) / _gridV;
 357                 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)
 358                 if (Dist < D_Ref) {                                  // If smaller than ref. distance  ->  update ref. distance
 359                   D_Ref = Dist;
 360                   Dist_changed = true;
 361                 }
 362               }
 363             }
 364           } // End of the loop on neighbours --------------------------------------------------------------------------------------------------------------
 365
 366           if (Dist_changed) {                              // If distance has been updated, update _trial
 367             found=_trial.find(TPnt);
 368             if (found != _trial.end()){
 369               _trial.erase(found);                         // Erase the point if it was already in _trial
 370             }
 371             TPnt[0] = D_Ref;
 372             TPnt[1] = ip;
 373             TPnt[2] = jp;
 374             _DMap[ip][jp] = D_Ref;                         // Set it distance to the minimum distance found during the loop above
 375             _trial.insert(TPnt);                           // Insert it (or reinsert it) in _trial
 376           }
 377         } // end if (!_known[ip][jp])
 378       } // for
 379     } // for
 380   } // while (_trial)
 381   _known.clear();
 382   _trial.clear();
 383   _isMapBuilt = true;
 384 } // end of BuildMap()
 385
 386
 387