#include <GCPnts_AbscissaPoint.hxx>
#include <BRepBuilderAPI_MakeEdge.hxx>
#include <limits>
+#include <Bnd_Box2d.hxx>
#include <BRepLib_MakeEdge.hxx>
#include <BRepLib_MakeWire.hxx>
#include <TopTools_ListIteratorOfListOfShape.hxx>
#include <TopTools_SequenceOfShape.hxx>
#include <assert.h>
+#include <NCollection_DataMap.hxx>
+#include <QSet>
+#include <QString>
#include <float.h>
+double PREC = 0.001;
IMPLEMENT_STANDARD_RTTIEXT( HYDROData_DTM, HYDROData_Bathymetry )
HYDROData_DTM::CurveUZ::CurveUZ( double theXCurv, const gp_Vec2d& theProfileDir, double theDeltaZ, double theMaxZ )
double step = GetSpatialStep();
std::set<int> InvInd;
bool WireIntersections; //__TODO
- CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, true, true, InvInd, -1, WireIntersections );
+ myWarnings.Clear();
+ bool ToEstimateWarnings; //NOT NEEDED here
+ CreateProfilesFromDTM( objs, ddz, step, points, Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet,
+ true, true, InvInd, -1, WireIntersections, myWarnings, ToEstimateWarnings );
SetAltitudePoints( points );
SetShape( DataTag_LeftBankShape, OutLeftB);
bool Create2dPres,
std::set<int>& InvInd,
int thePntsLimit,
- bool& WireIntersections)
+ bool& WireIntersections,
+ NCollection_DataMap<Handle(HYDROData_Profile), QSet<QString>>& warnings,
+ bool& ToEstimateWarnings)
{
int aLower = InpProfiles.Lower(), anUpper = InpProfiles.Upper();
size_t n = anUpper - aLower + 1;
AltitudePoints right;
std::vector<AltitudePoints> main_profiles;
+ bool ToEstimateWarningsOnly = false;
if( thePntsLimit > 0 )
{
int aNbPoints = EstimateNbPoints( profiles, ddz, step );
if( aNbPoints < 0 || aNbPoints > thePntsLimit )
- return;
+ {
+ ToEstimateWarningsOnly = true;
+ //return;
+ }
}
+ if( profiles.size() < 2 )
+ return;
+
if( ddz>EPS && step>EPS )
CreateProfiles(profiles, ddz, step, left, right, points, main_profiles,
- Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres, Create2dPres, InvInd, WireIntersections );
+ Out3dPres, Out2dPres, OutLeftB, OutRightB, OutInlet, OutOutlet, Create3dPres,
+ Create2dPres, InvInd, WireIntersections, warnings, ToEstimateWarningsOnly );
+
+ ToEstimateWarnings = ToEstimateWarningsOnly;
}
bool HYDROData_DTM::GetPlanarFaceFromBanks( const TopoDS_Edge& LB, const TopoDS_Edge& RB, TopoDS_Face& outF,
bool Create3dPres,
bool Create2dPres,
std::set<int>& InvInd,
- bool& ProjStat)
+ bool& ProjStat,
+ NCollection_DataMap<Handle(HYDROData_Profile), QSet<QString>>& warnings,
+ bool ToEstimateWarningsOnly)
{
if (theProfiles.empty())
return;
- theOutPoints = Interpolate( theProfiles, theDDZ, theSpatialStep, theOutLeft, theOutRight, theOutMainProfiles, InvInd );
+ theOutPoints = Interpolate( theProfiles, theDDZ, theSpatialStep, theOutLeft, theOutRight,
+ theOutMainProfiles, InvInd, warnings, ToEstimateWarningsOnly );
//note that if Create3dPres is false => Create2dPres flag is meaningless!
if( theOutPoints.empty() )
return;
}
-bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
+//bool CalcMidWidth( const std::set<double>& intersections, double& theMid, double& theWid )
+//{
+// double umin = std::numeric_limits<double>::max(),
+// umax = -umin;
+//
+// size_t n = intersections.size();
+// if( n <= 0 )
+// return false;
+//
+// std::set<double>::const_iterator it = intersections.begin(), last = intersections.end();
+// for( ; it!=last; it++ )
+// {
+// double u = *it;
+// if( u<umin )
+// umin = u;
+// if( u>umax )
+// umax = u;
+// }
+// theMid = ( umin+umax )/2;
+// theWid = umax-umin;
+// return true;
+//}
+
+
+
+bool CalcMidWidth( const Bnd_Box2d& inters_bnd, double& theMid, double& theWid )
{
double umin = DBL_MAX,
umax = -umin;
- size_t n = intersections.size();
- if( n <= 0 )
+ if (inters_bnd.IsVoid())
return false;
+ //size_t n = intersections.size();
+ //if( n <= 0 )
+ // return false;
- std::set<double>::const_iterator it = intersections.begin(), last = intersections.end();
- for( ; it!=last; it++ )
- {
- double u = *it;
- if( u<umin )
- umin = u;
- if( u>umax )
- umax = u;
- }
- theMid = ( umin+umax )/2;
- theWid = umax-umin;
+ double xmin, ymin, xmax, ymax;
+ inters_bnd.Get(xmin, ymin, xmax, ymax);
+ if (Abs(ymax-ymin)>PREC)
+ return false;
+
+ theMid = ( xmax+xmin )/2;
+ theWid = xmax-xmin;
return true;
}
+
+
void HYDROData_DTM::ProfileDiscretization( const Handle(HYDROData_Profile)& theProfile,
double theXCurv, double theMinZ, double theMaxZ, double theTopZ, double theDDZ,
CurveUZ& theMidPointCurve,
CurveUZ& theWidthCurve,
int& intersection_nb,
- double theTolerance)
+ double theTolerance,
+ QSet<QString>& warnings)
{
- double aDblMax = DBL_MAX,
+ warnings.clear();
+ double aDblMax = std::numeric_limits<double>::max(),
aUMin = aDblMax,
aUMax = -aUMin,
aVMax = 1000000;
size_t n = curves.size();
if( n==0 )
+ {
+ warnings.insert("no curves for discretization; skipped");
return;
+ }
// we add the "virtual" vertical lines to simulate the intersection with profile
gp_Pnt2d aFirst, aLast;
theWidthCurve = CurveUZ( theXCurv, aProfileDir, theMinZ, theTopZ );
theWidthCurve.reserve( psize );
- n = curves.size();
+ {
+ bool PlatoCase = false;
+ if (Abs(aFirst.Y() - aLast.Y() > PREC))
+ {
+ warnings.insert("One of the extreme points is higher than another");
+ }
+ double ZminExtr = Min(aFirst.Y(), aLast.Y());
+ double ZmaxExtr = Max(aFirst.Y(), aLast.Y());
+ Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, ZminExtr ), gp_Dir2d( 1, 0 ) );
+ std::vector<gp_Pnt2d> intersections;
+ for( int i = 0; i < n; i++ )
+ {
+ Handle(Geom2d_Curve) aCurve = curves[i];
+ Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
+ for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
+ intersections.push_back(anIntersect.Point(k));
+ if (anIntersect.NbSegments() > 0 )
+ PlatoCase = true;
+ }
+ std::vector<gp_Pnt2d> interm_intersections;
+ for (int i=0;i<intersections.size();i++)
+ {
+ gp_Pnt2d int_p2d = intersections[i];
+ //check for intermid. points: shoudl be higher than ZminExtr and not in intersection with first,last points
+ if (aFirst.Distance(int_p2d) > PREC && aLast.Distance(int_p2d) > PREC && int_p2d.Y() > ZminExtr)
+ interm_intersections.push_back(int_p2d);
+ }
+ if (!interm_intersections.empty())
+ warnings.insert("One of the internal points is higher than extreme points");
+ if (ZminExtr != ZmaxExtr && !PlatoCase)
+ {
+ //additional check of plato for zmax
+ Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, ZmaxExtr ), gp_Dir2d( 1, 0 ) );
+ std::vector<gp_Pnt2d> intersections;
+ for( int i = 0; i < n; i++ )
+ {
+ Handle(Geom2d_Curve) aCurve = curves[i];
+ Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
+ if (anIntersect.NbSegments() > 0 )
+ {
+ PlatoCase = true;
+ break;
+ }
+ }
+ }
+ if (PlatoCase)
+ warnings.insert("Plato case on extremes");
+ }
+
+ n = curves.size();
+
// for each discrete value of z we search intersection with profile
for( double z1 = theMinZ; z1 <= theMaxZ; z1 += theDDZ )
{
Handle(Geom2d_Line) aLine = new Geom2d_Line( gp_Pnt2d( 0, z1 ), gp_Dir2d( 1, 0 ) );
std::set<double> intersections;
+ Bnd_Box2d intersect_bndbox;
for( size_t i = 0; i < n; i++ )
{
Handle(Geom2d_Curve) aCurve = curves[i];
Geom2dAPI_InterCurveCurve anIntersect( aCurve, aLine, theTolerance );
for( int k=1, m=anIntersect.NbPoints(); k<=m; k++ )
+ {
intersections.insert( anIntersect.Point( k ).X() );
+ intersect_bndbox.Add( anIntersect.Point(k));
+ }
+ //
+ for( int k=1, m=anIntersect.NbSegments(); k<=m; k++ )
+ {
+ Handle(Geom2d_Curve) Curve1,Curve2;
+ anIntersect.Segment(k, Curve1, Curve2 );
+ double f = Curve2->FirstParameter();
+ double l = Curve2->LastParameter();
+ gp_Pnt2d Pf, Pl;
+ Curve2->D0(f, Pf);
+ Curve2->D0(l, Pl);
+ intersect_bndbox.Add( Pf );
+ intersect_bndbox.Add( Pl );
+ intersections.insert( Pf.X() );
+ intersections.insert( Pl.X() );
+ }
}
intersection_nb = intersections.size();
- if( intersection_nb >= 1 )
+ if (intersection_nb == 0)
+ {
+ warnings.insert("No intersections between profile & altitude Z-lines found; skipped");
+ return;
+ }
+ else if (intersection_nb > 2)
+ {
+ warnings.insert("More than 2 intersections between profile & altitude Z-lines found");
+ }
+ double xmin, ymin, xmax, ymax;
+ intersect_bndbox.Get(xmin, ymin, xmax, ymax);
+ if (Abs(xmax-xmin)>PREC)
{
double u_mid, u_wid;
- if( !CalcMidWidth( intersections, u_mid, u_wid ) )
+ if( !CalcMidWidth( intersect_bndbox, u_mid, u_wid ) )
continue;
double z = z1 - theMinZ;
p_wid.Z = z;
theWidthCurve.push_back( p_wid );
}
+
+
+ //intersection_nb = intersections.size();
+ //if( intersection_nb >= 1 )
+ //{
+ // double u_mid, u_wid;
+ // if( !CalcMidWidth( intersections, u_mid, u_wid ) )
+ // continue;
+ //
+ // double z = z1 - theMinZ;
+ // PointUZ p_mid;
+ // p_mid.U = u_mid;
+ // p_mid.Z = z;
+ // theMidPointCurve.push_back( p_mid );
+ //
+ // PointUZ p_wid;
+ // p_wid.U = u_wid;
+ // p_wid.Z = z;
+ // theWidthCurve.push_back( p_wid );
+ //}
}
}
const Handle(HYDROData_Profile)& theProfileB,
double theXCurvB,
double theDDZ, int theNbSteps, bool isAddSecond,
- int& inter_nb_1, int& inter_nb_2)
+ int& inter_nb_1, int& inter_nb_2,
+ NCollection_DataMap<Handle(HYDROData_Profile), QSet<QString>>& warnings,
+ bool ToEstimateWarningsOnly)
{
double zminA, zmaxA, zminB, zmaxB;
gp_Pnt lowestA, lowestB;
CurveUZ midA(0, gp_Vec2d(), 0, 0), midB(0, gp_Vec2d(), 0, 0);
CurveUZ widA(0, gp_Vec2d(), 0, 0), widB(0, gp_Vec2d(), 0, 0);
- ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, zmaxA-zminA, theDDZ, midA, widA, inter_nb_1 );
- ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, zmaxB-zminB, theDDZ, midB, widB, inter_nb_2 );
+ QSet<QString> warnings_per_profileA, warnings_per_profileB;
+
+ ProfileDiscretization( theProfileA, theXCurvA, zminA, zminA+hmax, zmaxA-zminA,
+ theDDZ, midA, widA, inter_nb_1, 1E-6, warnings_per_profileA );
+ ProfileDiscretization( theProfileB, theXCurvB, zminB, zminB+hmax, zmaxB-zminB,
+ theDDZ, midB, widB, inter_nb_2, 1E-6, warnings_per_profileB );
+
+ //process warnings
+ if (warnings.IsBound(theProfileA))
+ {
+ QSet<QString>& warnings_per_profileA_old = warnings.ChangeFind(theProfileA);
+ warnings_per_profileA_old+=warnings_per_profileA;
+ }
+ else
+ warnings.Bind(theProfileA, warnings_per_profileA);
+
+ if (warnings.IsBound(theProfileB))
+ {
+ QSet<QString>& warnings_per_profileB_old = warnings.ChangeFind(theProfileB);
+ warnings_per_profileB_old+=warnings_per_profileB;
+ }
+ else
+ warnings.Bind(theProfileB, warnings_per_profileB);
+ //
+
+ if (ToEstimateWarningsOnly)
+ return std::vector<AltitudePoints>();
std::vector<CurveUZ> mid, wid;
Interpolate( midA, midB, theNbSteps, mid, isAddSecond );
AltitudePoints& theLeft,
AltitudePoints& theRight,
std::vector<AltitudePoints>& theMainProfiles,
- std::set<int>& invalInd)
+ std::set<int>& invalInd,
+ NCollection_DataMap<Handle(HYDROData_Profile), QSet<QString>>& warnings,
+ bool ToEstimateWarningsOnly)
{
AltitudePoints points;
size_t n = theProfiles.size();
return points;
std::vector<double> distances;
- Handle(Geom2d_BSplineCurve) aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
- if( aHydraulicAxis.IsNull() )
- return points;
+ Handle(Geom2d_BSplineCurve) aHydraulicAxis;
+
+ if (!ToEstimateWarningsOnly)
+ {
+ aHydraulicAxis = CreateHydraulicAxis( theProfiles, distances );
+ if( aHydraulicAxis.IsNull() )
+ return points;
+ }
theMainProfiles.reserve( n );
for( size_t i=0, n1=n-1; i<n1; i++ )
{
- double aDistance = distances[i+1]-distances[i];
- int aNbSteps = int(aDistance/theSpatialStep);
- bool isAddSecond = i==n1-1;
+ double aDistance = 0;
+ int aNbSteps = -1;
+ bool isAddSecond = false;
+
+ if (!ToEstimateWarningsOnly)
+ {
+ aDistance = distances[i+1]-distances[i];
+ aNbSteps = int(aDistance/theSpatialStep);
+ isAddSecond = i==n1-1;
+ }
// 1. Calculate interpolated profiles
int inter_nb_1, inter_nb_2;
- std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i], distances[i],
- theProfiles[i+1], distances[i+1], theDDZ, aNbSteps, isAddSecond, inter_nb_1, inter_nb_2 );
+ std::vector<AltitudePoints> local_points = Interpolate( aHydraulicAxis, theProfiles[i],
+ ToEstimateWarningsOnly ? 0 : distances[i],
+ theProfiles[i+1],
+ ToEstimateWarningsOnly ? 0 : distances[i+1],
+ theDDZ, aNbSteps, isAddSecond, inter_nb_1, inter_nb_2, warnings,
+ ToEstimateWarningsOnly);
int lps = local_points.size();
+ if (lps == 0)
+ continue;
+
if (inter_nb_1 > 2)
invalInd.insert(i);
return aNbSteps * aNbZSteps;
}
+
+void HYDROData_DTM::GetWarnings(NCollection_DataMap<Handle(HYDROData_Profile), QSet<QString>>& warnings)
+{
+ warnings = myWarnings;
+}
