+// Copyright (C) 2007-2010 CEA/DEN, EDF R&D
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
+//
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+
+// ---
+// File : BLSURFPlugin_Attractor.cxx
+// Authors : Renaud Nédélec (OCC)
+// ---
+//
+// The idea of the algorithm used to calculate the distance on a
+// non-euclidian parametric surface has been found in the ref. below:
+//
+// Ref:"Accurate Anisotropic Fast Marching for Diffusion-Based Geodesic Tractography"
+// S. Jbabdi, P. Bellec, R. Toro, Daunizeau, M. Pélégrini-Issac, and H. Benali1
+//
+
#include "BLSURFPlugin_Attractor.hxx"
#include <utilities.h>
#include <algorithm>
_isMapBuilt(false),
_isEmpty(true){ MESSAGE("construction of a void attractor"); }
-BLSURFPlugin_Attractor::BLSURFPlugin_Attractor (const TopoDS_Face& Face, const TopoDS_Shape& Attractor, const std::string& attEntry, double Step)
+BLSURFPlugin_Attractor::BLSURFPlugin_Attractor (const TopoDS_Face& Face, const TopoDS_Shape& Attractor, const std::string& attEntry, double Step) // TODO Step is now unused -> remove it if testing is OK
: _face(),
_attractorShape(),
_attEntry(attEntry),
- _step(Step),
+ _step(),
_gridU(),
_gridV(),
_vectU(),
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;
+// _gridU = floor (_u2 - _u1) / _step;
+// _gridV = floor (_v2 - _v1) / _step;
+ // TEST
+ _gridU = 300;
+ _gridV = 300;
+ _step = std::min((_u2-_u1)/_gridU,(_v2-_v1)/_gridV);
+
for (i=0; i<=_gridU; i++){
_vectU.push_back(_u1+i*(_u2-_u1)/_gridU) ;
}
ShapeConstruct_ProjectCurveOnSurface curveProjector;
curveProjector.Init(aSurf, Precision::Confusion());
curveProjector.PerformAdvanced (aCurve3d, first, last, aCurve2d);
- 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
+ // TEST
+ //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
+ int N = 1200;
+ MESSAGE("Initialisation des points de départ")
for (i=0; i<=N; i++){
P2 = aCurve2d->Value(first + i * (last-first) / N);
i0 = floor( (P2.X() - _u1) * _gridU / (_u2 - _u1) + 0.5 );
//
// double u_inf = _vectU[i_inf];
// double v_inf = _vectV[j_inf];
-//
-// // MESSAGE("i_inf , i_sup, j_inf, j_sup = "<<i_inf<<" , "<<i_sup<<" , "<<j_inf<<" , "<<j_sup)
-// // MESSAGE("u = "<<u<<", _u1 ="<<_u1)
-// // MESSAGE("(u - _u1) / stepU = "<< (u - _u1) / stepU)
+// //
+// // // MESSAGE("i_inf , i_sup, j_inf, j_sup = "<<i_inf<<" , "<<i_sup<<" , "<<j_inf<<" , "<<j_sup)
+// // // MESSAGE("u = "<<u<<", _u1 ="<<_u1)
+// // // MESSAGE("(u - _u1) / stepU = "<< (u - _u1) / stepU)
// double d1 = _DMap[i_sup][j_sup];
// double d2 = _DMap[i_sup][j_inf];
// double d3 = _DMap[i_inf][j_sup];
// double d4 = _DMap[i_inf][j_inf];
-//
-// //double d = 0.25 * (d1 + d2 + d3 + d4);
-// //double d = d1;
-// //
-// // if (fabs(v_inf-v_sup) < 1e-14 || fabs(u_inf-u_sup) < 1e-14 ){
-// // MESSAGE("division par zero v_inf-v_sup = "<< fabs(v_inf-v_sup)<<" , u_inf-u_sup"<<fabs(u_inf-u_sup))
-// // MESSAGE("v_inf = "<< v_inf<<" , v_sup"<<v_sup)
-// // MESSAGE("u_inf = "<< u_inf<<" , u_sup"<<u_sup)
-// // }
+// //
+// // double d = 0.25 * (d1 + d2 + d3 + d4);
+// // //double d = d1;
+// // //
+// // // if (fabs(v_inf-v_sup) < 1e-14 || fabs(u_inf-u_sup) < 1e-14 ){
+// // // MESSAGE("division par zero v_inf-v_sup = "<< fabs(v_inf-v_sup)<<" , u_inf-u_sup"<<fabs(u_inf-u_sup))
+// // // MESSAGE("v_inf = "<< v_inf<<" , v_sup"<<v_sup)
+// // // MESSAGE("u_inf = "<< u_inf<<" , u_sup"<<u_sup)
+// // // }
// double P_inf = d4 * ( 1 + (u_inf - u) / stepU + (v_inf - v) / stepV )
// + d3 * ( (v - v_inf) / stepV )
// + d2 * ( (u - u_inf) / stepU ) ;
//
// return d;
- // BLSURF seems to perform a linear interpolation so it's sufficient to give it a non-continuous distance map
+ // 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 );