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