From 8c841070095678b9857dcbfb053e1b77486b3f48 Mon Sep 17 00:00:00 2001 From: Paul RASCLE Date: Fri, 20 Nov 2020 18:35:15 +0100 Subject: [PATCH] lot 4: improve stream linear interpolation, first working version --- doc/salome/examples/profilsStreamLin.sx | 45 ++ doc/salome/examples/profilsStreamLin2.sx | 45 ++ .../HYDROData_StreamLinearInterpolation.cxx | 687 ++++++++++++------ 3 files changed, 544 insertions(+), 233 deletions(-) create mode 100644 doc/salome/examples/profilsStreamLin.sx create mode 100644 doc/salome/examples/profilsStreamLin2.sx diff --git a/doc/salome/examples/profilsStreamLin.sx b/doc/salome/examples/profilsStreamLin.sx new file mode 100644 index 00000000..0316d28a --- /dev/null +++ b/doc/salome/examples/profilsStreamLin.sx @@ -0,0 +1,45 @@ +C Generated by HYDRO Module +C +B P +CN Profile_1 +CP 0 0 +CP 20.000 5.000 10.000 15.000 +CP 2 + 0.000 0.000 11.000 + 2.000 0.000 9.000 + 4.000 0.000 9.000 + 5.000 0.000 10.000 + 6.000 0.000 10.000 + 7.000 0.000 6.000 + 8.000 0.000 6.000 + 10.000 0.000 8.000 + 13.000 0.000 8.000 + 14.000 0.000 9.000 +B P +CN Profile_2 +CP 0 0 +CP 25.000 7.000 25.000 27.000 +CP 2 + 0.000 0.000 11.000 + 3.000 0.000 7.000 + 5.000 0.000 6.000 + 8.000 0.000 7.000 + 12.000 0.000 6.000 + 13.000 0.000 6.000 + 15.000 0.000 5.000 + 20.000 0.000 10.000 +B P +CN Profile_3 +CP 0 0 +CP 30.000 5.000 40.000 15.000 +CP 2 + 0.000 0.000 10.000 + 1.000 0.000 6.000 + 3.000 0.000 6.000 + 4.000 0.000 5.000 + 6.000 0.000 5.000 + 7.000 0.000 8.000 + 9.000 0.000 8.000 + 10.000 0.000 7.000 + 12.000 0.000 7.000 + 14.000 0.000 11.000 diff --git a/doc/salome/examples/profilsStreamLin2.sx b/doc/salome/examples/profilsStreamLin2.sx new file mode 100644 index 00000000..085880c8 --- /dev/null +++ b/doc/salome/examples/profilsStreamLin2.sx @@ -0,0 +1,45 @@ +C Generated by HYDRO Module +C +B P +CN Profile_1 +CP 0 0 +CP 20.000 5.000 6.000 5.000 +CP 2 + 0.000 0.000 11.000 + 2.000 0.000 9.000 + 4.000 0.000 9.000 + 5.000 0.000 10.000 + 6.000 0.000 10.000 + 7.000 0.000 6.000 + 8.000 0.000 6.000 + 10.000 0.000 8.000 + 13.000 0.000 8.000 + 14.000 0.000 9.000 +B P +CN Profile_2 +CP 0 0 +CP 25.000 7.000 25.000 27.000 +CP 2 + 0.000 0.000 11.000 + 3.000 0.000 7.000 + 5.000 0.000 6.000 + 8.000 0.000 7.000 + 12.000 0.000 6.000 + 13.000 0.000 6.000 + 15.000 0.000 5.000 + 20.000 0.000 10.000 +B P +CN Profile_3 +CP 0 0 +CP 30.000 5.000 44.000 5.000 +CP 2 + 0.000 0.000 10.000 + 1.000 0.000 6.000 + 3.000 0.000 6.000 + 4.000 0.000 5.000 + 6.000 0.000 5.000 + 7.000 0.000 8.000 + 9.000 0.000 8.000 + 10.000 0.000 7.000 + 12.000 0.000 7.000 + 14.000 0.000 11.000 diff --git a/src/HYDROData/HYDROData_StreamLinearInterpolation.cxx b/src/HYDROData/HYDROData_StreamLinearInterpolation.cxx index 23eab97e..758fa277 100644 --- a/src/HYDROData/HYDROData_StreamLinearInterpolation.cxx +++ b/src/HYDROData/HYDROData_StreamLinearInterpolation.cxx @@ -42,6 +42,7 @@ #include #include #include +#include #include #include #include @@ -58,6 +59,9 @@ #include #include +#define _DEVDEBUG_ +#include "HYDRO_trace.hxx" + static std::vector GetProfileUZPoints(const Handle(HYDROData_Profile)& theProfile) { Handle(HYDROData_ProfileUZ) aProfileUZ = theProfile->GetProfileUZ( false ); @@ -72,10 +76,11 @@ static std::vector GetProfileUZPoints(const Handle(HYDROData_Profile)& return points; } -static void GetMaxDist(const std::vector& PNTS, - double& dmax, - double& dmax2, - int& imax) +static void GetMaxDist( + const std::vector& PNTS, + double& dmax, + double& dmax2, + int& imax) { dmax = 0; dmax2 = 0; @@ -104,7 +109,8 @@ static void GetMaxDist(const std::vector& PNTS, } } -static void InsertPoints(std::vector& points, int nbpoints) +static void InsertPoints(std::vector& points, // existing points of the profile, completed + int nbpoints) // number of points to add { points.reserve(points.size() + nbpoints); while (nbpoints>0) @@ -171,6 +177,9 @@ static void PolyToCurve2d(const Handle(HYDROData_PolylineXY)& poly, Handle(Geom2 static void InterProfilesAndHAX(const HYDROData_SequenceOfObjects& profiles, const Handle(Geom2d_Curve)& Hax2d, std::map& profToInterParam, + std::map& interParamToPnt, + std::map& interParamToVec, + std::map& interParamToLong, std::vector& paramToSwapFlag, std::vector* warnings) { @@ -223,6 +232,9 @@ static void InterProfilesAndHAX(const HYDROData_SequenceOfObjects& profiles, else paramToSwapFlag.push_back(false); profToInterParam[param] = aProfile; + interParamToPnt[param] = int_p2d; + interParamToVec[param] = prof_dir2d; + interParamToLong[param] = d; } } @@ -239,28 +251,31 @@ static bool EquidParamsOnHAX(const Handle(Geom2d_Curve)& Hax2d, double step, dou return true; } -static void GetPointsOnBanks(const std::vector& gua_params, - const Handle(Geom2d_Curve)& Hax2d, - const Handle(Geom2d_Curve)& LB2d, - const Handle(Geom2d_Curve)& RB2d, - std::map>& parToBankPoints, - std::vector* warnings, - const std::map>* intermParamToProfPair) +static void GetPointsOnBanks( + const std::vector& gua_params, // params on hydr.axis for intermediate profiles + const Handle(Geom2d_Curve)& Hax2d, // hydr.axis + const Handle(Geom2d_Curve)& LB2d, // left bank + const Handle(Geom2d_Curve)& RB2d, // right bank + const std::map>* intermParamToProfPair, // param --> + std::map>& parToBankPoints, // return param on hydr.axis --> + std::vector* warnings) // return warnings { for(int i = 0; i < gua_params.size(); i++) { double par = gua_params[i]; gp_Pnt2d p; gp_Vec2d v; - Hax2d->D1(par, p, v); //if non-diff? - gp_Pnt2d outpLB, outpRB; - - gp_Dir2d n2d(v.X(), v.Y()); //ignore Z (Z==0) - n2d.Rotate(M_PI/2); - //nLine2d is an orthogonal line to hydr. axis + Hax2d->D1(par, p, v); // get point and tangent on hydr.axis (what if non-diff?) + //gp_Pnt2d outpLB, outpRB; + gp_Dir2d n2d(v.X(), v.Y()); // ignore Z (Z==0) + n2d.Rotate(M_PI/2); // normal to hydr.axis gp_Pnt2d p2d(gp_Pnt2d(p.X(), p.Y())); - Handle(Geom2d_Line) nLine2d = new Geom2d_Line(p2d, n2d); - //intersect nLine2d with banks + Handle(Geom2d_Line) nLine2d = new Geom2d_Line(p2d, n2d); // nLine2d is an orthogonal line to hydr. axis + + // --- intersect nLine2d with banks + Geom2dAPI_InterCurveCurve intL(LB2d, nLine2d); Geom2dAPI_InterCurveCurve intR(RB2d, nLine2d); @@ -303,7 +318,6 @@ static void GetPointsOnBanks(const std::vector& gua_params, if (min_sq_dist > sq_dist) { min_sq_dist = sq_dist; - //aNearSolL = intL.Point(j); aNearSolL = intL.Intersector().Point(j); } } @@ -314,7 +328,6 @@ static void GetPointsOnBanks(const std::vector& gua_params, if (min_sq_dist > sq_dist) { min_sq_dist = sq_dist; - //aNearSolR = intR.Point(j); aNearSolR = intR.Intersector().Point(j); } } @@ -324,11 +337,158 @@ static void GetPointsOnBanks(const std::vector& gua_params, } } +static void GetPointsOnBanks2( + const std::map &InterParamToProf, // hydr.axis param --> profile + const std::vector> &ParamsPerSegm, // between profiles, list of intermediate param on hydr. axis: used here for number of profiles + const Handle(Geom2d_Curve) &Hax2d, // hydr.axis + const Handle(Geom2d_Curve) &LB2d, // left bank + const Handle(Geom2d_Curve) &RB2d, // right bank + std::map> &parToBankParams, // return param on hydr.axis --> + std::vector *warnings) // return warnings +{ + DEBTRACE("GetPointsOnBanks2"); + std::vector mainProfiles; + std::vector mainParams; + mainProfiles.reserve(InterParamToProf.size()); + for (std::map::const_iterator it = + InterParamToProf.begin(); it != InterParamToProf.end(); ++it) + { + mainProfiles.push_back(it->second); + mainParams.push_back(it->first); + } + + std::vector leftInters; + std::vector rightInters; + + for (int i = 0; i < mainProfiles.size(); i++) + { + gp_XY lptxy; + gp_XY rptxy; + mainProfiles[i]->GetLeftPoint(lptxy); + mainProfiles[i]->GetRightPoint(rptxy); + gp_Pnt2d lpt(lptxy); + gp_Pnt2d rpt(rptxy); + gp_Vec2d vec(lpt, rpt); + gp_Dir2d dir(vec); + Handle(Geom2d_Line) nLine2d = new Geom2d_Line(lpt, dir); + + // --- intersect nLine2d with banks + + Geom2dAPI_InterCurveCurve intL(LB2d, nLine2d); + Geom2dAPI_InterCurveCurve intR(RB2d, nLine2d); + + if (warnings && intL.NbPoints() == 0) + { + std::string pn = mainProfiles[i]->GetName().toStdString(); + warnings->push_back( + "no intersection between profile " + pn + + " and left bank found; skipped"); + } + + if (warnings && intR.NbPoints() == 0) + { + std::string pn = mainProfiles[i]->GetName().toStdString(); + warnings->push_back( + "no intersection between profile " + pn + + " and right bank found; skipped"); + } + + if (intL.NbPoints() == 0 || intR.NbPoints() == 0) + continue; + + IntRes2d_IntersectionPoint aNearSolL, aNearSolR; + double min_sq_dist = Precision::Infinite(); + for (int j = 1; j <= intL.NbPoints(); j++) + { + double sq_dist = lpt.SquareDistance(intL.Point(j)); + if (min_sq_dist > sq_dist) + { + min_sq_dist = sq_dist; + aNearSolL = intL.Intersector().Point(j); + } + } + min_sq_dist = Precision::Infinite(); + for (int j = 1; j <= intR.NbPoints(); j++) + { + double sq_dist = rpt.SquareDistance(intR.Point(j)); + if (min_sq_dist > sq_dist) + { + min_sq_dist = sq_dist; + aNearSolR = intR.Intersector().Point(j); + } + } + leftInters.push_back(aNearSolL); + rightInters.push_back(aNearSolR); + } + + if (leftInters.size() != mainProfiles.size()) + { + DEBTRACE("profiles skipped: " << mainProfiles.size() << " " << leftInters.size() << " " << rightInters.size()); + return; + } + + // --- for each segment between 2 profiles + // compute intermediate points regularly spaced on each bank + // store parameters in parToBankParams + + for (int i = 0; i < ParamsPerSegm.size(); i++) + { + DEBTRACE("Segment " << i << " between profiles"); + int nbProfs = ParamsPerSegm[i].size(); + double fpL = leftInters[i].ParamOnFirst(); + double fpR = rightInters[i].ParamOnFirst(); + double lpL = leftInters[i + 1].ParamOnFirst(); + double lpR = rightInters[i + 1].ParamOnFirst(); + DEBTRACE(nbProfs << " " << fpL << " " << fpR << " " << lpL << " " << lpR); + Geom2dAdaptor_Curve lad(LB2d); + GCPnts_UniformAbscissa lgua; + lgua.Initialize(lad, nbProfs, fpL, lpL); + if (!lgua.IsDone()) + { + DEBTRACE("Pb intermediate points on left bank") + return; + } + std::vector lparams; + lparams.reserve(lgua.NbPoints()); + DEBTRACE("lgua.NbPoints() "<< lgua.NbPoints()); + for (int j = 1; j <= lgua.NbPoints(); j++) + lparams[j - 1] = lgua.Parameter(j); + Geom2dAdaptor_Curve rad(RB2d); + GCPnts_UniformAbscissa rgua; + rgua.Initialize(rad, nbProfs, fpR, lpR); + if (!rgua.IsDone()) + { + DEBTRACE("Pb intermediate points on right bank") + return; + } + std::vector rparams; + rparams.reserve(rgua.NbPoints()); + DEBTRACE("rgua.NbPoints() "<< rgua.NbPoints()); + for (int j = 1; j <= rgua.NbPoints(); j++) + rparams[j - 1] = rgua.Parameter(j); + if (lparams.size() != rparams.size()) + { + DEBTRACE("Number of intermediate points left and right unequal " << lparams.size() << " " << rparams.size()); + return; + } + for (int j = 0; j < nbProfs; j++) + { + std::pair paramsPair(lparams[j], rparams[j]); + parToBankParams[ParamsPerSegm[i][j]] = paramsPair; + DEBTRACE(" parToBankParams " << ParamsPerSegm[i][j] << " -->(" << lparams[j] << ", "<< rparams[j] << ")"); + } + } +} + -void BuildFace(const std::map>& parToBankPoints, - const Handle(Geom2d_Curve)& LB2d, - const Handle(Geom2d_Curve)& RB2d, - HYDROData_Stream::PrsDefinition& prsDef) +void BuildFace( + bool mode, // true: intermediate point profiles based on hydraulic axis normals + const std::map>& parToBankPoints, // param on hydr.axis --> + const std::map> parToBankParams, // param on hydr.axis --> + const Handle(Geom2d_Curve)& LB2d, // left bank + const Handle(Geom2d_Curve)& RB2d, // right bank + HYDROData_Stream::PrsDefinition& prsDef) // all the resulting shapes { prsDef.myPrs3D.Nullify(); prsDef.myPrs2D.Nullify(); @@ -338,23 +498,29 @@ void BuildFace(const std::mapsecond.first.ParamOnFirst(); - rb_1 = parToBankPoints.begin()->second.second.ParamOnFirst(); - - lb_2 = parToBankPoints.rbegin()->second.first.ParamOnFirst(); - rb_2 = parToBankPoints.rbegin()->second.second.ParamOnFirst(); - - //stat = GetItersectionParam(F_prof, LB2d, lb_1) && GetItersectionParam(F_prof, RB2d, rb_1) - // && GetItersectionParam(S_prof, LB2d, lb_2) && GetItersectionParam(S_prof, RB2d, rb_2); + if (mode) + { + lb_1 = parToBankPoints.begin()->second.first.ParamOnFirst(); // param of first profile on left bank + rb_1 = parToBankPoints.begin()->second.second.ParamOnFirst(); // param of first profile on right bank + lb_2 = parToBankPoints.rbegin()->second.first.ParamOnFirst(); // param of last profile on left bank + rb_2 = parToBankPoints.rbegin()->second.second.ParamOnFirst(); // param of last profile on right bank + } + else + { + lb_1 = parToBankParams.begin()->second.first; // param of first profile on left bank + rb_1 = parToBankParams.begin()->second.second; // param of first profile on right bank + lb_2 = parToBankParams.rbegin()->second.first; // param of last profile on left bank + rb_2 = parToBankParams.rbegin()->second.second; // param of last profile on right bank + } - BRepBuilderAPI_MakeEdge2d LEM(LB2d, lb_1, lb_2); + BRepBuilderAPI_MakeEdge2d LEM(LB2d, lb_1, lb_2); // left bank limited --> edge if (!LEM.IsDone()) return; - BRepBuilderAPI_MakeEdge2d REM(RB2d, rb_1, rb_2); + BRepBuilderAPI_MakeEdge2d REM(RB2d, rb_1, rb_2); // right bank limited --> edge if (!REM.IsDone()) return; - TopoDS_Edge LBE = LEM.Edge(); - TopoDS_Edge RBE = REM.Edge(); + TopoDS_Edge LBE = LEM.Edge(); // left edge + TopoDS_Edge RBE = REM.Edge(); // right edge if (LBE.IsNull() || RBE.IsNull()) return; BRepBuilderAPI_MakeEdge PFEM(BRep_Tool::Pnt(LEM.Vertex1()), BRep_Tool::Pnt(REM.Vertex1())); @@ -363,14 +529,14 @@ void BuildFace(const std::map& points) std::reverse(points.begin(), points.end()); } -void HYDROData_StreamLinearInterpolation::Perform(const HYDROData_SequenceOfObjects& profiles, - int pointsToInsert, - double stepOnHA, - const Handle(HYDROData_PolylineXY)& hax, - const Handle(HYDROData_PolylineXY)& LB, - const Handle(HYDROData_PolylineXY)& RB, - HYDROData_Bathymetry::AltitudePoints& outBathypoints, - bool buildPresentationShapes, - bool estimateWarnOnly, - HYDROData_Stream::PrsDefinition& prsDef, - std::vector* warnings) +void HYDROData_StreamLinearInterpolation::Perform( + const HYDROData_SequenceOfObjects &profiles, + int pointsToInsert, + double stepOnHA, + const Handle(HYDROData_PolylineXY) &hax, + const Handle(HYDROData_PolylineXY) &LB, + const Handle(HYDROData_PolylineXY) &RB, + HYDROData_Bathymetry::AltitudePoints &outBathypoints, + bool buildPresentationShapes, + bool estimateWarnOnly, + HYDROData_Stream::PrsDefinition &prsDef, + std::vector *warnings) { - //intersect profiles with given hydr.axis - if (hax.IsNull() || LB.IsNull() || RB.IsNull()) - return; + bool isHaxNormal = false; // build intermediate point profiles based on hydraulic axis normals or on regular points on left and right banks + DEBTRACE("Perform, mode " << isHaxNormal); - Handle(Geom2d_Curve) Hax2d; - PolyToCurve2d(hax, Hax2d); - if (hax->IsClosed() || Hax2d->IsClosed()) //can't be closed - { - if (warnings) - warnings->push_back(hax->GetName().toStdString() + " is closed; abort"); - return; - } + if (hax.IsNull() || LB.IsNull() || RB.IsNull()) + return; - std::map InterParamToProf; - std::vector paramToSwapFlag; - //profToInterParam is output param map: profile to intersection param on hydr.axis - InterProfilesAndHAX(profiles, Hax2d, InterParamToProf, paramToSwapFlag, warnings); - // - std::vector> profilesPoints; - profilesPoints.reserve(InterParamToProf.size()); - int maxNbPoints = 0; - - //for (int i=1;i<=InterParamToProf.Extent();i++) - for( std::map::iterator it = InterParamToProf.begin(); it != InterParamToProf.end(); ++it ) - { - //ordered by param on hydr axis - Handle(HYDROData_Profile) aProfile = it->second; - const std::vector& profile_points = GetProfileUZPoints(aProfile); - profilesPoints.push_back(profile_points); - if (profile_points.size() > maxNbPoints) - maxNbPoints = profile_points.size(); - } - std::set paramsOnHAX; + Handle(Geom2d_Curve) Hax2d; + PolyToCurve2d(hax, Hax2d); + if (hax->IsClosed() || Hax2d->IsClosed()) //can't be closed + { + if (warnings) + warnings->push_back( + hax->GetName().toStdString() + " is closed; abort"); + return; + } - if (InterParamToProf.size() <= 1) - { - if (warnings) - warnings->push_back("Insufficient count of correct input profiles; abort"); - return; - } + // --- intersect profiles with given hydr.axis + + std::map InterParamToProf; // hydr.axis param --> profile + std::map InterParamToPnt; // hydr.axis param --> intersection point + std::map InterParamToVec; // hydr.axis param --> profile direction + std::map InterParamToLong; // hydr.axis param --> profile size + std::vector paramToSwapFlag; + // --- profiles ordered and oriented by hydraulic axis: order = param on hydr.axis + InterProfilesAndHAX(profiles, Hax2d, InterParamToProf, InterParamToPnt, + InterParamToVec, InterParamToLong, paramToSwapFlag, warnings); + + // --- vectors of points (U, Z) for profiles + + std::vector> profilesPoints; // vector of vector of points per profile (U, Z) + profilesPoints.reserve(InterParamToProf.size()); + int maxNbPoints = 0; + std::vector profileParamsHax; // vector of hydr.axis params for main profiles + profileParamsHax.reserve(InterParamToProf.size()); + for (std::map::iterator it = + InterParamToProf.begin(); it != InterParamToProf.end(); ++it) + { + profileParamsHax.push_back(it->first); + Handle(HYDROData_Profile) aProfile = it->second; + const std::vector &profile_points = GetProfileUZPoints( + aProfile); + profilesPoints.push_back(profile_points); + if (profile_points.size() > maxNbPoints) + maxNbPoints = profile_points.size(); + } - EquidParamsOnHAX(Hax2d, stepOnHA, InterParamToProf.begin()->first, - InterParamToProf.rbegin()->first, paramsOnHAX); - - //prepare ParamsPerSegm - number of interm.profiles for profile_points - std::vector> ParamsPerSegm; - std::set::iterator it_params_hax = paramsOnHAX.begin(); - std::map::iterator it_p = InterParamToProf.begin(); - // - std::map> intermParamToProfPair; //for warnings only - // - it_p++; - ParamsPerSegm.resize(profilesPoints.size()-1); - for( int k=0; it_p != InterParamToProf.end(); ) - { - if (*it_params_hax < it_p->first) + if (InterParamToProf.size() <= 1) { - double val = *it_params_hax; - ParamsPerSegm[k].push_back(val); - it_params_hax++; - - it_p--; - const Handle(HYDROData_Profile)& cp = it_p->second; - it_p++; - const Handle(HYDROData_Profile)& np = it_p->second; - std::pair profPair(cp, np); - intermParamToProfPair[val] = profPair; + if (warnings) + warnings->push_back( + "Insufficient count of correct input profiles; abort"); + return; } - else + + // --- ordered params on hydr.axis for intermediate profile_points + + std::set paramsOnHAX; // ordered params on hydr.axis for intermediate profile_points + EquidParamsOnHAX(Hax2d, stepOnHA, InterParamToProf.begin()->first, + InterParamToProf.rbegin()->first, paramsOnHAX); + + // --- prepare ParamsPerSegm - number of interm.profiles for profile_points between 2 profiles + + std::vector> ParamsPerSegm; + std::set::iterator it_params_hax = paramsOnHAX.begin(); + std::map::iterator it_p = + InterParamToProf.begin(); + std::map> intermParamToProfPair; //for warnings only + + it_p++; + ParamsPerSegm.resize(profilesPoints.size() - 1); + for (int k = 0; it_p != InterParamToProf.end();) { - it_p++; - k++; + if (*it_params_hax < it_p->first) // intermediate profile_point before end of segment + { + double val = *it_params_hax; + ParamsPerSegm[k].push_back(val); // param for intermediate profile_point + it_params_hax++; + + it_p--; + const Handle(HYDROData_Profile) &cp = it_p->second; + it_p++; + const Handle(HYDROData_Profile) &np = it_p->second; + std::pair profPair( + cp, np); + intermParamToProfPair[val] = profPair; // param of prof.point --> + } + else + { + it_p++; + k++; + } } - } - ParamsPerSegm.back().push_back(*paramsOnHAX.rbegin()); - // - std::map> parToBankPoints; - std::vector paramHAXVec; - paramHAXVec.reserve(paramsOnHAX.size()); - for( std::set::iterator it = paramsOnHAX.begin(); it != paramsOnHAX.end(); ++it ) - paramHAXVec.push_back(*it); - // - Handle(Geom2d_Curve) LB2d, RB2d; - PolyToCurve2d(LB, LB2d); - PolyToCurve2d(RB, RB2d); - GetPointsOnBanks( paramHAXVec, Hax2d, LB2d, RB2d, parToBankPoints, warnings, &intermParamToProfPair); - - if (buildPresentationShapes) - { - //parToBankPoints.begin()->second.; - //double lp1 = parToBankPoints.rbegin()->first; + ParamsPerSegm.back().push_back(*paramsOnHAX.rbegin()); - //Handle(HYDROData_Profile) F_prof = InterParamToProf.begin()->second; - //Handle(HYDROData_Profile) S_prof = InterParamToProf.rbegin()->second; - BuildFace(parToBankPoints, LB2d, RB2d, prsDef); - } - // + // --- get points on banks - if (estimateWarnOnly) - return; + Handle(Geom2d_Curve) LB2d, RB2d; + PolyToCurve2d(LB, LB2d); + PolyToCurve2d(RB, RB2d); - maxNbPoints = Max(pointsToInsert, maxNbPoints); + std::vector paramHAXVec; // replace set by vector, for params on hydr.axis for intermediate profile_points + std::map> parToBankPoints; - //insert points to profiles - for (int i=0;i> parToBankParams; - for (int i=0;i& prof = profilesPoints[i]; - reversePoints(prof); + paramHAXVec.reserve(paramsOnHAX.size()); + for (std::set::iterator it = paramsOnHAX.begin(); + it != paramsOnHAX.end(); ++it) + paramHAXVec.push_back(*it); + GetPointsOnBanks(paramHAXVec, Hax2d, LB2d, RB2d,&intermParamToProfPair, + parToBankPoints, warnings); // points on banks for extremities of intermediates profiles normal to axis + } + else + { + GetPointsOnBanks2(InterParamToProf, ParamsPerSegm, Hax2d, LB2d, RB2d, + parToBankParams, warnings); // points on banks regularly spaced, for extremities of intermediates profiles } - } - TopoDS_Compound cmp2d; - if ( buildPresentationShapes) - { - BRep_Builder().MakeCompound(cmp2d); - BRep_Builder().Add(cmp2d, prsDef.myLeftBank); - BRep_Builder().Add(cmp2d, prsDef.myRightBank); - } + if (buildPresentationShapes) + { + BuildFace(isHaxNormal, parToBankPoints, parToBankParams, LB2d, RB2d, prsDef); + } - for (int i=0;i& prof1 = profilesPoints[i]; - const std::vector& prof2 = profilesPoints[i+1]; - std::vector> IntermProf; + if (estimateWarnOnly) + return; + + maxNbPoints = Max(pointsToInsert, maxNbPoints); + + // --- insert points to profiles - std::map indToParam; //index of im.profile to param on hax. () - for (int j=0;j 0) //not intersected with banks; skip this - indToParam[j]=param; + int nbPointsToInsert = maxNbPoints - profilesPoints[i].size(); + InsertPoints(profilesPoints[i], nbPointsToInsert); } - int NbIntermProf = indToParam.size(); - IntermProf.resize(NbIntermProf); - - for (int l=0;l::iterator it = indToParam.begin(); - for ( int m=0; it != indToParam.end(); it++, m++ ) - { - int ind = it->first; - double t = ((double)ind+1)/((double)NbIntermProf+1); - double Xp = X0 + (X1-X0)*t; - double Yp = Y0 + (Y1-Y0)*t; - gp_Pnt2d P(Xp, Yp); - IntermProf[m][j] = P; - } + bool toSwap = paramToSwapFlag[i]; + if (toSwap) + { + std::vector &prof = profilesPoints[i]; + reversePoints(prof); + } } - std::map::iterator it = indToParam.begin(); - for ( int m=0; it != indToParam.end(); it++, m++ ) + + TopoDS_Compound cmp2d; + if (buildPresentationShapes) { - const std::vector& im_prof = IntermProf[m]; - double param = it->second; - const std::pair& BB = parToBankPoints[param]; //param is included in map; no checks - gp_Pnt2d LP = BB.first.Value(); - gp_Pnt2d RP = BB.second.Value(); - HYDROData_ProfileUZ::PointsList pl; - for (int k=0;ksecond) //if last prof1/prof2 segment => take a last profile too + BRep_Builder().MakeCompound(cmp2d); + BRep_Builder().Add(cmp2d, prsDef.myLeftBank); + BRep_Builder().Add(cmp2d, prsDef.myRightBank); + } + + for (int i = 0; i < profilesPoints.size() - 1; i++) + { + const std::vector &prof1 = profilesPoints[i]; // profile points (U, Z), begin of segment + const std::vector &prof2 = profilesPoints[i + 1]; // profile points (U, Z), end of segment + if (prof1.size() != prof2.size()) + throw std::range_error( + "internal error: profiles are not adjusted with the same number of points"); + + std::map indToParam; //index of intermediate profile to param on hydr.axis + for (int j = 0; j < ParamsPerSegm[i].size(); j++) { - BRepLib_MakePolygon PM; - for (int k=1;k<=profile_points3d.Length();k++) - PM.Add(gp_Pnt(profile_points3d(k))); - if (PM.IsDone()) - { - const TopoDS_Wire& resW = PM.Wire(); - BRep_Builder().Add(cmp2d, resW); - } + double param = ParamsPerSegm[i][j]; + if (isHaxNormal) + { + if (parToBankPoints.count(param) > 0) // keep only when intermediate profiles normal to axis intersected with banks + indToParam[j] = param; + } + else + indToParam[j] = param; // always keep when when intermediate profiles based on left and right banks nodes + } + int NbIntermProf = indToParam.size(); + std::vector> IntermProf; // vector of intermediate profiles points to compute + IntermProf.resize(NbIntermProf); + + // --- calculate intermediate profile (U, Z) between prof1 & prof2 + + for (int l = 0; l < NbIntermProf; l++) + IntermProf[l].resize(prof1.size()); + for (int j = 0; j < prof1.size(); j++) + { + gp_Pnt2d P1 = prof1[j]; // point (u,z) of rank j in first profile of segment + gp_Pnt2d P2 = prof2[j]; // point (u,z) of rank j in last profile of segment + double X0 = P1.X(); + double Y0 = P1.Y(); + double X1 = P2.X(); + double Y1 = P2.Y(); + std::map::iterator it = indToParam.begin(); + for (int m = 0; it != indToParam.end(); it++, m++) // iterate on intermediate profiles to build points of rank j + { + int ind = it->first; + double t = ((double) ind + 1) / ((double) NbIntermProf + 1); + double Xp = X0 + (X1 - X0) * t; // linear interpolation on u + double Yp = Y0 + (Y1 - Y0) * t; // linear interpolation on z + gp_Pnt2d P(Xp, Yp); + IntermProf[m][j] = P; // interpolated point (u,z) + } + } + std::map::iterator it = indToParam.begin(); + for (int m = 0; it != indToParam.end(); it++, m++) // iterate on intermediate profiles + { + const std::vector &im_prof = IntermProf[m]; + double param = it->second; + gp_Pnt2d LP; + gp_Pnt2d RP; + if (isHaxNormal) + { + const std::pair &BB = parToBankPoints[param]; // intersection points with banks + LP = BB.first.Value(); // left bank point + RP = BB.second.Value(); // right bank point + } + else + { + double lparam = parToBankParams[param].first; + double rparam = parToBankParams[param].second; + LB2d->D0(lparam, LP); + RB2d->D0(rparam, RP); + } + HYDROData_ProfileUZ::PointsList pl; + for (int k = 0; k < im_prof.size(); k++) + pl.Append(im_prof[k].XY()); + HYDROData_Profile::ProfilePoints profile_points3d = + HYDROData_Profile::CalculateProfilePoints(pl, LP.XY(), + RP.XY()); // intermediate profile geolocalization 3D + for (int k = 1; k <= profile_points3d.Length(); k++) + { + HYDROData_Bathymetry::AltitudePoint AP(profile_points3d(k).X(), + profile_points3d(k).Y(), profile_points3d(k).Z()); + outBathypoints.push_back(AP); + } + if (buildPresentationShapes) + { + if (m == 0 + || i == profilesPoints.size() - 2 + && param == indToParam.rbegin()->second) //if last prof1/prof2 segment => take a last profile too + { + BRepLib_MakePolygon PM; + for (int k = 1; k <= profile_points3d.Length(); k++) + PM.Add(gp_Pnt(profile_points3d(k))); + if (PM.IsDone()) + { + const TopoDS_Wire &resW = PM.Wire(); + BRep_Builder().Add(cmp2d, resW); + } + } + } } - } } - } - if (buildPresentationShapes) - { - prsDef.myPrs3D = cmp2d; - } - //outBathy->SetAltitudePoints(bathypoints); + if (buildPresentationShapes) + { + prsDef.myPrs3D = cmp2d; + } } -- 2.39.2