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
  30 #define _DEVDEBUG_
  31 #include "HYDRO_trace.hxx"
  32 #include <QString>
  33
  34 #include <BRepBuilderAPI_MakeVertex.hxx>
  35 #include <BRepExtrema_DistShapeShape.hxx>
  36
  37 #include <Precision.hxx>
  38 #include <TopAbs.hxx>
  39 #include <TopExp.hxx>
  40 #include <TopoDS.hxx>
  41 #include <TopoDS_Shape.hxx>
  42 #include <TopoDS_Edge.hxx>
  43 #include <TopoDS_Face.hxx>
  44 #include <TopoDS_Vertex.hxx>
  45 #include <TopoDS_Wire.hxx>
  46 #include <TopExp_Explorer.hxx>
  47 #include <TopAbs_ShapeEnum.hxx>
  48 #include <BRep_Tool.hxx>
  49 #include <BRepTools.hxx>
  50 #include <Geom_Curve.hxx>
  51
  52 #include <gp_Trsf.hxx>
  53 #include <gp_Pnt.hxx>
  54 #include <gp_Vec.hxx>
  55
  56 IMPLEMENT_STANDARD_HANDLE(HYDROData_ChannelAltitude, HYDROData_IAltitudeObject)
  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   TopoDS_Shape aShape =  aGuideXY->GetShape();
 121   gp_Pnt P1(thePoint.X(), thePoint.Y(), 0);
 122   TopoDS_Shape aPoint = BRepBuilderAPI_MakeVertex(P1).Shape();
 123
 124   if (aPoint.IsNull() || aShape.IsNull())
 125     {
 126       DEBTRACE("aPoint.IsNull() || aShape.IsNull()");
 127       return aResAltitude;
 128     }
 129
 130   if (aShape.ShapeType() != TopAbs_EDGE && aShape.ShapeType() != TopAbs_WIRE)
 131     {
 132       DEBTRACE("Projection aborted : the shape is neither an edge nor a wire");
 133       return aResAltitude;
 134     }
 135
 136   // Perform projection.
 137   BRepExtrema_DistShapeShape aDistShSh(aPoint, aShape, Extrema_ExtFlag_MIN);
 138
 139   if (aDistShSh.IsDone() == Standard_False)
 140     {
 141       DEBTRACE("Projection not done");
 142       return aResAltitude;
 143     }
 144
 145   Standard_Boolean hasValidSolution = Standard_False;
 146   Standard_Integer aNbSolutions = aDistShSh.NbSolution();
 147   Standard_Integer i;
 148   double aParam = 0.;
 149   Standard_Real aTolConf = BRep_Tool::Tolerance(TopoDS::Vertex(aPoint));
 150   Standard_Real aTolAng = 1.e-4;
 151
 152   for (i = 1; i <= aNbSolutions; i++)
 153     {
 154       Standard_Boolean isValid = Standard_False;
 155       BRepExtrema_SupportType aSupportType = aDistShSh.SupportTypeShape2(i);
 156       TopoDS_Shape aSupportShape = aDistShSh.SupportOnShape2(i);
 157
 158       if (aSupportType == BRepExtrema_IsOnEdge)
 159         {
 160           // Minimal distance inside edge is really a projection.
 161           isValid = Standard_True;
 162           aDistShSh.ParOnEdgeS2(i, aParam);
 163         }
 164       else if (aSupportType == BRepExtrema_IsVertex)
 165         {
 166           TopExp_Explorer anExp(aShape, TopAbs_EDGE);
 167
 168           if (aDistShSh.Value() <= aTolConf)
 169             {
 170               // The point lies on the shape. This means this point
 171               // is really a projection.
 172               for (; anExp.More() && !isValid; anExp.Next())
 173                 {
 174                   TopoDS_Edge aCurEdge = TopoDS::Edge(anExp.Current());
 175
 176                   if (aCurEdge.IsNull() == Standard_False)
 177                     {
 178                       TopoDS_Vertex aVtx[2];
 179
 180                       TopExp::Vertices(aCurEdge, aVtx[0], aVtx[1]);
 181
 182                       for (int j = 0; j < 2; j++)
 183                         {
 184                           if (aSupportShape.IsSame(aVtx[j]))
 185                             {
 186                               // The current edge is a projection edge.
 187                               isValid = Standard_True;
 188                               aSupportShape = aCurEdge;
 189                               aParam = BRep_Tool::Parameter(aVtx[j], aCurEdge);
 190                               break;
 191                             }
 192                         }
 193                     }
 194                 }
 195             }
 196           else
 197             {
 198               // Minimal distance to vertex is not always a real projection.
 199               gp_Pnt aPnt = BRep_Tool::Pnt(TopoDS::Vertex(aPoint));
 200               gp_Pnt aPrjPnt = BRep_Tool::Pnt(TopoDS::Vertex(aSupportShape));
 201               gp_Vec aDProjP(aPrjPnt, aPnt);
 202
 203               for (; anExp.More() && !isValid; anExp.Next())
 204                 {
 205                   TopoDS_Edge aCurEdge = TopoDS::Edge(anExp.Current());
 206
 207                   if (aCurEdge.IsNull() == Standard_False)
 208                     {
 209                       TopoDS_Vertex aVtx[2];
 210
 211                       TopExp::Vertices(aCurEdge, aVtx[0], aVtx[1]);
 212
 213                       for (int j = 0; j < 2; j++)
 214                         {
 215                           if (aSupportShape.IsSame(aVtx[j]))
 216                             {
 217                               // Check if the point is a projection to the current edge.
 218                               Standard_Real anEdgePars[2];
 219                               Handle(Geom_Curve) aCurve = BRep_Tool::Curve(aCurEdge, anEdgePars[0], anEdgePars[1]);
 220                               gp_Pnt aVal;
 221                               gp_Vec aD1;
 222
 223                               aParam = BRep_Tool::Parameter(aVtx[j], aCurEdge);
 224                               aCurve->D1(aParam, aVal, aD1);
 225
 226                               if (Abs(aD1.Dot(aDProjP)) <= aTolAng)
 227                                 {
 228                                   // The current edge is a projection edge.
 229                                   isValid = Standard_True;
 230                                   aSupportShape = aCurEdge;
 231                                   break;
 232                                 }
 233                             }
 234                         }
 235                     }
 236                 }
 237             }
 238         }
 239
 240       if (isValid)
 241         {
 242           if (hasValidSolution)
 243             {
 244               DEBTRACE("Projection aborted : multiple solutions");
 245               return aResAltitude;
 246             }
 247
 248           // Store the valid solution.
 249           hasValidSolution = Standard_True;
 250
 251           // profile altitude at projection point
 252           HYDROData_IPolyline::PointsList aProfilePoints = aGuideUZ->GetPoints();
 253           if ( aProfilePoints.IsEmpty() )
 254             {
 255               DEBTRACE("empty profile UZ");
 256               return aResAltitude;
 257             }
 258           double aDepth = HYDROData_ProfileUZ::GetDepthFromDistance( aProfilePoints, aParam );
 259           //DEBTRACE("profile altitude: " << aDepth);
 260
 261           // Compute edge index.
 262           TopExp_Explorer anExp(aShape, TopAbs_EDGE);
 263           int anIndex = 0;
 264           for (; anExp.More(); anExp.Next(), anIndex++)
 265             {
 266               if (aSupportShape.IsSame(anExp.Current()))
 267                 {
 268                   break;
 269                 }
 270             }
 271
 272           // get the XY distance from point to guideline
 273           const gp_Pnt &aPntProj = aDistShSh.PointOnShape2(i);
 274           //DEBTRACE("projection: (" << aPntProj.X() << ", " << aPntProj.Y() << ", " << aPntProj.Z() << ")");
 275           gp_XY aProjXY = gp_XY(aPntProj.X(), aPntProj.Y());
 276           aProjXY.Subtract(thePoint);
 277           double distance = aProjXY.Modulus();
 278           //DEBTRACE("distance to guideline " << distance);
 279
 280           // get delta altitude on section (supposed symmetric) from guideline distance (aParam)
 281           double delta = 0;
 282           int i1 = 0;
 283           gp_XY pt1 = gp_XY();
 284           gp_XY pt2 = gp_XY();
 285           HYDROData_ProfileUZ::PointsList aSectionPoints = aProfile->GetParametricPoints();
 286           for ( int i = 1, aNbPoints = aSectionPoints.Size(); i <= aNbPoints; ++i )
 287             {
 288               const HYDROData_IPolyline::Point& aPolylinePoint = aSectionPoints.Value( i );
 289               //DEBTRACE("  profile point: " << aPolylinePoint.X() << " " << aPolylinePoint.Y());
 290               if (aPolylinePoint.X() < distance)
 291                 {
 292                   i1 = i;
 293                   pt1 = aPolylinePoint;
 294                 }
 295               if (aPolylinePoint.X() >= distance)
 296                 {
 297                   pt2 = aPolylinePoint;
 298                   break;
 299                 }
 300             }
 301           if ((i1 == 0) && (distance > 0))
 302             {
 303               DEBTRACE("Projection aborted : non centered profile");
 304               return aResAltitude;
 305             }
 306           if (i1 == aProfilePoints.Size()) // distance >= profile width
 307             {
 308               delta = pt1.Y();
 309             }
 310           else
 311             {
 312               delta = pt1.Y() + (pt2.Y() - pt1.Y())*(distance -pt1.X())/(pt2.X()-pt1.X());
 313             }
 314           aResAltitude = delta + aDepth;
 315           DEBTRACE("distance XY: "<< aParam << " distance to guideline: " << distance << " final altitude: " << aResAltitude << " delta: " << delta);
 316           return aResAltitude;
 317         }
 318     }
 319
 320   if (!hasValidSolution)
 321     {
 322       DEBTRACE("Projection aborted : no projection");
 323       return aResAltitude;
 324     }
 325
 326   return aResAltitude;
 327 }
 328
 329
 330
 331