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