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 <stdexcept>
  46 #include <BRepBuilderAPI_MakeEdge2d.hxx>
  47 #include <BRepBuilderAPI_MakeEdge.hxx>
  48 #include <BRep_Tool.hxx>
  49 #include <BRepBuilderAPI_MakeWire.hxx>
  50 #include <TopoDS_Face.hxx>
  51 #include <gp_Pln.hxx>
  52 #include <BRepBuilderAPI_MakeFace.hxx>
  53 #include <gp.hxx>
  54 #include <gp_Ax2.hxx>
  55 #include <Geom_Plane.hxx>
  56 #include <TopoDS_Wire.hxx>
  57 #include <BRepLib.hxx>
  58 #include <BRepLib_MakePolygon.hxx>
  59 #include <TopoDS_Compound.hxx>
  60 #include <BRep_Builder.hxx>
  61
  62 //#define _DEVDEBUG_
  63 #include "HYDRO_trace.hxx"
  64
  65 static std::vector<gp_Pnt2d> GetProfileUZPoints(const Handle(HYDROData_Profile)& theProfile)
  66 {
  67   Handle(HYDROData_ProfileUZ) aProfileUZ = theProfile->GetProfileUZ( false );
  68   HYDROData_ProfileUZ::PointsList aSectPointsList = aProfileUZ->GetPoints();
  69   std::vector<gp_Pnt2d> points;
  70   points.reserve(aSectPointsList.Size());
  71   for ( int i = 1; i <= aSectPointsList.Size(); i++ )
  72   {
  73     const HYDROData_ProfileUZ::Point& aSectPoint = aSectPointsList(i);
  74     points.push_back( gp_Pnt2d(aSectPoint.X(), aSectPoint.Y()));
  75   }
  76   return points;
  77 }
  78
  79 static void GetMaxDist(
  80     const std::vector<gp_Pnt2d>& PNTS,
  81     double& dmax,
  82     double& dmax2,
  83     int& imax)
  84 {
  85   dmax = 0;
  86   dmax2 = 0;
  87   imax = -1;
  88   std::vector<double> dist;
  89   dist.reserve(PNTS.size()-1);
  90   for (int i=0; i<PNTS.size()-1; i++)
  91   {
  92     //dist.push_back(PNTS[i].Distance(PNTS[i+1])); //via distance between p1/p2
  93     dist.push_back(std::abs(PNTS[i].X()-PNTS[i+1].X())); //via U-diff of p1/p2
  94   }
  95   for (int i=0; i<dist.size(); i++)
  96   {
  97     double d = dist[i];
  98     if (d > dmax)
  99     {
 100       dmax = d;
 101       imax = i;
 102     }
 103   }
 104   for (int i=0; i<dist.size(); i++)
 105   {
 106     double d = dist[i];
 107     if (d > dmax2 && d < dmax)
 108       dmax2 = d;
 109   }
 110 }
 111
 112 static void InsertPoints(std::vector<gp_Pnt2d>& points,    // existing points of the profile, completed
 113                          int nbpoints)                     // number of points to add
 114 {
 115   points.reserve(points.size() + nbpoints);
 116   while (nbpoints>0)
 117   {
 118     double dmax=0, dmax2=0;
 119     int imax=-1;
 120     GetMaxDist(points, dmax, dmax2, imax);
 121     int nbPins = 0;
 122     if (dmax2 != 0)
 123       nbPins = floor(dmax/dmax2);
 124     else
 125       nbPins = nbpoints; //one segment, put all points here
 126     //insert nbPins points in [imax, imax+1] segment
 127     gp_Pnt2d p1 = points[imax];
 128     gp_Pnt2d p2 = points[imax+1];
 129     double X0 = p1.X();
 130     double Y0 = p1.Y();
 131     double X1 = p2.X();
 132     double Y1 = p2.Y();
 133     if (nbPins>nbpoints) //deny to add more points than needed
 134       nbPins = nbpoints;
 135     for (int i=0;i<nbPins;i++)
 136     {
 137       double t = ((double)i+1)/((double)nbPins+1);
 138       double Xp = X0 + (X1-X0)*t;
 139       double Yp = Y0 + (Y1-Y0)*t;
 140       gp_Pnt2d P(Xp, Yp);
 141       points.insert(points.begin() + imax + 1, P);
 142     }
 143     nbpoints-=nbPins;
 144   }
 145   //assert (nbpoints == 0);
 146 }
 147
 148 static void PolyToCurve2d(const Handle(HYDROData_PolylineXY)& poly, Handle(Geom2d_Curve)& c2d)
 149 {
 150   if (poly.IsNull())
 151     return;
 152   TopoDS_Shape sh = poly->GetShape();
 153   if (sh.ShapeType() == TopAbs_COMPOUND)
 154   {
 155     TopoDS_Iterator it(sh);
 156     if (it.More())
 157       sh = it.Value();
 158   }
 159
 160   if (sh.IsNull())
 161     return;
 162
 163   TopoDS_Edge ee;
 164   if (sh.ShapeType() == TopAbs_EDGE)
 165     ee = TopoDS::Edge(sh);
 166   else if (sh.ShapeType() == TopAbs_WIRE)
 167     ee = BRepAlgo::ConcatenateWireC0(TopoDS::Wire(sh)); //convert to bspline
 168
 169   if (ee.IsNull())
 170     return;
 171
 172   BRepAdaptor_Curve Ad(ee);
 173
 174   c2d = HYDROData_Tool::BRepAdaptorTo2DCurve(Ad);
 175 }
 176
 177 static void InterProfilesAndHAX(const HYDROData_SequenceOfObjects& profiles,
 178   const Handle(Geom2d_Curve)& Hax2d,
 179   std::map<double, Handle(HYDROData_Profile)>& profToInterParam,
 180   std::map<double, gp_Pnt2d>& interParamToPnt,
 181   std::map<double, gp_Dir2d>& interParamToVec,
 182   std::map<double, double>& interParamToLong,
 183   std::vector<bool>& paramToSwapFlag,
 184   std::vector<std::string>* warnings)
 185 {
 186   for (int i=1;i<=profiles.Size();i++)
 187   {
 188     Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast(profiles(i));
 189     gp_XY LP, RP;
 190     aProfile->GetLeftPoint( LP, false, true );
 191     aProfile->GetRightPoint( RP, false, true );
 192     gp_Pnt2d P1(LP), P2(RP);
 193     double d = P2.Distance(P1);
 194     if (d < gp::Resolution())
 195     {
 196       if (warnings)
 197         warnings->push_back(aProfile->GetName().toStdString() + " is skipped: Left and Right points is the same");
 198       continue;
 199     }
 200     Handle(Geom2d_Line) lin2d = new Geom2d_Line(P1,gp_Dir2d(P2.XY()-P1.XY()));
 201     gp_Dir2d prof_dir2d = lin2d->Direction();
 202
 203     Geom2dAdaptor_Curve linAd(lin2d, 0, d);
 204     Geom2dAdaptor_Curve haxAd(Hax2d);
 205     Geom2dInt_GInter isec( linAd, haxAd, 1.0e-6, 1.0e-6);
 206     if (!isec.IsDone())
 207     {
 208       if (warnings)
 209         warnings->push_back(aProfile->GetName().toStdString() + " is skipped: intersection between hydraulic axis & profile is failed");
 210       continue;
 211     }
 212     if (isec.NbPoints() == 0)
 213     {
 214       if (warnings)
 215         warnings->push_back(aProfile->GetName().toStdString() + " is skipped: no intersection between hydraulic axis & profile");
 216       continue;
 217     }
 218     if (isec.NbPoints() > 1)
 219     {
 220       if (warnings)
 221         warnings->push_back(aProfile->GetName().toStdString() + " is skipped; intersection between hydraulic axis & profile produces more than 1 point");
 222       continue;
 223     }
 224     double param = isec.Point(1).ParamOnSecond();
 225     gp_Pnt2d int_p2d;
 226     gp_Vec2d vec_on_hax_intersec;
 227     haxAd.D1(param, int_p2d, vec_on_hax_intersec);
 228     gp_Dir2d d_hax(vec_on_hax_intersec);
 229     double ang = d_hax.Angle(prof_dir2d);
 230     if (ang>0)
 231       paramToSwapFlag.push_back(true);
 232     else
 233       paramToSwapFlag.push_back(false);
 234     profToInterParam[param] = aProfile;
 235     interParamToPnt[param] = int_p2d;
 236     interParamToVec[param] = prof_dir2d;
 237     interParamToLong[param] = d;
 238   }
 239 }
 240
 241 static bool EquidParamsOnHAX(const Handle(Geom2d_Curve)& Hax2d, double step, double fp, double lp, std::set<double>& params)
 242 {
 243   Geom2dAdaptor_Curve ad(Hax2d);
 244   GCPnts_UniformAbscissa GUA(ad, step, fp, lp);
 245   if (!GUA.IsDone())
 246     return false;
 247
 248   for(int i = 1; i<= GUA.NbPoints(); i++)
 249     params.insert(GUA.Parameter(i));
 250
 251   return true;
 252 }
 253
 254 static void GetPointsOnBanks(
 255     const std::vector<double>& gua_params,                                                // params on hydr.axis for intermediate profiles
 256     const Handle(Geom2d_Curve)& Hax2d,                                                    // hydr.axis
 257     const Handle(Geom2d_Curve)& LB2d,                                                     // left bank
 258     const Handle(Geom2d_Curve)& RB2d,                                                     // right bank
 259     const std::map<double, std::pair<Handle(HYDROData_Profile),
 260                                      Handle(HYDROData_Profile)>>* intermParamToProfPair,  // param --> <profile before, profile after>
 261     std::map<double, std::pair<IntRes2d_IntersectionPoint,
 262                                IntRes2d_IntersectionPoint>>& parToBankPoints,             // return param on hydr.axis --> <intersect.pt left, intersect.pt right>
 263     std::vector<std::string>* warnings)                                                   // return warnings
 264 {
 265   for(int i = 0; i < gua_params.size(); i++)
 266   {
 267     double par = gua_params[i];
 268     gp_Pnt2d p;
 269     gp_Vec2d v;
 270     Hax2d->D1(par, p, v);                                    // get point and tangent on hydr.axis (what if non-diff?)
 271     //gp_Pnt2d outpLB, outpRB;
 272     gp_Dir2d n2d(v.X(), v.Y());                              // ignore Z (Z==0)
 273     n2d.Rotate(M_PI/2);                                      // normal to hydr.axis
 274     gp_Pnt2d p2d(gp_Pnt2d(p.X(), p.Y()));
 275     Handle(Geom2d_Line) nLine2d = new Geom2d_Line(p2d, n2d); // nLine2d is an orthogonal line to hydr. axis
 276
 277     // --- intersect nLine2d with banks
 278
 279     Geom2dAPI_InterCurveCurve intL(LB2d, nLine2d);
 280     Geom2dAPI_InterCurveCurve intR(RB2d, nLine2d);
 281
 282     if (warnings && intL.NbPoints()==0 )
 283     {
 284       std::string coord = (QString::number(p2d.X()) + ", "+QString::number(p2d.Y())).toStdString();
 285       if (intermParamToProfPair && intermParamToProfPair->count(par) != 0)
 286       {
 287         std::string fpn = intermParamToProfPair->at(par).first->GetName().toStdString();
 288         std::string spn = intermParamToProfPair->at(par).second->GetName().toStdString();
 289         warnings->push_back("no intersection between intermediate profile (point on h.axis:("+
 290           coord +"), located between " + fpn + " & " + spn + ") and left bank found; skipped");
 291       }
 292       else
 293         warnings->push_back("no intersection between intermediate profile (point on h.axis:("+ coord +")) and left bank found; skipped");
 294     }
 295
 296     if (warnings && intR.NbPoints()==0)
 297     {
 298       std::string coord = (QString::number(p2d.X()) + ", "+QString::number(p2d.Y())).toStdString();
 299       if (intermParamToProfPair && intermParamToProfPair->count(par) != 0)
 300       {
 301         std::string fpn = intermParamToProfPair->at(par).first->GetName().toStdString();
 302         std::string spn = intermParamToProfPair->at(par).second->GetName().toStdString();
 303         warnings->push_back("no intersection between intermediate profile (point on h.axis:("+
 304           coord +"), located between " + fpn + " & " + spn + ") and right bank found; skipped");
 305       }
 306       else
 307         warnings->push_back("no intersection between intermediate profile (point on h.axis:("+ coord +")) and right bank found; skipped");
 308     }
 309
 310     if ( intL.NbPoints()==0 || intR.NbPoints()==0)
 311       continue;
 312
 313     IntRes2d_IntersectionPoint aNearSolL, aNearSolR;
 314     double min_sq_dist = Precision::Infinite();
 315     for (int j=1;j<=intL.NbPoints();j++)
 316     {
 317       double sq_dist = p2d.SquareDistance(intL.Point(j));
 318       if (min_sq_dist > sq_dist)
 319       {
 320         min_sq_dist = sq_dist;
 321         aNearSolL = intL.Intersector().Point(j);
 322       }
 323     }
 324     min_sq_dist = Precision::Infinite();
 325     for (int j=1;j<=intR.NbPoints();j++)
 326     {
 327       double sq_dist = p2d.SquareDistance(intR.Point(j));
 328       if (min_sq_dist > sq_dist)
 329       {
 330         min_sq_dist = sq_dist;
 331         aNearSolR =  intR.Intersector().Point(j);
 332       }
 333     }
 334
 335     std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint> int_pair(aNearSolL, aNearSolR);
 336     parToBankPoints[par]=int_pair;
 337   }
 338 }
 339
 340 static void GetPointsOnBanks2(
 341     const std::map<double, Handle(HYDROData_Profile)> &InterParamToProf, // hydr.axis param --> profile
 342     const std::vector<std::vector<double>> &ParamsPerSegm,               // between profiles, list of intermediate param on hydr. axis: used here for number of profiles
 343     const Handle(Geom2d_Curve) &Hax2d,                                   // hydr.axis
 344     const Handle(Geom2d_Curve) &LB2d,                                    // left bank
 345     const Handle(Geom2d_Curve) &RB2d,                                    // right bank
 346     std::map<double, std::pair<double, double>> &parToBankParams,        // return param on hydr.axis --> <param left, param right>
 347     std::vector<std::string> *warnings)                                  // return warnings
 348 {
 349     DEBTRACE("GetPointsOnBanks2");
 350     std::vector<Handle(HYDROData_Profile)> mainProfiles;
 351     std::vector<double> mainParams;
 352     mainProfiles.reserve(InterParamToProf.size());
 353     for (std::map<double, Handle(HYDROData_Profile)>::const_iterator it =
 354             InterParamToProf.begin(); it != InterParamToProf.end(); ++it)
 355     {
 356         mainProfiles.push_back(it->second);
 357         mainParams.push_back(it->first);
 358     }
 359
 360     std::vector<IntRes2d_IntersectionPoint> leftInters;
 361     std::vector<IntRes2d_IntersectionPoint> rightInters;
 362
 363     for (int i = 0; i < mainProfiles.size(); i++)
 364     {
 365         gp_XY lptxy;
 366         gp_XY rptxy;
 367         mainProfiles[i]->GetLeftPoint(lptxy);
 368         mainProfiles[i]->GetRightPoint(rptxy);
 369         gp_Pnt2d lpt(lptxy);
 370         gp_Pnt2d rpt(rptxy);
 371         gp_Vec2d vec(lpt, rpt);
 372         gp_Dir2d dir(vec);
 373         Handle(Geom2d_Line) nLine2d = new Geom2d_Line(lpt, dir);
 374
 375         // --- intersect nLine2d with banks
 376
 377         Geom2dAPI_InterCurveCurve intL(LB2d, nLine2d);
 378         Geom2dAPI_InterCurveCurve intR(RB2d, nLine2d);
 379
 380         if (warnings && intL.NbPoints() == 0)
 381         {
 382             std::string pn = mainProfiles[i]->GetName().toStdString();
 383             warnings->push_back(
 384                     "no intersection between profile " + pn
 385                             + " and left bank found; skipped");
 386         }
 387
 388         if (warnings && intR.NbPoints() == 0)
 389         {
 390             std::string pn = mainProfiles[i]->GetName().toStdString();
 391             warnings->push_back(
 392                     "no intersection between profile " + pn
 393                             + " and right bank found; skipped");
 394         }
 395
 396         if (intL.NbPoints() == 0 || intR.NbPoints() == 0)
 397             continue;
 398
 399         IntRes2d_IntersectionPoint aNearSolL, aNearSolR;
 400         double min_sq_dist = Precision::Infinite();
 401         for (int j = 1; j <= intL.NbPoints(); j++)
 402         {
 403             double sq_dist = lpt.SquareDistance(intL.Point(j));
 404             if (min_sq_dist > sq_dist)
 405             {
 406                 min_sq_dist = sq_dist;
 407                 aNearSolL = intL.Intersector().Point(j);
 408             }
 409         }
 410         min_sq_dist = Precision::Infinite();
 411         for (int j = 1; j <= intR.NbPoints(); j++)
 412         {
 413             double sq_dist = rpt.SquareDistance(intR.Point(j));
 414             if (min_sq_dist > sq_dist)
 415             {
 416                 min_sq_dist = sq_dist;
 417                 aNearSolR = intR.Intersector().Point(j);
 418             }
 419         }
 420         leftInters.push_back(aNearSolL);
 421         rightInters.push_back(aNearSolR);
 422     }
 423
 424     if (leftInters.size() != mainProfiles.size())
 425     {
 426         DEBTRACE("profiles skipped: " << mainProfiles.size() << " " << leftInters.size() << " " << rightInters.size());
 427         return;
 428     }
 429
 430     // --- for each segment between 2 profiles
 431     //     compute intermediate points regularly spaced on each bank
 432     //     store parameters in parToBankParams
 433
 434     for (int i = 0; i < ParamsPerSegm.size(); i++)
 435     {
 436         DEBTRACE("Segment " << i << " between profiles");
 437         int nbProfs = ParamsPerSegm[i].size();
 438         double fpL = leftInters[i].ParamOnFirst();
 439         double fpR = rightInters[i].ParamOnFirst();
 440         double lpL = leftInters[i + 1].ParamOnFirst();
 441         double lpR = rightInters[i + 1].ParamOnFirst();
 442         DEBTRACE(nbProfs << " " << fpL << " " << fpR << " " << lpL << " " << lpR);
 443         Geom2dAdaptor_Curve lad(LB2d);
 444         GCPnts_UniformAbscissa lgua;
 445         lgua.Initialize(lad, nbProfs, fpL, lpL);
 446         if (!lgua.IsDone())
 447         {
 448             DEBTRACE("Pb intermediate points on left bank")
 449             return;
 450         }
 451         std::vector<double> lparams;
 452         lparams.reserve(lgua.NbPoints());
 453         DEBTRACE("lgua.NbPoints() "<< lgua.NbPoints());
 454         for (int j = 1; j <= lgua.NbPoints(); j++)
 455             lparams[j - 1] = lgua.Parameter(j);
 456         Geom2dAdaptor_Curve rad(RB2d);
 457         GCPnts_UniformAbscissa rgua;
 458         rgua.Initialize(rad, nbProfs, fpR, lpR);
 459         if (!rgua.IsDone())
 460         {
 461             DEBTRACE("Pb intermediate points on right bank")
 462             return;
 463         }
 464         std::vector<double> rparams;
 465         rparams.reserve(rgua.NbPoints());
 466         DEBTRACE("rgua.NbPoints() "<< rgua.NbPoints());
 467         for (int j = 1; j <= rgua.NbPoints(); j++)
 468             rparams[j - 1] = rgua.Parameter(j);
 469         if (lparams.size() != rparams.size())
 470         {
 471             DEBTRACE("Number of intermediate points left and right unequal " << lparams.size() << " " << rparams.size());
 472             return;
 473         }
 474         for (int j = 0; j < nbProfs; j++)
 475         {
 476             std::pair<double, double> paramsPair(lparams[j], rparams[j]);
 477             parToBankParams[ParamsPerSegm[i][j]] = paramsPair;
 478             DEBTRACE("  parToBankParams " << ParamsPerSegm[i][j] << " -->(" << lparams[j] << ", "<< rparams[j] << ")");
 479         }
 480     }
 481 }
 482
 483
 484 void BuildFace(
 485     bool mode,                                                                      // true: intermediate point profiles based on hydraulic axis normals
 486     const std::map<double, std::pair<IntRes2d_IntersectionPoint,
 487                                      IntRes2d_IntersectionPoint>>& parToBankPoints, // param on hydr.axis --> <intersect.pt left, intersect.pt right>
 488     const std::map<double, std::pair<double, double>> parToBankParams,              // param on hydr.axis --> <param left, param right>
 489     const Handle(Geom2d_Curve)& LB2d,                                               // left bank
 490     const Handle(Geom2d_Curve)& RB2d,                                               // right bank
 491     HYDROData_Stream::PrsDefinition& prsDef)                                        // all the resulting shapes
 492 {
 493   prsDef.myPrs3D.Nullify();
 494   prsDef.myPrs2D.Nullify();
 495   prsDef.myLeftBank.Nullify();
 496   prsDef.myRightBank.Nullify();
 497   prsDef.myInlet.Nullify();
 498   prsDef.myOutlet.Nullify();
 499
 500   double lb_1, lb_2, rb_1, rb_2;
 501   if (mode)
 502   {
 503       lb_1 = parToBankPoints.begin()->second.first.ParamOnFirst();                      // param of first profile on left bank
 504       rb_1 = parToBankPoints.begin()->second.second.ParamOnFirst();                     // param of first profile on right bank
 505       lb_2 = parToBankPoints.rbegin()->second.first.ParamOnFirst();                     // param of last profile on left bank
 506       rb_2 = parToBankPoints.rbegin()->second.second.ParamOnFirst();                    // param of last profile on right bank
 507   }
 508   else
 509   {
 510       lb_1 = parToBankParams.begin()->second.first;                      // param of first profile on left bank
 511       rb_1 = parToBankParams.begin()->second.second;                     // param of first profile on right bank
 512       lb_2 = parToBankParams.rbegin()->second.first;                     // param of last profile on left bank
 513       rb_2 = parToBankParams.rbegin()->second.second;                    // param of last profile on right bank
 514   }
 515
 516   BRepBuilderAPI_MakeEdge2d LEM(LB2d, lb_1, lb_2);                                  // left bank limited --> edge
 517   if (!LEM.IsDone())
 518     return;
 519   BRepBuilderAPI_MakeEdge2d REM(RB2d, rb_1, rb_2);                                  // right bank limited --> edge
 520   if (!REM.IsDone())
 521     return;
 522   TopoDS_Edge LBE = LEM.Edge();                                                     // left edge
 523   TopoDS_Edge RBE = REM.Edge();                                                     // right edge
 524   if (LBE.IsNull() || RBE.IsNull())
 525     return;
 526   BRepBuilderAPI_MakeEdge PFEM(BRep_Tool::Pnt(LEM.Vertex1()), BRep_Tool::Pnt(REM.Vertex1()));
 527   BRepBuilderAPI_MakeEdge PLEM(BRep_Tool::Pnt(LEM.Vertex2()), BRep_Tool::Pnt(REM.Vertex2()));
 528   if (!PFEM.IsDone())
 529     return;
 530   if (!PLEM.IsDone())
 531     return;
 532   TopoDS_Edge FProfE = PFEM.Edge();                                                 // edge first section
 533   TopoDS_Edge SProfE = PLEM.Edge();                                                 // edge last section
 534   BRepBuilderAPI_MakeWire WM(FProfE, LBE, SProfE, RBE);                             // wire with 4 edges
 535   if (WM.IsDone())
 536   {
 537     TopoDS_Wire W = WM.Wire();
 538     BRepLib::BuildCurves3d(W);
 539     BRepBuilderAPI_MakeFace FM(Geom_Plane(gp::XOY()).Pln(), WM.Wire());             // face on wire
 540     if (FM.IsDone())
 541     {
 542       prsDef.myPrs2D = FM.Face();
 543       prsDef.myLeftBank = LBE;
 544       prsDef.myRightBank = RBE;
 545       prsDef.myInlet = FProfE;
 546       prsDef.myOutlet = SProfE;
 547     }
 548   }
 549 }
 550
 551 void reversePoints(std::vector<gp_Pnt2d>& points)
 552 {
 553   std::reverse(points.begin(), points.end());
 554 }
 555
 556 void HYDROData_StreamLinearInterpolation::Perform(
 557     const HYDROData_SequenceOfObjects &profiles,
 558     int pointsToInsert,
 559     double stepOnHA,
 560     const Handle(HYDROData_PolylineXY) &hax,
 561     const Handle(HYDROData_PolylineXY) &LB,
 562     const Handle(HYDROData_PolylineXY) &RB,
 563     HYDROData_Bathymetry::AltitudePoints &outBathypoints,
 564     bool buildPresentationShapes,
 565     bool estimateWarnOnly,
 566     HYDROData_Stream::PrsDefinition &prsDef,
 567     std::vector<std::string> *warnings)
 568 {
 569     bool isHaxNormal = false; // build intermediate point profiles based on hydraulic axis normals or on regular points on left and right banks
 570     DEBTRACE("Perform, mode " << isHaxNormal);
 571
 572     if (hax.IsNull() || LB.IsNull() || RB.IsNull())
 573         return;
 574
 575     Handle(Geom2d_Curve) Hax2d;
 576     PolyToCurve2d(hax, Hax2d);
 577     if (hax->IsClosed() || Hax2d->IsClosed()) //can't be closed
 578     {
 579         if (warnings)
 580             warnings->push_back(
 581                     hax->GetName().toStdString() + " is closed; abort");
 582         return;
 583     }
 584
 585     // --- intersect profiles with given hydr.axis
 586
 587     std::map<double, Handle(HYDROData_Profile)> InterParamToProf; // hydr.axis param --> profile
 588     std::map<double, gp_Pnt2d> InterParamToPnt;                   // hydr.axis param --> intersection point
 589     std::map<double, gp_Dir2d> InterParamToVec;                   // hydr.axis param --> profile direction
 590     std::map<double, double> InterParamToLong;                    // hydr.axis param --> profile size
 591     std::vector<bool> paramToSwapFlag;
 592     // --- profiles ordered and oriented by hydraulic axis: order = param on hydr.axis
 593     InterProfilesAndHAX(profiles, Hax2d, InterParamToProf, InterParamToPnt,
 594             InterParamToVec, InterParamToLong, paramToSwapFlag, warnings);
 595
 596     // --- vectors of points (U, Z) for profiles
 597
 598     std::vector<std::vector<gp_Pnt2d>> profilesPoints;            // vector of vector of points per profile (U, Z)
 599     profilesPoints.reserve(InterParamToProf.size());
 600     int maxNbPoints = 0;
 601     std::vector<double> profileParamsHax;                         // vector of hydr.axis params for main profiles
 602     profileParamsHax.reserve(InterParamToProf.size());
 603     for (std::map<double, Handle(HYDROData_Profile)>::iterator it =
 604             InterParamToProf.begin(); it != InterParamToProf.end(); ++it)
 605     {
 606         profileParamsHax.push_back(it->first);
 607         Handle(HYDROData_Profile) aProfile = it->second;
 608         const std::vector<gp_Pnt2d> &profile_points = GetProfileUZPoints(
 609                 aProfile);
 610         profilesPoints.push_back(profile_points);
 611         if (profile_points.size() > maxNbPoints)
 612             maxNbPoints = profile_points.size();
 613     }
 614
 615     if (InterParamToProf.size() <= 1)
 616     {
 617         if (warnings)
 618             warnings->push_back(
 619                     "Insufficient count of correct input profiles; abort");
 620         return;
 621     }
 622
 623     // --- ordered params on hydr.axis for intermediate profile_points
 624
 625     std::set<double> paramsOnHAX;                                 // ordered params on hydr.axis for intermediate profile_points
 626     EquidParamsOnHAX(Hax2d, stepOnHA, InterParamToProf.begin()->first,
 627             InterParamToProf.rbegin()->first, paramsOnHAX);
 628
 629     // --- prepare ParamsPerSegm - number of interm.profiles for profile_points between 2 profiles
 630
 631     std::vector<std::vector<double>> ParamsPerSegm;
 632     std::set<double>::iterator it_params_hax = paramsOnHAX.begin();
 633     std::map<double, Handle(HYDROData_Profile)>::iterator it_p =
 634             InterParamToProf.begin();
 635     std::map<double, std::pair<Handle(HYDROData_Profile),
 636     Handle(HYDROData_Profile)>> intermParamToProfPair; //for warnings only
 637
 638     it_p++;
 639     ParamsPerSegm.resize(profilesPoints.size() - 1);
 640     for (int k = 0; it_p != InterParamToProf.end();)
 641     {
 642         if (*it_params_hax < it_p->first)          // intermediate profile_point before end of segment
 643         {
 644             double val = *it_params_hax;
 645             ParamsPerSegm[k].push_back(val);       // param for intermediate profile_point
 646             it_params_hax++;
 647
 648             it_p--;
 649             const Handle(HYDROData_Profile) &cp = it_p->second;
 650             it_p++;
 651             const Handle(HYDROData_Profile) &np = it_p->second;
 652             std::pair<Handle(HYDROData_Profile), Handle(HYDROData_Profile)> profPair(
 653                     cp, np);
 654             intermParamToProfPair[val] = profPair; // param of prof.point --> <begin prof., end prof.>
 655         }
 656         else
 657         {
 658             it_p++;
 659             k++;
 660         }
 661     }
 662     ParamsPerSegm.back().push_back(*paramsOnHAX.rbegin());
 663
 664     // --- get points on banks
 665
 666     Handle(Geom2d_Curve) LB2d, RB2d;
 667     PolyToCurve2d(LB, LB2d);
 668     PolyToCurve2d(RB, RB2d);
 669
 670     std::vector<double> paramHAXVec;               // replace set by vector, for params on hydr.axis for intermediate profile_points
 671     std::map<double,
 672             std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>> parToBankPoints;
 673
 674     std::map<double, std::pair<double, double>>  parToBankParams;
 675
 676     if (isHaxNormal)
 677     {
 678         paramHAXVec.reserve(paramsOnHAX.size());
 679         for (std::set<double>::iterator it = paramsOnHAX.begin();
 680                 it != paramsOnHAX.end(); ++it)
 681             paramHAXVec.push_back(*it);
 682         GetPointsOnBanks(paramHAXVec, Hax2d, LB2d, RB2d,&intermParamToProfPair,
 683                          parToBankPoints, warnings);                         // points on banks for extremities of intermediates profiles normal to axis
 684     }
 685     else
 686     {
 687         GetPointsOnBanks2(InterParamToProf, ParamsPerSegm, Hax2d, LB2d, RB2d,
 688                           parToBankParams, warnings);                         // points on banks regularly spaced, for extremities of intermediates profiles
 689     }
 690
 691     if (buildPresentationShapes)
 692     {
 693         BuildFace(isHaxNormal, parToBankPoints, parToBankParams, LB2d, RB2d, prsDef);
 694     }
 695
 696     if (estimateWarnOnly)
 697         return;
 698
 699     maxNbPoints = Max(pointsToInsert, maxNbPoints);
 700
 701     // --- insert points to profiles
 702
 703     for (int i = 0; i < profilesPoints.size(); i++)
 704     {
 705         int nbPointsToInsert = maxNbPoints - profilesPoints[i].size();
 706         InsertPoints(profilesPoints[i], nbPointsToInsert);
 707     }
 708
 709     for (int i = 0; i < profilesPoints.size(); i++)
 710     {
 711         bool toSwap = paramToSwapFlag[i];
 712         if (toSwap)
 713         {
 714             std::vector<gp_Pnt2d> &prof = profilesPoints[i];
 715             reversePoints(prof);
 716         }
 717     }
 718
 719     TopoDS_Compound cmp2d;
 720     if (buildPresentationShapes)
 721     {
 722         BRep_Builder().MakeCompound(cmp2d);
 723         BRep_Builder().Add(cmp2d, prsDef.myLeftBank);
 724         BRep_Builder().Add(cmp2d, prsDef.myRightBank);
 725     }
 726
 727     for (int i = 0; i < profilesPoints.size() - 1; i++)
 728     {
 729         const std::vector<gp_Pnt2d> &prof1 = profilesPoints[i];      // profile points (U, Z), begin of segment
 730         const std::vector<gp_Pnt2d> &prof2 = profilesPoints[i + 1];  // profile points (U, Z), end of segment
 731         if (prof1.size() != prof2.size())
 732             throw std::range_error(
 733                     "internal error: profiles are not adjusted with the same number of points");
 734
 735         std::map<int, double> indToParam;                            //index of intermediate profile to param on hydr.axis
 736         for (int j = 0; j < ParamsPerSegm[i].size(); j++)
 737         {
 738             double param = ParamsPerSegm[i][j];
 739             if (isHaxNormal)
 740             {
 741                 if (parToBankPoints.count(param) > 0)                         // keep only when intermediate profiles normal to axis intersected with banks
 742                     indToParam[j] = param;
 743             }
 744             else
 745                 indToParam[j] = param;                                        // always keep when when intermediate profiles based on left and right banks nodes
 746         }
 747         int NbIntermProf = indToParam.size();
 748         std::vector<std::vector<gp_Pnt2d>> IntermProf;                        // vector of intermediate profiles points to compute
 749         IntermProf.resize(NbIntermProf);
 750
 751         // --- calculate intermediate profile (U, Z) between prof1 & prof2
 752
 753         for (int l = 0; l < NbIntermProf; l++)
 754             IntermProf[l].resize(prof1.size());
 755         for (int j = 0; j < prof1.size(); j++)
 756         {
 757             gp_Pnt2d P1 = prof1[j];                                           // point (u,z) of rank j in first profile of segment
 758             gp_Pnt2d P2 = prof2[j];                                           // point (u,z) of rank j in last profile of segment
 759             double X0 = P1.X();
 760             double Y0 = P1.Y();
 761             double X1 = P2.X();
 762             double Y1 = P2.Y();
 763             std::map<int, double>::iterator it = indToParam.begin();
 764             for (int m = 0; it != indToParam.end(); it++, m++)                // iterate on intermediate profiles to build points of rank j
 765             {
 766                 int ind = it->first;
 767                 double t = ((double) ind + 1) / ((double) NbIntermProf + 1);
 768                 double Xp = X0 + (X1 - X0) * t;                               // linear interpolation on u
 769                 double Yp = Y0 + (Y1 - Y0) * t;                               // linear interpolation on z
 770                 gp_Pnt2d P(Xp, Yp);
 771                 IntermProf[m][j] = P;                                         // interpolated point (u,z)
 772             }
 773         }
 774         std::map<int, double>::iterator it = indToParam.begin();
 775         for (int m = 0; it != indToParam.end(); it++, m++)                    // iterate on intermediate profiles
 776         {
 777             const std::vector<gp_Pnt2d> &im_prof = IntermProf[m];
 778             double param = it->second;
 779             gp_Pnt2d LP;
 780             gp_Pnt2d RP;
 781             if (isHaxNormal)
 782             {
 783                 const std::pair<IntRes2d_IntersectionPoint,
 784                         IntRes2d_IntersectionPoint> &BB = parToBankPoints[param]; // intersection points with banks
 785                 LP = BB.first.Value();                                            // left bank point
 786                 RP = BB.second.Value();                                           // right bank point
 787             }
 788             else
 789             {
 790                 double lparam = parToBankParams[param].first;
 791                 double rparam = parToBankParams[param].second;
 792                 LB2d->D0(lparam, LP);
 793                 RB2d->D0(rparam, RP);
 794             }
 795             HYDROData_ProfileUZ::PointsList pl;
 796             for (int k = 0; k < im_prof.size(); k++)
 797                 pl.Append(im_prof[k].XY());
 798             HYDROData_Profile::ProfilePoints profile_points3d =
 799                     HYDROData_Profile::CalculateProfilePoints(pl, LP.XY(),
 800                             RP.XY());                                         // intermediate profile geolocalization 3D
 801             DEBTRACE(" ---------------- fill outBathypoints");
 802             for (int k = 1; k <= profile_points3d.Length(); k++)
 803             {
 804                 HYDROData_Bathymetry::AltitudePoint AP(profile_points3d(k).X(),
 805                         profile_points3d(k).Y(), profile_points3d(k).Z());
 806                 outBathypoints.push_back(AP);
 807             }
 808             if (buildPresentationShapes)
 809             {
 810                 if (m == 0
 811                         || i == profilesPoints.size() - 2
 812                                 && param == indToParam.rbegin()->second)      //if last prof1/prof2 segment => take a last profile too
 813                 {
 814                     BRepLib_MakePolygon PM;
 815                     for (int k = 1; k <= profile_points3d.Length(); k++)
 816                         PM.Add(gp_Pnt(profile_points3d(k)));
 817                     if (PM.IsDone())
 818                     {
 819                         const TopoDS_Wire &resW = PM.Wire();
 820                         BRep_Builder().Add(cmp2d, resW);
 821                     }
 822                 }
 823             }
 824         }
 825     }
 826     if (buildPresentationShapes)
 827     {
 828         prsDef.myPrs3D = cmp2d;
 829     }
 830 }
 831
 832
 833