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 <GeomAPI_Interpolate.hxx>
  22 #include <TColStd_Array1OfReal.hxx>
  23 #include <TColStd_Array1OfInteger.hxx>
  24 #include <TColgp_Array1OfPnt.hxx>
  25 #include <TColgp_Array1OfVec.hxx>
  26 #include <TColgp_HArray1OfPnt.hxx>
  27 #include <Geom2dAPI_InterCurveCurve.hxx>
  28 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
  29 #include <Geom2dAdaptor_Curve.hxx>
  30 #include <GCPnts_AbscissaPoint.hxx>
  31 #include <BRepBuilderAPI_MakeEdge.hxx>
  32 #include <limits>
  33
  34 #include <BRepLib_MakeEdge.hxx>
  35 #include <BRepLib_MakeWire.hxx>
  36 #include <BRep_Builder.hxx>
  37 #include <GeomProjLib.hxx>
  38 #include <Geom_TrimmedCurve.hxx>
  39 #include <BRepTools_WireExplorer.hxx>
  40 #include <TopTools_IndexedMapOfShape.hxx>
  41 #include <BRepBuilderAPI_MakeFace.hxx>
  42 #include <TopExp.hxx>
  43 #include <TopTools_IndexedMapOfOrientedShape.hxx>
  44 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  45
  46 #include <BRepLib_MakeEdge.hxx>
  47 #include <BRepLib_MakeWire.hxx>
  48 #include <BRep_Builder.hxx>
  49 #include <ShapeAnalysis_Wire.hxx>
  50 #include <BRepAlgo_NormalProjection.hxx>
  51 #include <ShapeUpgrade_UnifySameDomain.hxx>
  52 #include <BRepBuilderAPI_MakePolygon.hxx>
  53 #include <BOPAlgo_Builder.hxx>
  54 #include <BRepAdaptor_Curve.hxx>
  55 #include <GeomProjLib.hxx>
  56 #include <gp_Pln.hxx>
  57 #include <TopTools_IndexedDataMapOfShapeShape.hxx>
  58 #include <TopTools_ListIteratorOfListOfShape.hxx>
  59 #include <TopTools_SequenceOfShape.hxx>
  60 #include <assert.h>
  61 #include <float.h>
  62
  63 IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
  64
  65 HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ, double theMaxZ )
  66   : myXcurv( theXCurv ), myProfileDir( theProfileDir ), myDeltaZ( theDeltaZ ), myMaxZ (theMaxZ)
  67 {
  68 }
  69
  70 HYDROData_DTM::CurveUZ::~CurveUZ()
  71 {
  72 }
  73
  74 double HYDROData_DTM::CurveUZ::Xcurv() const
  75 {
  76   return myXcurv;
  77 }
  78
  79 gp_Vec2d HYDROData_DTM::CurveUZ::ProfileDir() const
  80 {
  81   return myProfileDir;
  82 }
  83
  84 double HYDROData_DTM::CurveUZ::DeltaZ() const
  85 {
  86   return myDeltaZ;
  87 }
  88
  89 double HYDROData_DTM::CurveUZ::MaxZ() const
  90 {
  91   return myMaxZ;
  92 }
  93
  94 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator + ( const CurveUZ& c ) const
  95 {
  96   HYDROData_DTM::CurveUZ res( Xcurv() + c.Xcurv(), ProfileDir() + c.ProfileDir(), DeltaZ() + c.DeltaZ(), MaxZ() + c.MaxZ() );
  97   size_t n = size(), n1 = c.size();
  98   if( n!=n1 )
  99   {
 100     std::cout << "Warning: different number of points in curves: " << n << ", " << n1 << std::endl;
 101   }
 102   int q = min(n, n1);
 103   res.reserve( q );
 104   for( int i=0; i<q; i++ )
 105   {
 106     PointUZ p;
 107     p.U = operator[]( i ).U + c[i].U;
 108     p.Z = operator[]( i ).Z;
 109     res.push_back( p );
 110   }
 111   return res;
 112 }
 113
 114 HYDROData_DTM::CurveUZ HYDROData_DTM::CurveUZ::operator * ( double d ) const
 115 {
 116   HYDROData_DTM::CurveUZ res( Xcurv()*d, ProfileDir()*d, DeltaZ()*d, MaxZ()*d );
 117   size_t n = size();
 118   res.reserve( n );
 119   for( int i=0; i<n; i++ )
 120   {
 121     PointUZ p;
 122     p.U = operator[]( i ).U * d;
 123     p.Z = operator[]( i ).Z;
 124     res.push_back( p );
 125   }
 126   return res;
 127 }
 128
 129
 130
 131
 132 HYDROData_DTM::HYDROData_DTM()
 133 {
 134 }
 135
 136 HYDROData_DTM::~HYDROData_DTM()
 137 {
 138 }
 139
 140 HYDROData_SequenceOfObjects HYDROData_DTM::GetProfiles() const
 141 {
 142   return GetReferenceObjects( DataTag_Profiles );
 143 }
 144
 145 void HYDROData_DTM::SetProfiles( const HYDROData_SequenceOfObjects& theProfiles )
 146 {
 147   SetReferenceObjects( theProfiles, DataTag_Profiles );
 148   Changed( Geom_3d );
 149 }
 150
 151 double HYDROData_DTM::GetDDZ() const
 152 {
 153   return GetDouble( DataTag_DDZ );
 154 }
 155
 156 void HYDROData_DTM::SetDDZ( double theDDZ )
 157 {
 158   SetDouble( DataTag_DDZ, theDDZ );
 159   Changed( Geom_3d );
 160 }
 161
 162 double HYDROData_DTM::GetSpatialStep() const
 163 {
 164   return GetDouble( DataTag_SpatialStep );
 165 }
 166
 167 void HYDROData_DTM::SetSpatialStep( double theSpatialStep )
 168 {
 169   SetDouble( DataTag_SpatialStep, theSpatialStep );
 170   Changed( Geom_3d );
 171 }
 172
 173 void HYDROData_DTM::PointsToWire(const AltitudePoints& pnts, TopoDS_Wire& W )
 174 {
 175
 176   BRepBuilderAPI_MakePolygon PM;
 177   for (int i = 0; i < pnts.size(); i++)
 178     PM.Add(gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
 179
 180   W = PM.Wire();
 181 }
 182
 183 void HYDROData_DTM::PointsToEdge(const AltitudePoints& pnts, TopoDS_Edge& E )
 184 {
 185   Handle(TColgp_HArray1OfPnt) gpPoints = new TColgp_HArray1OfPnt(1, (int)pnts.size());
 186
 187   for (int i = 0; i < pnts.size(); i++)
 188     gpPoints->SetValue(i+1, gp_Pnt(pnts[i].X, pnts[i].Y, pnts[i].Z));
 189
 190   GeomAPI_Interpolate anInterpolator(gpPoints, Standard_False,1.0e-6);
 191   anInterpolator.Perform() ;
 192   if (anInterpolator.IsDone())
 193   {
 194     Handle(Geom_BSplineCurve) C = anInterpolator.Curve();
 195     E = BRepBuilderAPI_MakeEdge(C).Edge();
 196   }
 197 }
 198
 199 TopTools_IndexedMapOfOrientedShape HYDROData_DTM::Create3DShape(const AltitudePoints& left,
 200                                                                 const AltitudePoints& right,
 201                                                                 const std::vector<AltitudePoints>& main_profiles)
 202 {
 203   TopTools_IndexedMapOfOrientedShape ll;
 204   //TopoDS_Wire LWire, RWire;
 205   //PointsToWire(left, LWire);
 206   //PointsToWire(right, RWire);
 207   TopoDS_Edge LEdge, REdge;
 208   PointsToEdge(left, LEdge);
 209   PointsToEdge(right, REdge);
 210   if (!LEdge.IsNull())
 211     ll.Add(LEdge.Oriented(TopAbs_FORWARD));
 212
 213   for (int k = 0; k < main_profiles.size(); k++)
 214   {
 215     TopoDS_Wire W;
 216     PointsToWire(main_profiles[k], W);
 217     TopAbs_Orientation Ori = TopAbs_INTERNAL;
 218     if (k == 0 || k == main_profiles.size() - 1)
 219       Ori = TopAbs_FORWARD;
 220     ll.Add(W.Oriented(Ori));
 221   }
 222
 223   if (!REdge.IsNull())
 224     ll.Add(REdge.Oriented(TopAbs_FORWARD));
 225   //yes, add subshapes in this order (left + profiles + right)
 226   //otherwise the projected wire will be non-manifold
 227
 228   return ll;
 229 }
 230
 231
 232 void HYDROData_DTM::Update()
 233 {
 234   AltitudePoints points;
 235   TopoDS_Shape Out3dPres;
 236   TopoDS_Shape Out2dPres;
 237   TopoDS_Shape OutLeftB;
 238   TopoDS_Shape OutRightB;
 239   TopoDS_Shape OutInlet;
 240   TopoDS_Shape OutOutlet;
 241
 242   HYDROData_SequenceOfObjects objs = GetProfiles();
 243   double ddz = GetDDZ();
 244   double step = GetSpatialStep();
 245   std::set<int> InvInd;
 246   bool WireIntersections; //__TODO
 247   CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, InvInd, -1, WireIntersections );
 248   SetAltitudePoints( points );
 249
 250   SetShape( DataTag_LeftBankShape, OutLeftB);
 251   SetShape( DataTag_RightBankShape, OutRightB);
 252   SetShape( DataTag_InletShape, OutInlet);
 253   SetShape( DataTag_OutletShape, OutOutlet );
 254   SetShape( DataTag_3DShape, Out3dPres );
 255   SetShape( DataTag_2DShape, Out2dPres );
 256
 257   HYDROData_Bathymetry::Update();
 258 }
 259
 260 void HYDROData_DTM::GetPresentationShapes( TopoDS_Shape& Out3dPres,
 261                                            TopoDS_Shape& Out2dPres,
 262                                            TopoDS_Shape& OutLeftB,
 263                                            TopoDS_Shape& OutRightB,
 264                                            TopoDS_Shape& OutInlet,
 265                                            TopoDS_Shape& OutOutlet )
 266 {
 267   //without update!
 268   OutLeftB = GetShape( DataTag_LeftBankShape);
 269   OutRightB = GetShape( DataTag_RightBankShape);
 270   OutInlet = GetShape( DataTag_InletShape);
 271   OutOutlet = GetShape( DataTag_OutletShape );
 272   Out3dPres = GetShape( DataTag_3DShape );
 273   Out2dPres = GetShape( DataTag_2DShape );
 274 }
 275 void HYDROData_DTM::CreateProfilesFromDTM (const HYDROData_SequenceOfObjects& InpProfiles,
 276                                            double ddz,
 277                                            double step,
 278                                            AltitudePoints& points,
 279                                            TopoDS_Shape& Out3dPres,
 280                                            TopoDS_Shape& Out2dPres,
 281                                            TopoDS_Shape& OutLeftB,
 282                                            TopoDS_Shape& OutRightB,
 283                                            TopoDS_Shape& OutInlet,
 284                                            TopoDS_Shape& OutOutlet,
 285                                            bool Create3dPres,
 286                                            bool Create2dPres,
 287                                            std::set<int>& InvInd,
 288                                            int thePntsLimit,
 289                                            bool& WireIntersections)
 290 {
 291   int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper();
 292   size_t n = anUpper - aLower + 1;
 293
 294   std::vector<Handle(HYDROData_Profile)> profiles;
 295   profiles.reserve( n );
 296   for( int i=aLower; i<=anUpper; i++ )
 297   {
 298     Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast( InpProfiles.Value( i ) );
 299     if( !aProfile.IsNull() )
 300       profiles.push_back( aProfile );
 301   }
 302   const double EPS = 1E-3;
 303   AltitudePoints left;
 304   AltitudePoints right;
 305   std::vector<AltitudePoints> main_profiles;
 306
 307   if( thePntsLimit > 0 )
 308   {
 309     int aNbPoints = EstimateNbPoints( profiles, ddz, step );
 310     if( aNbPoints < 0 || aNbPoints > thePntsLimit )
 311       return;
 312   }
 313
 314   if( ddz>EPS && step>EPS )
 315     CreateProfiles(profiles, ddz, step, left, right, points, main_profiles,
 316     Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections );
 317 }
 318
 319 bool HYDROData_DTM::GetPlanarFaceFromBanks( const TopoDS_Edge& LB, const TopoDS_Edge& RB, TopoDS_Face& outF,
 320   TopTools_SequenceOfShape* Boundr)
 321 {
 322   BRep_Builder BB;
 323   TopoDS_Face F;
 324   Handle_Geom_Plane refpl = new Geom_Plane(gp_Pnt(0,0,0), gp_Dir(0,0,1));
 325
 326   TopoDS_Vertex VFI, VLI, VFO, VLO;
 327   TopoDS_Edge prLB;
 328   TopoDS_Edge prRB;
 329
 330
 331   BRepAdaptor_Curve LBAD(LB);
 332   Handle_Geom_Curve LBPC = GeomProjLib::ProjectOnPlane(LBAD.Curve().Curve(), refpl, gp_Dir(0, 0, -1), 1 );
 333   prLB = BRepLib_MakeEdge(LBPC).Edge();
 334
 335   BRepAdaptor_Curve RBAD(RB);
 336   Handle_Geom_Curve RBPC = GeomProjLib::ProjectOnPlane(RBAD.Curve().Curve(), refpl, gp_Dir(0, 0, -1), 1 );
 337   prRB = BRepLib_MakeEdge(RBPC).Edge();
 338
 339   TopExp::Vertices(prLB, VFI, VFO, 1);
 340   TopExp::Vertices(prRB, VLI, VLO, 1);
 341   TopoDS_Edge prIL = BRepLib_MakeEdge(VFI, VLI).Edge();
 342   TopoDS_Edge prOL = BRepLib_MakeEdge(VFO, VLO).Edge();
 343   TopoDS_Wire prW = BRepLib_MakeWire(prLB, prIL, prOL, prRB).Wire();
 344   outF = BRepBuilderAPI_MakeFace(refpl->Pln(), prW, 1).Face();
 345
 346   if (Boundr)
 347   {
 348     Boundr->Append(prLB);
 349     Boundr->Append(prIL);
 350     Boundr->Append(prOL);
 351     Boundr->Append(prRB);
 352   }
 353
 354   ShapeAnalysis_Wire WA(prW, outF, Precision::Confusion());
 355   bool res = WA.CheckSelfIntersection();
 356   return !res;
 357 }
 358
 359 void HYDROData_DTM::CreateProfiles(const std::vector<Handle(HYDROData_Profile)>& theProfiles,
 360                                    double theDDZ,
 361                                    double theSpatialStep,
 362                                    AltitudePoints& theOutLeft,
 363                                    AltitudePoints& theOutRight,
 364                                    AltitudePoints& theOutPoints,
 365                                    std::vector<AltitudePoints>& theOutMainProfiles,
 366                                    TopoDS_Shape& Out3dPres,
 367                                    TopoDS_Shape& Out2dPres,
 368                                    TopoDS_Shape& OutLeftB,
 369                                    TopoDS_Shape& OutRightB,
 370                                    TopoDS_Shape& OutInlet,
 371                                    TopoDS_Shape& OutOutlet,
 372                                    bool Create3dPres,
 373                                    bool Create2dPres,
 374                                    std::set<int>& InvInd,
 375                                    bool& ProjStat)
 376 {
 377   if (theProfiles.empty())
 378     return;
 379   theOutPoints = Interpolate( theProfiles, theDDZ, theSpatialStep, theOutLeft, theOutRight, theOutMainProfiles, InvInd );
 380   //note that if Create3dPres is false => Create2dPres flag is meaningless!
 381   if( theOutPoints.empty() )
 382     return;
 383
 384   if (Create3dPres)
 385   {
 386     TopTools_IndexedMapOfOrientedShape ll = Create3DShape( theOutLeft, theOutRight, theOutMainProfiles);
 387
 388     if (ll.IsEmpty())
 389       return;
 390     BRep_Builder BB;
 391     TopoDS_Compound cmp;
 392     BB.MakeCompound(cmp);
 393     for (int i = 1; i <= ll.Extent(); i++)
 394       BB.Add(cmp, ll(i));
 395
 396     Out3dPres = cmp;
 397
 398     //same order as in HYDROData_DTM::Update()
 399     OutLeftB = ll(1);
 400     OutRightB = ll(ll.Extent());
 401     OutInlet = ll(2);
 402     OutOutlet = ll(ll.Extent() - 1);
 403
 404     if (Create2dPres)
 405     {
 406       TopoDS_Face outF;
 407       ProjStat = GetPlanarFaceFromBanks(TopoDS::Edge(OutLeftB), TopoDS::Edge(OutRightB), outF, NULL);
 408       Out2dPres = outF;
 409     };
 410   }
 411 }
 412
 413
 414
 415
 416 void HYDROData_DTM::GetProperties( const Handle(HYDROData_Profile)& theProfile,
 417                     gp_Pnt& theLowestPoint, gp_Vec2d& theDir,
 418                     double& theZMin, double& theZMax )
 419 {
 420   theLowestPoint = theProfile->GetBottomPoint();
 421
 422   gp_XY aLeft, aRight;
 423   theProfile->GetLeftPoint( aLeft, true, true );
 424   theProfile->GetRightPoint( aRight, true, true );
 425   double x = aRight.X()-aLeft.X();
 426   double y = aRight.Y()-aLeft.Y();
 427   theDir = gp_Vec2d( x, y );
 428
 429   HYDROData_Profile::ProfilePoints points = theProfile->GetProfilePoints();
 430   int lo = points.Lower();
 431   int up = points.Upper();
 432   theZMin = DBL_MAX;
 433   theZMax = -theZMin;
 434   for( int i=lo; i<=up; i++ )
 435   {
 436     double z = points.Value( i ).Z();
 437     if( z>theZMax )
 438       theZMax = z;
 439     if( z<theZMin )
 440       theZMin = z;
 441   }
 442 }
 443
 444 inline gp_Pnt2d To2D( const gp_Pnt& thePnt, const gp_Trsf& theTr,
 445                       double& theUMin, double& theUMax )
 446 {
 447   gp_Pnt p = thePnt.Transformed( theTr );
 448   double u = p.X();
 449   double z = p.Z();
 450   if( u<theUMin )
 451     theUMin = u;
 452   if( u>theUMax )
 453     theUMax = u;
 454   return gp_Pnt2d( u, z );
 455 }
 456
 457 Handle(TColgp_HArray1OfPnt2d) To2D( const TColgp_Array1OfPnt& thePoints,
 458                                     const gp_Trsf& theTr,
 459                                     double& theUMin, double& theUMax )
 460 {
 461   int low = thePoints.Lower(), up = thePoints.Upper();
 462   Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( low, up );
 463   for( int i=low; i<=up; i++ )
 464     points->SetValue( i, To2D( thePoints.Value( i ), theTr, theUMin, theUMax ) );
 465   return points;
 466 }
 467
 468 Handle(Geom2d_Curve) CurveTo2D( const Handle(Geom_Curve)& theCurve,
 469                                 Standard_Real theFirst, Standard_Real theLast,
 470                                 const gp_Trsf& theTr,
 471                                 double& theUMin, double& theUMax )
 472 {
 473   if( theCurve->IsKind( STANDARD_TYPE( Geom_Line ) ) )
 474   {
 475     gp_Pnt aFirstPnt, aLastPnt;
 476     theCurve->D0( theFirst, aFirstPnt );
 477     theCurve->D0( theLast, aLastPnt );
 478
 479     gp_Pnt2d
 480       aFirst2d = To2D( aFirstPnt, theTr, theUMin, theUMax ),
 481       aLast2d = To2D( aLastPnt, theTr, theUMin, theUMax );
 482
 483     gp_Vec2d dir( aFirst2d, aLast2d );
 484     Handle(Geom2d_Line) aLine2d = new Geom2d_Line( aFirst2d, gp_Dir2d( dir.X(), dir.Y() ) );
 485     return new Geom2d_TrimmedCurve( aLine2d, 0, aLast2d.Distance( aFirst2d ) );
 486   }
 487
 488   if( theCurve->IsKind( STANDARD_TYPE( Geom_BSplineCurve ) ) )
 489   {
 490     Handle(Geom_BSplineCurve) aSpline = Handle(Geom_BSplineCurve)::DownCast( theCurve );
 491
 492     Handle(TColgp_HArray1OfPnt2d) poles = To2D( aSpline->Poles(), theTr, theUMin, theUMax );
 493     const TColStd_Array1OfReal& knots = aSpline->Knots();
 494     const TColStd_Array1OfInteger& multiplicities = aSpline->Multiplicities();
 495     int aDegree = aSpline->Degree();
 496
 497     return new Geom2d_BSplineCurve( poles->Array1(), knots, multiplicities, aDegree );
 498   }
 499
 500   return Handle(Geom2d_Curve)();
 501 }
 502
 503 #include <GCE2d_MakeSegment.hxx>
 504 #include <Geom2dAPI_InterCurveCurve.hxx>
 505 bool IsCooriented( const Handle(HYDROData_Profile)& theProfile1,
 506                    const Handle(HYDROData_Profile)& theProfile2 )
 507 {
 508   if( theProfile1==theProfile2 )
 509     return true;
 510
 511   gp_XY lp1, rp1, lp2, rp2;
 512   theProfile1->GetLeftPoint(lp1);
 513   theProfile1->GetRightPoint(rp1);
 514   theProfile2->GetLeftPoint(lp2);
 515   theProfile2->GetRightPoint(rp2);
 516
 517   GCE2d_MakeSegment s1(lp1, lp2);
 518   GCE2d_MakeSegment s2(rp1, rp2);
 519
 520   Geom2dAPI_InterCurveCurve inter;
 521   inter.Init(s1.Value(), s2.Value());
 522   if (inter.NbPoints() == 0)
 523     return true;
 524   else
 525     return false;
 526 }
 527
 528 Handle(Geom2d_BSplineCurve) HYDROData_DTM::CreateHydraulicAxis(
 529   const std::vector<Handle(HYDROData_Profile)>& theProfiles,
 530   std::vector<double>& theDistances )
 531 {
 532   size_t n = theProfiles.size();
 533   if( n==1 )
 534     return Handle_Geom2d_BSplineCurve();
 535
 536   Handle_Geom2d_BSplineCurve aResult;
 537
 538   Handle(TColgp_HArray1OfPnt2d) points = new TColgp_HArray1OfPnt2d( 1, (int)n );
 539   TColgp_Array1OfVec2d tangents( 1, (int)n );
 540   Handle(TColStd_HArray1OfBoolean) flags = new TColStd_HArray1OfBoolean( 1, (int)n );
 541
 542   // Stage 1. Orient all profiles to be co-oriented with the first profile
 543   theProfiles[0]->Update();
 544   for( size_t i = 1; i < n; i++ )
 545   {
 546     Handle(HYDROData_Profile) aProfile = theProfiles[i];
 547     Handle(HYDROData_Profile) aPrevProfile = theProfiles[i-1];
 548
 549     if( !IsCooriented( aProfile, aPrevProfile ) )
 550     {
 551       gp_XY lp, rp;
 552       aProfile->GetLeftPoint( lp, true );
 553       aProfile->GetRightPoint( rp, true );
 554       aProfile->SetLeftPoint( rp, true );
 555       aProfile->SetRightPoint( lp, true );
 556     }
 557     aProfile->Update();
 558   }
 559
 560   // Stage 2. Calculate normals so that each normal "points" to the next profile
 561   for( size_t i = 0; i < n; i++ )
 562   {
 563     Handle(HYDROData_Profile) aProfile = theProfiles[i];
 564     Handle(HYDROData_Profile) aNextProfile = i==n-1 ? theProfiles[i-1] : theProfiles[i+1];
 565
 566     gp_Pnt aLowest;
 567     gp_Vec2d aNormal;
 568     double zmin, zmax;
 569
 570     gp_XYZ curP = aProfile->GetBottomPoint(true);
 571     gp_XY curP2d = gp_XY(curP.X(), curP.Y());
 572
 573     gp_XYZ nextP = aNextProfile->GetBottomPoint(true);
 574     gp_XY nextP2d = gp_XY(nextP.X(), nextP.Y());
 575
 576     gp_Vec2d aPrTangent;
 577     GetProperties( aProfile, aLowest, aPrTangent, zmin, zmax );
 578     aNormal.SetCoord( -aPrTangent.Y(), aPrTangent.X() );
 579
 580     gp_Vec2d aDirToNextProfile(nextP2d.X() - curP2d.X(), nextP2d.Y() - curP2d.Y() );
 581     if( i==n-1 )
 582       aDirToNextProfile.Reverse();
 583
 584     if (aNormal.Dot(aDirToNextProfile) < 0)
 585       aNormal.Reverse();
 586
 587     aNormal.Normalize();
 588
 589     points->SetValue( (int)(i+1), gp_Pnt2d( aLowest.X(), aLowest.Y() ) );
 590     tangents.SetValue( (int)(i+1), aNormal );
 591     flags->SetValue( (int)(i+1), Standard_True );
 592   }
 593
 594   Geom2dAPI_Interpolate anInterpolator( points, Standard_False, Standard_False );
 595   anInterpolator.Load( tangents, flags );
 596   anInterpolator.Perform();
 597   if( anInterpolator.IsDone() )
 598   {
 599     aResult = anInterpolator.Curve();
 600
 601     //fill the distances vector
 602     Geom2dAdaptor_Curve anAdaptor( aResult );
 603
 604     theDistances.clear();
 605     theDistances.reserve( n );
 606     Standard_Real aParamFirst = anAdaptor.FirstParameter(), aParamLast = anAdaptor.LastParameter();
 607     for( size_t i = 1; i <= n; i++ )
 608     {
 609       gp_Pnt2d aPnt = points->Value( (Standard_Integer)i );
 610       Geom2dAPI_ProjectPointOnCurve aProject( aPnt, aResult );
 611       Standard_Real aParam = aProject.LowerDistanceParameter();
 612       double aDistance = GCPnts_AbscissaPoint::Length( anAdaptor, aParamFirst, aParam );
 613       theDistances.push_back( aDistance );
 614     }
 615   }
 616   return aResult;
 617 }
 618
 619 std::vector<Handle(Geom2d_Curve)> HYDROData_DTM::ProfileToParametric(
 620   const Handle(HYDROData_Profile)& theProfile,
 621   double& theUMin, double& theUMax, gp_Vec2d& theDir )
 622 {
 623   std::vector<Handle(Geom2d_Curve)> curves;
 624
 625   // Transformation of the coordinate systems
 626   gp_Pnt aLowest;
 627   double zmin, zmax;
 628   GetProperties( theProfile, aLowest, theDir, zmin, zmax );
 629
 630   gp_Ax3 aStd3d( gp_Pnt( 0, 0, 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( 1, 0, 0 ) );
 631   gp_Ax3 aLocal( gp_Pnt( aLowest.X(), aLowest.Y(), 0 ), gp_Dir( 0, 0, 1 ), gp_Dir( theDir.X(), theDir.Y(), 0 ) );
 632
 633   gp_Trsf aTransf;
 634   aTransf.SetTransformation( aStd3d, aLocal );
 635
 636   // Iteration via edges
 637   TopoDS_Wire aWire = TopoDS::Wire( theProfile->GetShape3D() );
 638   TopExp_Explorer anExp( aWire, TopAbs_EDGE );
 639   for( ; anExp.More(); anExp.Next() )
 640   {
 641     // Extract an edge from wire
 642     TopoDS_Edge anEdge = TopoDS::Edge( anExp.Current() );
 643
 644     // Extract a curve corresponding to the edge
 645     TopLoc_Location aLoc;
 646     Standard_Real aFirst, aLast;
 647     Handle(Geom_Curve) aCurve = BRep_Tool::Curve( anEdge, aLoc, aFirst, aLast );
 648
 649     // Convert the curve to 2d CS
 650     Handle(Geom2d_Curve) aCurve2d = CurveTo2D( aCurve, aFirst, aLast, aTransf, theUMin, theUMax );
 651     if( !aCurve2d.IsNull() )
 652       curves.push_back( aCurve2d );
 653   }
 654   return curves;
 655 }
 656
 657
 658 bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
 659 {
 660   double umin = DBL_MAX,
 661          umax = -umin;
 662
 663   size_t n = intersections.size();
 664   if( n <= 0 )
 665     return false;
 666
 667   std::set<double>::const_iterator it = intersections.begin(), last = intersections.end();
 668   for( ; it!=last; it++ )
 669   {
 670     double u = *it;
 671     if( u<umin )
 672       umin = u;
 673     if( u>umax )
 674       umax = u;
 675   }
 676   theMid = ( umin+umax )/2;
 677   theWid = umax-umin;
 678   return true;
 679 }
 680
 681 void HYDROData_DTM::ProfileDiscretization( const Handle(HYDROData_Profile)& theProfile,
 682                                            double theXCurv, double theMinZ, double theMaxZ, double theTopZ, double theDDZ,
 683                                            CurveUZ& theMidPointCurve,
 684                                            CurveUZ& theWidthCurve,
 685                                            int& intersection_nb,
 686                                            double theTolerance)
 687 {
 688   double aDblMax = DBL_MAX,
 689          aUMin = aDblMax,
 690          aUMax = -aUMin,
 691          aVMax = 1000000;
 692
 693   gp_Vec2d aProfileDir;
 694   std::vector<Handle(Geom2d_Curve)> curves = ProfileToParametric( theProfile, aUMin, aUMax, aProfileDir );
 695   size_t n = curves.size();
 696
 697   if( n==0 )
 698     return;
 699
 700   // we add the "virtual" vertical lines to simulate the intersection with profile
 701   gp_Pnt2d aFirst, aLast;
 702   curves[0]->D0( curves[0]->FirstParameter(), aFirst );
 703   curves[n-1]->D0( curves[n-1]->LastParameter(), aLast );
 704   Handle(Geom2d_Line) aV1 = new Geom2d_Line( aFirst, gp_Dir2d( 0, 1 ) );
 705   Handle(Geom2d_TrimmedCurve) aT1 = new Geom2d_TrimmedCurve( aV1, 0.0, aVMax );
 706
 707   Handle(Geom2d_Line) aV2 = new Geom2d_Line( aLast, gp_Dir2d( 0, 1 ) );
 708   Handle(Geom2d_TrimmedCurve) aT2 = new Geom2d_TrimmedCurve( aV2, 0.0, aVMax );
 709
 710   curves.push_back( aT1 );
 711   curves.push_back( aT2 );
 712
 713   int psize = ( int )( ( theMaxZ-theMinZ ) / theDDZ + 1 );
 714   theMidPointCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ);
 715   theMidPointCurve.reserve( psize );
 716   theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ );
 717   theWidthCurve.reserve( psize );
 718
 719   n = curves.size();
 720   // for each discrete value of z we search intersection with profile
 721   for( double z1 = theMinZ; z1 <= theMaxZ; z1 += theDDZ )
 722   {
 723     Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z1 ), gp_Dir2d( 1, 0 ) );
 724     std::set<double> intersections;
 725     for( size_t i = 0; i < n; i++ )
 726     {
 727       Handle(Geom2d_Curve) aCurve = curves[i];
 728       Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
 729       for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
 730         intersections.insert( anIntersect.Point( k ).X() );
 731     }
 732
 733     intersection_nb = intersections.size();
 734     if( intersection_nb >= 1 )
 735     {
 736       double u_mid, u_wid;
 737       if( !CalcMidWidth( intersections, u_mid, u_wid ) )
 738         continue;
 739
 740       double z = z1 - theMinZ;
 741       PointUZ p_mid;
 742       p_mid.U = u_mid;
 743       p_mid.Z = z;
 744       theMidPointCurve.push_back( p_mid );
 745
 746       PointUZ p_wid;
 747       p_wid.U = u_wid;
 748       p_wid.Z = z;
 749       theWidthCurve.push_back( p_wid );
 750     }
 751   }
 752 }
 753
 754 void HYDROData_DTM::Interpolate( const CurveUZ& theCurveA, const CurveUZ& theCurveB,
 755                                  int theNbSteps, std::vector<CurveUZ>& theInterpolation,
 756                                  bool isAddSecond )
 757 {
 758   theInterpolation.clear();
 759   int d = isAddSecond ? 2 : 1;
 760   theInterpolation.reserve( theNbSteps+d );
 761   double dt = 1.0 / double( theNbSteps + 1 );
 762   double t = dt;
 763   theInterpolation.push_back( theCurveA );
 764   for( int i=0; i<theNbSteps; i++, t+=dt )
 765   {
 766     CurveUZ anInterp = theCurveA*(1-t) + theCurveB*t;
 767     theInterpolation.push_back( anInterp );
 768   }
 769   if( isAddSecond )
 770     theInterpolation.push_back( theCurveB );
 771 }
 772 #include <BRepLib_MakeEdge2d.hxx>
 773 void HYDROData_DTM::CurveTo3D( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
 774                                const CurveUZ& theMidCurve, const CurveUZ& theWidthCurve,
 775                                AltitudePoints& thePoints )
 776 {
 777   Geom2dAdaptor_Curve anAdaptor( theHydraulicAxis );
 778   TopoDS_Edge E2d = BRepLib_MakeEdge2d(theHydraulicAxis).Edge();
 779   GCPnts_AbscissaPoint ap( anAdaptor, theMidCurve.Xcurv(), anAdaptor.FirstParameter() );
 780   double aParam = ap.Parameter();
 781
 782   gp_Pnt2d point;
 783   anAdaptor.D0( aParam, point );
 784   gp_Vec2d profile_dir = theMidCurve.ProfileDir();
 785   //gp_Dir tangent_n( -profile_dir.Y(), profile_dir.X(), dz );
 786   profile_dir.Normalize();
 787
 788   size_t n = theMidCurve.size();
 789   std::map<double, AltitudePoint> sorted_points;
 790   bool isOnTop = false;
 791   for( size_t i=0; i<n; i++ ) // build the two banks of the interpolated profile, from bottom to top
 792   {
 793     double param1 = theMidCurve[i].U - theWidthCurve[i].U / 2;
 794     double param2 = theMidCurve[i].U + theWidthCurve[i].U / 2;
 795
 796     gp_Pnt2d p1 = point.Translated( param1 * profile_dir);
 797     gp_Pnt2d p2 = point.Translated( param2 * profile_dir);
 798
 799     bool arrivedOnTop = false;
 800     double z = 0;
 801     if (theMidCurve[i].Z <= theMidCurve.MaxZ())
 802       z = theMidCurve[i].Z + theMidCurve.DeltaZ();
 803     else
 804       {
 805         z = theMidCurve.MaxZ() + theMidCurve.DeltaZ(); // limit z to linear interpolation between maxima on extremity profiles
 806         arrivedOnTop = true; // do not keep points after this one
 807       }
 808     if (!isOnTop)
 809       {
 810         AltitudePoint p3d_1( p1.X(), p1.Y(), z ), p3d_2( p2.X(), p2.Y(), z );
 811
 812         sorted_points[param1] = p3d_1;
 813         sorted_points[param2] = p3d_2;
 814       }
 815     //if (arrivedOnTop)
 816     //  isOnTop =true; // do not keep points after this one (commented: leads to strange limits of 2D shape)
 817   }
 818
 819   thePoints.reserve( sorted_points.size() );
 820   const double EPS = 1E-12;
 821   std::map<double, AltitudePoint>::const_iterator it = sorted_points.begin(), last = sorted_points.end();
 822   for( ; it!=last; it++ )
 823     if( thePoints.empty() || thePoints.back().SquareDistance( it->second ) > EPS )
 824       thePoints.push_back( it->second );
 825 }
 826
 827 //inline double max( double a, double b )
 828 //{
 829 //  if( a>b )
 830 //    return a;
 831 //  else
 832 //    return b;
 833 //}
 834 //
 835 //inline double min( double a, double b )
 836 //{
 837 //  if( a<b )
 838 //    return a;
 839 //  else
 840 //    return b;
 841 //}
 842
 843 #include <BRepLib_MakeWire.hxx>
 844
 845 std::vector<HYDROData_Bathymetry::AltitudePoints> HYDROData_DTM::Interpolate
 846   ( const Handle(Geom2d_BSplineCurve)& theHydraulicAxis,
 847     const Handle(HYDROData_Profile)& theProfileA,
 848     double theXCurvA,
 849     const Handle(HYDROData_Profile)& theProfileB,
 850     double theXCurvB,
 851     double theDDZ, int theNbSteps, bool isAddSecond,
 852     int& inter_nb_1, int& inter_nb_2)
 853 {
 854   double zminA, zmaxA, zminB, zmaxB;
 855   gp_Pnt lowestA, lowestB;
 856   gp_Vec2d dirA, dirB;
 857
 858   GetProperties( theProfileA, lowestA, dirA, zminA, zmaxA );
 859   GetProperties( theProfileB, lowestB, dirB, zminB, zmaxB );
 860
 861
 862   double hmax = max( zmaxA-zminA, zmaxB-zminB );
 863
 864   //double dz = zminB - zminA;
 865   //double zmin = min( zminA, zminB );
 866   //double zmax = max( zmaxA, zmaxB );
 867
 868   CurveUZ midA(0, gp_Vec2d(), 0, 0), midB(0, gp_Vec2d(), 0, 0);
 869   CurveUZ widA(0, gp_Vec2d(), 0, 0), widB(0, gp_Vec2d(), 0, 0);
 870
 871   ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, zmaxA-zminA, theDDZ, midA, widA, inter_nb_1 );
 872   ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, zmaxB-zminB, theDDZ, midB, widB, inter_nb_2 );
 873
 874   std::vector<CurveUZ> mid, wid;
 875   Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
 876   Interpolate( widA, widB, theNbSteps, wid, isAddSecond );
 877
 878   size_t p = mid.size();
 879   size_t q = p>0 ? 2*mid[0].size() : 1;
 880   std::vector<AltitudePoints> points;
 881   points.resize( p );
 882
 883   for( size_t i=0; i<p; i++ )
 884   {
 885     points[i].reserve( q );
 886     CurveTo3D( theHydraulicAxis, mid[i], wid[i], points[i] );
 887   }
 888
 889   return points;
 890 }
 891
 892 HYDROData_Bathymetry::AltitudePoints HYDROData_DTM::Interpolate
 893   ( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
 894     double theDDZ, double theSpatialStep,
 895     AltitudePoints& theLeft,
 896     AltitudePoints& theRight,
 897     std::vector<AltitudePoints>& theMainProfiles,
 898     std::set<int>& invalInd)
 899 {
 900   AltitudePoints points;
 901   size_t n = theProfiles.size();
 902   if( n<=1 )
 903     return points;
 904
 905   std::vector<double> distances;
 906   Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
 907   if( aHydraulicAxis.IsNull() )
 908     return points;
 909
 910   theMainProfiles.reserve( n );
 911
 912   for( size_t i=0, n1=n-1; i<n1; i++ )
 913   {
 914     double aDistance = distances[i+1]-distances[i];
 915     int aNbSteps = int(aDistance/theSpatialStep);
 916     bool isAddSecond = i==n1-1;
 917
 918     // 1. Calculate interpolated profiles
 919     int inter_nb_1, inter_nb_2;
 920     std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
 921       theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, inter_nb_1, inter_nb_2 );
 922     int lps = local_points.size();
 923
 924     if (inter_nb_1 > 2)
 925       invalInd.insert(i);
 926
 927     if (inter_nb_2 > 2)
 928       invalInd.insert(i+1);
 929
 930     // 2. Put all points into the global container
 931     for( size_t j=0; j<lps; j++ )
 932     {
 933       const AltitudePoints& lp = local_points[j];
 934       if( i==0 && j==0 )
 935         points.reserve( lp.size() * n );
 936       for( size_t k=0, ks=lp.size(); k<ks; k++ )
 937         points.push_back( lp[k] );
 938     }
 939
 940     // 3. Get left/right banks' points
 941     if( i==0 )
 942     {
 943       theLeft.reserve( lps * n );
 944       theRight.reserve( lps * n );
 945     }
 946     for( size_t j=0; j<lps; j++ )
 947     {
 948       const AltitudePoints& lp = local_points[j];
 949       theLeft.push_back( lp[0] );
 950       theRight.push_back( lp[lp.size()-1] );
 951     }
 952
 953     // 4. Get main profiles points
 954     theMainProfiles.push_back( local_points[0] );
 955     if( isAddSecond )
 956       theMainProfiles.push_back( local_points[lps-1] );
 957   }
 958   return points;
 959 }
 960
 961 int HYDROData_DTM::EstimateNbPoints( const std::vector<Handle(HYDROData_Profile)>& theProfiles,
 962                                      double theDDZ, double theSpatialStep )
 963 {
 964   size_t n = theProfiles.size();
 965   if( n<=1 )
 966     return 0;
 967   if( theDDZ<1E-6 || theSpatialStep<1E-6 )
 968     return 1 << 20;
 969
 970   std::vector<double> distances;
 971   Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
 972   if( aHydraulicAxis.IsNull() )
 973     return 0;
 974
 975   double aCompleteDistance = distances[n-1];
 976   int aNbSteps = int( aCompleteDistance / theSpatialStep ) + 1;
 977   gp_Pnt aLowest;
 978   gp_Vec2d aDir;
 979   double aZMin, aZMax;
 980   GetProperties( theProfiles[0], aLowest, aDir, aZMin, aZMax );
 981   int aNbZSteps = (aZMax-aZMin)/theDDZ;
 982
 983   if( aNbSteps > ( 1<<16 ) || aNbZSteps > ( 1<<16 ) )
 984     return 1 << 20;
 985
 986   return aNbSteps * aNbZSteps;
 987 }