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