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 #include <BRepLib_MakePolygon.hxx>
  58 #include <TopoDS_Compound.hxx>
  59 #include <BRep_Builder.hxx>
  60
  61 static std::vector<gp_Pnt2d> GetProfileUZPoints(const Handle(HYDROData_Profile)& theProfile)
  62 {
  63   Handle(HYDROData_ProfileUZ) aProfileUZ = theProfile->GetProfileUZ( false );
  64   HYDROData_ProfileUZ::PointsList aSectPointsList = aProfileUZ->GetPoints();
  65   std::vector<gp_Pnt2d> points;
  66   points.reserve(aSectPointsList.Size());
  67   for ( int i = 1; i <= aSectPointsList.Size(); i++ )
  68   {
  69     const HYDROData_ProfileUZ::Point& aSectPoint = aSectPointsList(i);
  70     points.push_back( gp_Pnt2d(aSectPoint.X(), aSectPoint.Y()));
  71   }
  72   return points;
  73 }
  74
  75 static void GetMaxDist(const std::vector<gp_Pnt2d>& PNTS,
  76   double& dmax,
  77   double& dmax2,
  78   int& imax)
  79 {
  80   dmax = 0;
  81   dmax2 = 0;
  82   imax = -1;
  83   std::vector<double> dist;
  84   dist.reserve(PNTS.size()-1);
  85   for (int i=0; i<PNTS.size()-1; i++)
  86   {
  87     //dist.push_back(PNTS[i].Distance(PNTS[i+1])); //via distance between p1/p2
  88     dist.push_back(std::abs(PNTS[i].X()-PNTS[i+1].X())); //via U-diff of p1/p2
  89   }
  90   for (int i=0; i<dist.size(); i++)
  91   {
  92     double d = dist[i];
  93     if (d > dmax)
  94     {
  95       dmax = d;
  96       imax = i;
  97     }
  98   }
  99   for (int i=0; i<dist.size(); i++)
 100   {
 101     double d = dist[i];
 102     if (d > dmax2 && d < dmax)
 103       dmax2 = d;
 104   }
 105 }
 106
 107 static void InsertPoints(std::vector<gp_Pnt2d>& points, int nbpoints)
 108 {
 109   points.reserve(points.size() + nbpoints);
 110   while (nbpoints>0)
 111   {
 112     double dmax=0, dmax2=0;
 113     int imax=-1;
 114     GetMaxDist(points, dmax, dmax2, imax);
 115     int nbPins = 0;
 116     if (dmax2 != 0)
 117       nbPins = floor(dmax/dmax2);
 118     else
 119       nbPins = nbpoints; //one segment, put all points here
 120     //insert nbPins points in [imax, imax+1] segment
 121     gp_Pnt2d p1 = points[imax];
 122     gp_Pnt2d p2 = points[imax+1];
 123     double X0 = p1.X();
 124     double Y0 = p1.Y();
 125     double X1 = p2.X();
 126     double Y1 = p2.Y();
 127     if (nbPins>nbpoints) //deny to add more points than needed
 128       nbPins = nbpoints;
 129     for (int i=0;i<nbPins;i++)
 130     {
 131       double t = ((double)i+1)/((double)nbPins+1);
 132       double Xp = X0 + (X1-X0)*t;
 133       double Yp = Y0 + (Y1-Y0)*t;
 134       gp_Pnt2d P(Xp, Yp);
 135       points.insert(points.begin() + imax + 1, P);
 136     }
 137     nbpoints-=nbPins;
 138   }
 139   assert (nbpoints == 0);
 140 }
 141
 142 static void PolyToCurve2d(const Handle(HYDROData_PolylineXY)& poly, Handle(Geom2d_Curve)& c2d)
 143 {
 144   if (poly.IsNull())
 145     return;
 146   TopoDS_Shape sh = poly->GetShape();
 147   if (sh.ShapeType() == TopAbs_COMPOUND)
 148   {
 149     TopoDS_Iterator it(sh);
 150     if (it.More())
 151       sh = it.Value();
 152   }
 153
 154   if (sh.IsNull())
 155     return;
 156
 157   TopoDS_Edge ee;
 158   if (sh.ShapeType() == TopAbs_EDGE)
 159     ee = TopoDS::Edge(sh);
 160   else if (sh.ShapeType() == TopAbs_WIRE)
 161     ee = BRepAlgo::ConcatenateWireC0(TopoDS::Wire(sh)); //convert to bspline
 162
 163   if (ee.IsNull())
 164     return;
 165
 166   BRepAdaptor_Curve Ad(ee);
 167
 168   c2d = HYDROData_Tool::BRepAdaptorTo2DCurve(Ad);
 169 }
 170
 171 static void InterProfilesAndHAX(const HYDROData_SequenceOfObjects& profiles,
 172   const Handle(Geom2d_Curve)& Hax2d,
 173   std::map<double, Handle(HYDROData_Profile)>& profToInterParam,
 174   std::vector<bool>& paramToSwapFlag,
 175   std::vector<std::string>* warnings)
 176 {
 177   for (int i=1;i<=profiles.Size();i++)
 178   {
 179     Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast(profiles(i));
 180     gp_XY LP, RP;
 181     aProfile->GetLeftPoint( LP, false, true );
 182     aProfile->GetRightPoint( RP, false, true );
 183     gp_Pnt2d P1(LP), P2(RP);
 184     double d = P2.Distance(P1);
 185     if (d < gp::Resolution())
 186     {
 187       if (warnings)
 188         warnings->push_back(aProfile->GetName().toStdString() + " is skipped: Left and Right points is the same");
 189       continue;
 190     }
 191     Handle(Geom2d_Line) lin2d = new Geom2d_Line(P1,gp_Dir2d(P2.XY()-P1.XY()));
 192     gp_Dir2d prof_dir2d = lin2d->Direction();
 193
 194     Geom2dAdaptor_Curve linAd(lin2d, 0, d);
 195     Geom2dAdaptor_Curve haxAd(Hax2d);
 196     Geom2dInt_GInter isec( linAd, haxAd, 1.0e-6, 1.0e-6);
 197     if (!isec.IsDone())
 198     {
 199       if (warnings)
 200         warnings->push_back(aProfile->GetName().toStdString() + " is skipped: intersection between hydraulic axis & profile is failed");
 201       continue;
 202     }
 203     if (isec.NbPoints() == 0)
 204     {
 205       if (warnings)
 206         warnings->push_back(aProfile->GetName().toStdString() + " is skipped: no intersection between hydraulic axis & profile");
 207       continue;
 208     }
 209     if (isec.NbPoints() > 1)
 210     {
 211       if (warnings)
 212         warnings->push_back(aProfile->GetName().toStdString() + " is skipped; intersection between hydraulic axis & profile produces more than 1 point");
 213       continue;
 214     }
 215     double param = isec.Point(1).ParamOnSecond();
 216     gp_Pnt2d int_p2d;
 217     gp_Vec2d vec_on_hax_intersec;
 218     haxAd.D1(param, int_p2d, vec_on_hax_intersec);
 219     gp_Dir2d d_hax(vec_on_hax_intersec);
 220     double ang = d_hax.Angle(prof_dir2d);
 221     if (ang>0)
 222       paramToSwapFlag.push_back(true);
 223     else
 224       paramToSwapFlag.push_back(false);
 225     profToInterParam[param] = aProfile;
 226   }
 227 }
 228
 229 static bool EquidParamsOnHAX(const Handle(Geom2d_Curve)& Hax2d, double step, double fp, double lp, std::set<double>& params)
 230 {
 231   Geom2dAdaptor_Curve ad(Hax2d);
 232   GCPnts_UniformAbscissa GUA(ad, step, fp, lp);
 233   if (!GUA.IsDone())
 234     return false;
 235
 236   for(int i = 1; i<= GUA.NbPoints(); i++)
 237     params.insert(GUA.Parameter(i));
 238
 239   return true;
 240 }
 241
 242 static void GetPointsOnBanks(const std::vector<double>& gua_params,
 243   const Handle(Geom2d_Curve)& Hax2d,
 244   const Handle(Geom2d_Curve)& LB2d,
 245   const Handle(Geom2d_Curve)& RB2d,
 246   std::map<double, std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>>& parToBankPoints,
 247   std::vector<std::string>* warnings,
 248   const std::map<double, std::pair<Handle(HYDROData_Profile), Handle(HYDROData_Profile)>>* intermParamToProfPair)
 249 {
 250   for(int i = 0; i < gua_params.size(); i++)
 251   {
 252     double par = gua_params[i];
 253     gp_Pnt2d p;
 254     gp_Vec2d v;
 255     Hax2d->D1(par, p, v); //if non-diff?
 256     gp_Pnt2d outpLB, outpRB;
 257
 258     gp_Dir2d n2d(v.X(), v.Y()); //ignore Z (Z==0)
 259     n2d.Rotate(M_PI/2);
 260     //nLine2d is an orthogonal line to hydr. axis
 261     gp_Pnt2d p2d(gp_Pnt2d(p.X(), p.Y()));
 262     Handle(Geom2d_Line) nLine2d = new Geom2d_Line(p2d, n2d);
 263     //intersect nLine2d with banks
 264     Geom2dAPI_InterCurveCurve intL(LB2d, nLine2d);
 265     Geom2dAPI_InterCurveCurve intR(RB2d, nLine2d);
 266
 267     if (warnings && intL.NbPoints()==0 )
 268     {
 269       std::string coord = (QString::number(p2d.X()) + ", "+QString::number(p2d.Y())).toStdString();
 270       if (intermParamToProfPair && intermParamToProfPair->count(par) != 0)
 271       {
 272         std::string fpn = intermParamToProfPair->at(par).first->GetName().toStdString();
 273         std::string spn = intermParamToProfPair->at(par).second->GetName().toStdString();
 274         warnings->push_back("no intersection between intermediate profile (point on h.axis:("+
 275           coord +"), located between " + fpn + " & " + spn + ") and left bank found; skipped");
 276       }
 277       else
 278         warnings->push_back("no intersection between intermediate profile (point on h.axis:("+ coord +")) and left bank found; skipped");
 279     }
 280
 281     if (warnings && intR.NbPoints()==0)
 282     {
 283       std::string coord = (QString::number(p2d.X()) + ", "+QString::number(p2d.Y())).toStdString();
 284       if (intermParamToProfPair && intermParamToProfPair->count(par) != 0)
 285       {
 286         std::string fpn = intermParamToProfPair->at(par).first->GetName().toStdString();
 287         std::string spn = intermParamToProfPair->at(par).second->GetName().toStdString();
 288         warnings->push_back("no intersection between intermediate profile (point on h.axis:("+
 289           coord +"), located between " + fpn + " & " + spn + ") and right bank found; skipped");
 290       }
 291       else
 292         warnings->push_back("no intersection between intermediate profile (point on h.axis:("+ coord +")) and right bank found; skipped");
 293     }
 294
 295     if ( intL.NbPoints()==0 || intR.NbPoints()==0)
 296       continue;
 297
 298     IntRes2d_IntersectionPoint aNearSolL, aNearSolR;
 299     double min_sq_dist = Precision::Infinite();
 300     for (int j=1;j<=intL.NbPoints();j++)
 301     {
 302       double sq_dist = p2d.SquareDistance(intL.Point(j));
 303       if (min_sq_dist > sq_dist)
 304       {
 305         min_sq_dist = sq_dist;
 306         //aNearSolL = intL.Point(j);
 307         aNearSolL = intL.Intersector().Point(j);
 308       }
 309     }
 310     min_sq_dist = Precision::Infinite();
 311     for (int j=1;j<=intR.NbPoints();j++)
 312     {
 313       double sq_dist = p2d.SquareDistance(intR.Point(j));
 314       if (min_sq_dist > sq_dist)
 315       {
 316         min_sq_dist = sq_dist;
 317         //aNearSolR = intR.Point(j);
 318         aNearSolR =  intR.Intersector().Point(j);
 319       }
 320     }
 321
 322     std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint> int_pair(aNearSolL, aNearSolR);
 323     parToBankPoints[par]=int_pair;
 324   }
 325 }
 326
 327
 328 void BuildFace(const std::map<double, std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>>& parToBankPoints,
 329   const Handle(Geom2d_Curve)& LB2d,
 330   const Handle(Geom2d_Curve)& RB2d,
 331   HYDROData_Stream::PrsDefinition& prsDef)
 332 {
 333   prsDef.myPrs3D.Nullify();
 334   prsDef.myPrs2D.Nullify();
 335   prsDef.myLeftBank.Nullify();
 336   prsDef.myRightBank.Nullify();
 337   prsDef.myInlet.Nullify();
 338   prsDef.myOutlet.Nullify();
 339
 340   double lb_1, lb_2, rb_1, rb_2;
 341   lb_1 = parToBankPoints.begin()->second.first.ParamOnFirst();
 342   rb_1 = parToBankPoints.begin()->second.second.ParamOnFirst();
 343
 344   lb_2 = parToBankPoints.rbegin()->second.first.ParamOnFirst();
 345   rb_2 = parToBankPoints.rbegin()->second.second.ParamOnFirst();
 346
 347   //stat = GetItersectionParam(F_prof, LB2d, lb_1) && GetItersectionParam(F_prof, RB2d, rb_1)
 348   //  && GetItersectionParam(S_prof, LB2d, lb_2) && GetItersectionParam(S_prof, RB2d, rb_2);
 349
 350   BRepBuilderAPI_MakeEdge2d LEM(LB2d, lb_1, lb_2);
 351   if (!LEM.IsDone())
 352     return;
 353   BRepBuilderAPI_MakeEdge2d REM(RB2d, rb_1, rb_2);
 354   if (!REM.IsDone())
 355     return;
 356   TopoDS_Edge LBE = LEM.Edge();
 357   TopoDS_Edge RBE = REM.Edge();
 358   if (LBE.IsNull() || RBE.IsNull())
 359     return;
 360   BRepBuilderAPI_MakeEdge PFEM(BRep_Tool::Pnt(LEM.Vertex1()), BRep_Tool::Pnt(REM.Vertex1()));
 361   BRepBuilderAPI_MakeEdge PLEM(BRep_Tool::Pnt(LEM.Vertex2()), BRep_Tool::Pnt(REM.Vertex2()));
 362   if (!PFEM.IsDone())
 363     return;
 364   if (!PLEM.IsDone())
 365     return;
 366   TopoDS_Edge FProfE = PFEM.Edge();
 367   TopoDS_Edge SProfE = PLEM.Edge();
 368   BRepBuilderAPI_MakeWire WM(FProfE, LBE, SProfE, RBE);
 369   if (WM.IsDone())
 370   {
 371     TopoDS_Wire W = WM.Wire();
 372     BRepLib::BuildCurves3d(W);
 373     BRepBuilderAPI_MakeFace FM(Geom_Plane(gp::XOY()).Pln(), WM.Wire());
 374     if (FM.IsDone())
 375     {
 376       prsDef.myPrs2D = FM.Face();
 377       prsDef.myLeftBank = LBE;
 378       prsDef.myRightBank = RBE;
 379       prsDef.myInlet = FProfE;
 380       prsDef.myOutlet = SProfE;
 381     }
 382   }
 383 }
 384
 385 void reversePoints(std::vector<gp_Pnt2d>& points)
 386 {
 387   std::reverse(points.begin(), points.end());
 388 }
 389
 390 void HYDROData_StreamLinearInterpolation::Perform(const HYDROData_SequenceOfObjects& profiles,
 391   int pointsToInsert,
 392   double stepOnHA,
 393   const Handle(HYDROData_PolylineXY)& hax,
 394   const Handle(HYDROData_PolylineXY)& LB,
 395   const Handle(HYDROData_PolylineXY)& RB,
 396   HYDROData_Bathymetry::AltitudePoints& outBathypoints,
 397   bool buildPresentationShapes,
 398   bool estimateWarnOnly,
 399   HYDROData_Stream::PrsDefinition& prsDef,
 400   std::vector<std::string>* warnings)
 401 {
 402   //intersect profiles with given hydr.axis
 403   if (hax.IsNull() || LB.IsNull() || RB.IsNull())
 404     return;
 405
 406   Handle(Geom2d_Curve) Hax2d;
 407   PolyToCurve2d(hax, Hax2d);
 408   if (hax->IsClosed() || Hax2d->IsClosed()) //can't be closed
 409   {
 410     if (warnings)
 411       warnings->push_back(hax->GetName().toStdString() + " is closed; abort");
 412     return;
 413   }
 414
 415   std::map<double, Handle(HYDROData_Profile)> InterParamToProf;
 416   std::vector<bool> paramToSwapFlag;
 417   //profToInterParam is output param map: profile to intersection param on hydr.axis
 418   InterProfilesAndHAX(profiles, Hax2d, InterParamToProf, paramToSwapFlag, warnings);
 419   //
 420   std::vector<std::vector<gp_Pnt2d>> profilesPoints;
 421   profilesPoints.reserve(InterParamToProf.size());
 422   int maxNbPoints = 0;
 423
 424   //for (int i=1;i<=InterParamToProf.Extent();i++)
 425   for( std::map<double, Handle(HYDROData_Profile)>::iterator it = InterParamToProf.begin(); it != InterParamToProf.end(); ++it )
 426   {
 427     //ordered by param on hydr axis
 428     Handle(HYDROData_Profile) aProfile = it->second;
 429     const std::vector<gp_Pnt2d>& profile_points = GetProfileUZPoints(aProfile);
 430     profilesPoints.push_back(profile_points);
 431     if (profile_points.size() > maxNbPoints)
 432       maxNbPoints = profile_points.size();
 433   }
 434   std::set<double> paramsOnHAX;
 435
 436   if (InterParamToProf.size() <= 1)
 437   {
 438     if (warnings)
 439       warnings->push_back("Insufficient count of correct input profiles; abort");
 440     return;
 441   }
 442
 443   EquidParamsOnHAX(Hax2d, stepOnHA, InterParamToProf.begin()->first,
 444     InterParamToProf.rbegin()->first, paramsOnHAX);
 445
 446   //prepare ParamsPerSegm - number of interm.profiles for profile_points
 447   std::vector<std::vector<double>> ParamsPerSegm;
 448   std::set<double>::iterator it_params_hax = paramsOnHAX.begin();
 449   std::map<double, Handle(HYDROData_Profile)>::iterator it_p = InterParamToProf.begin();
 450   //
 451   std::map<double, std::pair<Handle(HYDROData_Profile), Handle(HYDROData_Profile)>> intermParamToProfPair; //for warnings only
 452   //
 453   it_p++;
 454   ParamsPerSegm.resize(profilesPoints.size()-1);
 455   for( int k=0; it_p != InterParamToProf.end(); )
 456   {
 457     if (*it_params_hax < it_p->first)
 458     {
 459       double val = *it_params_hax;
 460       ParamsPerSegm[k].push_back(val);
 461       it_params_hax++;
 462
 463       it_p--;
 464       const Handle(HYDROData_Profile)& cp = it_p->second;
 465       it_p++;
 466       const Handle(HYDROData_Profile)& np = it_p->second;
 467       std::pair<Handle(HYDROData_Profile), Handle(HYDROData_Profile)> profPair(cp, np);
 468       intermParamToProfPair[val] = profPair;
 469     }
 470     else
 471     {
 472       it_p++;
 473       k++;
 474     }
 475   }
 476   ParamsPerSegm.back().push_back(*paramsOnHAX.rbegin());
 477   //
 478   std::map<double, std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>> parToBankPoints;
 479   std::vector<double> paramHAXVec;
 480   paramHAXVec.reserve(paramsOnHAX.size());
 481   for( std::set<double>::iterator it = paramsOnHAX.begin(); it != paramsOnHAX.end(); ++it )
 482     paramHAXVec.push_back(*it);
 483   //
 484   Handle(Geom2d_Curve) LB2d, RB2d;
 485   PolyToCurve2d(LB, LB2d);
 486   PolyToCurve2d(RB, RB2d);
 487   GetPointsOnBanks( paramHAXVec, Hax2d, LB2d, RB2d, parToBankPoints, warnings, &intermParamToProfPair);
 488
 489   if (buildPresentationShapes)
 490   {
 491     //parToBankPoints.begin()->second.;
 492     //double lp1 = parToBankPoints.rbegin()->first;
 493
 494     //Handle(HYDROData_Profile) F_prof = InterParamToProf.begin()->second;
 495     //Handle(HYDROData_Profile) S_prof = InterParamToProf.rbegin()->second;
 496     BuildFace(parToBankPoints, LB2d, RB2d, prsDef);
 497   }
 498   //
 499
 500   if (estimateWarnOnly)
 501     return;
 502
 503   maxNbPoints = Max(pointsToInsert, maxNbPoints);
 504
 505   //insert points to profiles
 506   for (int i=0;i<profilesPoints.size();i++)
 507   {
 508     int nbPointsToInsert = maxNbPoints - profilesPoints[i].size();
 509     //insert nbPointsToInsert points to current profile
 510     InsertPoints(profilesPoints[i], nbPointsToInsert);
 511   }
 512
 513   for (int i=0;i<profilesPoints.size();i++)
 514   {
 515     bool toSwap = paramToSwapFlag[i];
 516     if (toSwap)
 517     {
 518       std::vector<gp_Pnt2d>& prof = profilesPoints[i];
 519       reversePoints(prof);
 520     }
 521   }
 522
 523   TopoDS_Compound cmp2d;
 524   if ( buildPresentationShapes)
 525   {
 526     BRep_Builder().MakeCompound(cmp2d);
 527     BRep_Builder().Add(cmp2d, prsDef.myLeftBank);
 528     BRep_Builder().Add(cmp2d, prsDef.myRightBank);
 529   }
 530
 531   for (int i=0;i<profilesPoints.size()-1;i++)
 532   {
 533     const std::vector<gp_Pnt2d>& prof1 = profilesPoints[i];
 534     const std::vector<gp_Pnt2d>& prof2 = profilesPoints[i+1];
 535     std::vector<std::vector<gp_Pnt2d>> IntermProf;
 536
 537     std::map<int, double> indToParam; //index of im.profile to param on hax. ()
 538     for (int j=0;j<ParamsPerSegm[i].size();j++)
 539     {
 540       double param = ParamsPerSegm[i][j];
 541       if (parToBankPoints.count(param) > 0) //not intersected with banks; skip this
 542         indToParam[j]=param;
 543     }
 544     int NbIntermProf = indToParam.size();
 545     IntermProf.resize(NbIntermProf);
 546
 547     for (int l=0;l<NbIntermProf;l++)
 548       IntermProf[l].resize(prof1.size());
 549     //calc interm. between prof1 & prof2
 550     assert (prof1.size() == prof2.size());
 551     for (int j=0;j<prof1.size();j++)
 552     {
 553       gp_Pnt2d P1 = prof1[j];
 554       gp_Pnt2d P2 = prof2[j];
 555       double X0 = P1.X();
 556       double Y0 = P1.Y();
 557       double X1 = P2.X();
 558       double Y1 = P2.Y();
 559       std::map<int, double>::iterator it = indToParam.begin();
 560       for ( int m=0; it != indToParam.end(); it++, m++ )
 561       {
 562         int ind = it->first;
 563         double t = ((double)ind+1)/((double)NbIntermProf+1);
 564         double Xp = X0 + (X1-X0)*t;
 565         double Yp = Y0 + (Y1-Y0)*t;
 566         gp_Pnt2d P(Xp, Yp);
 567         IntermProf[m][j] = P;
 568       }
 569     }
 570     std::map<int, double>::iterator it = indToParam.begin();
 571     for ( int m=0; it != indToParam.end(); it++, m++ )
 572     {
 573       const std::vector<gp_Pnt2d>& im_prof = IntermProf[m];
 574       double param = it->second;
 575       const std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>& BB = parToBankPoints[param]; //param is included in map; no checks
 576       gp_Pnt2d LP = BB.first.Value();
 577       gp_Pnt2d RP = BB.second.Value();
 578       HYDROData_ProfileUZ::PointsList pl;
 579       for (int k=0;k<im_prof.size();k++)
 580         pl.Append(im_prof[k].XY());
 581       HYDROData_Profile::ProfilePoints profile_points3d = HYDROData_Profile::CalculateProfilePoints(pl, LP.XY(), RP.XY());
 582       for (int k=1;k<=profile_points3d.Length();k++)
 583       {
 584         HYDROData_Bathymetry::AltitudePoint AP(profile_points3d(k).X(), profile_points3d(k).Y(), profile_points3d(k).Z());
 585         outBathypoints.push_back(AP);
 586       }
 587       if (buildPresentationShapes)
 588       {
 589         if (m==0 || i == profilesPoints.size()-2 && param == indToParam.rbegin()->second) //if last prof1/prof2 segment => take a last profile too
 590         {
 591           BRepLib_MakePolygon PM;
 592           for (int k=1;k<=profile_points3d.Length();k++)
 593             PM.Add(gp_Pnt(profile_points3d(k)));
 594           if (PM.IsDone())
 595           {
 596             const TopoDS_Wire& resW = PM.Wire();
 597             BRep_Builder().Add(cmp2d, resW);
 598           }
 599         }
 600       }
 601     }
 602   }
 603   if (buildPresentationShapes)
 604   {
 605     prsDef.myPrs3D = cmp2d;
 606   }
 607   //outBathy->SetAltitudePoints(bathypoints);
 608 }
 609
 610
 611