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
  35 #define PYTHON_STREAM_ALTITUDE_ID "KIND_STREAM_ALTITUDE"
  36
  37
  38 IMPLEMENT_STANDARD_HANDLE(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
  39 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
  40
  41 HYDROData_StreamAltitude::HYDROData_StreamAltitude()
  42 : HYDROData_IAltitudeObject()
  43 {
  44 }
  45
  46 HYDROData_StreamAltitude::~HYDROData_StreamAltitude()
  47 {
  48 }
  49
  50 QStringList HYDROData_StreamAltitude::DumpToPython( MapOfTreatedObjects& theTreatedObjects ) const
  51 {
  52   QStringList aResList;
  53
  54   Handle(HYDROData_Document) aDocument = HYDROData_Document::Document( myLab );
  55   if ( aDocument.IsNull() )
  56     return aResList;
  57
  58   QString aDocName = aDocument->GetDocPyName();
  59   QString anAltitudeName = GetName();
  60
  61   aResList << QString( "%1 = %2.CreateObject( %3 );" )
  62               .arg( anAltitudeName ).arg( aDocName ).arg( PYTHON_STREAM_ALTITUDE_ID );
  63   aResList << QString( "%1.SetName( \"%2\" );" )
  64               .arg( anAltitudeName ).arg( anAltitudeName );
  65
  66   // TODO
  67
  68   return aResList;
  69 }
  70
  71 Standard_Real getAltitudeFromProfile( const Handle(HYDROData_Profile)& theProfile,
  72                                       const Standard_Real&             theLeftDist,
  73                                       const Standard_Real&             theRightDist )
  74 {
  75   Standard_Real aResAlt = 0.0;
  76
  77   gp_XY aFirstPoint, aLastPoint;
  78   if ( !theProfile->GetLeftPoint( aFirstPoint ) ||
  79        !theProfile->GetRightPoint( aLastPoint ) )
  80     return aResAlt;
  81
  82   gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
  83   gp_Pnt aPnt2( aLastPoint.X(),  aLastPoint.Y(),  0 );
  84
  85   Standard_Real aProfileDist = aPnt1.Distance( aPnt2 );
  86
  87   Standard_Real aCoeff = aProfileDist / ( theLeftDist + theRightDist );
  88
  89   gp_Pnt anIntPoint( aPnt1.XYZ() + ( aCoeff * theLeftDist ) * gp_Dir( gp_Vec( aPnt1, aPnt2 ) ).XYZ() );
  90
  91   gp_Lin aPointLine( anIntPoint, gp::DZ() );
  92
  93   gp_Pnt aPrevPoint;
  94   gp_Lin aPrevNormal;
  95   HYDROData_Profile::ProfilePoints aProfilePoints = theProfile->GetProfilePoints();
  96   for ( int i = 1, n = aProfilePoints.Length(); i <= n; ++i )
  97   {
  98     gp_Pnt aProfPoint( aProfilePoints.Value( i ) );
  99
 100     Standard_Real aDist = aPointLine.Distance( aProfPoint );
 101     if ( aDist <= gp::Resolution() )
 102     {
 103       // We found the intersected point
 104       aResAlt = aProfPoint.Z();
 105       break;
 106     }
 107
 108     gp_Lin aNormal = aPointLine.Normal( aProfPoint );
 109     if ( i == 1 )
 110     {
 111       aPrevNormal = aNormal;
 112       aPrevPoint = aProfPoint;
 113       continue;
 114     }
 115
 116     if ( aPrevNormal.Direction().Dot( aNormal.Direction() ) < 0 )
 117     {
 118       // We found the intersected edge
 119       gp_Lin anEdgeLine( aPrevPoint, gp_Dir( gp_Vec( aPrevPoint, aProfPoint ) ) );
 120
 121       Extrema_ExtElC anExtrema( aPointLine, anEdgeLine, Precision::Angular() );
 122       if ( !anExtrema.IsParallel() )
 123       {
 124         Extrema_POnCurv aFirstPnt, aSecPnt;
 125         anExtrema.Points( 1, aFirstPnt, aSecPnt );
 126
 127         const gp_Pnt& anIntPnt = aSecPnt.Value();
 128         aResAlt = anIntPnt.Z();
 129
 130         break;
 131       }
 132     }
 133
 134     aPrevNormal = aNormal;
 135     aPrevPoint = aProfPoint;
 136   }
 137
 138   return aResAlt;
 139 }
 140
 141 bool HYDROData_StreamAltitude::getBoundaryProfilesForPoint(
 142   const gp_XY&               thePoint,
 143   Handle(HYDROData_Profile)& theLeftProfile,
 144   Handle(HYDROData_Profile)& theRightProfile ) const
 145 {
 146   Handle(HYDROData_Stream) aStream =
 147     Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
 148   if ( aStream.IsNull() )
 149     return false;
 150
 151   HYDROData_SequenceOfObjects aStreamProfiles = aStream->GetProfiles();
 152   if ( aStreamProfiles.Length() < 2 )
 153     return false;
 154
 155   Handle(HYDROData_Profile) aPrevProfile;
 156   gp_Pnt aPrevPnt1, aPrevPnt2;
 157   for ( int i = 1, n = aStreamProfiles.Length(); i <= n; ++i )
 158   {
 159     Handle(HYDROData_Profile) aProfile =
 160       Handle(HYDROData_Profile)::DownCast( aStreamProfiles.Value( i ) );
 161     if ( aProfile.IsNull() )
 162       continue;
 163
 164     gp_XY aFirstPoint, aLastPoint;
 165     if ( !aProfile->GetLeftPoint( aFirstPoint ) || !aProfile->GetRightPoint( aLastPoint ) )
 166       continue;
 167
 168     gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
 169     gp_Pnt aPnt2( aLastPoint.X(),  aLastPoint.Y(),  0 );
 170
 171     if ( !aPrevProfile.IsNull() )
 172     {
 173       BRepBuilderAPI_MakeEdge aLeftMakeEdge( aPrevPnt1, aPrevPnt2 );
 174       BRepBuilderAPI_MakeEdge aBotMakeEdge( aPrevPnt2, aPnt2 );
 175       BRepBuilderAPI_MakeEdge aRightMakeEdge( aPnt2, aPnt1 );
 176       BRepBuilderAPI_MakeEdge aTopMakeEdge( aPnt1, aPrevPnt1 );
 177
 178       BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(),
 179                                          aRightMakeEdge.Edge(), aTopMakeEdge.Edge() );
 180
 181       BRepBuilderAPI_MakeFace aMakeFace( aMakeWire.Wire() );
 182
 183       TopoDS_Face aProfilesFace = aMakeFace.Face();
 184
 185       BRepTopAdaptor_FClass2d aClassifier( aProfilesFace, Precision::Confusion() );
 186       TopAbs_State aPointState = aClassifier.Perform( gp_Pnt2d( thePoint ), Standard_False );
 187       if ( aPointState != TopAbs_OUT )
 188       {
 189         theLeftProfile = aPrevProfile;
 190         theRightProfile = aProfile;
 191         break;
 192       }
 193     }
 194
 195     aPrevProfile = aProfile;
 196     aPrevPnt1 = aPnt1;
 197     aPrevPnt2 = aPnt2;
 198   }
 199
 200   return !theLeftProfile.IsNull() && !theRightProfile.IsNull();
 201 }
 202
 203 double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) const
 204 {
 205   double aResAltitude = GetInvalidAltitude();
 206
 207   Handle(HYDROData_Stream) aStream =
 208     Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
 209   if ( aStream.IsNull() )
 210     return aResAltitude;
 211
 212   TopoDS_Shape aStreamShape = aStream->GetTopShape();
 213   if ( aStreamShape.IsNull() )
 214     return aResAltitude;
 215
 216   TopExp_Explorer aStreamFaceExp( aStreamShape, TopAbs_FACE );
 217   if ( !aStreamFaceExp.More() )
 218     return aResAltitude;
 219
 220   // Get only face because of 2d profile wires is in compound
 221   TopoDS_Face aStreamFace = TopoDS::Face( aStreamFaceExp.Current() );
 222
 223   // Check if point is inside of stream presentation
 224   BRepTopAdaptor_FClass2d aClassifier( aStreamFace, Precision::Confusion() );
 225   TopAbs_State aPointState = aClassifier.Perform( gp_Pnt2d( thePoint ), Standard_False );
 226   if ( aPointState == TopAbs_OUT )
 227     return aResAltitude;
 228
 229   // Find the two profiles between which the point is lies
 230   Handle(HYDROData_Profile) aLeftProfile, aRightProfile;
 231   if ( !getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile ) )
 232     return aResAltitude;
 233
 234   // Find the projections of point to borders of stream
 235   gp_XYZ aPointToTest( thePoint.X(), thePoint.Y(), 0.0 );
 236
 237   TopoDS_Edge aStreamLeftEdge = TopoDS::Edge( aStream->GetLeftShape() );
 238   TopoDS_Edge aStreamRightEdge = TopoDS::Edge( aStream->GetRightShape() );
 239
 240   Standard_Real aFirst = 0.0, aLast = 0.0;
 241   Handle(Geom_Curve) anEdgeLeftCurve = BRep_Tool::Curve( aStreamLeftEdge, aFirst, aLast );
 242   Handle(Geom_Curve) anEdgeRightCurve = BRep_Tool::Curve( aStreamRightEdge, aFirst, aLast );
 243
 244   GeomAPI_ProjectPointOnCurve aLeftProject( aPointToTest, anEdgeLeftCurve );
 245   GeomAPI_ProjectPointOnCurve aRightProject( aPointToTest, anEdgeRightCurve );
 246
 247   Standard_Real aLeftDist = aLeftProject.LowerDistance();
 248   Standard_Real aRightDist = aRightProject.LowerDistance();
 249
 250   // Find the altitude in profiles
 251   Standard_Real aLeftAlt = getAltitudeFromProfile( aLeftProfile, aLeftDist, aRightDist );
 252   Standard_Real aRightAlt = getAltitudeFromProfile( aRightProfile, aLeftDist, aRightDist );
 253
 254   // Interpolate altitudes
 255   Standard_Real aFirstCoeff = aLeftDist / ( aLeftDist + aRightDist );
 256   Standard_Real aSecCoeff = aRightDist / ( aLeftDist + aRightDist );
 257
 258   aResAltitude = aLeftAlt * aFirstCoeff + aRightAlt * aSecCoeff;
 259
 260   return aResAltitude;
 261 }
 262
 263
 264
 265