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 <Geom2dAPI_ProjectPointOnCurve.hxx>
  26 #include <Geom2dAdaptor_Curve.hxx>
  27 #include <GCPnts_AbscissaPoint.hxx>
  28 #include <limits>
  29
  30 IMPLEMENT_STANDARD_HANDLE( HYDROData_DTM, HYDROData_Bathymetry )
  31 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
  32
  33
  34 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv )
  35   : myXcurv( theXCurv )
  36 {
  37 }
  38
  39 HYDROData_DTM::CurveUZ::~CurveUZ()
  40 {
  41 }
  42
  43 double HYDROData_DTM::CurveUZ::Xcurv() const
  44 {
  45   return myXcurv;
  46 }
  47
  48 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
  49 {
  50   HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv() );
  51   size_t n = size();
  52   res.reserve( n );
  53   for( int i=0; i<n; i++ )
  54   {
  55     PointUZ p;
  56     p.U = operator[]( i ).U + c[i].U;
  57     p.Z = operator[]( i ).Z;
  58     res.push_back( p );
  59   }
  60   return res;
  61 }
  62
  63 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
  64 {
  65   HYDROData_DTM::CurveUZ res( Xcurv()*d );
  66   size_t n = size();
  67   res.reserve( n );
  68   for( int i=0; i<n; i++ )
  69   {
  70     PointUZ p;
  71     p.U = operator[]( i ).U * d;
  72     p.Z = operator[]( i ).Z;
  73     res.push_back( p );
  74   }
  75   return res;
  76 }
  77
  78
  79
  80 HYDROData_DTM::HYDROData_DTM()
  81 {
  82 }
  83
  84 HYDROData_DTM::~HYDROData_DTM()
  85 {
  86 }
  87
  88 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
  89 {
  90   return GetReferenceObjects( DataTag_Profiles );
  91 }
  92
  93 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
  94 {
  95   SetReferenceObjects( theProfiles, DataTag_Profiles );
  96 }
  97
  98 double HYDROData_DTM::GetDDZ() const
  99 {
 100   return GetDouble( DataTag_DDZ );
 101 }
 102
 103 void HYDROData_DTM::SetDDZ( double theDDZ )
 104 {
 105   SetDouble( DataTag_DDZ, theDDZ );
 106 }
 107
 108 double HYDROData_DTM::GetSpatialStep() const
 109 {
 110   return GetDouble( DataTag_SpatialStep );
 111 }
 112
 113 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
 114 {
 115   SetDouble( DataTag_SpatialStep, theSpatialStep );
 116 }
 117
 118 void HYDROData_DTM::Update()
 119 {
 120   HYDROData_SequenceOfObjects objs = GetProfiles();
 121   int aLower = objs.Lower(), anUpper = objs.Upper();
 122   size_t n = anUpper-aLower+1;
 123
 124   std::vector<Handle_HYDROData_Profile> profiles;
 125   profiles.reserve( n );
 126   for( int i=aLower; i<=anUpper; i++ )
 127   {
 128     Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( objs.Value( i ) );
 129     if( !aProfile.IsNull() )
 130       profiles.push_back( aProfile );
 131   }
 132
 133   double ddz = GetDDZ();
 134   double step = GetSpatialStep();
 135   const double EPS = 1E-3;
 136   AltitudePoints points;
 137   if( ddz>EPS && step>EPS )
 138     points = Interpolate( profiles, ddz, step );
 139   SetAltitudePoints( points );
 140 }
 141
 142
 143
 144
 145
 146
 147
 148
 149
 150
 151 void HYDROData_DTM::GetProperties( const Handle_HYDROData_Profile& theProfile,
 152                     gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
 153                     bool isNormalDir,
 154                     double& theZMin, double& theZMax )
 155 {
 156   theLowestPoint = theProfile->GetBottomPoint();
 157
 158   gp_XY aLeft, aRight;
 159   theProfile->GetLeftPoint( aLeft, true, true );
 160   theProfile->GetRightPoint( aRight, true, true );
 161   double x = aRight.X()-aLeft.X();
 162   double y = aRight.Y()-aLeft.Y();
 163   if( isNormalDir )
 164     theDir = gp_Vec2d( -y, x );
 165   else
 166     theDir = gp_Vec2d( x, y );
 167
 168   HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
 169   int lo = points.Lower();
 170   int up = points.Upper();
 171   theZMin = std::numeric_limits<double>::max();
 172   theZMax = -theZMin;
 173   for( int i=lo; i<=up; i++ )
 174   {
 175     double z = points.Value( i ).Z();
 176     if( z>theZMax )
 177       theZMax = z;
 178     if( z<theZMin )
 179       theZMin = z;
 180   }
 181 }
 182
 183 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
 184                       double& theUMin, double& theUMax )
 185 {
 186   gp_Pnt p = thePnt.Transformed( theTr );
 187   double u = p.X();
 188   double z = p.Z();
 189   if( u<theUMin )
 190     theUMin = u;
 191   if( u>theUMax )
 192     theUMax = u;
 193   return gp_Pnt2d( u, z );
 194 }
 195
 196 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
 197                                     const gp_Trsf& theTr,
 198                                     double& theUMin, double& theUMax )
 199 {
 200   int low = thePoints.Lower(), up = thePoints.Upper();
 201   Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
 202   for( int i=low; i<=up; i++ )
 203     points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
 204   return points;
 205 }
 206
 207 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
 208                                 Standard_Real theFirst, Standard_Real theLast,
 209                                 const gp_Trsf& theTr,
 210                                 double& theUMin, double& theUMax )
 211 {
 212   if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
 213   {
 214     gp_Pnt aFirstPnt, aLastPnt;
 215     theCurve->D0( theFirst, aFirstPnt );
 216     theCurve->D0( theLast, aLastPnt );
 217
 218     gp_Pnt2d
 219       aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
 220       aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
 221
 222     gp_Vec2d dir( aFirst2d, aLast2d );
 223     Handle_Geom2d_Line aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
 224     return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
 225   }
 226
 227   if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
 228   {
 229     Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
 230
 231     Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
 232     const TColStd_Array1OfReal& knots = aSpline->Knots();
 233     const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
 234     int aDegree = aSpline->Degree();
 235
 236     return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
 237   }
 238
 239   return Handle(Geom2d_Curve)();
 240 }
 241
 242 Handle_Geom2d_BSplineCurve HYDROData_DTM::CreateHydraulicAxis(
 243   const std::vector<Handle_HYDROData_Profile>& theProfiles,
 244   std::vector<double>& theDistances )
 245 {
 246   size_t n = theProfiles.size();
 247   Handle_Geom2d_BSplineCurve aResult;
 248
 249   Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
 250   TColgp_Array1OfVec2d tangents( 1, (int)n );
 251   Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
 252
 253   for( size_t i = 1; i <= n; i++ )
 254   {
 255     Handle_HYDROData_Profile aProfile = theProfiles[i-1];
 256     aProfile->Update();
 257
 258     gp_Pnt aLowest;
 259     gp_Vec2d aTangent;
 260     double zmin, zmax;
 261     GetProperties( aProfile, aLowest, aTangent, true, zmin, zmax );
 262     aTangent.Normalize();
 263
 264     points->SetValue( (int)i, gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
 265     tangents.SetValue( (int)i, aTangent );
 266     flags->SetValue( (int)i, Standard_True );
 267   }
 268
 269   Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
 270   anInterpolator.Load( tangents, flags );
 271   anInterpolator.Perform();
 272   if( anInterpolator.IsDone() )
 273   {
 274     aResult = anInterpolator.Curve();
 275
 276     //fill the distances vector
 277     Geom2dAdaptor_Curve anAdaptor( aResult );
 278
 279     theDistances.clear();
 280     theDistances.reserve( n );
 281     Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
 282     for( size_t i = 1; i <= n; i++ )
 283     {
 284       gp_Pnt2d aPnt = points->Value( i );
 285       Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
 286       Standard_Real aParam = aProject.LowerDistanceParameter();
 287       double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );
 288       theDistances.push_back( aDistance );
 289     }
 290   }
 291   return aResult;
 292 }
 293
 294 std::vector<Handle_Geom2d_Curve> HYDROData_DTM::ProfileToParametric(
 295   const Handle_HYDROData_Profile& theProfile,
 296   double& theUMin, double& theUMax )
 297 {
 298   std::vector<Handle_Geom2d_Curve> curves;
 299
 300   // Transformation of the coordinate systems
 301   gp_Pnt aLowest;
 302   gp_Vec2d aDir;
 303   double zmin, zmax;
 304   GetProperties( theProfile, aLowest, aDir, false, zmin, zmax );
 305
 306   gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
 307   gp_Ax3 aLocal( aLowest, gp_Dir( 0, 0, 1 ), gp_Dir( aDir.X(), aDir.Y(), 0 ) );
 308
 309   gp_Trsf aTransf;
 310   aTransf.SetTransformation( aStd3d, aLocal );
 311
 312   // Iteration via edges
 313   TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
 314   TopExp_Explorer anExp( aWire, TopAbs_EDGE );
 315   for( ; anExp.More(); anExp.Next() )
 316   {
 317     // Extract an edge from wire
 318     TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
 319
 320     // Extract a curve corresponding to the edge
 321     TopLoc_Location aLoc;
 322     Standard_Real aFirst, aLast;
 323     Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
 324
 325     // Convert the curve to 2d CS
 326     Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
 327     if( !aCurve2d.IsNull() )
 328       curves.push_back( aCurve2d );
 329   }
 330   return curves;
 331 }
 332
 333
 334 double CalcGC( const std::vector<gp_Pnt2d>& intersections )
 335 {
 336   double u = 0;
 337   size_t n = intersections.size();
 338   for( size_t i = 0; i < n; i++ )
 339   {
 340     u += intersections[i].X();
 341   }
 342   u /= n;
 343   return u;
 344 }
 345
 346 double CalcWidth( const std::vector<gp_Pnt2d>& intersections )
 347 {
 348   double umin = std::numeric_limits<double>::max(),
 349          umax = -umin;
 350
 351   size_t n = intersections.size();
 352   if( n <= 1 )
 353     return 0;
 354
 355   for( size_t i = 0; i < n; i++ )
 356   {
 357     double u = intersections[i].X();
 358     if( u<umin )
 359       umin = u;
 360     if( u>umax )
 361       umax = u;
 362   }
 363   return umax-umin;
 364 }
 365
 366 void HYDROData_DTM::ProfileDiscretization( const Handle_HYDROData_Profile& theProfile,
 367                                            double theXCurv, double theMinZ, double theMaxZ, double theDDZ,
 368                                            CurveUZ& theMidPointCurve,
 369                                            CurveUZ& theWidthCurve,
 370                                            double theTolerance )
 371 {
 372   double aDblMax = std::numeric_limits<double>::max(),
 373          aUMin = aDblMax,
 374          aUMax = -aUMin,
 375          aVMax = 1000000;
 376
 377   std::vector<Handle_Geom2d_Curve> curves = ProfileToParametric( theProfile, aUMin, aUMax );
 378   size_t n = curves.size();
 379
 380   if( n==0 )
 381     return;
 382
 383   // we add the "virtual" vertical lines to simulate the intersection with profile
 384   gp_Pnt2d aFirst, aLast;
 385   curves[0]->D0( curves[0]->FirstParameter(), aFirst );
 386   curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
 387   Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
 388   Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
 389
 390   Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
 391   Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
 392
 393   curves.push_back( aT1 );
 394   curves.push_back( aT2 );
 395
 396   int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
 397   theMidPointCurve = CurveUZ( theXCurv );
 398   theMidPointCurve.reserve( psize );
 399   theWidthCurve = CurveUZ( theXCurv );
 400   theWidthCurve.reserve( psize );
 401
 402   n = curves.size();
 403   // for each discrete value of z we search intersection with profile
 404   for( double z = theMinZ; z <= theMaxZ; z += theDDZ )
 405   {
 406     Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z ), gp_Dir2d( 1, 0 ) );
 407     std::vector<gp_Pnt2d> intersections;
 408     for( size_t i = 0; i < n; i++ )
 409     {
 410       Handle_Geom2d_Curve aCurve = curves[i];
 411       Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
 412       for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
 413         intersections.push_back( anIntersect.Point( k ) );
 414     }
 415
 416     if( intersections.size() >= 2 )
 417     {
 418       double u_mid = CalcGC( intersections );
 419       PointUZ p_mid;
 420       p_mid.U = u_mid;
 421       p_mid.Z = z;
 422       theMidPointCurve.push_back( p_mid );
 423
 424       double u_wid = CalcWidth( intersections );
 425       PointUZ p_wid;
 426       p_wid.U = u_wid;
 427       p_wid.Z = z;
 428       theWidthCurve.push_back( p_wid );
 429     }
 430   }
 431 }
 432
 433 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
 434                                  int theNbSteps, std::vector<CurveUZ>& theInterpolation,
 435                                  bool isAddSecond )
 436 {
 437   theInterpolation.clear();
 438   int d = isAddSecond ? 2 : 1;
 439   theInterpolation.reserve( theNbSteps+d );
 440   double dt = 1.0 / double( theNbSteps + 1 );
 441   double t = dt;
 442   theInterpolation.push_back( theCurveA );
 443   for( int i=0; i<theNbSteps; i++, t+=dt )
 444   {
 445     CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
 446     theInterpolation.push_back( anInterp );
 447   }
 448   if( isAddSecond )
 449     theInterpolation.push_back( theCurveB );
 450 }
 451
 452 void HYDROData_DTM::CurveTo3D( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
 453                                const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
 454                                AltitudePoints& thePoints )
 455 {
 456   Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
 457   GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
 458   double aParam = ap.Parameter();
 459
 460   gp_Pnt2d point;
 461   gp_Vec2d tangent, profile_dir;
 462   anAdaptor.D1( aParam, point, tangent );
 463   profile_dir.SetCoord( tangent.Y(), -tangent.X() );
 464   profile_dir.Normalize();
 465
 466   size_t n = theMidCurve.size();
 467   for( size_t i=0; i<n; i++ )
 468   {
 469     double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
 470     double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
 471
 472     gp_Pnt2d p1 = point.Translated( param1 * profile_dir / 2 );
 473     gp_Pnt2d p2 = point.Translated( param2 * profile_dir / 2 );
 474
 475     double z = theMidCurve[i].Z;
 476
 477     AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
 478     thePoints.push_back( p3d_1 );
 479     thePoints.push_back( p3d_2 );
 480   }
 481 }
 482
 483 inline double max( double a, double b )
 484 {
 485   if( a>b )
 486     return a;
 487   else
 488     return b;
 489 }
 490
 491 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
 492   ( const Handle_Geom2d_BSplineCurve& theHydraulicAxis,
 493     const Handle_HYDROData_Profile& theProfileA,
 494     double theXCurvA,
 495     const Handle_HYDROData_Profile& theProfileB,
 496     double theXCurvB,
 497     double theDDZ, int theNbSteps, bool isAddSecond )
 498 {
 499   double zminA, zmaxA, zminB, zmaxB;
 500   gp_Pnt lowestA, lowestB;
 501   gp_Vec2d dirA, dirB;
 502
 503   GetProperties( theProfileA, lowestA, dirA, false, zminA, zmaxA );
 504   GetProperties( theProfileB, lowestB, dirB, false, zminB, zmaxB );
 505
 506   double zmin = max( zminA, zminB );
 507   double zmax = max( zmaxA, zmaxB );
 508
 509   CurveUZ midA(0), midB(0);
 510   CurveUZ widA(0), widB(0);
 511
 512   ProfileDiscretization( theProfileA, theXCurvA, zmin, zmax, theDDZ, midA, widA );
 513   ProfileDiscretization( theProfileB, theXCurvB, zmin, zmax, theDDZ, midB, widB );
 514
 515   std::vector<CurveUZ> mid, wid;
 516   Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
 517   Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
 518
 519   size_t p = mid.size();
 520   size_t q = p>0 ? 2*mid[0].size() : 1;
 521   AltitudePoints points;
 522   points.reserve( p*q );
 523   for( size_t i=0; i<p; i++ )
 524     CurveTo3D( theHydraulicAxis, mid[i], wid[i], points );
 525   return points;
 526 }
 527
 528 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
 529   ( const std::vector<Handle_HYDROData_Profile>& theProfiles,
 530     double theDDZ, double theSpatialStep )
 531 {
 532   AltitudePoints points;
 533   size_t n = theProfiles.size();
 534   if( n<=1 )
 535     return points;
 536
 537   std::vector<double> distances;
 538   Handle_Geom2d_BSplineCurve aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
 539   if( aHydraulicAxis.IsNull() )
 540     return points;
 541
 542   for( size_t i=0, n1=n-1; i<n1; i++ )
 543   {
 544     double aDistance = distances[i+1]-distances[i];
 545     int aNbSteps = int(aDistance/theSpatialStep) + 1;
 546     bool isAddSecond = i==n1-1;
 547     AltitudePoints local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
 548       theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond );
 549     if( i==0 )
 550       points.reserve( local_points.size() * ( n-1 ) );
 551     for( size_t j=0, m=local_points.size(); j<m; j++ )
 552       points.push_back( local_points[j] );
 553   }
 554   return points;
 555 }