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