From a02d4651d47c998d20ab26462c12a4eaba31fd2d Mon Sep 17 00:00:00 2001 From: gdd Date: Fri, 1 Apr 2011 14:14:05 +0000 Subject: [PATCH] Merge from BR_ATTRACTOR_GEO --- src/BLSURFPlugin/BLSURFPlugin_Attractor.cxx | 387 ++++++++++++++++++++ src/BLSURFPlugin/BLSURFPlugin_Attractor.hxx | 144 ++++++++ 2 files changed, 531 insertions(+) create mode 100644 src/BLSURFPlugin/BLSURFPlugin_Attractor.cxx create mode 100644 src/BLSURFPlugin/BLSURFPlugin_Attractor.hxx diff --git a/src/BLSURFPlugin/BLSURFPlugin_Attractor.cxx b/src/BLSURFPlugin/BLSURFPlugin_Attractor.cxx new file mode 100644 index 0000000..0f9478f --- /dev/null +++ b/src/BLSURFPlugin/BLSURFPlugin_Attractor.cxx @@ -0,0 +1,387 @@ +// 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 +#include +#include + +// cascade include +#include "ShapeAnalysis.hxx" +#include "ShapeConstruct_ProjectCurveOnSurface.hxx" +#include + +BLSURFPlugin_Attractor::BLSURFPlugin_Attractor () + : _face(), + _attractorShape(), + _attEntry(), + _step(0), + _gridU(0), + _gridV(0), + _vectU(), + _vectV(), + _DMap(), + _known(), + _trial(), + _u1 (0.), + _u2 (0.), + _v1 (0.), + _v2 (0.), + _startSize(-1), + _endSize(-1), + _actionRadius(-1), + _constantRadius(-1), + _type(-1), + _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) // TODO Step is now unused -> remove it if testing is OK + : _face(), + _attractorShape(), + _attEntry(attEntry), + _step(), + _gridU(), + _gridV(), + _vectU(), + _vectV(), + _DMap(), + _known(), + _trial(), + _u1 (0.), + _u2 (0.), + _v1 (0.), + _v2 (0.), + _startSize(-1), + _endSize(-1), + _actionRadius(-1), + _constantRadius(-1), + _type(0), + _isMapBuilt(false), + _isEmpty(false) +{ + _face = Face; + _attractorShape = Attractor; + + init(); +} + +bool BLSURFPlugin_Attractor::init(){ + Standard_Real u0,v0; + int i,j,i0,j0 ; + _known.clear(); + _trial.clear(); + Handle(Geom_Surface) aSurf = BRep_Tool::Surface(_face); + Trial_Pnt TPnt(3,0); + + // Calculation of the bounds of the face + 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; + // 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) ; + } + for (j=0; j<=_gridV; j++){ + _vectV.push_back(_v1+j*(_v2-_v1)/_gridV) ; + } + + // Initialization of _DMap and _known + std::vector temp(_gridV+1,std::numeric_limits::infinity()); // Set distance of all "far" points to Infinity + for (i=0; i<=_gridU; i++){ + _DMap.push_back(temp); + } + std::vector temp2(_gridV+1,false); + for (i=0; i<=_gridU; i++){ + _known.push_back(temp2); + } + + // Determination of the starting points + if (_attractorShape.ShapeType() == TopAbs_VERTEX){ + gp_Pnt P = BRep_Tool::Pnt(TopoDS::Vertex(_attractorShape)); + GeomAPI_ProjectPointOnSurf projector( P, aSurf ); + projector.LowerDistanceParameters(u0,v0); + MESSAGE("u0 = "<Value(first + i * (last-first) / N); + i0 = floor( (P2.X() - _u1) * _gridU / (_u2 - _u1) + 0.5 ); + j0 = floor( (P2.Y() - _v1) * _gridV / (_v2 - _v1) + 0.5 ); + //MESSAGE("i0 = "< u_sup+v_sup); +// double d; +// if (P_switch){ +// d = P_inf; +// } +// else { +// d = P_sup; +// } +// +// return d; + + // 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 ); + + return _DMap[i][j]; + +} + +double BLSURFPlugin_Attractor::GetSize(double u, double v){ + double myDist = 0.5 * (_distance(u,v) - _constantRadius + fabs(_distance(u,v) - _constantRadius)); + if (myDist<0){ + MESSAGE("Warning myDist<0 : myDist= "< 0 ){ + min = _trial.begin(); // Get trial point with min distance from start + i0 = (*min)[1]; + j0 = (*min)[2]; + _known[i0][j0] = true; // Move it to "Known" + _trial.erase(min); // Remove it from "Trial" + // Loop on neighbours of the trial min -------------------------------------------------------------------------------------------------------------- + for (i=i0 - 1 ; i <= i0 + 1 ; i++){ + if (!aSurf->IsUPeriodic()){ // Periodic conditions in U + if (i > _gridU ){ + break; } + else if (i < 0){ + i++; } + } + ip = (i + _gridU + 1) % (_gridU+1); // We get a periodic index : + 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; + if (!aSurf->IsVPeriodic()){ // Periodic conditions in V . + if (j > _gridV ){ + break; } + else if (j < 0){ + j++; + } + } + jp = (j + _gridV + 1) % (_gridV+1); + + if (!_known[ip][jp]){ // If the distance is not known yet + aSurf->D1(_vectU[ip],_vectV[jp],P,D1U,D1V); // Calculate the metric at (i,j) + // G(i,j) = | ||dS/du||**2 * | + // | ||dS/dv||**2 | + Guu = D1U.X()*D1U.X() + D1U.Y()*D1U.Y() + D1U.Z()*D1U.Z(); // Guu = ||dS/du||**2 + Gvv = D1V.X()*D1V.X() + D1V.Y()*D1V.Y() + D1V.Z()*D1V.Z(); // Gvv = ||dS/dv||**2 + Guv = D1U.X()*D1V.X() + D1U.Y()*D1V.Y() + D1U.Z()*D1V.Z(); // Guv = Gvu = < dS/du,dS/dv > + D_Ref = _DMap[ip][jp]; // Set a ref. distance of the point to its value in _DMap + TPnt[0] = D_Ref; // (may be infinite or uncertain) + TPnt[1] = ip; + TPnt[2] = jp; + Dist_changed = false; + // Loop on neighbours to calculate the min distance from them --------------------------------------------------------------------------------- + for (k=i - 1 ; k <= i + 1 ; k++){ + if (!aSurf->IsUPeriodic()){ // Periodic conditions in U + if(k > _gridU ){ + break; + } + else if (k < 0){ + k++; } + } + kp = (k + _gridU + 1) % (_gridU+1); // periodic index + for (n=j - 1 ; n <= j + 1 ; n++){ + if (!aSurf->IsVPeriodic()){ // Periodic conditions in V + if(n > _gridV){ + break; + } + else if (n < 0){ + n++; } + } + np = (n + _gridV + 1) % (_gridV+1); + if (_known[kp][np]){ // If the distance of the neighbour is known + // Calculate the distance from (k,n) + du = (k-i) * (_u2 - _u1) / _gridU; + dv = (n-j) * (_v2 - _v1) / _gridV; + 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) + if (Dist < D_Ref) { // If smaller than ref. distance -> update ref. distance + D_Ref = Dist; + Dist_changed = true; + } + } + } + } // End of the loop on neighbours -------------------------------------------------------------------------------------------------------------- + if (Dist_changed) { // If distance has been updated, update _trial + found=_trial.find(TPnt); + if (found != _trial.end()){ + _trial.erase(found); // Erase the point if it was already in _trial + } + TPnt[0] = D_Ref; + TPnt[1] = ip; + TPnt[2] = jp; + _DMap[ip][jp] = D_Ref; // Set it distance to the minimum distance found during the loop above + _trial.insert(TPnt); // Insert it (or reinsert it) in _trial + } + } // if + } // for + } // for + } // while (_trial) + _known.clear(); + _trial.clear(); + _isMapBuilt = true; + MESSAGE("_gridU = "<<_gridU<<" , _gridV = "<<_gridV) +} // end of BuildMap() + + + diff --git a/src/BLSURFPlugin/BLSURFPlugin_Attractor.hxx b/src/BLSURFPlugin/BLSURFPlugin_Attractor.hxx new file mode 100644 index 0000000..74260d6 --- /dev/null +++ b/src/BLSURFPlugin/BLSURFPlugin_Attractor.hxx @@ -0,0 +1,144 @@ +// 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.hxx +// 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 +// + +#ifndef _BLSURFPlugin_Attractor_HXX_ +#define _BLSURFPlugin_Attractor_HXX_ + +#include +#include +#include +#include +#include +#include +#include + +// OPENCASCADE includes +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +#ifndef WNT +#include +#endif + +#include +#include +#include +#include +#include +#include + +#define TYPE_EXP 0 +#define TYPE_LIN 1 + +class BLSURFPlugin_Attractor { + + public: + + BLSURFPlugin_Attractor (); + BLSURFPlugin_Attractor (const TopoDS_Face& Face, const TopoDS_Shape& Attractor, const std::string& attEntry, double Step = 0.015); + bool init(); // Calculates the discrete points correponding to attractor + // and intialises the map of distances + double GetSize (double u, double v); + TopoDS_Face GetFace() const { return _face; } + TopoDS_Shape GetAttractorShape() const { return _attractorShape; } + std::string GetAttractorEntry() const { return _attEntry; } + double GetStep() const { return _step; } + std::vector GetParameters() const + { + double tab_params[] = {_startSize, _endSize, _actionRadius, _constantRadius}; + std::vector params (tab_params, tab_params + sizeof(tab_params) / sizeof(double) ); + return params; + } + + void SetParameters(double Start_Size, double End_Size, double Action_Radius, double Constant_Radius); + void SetType(int type){ _type = type; } + + void BuildMap(); // Builds the map of distances between source point and any point P(u,v) + bool IsMapBuilt() const { return _isMapBuilt; } // Controls if the map has been built + bool Empty() const { return _isEmpty; } + + typedef std::vector TDiscreteParam; + typedef std::vector< std::vector > TDistMap; + typedef std::vector< std::vector > TPointSet; + typedef std::set< std::vector > TTrialSet; + typedef std::vector Trial_Pnt; + typedef std::vector IJ_Pnt; + + private: + + TopoDS_Face _face; + TopoDS_Shape _attractorShape; + std::string _attEntry; + TDiscreteParam _vectU; + TDiscreteParam _vectV; + TDistMap _DMap; + TPointSet _known; + TTrialSet _trial; + int _type; // Type of function used to calculate the size from the distance (unused for now) + double _step; // Step between to value of the discretized parametric space in U or V direction + int _gridU; // Number of grid points in U direction + int _gridV; // Number of grid points in V direction + double _u1, _u2, _v1, _v2; // Bounds of the parametric space of the face + double _startSize, _endSize; // User parameters + double _actionRadius, _constantRadius; // + + bool _isMapBuilt; + bool _isEmpty; + + double _distance(double u, double v); +}; + +#endif -- 2.39.2