#include <vector>
#include <map>
#include <set>
+#include <stdexcept>
#include <BRepBuilderAPI_MakeEdge2d.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <BRep_Tool.hxx>
#include <TopoDS_Compound.hxx>
#include <BRep_Builder.hxx>
+#define _DEVDEBUG_
+#include "HYDRO_trace.hxx"
+
static std::vector<gp_Pnt2d> GetProfileUZPoints(const Handle(HYDROData_Profile)& theProfile)
{
Handle(HYDROData_ProfileUZ) aProfileUZ = theProfile->GetProfileUZ( false );
return points;
}
-static void GetMaxDist(const std::vector<gp_Pnt2d>& PNTS,
- double& dmax,
- double& dmax2,
- int& imax)
+static void GetMaxDist(
+ const std::vector<gp_Pnt2d>& PNTS,
+ double& dmax,
+ double& dmax2,
+ int& imax)
{
dmax = 0;
dmax2 = 0;
}
}
-static void InsertPoints(std::vector<gp_Pnt2d>& points, int nbpoints)
+static void InsertPoints(std::vector<gp_Pnt2d>& points, // existing points of the profile, completed
+ int nbpoints) // number of points to add
{
points.reserve(points.size() + nbpoints);
while (nbpoints>0)
static void InterProfilesAndHAX(const HYDROData_SequenceOfObjects& profiles,
const Handle(Geom2d_Curve)& Hax2d,
std::map<double, Handle(HYDROData_Profile)>& profToInterParam,
+ std::map<double, gp_Pnt2d>& interParamToPnt,
+ std::map<double, gp_Dir2d>& interParamToVec,
+ std::map<double, double>& interParamToLong,
std::vector<bool>& paramToSwapFlag,
std::vector<std::string>* warnings)
{
else
paramToSwapFlag.push_back(false);
profToInterParam[param] = aProfile;
+ interParamToPnt[param] = int_p2d;
+ interParamToVec[param] = prof_dir2d;
+ interParamToLong[param] = d;
}
}
return true;
}
-static void GetPointsOnBanks(const std::vector<double>& gua_params,
- const Handle(Geom2d_Curve)& Hax2d,
- const Handle(Geom2d_Curve)& LB2d,
- const Handle(Geom2d_Curve)& RB2d,
- std::map<double, std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>>& parToBankPoints,
- std::vector<std::string>* warnings,
- const std::map<double, std::pair<Handle(HYDROData_Profile), Handle(HYDROData_Profile)>>* intermParamToProfPair)
+static void GetPointsOnBanks(
+ const std::vector<double>& 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<double, std::pair<Handle(HYDROData_Profile),
+ Handle(HYDROData_Profile)>>* intermParamToProfPair, // param --> <profile before, profile after>
+ std::map<double, std::pair<IntRes2d_IntersectionPoint,
+ IntRes2d_IntersectionPoint>>& parToBankPoints, // return param on hydr.axis --> <intersect.pt left, intersect.pt right>
+ std::vector<std::string>* 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);
if (min_sq_dist > sq_dist)
{
min_sq_dist = sq_dist;
- //aNearSolL = intL.Point(j);
aNearSolL = intL.Intersector().Point(j);
}
}
if (min_sq_dist > sq_dist)
{
min_sq_dist = sq_dist;
- //aNearSolR = intR.Point(j);
aNearSolR = intR.Intersector().Point(j);
}
}
}
}
+static void GetPointsOnBanks2(
+ const std::map<double, Handle(HYDROData_Profile)> &InterParamToProf, // hydr.axis param --> profile
+ const std::vector<std::vector<double>> &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<double, std::pair<double, double>> &parToBankParams, // return param on hydr.axis --> <param left, param right>
+ std::vector<std::string> *warnings) // return warnings
+{
+ DEBTRACE("GetPointsOnBanks2");
+ std::vector<Handle(HYDROData_Profile)> mainProfiles;
+ std::vector<double> mainParams;
+ mainProfiles.reserve(InterParamToProf.size());
+ for (std::map<double, Handle(HYDROData_Profile)>::const_iterator it =
+ InterParamToProf.begin(); it != InterParamToProf.end(); ++it)
+ {
+ mainProfiles.push_back(it->second);
+ mainParams.push_back(it->first);
+ }
+
+ std::vector<IntRes2d_IntersectionPoint> leftInters;
+ std::vector<IntRes2d_IntersectionPoint> 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<double> 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<double> 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<double, double> paramsPair(lparams[j], rparams[j]);
+ parToBankParams[ParamsPerSegm[i][j]] = paramsPair;
+ DEBTRACE(" parToBankParams " << ParamsPerSegm[i][j] << " -->(" << lparams[j] << ", "<< rparams[j] << ")");
+ }
+ }
+}
+
-void BuildFace(const std::map<double, std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>>& 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<double, std::pair<IntRes2d_IntersectionPoint,
+ IntRes2d_IntersectionPoint>>& parToBankPoints, // param on hydr.axis --> <intersect.pt left, intersect.pt right>
+ const std::map<double, std::pair<double, double>> parToBankParams, // param on hydr.axis --> <param left, param right>
+ 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();
prsDef.myOutlet.Nullify();
double lb_1, lb_2, rb_1, rb_2;
- lb_1 = parToBankPoints.begin()->second.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()));
return;
if (!PLEM.IsDone())
return;
- TopoDS_Edge FProfE = PFEM.Edge();
- TopoDS_Edge SProfE = PLEM.Edge();
- BRepBuilderAPI_MakeWire WM(FProfE, LBE, SProfE, RBE);
+ TopoDS_Edge FProfE = PFEM.Edge(); // edge first section
+ TopoDS_Edge SProfE = PLEM.Edge(); // edge last section
+ BRepBuilderAPI_MakeWire WM(FProfE, LBE, SProfE, RBE); // wire with 4 edges
if (WM.IsDone())
{
TopoDS_Wire W = WM.Wire();
BRepLib::BuildCurves3d(W);
- BRepBuilderAPI_MakeFace FM(Geom_Plane(gp::XOY()).Pln(), WM.Wire());
+ BRepBuilderAPI_MakeFace FM(Geom_Plane(gp::XOY()).Pln(), WM.Wire()); // face on wire
if (FM.IsDone())
{
prsDef.myPrs2D = FM.Face();
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<std::string>* 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<std::string> *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<double, Handle(HYDROData_Profile)> InterParamToProf;
- std::vector<bool> paramToSwapFlag;
- //profToInterParam is output param map: profile to intersection param on hydr.axis
- InterProfilesAndHAX(profiles, Hax2d, InterParamToProf, paramToSwapFlag, warnings);
- //
- std::vector<std::vector<gp_Pnt2d>> profilesPoints;
- profilesPoints.reserve(InterParamToProf.size());
- int maxNbPoints = 0;
-
- //for (int i=1;i<=InterParamToProf.Extent();i++)
- for( std::map<double, Handle(HYDROData_Profile)>::iterator it = InterParamToProf.begin(); it != InterParamToProf.end(); ++it )
- {
- //ordered by param on hydr axis
- Handle(HYDROData_Profile) aProfile = it->second;
- const std::vector<gp_Pnt2d>& profile_points = GetProfileUZPoints(aProfile);
- profilesPoints.push_back(profile_points);
- if (profile_points.size() > maxNbPoints)
- maxNbPoints = profile_points.size();
- }
- std::set<double> 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<double, Handle(HYDROData_Profile)> InterParamToProf; // hydr.axis param --> profile
+ std::map<double, gp_Pnt2d> InterParamToPnt; // hydr.axis param --> intersection point
+ std::map<double, gp_Dir2d> InterParamToVec; // hydr.axis param --> profile direction
+ std::map<double, double> InterParamToLong; // hydr.axis param --> profile size
+ std::vector<bool> 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<std::vector<gp_Pnt2d>> profilesPoints; // vector of vector of points per profile (U, Z)
+ profilesPoints.reserve(InterParamToProf.size());
+ int maxNbPoints = 0;
+ std::vector<double> profileParamsHax; // vector of hydr.axis params for main profiles
+ profileParamsHax.reserve(InterParamToProf.size());
+ for (std::map<double, Handle(HYDROData_Profile)>::iterator it =
+ InterParamToProf.begin(); it != InterParamToProf.end(); ++it)
+ {
+ profileParamsHax.push_back(it->first);
+ Handle(HYDROData_Profile) aProfile = it->second;
+ const std::vector<gp_Pnt2d> &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<std::vector<double>> ParamsPerSegm;
- std::set<double>::iterator it_params_hax = paramsOnHAX.begin();
- std::map<double, Handle(HYDROData_Profile)>::iterator it_p = InterParamToProf.begin();
- //
- std::map<double, std::pair<Handle(HYDROData_Profile), Handle(HYDROData_Profile)>> 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<Handle(HYDROData_Profile), Handle(HYDROData_Profile)> 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<double> 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<std::vector<double>> ParamsPerSegm;
+ std::set<double>::iterator it_params_hax = paramsOnHAX.begin();
+ std::map<double, Handle(HYDROData_Profile)>::iterator it_p =
+ InterParamToProf.begin();
+ std::map<double, std::pair<Handle(HYDROData_Profile),
+ Handle(HYDROData_Profile)>> 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<Handle(HYDROData_Profile), Handle(HYDROData_Profile)> profPair(
+ cp, np);
+ intermParamToProfPair[val] = profPair; // param of prof.point --> <begin prof., end prof.>
+ }
+ else
+ {
+ it_p++;
+ k++;
+ }
}
- }
- ParamsPerSegm.back().push_back(*paramsOnHAX.rbegin());
- //
- std::map<double, std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>> parToBankPoints;
- std::vector<double> paramHAXVec;
- paramHAXVec.reserve(paramsOnHAX.size());
- for( std::set<double>::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<double> paramHAXVec; // replace set by vector, for params on hydr.axis for intermediate profile_points
+ std::map<double,
+ std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>> parToBankPoints;
- //insert points to profiles
- for (int i=0;i<profilesPoints.size();i++)
- {
- int nbPointsToInsert = maxNbPoints - profilesPoints[i].size();
- //insert nbPointsToInsert points to current profile
- InsertPoints(profilesPoints[i], nbPointsToInsert);
- }
+ std::map<double, std::pair<double, double>> parToBankParams;
- for (int i=0;i<profilesPoints.size();i++)
- {
- bool toSwap = paramToSwapFlag[i];
- if (toSwap)
+ if (isHaxNormal)
{
- std::vector<gp_Pnt2d>& prof = profilesPoints[i];
- reversePoints(prof);
+ paramHAXVec.reserve(paramsOnHAX.size());
+ for (std::set<double>::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<profilesPoints.size()-1;i++)
- {
- const std::vector<gp_Pnt2d>& prof1 = profilesPoints[i];
- const std::vector<gp_Pnt2d>& prof2 = profilesPoints[i+1];
- std::vector<std::vector<gp_Pnt2d>> IntermProf;
+ if (estimateWarnOnly)
+ return;
+
+ maxNbPoints = Max(pointsToInsert, maxNbPoints);
+
+ // --- insert points to profiles
- std::map<int, double> indToParam; //index of im.profile to param on hax. ()
- for (int j=0;j<ParamsPerSegm[i].size();j++)
+ for (int i = 0; i < profilesPoints.size(); i++)
{
- double param = ParamsPerSegm[i][j];
- if (parToBankPoints.count(param) > 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<NbIntermProf;l++)
- IntermProf[l].resize(prof1.size());
- //calc interm. between prof1 & prof2
- assert (prof1.size() == prof2.size());
- for (int j=0;j<prof1.size();j++)
+
+ for (int i = 0; i < profilesPoints.size(); i++)
{
- gp_Pnt2d P1 = prof1[j];
- gp_Pnt2d P2 = prof2[j];
- double X0 = P1.X();
- double Y0 = P1.Y();
- double X1 = P2.X();
- double Y1 = P2.Y();
- std::map<int, double>::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<gp_Pnt2d> &prof = profilesPoints[i];
+ reversePoints(prof);
+ }
}
- std::map<int, double>::iterator it = indToParam.begin();
- for ( int m=0; it != indToParam.end(); it++, m++ )
+
+ TopoDS_Compound cmp2d;
+ if (buildPresentationShapes)
{
- const std::vector<gp_Pnt2d>& im_prof = IntermProf[m];
- double param = it->second;
- const std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>& 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;k<im_prof.size();k++)
- pl.Append(im_prof[k].XY());
- HYDROData_Profile::ProfilePoints profile_points3d = HYDROData_Profile::CalculateProfilePoints(pl, LP.XY(), RP.XY());
- 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
+ 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<gp_Pnt2d> &prof1 = profilesPoints[i]; // profile points (U, Z), begin of segment
+ const std::vector<gp_Pnt2d> &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<int, double> 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<std::vector<gp_Pnt2d>> 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<int, double>::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<int, double>::iterator it = indToParam.begin();
+ for (int m = 0; it != indToParam.end(); it++, m++) // iterate on intermediate profiles
+ {
+ const std::vector<gp_Pnt2d> &im_prof = IntermProf[m];
+ double param = it->second;
+ gp_Pnt2d LP;
+ gp_Pnt2d RP;
+ if (isHaxNormal)
+ {
+ const std::pair<IntRes2d_IntersectionPoint,
+ IntRes2d_IntersectionPoint> &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;
+ }
}
