1 // Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE
3 // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
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.
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.
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
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
23 #include "HYDROData_StreamLinearInterpolation.h"
24 #include <HYDROData_Profile.h>
25 #include <HYDROData_ProfileUZ.h>
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>
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>
45 #include <BRepBuilderAPI_MakeEdge2d.hxx>
46 #include <BRepBuilderAPI_MakeEdge.hxx>
47 #include <BRep_Tool.hxx>
48 #include <BRepBuilderAPI_MakeWire.hxx>
49 #include <TopoDS_Face.hxx>
51 #include <BRepBuilderAPI_MakeFace.hxx>
54 #include <Geom_Plane.hxx>
55 #include <TopoDS_Wire.hxx>
56 #include <BRepLib.hxx>
58 static std::vector<gp_Pnt2d> GetProfileUZPoints(const Handle(HYDROData_Profile)& theProfile)
60 Handle(HYDROData_ProfileUZ) aProfileUZ = theProfile->GetProfileUZ( false );
61 HYDROData_ProfileUZ::PointsList aSectPointsList = aProfileUZ->GetPoints();
62 std::vector<gp_Pnt2d> points;
63 points.reserve(aSectPointsList.Size());
64 for ( int i = 1; i <= aSectPointsList.Size(); i++ )
66 const HYDROData_ProfileUZ::Point& aSectPoint = aSectPointsList(i);
67 points.push_back( gp_Pnt2d(aSectPoint.X(), aSectPoint.Y()));
72 static void GetMaxDist(const std::vector<gp_Pnt2d>& PNTS,
80 std::vector<double> dist;
81 dist.reserve(PNTS.size()-1);
82 for (int i=0; i<PNTS.size()-1; i++)
84 //dist.push_back(PNTS[i].Distance(PNTS[i+1])); //via distance between p1/p2
85 dist.push_back(std::abs(PNTS[i].X()-PNTS[i+1].X())); //via U-diff of p1/p2
87 for (int i=0; i<dist.size(); i++)
96 for (int i=0; i<dist.size(); i++)
99 if (d > dmax2 && d < dmax)
104 static void InsertPoints(std::vector<gp_Pnt2d>& points, int nbpoints)
106 points.reserve(points.size() + nbpoints);
109 double dmax=0, dmax2=0;
111 GetMaxDist(points, dmax, dmax2, imax);
114 nbPins = floor(dmax/dmax2);
116 nbPins = nbpoints; //one segment, put all points here
117 //insert nbPins points in [imax, imax+1] segment
118 gp_Pnt2d p1 = points[imax];
119 gp_Pnt2d p2 = points[imax+1];
124 for (int i=0;i<nbPins;i++)
126 double t = ((double)i+1)/((double)nbPins+1);
127 double Xp = X0 + (X1-X0)*t;
128 double Yp = Y0 + (Y1-Y0)*t;
130 points.insert(points.begin() + imax + 1, P);
136 static void PolyToCurve2d(const Handle(HYDROData_PolylineXY)& poly, Handle(Geom2d_Curve)& c2d)
140 TopoDS_Shape sh = poly->GetShape();
141 if (sh.ShapeType() == TopAbs_COMPOUND)
143 TopoDS_Iterator it(sh);
152 if (sh.ShapeType() == TopAbs_EDGE)
153 ee = TopoDS::Edge(sh);
154 else if (sh.ShapeType() == TopAbs_WIRE)
155 ee = BRepAlgo::ConcatenateWireC0(TopoDS::Wire(sh)); //convert to bspline
160 BRepAdaptor_Curve Ad(ee);
162 c2d = HYDROData_Tool::BRepAdaptorTo2DCurve(Ad);
165 static void InterProfilesAndHAX(const HYDROData_SequenceOfObjects& profiles,
166 const Handle(Geom2d_Curve)& Hax2d,
167 std::map<double, Handle(HYDROData_Profile)>& profToInterParam,
168 std::vector<std::string>* warnings)
170 for (int i=1;i<=profiles.Size();i++)
172 Handle(HYDROData_Profile) aProfile = Handle(HYDROData_Profile)::DownCast(profiles(i));
174 aProfile->GetLeftPoint( LP, true, true );
175 aProfile->GetRightPoint( RP, true, true );
176 gp_Pnt2d P1(LP), P2(RP);
177 double d = P2.Distance(P1);
178 if (d < gp::Resolution())
181 warnings->push_back(aProfile->GetName().toStdString() + " is skipped: Left and Right points is the same");
184 Handle(Geom2d_Line) lin2d = new Geom2d_Line(P1,gp_Dir2d(P2.XY()-P1.XY()));
186 Geom2dAdaptor_Curve linAd(lin2d, 0, d);
187 Geom2dAdaptor_Curve haxAd(Hax2d);
188 Geom2dInt_GInter isec( linAd, haxAd, 1.0e-6, 1.0e-6);
192 warnings->push_back(aProfile->GetName().toStdString() + " is skipped: intersection between hydraulic axis & profile is failed");
195 if (isec.NbPoints() == 0)
198 warnings->push_back(aProfile->GetName().toStdString() + " is skipped: no intersection between hydraulic axis & profile");
201 if (isec.NbPoints() > 1)
204 warnings->push_back(aProfile->GetName().toStdString() + " is skipped; intersection between hydraulic axis & profile produces more than 1 point");
207 double param = isec.Point(1).ParamOnSecond();
208 profToInterParam[param] = aProfile;
212 static bool EquidParamsOnHAX(const Handle(Geom2d_Curve)& Hax2d, double step, double fp, double lp, std::set<double>& params)
214 Geom2dAdaptor_Curve ad(Hax2d);
215 GCPnts_UniformAbscissa GUA(ad, step, fp, lp);
219 for(int i = 1; i<= GUA.NbPoints(); i++)
220 params.insert(GUA.Parameter(i));
225 static void GetPointsOnBanks(const std::vector<double>& gua_params,
226 const Handle(Geom2d_Curve)& Hax2d,
227 const Handle(Geom2d_Curve)& LB2d,
228 const Handle(Geom2d_Curve)& RB2d,
229 std::map<double, std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>>& parToBankPoints,
230 std::vector<std::string>* warnings,
231 const std::map<double, std::pair<Handle(HYDROData_Profile), Handle(HYDROData_Profile)>>* intermParamToProfPair)
233 for(int i = 0; i < gua_params.size(); i++)
235 double par = gua_params[i];
238 Hax2d->D1(par, p, v); //if non-diff?
239 gp_Pnt2d outpLB, outpRB;
241 gp_Dir2d n2d(v.X(), v.Y()); //ignore Z (Z==0)
243 //nLine2d is an orthogonal line to hydr. axis
244 gp_Pnt2d p2d(gp_Pnt2d(p.X(), p.Y()));
245 Handle(Geom2d_Line) nLine2d = new Geom2d_Line(p2d, n2d);
246 //intersect nLine2d with banks
247 Geom2dAPI_InterCurveCurve intL(LB2d, nLine2d);
248 Geom2dAPI_InterCurveCurve intR(RB2d, nLine2d);
250 if (warnings && intL.NbPoints()==0 )
252 std::string coord = (QString::number(p2d.X()) + ", "+QString::number(p2d.Y())).toStdString();
253 if (intermParamToProfPair && intermParamToProfPair->count(par) != 0)
255 std::string fpn = intermParamToProfPair->at(par).first->GetName().toStdString();
256 std::string spn = intermParamToProfPair->at(par).second->GetName().toStdString();
257 warnings->push_back("no intersection between intermediate profile (point on h.axis:("+
258 coord +"), located between " + fpn + " & " + spn + ") and left bank found; skipped");
261 warnings->push_back("no intersection between intermediate profile (point on h.axis:("+ coord +")) and left bank found; skipped");
264 if (warnings && intR.NbPoints()==0)
266 std::string coord = (QString::number(p2d.X()) + ", "+QString::number(p2d.Y())).toStdString();
267 if (intermParamToProfPair && intermParamToProfPair->count(par) != 0)
269 std::string fpn = intermParamToProfPair->at(par).first->GetName().toStdString();
270 std::string spn = intermParamToProfPair->at(par).second->GetName().toStdString();
271 warnings->push_back("no intersection between intermediate profile (point on h.axis:("+
272 coord +"), located between " + fpn + " & " + spn + ") and right bank found; skipped");
275 warnings->push_back("no intersection between intermediate profile (point on h.axis:("+ coord +")) and right bank found; skipped");
278 if ( intL.NbPoints()==0 || intR.NbPoints()==0)
281 IntRes2d_IntersectionPoint aNearSolL, aNearSolR;
282 double min_sq_dist = Precision::Infinite();
283 for (int j=1;j<=intL.NbPoints();j++)
285 double sq_dist = p2d.SquareDistance(intL.Point(j));
286 if (min_sq_dist > sq_dist)
288 min_sq_dist = sq_dist;
289 //aNearSolL = intL.Point(j);
290 aNearSolL = intL.Intersector().Point(j);
293 min_sq_dist = Precision::Infinite();
294 for (int j=1;j<=intR.NbPoints();j++)
296 double sq_dist = p2d.SquareDistance(intR.Point(j));
297 if (min_sq_dist > sq_dist)
299 min_sq_dist = sq_dist;
300 //aNearSolR = intR.Point(j);
301 aNearSolR = intR.Intersector().Point(j);
305 std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint> int_pair(aNearSolL, aNearSolR);
306 parToBankPoints[par]=int_pair;
311 void BuildFace(const std::map<double, std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>>& parToBankPoints,
312 const Handle(Geom2d_Curve)& LB2d,
313 const Handle(Geom2d_Curve)& RB2d,
314 HYDROData_Stream::PrsDefinition& prsDef)
316 prsDef.myPrs3D.Nullify();
317 prsDef.myPrs2D.Nullify();
318 prsDef.myLeftBank.Nullify();
319 prsDef.myRightBank.Nullify();
320 prsDef.myInlet.Nullify();
321 prsDef.myOutlet.Nullify();
323 double lb_1, lb_2, rb_1, rb_2;
324 lb_1 = parToBankPoints.begin()->second.first.ParamOnFirst();
325 rb_1 = parToBankPoints.begin()->second.second.ParamOnFirst();
327 lb_2 = parToBankPoints.rbegin()->second.first.ParamOnFirst();
328 rb_2 = parToBankPoints.rbegin()->second.second.ParamOnFirst();
330 //stat = GetItersectionParam(F_prof, LB2d, lb_1) && GetItersectionParam(F_prof, RB2d, rb_1)
331 // && GetItersectionParam(S_prof, LB2d, lb_2) && GetItersectionParam(S_prof, RB2d, rb_2);
333 BRepBuilderAPI_MakeEdge2d LEM(LB2d, lb_1, lb_2);
336 BRepBuilderAPI_MakeEdge2d REM(RB2d, rb_1, rb_2);
339 TopoDS_Edge LBE = LEM.Edge();
340 TopoDS_Edge RBE = REM.Edge();
341 if (LBE.IsNull() || RBE.IsNull())
343 BRepBuilderAPI_MakeEdge PFEM(BRep_Tool::Pnt(LEM.Vertex1()), BRep_Tool::Pnt(REM.Vertex1()));
344 BRepBuilderAPI_MakeEdge PLEM(BRep_Tool::Pnt(LEM.Vertex2()), BRep_Tool::Pnt(REM.Vertex2()));
349 TopoDS_Edge FProfE = PFEM.Edge();
350 TopoDS_Edge SProfE = PLEM.Edge();
351 BRepBuilderAPI_MakeWire WM(FProfE, LBE, SProfE, RBE);
354 TopoDS_Wire W = WM.Wire();
355 BRepLib::BuildCurves3d(W);
356 BRepBuilderAPI_MakeFace FM(Geom_Plane(gp::XOY()).Pln(), WM.Wire());
359 prsDef.myPrs2D = FM.Face();
360 prsDef.myLeftBank = LBE;
361 prsDef.myRightBank = RBE;
362 prsDef.myInlet = FProfE;
363 prsDef.myOutlet = SProfE;
368 void HYDROData_StreamLinearInterpolation::Perform(const HYDROData_SequenceOfObjects& profiles,
371 const Handle(HYDROData_PolylineXY)& hax,
372 const Handle(HYDROData_PolylineXY)& LB,
373 const Handle(HYDROData_PolylineXY)& RB,
374 HYDROData_Bathymetry::AltitudePoints& outBathypoints,
375 bool buildPresentationShapes,
376 bool estimateWarnOnly,
377 HYDROData_Stream::PrsDefinition& prsDef,
378 std::vector<std::string>* warnings)
380 //intersect profiles with given hydr.axis
381 if (hax.IsNull() || LB.IsNull() || RB.IsNull())
384 Handle(Geom2d_Curve) Hax2d;
385 PolyToCurve2d(hax, Hax2d);
386 if (hax->IsClosed() || Hax2d->IsClosed()) //can't be closed
389 warnings->push_back(hax->GetName().toStdString() + " is closed; abort");
393 std::map<double, Handle(HYDROData_Profile)> InterParamToProf;
394 //profToInterParam is output param map: profile to intersection param on hydr.axis
395 InterProfilesAndHAX(profiles, Hax2d, InterParamToProf, warnings);
397 std::vector<std::vector<gp_Pnt2d>> profilesPoints;
398 profilesPoints.reserve(InterParamToProf.size());
401 //for (int i=1;i<=InterParamToProf.Extent();i++)
402 for( std::map<double, Handle(HYDROData_Profile)>::iterator it = InterParamToProf.begin(); it != InterParamToProf.end(); ++it )
404 //ordered by param on hydr axis
405 Handle(HYDROData_Profile) aProfile = it->second;
406 const std::vector<gp_Pnt2d>& profile_points = GetProfileUZPoints(aProfile);
407 profilesPoints.push_back(profile_points);
408 if (profile_points.size() > maxNbPoints)
409 maxNbPoints = profile_points.size();
411 std::set<double> paramsOnHAX;
413 if (InterParamToProf.size() <= 1)
416 warnings->push_back("Insufficient count of correct input profiles; abort");
420 EquidParamsOnHAX(Hax2d, stepOnHA, InterParamToProf.begin()->first,
421 InterParamToProf.rbegin()->first, paramsOnHAX);
423 //prepare ParamsPerSegm - number of interm.profiles for profile_points
424 std::vector<std::vector<double>> ParamsPerSegm;
425 std::set<double>::iterator it_params_hax = paramsOnHAX.begin();
426 std::map<double, Handle(HYDROData_Profile)>::iterator it_p = InterParamToProf.begin();
428 std::map<double, std::pair<Handle(HYDROData_Profile), Handle(HYDROData_Profile)>> intermParamToProfPair; //for warnings only
431 ParamsPerSegm.resize(profilesPoints.size()-1);
432 for( int k=0; it_p != InterParamToProf.end(); )
434 if (*it_params_hax < it_p->first)
436 double val = *it_params_hax;
437 ParamsPerSegm[k].push_back(val);
441 const Handle(HYDROData_Profile)& cp = it_p->second;
443 const Handle(HYDROData_Profile)& np = it_p->second;
444 std::pair<Handle(HYDROData_Profile), Handle(HYDROData_Profile)> profPair(cp, np);
445 intermParamToProfPair[val] = profPair;
453 ParamsPerSegm.back().push_back(*paramsOnHAX.rbegin());
455 std::map<double, std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>> parToBankPoints;
456 std::vector<double> paramHAXVec;
457 paramHAXVec.reserve(paramsOnHAX.size());
458 for( std::set<double>::iterator it = paramsOnHAX.begin(); it != paramsOnHAX.end(); ++it )
459 paramHAXVec.push_back(*it);
461 Handle(Geom2d_Curve) LB2d, RB2d;
462 PolyToCurve2d(LB, LB2d);
463 PolyToCurve2d(RB, RB2d);
464 GetPointsOnBanks( paramHAXVec, Hax2d, LB2d, RB2d, parToBankPoints, warnings, &intermParamToProfPair);
466 if (buildPresentationShapes)
468 //parToBankPoints.begin()->second.;
469 //double lp1 = parToBankPoints.rbegin()->first;
471 //Handle(HYDROData_Profile) F_prof = InterParamToProf.begin()->second;
472 //Handle(HYDROData_Profile) S_prof = InterParamToProf.rbegin()->second;
473 BuildFace(parToBankPoints, LB2d, RB2d, prsDef);
477 if (estimateWarnOnly)
480 maxNbPoints = Max(pointsToInsert, maxNbPoints);
482 //insert points to profiles
483 for (int i=0;i<profilesPoints.size();i++)
485 int nbPointsToInsert = maxNbPoints - profilesPoints[i].size();
486 //insert nbPointsToInsert points to current profile
487 InsertPoints(profilesPoints[i], nbPointsToInsert);
490 for (int i=0;i<profilesPoints.size()-1;i++)
492 const std::vector<gp_Pnt2d>& prof1 = profilesPoints[i];
493 const std::vector<gp_Pnt2d>& prof2 = profilesPoints[i+1];
494 std::vector<std::vector<gp_Pnt2d>> IntermProf;
496 std::map<int, double> indToParam; //index of im.profile to param on hax. ()
497 for (int j=0;j<ParamsPerSegm[i].size();j++)
499 double param = ParamsPerSegm[i][j];
500 if (parToBankPoints.count(param) > 0) //not intersected with banks; skip this
503 int NbIntermProf = indToParam.size();
504 IntermProf.resize(NbIntermProf);
506 for (int l=0;l<NbIntermProf;l++)
507 IntermProf[l].resize(prof1.size());
508 //calc interm. between prof1 & prof2
509 assert (prof1.size() == prof2.size());
510 for (int j=0;j<prof1.size();j++)
512 gp_Pnt2d P1 = prof1[j];
513 gp_Pnt2d P2 = prof2[j];
518 std::map<int, double>::iterator it = indToParam.begin();
519 for ( int m=0; it != indToParam.end(); it++, m++ )
522 double t = ((double)ind+1)/((double)NbIntermProf+1);
523 double Xp = X0 + (X1-X0)*t;
524 double Yp = Y0 + (Y1-Y0)*t;
526 IntermProf[m][j] = P;
529 std::map<int, double>::iterator it = indToParam.begin();
530 for ( int m=0; it != indToParam.end(); it++, m++ )
532 const std::vector<gp_Pnt2d>& im_prof = IntermProf[m];
533 double param = it->second;
534 const std::pair<IntRes2d_IntersectionPoint, IntRes2d_IntersectionPoint>& BB = parToBankPoints[param]; //param is included in map; no checks
535 gp_Pnt2d LP = BB.first.Value();
536 gp_Pnt2d RP = BB.second.Value();
537 HYDROData_ProfileUZ::PointsList pl;
538 for (int k=0;k<im_prof.size();k++)
539 pl.Append(im_prof[k].XY());
540 HYDROData_Profile::ProfilePoints profile_points3d = HYDROData_Profile::CalculateProfilePoints(pl, LP.XY(), RP.XY());
541 for (int k=1;k<=profile_points3d.Length();k++)
543 HYDROData_Bathymetry::AltitudePoint AP(profile_points3d(k).X(), profile_points3d(k).Y(), profile_points3d(k).Z());
544 outBathypoints.push_back(AP);
548 //outBathy->SetAltitudePoints(bathypoints);