src/HYDROData/HYDROData_ChannelAltitude.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_ChannelAltitude.h"
  20
  21 #include "HYDROData_Document.h"
  22 #include "HYDROData_Object.h"
  23 #include "HYDROData_Channel.h"
  24 #include "HYDROData_Projection.h"
  25 #include "HYDROData_Polyline3D.h"
  26 #include "HYDROData_PolylineXY.h"
  27 #include "HYDROData_ProfileUZ.h"
  28 #include "HYDROData_Profile.h"
  29 #include <TDataStd_Integer.hxx>
  30
  31 //#define _DEVDEBUG_
  32 #include "HYDRO_trace.hxx"
  33 #include <QString>
  34
  35 #include <BRepBuilderAPI_MakeVertex.hxx>
  36 #include <BRepExtrema_DistShapeShape.hxx>
  37
  38 #include <Precision.hxx>
  39 #include <TopAbs.hxx>
  40 #include <TopExp.hxx>
  41 #include <TopoDS.hxx>
  42 #include <TopoDS_Shape.hxx>
  43 #include <TopoDS_Edge.hxx>
  44 #include <TopoDS_Face.hxx>
  45 #include <TopoDS_Vertex.hxx>
  46 #include <TopoDS_Wire.hxx>
  47 #include <TopExp_Explorer.hxx>
  48 #include <TopAbs_ShapeEnum.hxx>
  49 #include <BRep_Tool.hxx>
  50 #include <BRepTools.hxx>
  51 #include <Geom_Curve.hxx>
  52
  53 #include <gp_Trsf.hxx>
  54 #include <gp_Pnt.hxx>
  55 #include <gp_Vec.hxx>
  56
  57 IMPLEMENT_STANDARD_RTTIEXT(HYDROData_ChannelAltitude, HYDROData_IAltitudeObject)
  58
  59 HYDROData_ChannelAltitude::HYDROData_ChannelAltitude()
  60 : HYDROData_IAltitudeObject()
  61 {
  62 }
  63
  64 HYDROData_ChannelAltitude::~HYDROData_ChannelAltitude()
  65 {
  66 }
  67
  68 double HYDROData_ChannelAltitude::GetAltitudeForPoint( const gp_XY& thePoint,
  69                                                        int theMethod) const
  70 {
  71   DEBTRACE("GetAltitudeForPoint p(" << thePoint.X() << ", " << thePoint.Y() << ")");
  72   double aResAltitude = GetInvalidAltitude();
  73
  74   Handle(HYDROData_Channel) aChannel =
  75   Handle(HYDROData_Channel)::DownCast(GetFatherObject());
  76   if (aChannel.IsNull())
  77     {
  78       DEBTRACE("aChannel.IsNull()");
  79       return aResAltitude;
  80     }
  81   DEBTRACE("aChannel: " << aChannel->GetName().toStdString());
  82
  83   Handle(HYDROData_Polyline3D) aGuideLine = aChannel->GetGuideLine();
  84   if (aGuideLine.IsNull())
  85     {
  86       DEBTRACE("aGuideLine.IsNull()");
  87       return aResAltitude;
  88     }
  89   //DEBTRACE("aGuideLine: " << aGuideLine->GetName().toStdString());
  90
  91   Handle(HYDROData_PolylineXY) aGuideXY = aGuideLine->GetPolylineXY();
  92   if (aGuideXY.IsNull())
  93     {
  94       DEBTRACE("aGuideXY.IsNull()");
  95       return aResAltitude;
  96     }
  97   //DEBTRACE("aGuideXY: " << aGuideXY->GetName().toStdString());
  98
  99   Handle(HYDROData_ProfileUZ) aGuideUZ = aGuideLine->GetProfileUZ();
 100   if (aGuideUZ.IsNull())
 101     {
 102       aGuideUZ = aGuideLine->GetChildProfileUZ(); // profile obtained from bathymetry
 103     }
 104   if (aGuideUZ.IsNull())
 105     {
 106       DEBTRACE("aGuideUZ.IsNull()");
 107       return aResAltitude;
 108     }
 109   //DEBTRACE("aGuideUZ: " << aGuideUZ->GetName().toStdString());
 110
 111   Handle (HYDROData_Profile) aProfile = aChannel->GetProfile();
 112   if (aProfile.IsNull())
 113     {
 114       return aResAltitude;
 115     }
 116   //DEBTRACE("aProfile: " << aProfile->GetName().toStdString());
 117
 118   // --- See GEOMImpl_ProjectionDriver.cxx
 119
 120   gp_XY LP, RP;
 121   aProfile->GetLeftPoint(LP);
 122   aProfile->GetRightPoint(RP);
 123
 124   TopoDS_Shape aShape =  aGuideXY->GetShape();
 125   gp_Pnt P1(thePoint.X(), thePoint.Y(), 0);
 126   TopoDS_Shape aPoint = BRepBuilderAPI_MakeVertex(P1).Shape();
 127
 128   if (aPoint.IsNull() || aShape.IsNull())
 129     {
 130       DEBTRACE("aPoint.IsNull() || aShape.IsNull()");
 131       return aResAltitude;
 132     }
 133
 134   if (aShape.ShapeType() != TopAbs_EDGE && aShape.ShapeType() != TopAbs_WIRE)
 135     {
 136       DEBTRACE("Projection aborted : the shape is neither an edge nor a wire");
 137       return aResAltitude;
 138     }
 139
 140   // Perform projection.
 141   BRepExtrema_DistShapeShape aDistShSh(aPoint, aShape, Extrema_ExtFlag_MIN);
 142
 143   if (aDistShSh.IsDone() == Standard_False)
 144     {
 145       DEBTRACE("Projection not done");
 146       return aResAltitude;
 147     }
 148
 149   Standard_Boolean hasValidSolution = Standard_False;
 150   Standard_Integer aNbSolutions = aDistShSh.NbSolution();
 151   Standard_Integer i;
 152   double aParam = 0.;
 153   Standard_Real aTolConf = BRep_Tool::Tolerance(TopoDS::Vertex(aPoint));
 154   Standard_Real aTolAng = 1.e-4;
 155
 156   for (i = 1; i <= aNbSolutions; i++)
 157     {
 158       Standard_Boolean isValid = Standard_False;
 159       BRepExtrema_SupportType aSupportType = aDistShSh.SupportTypeShape2(i);
 160       TopoDS_Shape aSupportShape = aDistShSh.SupportOnShape2(i);
 161
 162       if (aSupportType == BRepExtrema_IsOnEdge)
 163         {
 164           // Minimal distance inside edge is really a projection.
 165           isValid = Standard_True;
 166           aDistShSh.ParOnEdgeS2(i, aParam);
 167         }
 168       else if (aSupportType == BRepExtrema_IsVertex)
 169         {
 170           TopExp_Explorer anExp(aShape, TopAbs_EDGE);
 171
 172           if (aDistShSh.Value() <= aTolConf)
 173             {
 174               // The point lies on the shape. This means this point
 175               // is really a projection.
 176               for (; anExp.More() && !isValid; anExp.Next())
 177                 {
 178                   TopoDS_Edge aCurEdge = TopoDS::Edge(anExp.Current());
 179
 180                   if (aCurEdge.IsNull() == Standard_False)
 181                     {
 182                       TopoDS_Vertex aVtx[2];
 183
 184                       TopExp::Vertices(aCurEdge, aVtx[0], aVtx[1]);
 185
 186                       for (int j = 0; j < 2; j++)
 187                         {
 188                           if (aSupportShape.IsSame(aVtx[j]))
 189                             {
 190                               // The current edge is a projection edge.
 191                               isValid = Standard_True;
 192                               aSupportShape = aCurEdge;
 193                               aParam = BRep_Tool::Parameter(aVtx[j], aCurEdge);
 194                               break;
 195                             }
 196                         }
 197                     }
 198                 }
 199             }
 200           else
 201             {
 202               // Minimal distance to vertex is not always a real projection.
 203               gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(aPoint));
 204               gp_Pnt aPrjPnt = BRep_Tool::Pnt(TopoDS::Vertex(aSupportShape));
 205               gp_Vec aDProjP(aPrjPnt, aPnt);
 206
 207               for (; anExp.More() && !isValid; anExp.Next())
 208                 {
 209                   TopoDS_Edge aCurEdge = TopoDS::Edge(anExp.Current());
 210
 211                   if (aCurEdge.IsNull() == Standard_False)
 212                     {
 213                       TopoDS_Vertex aVtx[2];
 214
 215                       TopExp::Vertices(aCurEdge, aVtx[0], aVtx[1]);
 216
 217                       for (int j = 0; j < 2; j++)
 218                         {
 219                           if (aSupportShape.IsSame(aVtx[j]))
 220                             {
 221                               // Check if the point is a projection to the current edge.
 222                               Standard_Real anEdgePars[2];
 223                               Handle(Geom_Curve) aCurve = BRep_Tool::Curve(aCurEdge, anEdgePars[0], anEdgePars[1]);
 224                               gp_Pnt aVal;
 225                               gp_Vec aD1;
 226
 227                               aParam = BRep_Tool::Parameter(aVtx[j], aCurEdge);
 228                               aCurve->D1(aParam, aVal, aD1);
 229
 230                               if (Abs(aD1.Dot(aDProjP)) <= aTolAng)
 231                                 {
 232                                   // The current edge is a projection edge.
 233                                   isValid = Standard_True;
 234                                   aSupportShape = aCurEdge;
 235                                   break;
 236                                 }
 237                             }
 238                         }
 239                     }
 240                 }
 241             }
 242         }
 243
 244       if (isValid)
 245         {
 246           if (hasValidSolution)
 247             {
 248               DEBTRACE("Projection aborted : multiple solutions");
 249               return aResAltitude;
 250             }
 251
 252           // Store the valid solution.
 253           hasValidSolution = Standard_True;
 254
 255           // profile altitude at projection point
 256           HYDROData_IPolyline::PointsList aProfilePoints = aGuideUZ->GetPoints();
 257           if ( aProfilePoints.IsEmpty() )
 258             {
 259               DEBTRACE("empty profile UZ");
 260               return aResAltitude;
 261             }
 262           double aDepth = HYDROData_ProfileUZ::GetDepthFromDistance( aProfilePoints, aParam );
 263           DEBTRACE("profile altitude: " << aDepth);
 264
 265           // Compute edge index.
 266           TopExp_Explorer anExp(aShape, TopAbs_EDGE);
 267           int anIndex = 0;
 268           for (; anExp.More(); anExp.Next(), anIndex++)
 269             {
 270               if (aSupportShape.IsSame(anExp.Current()))
 271                 {
 272                   break;
 273                 }
 274             }
 275
 276           // get the XY distance from point to guideline
 277           const gp_Pnt &aPntProj = aDistShSh.PointOnShape2(i);
 278           DEBTRACE("projection: (" << aPntProj.X() << ", " << aPntProj.Y() << ", " << aPntProj.Z() << ")");
 279           gp_XY aProjXY = gp_XY(aPntProj.X(), aPntProj.Y());
 280           aProjXY.Subtract(thePoint);
 281           double distance = aProjXY.Modulus();
 282
 283           gp_Vec2d aProjDir = aProjXY;
 284           gp_Vec2d aProfileDir(LP, RP);
 285           double aSign = 1.0;
 286           if( aProfileDir.Dot(aProjDir) < 0 )
 287             aSign = -1.0;
 288
 289           distance *= aSign;
 290           DEBTRACE("distance to guideline " << distance);
 291
 292           if (GetInvertDirection())
 293             distance*=-1;
 294
 295           // get delta altitude on section (supposed symmetric) from guideline distance (aParam)
 296           double delta = 0;
 297           int i1 = 0;
 298           int i2 = 0;
 299           gp_XY pt1 = gp_XY();
 300           gp_XY pt2 = gp_XY();
 301           HYDROData_ProfileUZ::PointsList aSectionPoints = aProfile->GetParametricPoints();
 302           for ( int i = 1, aNbPoints = aSectionPoints.Size(); i <= aNbPoints; ++i )
 303             {
 304               const HYDROData_IPolyline::Point& aPolylinePoint = aSectionPoints.Value( i );
 305               DEBTRACE(" i, size, profile point: " << i <<  " " << aSectionPoints.Size() << " " << aPolylinePoint.X() << " " << aPolylinePoint.Y());
 306               if (aPolylinePoint.X() < distance)
 307                 {
 308                   i1 = i;
 309                   pt1 = aPolylinePoint;
 310                 }
 311               if (aPolylinePoint.X() >= distance)
 312                 {
 313                   i2 = i;
 314                   pt2 = aPolylinePoint;
 315                   break;
 316                 }
 317             }
 318           if ((i1 == 0) && (i2 == 0))
 319             {
 320               DEBTRACE("Projection aborted : non centered profile");
 321               return aResAltitude;
 322             }
 323           else if (i1 == 0)  // distance <  - profile width
 324             {
 325               DEBTRACE("distance < - profile width");
 326               delta = pt2.Y();
 327             }
 328           else if (i1 == aSectionPoints.Size()) // distance >= profile width
 329             {
 330               DEBTRACE("distance >= profile width");
 331               delta = pt1.Y();
 332             }
 333           else
 334             {
 335               delta = pt1.Y() + (pt2.Y() - pt1.Y())*(distance -pt1.X())/(pt2.X()-pt1.X());
 336             }
 337           aResAltitude = delta + aDepth;
 338           DEBTRACE("distance XY: "<< aParam << " distance to guideline: " << distance << " final altitude: " << aResAltitude << " delta: " << delta);
 339           return aResAltitude;
 340         }
 341     }
 342
 343   if (!hasValidSolution)
 344     {
 345       DEBTRACE("Projection aborted : no projection");
 346       return aResAltitude;
 347     }
 348
 349   return aResAltitude;
 350 }
 351
 352
 353 void HYDROData_ChannelAltitude::SetInvertDirection(bool IsInverted)
 354 {
 355   if ( GetInvertDirection() == IsInverted )
 356     return;
 357
 358   TDataStd_Integer::Set( myLab.FindChild( DataTag_InvertDirection ), (Standard_Integer)IsInverted );
 359 }
 360
 361 bool HYDROData_ChannelAltitude::GetInvertDirection() const
 362 {
 363   bool aRes = false;
 364
 365   TDF_Label aLabel = myLab.FindChild( DataTag_InvertDirection, false );
 366   if ( !aLabel.IsNull() )
 367   {
 368     Handle(TDataStd_Integer) anIntVal;
 369     if ( aLabel.FindAttribute( TDataStd_Integer::GetID(), anIntVal ) )
 370       aRes = (bool)anIntVal->Get();
 371   }
 372
 373   return aRes;
 374 }
 375
 376