src/HYDROData/HYDROData_StreamAltitude.cxx

   1
   2 #include "HYDROData_StreamAltitude.h"
   3
   4 #include "HYDROData_Document.h"
   5 #include "HYDROData_Profile.h"
   6 #include "HYDROData_Stream.h"
   7 #include "HYDROData_ShapesTool.h"
   8
   9 #include <BRep_Tool.hxx>
  10
  11 #include <BRepTopAdaptor_FClass2d.hxx>
  12
  13 #include <BRepBuilderAPI_MakeEdge.hxx>
  14 #include <BRepBuilderAPI_MakeFace.hxx>
  15 #include <BRepBuilderAPI_MakeWire.hxx>
  16
  17 #include <Extrema_ExtElC.hxx>
  18
  19 #include <GeomAPI_ProjectPointOnCurve.hxx>
  20
  21 #include <gp_Lin.hxx>
  22
  23 #include <Precision.hxx>
  24
  25 #include <TopExp_Explorer.hxx>
  26
  27 #include <TopoDS.hxx>
  28 #include <TopoDS_Wire.hxx>
  29 #include <TopoDS_Edge.hxx>
  30
  31 #include <TopTools_SequenceOfShape.hxx>
  32
  33 #include <QStringList>
  34 #define CLASS2D 1
  35 #ifdef CLASS2D
  36 #include <gp_Pln.hxx>
  37 #include <TopExp.hxx>
  38 #include <Geom_Plane.hxx>
  39 #include <BRepBuilderAPI_FindPlane.hxx>
  40 #include <TopoDS_Vertex.hxx>
  41 #include <TColgp_SequenceOfVec.hxx>
  42 #include <TopTools_ShapeMapHasher.hxx>
  43 #undef _NCollection_MapHasher
  44 #endif
  45 #ifdef DEB_CLASS2D
  46 #include <BRepTools.hxx>
  47 #include <BRep_Builder.hxx>
  48 #include <BRepBuilderAPI_MakeVertex.hxx>
  49 #endif
  50 IMPLEMENT_STANDARD_HANDLE(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
  51 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
  52
  53 HYDROData_StreamAltitude::HYDROData_StreamAltitude()
  54 : HYDROData_IAltitudeObject()
  55 {
  56 }
  57
  58 HYDROData_StreamAltitude::~HYDROData_StreamAltitude()
  59 {
  60 }
  61
  62 QStringList HYDROData_StreamAltitude::DumpToPython( MapOfTreatedObjects& theTreatedObjects ) const
  63 {
  64   QStringList aResList = dumpObjectCreation( theTreatedObjects );
  65   QString aName = GetObjPyName();
  66
  67   // TODO
  68
  69   aResList << QString( "" );
  70   aResList << QString( "%1.Update();" ).arg( aName );
  71   aResList << QString( "" );
  72
  73   return aResList;
  74 }
  75
  76 Standard_Real getAltitudeFromProfile( const Handle(HYDROData_Profile)& theProfile,
  77                                       const Standard_Real&             theLeftDist,
  78                                       const Standard_Real&             theRightDist )
  79 {
  80   Standard_Real aResAlt = 0.0;
  81
  82   gp_XY aFirstPoint, aLastPoint;
  83   if ( !theProfile->GetLeftPoint( aFirstPoint ) ||
  84        !theProfile->GetRightPoint( aLastPoint ) )
  85     return aResAlt;
  86
  87   gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
  88   gp_Pnt aPnt2( aLastPoint.X(),  aLastPoint.Y(),  0 );
  89
  90   Standard_Real aProfileDist = aPnt1.Distance( aPnt2 );
  91
  92   Standard_Real aCoeff = aProfileDist / ( theLeftDist + theRightDist );
  93
  94   gp_Pnt anIntPoint( aPnt1.XYZ() + ( aCoeff * theLeftDist ) * gp_Dir( gp_Vec( aPnt1, aPnt2 ) ).XYZ() );
  95
  96   gp_Lin aPointLine( anIntPoint, gp::DZ() );
  97
  98   gp_Pnt aPrevPoint;
  99   gp_Lin aPrevNormal;
 100   HYDROData_Profile::ProfilePoints aProfilePoints = theProfile->GetProfilePoints();
 101   for ( int i = 1, n = aProfilePoints.Length(); i <= n; ++i )
 102   {
 103     gp_Pnt aProfPoint( aProfilePoints.Value( i ) );
 104
 105     Standard_Real aDist = aPointLine.Distance( aProfPoint );
 106     if ( aDist <= gp::Resolution() )
 107     {
 108       // We found the intersected point
 109       aResAlt = aProfPoint.Z();
 110       break;
 111     }
 112
 113     gp_Lin aNormal = aPointLine.Normal( aProfPoint );
 114     if ( i == 1 )
 115     {
 116       aPrevNormal = aNormal;
 117       aPrevPoint = aProfPoint;
 118       continue;
 119     }
 120
 121     if ( aPrevNormal.Direction().Dot( aNormal.Direction() ) < 0 )
 122     {
 123       // We found the intersected edge
 124       gp_Lin anEdgeLine( aPrevPoint, gp_Dir( gp_Vec( aPrevPoint, aProfPoint ) ) );
 125
 126       Extrema_ExtElC anExtrema( aPointLine, anEdgeLine, Precision::Angular() );
 127       if ( !anExtrema.IsParallel() )
 128       {
 129         Extrema_POnCurv aFirstPnt, aSecPnt;
 130         anExtrema.Points( 1, aFirstPnt, aSecPnt );
 131
 132         const gp_Pnt& anIntPnt = aSecPnt.Value();
 133         aResAlt = anIntPnt.Z();
 134
 135         break;
 136       }
 137     }
 138
 139     aPrevNormal = aNormal;
 140     aPrevPoint = aProfPoint;
 141   }
 142
 143   return aResAlt;
 144 }
 145
 146 // =================================== workaround ==============================
 147 // till fix of Classifier2d bug
 148 #ifdef CLASS2D
 149 TopAbs_State ComputePointState(const gp_XY& thePnt2d, const TopoDS_Face& theFace)
 150 {
 151   TopAbs_State aStatus(TopAbs_UNKNOWN);
 152   TopoDS_Wire aWire;
 153   int nb(0);
 154   TopoDS_Iterator it(theFace);
 155   for(;it.More();it.Next()) {
 156           aWire = TopoDS::Wire(it.Value());
 157           nb++;
 158   }
 159   if(nb > 1 || aWire.IsNull()) return aStatus;
 160   gp_Pln aPlane;
 161   gp_Vec aNormal;
 162   BRepBuilderAPI_FindPlane fndPlane (theFace, Precision::Confusion());
 163   if(fndPlane.Found())
 164     aPlane = fndPlane.Plane()->Pln();
 165   else
 166     return aStatus;
 167   aNormal = gp_Vec(aPlane.Axis().Direction());
 168   if(theFace.Orientation() == TopAbs_REVERSED)
 169     aNormal.Reverse();
 170   gp_Pnt aPoint = gp_Pnt (thePnt2d.X(), thePnt2d.Y(), 0);
 171   TColgp_SequenceOfVec aSeq;
 172   TopTools_MapOfShape aMap;
 173   it.Initialize(aWire);
 174   for (;it.More(); it.Next()) {
 175         const TopoDS_Vertex& aVx = TopExp::FirstVertex(TopoDS::Edge(it.Value()), Standard_True);
 176         if(!aMap.Add(aVx))
 177           continue;
 178         const gp_Pnt& aCurPnt = BRep_Tool::Pnt(TopoDS::Vertex(aVx));
 179         if(aPoint.IsEqual(aCurPnt,Precision::Confusion())) {
 180       aStatus = TopAbs_ON;
 181           return aStatus;
 182         }
 183         gp_Vec aVec (aPoint, aCurPnt);
 184         aSeq.Append(aVec);
 185   }
 186   Standard_Real anAng(0.0);
 187   for(int i = 1;i < aSeq.Length();i++) {
 188     const gp_Vec& aV1 = aSeq.Value(i);
 189         const gp_Vec& aV2 = aSeq.Value(i+1);
 190          anAng += aV1.AngleWithRef(aV2, aNormal);
 191   }
 192   anAng += aSeq.Value(aSeq.Length()).AngleWithRef(aSeq.Value(1), aNormal);
 193   if(abs(anAng) > Precision::Angular())
 194     aStatus = TopAbs_IN;
 195   else
 196          aStatus = TopAbs_OUT;
 197   return aStatus;
 198 }
 199 #endif
 200 bool HYDROData_StreamAltitude::getBoundaryProfilesForPoint(
 201   const gp_XY&               thePoint,
 202   Handle(HYDROData_Profile)& theLeftProfile,
 203   Handle(HYDROData_Profile)& theRightProfile ) const
 204 {
 205   Handle(HYDROData_Stream) aStream =
 206     Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
 207   if ( aStream.IsNull() )
 208     return false;
 209
 210   HYDROData_SequenceOfObjects aStreamProfiles = aStream->GetProfiles();
 211   if ( aStreamProfiles.Length() < 2 )
 212     return false;
 213
 214   Handle(HYDROData_Profile) aPrevProfile;
 215   gp_Pnt aPrevPnt1, aPrevPnt2;
 216   for ( int i = 1, n = aStreamProfiles.Length(); i <= n; ++i )
 217   {
 218     Handle(HYDROData_Profile) aProfile =
 219       Handle(HYDROData_Profile)::DownCast( aStreamProfiles.Value( i ) );
 220     if ( aProfile.IsNull() )
 221       continue;
 222
 223     gp_XY aFirstPoint, aLastPoint;
 224     if ( !aProfile->GetLeftPoint( aFirstPoint ) || !aProfile->GetRightPoint( aLastPoint ) )
 225       continue;
 226
 227     gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
 228     gp_Pnt aPnt2( aLastPoint.X(),  aLastPoint.Y(),  0 );
 229
 230     if ( !aPrevProfile.IsNull() )
 231     {
 232       BRepBuilderAPI_MakeEdge aLeftMakeEdge( aPrevPnt1, aPrevPnt2 );
 233       BRepBuilderAPI_MakeEdge aBotMakeEdge( aPrevPnt2, aPnt2 );
 234       BRepBuilderAPI_MakeEdge aRightMakeEdge( aPnt2, aPnt1 );
 235       BRepBuilderAPI_MakeEdge aTopMakeEdge( aPnt1, aPrevPnt1 );
 236
 237       BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(),
 238                                          aRightMakeEdge.Edge(), aTopMakeEdge.Edge() );
 239
 240       BRepBuilderAPI_MakeFace aMakeFace( aMakeWire.Wire() );
 241
 242       TopoDS_Face aProfilesFace = aMakeFace.Face();
 243 #ifdef DEB_CLASS2D
 244           TopoDS_Compound aCmp;
 245       BRep_Builder aBB;
 246       aBB.MakeCompound(aCmp);
 247           aBB.Add(aCmp, aProfilesFace);
 248           gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
 249           BRepBuilderAPI_MakeVertex aMk(aPnt);
 250           aBB.Add(aCmp, aMk.Vertex());
 251           BRepTools::Write(aCmp, "ProfileFace.brep");
 252 #endif
 253       BRepTopAdaptor_FClass2d aClassifier( aProfilesFace, Precision::Confusion() );
 254           TopAbs_State aPointState =
 255 #ifdef CLASS2D
 256          ComputePointState(thePoint, aProfilesFace);
 257 #else
 258           aClassifier.Perform( gp_Pnt2d( thePoint ), Standard_False );
 259 #endif
 260 #ifdef DEB_CLASS2D
 261           cout << "Point status is = " << aPointState <<endl;
 262 #endif
 263       if ( aPointState != TopAbs_OUT )
 264       {
 265         theLeftProfile = aPrevProfile;
 266         theRightProfile = aProfile;
 267         break;
 268       }
 269     }
 270
 271     aPrevProfile = aProfile;
 272     aPrevPnt1 = aPnt1;
 273     aPrevPnt2 = aPnt2;
 274   }
 275
 276   return !theLeftProfile.IsNull() && !theRightProfile.IsNull();
 277 }
 278
 279 double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) const
 280 {
 281   double aResAltitude = GetInvalidAltitude();
 282
 283   Handle(HYDROData_Stream) aStream =
 284     Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
 285   if ( aStream.IsNull() )
 286     return aResAltitude;
 287
 288   TopoDS_Shape aStreamShape = aStream->GetTopShape();
 289   if ( aStreamShape.IsNull() )
 290     return aResAltitude;
 291
 292   TopExp_Explorer aStreamFaceExp( aStreamShape, TopAbs_FACE );
 293   if ( !aStreamFaceExp.More() )
 294     return aResAltitude;
 295
 296   // Get only face because of 2d profile wires is in compound
 297   TopoDS_Face aStreamFace = TopoDS::Face( aStreamFaceExp.Current() );
 298
 299   // Check if point is inside of stream presentation
 300   BRepTopAdaptor_FClass2d aClassifier( aStreamFace, Precision::Confusion() );
 301   TopAbs_State aPointState =
 302 #ifdef CLASS2D
 303     ComputePointState(thePoint, aStreamFace);
 304 #else
 305     aClassifier.Perform( gp_Pnt2d( thePoint ), Standard_False );
 306 #endif
 307
 308 #ifdef DEB_CLASS2D
 309       cout << "Point status is = " << aPointState <<endl;
 310           TopoDS_Compound aCmp;
 311       BRep_Builder aBB;
 312       aBB.MakeCompound(aCmp);
 313           aBB.Add(aCmp, aStreamFace);
 314           gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
 315           BRepBuilderAPI_MakeVertex aMk(aPnt);
 316           aBB.Add(aCmp, aMk.Vertex());
 317           BRepTools::Write(aCmp, "FCL2d.brep");
 318 #endif
 319   if ( aPointState == TopAbs_OUT )
 320     return aResAltitude;
 321
 322   // Find the two profiles between which the point is lies
 323   Handle(HYDROData_Profile) aLeftProfile, aRightProfile;
 324   if ( !getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile ) )
 325     return aResAltitude;
 326
 327   // Find the projections of point to borders of stream
 328   gp_XYZ aPointToTest( thePoint.X(), thePoint.Y(), 0.0 );
 329
 330   TopoDS_Edge aStreamLeftEdge = TopoDS::Edge( aStream->GetLeftShape() );
 331   TopoDS_Edge aStreamRightEdge = TopoDS::Edge( aStream->GetRightShape() );
 332
 333   Standard_Real aFirst = 0.0, aLast = 0.0;
 334   Handle(Geom_Curve) anEdgeLeftCurve = BRep_Tool::Curve( aStreamLeftEdge, aFirst, aLast );
 335   Handle(Geom_Curve) anEdgeRightCurve = BRep_Tool::Curve( aStreamRightEdge, aFirst, aLast );
 336
 337   GeomAPI_ProjectPointOnCurve aLeftProject( aPointToTest, anEdgeLeftCurve );
 338   GeomAPI_ProjectPointOnCurve aRightProject( aPointToTest, anEdgeRightCurve );
 339
 340   Standard_Real aLeftDist = aLeftProject.LowerDistance();
 341   Standard_Real aRightDist = aRightProject.LowerDistance();
 342
 343   // Find the altitude in profiles
 344   Standard_Real aLeftAlt = getAltitudeFromProfile( aLeftProfile, aLeftDist, aRightDist );
 345   Standard_Real aRightAlt = getAltitudeFromProfile( aRightProfile, aLeftDist, aRightDist );
 346
 347   // Interpolate altitudes
 348   Standard_Real aFirstCoeff = aLeftDist / ( aLeftDist + aRightDist );
 349   Standard_Real aSecCoeff = aRightDist / ( aLeftDist + aRightDist );
 350
 351   aResAltitude = aLeftAlt * aFirstCoeff + aRightAlt * aSecCoeff;
 352
 353   return aResAltitude;
 354 }
 355
 356
 357
 358