src/HYDROData/HYDROData_StreamAltitude.cxx

   1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "HYDROData_StreamAltitude.h"
  24
  25 #include "HYDROData_Document.h"
  26 #include "HYDROData_Profile.h"
  27 #include "HYDROData_Stream.h"
  28 #include "HYDROData_ShapesTool.h"
  29 #include <HYDROData_Tool.h>
  30 #include <BRep_Tool.hxx>
  31
  32 #include <BRepBuilderAPI_MakeEdge.hxx>
  33 #include <BRepBuilderAPI_MakeFace.hxx>
  34 #include <BRepBuilderAPI_MakeWire.hxx>
  35
  36 #include <Extrema_ExtElC.hxx>
  37
  38 #include <GeomAPI_ProjectPointOnCurve.hxx>
  39
  40 #include <gp_Lin.hxx>
  41
  42 #include <Precision.hxx>
  43
  44 #include <TopExp_Explorer.hxx>
  45
  46 #include <TopoDS.hxx>
  47 #include <TopoDS_Wire.hxx>
  48 #include <TopoDS_Edge.hxx>
  49
  50 #include <TopTools_SequenceOfShape.hxx>
  51
  52 #include <Geom_Line.hxx>
  53
  54 #include <QStringList>
  55
  56 #ifdef DEB_CLASS2D
  57 #include <BRepTools.hxx>
  58 #include <BRep_Builder.hxx>
  59 #include <BRepBuilderAPI_MakeVertex.hxx>
  60 #endif
  61 IMPLEMENT_STANDARD_HANDLE(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
  62 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_StreamAltitude, HYDROData_IAltitudeObject)
  63
  64 HYDROData_StreamAltitude::HYDROData_StreamAltitude()
  65 : HYDROData_IAltitudeObject()
  66 {
  67 }
  68
  69 HYDROData_StreamAltitude::~HYDROData_StreamAltitude()
  70 {
  71 }
  72
  73 Standard_Real getAltitudeFromProfile( const Handle(HYDROData_Profile)& theProfile,
  74                                       const Standard_Real&             theLeftDist,
  75                                       const Standard_Real&             theRightDist )
  76 {
  77   Standard_Real aResAlt = 0.0;
  78
  79   gp_XY aFirstPoint, aLastPoint;
  80   if ( !theProfile->GetLeftPoint( aFirstPoint, false ) ||
  81        !theProfile->GetRightPoint( aLastPoint, false ) )
  82     return aResAlt;
  83
  84   gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
  85   gp_Pnt aPnt2( aLastPoint.X(),  aLastPoint.Y(),  0 );
  86
  87   Standard_Real aProfileDist = aPnt1.Distance( aPnt2 );
  88
  89   Standard_Real aCoeff = aProfileDist / ( theLeftDist + theRightDist );
  90
  91   gp_Pnt anIntPoint( aPnt1.XYZ() + ( aCoeff * theLeftDist ) * gp_Dir( gp_Vec( aPnt1, aPnt2 ) ).XYZ() );
  92
  93   gp_Lin aPointLine( anIntPoint, gp::DZ() );
  94
  95   gp_Pnt aPrevPoint;
  96   gp_Lin aPrevNormal;
  97   HYDROData_Profile::ProfilePoints aProfilePoints = theProfile->GetProfilePoints( false );
  98   for ( int i = 1, n = aProfilePoints.Length(); i <= n; ++i )
  99   {
 100     gp_Pnt aProfPoint( aProfilePoints.Value( i ) );
 101
 102     Standard_Real aDist = aPointLine.Distance( aProfPoint );
 103     if ( aDist <= gp::Resolution() )
 104     {
 105       // We found the intersected point
 106       aResAlt = aProfPoint.Z();
 107       break;
 108     }
 109
 110     gp_Lin aNormal = aPointLine.Normal( aProfPoint );
 111     if ( i == 1 )
 112     {
 113       aPrevNormal = aNormal;
 114       aPrevPoint = aProfPoint;
 115       continue;
 116     }
 117
 118     if ( aPrevNormal.Direction().Dot( aNormal.Direction() ) < 0 )
 119     {
 120       // We found the intersected edge
 121       gp_Lin anEdgeLine( aPrevPoint, gp_Dir( gp_Vec( aPrevPoint, aProfPoint ) ) );
 122
 123       Extrema_ExtElC anExtrema( aPointLine, anEdgeLine, Precision::Angular() );
 124       if ( !anExtrema.IsParallel() )
 125       {
 126         Extrema_POnCurv aFirstPnt, aSecPnt;
 127         anExtrema.Points( 1, aFirstPnt, aSecPnt );
 128
 129         const gp_Pnt& anIntPnt = aSecPnt.Value();
 130         aResAlt = anIntPnt.Z();
 131
 132         break;
 133       }
 134     }
 135
 136     aPrevNormal = aNormal;
 137     aPrevPoint = aProfPoint;
 138   }
 139
 140   return aResAlt;
 141 }
 142
 143 bool HYDROData_StreamAltitude::getBoundaryProfilesForPoint(
 144   const gp_XY&               thePoint,
 145   Handle(HYDROData_Profile)& theLeftProfile,
 146   Handle(HYDROData_Profile)& theRightProfile ) const
 147 {
 148   Handle(HYDROData_Stream) aStream =
 149     Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
 150   if ( aStream.IsNull() )
 151     return false;
 152
 153   HYDROData_SequenceOfObjects aStreamProfiles = aStream->GetProfiles();
 154   if ( aStreamProfiles.Length() < 2 )
 155     return false;
 156
 157   Handle(HYDROData_Profile) aPrevProfile;
 158   gp_Pnt aPrevPnt1, aPrevPnt2;
 159   for ( int i = 1, n = aStreamProfiles.Length(); i <= n; ++i )
 160   {
 161     Handle(HYDROData_Profile) aProfile =
 162       Handle(HYDROData_Profile)::DownCast( aStreamProfiles.Value( i ) );
 163     if ( aProfile.IsNull() )
 164       continue;
 165
 166     gp_XY aFirstPoint, aLastPoint;
 167     if ( !aProfile->GetLeftPoint( aFirstPoint, false ) ||
 168          !aProfile->GetRightPoint( aLastPoint, false ) )
 169       continue;
 170
 171     gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
 172     gp_Pnt aPnt2( aLastPoint.X(),  aLastPoint.Y(),  0 );
 173
 174     if ( !aPrevProfile.IsNull() )
 175     {
 176       BRepBuilderAPI_MakeEdge aLeftMakeEdge( aPrevPnt1, aPrevPnt2 );
 177       BRepBuilderAPI_MakeEdge aBotMakeEdge( aPrevPnt2, aPnt2 );
 178       BRepBuilderAPI_MakeEdge aRightMakeEdge( aPnt2, aPnt1 );
 179       BRepBuilderAPI_MakeEdge aTopMakeEdge( aPnt1, aPrevPnt1 );
 180
 181       BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(),
 182                                          aRightMakeEdge.Edge(), aTopMakeEdge.Edge() );
 183
 184       BRepBuilderAPI_MakeFace aMakeFace( aMakeWire.Wire() );
 185
 186       TopoDS_Face aProfilesFace = aMakeFace.Face();
 187 #ifdef DEB_CLASS2D
 188           TopoDS_Compound aCmp;
 189       BRep_Builder aBB;
 190       aBB.MakeCompound(aCmp);
 191           aBB.Add(aCmp, aProfilesFace);
 192           gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
 193           BRepBuilderAPI_MakeVertex aMk(aPnt);
 194           aBB.Add(aCmp, aMk.Vertex());
 195           BRepTools::Write(aCmp, "ProfileFace.brep");
 196 #endif
 197
 198           TopAbs_State aPointState =  HYDROData_Tool::ComputePointState(thePoint, aProfilesFace);
 199
 200 #ifdef DEB_CLASS2D
 201           cout << "Point status is = " << aPointState <<endl;
 202 #endif
 203       if ( aPointState != TopAbs_OUT )
 204       {
 205         theLeftProfile = aPrevProfile;
 206         theRightProfile = aProfile;
 207         break;
 208       }
 209     }
 210
 211     aPrevProfile = aProfile;
 212     aPrevPnt1 = aPnt1;
 213     aPrevPnt2 = aPnt2;
 214   }
 215
 216   return !theLeftProfile.IsNull() && !theRightProfile.IsNull();
 217 }
 218
 219 double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) const
 220 {
 221   double aResAltitude = GetInvalidAltitude();
 222
 223   Handle(HYDROData_Stream) aStream =
 224     Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
 225   if ( aStream.IsNull() )
 226     return aResAltitude;
 227
 228   TopoDS_Shape aStreamShape = aStream->GetTopShape();
 229   if ( aStreamShape.IsNull() )
 230     return aResAltitude;
 231
 232   TopExp_Explorer aStreamFaceExp( aStreamShape, TopAbs_FACE );
 233   if ( !aStreamFaceExp.More() )
 234     return aResAltitude;
 235
 236   // Get only face because of 2d profile wires is in compound
 237   TopoDS_Face aStreamFace = TopoDS::Face( aStreamFaceExp.Current() );
 238
 239   // Check if point is inside of stream presentation
 240   TopAbs_State aPointState = HYDROData_Tool::ComputePointState(thePoint, aStreamFace);
 241
 242 #ifdef DEB_CLASS2D
 243       cout << "Point status is = " << aPointState <<endl;
 244           TopoDS_Compound aCmp;
 245       BRep_Builder aBB;
 246       aBB.MakeCompound(aCmp);
 247           aBB.Add(aCmp, aStreamFace);
 248           gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
 249           BRepBuilderAPI_MakeVertex aMk(aPnt);
 250           aBB.Add(aCmp, aMk.Vertex());
 251           BRepTools::Write(aCmp, "FCL2d.brep");
 252 #endif
 253   if ( aPointState == TopAbs_OUT )
 254     return aResAltitude;
 255
 256   // Find the two profiles between which the point is lies
 257   Handle(HYDROData_Profile) aLeftProfile, aRightProfile;
 258   if ( !getBoundaryProfilesForPoint( thePoint, aLeftProfile, aRightProfile ) )
 259     return aResAltitude;
 260
 261   // Find the projections of point to borders of stream
 262   gp_XYZ aPointToTest( thePoint.X(), thePoint.Y(), 0.0 );
 263
 264   TopoDS_Edge aStreamLeftEdge = TopoDS::Edge( aStream->GetLeftShape() );
 265   TopoDS_Edge aStreamRightEdge = TopoDS::Edge( aStream->GetRightShape() );
 266
 267   Standard_Real aFirst = 0.0, aLast = 0.0;
 268   Handle(Geom_Curve) anEdgeLeftCurve = BRep_Tool::Curve( aStreamLeftEdge, aFirst, aLast );
 269   Handle(Geom_Curve) anEdgeRightCurve = BRep_Tool::Curve( aStreamRightEdge, aFirst, aLast );
 270
 271   GeomAPI_ProjectPointOnCurve aLeftProject( aPointToTest, anEdgeLeftCurve );
 272   GeomAPI_ProjectPointOnCurve aRightProject( aPointToTest, anEdgeRightCurve );
 273
 274   Standard_Real aLeftDist = aLeftProject.LowerDistance();
 275   Standard_Real aRightDist = aRightProject.LowerDistance();
 276
 277   // Find the altitude in profiles
 278   Standard_Real aLeftAlt = getAltitudeFromProfile( aLeftProfile, aLeftDist, aRightDist );
 279   Standard_Real aRightAlt = getAltitudeFromProfile( aRightProfile, aLeftDist, aRightDist );
 280
 281   // Interpolate altitudes
 282   // Left profile line ( the segment between the firts and the last profile point )
 283   HYDROData_Profile::ProfilePoints aLeftProfilePoints = aLeftProfile->GetProfilePoints( false );
 284   gp_Pnt aLeftProfileP1( aLeftProfilePoints.First() );
 285   aLeftProfileP1.SetZ( 0 );
 286   gp_Pnt aLeftProfileP2( aLeftProfilePoints.Last() );
 287   aLeftProfileP2.SetZ( 0 );
 288   gp_Vec aLeftProfileVec( aLeftProfileP1, aLeftProfileP2 );
 289   Handle(Geom_Line) aLeftProfileLine = new Geom_Line( gp_Ax1( aLeftProfileP1, aLeftProfileVec ) );
 290   // Right profile line
 291   HYDROData_Profile::ProfilePoints aRightProfilePoints = aRightProfile->GetProfilePoints( false );
 292   gp_Pnt aRightProfileP1( aRightProfilePoints.First() );
 293   aRightProfileP1.SetZ( 0 );
 294   gp_Pnt aRightProfileP2( aRightProfilePoints.Last() );
 295   aRightProfileP2.SetZ( 0 );
 296   gp_Vec aRightProfileVec( aRightProfileP1, aRightProfileP2 );
 297   Handle(Geom_Line) aRightProfileLine = new Geom_Line( gp_Ax1( aRightProfileP1, aRightProfileVec ) );
 298   // The point projections on the left and right profiles
 299   GeomAPI_ProjectPointOnCurve aLeftProfileProject( aPointToTest, aLeftProfileLine );
 300   GeomAPI_ProjectPointOnCurve aRightProfileProject( aPointToTest, aRightProfileLine );
 301   // The point distance to the left and right profiles
 302   Standard_Real aLeftProfileDist = aLeftProfileProject.LowerDistance();
 303   Standard_Real aRightProfileDist = aRightProfileProject.LowerDistance();
 304   // The coefficients
 305   Standard_Real aFirstCoeff = aLeftProfileDist / ( aLeftProfileDist + aRightProfileDist );
 306   Standard_Real aSecCoeff = aRightProfileDist / ( aLeftProfileDist + aRightProfileDist );
 307
 308   aResAltitude = aLeftAlt * ( 1 - aFirstCoeff ) + aRightAlt * ( 1 - aSecCoeff );
 309
 310   return aResAltitude;
 311 }
 312
 313
 314
 315