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_Channel.h"
  22 #include "HYDROData_Document.h"
  23 #include "HYDROData_Profile.h"
  24 #include "HYDROData_Stream.h"
  25 #include "HYDROData_ShapesTool.h"
  26 #include <HYDROData_Tool.h>
  27 #include <BRep_Tool.hxx>
  28
  29 #include <BRepBuilderAPI_MakeEdge.hxx>
  30 #include <BRepBuilderAPI_MakeFace.hxx>
  31 #include <BRepBuilderAPI_MakeWire.hxx>
  32
  33 #include <Extrema_ExtElC.hxx>
  34
  35 #include <GeomAPI_ProjectPointOnCurve.hxx>
  36
  37 #include <gp_Lin.hxx>
  38
  39 #include <Precision.hxx>
  40
  41 #include <TopExp.hxx>
  42 #include <TopExp_Explorer.hxx>
  43
  44 #include <TopoDS.hxx>
  45 #include <TopoDS_Wire.hxx>
  46 #include <TopoDS_Edge.hxx>
  47
  48 #include <TopTools_SequenceOfShape.hxx>
  49 #include <TopTools_ListIteratorOfListOfShape.hxx>
  50
  51 #include <Geom_Line.hxx>
  52
  53 #include <QStringList>
  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 bool IsPointBetweenEdges( const gp_Pnt& aFirstPnt1, const gp_Pnt& aLastPnt1,
  73                           const gp_Pnt& aFirstPnt2, const gp_Pnt& aLastPnt2,
  74                           const gp_Pnt& thePoint) {
  75   BRepBuilderAPI_MakeEdge aLeftMakeEdge( aFirstPnt1, aLastPnt1 );
  76   BRepBuilderAPI_MakeEdge aBotMakeEdge( aLastPnt1, aLastPnt2 );
  77   BRepBuilderAPI_MakeEdge aRightMakeEdge( aLastPnt2, aFirstPnt2 );
  78   BRepBuilderAPI_MakeEdge aTopMakeEdge( aFirstPnt2, aFirstPnt1 );
  79
  80   BRepBuilderAPI_MakeWire aMakeWire( aLeftMakeEdge.Edge(), aBotMakeEdge.Edge(),
  81                                      aRightMakeEdge.Edge(), aTopMakeEdge.Edge() );
  82
  83   BRepBuilderAPI_MakeFace aMakeFace( aMakeWire.Wire() );
  84
  85   TopoDS_Face aFace = aMakeFace.Face();
  86 #ifdef DEB_CLASS2D
  87   TopoDS_Compound aCmp;
  88   BRep_Builder aBB;
  89   aBB.MakeCompound(aCmp);
  90   aBB.Add(aCmp, aFace);
  91   BRepBuilderAPI_MakeVertex aMk(thePoint);
  92   aBB.Add(aCmp, aMk.Vertex());
  93   BRepTools::Write(aCmp, "ProfileFace.brep");
  94 #endif
  95
  96   gp_XY anXY( thePoint.X(), thePoint.Y() );
  97   TopAbs_State aPointState =  HYDROData_Tool::ComputePointState(anXY, aFace);
  98
  99 #ifdef DEB_CLASS2D
 100           cout << "Point status is = " << aPointState <<endl;
 101 #endif
 102   return aPointState != TopAbs_OUT;
 103 }
 104
 105 Standard_Real getAltitudeFromWire( const TopoDS_Wire& theWire,
 106                                    const Standard_Real& theLeftDist,
 107                                    const Standard_Real& theRightDist )
 108 {
 109   Standard_Real aResAlt = 0.0;
 110
 111   gp_XY aFirstPoint, aLastPoint;
 112   TopoDS_Vertex aFirstVertex, aLastVertex;
 113   TopExp::Vertices( theWire, aFirstVertex, aLastVertex );
 114
 115   gp_Pnt aPnt1( BRep_Tool::Pnt( aFirstVertex ) );
 116   aPnt1.SetZ( 0 );
 117   gp_Pnt aPnt2( BRep_Tool::Pnt( aLastVertex ) );
 118   aPnt2.SetZ( 0 );
 119
 120   Standard_Real aProfileDist = aPnt1.Distance( aPnt2 );
 121
 122   Standard_Real aCoeff = aProfileDist / ( theLeftDist + theRightDist );
 123
 124   gp_Pnt anIntPoint( aPnt1.XYZ() + ( aCoeff * theLeftDist ) * gp_Dir( gp_Vec( aPnt1, aPnt2 ) ).XYZ() );
 125
 126   return HYDROData_Tool::GetAltitudeForWire( theWire,
 127                                              gp_XY(anIntPoint.X(), anIntPoint.Y()),
 128                                              1E-2,
 129                                              1E-2,
 130                                              HYDROData_IAltitudeObject::GetInvalidAltitude() );
 131 }
 132
 133 bool HYDROData_StreamAltitude::getBoundaryWiresForPoint(
 134   const gp_XY& thePoint,
 135   TopoDS_Wire& theLeftWire,
 136   TopoDS_Wire& theRightWire ) const
 137 {
 138   gp_Pnt aTestPnt( thePoint.X(),  thePoint.Y(), 0 );
 139
 140   Handle(HYDROData_Object) anObject =
 141     Handle(HYDROData_Object)::DownCast( GetFatherObject() );
 142   if ( anObject.IsNull() ) {
 143     return false;
 144   }
 145
 146   if ( anObject->GetKind() == KIND_STREAM ) {
 147     Handle(HYDROData_Stream) aStream = Handle(HYDROData_Stream)::DownCast( anObject );
 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, false ) ||
 166            !aProfile->GetRightPoint( aLastPoint, false ) )
 167         continue;
 168
 169       gp_Pnt aPnt1( aFirstPoint.X(), aFirstPoint.Y(), 0 );
 170       gp_Pnt aPnt2( aLastPoint.X(),  aLastPoint.Y(),  0 );
 171
 172       if ( !aPrevProfile.IsNull() )
 173       {
 174         if ( IsPointBetweenEdges( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) )
 175         {
 176           theLeftWire = TopoDS::Wire( aPrevProfile->GetShape3D() );
 177           theRightWire = TopoDS::Wire( aProfile->GetShape3D() );
 178           break;
 179         }
 180       }
 181
 182       aPrevProfile = aProfile;
 183       aPrevPnt1 = aPnt1;
 184       aPrevPnt2 = aPnt2;
 185     }
 186   } else if ( anObject->GetKind() == KIND_CHANNEL ) {
 187     Handle(HYDROData_Channel) aChannel = Handle(HYDROData_Channel)::DownCast( anObject );
 188     if ( aChannel.IsNull() )
 189       return false;
 190
 191     TopTools_ListOfShape aWiresList;
 192     TopExp_Explorer anExp( aChannel->GetShape3D(), TopAbs_WIRE );
 193     for ( ; anExp.More(); anExp.Next() ) {
 194       if(!anExp.Current().IsNull()) {
 195         const TopoDS_Wire& aWire = TopoDS::Wire( anExp.Current() );
 196         aWiresList.Append( aWire );
 197       }
 198     }
 199
 200     if ( aWiresList.Extent() < 2 )
 201       return false;
 202
 203     TopoDS_Wire aPrevWire;
 204     gp_Pnt aPrevPnt1, aPrevPnt2;
 205
 206     TopTools_ListIteratorOfListOfShape anIt( aWiresList );
 207     for ( ; anIt.More(); anIt.Next() ) {
 208       TopoDS_Wire& aWire = TopoDS::Wire( anIt.Value() );
 209       if ( aWire.IsNull() )
 210         continue;
 211
 212       TopoDS_Vertex aFirstVertex, aLastVertex;
 213       TopExp::Vertices( aWire, aFirstVertex, aLastVertex );
 214
 215       gp_Pnt aPnt1( BRep_Tool::Pnt( aFirstVertex ) );
 216       aPnt1.SetZ( 0 );
 217       gp_Pnt aPnt2( BRep_Tool::Pnt( aLastVertex ) );
 218       aPnt2.SetZ( 0 );
 219
 220       if ( !aPrevWire.IsNull() ) {
 221         if ( IsPointBetweenEdges( aPrevPnt1, aPrevPnt2, aPnt1, aPnt2, aTestPnt ) ) {
 222           theLeftWire = aPrevWire;
 223           theRightWire = aWire;
 224           break;
 225         }
 226       }
 227
 228       aPrevWire = aWire;
 229       aPrevPnt1 = aPnt1;
 230       aPrevPnt2 = aPnt2;
 231     }
 232   }
 233
 234   return !theLeftWire.IsNull() && !theRightWire.IsNull();
 235 }
 236
 237 double HYDROData_StreamAltitude::GetAltitudeForPoint( const gp_XY& thePoint ) const
 238 {
 239   double aResAltitude = GetInvalidAltitude();
 240
 241   Handle(HYDROData_Object) anObject =
 242     Handle(HYDROData_Object)::DownCast( GetFatherObject() );
 243   if ( anObject.IsNull() )
 244     return aResAltitude;
 245
 246   TopoDS_Shape aTopShape = anObject->GetTopShape();
 247   if ( aTopShape.IsNull() )
 248     return aResAltitude;
 249
 250   TopExp_Explorer aFaceExp( aTopShape, TopAbs_FACE );
 251   if ( !aFaceExp.More() )
 252     return aResAltitude;
 253
 254   // Get only face because of 2d profile wires is in compound
 255   TopoDS_Face aFace = TopoDS::Face( aFaceExp.Current() );
 256
 257   // Check if point is inside of stream/channel presentation
 258   TopAbs_State aPointState = HYDROData_Tool::ComputePointState(thePoint, aFace);
 259
 260 #ifdef DEB_CLASS2D
 261       cout << "Point status is = " << aPointState <<endl;
 262           TopoDS_Compound aCmp;
 263       BRep_Builder aBB;
 264       aBB.MakeCompound(aCmp);
 265           aBB.Add(aCmp, aFace);
 266           gp_Pnt aPnt (thePoint.X(), thePoint.Y(), 0.);
 267           BRepBuilderAPI_MakeVertex aMk(aPnt);
 268           aBB.Add(aCmp, aMk.Vertex());
 269           BRepTools::Write(aCmp, "FCL2d.brep");
 270 #endif
 271   if ( aPointState == TopAbs_OUT )
 272     return aResAltitude;
 273
 274   TopoDS_Edge aLeftEdge, aRightEdge;
 275
 276   if ( anObject->GetKind() == KIND_STREAM ) {
 277     Handle(HYDROData_Stream) aStream =
 278       Handle(HYDROData_Stream)::DownCast( GetFatherObject() );
 279     if ( !aStream.IsNull() ) {
 280       aLeftEdge = TopoDS::Edge( aStream->GetLeftShape() );
 281       aRightEdge = TopoDS::Edge( aStream->GetRightShape() );
 282     }
 283   } else if ( anObject->GetKind() == KIND_CHANNEL ) {
 284     Handle(HYDROData_Channel) aChannel =
 285       Handle(HYDROData_Channel)::DownCast( GetFatherObject() );
 286     if ( !aChannel.IsNull() ) {
 287       aLeftEdge = TopoDS::Edge( aChannel->GetLeftShape() );
 288       aRightEdge = TopoDS::Edge( aChannel->GetRightShape() );
 289     }
 290   }
 291
 292   // Find the two profiles between which the point is lies
 293   TopoDS_Wire aLeftWire, aRightWire;
 294   if ( !getBoundaryWiresForPoint( thePoint, aLeftWire, aRightWire ) )
 295       return aResAltitude;
 296
 297   // Find the projections of point to borders of stream
 298   gp_XYZ aPointToTest( thePoint.X(), thePoint.Y(), 0.0 );
 299
 300   Standard_Real aFirst = 0.0, aLast = 0.0;
 301   Handle(Geom_Curve) anEdgeLeftCurve = BRep_Tool::Curve( aLeftEdge, aFirst, aLast );
 302   Handle(Geom_Curve) anEdgeRightCurve = BRep_Tool::Curve( aRightEdge, aFirst, aLast );
 303
 304   GeomAPI_ProjectPointOnCurve aLeftProject( aPointToTest, anEdgeLeftCurve );
 305   GeomAPI_ProjectPointOnCurve aRightProject( aPointToTest, anEdgeRightCurve );
 306
 307   int aNbPoints1 = aRightProject.NbPoints();
 308   int aNbPoints2 = aLeftProject.NbPoints();
 309   if ( aNbPoints1 < 1 || aNbPoints2 < 1)
 310      return aResAltitude;
 311   Standard_Real aLeftDist = aLeftProject.LowerDistance();
 312   Standard_Real aRightDist = aRightProject.LowerDistance();
 313
 314   // Find the altitude in profiles
 315   Standard_Real aLeftAlt, aRightAlt;
 316   gp_Pnt aLeftProfileP1, aLeftProfileP2, aRightProfileP1, aRightProfileP2;
 317   aLeftAlt = getAltitudeFromWire( aLeftWire, aLeftDist, aRightDist );
 318   aRightAlt = getAltitudeFromWire( aRightWire, aLeftDist, aRightDist );
 319
 320   TopoDS_Vertex aFirstVertex, aLastVertex;
 321   TopExp::Vertices( aLeftWire, aFirstVertex, aLastVertex );
 322   aLeftProfileP1 = BRep_Tool::Pnt( aFirstVertex );
 323   aLeftProfileP2 = BRep_Tool::Pnt( aLastVertex );
 324
 325   TopExp::Vertices( aRightWire, aFirstVertex, aLastVertex );
 326   aRightProfileP1 = BRep_Tool::Pnt( aFirstVertex );
 327   aRightProfileP2 = BRep_Tool::Pnt( aLastVertex );
 328
 329   // Interpolate altitudes
 330   // Left profile line ( the segment between the firts and the last profile point )
 331   aLeftProfileP1.SetZ( 0 );
 332   aLeftProfileP2.SetZ( 0 );
 333   gp_Vec aLeftProfileVec( aLeftProfileP1, aLeftProfileP2 );
 334   Handle(Geom_Line) aLeftProfileLine = new Geom_Line( gp_Ax1( aLeftProfileP1, aLeftProfileVec ) );
 335   // Right profile line
 336   aRightProfileP1.SetZ( 0 );
 337   aRightProfileP2.SetZ( 0 );
 338   gp_Vec aRightProfileVec( aRightProfileP1, aRightProfileP2 );
 339   Handle(Geom_Line) aRightProfileLine = new Geom_Line( gp_Ax1( aRightProfileP1, aRightProfileVec ) );
 340   // The point projections on the left and right profiles
 341   GeomAPI_ProjectPointOnCurve aLeftProfileProject( aPointToTest, aLeftProfileLine );
 342   GeomAPI_ProjectPointOnCurve aRightProfileProject( aPointToTest, aRightProfileLine );
 343   // The point distance to the left and right profiles
 344   Standard_Real aLeftProfileDist = aLeftProfileProject.LowerDistance();
 345   Standard_Real aRightProfileDist = aRightProfileProject.LowerDistance();
 346   // The coefficients
 347   Standard_Real aFirstCoeff = aLeftProfileDist / ( aLeftProfileDist + aRightProfileDist );
 348   Standard_Real aSecCoeff = aRightProfileDist / ( aLeftProfileDist + aRightProfileDist );
 349
 350   aResAltitude = aLeftAlt * ( 1 - aFirstCoeff ) + aRightAlt * ( 1 - aSecCoeff );
 351
 352   return aResAltitude;
 353 }
 354
 355
 356
 357