src/HYDROData/HYDROData_StreamLinearInterpolation.cxx

   1 // Copyright (C) 2007-2014  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include "HYDROData_StreamLinearInterpolation.h"
  24 #include <HYDROData_Profile.h>
  25 #include <HYDROData_ProfileUZ.h>
  26
  27 #include <gp_Lin2d.hxx>
  28 #include <BRepAdaptor_Curve.hxx>
  29 #include <HYDROData_ShapesTool.h>
  30 #include <GCPnts_UniformAbscissa.hxx>
  31 #include <Geom2dAPI_InterCurveCurve.hxx>
  32 #include <Geom2d_Line.hxx>
  33 #include <TopoDS_Shape.hxx>
  34 #include <TopoDS.hxx>
  35 #include <HYDROData_Tool.h>
  36 #include <GCE2d_MakeLine.hxx>
  37 #include <TopoDS_Iterator.hxx>
  38 #include <TopExp_Explorer.hxx>
  39 #include <BRepAlgo.hxx>
  40 #include <Geom2dAdaptor_Curve.hxx>
  41 #include <assert.h>
  42 #include <vector>
  43 #include <map>
  44 #include <set>
  45 #include <BRepBuilderAPI_MakeEdge2d.hxx>
  46 #include <BRepBuilderAPI_MakeEdge.hxx>
  47 #include <BRep_Tool.hxx>
  48 #include <BRepBuilderAPI_MakeWire.hxx>
  49 #include <TopoDS_Face.hxx>
  50 #include <gp_Pln.hxx>
  51 #include <BRepBuilderAPI_MakeFace.hxx>
  52 #include <gp.hxx>
  53 #include <gp_Ax2.hxx>
  54 #include <Geom_Plane.hxx>
  55 #include <TopoDS_Wire.hxx>
  56 #include <BRepLib.hxx>
  57
  58 static std::vector<gp_Pnt2d> GetProfileUZPoints(const Handle(HYDROData_Profile)& theProfile)
  59 {
  60   Handle(HYDROData_ProfileUZ) aProfileUZ = theProfile->GetProfileUZ( false );
  61   HYDROData_ProfileUZ::PointsList aSectPointsList = aProfileUZ->GetPoints();
  62   std::vector<gp_Pnt2d> points;
  63   points.reserve(aSectPointsList.Size());
  64   for ( int i = 1; i <= aSectPointsList.Size(); i++ )
  65   {
  66     const HYDROData_ProfileUZ::Point& aSectPoint = aSectPointsList(i);
  67     points.push_back( gp_Pnt2d(aSectPoint.X(), aSectPoint.Y()));
  68   }
  69   return points;
  70 }
  71
  72 static void GetMaxDist(const std::vector<gp_Pnt2d>& PNTS,
  73   double& dmax,
  74   double& dmax2,
  75   int& imax)
  76 {
  77   dmax = 0;
  78   dmax2 = 0;
  79   imax = -1;
  80   std::vector<double> dist;
  81   dist.reserve(PNTS.size()-1);
  82   for (int i=0; i<PNTS.size()-1; i++)
  83   {
  84     //dist.push_back(PNTS[i].Distance(PNTS[i+1])); //via distance between p1/p2
  85     dist.push_back(std::abs(PNTS[i].X()-PNTS[i+1].X())); //via U-diff of p1/p2
  86   }
  87   for (int i=0; i<dist.size(); i++)
  88   {
  89     double d = dist[i];
  90     if (d > dmax)
  91     {
  92       dmax = d;
  93       imax = i;
  94     }
  95   }
  96   for (int i=0; i<dist.size(); i++)
  97   {
  98     double d = dist[i];
  99     if (d > dmax2 && d < dmax)
 100       dmax2 = d;
 101   }
 102 }
 103
 104 static void InsertPoints(std::vector<gp_Pnt2d>& points, int nbpoints)
 105 {
 106   points.reserve(points.size() + nbpoints);
 107   while (nbpoints)
 108   {
 109     double dmax=0, dmax2=0;
 110     int imax=-1;
 111     GetMaxDist(points, dmax, dmax2, imax);
 112     int nbPins = 0;
 113     if (dmax2 != 0)
 114       nbPins = floor(dmax/dmax2);
 115     else
 116       nbPins = nbpoints; //one segment, put all points here
 117     //insert nbPins points in [imax, imax+1] segment
 118     gp_Pnt2d p1 = points[imax];
 119     gp_Pnt2d p2 = points[imax+1];
 120     double X0 = p1.X();
 121     double Y0 = p1.Y();
 122     double X1 = p2.X();
 123     double Y1 = p2.Y();
 124     for (int i=0;i<nbPins;i++)
 125     {
 126       double t = ((double)i+1)/((double)nbPins+1);
 127       double Xp = X0 + (X1-X0)*t;
 128       double Yp = Y0 + (Y1-Y0)*t;
 129       gp_Pnt2d P(Xp, Yp);
 130       points.insert(points.begin() + imax + 1, P);
 131     }
 132     nbpoints-=nbPins;
 133   }
 134 }
 135
 136 static void PolyToCurve2d(const Handle(HYDROData_PolylineXY)& poly, Handle(Geom2d_Curve)& c2d)
 137 {
 138   if (poly.IsNull())
 139     return;
 140   TopoDS_Shape sh = poly->GetShape();
 141   if (sh.ShapeType() == TopAbs_COMPOUND)
 142   {
 143     TopoDS_Iterator it(sh);
 144     if (it.More())
 145       sh = it.Value();
 146   }
 147
 148   if (sh.IsNull())
 149     return;
 150
 151   TopoDS_Edge ee;
 152   if (sh.ShapeType() == TopAbs_EDGE)
 153     ee = TopoDS::Edge(sh);
 154   else if (sh.ShapeType() == TopAbs_WIRE)
 155     ee = BRepAlgo::ConcatenateWireC0(TopoDS::Wire(sh)); //convert to bspline
 156
 157   if (ee.IsNull())
 158     return;
 159
 160   BRepAdaptor_Curve Ad(ee);
 161
 162   c2d = HYDROData_Tool::BRepAdaptorTo2DCurve(Ad);
 163 }
 164
 165 static void InterProfilesAndHAX(const HYDROData_SequenceOfObjects& profiles,
 166   const Handle(Geom2d_Curve)& Hax2d,
 167   std::map<double, Handle(HYDROData_Profile)>& profToInterParam,
 168   std::vector<std::string>* warnings)
 169 {
 170   for (int i=1;i<=profiles.Size();i++)
 171   {
 172     Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast(profiles(i));
 173     gp_XY LP, RP;
 174     aProfile->GetLeftPoint( LP, true, true );
 175     aProfile->GetRightPoint( RP, true, true );
 176     gp_Pnt2d P1(LP), P2(RP);
 177     double d = P2.Distance(P1);
 178     if (d < gp::Resolution())
 179     {
 180       if (warnings)
 181         warnings->push_back(aProfile->GetName().toStdString() + " is skipped: Left and Right points is the same");
 182       continue;
 183     }
 184     Handle(Geom2d_Line) lin2d = new Geom2d_Line(P1,gp_Dir2d(P2.XY()-P1.XY()));
 185
 186     Geom2dAdaptor_Curve linAd(lin2d, 0, d);
 187     Geom2dAdaptor_Curve haxAd(Hax2d);
 188     Geom2dInt_GInter isec( linAd, haxAd, 1.0e-6, 1.0e-6);
 189     if (!isec.IsDone())
 190     {
 191       if (warnings)
 192         warnings->push_back(aProfile->GetName().toStdString() + " is skipped: intersection between hydraulic axis & profile is failed");
 193       continue;
 194     }
 195     if (isec.NbPoints() == 0)
 196     {
 197       if (warnings)
 198         warnings->push_back(aProfile->GetName().toStdString() + " is skipped: no intersection between hydraulic axis & profile");
 199       continue;
 200     }
 201     if (isec.NbPoints() > 1)
 202     {
 203       if (warnings)
 204         warnings->push_back(aProfile->GetName().toStdString() + " is skipped; intersection between hydraulic axis & profile produces more than 1 point");
 205       continue;
 206     }
 207     double param = isec.Point(1).ParamOnSecond();
 208     profToInterParam[param] = aProfile;
 209   }
 210 }
 211
 212 static bool EquidParamsOnHAX(const Handle(Geom2d_Curve)& Hax2d, double step, double fp, double lp, std::set<double>& params)
 213 {
 214   Geom2dAdaptor_Curve ad(Hax2d);
 215   GCPnts_UniformAbscissa GUA(ad, step, fp, lp);
 216   if (!GUA.IsDone())
 217     return false;
 218
 219   for(int i = 1; i<= GUA.NbPoints(); i++)
 220     params.insert(GUA.Parameter(i));
 221
 222   return true;
 223 }
 224
 225 static void GetPointsOnBanks(const std::vector<double>& gua_params,
 226   const Handle(Geom2d_Curve)& Hax2d,
 227   const Handle(Geom2d_Curve)& LB2d,
 228   const Handle(Geom2d_Curve)& RB2d,
 229   std::map<double, std::pair<gp_Pnt2d, gp_Pnt2d>>& parToBankPoints,
 230   std::vector<std::string>* warnings)
 231 {
 232   for(int i = 0; i < gua_params.size(); i++)
 233   {
 234     double par = gua_params[i];
 235     gp_Pnt2d p;
 236     gp_Vec2d v;
 237     Hax2d->D1(par, p, v); //if non-diff?
 238     gp_Pnt2d outpLB, outpRB;
 239
 240     gp_Dir2d n2d(v.X(), v.Y()); //ignore Z (Z==0)
 241     n2d.Rotate(M_PI/2);
 242     //nLine2d is an orthogonal line to hydr. axis
 243     gp_Pnt2d p2d(gp_Pnt2d(p.X(), p.Y()));
 244     Handle(Geom2d_Line) nLine2d = new Geom2d_Line(p2d, n2d);
 245     //intersect nLine2d with banks
 246     Geom2dAPI_InterCurveCurve intL(LB2d, nLine2d);
 247     Geom2dAPI_InterCurveCurve intR(RB2d, nLine2d);
 248
 249     if (warnings && intL.NbPoints()==0 )
 250     {
 251       std::string coord = (QString::number(p2d.X()) + ", "+QString::number(p2d.Y())).toStdString();
 252       warnings->push_back("no intersection between intermediate profile (point on h.axis:("+ coord +")) and left bank found; skipped");
 253     }
 254
 255     if (warnings && intR.NbPoints()==0)
 256     {
 257       std::string coord = (QString::number(p2d.X()) + ", "+QString::number(p2d.Y())).toStdString();
 258       warnings->push_back("no intersection between intermediate profile (point on h.axis:("+ coord +")) and right bank found; skipped");
 259     }
 260
 261     if ( intL.NbPoints()==0 || intR.NbPoints()==0)
 262       continue;
 263
 264     gp_Pnt2d aNearSolL, aNearSolR;
 265     double min_sq_dist = Precision::Infinite();
 266     for (int j=1;j<=intL.NbPoints();j++)
 267     {
 268       double sq_dist = p2d.SquareDistance(intL.Point(j));
 269       if (min_sq_dist > sq_dist)
 270       {
 271         min_sq_dist = sq_dist;
 272         aNearSolL = intL.Point(j);
 273       }
 274     }
 275     min_sq_dist = Precision::Infinite();
 276     for (int j=1;j<=intR.NbPoints();j++)
 277     {
 278       double sq_dist = p2d.SquareDistance(intR.Point(j));
 279       if (min_sq_dist > sq_dist)
 280       {
 281         min_sq_dist = sq_dist;
 282         aNearSolR = intR.Point(j);
 283       }
 284     }
 285
 286     std::pair<gp_Pnt2d, gp_Pnt2d> int_pair(aNearSolL, aNearSolR);
 287     parToBankPoints[par]=int_pair;
 288   }
 289 }
 290
 291
 292 bool GetItersectionParam(const Handle(HYDROData_Profile)& aProfile, const Handle(Geom2d_Curve)& B2d, double& outIntersectionParam)
 293 {
 294   gp_XY LP, RP;
 295   aProfile->GetLeftPoint( LP, true, true );
 296   aProfile->GetRightPoint( RP, true, true );
 297   gp_Pnt2d P1(LP), P2(RP);
 298   double d = P2.Distance(P1);
 299   if (d < gp::Resolution())
 300     return false;
 301   Handle(Geom2d_Line) lin2d = new Geom2d_Line(P1,gp_Dir2d(P2.XY()-P1.XY()));
 302
 303   Geom2dAdaptor_Curve linAd(lin2d, 0, d);
 304   Geom2dAdaptor_Curve BAd(B2d);
 305   Geom2dInt_GInter isec( linAd, BAd, 1.0e-6, 1.0e-6);
 306   if (!isec.IsDone())
 307     return false;
 308   if (isec.NbPoints() == 0)
 309     return false;
 310   if (isec.NbPoints() > 1)
 311     return false;
 312   outIntersectionParam = isec.Point(1).ParamOnSecond();
 313   return true;
 314 }
 315
 316 void BuildFace(const Handle(HYDROData_Profile)& F_prof,
 317   const Handle(HYDROData_Profile)& S_prof,
 318   const Handle(Geom2d_Curve)& LB2d,
 319   const Handle(Geom2d_Curve)& RB2d,
 320   HYDROData_Stream::PrsDefinition& prsDef)
 321 {
 322   prsDef.myPrs3D.Nullify();
 323   prsDef.myPrs2D.Nullify();
 324   prsDef.myLeftBank.Nullify();
 325   prsDef.myRightBank.Nullify();
 326   prsDef.myInlet.Nullify();
 327   prsDef.myOutlet.Nullify();
 328
 329   bool stat = true;
 330   double lb_1, lb_2, rb_1, rb_2;
 331   stat = GetItersectionParam(F_prof, LB2d, lb_1) && GetItersectionParam(F_prof, RB2d, rb_1)
 332     && GetItersectionParam(S_prof, LB2d, lb_2) && GetItersectionParam(S_prof, RB2d, rb_2);
 333   if (stat)
 334   {
 335     BRepBuilderAPI_MakeEdge2d LEM(LB2d, lb_1, lb_2);
 336     if (!LEM.IsDone())
 337       return;
 338     BRepBuilderAPI_MakeEdge2d REM(RB2d, rb_1, rb_2);
 339     if (!REM.IsDone())
 340       return;
 341     TopoDS_Edge LBE = LEM.Edge();
 342     TopoDS_Edge RBE = REM.Edge();
 343     if (LBE.IsNull() || RBE.IsNull())
 344       return;
 345     BRepBuilderAPI_MakeEdge PFEM(BRep_Tool::Pnt(LEM.Vertex1()), BRep_Tool::Pnt(REM.Vertex1()));
 346     BRepBuilderAPI_MakeEdge PLEM(BRep_Tool::Pnt(LEM.Vertex2()), BRep_Tool::Pnt(REM.Vertex2()));
 347     if (!PFEM.IsDone())
 348       return;
 349     if (!PLEM.IsDone())
 350       return;
 351     TopoDS_Edge FProfE = PFEM.Edge();
 352     TopoDS_Edge SProfE = PLEM.Edge();
 353     BRepBuilderAPI_MakeWire WM(FProfE, LBE, SProfE, RBE);
 354     if (WM.IsDone())
 355     {
 356       TopoDS_Wire W = WM.Wire();
 357       BRepLib::BuildCurves3d(W);
 358       BRepBuilderAPI_MakeFace FM(Geom_Plane(gp::XOY()).Pln(), WM.Wire());
 359       if (FM.IsDone())
 360       {
 361         prsDef.myPrs2D = FM.Face();
 362         prsDef.myLeftBank = LBE;
 363         prsDef.myRightBank = RBE;
 364         prsDef.myInlet = FProfE;
 365         prsDef.myOutlet = SProfE;
 366       }
 367     }
 368   }
 369 }
 370
 371 void HYDROData_StreamLinearInterpolation::Perform(const HYDROData_SequenceOfObjects& profiles,
 372   int pointsToInsert,
 373   double stepOnHA,
 374   const Handle(HYDROData_PolylineXY)& hax,
 375   const Handle(HYDROData_PolylineXY)& LB,
 376   const Handle(HYDROData_PolylineXY)& RB,
 377   HYDROData_Bathymetry::AltitudePoints& outBathypoints,
 378   bool buildPresentationShapes,
 379   bool estimateWarnOnly,
 380   HYDROData_Stream::PrsDefinition& prsDef,
 381   std::vector<std::string>* warnings)
 382 {
 383   //intersect profiles with given hydr.axis
 384   if (hax.IsNull() || LB.IsNull() || RB.IsNull())
 385     return;
 386
 387   Handle(Geom2d_Curve) Hax2d;
 388   PolyToCurve2d(hax, Hax2d);
 389   if (Hax2d->IsClosed()) //can't be closed
 390   {
 391     if (warnings)
 392       warnings->push_back(hax->GetName().toStdString() + " is closed; abort");
 393     return;
 394   }
 395
 396   std::map<double, Handle(HYDROData_Profile)> InterParamToProf;
 397   //profToInterParam is output param map: profile to intersection param on hydr.axis
 398   InterProfilesAndHAX(profiles, Hax2d, InterParamToProf, warnings);
 399   //
 400   std::vector<std::vector<gp_Pnt2d>> profilesPoints;
 401   profilesPoints.reserve(InterParamToProf.size());
 402   int maxNbPoints = 0;
 403
 404   //for (int i=1;i<=InterParamToProf.Extent();i++)
 405   for( std::map<double, Handle(HYDROData_Profile)>::iterator it = InterParamToProf.begin(); it != InterParamToProf.end(); ++it )
 406   {
 407     //ordered by param on hydr axis
 408     Handle(HYDROData_Profile) aProfile = it->second;
 409     const std::vector<gp_Pnt2d>& profile_points = GetProfileUZPoints(aProfile);
 410     profilesPoints.push_back(profile_points);
 411     if (profile_points.size() > maxNbPoints)
 412       maxNbPoints = profile_points.size();
 413   }
 414   std::set<double> paramsOnHAX;
 415
 416   if (InterParamToProf.size() <= 1)
 417   {
 418     if (warnings)
 419       warnings->push_back("Insufficient count of correct input profiles; abort");
 420     return;
 421   }
 422
 423   EquidParamsOnHAX(Hax2d, stepOnHA, InterParamToProf.begin()->first,
 424     InterParamToProf.rbegin()->first, paramsOnHAX);
 425
 426   //prepare ParamsPerSegm - number of interm.profiles for profile_points
 427   std::vector<std::vector<double>> ParamsPerSegm;
 428   std::set<double>::iterator it_params_hax = paramsOnHAX.begin();
 429   std::map<double, Handle(HYDROData_Profile)>::iterator it_p = InterParamToProf.begin();
 430   it_p++;
 431   ParamsPerSegm.resize(profilesPoints.size()-1);
 432   for( int k=0; it_p != InterParamToProf.end(); )
 433   {
 434     if (*it_params_hax < it_p->first)
 435     {
 436       ParamsPerSegm[k].push_back(*it_params_hax);
 437       it_params_hax++;
 438     }
 439     else
 440     {
 441       it_p++;
 442       k++;
 443     }
 444   }
 445   ParamsPerSegm.back().push_back(*paramsOnHAX.rbegin());
 446   //
 447   std::map<double, std::pair<gp_Pnt2d, gp_Pnt2d>> parToBankPoints;
 448   std::vector<double> paramHAXVec;
 449   paramHAXVec.reserve(paramsOnHAX.size());
 450   for( std::set<double>::iterator it = paramsOnHAX.begin(); it != paramsOnHAX.end(); ++it )
 451     paramHAXVec.push_back(*it);
 452   //
 453   Handle(Geom2d_Curve) LB2d, RB2d;
 454   PolyToCurve2d(LB, LB2d);
 455   PolyToCurve2d(RB, RB2d);
 456   GetPointsOnBanks( paramHAXVec, Hax2d, LB2d, RB2d, parToBankPoints, warnings);
 457
 458   if (buildPresentationShapes)
 459   {
 460     Handle(HYDROData_Profile) F_prof = InterParamToProf.begin()->second;
 461     Handle(HYDROData_Profile) S_prof = InterParamToProf.rbegin()->second;
 462     BuildFace(F_prof, S_prof, LB2d, RB2d, prsDef);
 463   }
 464   //
 465
 466   if (estimateWarnOnly)
 467     return;
 468
 469   maxNbPoints = Max(pointsToInsert, maxNbPoints);
 470
 471   //insert points to profiles
 472   for (int i=0;i<profilesPoints.size();i++)
 473   {
 474     int nbPointsToInsert = maxNbPoints - profilesPoints[i].size();
 475     //insert nbPointsToInsert points to current profile
 476     InsertPoints(profilesPoints[i], nbPointsToInsert);
 477   }
 478
 479   for (int i=0;i<profilesPoints.size()-1;i++)
 480   {
 481     const std::vector<gp_Pnt2d>& prof1 = profilesPoints[i];
 482     const std::vector<gp_Pnt2d>& prof2 = profilesPoints[i+1];
 483     std::vector<std::vector<gp_Pnt2d>> IntermProf;
 484
 485     std::map<int, double> indToParam; //index of im.profile to param on hax. ()
 486     for (int j=0;j<ParamsPerSegm[i].size();j++)
 487     {
 488       double param = ParamsPerSegm[i][j];
 489       if (parToBankPoints.count(param) > 0) //not intersected with banks; skip this
 490         indToParam[j]=param;
 491     }
 492     int NbIntermProf = indToParam.size();
 493     IntermProf.resize(NbIntermProf);
 494
 495     for (int l=0;l<NbIntermProf;l++)
 496       IntermProf[l].resize(prof1.size());
 497     //calc interm. between prof1 & prof2
 498     assert (prof1.size() == prof2.size());
 499     for (int j=0;j<prof1.size();j++)
 500     {
 501       gp_Pnt2d P1 = prof1[j];
 502       gp_Pnt2d P2 = prof2[j];
 503       double X0 = P1.X();
 504       double Y0 = P1.Y();
 505       double X1 = P2.X();
 506       double Y1 = P2.Y();
 507       std::map<int, double>::iterator it = indToParam.begin();
 508       for ( int m=0; it != indToParam.end(); it++, m++ )
 509       {
 510         int ind = it->first;
 511         double t = ((double)ind+1)/((double)NbIntermProf+1);
 512         double Xp = X0 + (X1-X0)*t;
 513         double Yp = Y0 + (Y1-Y0)*t;
 514         gp_Pnt2d P(Xp, Yp);
 515         IntermProf[m][j] = P;
 516       }
 517     }
 518     std::map<int, double>::iterator it = indToParam.begin();
 519     for ( int m=0; it != indToParam.end(); it++, m++ )
 520     {
 521       const std::vector<gp_Pnt2d>& im_prof = IntermProf[m];
 522       double param = it->second;
 523       const std::pair<gp_Pnt2d, gp_Pnt2d>& BB = parToBankPoints[param]; //param is included in map; no checks
 524       gp_Pnt2d LP = BB.first;
 525       gp_Pnt2d RP = BB.second;
 526       HYDROData_ProfileUZ::PointsList pl;
 527       for (int k=0;k<im_prof.size();k++)
 528         pl.Append(im_prof[k].XY());
 529       HYDROData_Profile::ProfilePoints profile_points3d = HYDROData_Profile::CalculateProfilePoints(pl, LP.XY(), RP.XY());
 530       for (int k=1;k<=profile_points3d.Length();k++)
 531       {
 532         HYDROData_Bathymetry::AltitudePoint AP(profile_points3d(k).X(), profile_points3d(k).Y(), profile_points3d(k).Z());
 533         outBathypoints.push_back(AP);
 534       }
 535     }
 536   }
 537   //outBathy->SetAltitudePoints(bathypoints);
 538 }
 539
 540
 541