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<gp_Pnt2d, gp_Pnt2d>>& parToBankPoints,
230 std::vector<std::string>* warnings)
232 for(int i = 0; i < gua_params.size(); i++)
234 double par = gua_params[i];
237 Hax2d->D1(par, p, v); //if non-diff?
238 gp_Pnt2d outpLB, outpRB;
240 gp_Dir2d n2d(v.X(), v.Y()); //ignore Z (Z==0)
242 //nLine2d is an orthogonal line to hydr. axis
243 gp_Pnt2d p2d(gp_Pnt2d(p.X(), p.Y()));
244 Handle(Geom2d_Line) nLine2d = new Geom2d_Line(p2d, n2d);
245 //intersect nLine2d with banks
246 Geom2dAPI_InterCurveCurve intL(LB2d, nLine2d);
247 Geom2dAPI_InterCurveCurve intR(RB2d, nLine2d);
249 if (warnings && intL.NbPoints()==0 )
251 std::string coord = (QString::number(p2d.X()) + ", "+QString::number(p2d.Y())).toStdString();
252 warnings->push_back("no intersection between intermediate profile (point on h.axis:("+ coord +")) and left bank found; skipped");
255 if (warnings && intR.NbPoints()==0)
257 std::string coord = (QString::number(p2d.X()) + ", "+QString::number(p2d.Y())).toStdString();
258 warnings->push_back("no intersection between intermediate profile (point on h.axis:("+ coord +")) and right bank found; skipped");
261 if ( intL.NbPoints()==0 || intR.NbPoints()==0)
264 gp_Pnt2d aNearSolL, aNearSolR;
265 double min_sq_dist = Precision::Infinite();
266 for (int j=1;j<=intL.NbPoints();j++)
268 double sq_dist = p2d.SquareDistance(intL.Point(j));
269 if (min_sq_dist > sq_dist)
271 min_sq_dist = sq_dist;
272 aNearSolL = intL.Point(j);
275 min_sq_dist = Precision::Infinite();
276 for (int j=1;j<=intR.NbPoints();j++)
278 double sq_dist = p2d.SquareDistance(intR.Point(j));
279 if (min_sq_dist > sq_dist)
281 min_sq_dist = sq_dist;
282 aNearSolR = intR.Point(j);
286 std::pair<gp_Pnt2d, gp_Pnt2d> int_pair(aNearSolL, aNearSolR);
287 parToBankPoints[par]=int_pair;
292 bool GetItersectionParam(const Handle(HYDROData_Profile)& aProfile, const Handle(Geom2d_Curve)& B2d, double& outIntersectionParam)
295 aProfile->GetLeftPoint( LP, true, true );
296 aProfile->GetRightPoint( RP, true, true );
297 gp_Pnt2d P1(LP), P2(RP);
298 double d = P2.Distance(P1);
299 if (d < gp::Resolution())
301 Handle(Geom2d_Line) lin2d = new Geom2d_Line(P1,gp_Dir2d(P2.XY()-P1.XY()));
303 Geom2dAdaptor_Curve linAd(lin2d, 0, d);
304 Geom2dAdaptor_Curve BAd(B2d);
305 Geom2dInt_GInter isec( linAd, BAd, 1.0e-6, 1.0e-6);
308 if (isec.NbPoints() == 0)
310 if (isec.NbPoints() > 1)
312 outIntersectionParam = isec.Point(1).ParamOnSecond();
316 void BuildFace(const Handle(HYDROData_Profile)& F_prof,
317 const Handle(HYDROData_Profile)& S_prof,
318 const Handle(Geom2d_Curve)& LB2d,
319 const Handle(Geom2d_Curve)& RB2d,
320 HYDROData_Stream::PrsDefinition& prsDef)
322 prsDef.myPrs3D.Nullify();
323 prsDef.myPrs2D.Nullify();
324 prsDef.myLeftBank.Nullify();
325 prsDef.myRightBank.Nullify();
326 prsDef.myInlet.Nullify();
327 prsDef.myOutlet.Nullify();
330 double lb_1, lb_2, rb_1, rb_2;
331 stat = GetItersectionParam(F_prof, LB2d, lb_1) && GetItersectionParam(F_prof, RB2d, rb_1)
332 && GetItersectionParam(S_prof, LB2d, lb_2) && GetItersectionParam(S_prof, RB2d, rb_2);
335 BRepBuilderAPI_MakeEdge2d LEM(LB2d, lb_1, lb_2);
338 BRepBuilderAPI_MakeEdge2d REM(RB2d, rb_1, rb_2);
341 TopoDS_Edge LBE = LEM.Edge();
342 TopoDS_Edge RBE = REM.Edge();
343 if (LBE.IsNull() || RBE.IsNull())
345 BRepBuilderAPI_MakeEdge PFEM(BRep_Tool::Pnt(LEM.Vertex1()), BRep_Tool::Pnt(REM.Vertex1()));
346 BRepBuilderAPI_MakeEdge PLEM(BRep_Tool::Pnt(LEM.Vertex2()), BRep_Tool::Pnt(REM.Vertex2()));
351 TopoDS_Edge FProfE = PFEM.Edge();
352 TopoDS_Edge SProfE = PLEM.Edge();
353 BRepBuilderAPI_MakeWire WM(FProfE, LBE, SProfE, RBE);
356 TopoDS_Wire W = WM.Wire();
357 BRepLib::BuildCurves3d(W);
358 BRepBuilderAPI_MakeFace FM(Geom_Plane(gp::XOY()).Pln(), WM.Wire());
361 prsDef.myPrs2D = FM.Face();
362 prsDef.myLeftBank = LBE;
363 prsDef.myRightBank = RBE;
364 prsDef.myInlet = FProfE;
365 prsDef.myOutlet = SProfE;
371 void HYDROData_StreamLinearInterpolation::Perform(const HYDROData_SequenceOfObjects& profiles,
374 const Handle(HYDROData_PolylineXY)& hax,
375 const Handle(HYDROData_PolylineXY)& LB,
376 const Handle(HYDROData_PolylineXY)& RB,
377 HYDROData_Bathymetry::AltitudePoints& outBathypoints,
378 bool buildPresentationShapes,
379 bool estimateWarnOnly,
380 HYDROData_Stream::PrsDefinition& prsDef,
381 std::vector<std::string>* warnings)
383 //intersect profiles with given hydr.axis
384 if (hax.IsNull() || LB.IsNull() || RB.IsNull())
387 Handle(Geom2d_Curve) Hax2d;
388 PolyToCurve2d(hax, Hax2d);
389 if (Hax2d->IsClosed()) //can't be closed
392 warnings->push_back(hax->GetName().toStdString() + " is closed; abort");
396 std::map<double, Handle(HYDROData_Profile)> InterParamToProf;
397 //profToInterParam is output param map: profile to intersection param on hydr.axis
398 InterProfilesAndHAX(profiles, Hax2d, InterParamToProf, warnings);
400 std::vector<std::vector<gp_Pnt2d>> profilesPoints;
401 profilesPoints.reserve(InterParamToProf.size());
404 //for (int i=1;i<=InterParamToProf.Extent();i++)
405 for( std::map<double, Handle(HYDROData_Profile)>::iterator it = InterParamToProf.begin(); it != InterParamToProf.end(); ++it )
407 //ordered by param on hydr axis
408 Handle(HYDROData_Profile) aProfile = it->second;
409 const std::vector<gp_Pnt2d>& profile_points = GetProfileUZPoints(aProfile);
410 profilesPoints.push_back(profile_points);
411 if (profile_points.size() > maxNbPoints)
412 maxNbPoints = profile_points.size();
414 std::set<double> paramsOnHAX;
416 if (InterParamToProf.size() <= 1)
419 warnings->push_back("Insufficient count of correct input profiles; abort");
423 EquidParamsOnHAX(Hax2d, stepOnHA, InterParamToProf.begin()->first,
424 InterParamToProf.rbegin()->first, paramsOnHAX);
426 //prepare ParamsPerSegm - number of interm.profiles for profile_points
427 std::vector<std::vector<double>> ParamsPerSegm;
428 std::set<double>::iterator it_params_hax = paramsOnHAX.begin();
429 std::map<double, Handle(HYDROData_Profile)>::iterator it_p = InterParamToProf.begin();
431 ParamsPerSegm.resize(profilesPoints.size()-1);
432 for( int k=0; it_p != InterParamToProf.end(); )
434 if (*it_params_hax < it_p->first)
436 ParamsPerSegm[k].push_back(*it_params_hax);
445 ParamsPerSegm.back().push_back(*paramsOnHAX.rbegin());
447 std::map<double, std::pair<gp_Pnt2d, gp_Pnt2d>> parToBankPoints;
448 std::vector<double> paramHAXVec;
449 paramHAXVec.reserve(paramsOnHAX.size());
450 for( std::set<double>::iterator it = paramsOnHAX.begin(); it != paramsOnHAX.end(); ++it )
451 paramHAXVec.push_back(*it);
453 Handle(Geom2d_Curve) LB2d, RB2d;
454 PolyToCurve2d(LB, LB2d);
455 PolyToCurve2d(RB, RB2d);
456 GetPointsOnBanks( paramHAXVec, Hax2d, LB2d, RB2d, parToBankPoints, warnings);
458 if (buildPresentationShapes)
460 Handle(HYDROData_Profile) F_prof = InterParamToProf.begin()->second;
461 Handle(HYDROData_Profile) S_prof = InterParamToProf.rbegin()->second;
462 BuildFace(F_prof, S_prof, LB2d, RB2d, prsDef);
466 if (estimateWarnOnly)
469 maxNbPoints = Max(pointsToInsert, maxNbPoints);
471 //insert points to profiles
472 for (int i=0;i<profilesPoints.size();i++)
474 int nbPointsToInsert = maxNbPoints - profilesPoints[i].size();
475 //insert nbPointsToInsert points to current profile
476 InsertPoints(profilesPoints[i], nbPointsToInsert);
479 for (int i=0;i<profilesPoints.size()-1;i++)
481 const std::vector<gp_Pnt2d>& prof1 = profilesPoints[i];
482 const std::vector<gp_Pnt2d>& prof2 = profilesPoints[i+1];
483 std::vector<std::vector<gp_Pnt2d>> IntermProf;
485 std::map<int, double> indToParam; //index of im.profile to param on hax. ()
486 for (int j=0;j<ParamsPerSegm[i].size();j++)
488 double param = ParamsPerSegm[i][j];
489 if (parToBankPoints.count(param) > 0) //not intersected with banks; skip this
492 int NbIntermProf = indToParam.size();
493 IntermProf.resize(NbIntermProf);
495 for (int l=0;l<NbIntermProf;l++)
496 IntermProf[l].resize(prof1.size());
497 //calc interm. between prof1 & prof2
498 assert (prof1.size() == prof2.size());
499 for (int j=0;j<prof1.size();j++)
501 gp_Pnt2d P1 = prof1[j];
502 gp_Pnt2d P2 = prof2[j];
507 std::map<int, double>::iterator it = indToParam.begin();
508 for ( int m=0; it != indToParam.end(); it++, m++ )
511 double t = ((double)ind+1)/((double)NbIntermProf+1);
512 double Xp = X0 + (X1-X0)*t;
513 double Yp = Y0 + (Y1-Y0)*t;
515 IntermProf[m][j] = P;
518 std::map<int, double>::iterator it = indToParam.begin();
519 for ( int m=0; it != indToParam.end(); it++, m++ )
521 const std::vector<gp_Pnt2d>& im_prof = IntermProf[m];
522 double param = it->second;
523 const std::pair<gp_Pnt2d, gp_Pnt2d>& BB = parToBankPoints[param]; //param is included in map; no checks
524 gp_Pnt2d LP = BB.first;
525 gp_Pnt2d RP = BB.second;
526 HYDROData_ProfileUZ::PointsList pl;
527 for (int k=0;k<im_prof.size();k++)
528 pl.Append(im_prof[k].XY());
529 HYDROData_Profile::ProfilePoints profile_points3d = HYDROData_Profile::CalculateProfilePoints(pl, LP.XY(), RP.XY());
530 for (int k=1;k<=profile_points3d.Length();k++)
532 HYDROData_Bathymetry::AltitudePoint AP(profile_points3d(k).X(), profile_points3d(k).Y(), profile_points3d(k).Z());
533 outBathypoints.push_back(AP);
537 //outBathy->SetAltitudePoints(bathypoints);