src/HYDROData/HYDROData_DTM.cxx

   1
   2 #include <HYDROData_DTM.h>
   3 #include <HYDROData_Profile.h>
   4
   5 #include <Geom2d_BSplineCurve.hxx>
   6 #include <Geom2dAPI_Interpolate.hxx>
   7 #include <TColgp_HArray1OfPnt2d.hxx>
   8 #include <TColgp_Array1OfVec2d.hxx>
   9 #include <TColStd_HArray1OfBoolean.hxx>
  10 #include <TopoDS.hxx>
  11 #include <TopoDS_Edge.hxx>
  12 #include <TopoDS_Wire.hxx>
  13 #include <TopExp_Explorer.hxx>
  14 #include <BRep_Tool.hxx>
  15 #include <gp_Ax3.hxx>
  16 #include <Geom_Line.hxx>
  17 #include <Geom2d_Line.hxx>
  18 #include <Geom2d_TrimmedCurve.hxx>
  19 #include <Geom_BSplineCurve.hxx>
  20 #include <Geom2d_BSplineCurve.hxx>
  21 #include <TColStd_Array1OfReal.hxx>
  22 #include <TColStd_Array1OfInteger.hxx>
  23 #include <TColgp_Array1OfPnt.hxx>
  24 #include <Geom2dAPI_InterCurveCurve.hxx>
  25 #include <limits>
  26
  27 IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
  28 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
  29
  30
  31 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv )
  32   : myXcurv( theXCurv )
  33 {
  34 }
  35
  36 HYDROData_DTM::CurveUZ::~CurveUZ()
  37 {
  38 }
  39
  40 double HYDROData_DTM::CurveUZ::Xcurv() const
  41 {
  42   return myXcurv;
  43 }
  44
  45 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
  46 {
  47   HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv() );
  48   size_t n = size();
  49   res.reserve( n );
  50   for( int i=0; i<n; i++ )
  51   {
  52     PointUZ p;
  53     p.U = operator[]( i ).U + c[i].U;
  54     p.Z = operator[]( i ).Z;
  55     res.push_back( p );
  56   }
  57   return res;
  58 }
  59
  60 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
  61 {
  62   HYDROData_DTM::CurveUZ res( Xcurv()*d );
  63   size_t n = size();
  64   res.reserve( n );
  65   for( int i=0; i<n; i++ )
  66   {
  67     PointUZ p;
  68     p.U = operator[]( i ).U * d;
  69     p.Z = operator[]( i ).Z;
  70     res.push_back( p );
  71   }
  72   return res;
  73 }
  74
  75
  76
  77 HYDROData_DTM::HYDROData_DTM()
  78 {
  79 }
  80
  81 HYDROData_DTM::~HYDROData_DTM()
  82 {
  83 }
  84
  85
  86 void GetProperties( const Handle_HYDROData_Profile& theProfile,
  87                     gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
  88                     bool isNormalDir,
  89                     double& theZMin, double& theZMax )
  90 {
  91   theLowestPoint = theProfile->GetBottomPoint();
  92
  93   gp_XY aLeft, aRight;
  94   theProfile->GetLeftPoint( aLeft, true, true );
  95   theProfile->GetRightPoint( aRight, true, true );
  96   double x = aRight.X()-aLeft.X();
  97   double y = aRight.Y()-aLeft.Y();
  98   if( isNormalDir )
  99     theDir = gp_Vec2d( -y, x );
 100   else
 101     theDir = gp_Vec2d( x, y );
 102 }
 103
 104 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
 105                       double& theUMin, double& theUMax )
 106 {
 107   gp_Pnt p = thePnt.Transformed( theTr );
 108   double u = p.X();
 109   double z = p.Z();
 110   if( u<theUMin )
 111     theUMin = u;
 112   if( u>theUMax )
 113     theUMax = u;
 114   return gp_Pnt2d( u, z );
 115 }
 116
 117 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
 118                                     const gp_Trsf& theTr,
 119                                     double& theUMin, double& theUMax )
 120 {
 121   int low = thePoints.Lower(), up = thePoints.Upper();
 122   Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
 123   for( int i=low; i<=up; i++ )
 124     points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
 125   return points;
 126 }
 127
 128 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
 129                                 Standard_Real theFirst, Standard_Real theLast,
 130                                 const gp_Trsf& theTr,
 131                                 double& theUMin, double& theUMax )
 132 {
 133   if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
 134   {
 135     gp_Pnt aFirstPnt, aLastPnt;
 136     theCurve->D0( theFirst, aFirstPnt );
 137     theCurve->D0( theLast, aLastPnt );
 138
 139     gp_Pnt2d
 140       aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
 141       aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
 142
 143     gp_Vec2d dir( aFirst2d, aLast2d );
 144     Handle_Geom2d_Line aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
 145     return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
 146   }
 147
 148   if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
 149   {
 150     Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
 151
 152     Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
 153     const TColStd_Array1OfReal& knots = aSpline->Knots();
 154     const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
 155     int aDegree = aSpline->Degree();
 156
 157     return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
 158   }
 159
 160   return Handle(Geom2d_Curve)();
 161 }
 162
 163 Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
 164   const std::vector<Handle_HYDROData_Profile>& theProfiles,
 165   std::vector<double>& theDistances )
 166 {
 167   size_t n = theProfiles.size();
 168   Handle_Geom2d_BSplineCurve aResult;
 169
 170   Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
 171   TColgp_Array1OfVec2d tangents( 1, (int)n );
 172   Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
 173
 174   for( size_t i = 1; i <= n; i++ )
 175   {
 176     Handle_HYDROData_Profile aProfile = theProfiles[i-1];
 177     aProfile->Update();
 178
 179     gp_Pnt aLowest;
 180     gp_Vec2d aTangent;
 181     double zmin, zmax;
 182     GetProperties( aProfile, aLowest, aTangent, true, zmin, zmax );
 183     aTangent.Normalize();
 184
 185     points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
 186     tangents.SetValue( (int)i, aTangent );
 187     flags->SetValue( (int)i, Standard_True );
 188   }
 189
 190   Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
 191   anInterpolator.Load( tangents, flags );
 192   anInterpolator.Perform();
 193   if( anInterpolator.IsDone() )
 194   {
 195     aResult = anInterpolator.Curve();
 196     //TODO: fill the distances vector
 197   }
 198   return aResult;
 199 }
 200
 201 std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric(
 202   const Handle_HYDROData_Profile& theProfile,
 203   double& theUMin, double& theUMax )
 204 {
 205   std::vector<Handle_Geom2d_Curve> curves;
 206
 207   // Transformation of the coordinate systems
 208   gp_Pnt aLowest;
 209   gp_Vec2d aDir;
 210   double zmin, zmax;
 211   GetProperties( theProfile, aLowest, aDir, false, zmin, zmax );
 212
 213   gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
 214   gp_Ax3 aLocal( aLowest, gp_Dir( 0, 0, 1 ), gp_Dir( aDir.X(), aDir.Y(), 0 ) );
 215
 216   gp_Trsf aTransf;
 217   aTransf.SetTransformation( aStd3d, aLocal );
 218
 219   // Iteration via edges
 220   TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
 221   TopExp_Explorer anExp( aWire, TopAbs_EDGE );
 222   for( ; anExp.More(); anExp.Next() )
 223   {
 224     // Extract an edge from wire
 225     TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
 226
 227     // Extract a curve corresponding to the edge
 228     TopLoc_Location aLoc;
 229     Standard_Real aFirst, aLast;
 230     Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
 231
 232     // Convert the curve to 2d CS
 233     Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
 234     if( !aCurve2d.IsNull() )
 235       curves.push_back( aCurve2d );
 236   }
 237   return curves;
 238 }
 239
 240
 241 double CalcGC( const std::vector<gp_Pnt2d>& intersections )
 242 {
 243   double u = 0;
 244   size_t n = intersections.size();
 245   for( size_t i = 0; i < n; i++ )
 246   {
 247     u += intersections[i].X();
 248   }
 249   u /= n;
 250   return u;
 251 }
 252
 253 double CalcWidth( const std::vector<gp_Pnt2d>& intersections )
 254 {
 255   double umin = std::numeric_limits<double>::max(),
 256          umax = -umin;
 257
 258   size_t n = intersections.size();
 259   if( n <= 1 )
 260     return 0;
 261
 262   for( size_t i = 0; i < n; i++ )
 263   {
 264     double u = intersections[i].X();
 265     if( u<umin )
 266       umin = u;
 267     if( u>umax )
 268       umax = u;
 269   }
 270   return umax-umin;
 271 }
 272
 273 void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile,
 274                                            double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
 275                                            CurveUZ& theMidPointCurve,
 276                                            CurveUZ& theWidthCurve,
 277                                            double theTolerance )
 278 {
 279   double aDblMax = std::numeric_limits<double>::max(),
 280          aUMin = aDblMax,
 281          aUMax = -aUMin;
 282
 283   std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax );
 284   size_t n = curves.size();
 285
 286   if( n==0 )
 287     return;
 288
 289   // we add the "virtual" vertical lines to simulate the intersection with profile
 290   gp_Pnt2d aFirst, aLast;
 291   curves[0]->D0( curves[0]->FirstParameter(), aFirst );
 292   curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
 293   Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
 294   Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aDblMax );
 295
 296   Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
 297   Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aDblMax );
 298
 299   curves.push_back( aT1 );
 300   curves.push_back( aT2 );
 301
 302   int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
 303   theMidPointCurve = CurveUZ( theXCurv );
 304   theMidPointCurve.reserve( psize );
 305   theWidthCurve = CurveUZ( theXCurv );
 306   theWidthCurve.reserve( psize );
 307
 308   n = curves.size();
 309   // for each discrete value of z we search intersection with profile
 310   for( double z = theMinZ; z <= theMaxZ; z += theDDZ )
 311   {
 312     Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z ), gp_Dir2d( 1, 0 ) );
 313     std::vector<gp_Pnt2d> intersections;
 314     for( size_t i = 0; i < n; i++ )
 315     {
 316       Handle_Geom2d_Curve aCurve = curves[i];
 317       Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
 318       for( int k=0, m=anIntersect.NbPoints(); k<m; k++ )
 319         intersections.push_back( anIntersect.Point( k ) );
 320     }
 321
 322     if( intersections.size() >= 2 )
 323     {
 324       double u_mid = CalcGC( intersections );
 325       PointUZ p_mid;
 326       p_mid.U = u_mid;
 327       p_mid.Z = z;
 328       theMidPointCurve.push_back( p_mid );
 329
 330       double u_wid = CalcWidth( intersections );
 331       PointUZ p_wid;
 332       p_wid.U = u_wid;
 333       p_wid.Z = z;
 334       theWidthCurve.push_back( p_wid );
 335     }
 336   }
 337 }
 338
 339 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
 340                                  int theNbSteps, std::vector<CurveUZ>& theInterpolation,
 341                                  bool isAddSecond )
 342 {
 343   theInterpolation.clear();
 344   int d = isAddSecond ? 2 : 1;
 345   theInterpolation.reserve( theNbSteps+d );
 346   double dt = 1.0 / double( theNbSteps + 1 );
 347   double t = dt;
 348   theInterpolation.push_back( theCurveA );
 349   for( int i=0; i<theNbSteps; i++, t+=dt )
 350   {
 351     CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
 352     theInterpolation.push_back( anInterp );
 353   }
 354   if( isAddSecond )
 355     theInterpolation.push_back( theCurveB );
 356 }
 357
 358 void HYDROData_DTM::CurveTo3d( const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
 359                                AltitudePoints& thePoints )
 360 {
 361   //TODO
 362 }
 363
 364 inline double max( double a, double b )
 365 {
 366   if( a>b )
 367     return a;
 368   else
 369     return b;
 370 }
 371
 372 void HYDROData_DTM::Interpolate( const Handle_HYDROData_Profile& theProfileA,
 373                                  double theXCurvA,
 374                                  const Handle_HYDROData_Profile& theProfileB,
 375                                  double theXCurvB,
 376                                  double theDDZ, int theNbSteps, bool isAddSecond,
 377                                  AltitudePoints& thePoints )
 378 {
 379   double zminA, zmaxA, zminB, zmaxB;
 380   gp_Pnt lowestA, lowestB;
 381   gp_Vec2d dirA, dirB;
 382
 383   GetProperties( theProfileA, lowestA, dirA, false, zminA, zmaxA );
 384   GetProperties( theProfileB, lowestB, dirB, false, zminB, zmaxB );
 385
 386   double zmin = max( zminA, zminB );
 387   double zmax = max( zmaxA, zmaxB );
 388
 389   CurveUZ midA(0), midB(0);
 390   CurveUZ widA(0), widB(0);
 391
 392   ProfileDiscretization( theProfileA, theXCurvA, zmin, zmax, theDDZ, midA, widA );
 393   ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB );
 394
 395   std::vector<CurveUZ> mid, wid;
 396   Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
 397   Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
 398
 399   size_t q = mid.size();
 400   for( size_t i=0; i<q; i++ )
 401     CurveTo3d( mid[i], wid[i], thePoints );
 402 }
 403
 404 void HYDROData_DTM::Interpolate( const std::vector<Handle_HYDROData_Profile>& theProfiles,
 405                                  double theDDZ, double theSpatialStep, AltitudePoints& thePoints )
 406 {
 407   size_t n = theProfiles.size();
 408   if( n<=1 )
 409     return;
 410
 411   std::vector<double> distances;
 412   Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
 413   if( aHydraulicAxis.IsNull() )
 414     return;
 415
 416   for( size_t i=0, n1=n-1; i<n1; i++ )
 417   {
 418     double aDistance = distances[i+1]-distances[i];
 419     int aNbSteps = int(aDistance/theSpatialStep) + 1;
 420     bool isAddSecond = i==n1-1;
 421     Interpolate( theProfiles[i], distances[i],
 422       theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, thePoints );
 423   }
 424 }