src/HYDROData/HYDROData_StreamAltitude.cxx

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