diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp index 8d1b5d3530d860e5d713429ffa9565f9aa347927..3bcab2ef2728d1eb59ad7573649cd722c31eb96e 100644 --- a/GeoLib/AnalyticalGeometry.cpp +++ b/GeoLib/AnalyticalGeometry.cpp @@ -14,31 +14,27 @@ #include "AnalyticalGeometry.h" +#include <Eigen/Dense> #include <algorithm> #include <cmath> #include <limits> - -#include <Eigen/Dense> - #include "BaseLib/StringTools.h" - -#include "Polyline.h" -#include "PointVec.h" - #include "MathLib/GeometricBasics.h" +#include "PointVec.h" +#include "Polyline.h" -extern double orient2d(double *, double *, double *); +extern double orient2d(double*, double*, double*); extern double orient2dfast(double*, double*, double*); namespace ExactPredicates { -double getOrientation2d(MathLib::Point3d const& a, - MathLib::Point3d const& b, MathLib::Point3d const& c) +double getOrientation2d(MathLib::Point3d const& a, MathLib::Point3d const& b, + MathLib::Point3d const& c) { return orient2d(const_cast<double*>(a.getCoords()), - const_cast<double*>(b.getCoords()), - const_cast<double*>(c.getCoords())); + const_cast<double*>(b.getCoords()), + const_cast<double*>(c.getCoords())); } double getOrientation2dFast(MathLib::Point3d const& a, @@ -119,7 +115,8 @@ bool parallel(Eigen::Vector3d v, Eigen::Vector3d w) parallel = false; } - if (! parallel) { + if (!parallel) + { parallel = true; // change sense of direction of v_normalised v *= -1.0; @@ -185,14 +182,15 @@ bool lineSegmentIntersect(GeoLib::LineSegment const& s0, } auto isLineSegmentIntersectingAB = [&v](Eigen::Vector3d const& ap, - std::size_t i) - { + std::size_t i) { // check if p is located at v=(a,b): (ap = t*v, t in [0,1]) return 0.0 <= ap[i] / v[i] && ap[i] / v[i] <= 1.0; }; - if (parallel(v,w)) { // original line segments (a,b) and (c,d) are parallel - if (parallel(pq,v)) { // line segment (a,b) and (a,c) are also parallel + if (parallel(v, w)) + { // original line segments (a,b) and (c,d) are parallel + if (parallel(pq, v)) + { // line segment (a,b) and (a,c) are also parallel // Here it is already checked that the line segments (a,b) and (c,d) // are parallel. At this point it is also known that the line // segment (a,c) is also parallel to (a,b). In that case it is @@ -207,7 +205,8 @@ bool lineSegmentIntersect(GeoLib::LineSegment const& s0, // coordinate axis. std::size_t i_max(std::abs(v[0]) <= std::abs(v[1]) ? 1 : 0); i_max = std::abs(v[i_max]) <= std::abs(v[2]) ? 2 : i_max; - if (isLineSegmentIntersectingAB(qp, i_max)) { + if (isLineSegmentIntersectingAB(qp, i_max)) + { s = pc; return true; } @@ -227,10 +226,10 @@ bool lineSegmentIntersect(GeoLib::LineSegment const& s0, const double sqr_len_w(w.squaredNorm()); Eigen::Matrix2d mat; - mat(0,0) = sqr_len_v; - mat(0,1) = -v.dot(w); - mat(1,1) = sqr_len_w; - mat(1,0) = mat(0,1); + mat(0, 0) = sqr_len_v; + mat(0, 1) = -v.dot(w); + mat(1, 1) = sqr_len_w; + mat(1, 0) = mat(0, 1); Eigen::Vector2d rhs{v.dot(qp), w.dot(pq)}; @@ -238,21 +237,25 @@ bool lineSegmentIntersect(GeoLib::LineSegment const& s0, // no theory for the following tolerances, determined by testing // lower tolerance: little bit smaller than zero - const double l(-1.0*std::numeric_limits<float>::epsilon()); + const double l(-1.0 * std::numeric_limits<float>::epsilon()); // upper tolerance a little bit greater than one - const double u(1.0+std::numeric_limits<float>::epsilon()); - if (rhs[0] < l || u < rhs[0] || rhs[1] < l || u < rhs[1]) { + const double u(1.0 + std::numeric_limits<float>::epsilon()); + if (rhs[0] < l || u < rhs[0] || rhs[1] < l || u < rhs[1]) + { return false; } // compute points along line segments with minimal distance - GeoLib::Point const p0(a[0]+rhs[0]*v[0], a[1]+rhs[0]*v[1], a[2]+rhs[0]*v[2]); - GeoLib::Point const p1(c[0]+rhs[1]*w[0], c[1]+rhs[1]*w[1], c[2]+rhs[1]*w[2]); + GeoLib::Point const p0(a[0] + rhs[0] * v[0], a[1] + rhs[0] * v[1], + a[2] + rhs[0] * v[2]); + GeoLib::Point const p1(c[0] + rhs[1] * w[0], c[1] + rhs[1] * w[1], + c[2] + rhs[1] * w[2]); double const min_dist(std::sqrt(MathLib::sqrDist(p0, p1))); double const min_seg_len( std::min(std::sqrt(sqr_len_v), std::sqrt(sqr_len_w))); - if (min_dist < min_seg_len * 1e-6) { + if (min_dist < min_seg_len * 1e-6) + { s[0] = 0.5 * (p0[0] + p1[0]); s[1] = 0.5 * (p0[1] + p1[1]); s[2] = 0.5 * (p0[2] + p1[2]); @@ -263,8 +266,8 @@ bool lineSegmentIntersect(GeoLib::LineSegment const& s0, } bool lineSegmentsIntersect(const GeoLib::Polyline* ply, - GeoLib::Polyline::SegmentIterator &seg_it0, - GeoLib::Polyline::SegmentIterator &seg_it1, + GeoLib::Polyline::SegmentIterator& seg_it0, + GeoLib::Polyline::SegmentIterator& seg_it1, GeoLib::Point& intersection_pnt) { std::size_t const n_segs(ply->getNumberOfSegments()); @@ -272,13 +275,16 @@ bool lineSegmentsIntersect(const GeoLib::Polyline* ply, // checks for intersections of non-neighbouring segments. for (seg_it0 = ply->begin(); seg_it0 != ply->end() - 2; ++seg_it0) { - seg_it1 = seg_it0+2; + seg_it1 = seg_it0 + 2; std::size_t const seg_num_0 = seg_it0.getSegmentNumber(); - for ( ; seg_it1 != ply->end(); ++seg_it1) { + for (; seg_it1 != ply->end(); ++seg_it1) + { // Do not check first and last segment, because they are // neighboured. - if (!(seg_num_0 == 0 && seg_it1.getSegmentNumber() == n_segs - 1)) { - if (lineSegmentIntersect(*seg_it0, *seg_it1, intersection_pnt)) { + if (!(seg_num_0 == 0 && seg_it1.getSegmentNumber() == n_segs - 1)) + { + if (lineSegmentIntersect(*seg_it0, *seg_it1, intersection_pnt)) + { return true; } } @@ -402,12 +408,11 @@ std::unique_ptr<GeoLib::Point> triangleLineIntersection( u * a[2] + v * b[2] + w * c[2]); } -void computeAndInsertAllIntersectionPoints(GeoLib::PointVec &pnt_vec, - std::vector<GeoLib::Polyline*> & plys) +void computeAndInsertAllIntersectionPoints(GeoLib::PointVec& pnt_vec, + std::vector<GeoLib::Polyline*>& plys) { auto computeSegmentIntersections = [&pnt_vec](GeoLib::Polyline& poly0, - GeoLib::Polyline& poly1) - { + GeoLib::Polyline& poly1) { for (auto seg0_it(poly0.begin()); seg0_it != poly0.end(); ++seg0_it) { for (auto seg1_it(poly1.begin()); seg1_it != poly1.end(); ++seg1_it) @@ -424,10 +429,12 @@ void computeAndInsertAllIntersectionPoints(GeoLib::PointVec &pnt_vec, } }; - for (auto it0(plys.begin()); it0 != plys.end(); ++it0) { + for (auto it0(plys.begin()); it0 != plys.end(); ++it0) + { auto it1(it0); ++it1; - for (; it1 != plys.end(); ++it1) { + for (; it1 != plys.end(); ++it1) + { computeSegmentIntersections(*(*it0), *(*it1)); } } @@ -452,8 +459,7 @@ GeoLib::Polygon rotatePolygonToXY(GeoLib::Polygon const& polygon_in, // 3 set z coord to zero std::for_each(polygon_pnts->begin(), polygon_pnts->end(), - [] (GeoLib::Point* p) { (*p)[2] = 0.0; } - ); + [](GeoLib::Point* p) { (*p)[2] = 0.0; }); // 4 create new polygon GeoLib::Polyline rot_polyline(*polygon_pnts); @@ -477,12 +483,15 @@ std::vector<MathLib::Point3d> lineSegmentIntersect2d( double const orient_abd(getOrientation(a, b, d)); // check if the segment (cd) lies on the left or on the right of (ab) - if ((orient_abc > 0 && orient_abd > 0) || (orient_abc < 0 && orient_abd < 0)) { + if ((orient_abc > 0 && orient_abd > 0) || + (orient_abc < 0 && orient_abd < 0)) + { return std::vector<MathLib::Point3d>(); } // check: (cd) and (ab) are on the same line - if (orient_abc == 0.0 && orient_abd == 0.0) { + if (orient_abc == 0.0 && orient_abd == 0.0) + { double const eps(std::numeric_limits<double>::epsilon()); if (MathLib::sqrDist2d(a, c) < eps && MathLib::sqrDist2d(b, d) < eps) { @@ -495,28 +504,32 @@ std::vector<MathLib::Point3d> lineSegmentIntersect2d( // Since orient_ab and orient_abd vanish, a, b, c, d are on the same // line and for this reason it is enough to check the x-component. - auto isPointOnSegment = [](double q, double p0, double p1) - { + auto isPointOnSegment = [](double q, double p0, double p1) { double const t((q - p0) / (p1 - p0)); return 0 <= t && t <= 1; }; // check if c in (ab) - if (isPointOnSegment(c[0], a[0], b[0])) { + if (isPointOnSegment(c[0], a[0], b[0])) + { // check if a in (cd) - if (isPointOnSegment(a[0], c[0], d[0])) { + if (isPointOnSegment(a[0], c[0], d[0])) + { return {{a, c}}; } // check b == c - if (MathLib::sqrDist2d(b,c) < eps) { + if (MathLib::sqrDist2d(b, c) < eps) + { return {{b}}; } // check if b in (cd) - if (isPointOnSegment(b[0], c[0], d[0])) { + if (isPointOnSegment(b[0], c[0], d[0])) + { return {{b, c}}; } // check d in (ab) - if (isPointOnSegment(d[0], a[0], b[0])) { + if (isPointOnSegment(d[0], a[0], b[0])) + { return {{c, d}}; } std::stringstream err; @@ -529,21 +542,26 @@ std::vector<MathLib::Point3d> lineSegmentIntersect2d( } // check if d in (ab) - if (isPointOnSegment(d[0], a[0], b[0])) { + if (isPointOnSegment(d[0], a[0], b[0])) + { // check if a in (cd) - if (isPointOnSegment(a[0], c[0], d[0])) { + if (isPointOnSegment(a[0], c[0], d[0])) + { return {{a, d}}; } // check if b==d - if (MathLib::sqrDist2d(b, d) < eps) { + if (MathLib::sqrDist2d(b, d) < eps) + { return {{b}}; } // check if b in (cd) - if (isPointOnSegment(b[0], c[0], d[0])) { + if (isPointOnSegment(b[0], c[0], d[0])) + { return {{b, d}}; } // d in (ab), b not in (cd): check c in (ab) - if (isPointOnSegment(c[0], a[0], b[0])) { + if (isPointOnSegment(c[0], a[0], b[0])) + { return {{c, d}}; } @@ -581,7 +599,8 @@ std::vector<MathLib::Point3d> lineSegmentIntersect2d( return false; }; - if (orient_abc == 0.0) { + if (orient_abc == 0.0) + { if (isCollinearPointOntoLineSegment(a, b, c)) { return {{c}}; @@ -589,7 +608,8 @@ std::vector<MathLib::Point3d> lineSegmentIntersect2d( return std::vector<MathLib::Point3d>(); } - if (orient_abd == 0.0) { + if (orient_abd == 0.0) + { if (isCollinearPointOntoLineSegment(a, b, d)) { return {{d}}; @@ -600,7 +620,9 @@ std::vector<MathLib::Point3d> lineSegmentIntersect2d( // check if the segment (ab) lies on the left or on the right of (cd) double const orient_cda(getOrientation(c, d, a)); double const orient_cdb(getOrientation(c, d, b)); - if ((orient_cda > 0 && orient_cdb > 0) || (orient_cda < 0 && orient_cdb < 0)) { + if ((orient_cda > 0 && orient_cdb > 0) || + (orient_cda < 0 && orient_cdb < 0)) + { return std::vector<MathLib::Point3d>(); } @@ -608,63 +630,64 @@ std::vector<MathLib::Point3d> lineSegmentIntersect2d( // linear equations will be invertible // solve the two linear equations (b-a, c-d) (t, s)^T = (c-a) simultaneously Eigen::Matrix2d mat; - mat(0,0) = b[0]-a[0]; - mat(0,1) = c[0]-d[0]; - mat(1,0) = b[1]-a[1]; - mat(1,1) = c[1]-d[1]; + mat(0, 0) = b[0] - a[0]; + mat(0, 1) = c[0] - d[0]; + mat(1, 0) = b[1] - a[1]; + mat(1, 1) = c[1] - d[1]; Eigen::Vector2d rhs{c[0] - a[0], c[1] - a[1]}; rhs = mat.partialPivLu().solve(rhs); - if (0 <= rhs[1] && rhs[1] <= 1.0) { - return { MathLib::Point3d{std::array<double,3>{{ - c[0]+rhs[1]*(d[0]-c[0]), c[1]+rhs[1]*(d[1]-c[1]), - c[2]+rhs[1]*(d[2]-c[2])}} } }; + if (0 <= rhs[1] && rhs[1] <= 1.0) + { + return {MathLib::Point3d{std::array<double, 3>{ + {c[0] + rhs[1] * (d[0] - c[0]), c[1] + rhs[1] * (d[1] - c[1]), + c[2] + rhs[1] * (d[2] - c[2])}}}}; } return std::vector<MathLib::Point3d>(); // parameter s not in the valid // range } -void sortSegments( - MathLib::Point3d const& seg_beg_pnt, - std::vector<GeoLib::LineSegment>& sub_segments) +void sortSegments(MathLib::Point3d const& seg_beg_pnt, + std::vector<GeoLib::LineSegment>& sub_segments) { double const eps(std::numeric_limits<double>::epsilon()); - auto findNextSegment = [&eps]( - MathLib::Point3d const& seg_beg_pnt, - std::vector<GeoLib::LineSegment>& sub_segments, - std::vector<GeoLib::LineSegment>::iterator& - sub_seg_it) { - if (sub_seg_it == sub_segments.end()) - { - return; - } - // find appropriate segment for the given segment begin point - auto act_beg_seg_it = std::find_if( - sub_seg_it, sub_segments.end(), - [&seg_beg_pnt, &eps](GeoLib::LineSegment const& seg) + auto findNextSegment = + [&eps](MathLib::Point3d const& seg_beg_pnt, + std::vector<GeoLib::LineSegment>& sub_segments, + std::vector<GeoLib::LineSegment>::iterator& sub_seg_it) { + if (sub_seg_it == sub_segments.end()) { - return MathLib::sqrDist(seg_beg_pnt, seg.getBeginPoint()) < eps || - MathLib::sqrDist(seg_beg_pnt, seg.getEndPoint()) < eps; - }); - if (act_beg_seg_it == sub_segments.end()) - { - return; - } - // if necessary correct orientation of segment, i.e. swap beg and end - if (MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getEndPoint()) < - MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getBeginPoint())) - { - std::swap(act_beg_seg_it->getBeginPoint(), - act_beg_seg_it->getEndPoint()); - } - assert(sub_seg_it != sub_segments.end()); - // exchange segments within the container - if (sub_seg_it != act_beg_seg_it) - { - std::swap(*sub_seg_it, *act_beg_seg_it); - } - }; + return; + } + // find appropriate segment for the given segment begin point + auto act_beg_seg_it = std::find_if( + sub_seg_it, sub_segments.end(), + [&seg_beg_pnt, &eps](GeoLib::LineSegment const& seg) { + return MathLib::sqrDist(seg_beg_pnt, seg.getBeginPoint()) < + eps || + MathLib::sqrDist(seg_beg_pnt, seg.getEndPoint()) < + eps; + }); + if (act_beg_seg_it == sub_segments.end()) + { + return; + } + // if necessary correct orientation of segment, i.e. swap beg and + // end + if (MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getEndPoint()) < + MathLib::sqrDist(seg_beg_pnt, act_beg_seg_it->getBeginPoint())) + { + std::swap(act_beg_seg_it->getBeginPoint(), + act_beg_seg_it->getEndPoint()); + } + assert(sub_seg_it != sub_segments.end()); + // exchange segments within the container + if (sub_seg_it != act_beg_seg_it) + { + std::swap(*sub_seg_it, *act_beg_seg_it); + } + }; // find start segment auto seg_it = sub_segments.begin(); @@ -672,7 +695,7 @@ void sortSegments( while (seg_it != sub_segments.end()) { - MathLib::Point3d & new_seg_beg_pnt(seg_it->getEndPoint()); + MathLib::Point3d& new_seg_beg_pnt(seg_it->getEndPoint()); seg_it++; if (seg_it != sub_segments.end()) { @@ -686,10 +709,10 @@ Eigen::Matrix3d compute2DRotationMatrixToX(Eigen::Vector3d const& v) Eigen::Matrix3d rot_mat = Eigen::Matrix3d::Zero(); const double cos_theta = v[0]; const double sin_theta = v[1]; - rot_mat(0,0) = rot_mat(1,1) = cos_theta; - rot_mat(0,1) = sin_theta; - rot_mat(1,0) = -sin_theta; - rot_mat(2,2) = 1.0; + rot_mat(0, 0) = rot_mat(1, 1) = cos_theta; + rot_mat(0, 1) = sin_theta; + rot_mat(1, 0) = -sin_theta; + rot_mat(2, 2) = 1.0; return rot_mat; } @@ -733,4 +756,4 @@ Eigen::Matrix3d compute3DRotationMatrixToX(Eigen::Vector3d const& v) return rot_mat; } -} // end namespace GeoLib +} // end namespace GeoLib diff --git a/GeoLib/DuplicateGeometry.cpp b/GeoLib/DuplicateGeometry.cpp index 92313db2e49d88e07ade53a540e07202e080f14f..45d40cb609e7a1ba94447931998457eea57f18cf 100644 --- a/GeoLib/DuplicateGeometry.cpp +++ b/GeoLib/DuplicateGeometry.cpp @@ -12,15 +12,14 @@ #include <map> #include <utility> -#include "BaseLib/Logging.h" +#include "BaseLib/Logging.h" #include "GeoLib/GEOObjects.h" #include "GeoLib/Point.h" #include "GeoLib/Polyline.h" #include "GeoLib/Surface.h" #include "GeoLib/Triangle.h" - namespace GeoLib { DuplicateGeometry::DuplicateGeometry(GeoLib::GEOObjects& geo_objects, @@ -33,7 +32,8 @@ DuplicateGeometry::DuplicateGeometry(GeoLib::GEOObjects& geo_objects, void DuplicateGeometry::duplicate(std::string const& input_name) { - std::vector<GeoLib::Point*> const*const pnts (_geo_objects.getPointVec(input_name)); + std::vector<GeoLib::Point*> const* const pnts( + _geo_objects.getPointVec(input_name)); if (pnts == nullptr) { ERR("Geometry '{:s}' not found.", input_name); @@ -42,7 +42,8 @@ void DuplicateGeometry::duplicate(std::string const& input_name) auto new_pnts = std::make_unique<std::vector<GeoLib::Point*>>(); new_pnts->reserve(pnts->size()); - std::transform(pnts->cbegin(), pnts->cend(), std::back_inserter(*new_pnts), + std::transform( + pnts->cbegin(), pnts->cend(), std::back_inserter(*new_pnts), [](GeoLib::Point* point) { return new GeoLib::Point(*point); }); auto pnt_name_id_map = std::make_unique<std::map<std::string, std::size_t>>( _geo_objects.getPointVecObj(input_name)->getNameIDMapBegin(), @@ -50,7 +51,8 @@ void DuplicateGeometry::duplicate(std::string const& input_name) _geo_objects.addPointVec(std::move(new_pnts), _output_name, std::move(pnt_name_id_map)); - std::vector<GeoLib::Polyline*> const* plys (_geo_objects.getPolylineVec(input_name)); + std::vector<GeoLib::Polyline*> const* plys( + _geo_objects.getPolylineVec(input_name)); if (plys) { auto new_plys = copyPolylinesVector(*plys); @@ -62,7 +64,8 @@ void DuplicateGeometry::duplicate(std::string const& input_name) std::move(ply_name_id_map)); } - std::vector<GeoLib::Surface*> const* sfcs (_geo_objects.getSurfaceVec(input_name)); + std::vector<GeoLib::Surface*> const* sfcs( + _geo_objects.getSurfaceVec(input_name)); if (sfcs) { auto new_sfcs = copySurfacesVector(*sfcs); @@ -75,21 +78,23 @@ void DuplicateGeometry::duplicate(std::string const& input_name) } } -std::unique_ptr<std::vector<GeoLib::Polyline*>> DuplicateGeometry::copyPolylinesVector( +std::unique_ptr<std::vector<GeoLib::Polyline*>> +DuplicateGeometry::copyPolylinesVector( std::vector<GeoLib::Polyline*> const& polylines) const { std::size_t const n_plys = polylines.size(); auto new_lines = std::make_unique<std::vector<GeoLib::Polyline*>>(n_plys, nullptr); - for (std::size_t i=0; i<n_plys; ++i) + for (std::size_t i = 0; i < n_plys; ++i) { if (polylines[i] == nullptr) { continue; } - (*new_lines)[i] = new GeoLib::Polyline(*_geo_objects.getPointVec(_output_name)); - std::size_t const nLinePnts (polylines[i]->getNumberOfPoints()); + (*new_lines)[i] = + new GeoLib::Polyline(*_geo_objects.getPointVec(_output_name)); + std::size_t const nLinePnts(polylines[i]->getNumberOfPoints()); for (std::size_t j = 0; j < nLinePnts; ++j) { (*new_lines)[i]->addPoint(polylines[i]->getPointID(j)); @@ -105,20 +110,22 @@ std::unique_ptr<std::vector<Surface*>> DuplicateGeometry::copySurfacesVector( auto new_surfaces = std::make_unique<std::vector<GeoLib::Surface*>>(n_sfc, nullptr); - for (std::size_t i=0; i<n_sfc; ++i) + for (std::size_t i = 0; i < n_sfc; ++i) { if (surfaces[i] == nullptr) { continue; } - (*new_surfaces)[i] = new GeoLib::Surface(*_geo_objects.getPointVec(_output_name)); + (*new_surfaces)[i] = + new GeoLib::Surface(*_geo_objects.getPointVec(_output_name)); - std::size_t const n_tris (surfaces[i]->getNumberOfTriangles()); - for (std::size_t j=0; j<n_tris; ++j) + std::size_t const n_tris(surfaces[i]->getNumberOfTriangles()); + for (std::size_t j = 0; j < n_tris; ++j) { GeoLib::Triangle const* t = (*surfaces[i])[j]; - (*new_surfaces)[i]->addTriangle( - t->getPoint(0)->getID(), t->getPoint(1)->getID(), t->getPoint(2)->getID()); + (*new_surfaces)[i]->addTriangle(t->getPoint(0)->getID(), + t->getPoint(1)->getID(), + t->getPoint(2)->getID()); } } return new_surfaces; @@ -126,17 +133,20 @@ std::unique_ptr<std::vector<Surface*>> DuplicateGeometry::copySurfacesVector( std::vector<GeoLib::Point*>& DuplicateGeometry::getPointVectorCopy() { - return const_cast<std::vector<GeoLib::Point*>&>(*_geo_objects.getPointVec(_output_name)); + return const_cast<std::vector<GeoLib::Point*>&>( + *_geo_objects.getPointVec(_output_name)); } std::vector<GeoLib::Polyline*>& DuplicateGeometry::getPolylineVectorCopy() { - return const_cast<std::vector<GeoLib::Polyline*>&>(*_geo_objects.getPolylineVec(_output_name)); + return const_cast<std::vector<GeoLib::Polyline*>&>( + *_geo_objects.getPolylineVec(_output_name)); } std::vector<GeoLib::Surface*>& DuplicateGeometry::getSurfaceVectorCopy() { - return const_cast<std::vector<GeoLib::Surface*>&>(*_geo_objects.getSurfaceVec(_output_name)); + return const_cast<std::vector<GeoLib::Surface*>&>( + *_geo_objects.getSurfaceVec(_output_name)); } } // namespace GeoLib diff --git a/GeoLib/GEOObjects.cpp b/GeoLib/GEOObjects.cpp index d45a6eb65a5351fbac3d355e165ce6d4eb2ba524..d831f6b1a05ac6b839062fe439ceaaa97f7a8020 100644 --- a/GeoLib/GEOObjects.cpp +++ b/GeoLib/GEOObjects.cpp @@ -15,9 +15,9 @@ #include "GEOObjects.h" #include <fstream> -#include "BaseLib/StringTools.h" -#include "BaseLib/Logging.h" +#include "BaseLib/Logging.h" +#include "BaseLib/StringTools.h" #include "Triangle.h" namespace GeoLib @@ -50,8 +50,10 @@ void GEOObjects::addPointVec(std::unique_ptr<std::vector<Point*>> points, double eps) { isUniquePointVecName(name); - if (!points || points->empty()) { - DBUG("GEOObjects::addPointVec(): Failed to create PointVec, because " + if (!points || points->empty()) + { + DBUG( + "GEOObjects::addPointVec(): Failed to create PointVec, because " "there aren't any points in the given vector."); return; } @@ -61,7 +63,8 @@ void GEOObjects::addPointVec(std::unique_ptr<std::vector<Point*>> points, _callbacks->addPointVec(name); } -const std::vector<Point*>* GEOObjects::getPointVec(const std::string &name) const +const std::vector<Point*>* GEOObjects::getPointVec( + const std::string& name) const { std::size_t const idx = this->exists(name); if (idx != std::numeric_limits<std::size_t>::max()) @@ -73,7 +76,7 @@ const std::vector<Point*>* GEOObjects::getPointVec(const std::string &name) cons return nullptr; } -const PointVec* GEOObjects::getPointVecObj(const std::string &name) const +const PointVec* GEOObjects::getPointVecObj(const std::string& name) const { std::size_t const idx = this->exists(name); if (idx != std::numeric_limits<std::size_t>::max()) @@ -88,9 +91,11 @@ const PointVec* GEOObjects::getPointVecObj(const std::string &name) const bool GEOObjects::removePointVec(std::string const& name) { - if (isPntVecUsed (name)) + if (isPntVecUsed(name)) { - DBUG("GEOObjects::removePointVec() - There are still Polylines or Surfaces depending on these points."); + DBUG( + "GEOObjects::removePointVec() - There are still Polylines or " + "Surfaces depending on these points."); return false; } @@ -113,7 +118,8 @@ void GEOObjects::addStationVec(std::unique_ptr<std::vector<Point*>> stations, std::string& name) { isUniquePointVecName(name); - _pnt_vecs.push_back(new PointVec(name, std::move(stations), nullptr, PointVec::PointType::STATION)); + _pnt_vecs.push_back(new PointVec(name, std::move(stations), nullptr, + PointVec::PointType::STATION)); _callbacks->addStationVec(name); } @@ -146,7 +152,7 @@ void GEOObjects::addPolylineVec( if ((*it)->getNumberOfPoints() < 2) { auto it_erase(it); - it = lines->erase (it_erase); + it = lines->erase(it_erase); } else { @@ -168,8 +174,8 @@ bool GEOObjects::appendPolylineVec(const std::vector<Polyline*>& polylines, const std::string& name) { // search vector - std::size_t idx (0); - bool nfound (true); + std::size_t idx(0); + bool nfound(true); for (idx = 0; idx < _ply_vecs.size() && nfound; idx++) { if (_ply_vecs[idx]->getName() == name) @@ -181,7 +187,7 @@ bool GEOObjects::appendPolylineVec(const std::vector<Polyline*>& polylines, if (!nfound) { idx--; - std::size_t n_plys (polylines.size()); + std::size_t n_plys(polylines.size()); // append lines for (std::size_t k(0); k < n_plys; k++) { @@ -194,9 +200,10 @@ bool GEOObjects::appendPolylineVec(const std::vector<Polyline*>& polylines, return false; } -const std::vector<Polyline*>* GEOObjects::getPolylineVec(const std::string &name) const +const std::vector<Polyline*>* GEOObjects::getPolylineVec( + const std::string& name) const { - std::size_t size (_ply_vecs.size()); + std::size_t size(_ply_vecs.size()); for (std::size_t i = 0; i < size; i++) { if (_ply_vecs[i]->getName() == name) @@ -210,9 +217,9 @@ const std::vector<Polyline*>* GEOObjects::getPolylineVec(const std::string &name return nullptr; } -const PolylineVec* GEOObjects::getPolylineVecObj(const std::string &name) const +const PolylineVec* GEOObjects::getPolylineVecObj(const std::string& name) const { - std::size_t size (_ply_vecs.size()); + std::size_t size(_ply_vecs.size()); for (std::size_t i = 0; i < size; i++) { if (_ply_vecs[i]->getName() == name) @@ -259,8 +266,8 @@ bool GEOObjects::appendSurfaceVec(const std::vector<Surface*>& surfaces, const std::string& name) { // search vector - std::size_t idx (0); - bool nfound (true); + std::size_t idx(0); + bool nfound(true); for (idx = 0; idx < _sfc_vecs.size() && nfound; idx++) { if (_sfc_vecs[idx]->getName() == name) @@ -272,7 +279,7 @@ bool GEOObjects::appendSurfaceVec(const std::vector<Surface*>& surfaces, if (!nfound) { idx--; - std::size_t n_sfcs (surfaces.size()); + std::size_t n_sfcs(surfaces.size()); // append surfaces for (std::size_t k(0); k < n_sfcs; k++) { @@ -290,9 +297,10 @@ bool GEOObjects::appendSurfaceVec(const std::vector<Surface*>& surfaces, return false; } -const std::vector<Surface*>* GEOObjects::getSurfaceVec(const std::string &name) const +const std::vector<Surface*>* GEOObjects::getSurfaceVec( + const std::string& name) const { - std::size_t size (_sfc_vecs.size()); + std::size_t size(_sfc_vecs.size()); for (std::size_t i = 0; i < size; i++) { if (_sfc_vecs[i]->getName() == name) @@ -305,7 +313,7 @@ const std::vector<Surface*>* GEOObjects::getSurfaceVec(const std::string &name) return nullptr; } -bool GEOObjects::removeSurfaceVec(const std::string &name) +bool GEOObjects::removeSurfaceVec(const std::string& name) { _callbacks->removeSurfaceVec(name); for (auto it(_sfc_vecs.begin()); it != _sfc_vecs.end(); ++it) @@ -313,7 +321,7 @@ bool GEOObjects::removeSurfaceVec(const std::string &name) if ((*it)->getName() == name) { delete *it; - _sfc_vecs.erase (it); + _sfc_vecs.erase(it); return true; } } @@ -323,9 +331,9 @@ bool GEOObjects::removeSurfaceVec(const std::string &name) return false; } -const SurfaceVec* GEOObjects::getSurfaceVecObj(const std::string &name) const +const SurfaceVec* GEOObjects::getSurfaceVecObj(const std::string& name) const { - std::size_t size (_sfc_vecs.size()); + std::size_t size(_sfc_vecs.size()); for (std::size_t i = 0; i < size; i++) { if (_sfc_vecs[i]->getName() == name) @@ -351,7 +359,7 @@ bool GEOObjects::isUniquePointVecName(std::string& name) const return true; } -bool GEOObjects::isPntVecUsed (const std::string &name) const +bool GEOObjects::isPntVecUsed(const std::string& name) const { // search dependent data structures (Polyline) for (auto polyline : _ply_vecs) @@ -407,10 +415,12 @@ std::string GEOObjects::getElementNameByID(const std::string& geometry_name, this->getPointVecObj(geometry_name)->getNameOfElementByID(id, name); break; case GeoLib::GEOTYPE::POLYLINE: - this->getPolylineVecObj(geometry_name)->getNameOfElementByID(id, name); + this->getPolylineVecObj(geometry_name) + ->getNameOfElementByID(id, name); break; case GeoLib::GEOTYPE::SURFACE: - this->getSurfaceVecObj(geometry_name)->getNameOfElementByID(id, name); + this->getSurfaceVecObj(geometry_name) + ->getNameOfElementByID(id, name); } return name; } @@ -439,8 +449,9 @@ int GEOObjects::mergeGeometries(std::vector<std::string> const& geo_names, return 0; } -bool GEOObjects::mergePoints(std::vector<std::string> const & geo_names, - std::string & merged_geo_name, std::vector<std::size_t> &pnt_offsets) +bool GEOObjects::mergePoints(std::vector<std::string> const& geo_names, + std::string& merged_geo_name, + std::vector<std::size_t>& pnt_offsets) { const std::size_t n_geo_names(geo_names.size()); @@ -448,33 +459,40 @@ bool GEOObjects::mergePoints(std::vector<std::string> const & geo_names, auto merged_pnt_names = std::make_unique<std::map<std::string, std::size_t>>(); - for (std::size_t j(0); j < n_geo_names; ++j) { - GeoLib::PointVec const*const pnt_vec(this->getPointVecObj(geo_names[j])); + for (std::size_t j(0); j < n_geo_names; ++j) + { + GeoLib::PointVec const* const pnt_vec( + this->getPointVecObj(geo_names[j])); if (pnt_vec == nullptr) { continue; } const std::vector<GeoLib::Point*>* pnts(pnt_vec->getVector()); - if (pnts == nullptr) { + if (pnts == nullptr) + { return false; } // do not consider stations - if (dynamic_cast<GeoLib::Station*>((*pnts)[0])) { + if (dynamic_cast<GeoLib::Station*>((*pnts)[0])) + { continue; } std::size_t const n_pnts(pnts->size()); - for (std::size_t k(0); k < n_pnts; ++k) { + for (std::size_t k(0); k < n_pnts; ++k) + { merged_points->push_back( - new GeoLib::Point(*(*pnts)[k], pnt_offsets[j]+k)); + new GeoLib::Point(*(*pnts)[k], pnt_offsets[j] + k)); std::string const& item_name(pnt_vec->getItemNameByID(k)); - if (! item_name.empty()) { + if (!item_name.empty()) + { merged_pnt_names->insert( std::make_pair(item_name, pnt_offsets[j] + k)); } } - if (n_geo_names - 1 > j) { + if (n_geo_names - 1 > j) + { pnt_offsets[j + 1] = n_pnts + pnt_offsets[j]; } } @@ -495,34 +513,49 @@ void GEOObjects::mergePolylines(std::vector<std::string> const& geo_names, auto merged_ply_names = std::make_unique<std::map<std::string, std::size_t>>(); - std::vector<GeoLib::Point*> const* merged_points(this->getPointVecObj(merged_geo_name)->getVector()); - std::vector<std::size_t> const& id_map (this->getPointVecObj(merged_geo_name)->getIDMap ()); + std::vector<GeoLib::Point*> const* merged_points( + this->getPointVecObj(merged_geo_name)->getVector()); + std::vector<std::size_t> const& id_map( + this->getPointVecObj(merged_geo_name)->getIDMap()); - for (std::size_t j(0); j < n_geo_names; j++) { - const std::vector<GeoLib::Polyline*>* plys (this->getPolylineVec(geo_names[j])); - if (plys) { + for (std::size_t j(0); j < n_geo_names; j++) + { + const std::vector<GeoLib::Polyline*>* plys( + this->getPolylineVec(geo_names[j])); + if (plys) + { std::string tmp_name; - for (std::size_t k(0); k < plys->size(); k++) { + for (std::size_t k(0); k < plys->size(); k++) + { auto* kth_ply_new(new GeoLib::Polyline(*merged_points)); - GeoLib::Polyline const* const kth_ply_old ((*plys)[k]); - const std::size_t size_of_kth_ply (kth_ply_old->getNumberOfPoints()); - // copy point ids from old ply to new ply (considering the offset) - for (std::size_t i(0); i < size_of_kth_ply; i++) { - kth_ply_new->addPoint (id_map[pnt_offsets[j] + - kth_ply_old->getPointID(i)]); + GeoLib::Polyline const* const kth_ply_old((*plys)[k]); + const std::size_t size_of_kth_ply( + kth_ply_old->getNumberOfPoints()); + // copy point ids from old ply to new ply (considering the + // offset) + for (std::size_t i(0); i < size_of_kth_ply; i++) + { + kth_ply_new->addPoint( + id_map[pnt_offsets[j] + kth_ply_old->getPointID(i)]); } - merged_polylines->push_back (kth_ply_new); - if (this->getPolylineVecObj(geo_names[j])->getNameOfElementByID(k, tmp_name)) { - merged_ply_names->insert(std::pair<std::string, std::size_t>(tmp_name, ply_offsets[j] + k)); + merged_polylines->push_back(kth_ply_new); + if (this->getPolylineVecObj(geo_names[j]) + ->getNameOfElementByID(k, tmp_name)) + { + merged_ply_names->insert( + std::pair<std::string, std::size_t>( + tmp_name, ply_offsets[j] + k)); } } - if (n_geo_names - 1 > j) { + if (n_geo_names - 1 > j) + { ply_offsets[j + 1] = plys->size() + ply_offsets[j]; } } } - if (! merged_polylines->empty()) { + if (!merged_polylines->empty()) + { this->addPolylineVec(std::move(merged_polylines), merged_geo_name, std::move(merged_ply_names)); } @@ -532,42 +565,56 @@ void GEOObjects::mergeSurfaces(std::vector<std::string> const& geo_names, std::string const& merged_geo_name, std::vector<std::size_t> const& pnt_offsets) { - std::vector<GeoLib::Point*> const* merged_points(this->getPointVecObj(merged_geo_name)->getVector()); - std::vector<std::size_t> const& id_map (this->getPointVecObj(merged_geo_name)->getIDMap ()); + std::vector<GeoLib::Point*> const* merged_points( + this->getPointVecObj(merged_geo_name)->getVector()); + std::vector<std::size_t> const& id_map( + this->getPointVecObj(merged_geo_name)->getIDMap()); const std::size_t n_geo_names(geo_names.size()); std::vector<std::size_t> sfc_offsets(n_geo_names, 0); auto merged_sfcs = std::make_unique<std::vector<GeoLib::Surface*>>(); auto merged_sfc_names = std::make_unique<std::map<std::string, std::size_t>>(); - for (std::size_t j(0); j < n_geo_names; j++) { - const std::vector<GeoLib::Surface*>* sfcs (this->getSurfaceVec(geo_names[j])); - if (sfcs) { + for (std::size_t j(0); j < n_geo_names; j++) + { + const std::vector<GeoLib::Surface*>* sfcs( + this->getSurfaceVec(geo_names[j])); + if (sfcs) + { std::string tmp_name; - for (std::size_t k(0); k < sfcs->size(); k++) { + for (std::size_t k(0); k < sfcs->size(); k++) + { auto* kth_sfc_new(new GeoLib::Surface(*merged_points)); - GeoLib::Surface const* const kth_sfc_old ((*sfcs)[k]); - const std::size_t size_of_kth_sfc (kth_sfc_old->getNumberOfTriangles()); + GeoLib::Surface const* const kth_sfc_old((*sfcs)[k]); + const std::size_t size_of_kth_sfc( + kth_sfc_old->getNumberOfTriangles()); // clone surface elements using new ids - for (std::size_t i(0); i < size_of_kth_sfc; i++) { - const GeoLib::Triangle* tri ((*kth_sfc_old)[i]); - const std::size_t id0 (id_map[pnt_offsets[j] + (*tri)[0]]); - const std::size_t id1 (id_map[pnt_offsets[j] + (*tri)[1]]); - const std::size_t id2 (id_map[pnt_offsets[j] + (*tri)[2]]); - kth_sfc_new->addTriangle (id0, id1, id2); + for (std::size_t i(0); i < size_of_kth_sfc; i++) + { + const GeoLib::Triangle* tri((*kth_sfc_old)[i]); + const std::size_t id0(id_map[pnt_offsets[j] + (*tri)[0]]); + const std::size_t id1(id_map[pnt_offsets[j] + (*tri)[1]]); + const std::size_t id2(id_map[pnt_offsets[j] + (*tri)[2]]); + kth_sfc_new->addTriangle(id0, id1, id2); } - merged_sfcs->push_back (kth_sfc_new); - - if (this->getSurfaceVecObj(geo_names[j])->getNameOfElementByID(k, tmp_name)) { - merged_sfc_names->insert(std::pair<std::string, std::size_t>(tmp_name, sfc_offsets[j] + k)); + merged_sfcs->push_back(kth_sfc_new); + + if (this->getSurfaceVecObj(geo_names[j]) + ->getNameOfElementByID(k, tmp_name)) + { + merged_sfc_names->insert( + std::pair<std::string, std::size_t>( + tmp_name, sfc_offsets[j] + k)); } } - if (n_geo_names - 1 > j) { + if (n_geo_names - 1 > j) + { sfc_offsets[j + 1] = sfcs->size() + sfc_offsets[j]; } } } - if (! merged_sfcs->empty()) { + if (!merged_sfcs->empty()) + { this->addSurfaceVec(std::move(merged_sfcs), merged_geo_name, std::move(merged_sfc_names)); } @@ -692,44 +739,49 @@ const GeoLib::GeoObject* GEOObjects::getGeoObject( GeoLib::GEOTYPE type, const std::string& geo_obj_name) const { - GeoLib::GeoObject *geo_obj(nullptr); - switch (type) { - case GeoLib::GEOTYPE::POINT: { - GeoLib::PointVec const* pnt_vec(getPointVecObj(geo_name)); - if (pnt_vec) + GeoLib::GeoObject* geo_obj(nullptr); + switch (type) + { + case GeoLib::GEOTYPE::POINT: { - geo_obj = const_cast<GeoLib::GeoObject*>( - dynamic_cast<GeoLib::GeoObject const*>( - pnt_vec->getElementByName(geo_obj_name))); + GeoLib::PointVec const* pnt_vec(getPointVecObj(geo_name)); + if (pnt_vec) + { + geo_obj = const_cast<GeoLib::GeoObject*>( + dynamic_cast<GeoLib::GeoObject const*>( + pnt_vec->getElementByName(geo_obj_name))); + } + break; } - break; - } - case GeoLib::GEOTYPE::POLYLINE: { - GeoLib::PolylineVec const* ply_vec(getPolylineVecObj(geo_name)); - if (ply_vec) + case GeoLib::GEOTYPE::POLYLINE: { - geo_obj = const_cast<GeoLib::GeoObject*>( - dynamic_cast<GeoLib::GeoObject const*>( - ply_vec->getElementByName(geo_obj_name))); + GeoLib::PolylineVec const* ply_vec(getPolylineVecObj(geo_name)); + if (ply_vec) + { + geo_obj = const_cast<GeoLib::GeoObject*>( + dynamic_cast<GeoLib::GeoObject const*>( + ply_vec->getElementByName(geo_obj_name))); + } + break; } - break; - } - case GeoLib::GEOTYPE::SURFACE: { - GeoLib::SurfaceVec const* sfc_vec(getSurfaceVecObj(geo_name)); - if (sfc_vec) + case GeoLib::GEOTYPE::SURFACE: { - geo_obj = const_cast<GeoLib::GeoObject*>( - dynamic_cast<GeoLib::GeoObject const*>( - sfc_vec->getElementByName(geo_obj_name))); + GeoLib::SurfaceVec const* sfc_vec(getSurfaceVecObj(geo_name)); + if (sfc_vec) + { + geo_obj = const_cast<GeoLib::GeoObject*>( + dynamic_cast<GeoLib::GeoObject const*>( + sfc_vec->getElementByName(geo_obj_name))); + } + break; } - break; - } - default: - ERR("GEOObjects::getGeoObject(): geometric type not handled."); - return nullptr; + default: + ERR("GEOObjects::getGeoObject(): geometric type not handled."); + return nullptr; }; - if (!geo_obj) { + if (!geo_obj) + { DBUG( "GEOObjects::getGeoObject(): Could not find {:s} '{:s}' in " "geometry.", @@ -740,12 +792,10 @@ const GeoLib::GeoObject* GEOObjects::getGeoObject( } GeoLib::GeoObject const* GEOObjects::getGeoObject( - const std::string &geo_name, - const std::string &geo_obj_name) const + const std::string& geo_name, const std::string& geo_obj_name) const { GeoLib::GeoObject const* geo_obj( - getGeoObject(geo_name, GeoLib::GEOTYPE::POINT, geo_obj_name) - ); + getGeoObject(geo_name, GeoLib::GEOTYPE::POINT, geo_obj_name)); if (!geo_obj) { @@ -759,7 +809,8 @@ GeoLib::GeoObject const* GEOObjects::getGeoObject( getGeoObject(geo_name, GeoLib::GEOTYPE::SURFACE, geo_obj_name); } - if (!geo_obj) { + if (!geo_obj) + { DBUG( "GEOObjects::getGeoObject(): Could not find '{:s}' in geometry " "{:s}.", @@ -768,9 +819,9 @@ GeoLib::GeoObject const* GEOObjects::getGeoObject( return geo_obj; } -std::size_t GEOObjects::exists(const std::string &geometry_name) const +std::size_t GEOObjects::exists(const std::string& geometry_name) const { - std::size_t const size (_pnt_vecs.size()); + std::size_t const size(_pnt_vecs.size()); for (std::size_t i = 0; i < size; i++) { if (_pnt_vecs[i]->getName() == geometry_name) @@ -779,7 +830,8 @@ std::size_t GEOObjects::exists(const std::string &geometry_name) const } } - // HACK for enabling conversion of files without loading the associated geometry + // HACK for enabling conversion of files without loading the associated + // geometry if (size > 0 && _pnt_vecs[0]->getName() == "conversionTestRun#1") { return 1; diff --git a/GeoLib/GeoType.cpp b/GeoLib/GeoType.cpp index 21d82576304c597a6ed28a50f54ac1378f8f9f93..2d4bb326c056d61ac033000464e1164fa29c6f56 100644 --- a/GeoLib/GeoType.cpp +++ b/GeoLib/GeoType.cpp @@ -20,14 +20,16 @@ namespace GeoLib { - -std::string convertGeoTypeToString (GEOTYPE geo_type) +std::string convertGeoTypeToString(GEOTYPE geo_type) { switch (geo_type) { - case GEOTYPE::POINT: return "POINT"; - case GEOTYPE::POLYLINE: return "POLYLINE"; - case GEOTYPE::SURFACE: return "SURFACE"; + case GEOTYPE::POINT: + return "POINT"; + case GEOTYPE::POLYLINE: + return "POLYLINE"; + case GEOTYPE::SURFACE: + return "SURFACE"; } // Cannot happen, because switch covers all cases. @@ -35,4 +37,4 @@ std::string convertGeoTypeToString (GEOTYPE geo_type) OGS_FATAL("convertGeoTypeToString(): Given geo type is not supported"); } -} // end namespace GeoLib +} // end namespace GeoLib diff --git a/GeoLib/IO/TINInterface.cpp b/GeoLib/IO/TINInterface.cpp index 11535e56263d4107d867431e6a842b113579decb..dd89c6952889cb88d1dc91092471703902d3df4a 100644 --- a/GeoLib/IO/TINInterface.cpp +++ b/GeoLib/IO/TINInterface.cpp @@ -11,29 +11,26 @@ #include <fstream> #include <limits> -#include "BaseLib/Logging.h" - #include "BaseLib/FileTools.h" +#include "BaseLib/Logging.h" #include "BaseLib/StringTools.h" - #include "GeoLib/AnalyticalGeometry.h" #include "GeoLib/Surface.h" #include "GeoLib/Triangle.h" - #include "MathLib/GeometricBasics.h" namespace GeoLib { namespace IO { - GeoLib::Surface* TINInterface::readTIN(std::string const& fname, - GeoLib::PointVec &pnt_vec, - std::vector<std::string>* errors) + GeoLib::PointVec& pnt_vec, + std::vector<std::string>* errors) { // open file std::ifstream in(fname.c_str()); - if (!in) { + if (!in) + { WARN("readTIN(): could not open stream from {:s}.", fname); if (errors) { @@ -59,13 +56,15 @@ GeoLib::Surface* TINInterface::readTIN(std::string const& fname, // parse line std::stringstream input(line); // read id - if (!(input >> id)) { + if (!(input >> id)) + { in.close(); delete sfc; return nullptr; } // read first point - if (!(input >> p0[0] >> p0[1] >> p0[2])) { + if (!(input >> p0[0] >> p0[1] >> p0[2])) + { ERR("Could not read coords of 1st point of triangle {:d}.", id); if (errors) { @@ -80,7 +79,8 @@ GeoLib::Surface* TINInterface::readTIN(std::string const& fname, return nullptr; } // read second point - if (!(input >> p1[0] >> p1[1] >> p1[2])) { + if (!(input >> p1[0] >> p1[1] >> p1[2])) + { ERR("Could not read coords of 2nd point of triangle {:d}.", id); if (errors) { @@ -95,7 +95,8 @@ GeoLib::Surface* TINInterface::readTIN(std::string const& fname, return nullptr; } // read third point - if (!(input >> p2[0] >> p2[1] >> p2[2])) { + if (!(input >> p2[0] >> p2[1] >> p2[2])) + { ERR("Could not read coords of 3rd point of triangle {:d}.", id); if (errors) { @@ -112,7 +113,8 @@ GeoLib::Surface* TINInterface::readTIN(std::string const& fname, // check area of triangle double const d_eps(std::numeric_limits<double>::epsilon()); - if (MathLib::calcTriangleArea(p0, p1, p2) < d_eps) { + if (MathLib::calcTriangleArea(p0, p1, p2) < d_eps) + { ERR("readTIN: Triangle {:d} has zero area.", id); if (errors) { @@ -127,9 +129,12 @@ GeoLib::Surface* TINInterface::readTIN(std::string const& fname, // determine size pnt_vec to insert the correct ids std::size_t const s(pnt_vec.getVector()->size()); - std::size_t const pnt_pos_0(pnt_vec.push_back(new GeoLib::Point(p0, s))); - std::size_t const pnt_pos_1(pnt_vec.push_back(new GeoLib::Point(p1, s+1))); - std::size_t const pnt_pos_2(pnt_vec.push_back(new GeoLib::Point(p2, s+2))); + std::size_t const pnt_pos_0( + pnt_vec.push_back(new GeoLib::Point(p0, s))); + std::size_t const pnt_pos_1( + pnt_vec.push_back(new GeoLib::Point(p1, s + 1))); + std::size_t const pnt_pos_2( + pnt_vec.push_back(new GeoLib::Point(p2, s + 2))); // create new Triangle if (pnt_pos_0 != std::numeric_limits<std::size_t>::max() && pnt_pos_1 != std::numeric_limits<std::size_t>::max() && @@ -139,7 +144,8 @@ GeoLib::Surface* TINInterface::readTIN(std::string const& fname, } } - if (sfc->getNumberOfTriangles() == 0) { + if (sfc->getNumberOfTriangles() == 0) + { WARN("readTIN(): No triangle found.", fname); if (errors) { @@ -152,20 +158,24 @@ GeoLib::Surface* TINInterface::readTIN(std::string const& fname, return sfc; } -void TINInterface::writeSurfaceAsTIN(GeoLib::Surface const& surface, std::string const& file_name) +void TINInterface::writeSurfaceAsTIN(GeoLib::Surface const& surface, + std::string const& file_name) { - std::ofstream os (file_name.c_str()); - if (!os) { + std::ofstream os(file_name.c_str()); + if (!os) + { WARN("writeSurfaceAsTIN(): could not open stream to {:s}.", file_name); return; } os.precision(std::numeric_limits<double>::digits10); - const std::size_t n_tris (surface.getNumberOfTriangles()); - for (std::size_t l(0); l < n_tris; l++) { - GeoLib::Triangle const& tri (*(surface[l])); - os << l << " " << *(tri.getPoint(0)) << " " << *(tri.getPoint(1)) << " " << *(tri.getPoint(2)) << "\n"; + const std::size_t n_tris(surface.getNumberOfTriangles()); + for (std::size_t l(0); l < n_tris; l++) + { + GeoLib::Triangle const& tri(*(surface[l])); + os << l << " " << *(tri.getPoint(0)) << " " << *(tri.getPoint(1)) << " " + << *(tri.getPoint(2)) << "\n"; } os.close(); } -} // end namespace IO -} // end namespace GeoLib +} // end namespace IO +} // end namespace GeoLib diff --git a/GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.cpp b/GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.cpp index 93c4c7ccd3153d90d87b3c7a055c2c770f3ec09b..a47b433d80d56c1b1abece4b0399fcb750d2bedd 100644 --- a/GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.cpp +++ b/GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.cpp @@ -15,10 +15,10 @@ #include "BoostXmlGmlInterface.h" #include <boost/property_tree/xml_parser.hpp> -#include "BaseLib/Logging.h" #include "BaseLib/Algorithm.h" #include "BaseLib/ConfigTreeUtil.h" +#include "BaseLib/Logging.h" #include "GeoLib/GEOObjects.h" #include "GeoLib/Point.h" #include "GeoLib/PointVec.h" @@ -30,12 +30,12 @@ namespace GeoLib { namespace IO { +BoostXmlGmlInterface::BoostXmlGmlInterface(GeoLib::GEOObjects& geo_objs) + : _geo_objects(geo_objs) +{ +} -BoostXmlGmlInterface::BoostXmlGmlInterface(GeoLib::GEOObjects& geo_objs) : - _geo_objects(geo_objs) -{} - -bool BoostXmlGmlInterface::readFile(const std::string &fname) +bool BoostXmlGmlInterface::readFile(const std::string& fname) { //! \todo Reading geometries is always strict. auto doc = BaseLib::makeConfigTree(fname, true, "OpenGeoSysGLI"); @@ -89,12 +89,14 @@ bool BoostXmlGmlInterface::readFile(const std::string &fname) *sfc_names); } - if (!polylines->empty()) { + if (!polylines->empty()) + { _geo_objects.addPolylineVec(std::move(polylines), geo_name, std::move(ply_names)); } - if (!surfaces->empty()) { + if (!surfaces->empty()) + { _geo_objects.addSurfaceVec(std::move(surfaces), geo_name, std::move(sfc_names)); } @@ -102,9 +104,10 @@ bool BoostXmlGmlInterface::readFile(const std::string &fname) return true; } -void BoostXmlGmlInterface::readPoints(BaseLib::ConfigTree const& pointsRoot, - std::vector<GeoLib::Point*>& points, - std::map<std::string, std::size_t>& pnt_names ) +void BoostXmlGmlInterface::readPoints( + BaseLib::ConfigTree const& pointsRoot, + std::vector<GeoLib::Point*>& points, + std::map<std::string, std::size_t>& pnt_names) { //! \ogs_file_param{gml__points__point} for (auto const pt : pointsRoot.getConfigParameterList("point")) @@ -112,26 +115,28 @@ void BoostXmlGmlInterface::readPoints(BaseLib::ConfigTree const& pointsRoot, //! \ogs_file_attr{gml__points__point__id} auto const p_id = pt.getConfigAttribute<std::size_t>("id"); //! \ogs_file_attr{gml__points__point__x} - auto const p_x = pt.getConfigAttribute<double>("x"); + auto const p_x = pt.getConfigAttribute<double>("x"); //! \ogs_file_attr{gml__points__point__y} - auto const p_y = pt.getConfigAttribute<double>("y"); + auto const p_y = pt.getConfigAttribute<double>("y"); //! \ogs_file_attr{gml__points__point__z} - auto const p_z = pt.getConfigAttribute<double>("z"); + auto const p_z = pt.getConfigAttribute<double>("z"); auto const p_size = points.size(); BaseLib::insertIfKeyUniqueElseError(_idx_map, p_id, p_size, - "The point id is not unique."); + "The point id is not unique."); points.push_back(new GeoLib::Point(p_x, p_y, p_z, p_id)); //! \ogs_file_attr{gml__points__point__name} - if (auto const p_name = pt.getConfigAttributeOptional<std::string>("name")) + if (auto const p_name = + pt.getConfigAttributeOptional<std::string>("name")) { - if (p_name->empty()) { + if (p_name->empty()) + { OGS_FATAL("Empty point name found in geometry file."); } - BaseLib::insertIfKeyUniqueElseError(pnt_names, *p_name, p_size, - "The point name is not unique."); + BaseLib::insertIfKeyUniqueElseError( + pnt_names, *p_name, p_size, "The point name is not unique."); } } } @@ -150,18 +155,21 @@ void BoostXmlGmlInterface::readPolylines( auto const id = pl.getConfigAttribute<std::size_t>("id"); // The id is not used but must be present in the GML file. // That's why pl.ignore...() cannot be used. - (void) id; + (void)id; polylines.push_back(new GeoLib::Polyline(points)); //! \ogs_file_attr{gml__polylines__polyline__name} - if (auto const p_name = pl.getConfigAttributeOptional<std::string>("name")) + if (auto const p_name = + pl.getConfigAttributeOptional<std::string>("name")) { - if (p_name->empty()) { + if (p_name->empty()) + { OGS_FATAL("Empty polyline name found in geometry file."); } - BaseLib::insertIfKeyUniqueElseError(ply_names, *p_name, polylines.size()-1, + BaseLib::insertIfKeyUniqueElseError( + ply_names, *p_name, polylines.size() - 1, "The polyline name is not unique."); auto accessOrError = [this, &p_name](auto pt_idx) { @@ -191,7 +199,7 @@ void BoostXmlGmlInterface::readPolylines( } void BoostXmlGmlInterface::readSurfaces( - BaseLib::ConfigTree const& surfacesRoot, + BaseLib::ConfigTree const& surfacesRoot, std::vector<GeoLib::Surface*>& surfaces, std::vector<GeoLib::Point*> const& points, const std::vector<std::size_t>& pnt_id_map, @@ -204,27 +212,34 @@ void BoostXmlGmlInterface::readSurfaces( auto const id = sfc.getConfigAttribute<std::size_t>("id"); // The id is not used but must be present in the GML file. // That's why sfc.ignore...() cannot be used. - (void) id; + (void)id; surfaces.push_back(new GeoLib::Surface(points)); //! \ogs_file_attr{gml__surfaces__surface__name} - if (auto const s_name = sfc.getConfigAttributeOptional<std::string>("name")) + if (auto const s_name = + sfc.getConfigAttributeOptional<std::string>("name")) { - if (s_name->empty()) { + if (s_name->empty()) + { OGS_FATAL("Empty surface name found in geometry file."); } - BaseLib::insertIfKeyUniqueElseError(sfc_names, *s_name, surfaces.size()-1, + BaseLib::insertIfKeyUniqueElseError( + sfc_names, *s_name, surfaces.size() - 1, "The surface name is not unique."); //! \ogs_file_param{gml__surfaces__surface__element} - for (auto const& element : sfc.getConfigParameterList("element")) { + for (auto const& element : sfc.getConfigParameterList("element")) + { //! \ogs_file_attr{gml__surfaces__surface__element__p1} - auto const p1_attr = element.getConfigAttribute<std::size_t>("p1"); + auto const p1_attr = + element.getConfigAttribute<std::size_t>("p1"); //! \ogs_file_attr{gml__surfaces__surface__element__p2} - auto const p2_attr = element.getConfigAttribute<std::size_t>("p2"); + auto const p2_attr = + element.getConfigAttribute<std::size_t>("p2"); //! \ogs_file_attr{gml__surfaces__surface__element__p3} - auto const p3_attr = element.getConfigAttribute<std::size_t>("p3"); + auto const p3_attr = + element.getConfigAttribute<std::size_t>("p3"); auto accessOrError = [this, &s_name](std::size_t pt_idx) { auto search = _idx_map.find(pt_idx); @@ -242,7 +257,7 @@ void BoostXmlGmlInterface::readSurfaces( auto const p1 = pnt_id_map[accessOrError(p1_attr)]; auto const p2 = pnt_id_map[accessOrError(p2_attr)]; auto const p3 = pnt_id_map[accessOrError(p3_attr)]; - surfaces.back()->addTriangle(p1,p2,p3); + surfaces.back()->addTriangle(p1, p2, p3); } } else @@ -263,21 +278,24 @@ bool BoostXmlGmlInterface::write() GeoLib::PointVec const* const pnt_vec( _geo_objects.getPointVecObj(export_name)); - if (! pnt_vec) { + if (!pnt_vec) + { ERR("BoostXmlGmlInterface::write(): No PointVec within the geometry " "'{:s}'.", export_name); return false; } - std::vector<GeoLib::Point*> const*const pnts(pnt_vec->getVector()); - if (! pnts) { + std::vector<GeoLib::Point*> const* const pnts(pnt_vec->getVector()); + if (!pnts) + { ERR("BoostXmlGmlInterface::write(): No vector of points within the " "geometry '{:s}'.", export_name); return false; } - if (pnts->empty()) { + if (pnts->empty()) + { ERR("BoostXmlGmlInterface::write(): No points within the geometry " "'{:s}'.", export_name); @@ -294,7 +312,8 @@ bool BoostXmlGmlInterface::write() geometry_set.add("name", export_name); auto& pnts_tag = geometry_set.add("points", ""); - for (std::size_t k(0); k<pnts->size(); k++) { + for (std::size_t k(0); k < pnts->size(); k++) + { auto& pnt_tag = pnts_tag.add("point", ""); pnt_tag.put("<xmlattr>.id", k); pnt_tag.put("<xmlattr>.x", (*((*pnts)[k]))[0]); @@ -316,11 +335,12 @@ bool BoostXmlGmlInterface::write() } void BoostXmlGmlInterface::addSurfacesToPropertyTree( - boost::property_tree::ptree & geometry_set) + boost::property_tree::ptree& geometry_set) { GeoLib::SurfaceVec const* const sfc_vec( _geo_objects.getSurfaceVecObj(export_name)); - if (!sfc_vec) { + if (!sfc_vec) + { INFO( "BoostXmlGmlInterface::addSurfacesToPropertyTree(): " "No surfaces within the geometry '{:s}'.", @@ -328,7 +348,7 @@ void BoostXmlGmlInterface::addSurfacesToPropertyTree( return; } - std::vector<GeoLib::Surface*> const*const surfaces(sfc_vec->getVector()); + std::vector<GeoLib::Surface*> const* const surfaces(sfc_vec->getVector()); if (!surfaces || surfaces->empty()) { INFO( @@ -339,8 +359,9 @@ void BoostXmlGmlInterface::addSurfacesToPropertyTree( } auto& surfaces_tag = geometry_set.add("surfaces", ""); - for (std::size_t i=0; i<surfaces->size(); ++i) { - GeoLib::Surface const*const surface((*surfaces)[i]); + for (std::size_t i = 0; i < surfaces->size(); ++i) + { + GeoLib::Surface const* const surface((*surfaces)[i]); std::string sfc_name; sfc_vec->getNameOfElement(surface, sfc_name); auto& surface_tag = surfaces_tag.add("surface", ""); @@ -349,7 +370,8 @@ void BoostXmlGmlInterface::addSurfacesToPropertyTree( { surface_tag.put("<xmlattr>.name", sfc_name); } - for (std::size_t j=0; j<surface->getNumberOfTriangles(); ++j) { + for (std::size_t j = 0; j < surface->getNumberOfTriangles(); ++j) + { auto& element_tag = surface_tag.add("element", ""); element_tag.put("<xmlattr>.p1", (*(*surface)[j])[0]); element_tag.put("<xmlattr>.p2", (*(*surface)[j])[1]); @@ -359,11 +381,12 @@ void BoostXmlGmlInterface::addSurfacesToPropertyTree( } void BoostXmlGmlInterface::addPolylinesToPropertyTree( - boost::property_tree::ptree & geometry_set) + boost::property_tree::ptree& geometry_set) { GeoLib::PolylineVec const* const vec( _geo_objects.getPolylineVecObj(export_name)); - if (!vec) { + if (!vec) + { INFO( "BoostXmlGmlInterface::addPolylinesToPropertyTree(): " "No polylines within the geometry '{:s}'.", @@ -371,7 +394,7 @@ void BoostXmlGmlInterface::addPolylinesToPropertyTree( return; } - std::vector<GeoLib::Polyline*> const*const polylines(vec->getVector()); + std::vector<GeoLib::Polyline*> const* const polylines(vec->getVector()); if (!polylines || polylines->empty()) { INFO( @@ -382,8 +405,9 @@ void BoostXmlGmlInterface::addPolylinesToPropertyTree( } auto& polylines_tag = geometry_set.add("polylines", ""); - for (std::size_t i=0; i<polylines->size(); ++i) { - GeoLib::Polyline const*const polyline((*polylines)[i]); + for (std::size_t i = 0; i < polylines->size(); ++i) + { + GeoLib::Polyline const* const polyline((*polylines)[i]); std::string ply_name; vec->getNameOfElement(polyline, ply_name); auto& polyline_tag = polylines_tag.add("polyline", ""); @@ -392,11 +416,12 @@ void BoostXmlGmlInterface::addPolylinesToPropertyTree( { polyline_tag.put("<xmlattr>.name", ply_name); } - for (std::size_t j=0; j<polyline->getNumberOfPoints(); ++j) { + for (std::size_t j = 0; j < polyline->getNumberOfPoints(); ++j) + { polyline_tag.add("pnt", polyline->getPointID(j)); } } } -} // end namespace IO -} // end namespace GeoLib +} // end namespace IO +} // end namespace GeoLib diff --git a/GeoLib/IO/XmlIO/Qt/XmlGmlInterface.cpp b/GeoLib/IO/XmlIO/Qt/XmlGmlInterface.cpp index b0274052f6dc94895e05e44c8a01f6dc45ec94f5..5464f7e9209f5857c71e2ba1d4f049cb536cced7 100644 --- a/GeoLib/IO/XmlIO/Qt/XmlGmlInterface.cpp +++ b/GeoLib/IO/XmlIO/Qt/XmlGmlInterface.cpp @@ -19,7 +19,6 @@ #include <QtXml/QDomDocument> #include "BaseLib/Logging.h" - #include "GeoLib/Triangle.h" namespace @@ -63,7 +62,7 @@ XmlGmlInterface::XmlGmlInterface(GeoLib::GEOObjects& geo_objs) { } -int XmlGmlInterface::readFile(const QString &fileName) +int XmlGmlInterface::readFile(const QString& fileName) { if (XMLQtInterface::readFile(fileName) == 0) { @@ -72,7 +71,7 @@ int XmlGmlInterface::readFile(const QString &fileName) QDomDocument doc("OGS-GLI-DOM"); doc.setContent(getContent()); - QDomElement docElement = doc.documentElement(); //OpenGeoSysGLI + QDomElement docElement = doc.documentElement(); // OpenGeoSysGLI if (docElement.nodeName().compare("OpenGeoSysGLI")) { ERR("XmlGmlInterface::readFile() - Unexpected XML root."); @@ -105,7 +104,7 @@ int XmlGmlInterface::readFile(const QString &fileName) return 0; } - gliName = type_node.toElement().text().toStdString(); + gliName = type_node.toElement().text().toStdString(); } else if (nodeName.compare("points") == 0) { @@ -187,8 +186,9 @@ void XmlGmlInterface::readPoints(const QDomNode& pointsRoot, QDomElement point = pointsRoot.firstChildElement(); while (!point.isNull()) { - _idx_map.insert (std::pair<std::size_t, std::size_t>( - strtol((point.attribute("id")).toStdString().c_str(), &pEnd, 10), points->size())); + _idx_map.insert(std::pair<std::size_t, std::size_t>( + strtol((point.attribute("id")).toStdString().c_str(), &pEnd, 10), + points->size())); GeoLib::Point* p = new GeoLib::Point(point.attribute("x").toDouble(), point.attribute("y").toDouble(), point.attribute("z").toDouble(), @@ -217,16 +217,19 @@ void XmlGmlInterface::readPolylines( std::size_t const idx = polylines->size(); polylines->push_back(new GeoLib::Polyline(points)); - if (polyline.hasAttribute("name")) { + if (polyline.hasAttribute("name")) + { std::string const ply_name( - polyline.attribute("name").toStdString() - ); + polyline.attribute("name").toStdString()); std::map<std::string, std::size_t>::const_iterator it( - ply_names->find(ply_name) - ); - if (it == ply_names->end()) { - ply_names->insert(std::pair<std::string, std::size_t>(ply_name, idx)); - } else { + ply_names->find(ply_name)); + if (it == ply_names->end()) + { + ply_names->insert( + std::pair<std::string, std::size_t>(ply_name, idx)); + } + else + { WARN( "Polyline '{:s}' exists already. Polyline {:d} will be " "inserted without a name.", @@ -235,24 +238,22 @@ void XmlGmlInterface::readPolylines( } QDomElement point = polyline.firstChildElement(); - auto accessOrError = - [this, &polyline](auto pt_idx) { - auto search = _idx_map.find(pt_idx); - if (search == _idx_map.end()) + auto accessOrError = [this, &polyline](auto pt_idx) { + auto search = _idx_map.find(pt_idx); + if (search == _idx_map.end()) + { + std::string polyline_name; + if (polyline.hasAttribute("name")) { - std::string polyline_name; - if (polyline.hasAttribute("name")) - { - polyline_name = - polyline.attribute("name").toStdString(); - } - OGS_FATAL( - "Polyline `{:s}' contains the point id `{:d}' which is " - "not in the point list.", - polyline_name, pt_idx); + polyline_name = polyline.attribute("name").toStdString(); } - return search->second; - }; + OGS_FATAL( + "Polyline `{:s}' contains the point id `{:d}' which is " + "not in the point list.", + polyline_name, pt_idx); + } + return search->second; + }; while (!point.isNull()) { @@ -283,23 +284,22 @@ void XmlGmlInterface::readSurfaces( surface.attribute("name").toStdString(), surfaces->size() - 1)); } - auto accessOrError = - [this, &surface](auto pt_idx) { - auto search = _idx_map.find(pt_idx); - if (search == _idx_map.end()) + auto accessOrError = [this, &surface](auto pt_idx) { + auto search = _idx_map.find(pt_idx); + if (search == _idx_map.end()) + { + std::string surface_name; + if (surface.hasAttribute("name")) { - std::string surface_name; - if (surface.hasAttribute("name")) - { - surface_name = surface.attribute("name").toStdString(); - } - OGS_FATAL( - "Surface `{:s}' contains the point id `{:d}', which is " - "not in the point list.", - surface_name, pt_idx); + surface_name = surface.attribute("name").toStdString(); } - return search->second; - }; + OGS_FATAL( + "Surface `{:s}' contains the point id `{:d}', which is " + "not in the point list.", + surface_name, pt_idx); + } + return search->second; + }; QDomElement element = surface.firstChildElement(); while (!element.isNull()) @@ -310,7 +310,7 @@ void XmlGmlInterface::readSurfaces( pnt_id_map[accessOrError(element.attribute("p2").toInt())]; std::size_t p3 = pnt_id_map[accessOrError(element.attribute("p3").toInt())]; - surfaces->back()->addTriangle(p1,p2,p3); + surfaces->back()->addTriangle(p1, p2, p3); element = element.nextSiblingElement(); } @@ -330,8 +330,8 @@ bool XmlGmlInterface::write() QDomDocument doc("OGS-GML-DOM"); QDomElement root = doc.createElement("OpenGeoSysGLI"); - root.setAttribute( "xmlns:ogs", "http://www.opengeosys.org" ); - root.setAttribute( "xmlns:xsi", "http://www.w3.org/2001/XMLSchema-instance" ); + root.setAttribute("xmlns:ogs", "http://www.opengeosys.org"); + root.setAttribute("xmlns:xsi", "http://www.w3.org/2001/XMLSchema-instance"); doc.appendChild(root); @@ -348,7 +348,7 @@ bool XmlGmlInterface::write() const GeoLib::PointVec* pnt_vec(_geo_objs.getPointVecObj(export_name)); if (pnt_vec) { - const std::vector<GeoLib::Point*>* points (pnt_vec->getVector()); + const std::vector<GeoLib::Point*>* points(pnt_vec->getVector()); if (!points->empty()) { @@ -382,13 +382,15 @@ bool XmlGmlInterface::write() } else { - ERR("XmlGmlInterface::write(): Point vector is empty, abort writing geometry."); + ERR("XmlGmlInterface::write(): Point vector is empty, abort " + "writing geometry."); return false; } } else { - ERR("XmlGmlInterface::write(): No point vector found, abort writing geometry."); + ERR("XmlGmlInterface::write(): No point vector found, abort writing " + "geometry."); return false; } @@ -397,7 +399,7 @@ bool XmlGmlInterface::write() _geo_objs.getPolylineVecObj(export_name)); if (ply_vec) { - const std::vector<GeoLib::Polyline*>* polylines (ply_vec->getVector()); + const std::vector<GeoLib::Polyline*>* polylines(ply_vec->getVector()); if (polylines) { @@ -412,11 +414,16 @@ bool XmlGmlInterface::write() polylineTag.setAttribute("id", QString::number(i)); std::string ply_name; - if (ply_vec->getNameOfElementByID(i, ply_name)) { - polylineTag.setAttribute("name", QString::fromStdString(ply_name)); - } else { + if (ply_vec->getNameOfElementByID(i, ply_name)) + { + polylineTag.setAttribute( + "name", QString::fromStdString(ply_name)); + } + else + { ply_name = std::to_string(i); - polylineTag.setAttribute("name", QString::fromStdString(ply_name)); + polylineTag.setAttribute( + "name", QString::fromStdString(ply_name)); } plyListTag.appendChild(polylineTag); @@ -426,27 +433,32 @@ bool XmlGmlInterface::write() { QDomElement plyPointTag = doc.createElement("pnt"); polylineTag.appendChild(plyPointTag); - QDomText plyPointText = doc.createTextNode(QString::number(((*polylines)[i])->getPointID(j))); + QDomText plyPointText = doc.createTextNode( + QString::number(((*polylines)[i])->getPointID(j))); plyPointTag.appendChild(plyPointText); } } } else { - INFO("XmlGmlInterface::write(): Polyline vector is empty, no polylines written to file."); + INFO( + "XmlGmlInterface::write(): Polyline vector is empty, no " + "polylines written to file."); } } } else { - INFO("XmlGmlInterface::write(): Polyline vector is empty, no polylines written to file."); + INFO( + "XmlGmlInterface::write(): Polyline vector is empty, no polylines " + "written to file."); } // SURFACES const GeoLib::SurfaceVec* sfc_vec(_geo_objs.getSurfaceVecObj(export_name)); if (sfc_vec) { - const std::vector<GeoLib::Surface*>* surfaces (sfc_vec->getVector()); + const std::vector<GeoLib::Surface*>* surfaces(sfc_vec->getVector()); if (surfaces) { @@ -470,26 +482,34 @@ bool XmlGmlInterface::write() sfcListTag.appendChild(surfaceTag); // writing the elements compromising the surface - std::size_t nElements = ((*surfaces)[i])->getNumberOfTriangles(); + std::size_t nElements = + ((*surfaces)[i])->getNumberOfTriangles(); for (std::size_t j = 0; j < nElements; j++) { QDomElement elementTag = doc.createElement("element"); - elementTag.setAttribute("p1", QString::number((*(*(*surfaces)[i])[j])[0])); - elementTag.setAttribute("p2", QString::number((*(*(*surfaces)[i])[j])[1])); - elementTag.setAttribute("p3", QString::number((*(*(*surfaces)[i])[j])[2])); + elementTag.setAttribute( + "p1", QString::number((*(*(*surfaces)[i])[j])[0])); + elementTag.setAttribute( + "p2", QString::number((*(*(*surfaces)[i])[j])[1])); + elementTag.setAttribute( + "p3", QString::number((*(*(*surfaces)[i])[j])[2])); surfaceTag.appendChild(elementTag); } } } else { - INFO("XmlGmlInterface::write(): Surface vector is empty, no surfaces written to file."); + INFO( + "XmlGmlInterface::write(): Surface vector is empty, no " + "surfaces written to file."); } } } else { - INFO("XmlGmlInterface::write(): Surface vector is empty, no surfaces written to file."); + INFO( + "XmlGmlInterface::write(): Surface vector is empty, no surfaces " + "written to file."); } std::string xml = doc.toString().toStdString(); diff --git a/GeoLib/IO/XmlIO/Qt/XmlStnInterface.cpp b/GeoLib/IO/XmlIO/Qt/XmlStnInterface.cpp index 59e15aa4e112582d3e2d29478f9803149f9680b5..5321bb8351cc9cb0cc60c3c9d85ae54467eaced2 100644 --- a/GeoLib/IO/XmlIO/Qt/XmlStnInterface.cpp +++ b/GeoLib/IO/XmlIO/Qt/XmlStnInterface.cpp @@ -14,18 +14,16 @@ #include "XmlStnInterface.h" -#include <limits> - #include <QFile> #include <QtXml/QDomDocument> +#include <limits> #include <rapidxml.hpp> -#include "BaseLib/Logging.h" #include "BaseLib/DateTools.h" #include "BaseLib/FileTools.h" - -#include "GeoLib/StationBorehole.h" +#include "BaseLib/Logging.h" #include "GeoLib/GEOObjects.h" +#include "GeoLib/StationBorehole.h" namespace GeoLib { @@ -36,7 +34,7 @@ XmlStnInterface::XmlStnInterface(GeoLib::GEOObjects& geo_objs) { } -int XmlStnInterface::readFile(const QString &fileName) +int XmlStnInterface::readFile(const QString& fileName) { if (XMLQtInterface::readFile(fileName) == 0) { @@ -45,7 +43,8 @@ int XmlStnInterface::readFile(const QString &fileName) QDomDocument doc("OGS-STN-DOM"); doc.setContent(getContent()); - QDomElement docElement = doc.documentElement(); //root element, used for identifying file-type + QDomElement docElement = + doc.documentElement(); // root element, used for identifying file-type if (docElement.nodeName().compare("OpenGeoSysSTN")) { ERR("XmlStnInterface::readFile(): Unexpected XML root."); @@ -70,7 +69,8 @@ int XmlStnInterface::readFile(const QString &fileName) } else if (station_type.compare("stations") == 0) { - readStations(station_node, stations.get(), fileName.toStdString()); + readStations(station_node, stations.get(), + fileName.toStdString()); } else if (station_type.compare("boreholes") == 0) { @@ -88,9 +88,9 @@ int XmlStnInterface::readFile(const QString &fileName) return 1; } -void XmlStnInterface::readStations( const QDomNode &stationsRoot, - std::vector<GeoLib::Point*>* stations, - const std::string &station_file_name ) +void XmlStnInterface::readStations(const QDomNode& stationsRoot, + std::vector<GeoLib::Point*>* stations, + const std::string& station_file_name) { QDomElement station = stationsRoot.firstChildElement(); while (!station.isNull()) @@ -105,12 +105,12 @@ void XmlStnInterface::readStations( const QDomNode &stationsRoot, double stationValue(0.0); QDomNodeList stationFeatures = station.childNodes(); - for(int i = 0; i < stationFeatures.count(); i++) + for (int i = 0; i < stationFeatures.count(); i++) { // check for general station features - const QDomNode feature_node (stationFeatures.at(i)); - const QString feature_name (feature_node.nodeName()); - const QString element_text (feature_node.toElement().text()); + const QDomNode feature_node(stationFeatures.at(i)); + const QString feature_name(feature_node.nodeName()); + const QString element_text(feature_node.toElement().text()); if (feature_name.compare("name") == 0) { stationName = element_text.toStdString(); @@ -118,9 +118,9 @@ void XmlStnInterface::readStations( const QDomNode &stationsRoot, if (feature_name.compare("sensordata") == 0) { sensor_data_file_name = element_text.toStdString(); - /* add other station features here */ + /* add other station features here */ - // check for general borehole features + // check for general borehole features } else if (feature_name.compare("value") == 0) { @@ -137,14 +137,17 @@ void XmlStnInterface::readStations( const QDomNode &stationsRoot, /* add other borehole features here */ } - double zVal = (station.hasAttribute("z")) ? station.attribute("z").toDouble() : 0.0; + double zVal = (station.hasAttribute("z")) + ? station.attribute("z").toDouble() + : 0.0; if (station.nodeName().compare("station") == 0) { - GeoLib::Station* s = new GeoLib::Station(station.attribute("x").toDouble(), - station.attribute("y").toDouble(), - zVal, - stationName); + GeoLib::Station* s = + new GeoLib::Station(station.attribute("x").toDouble(), + station.attribute("y").toDouble(), + zVal, + stationName); s->setStationValue(stationValue); if (!sensor_data_file_name.empty()) { @@ -155,7 +158,8 @@ void XmlStnInterface::readStations( const QDomNode &stationsRoot, } else if (station.nodeName().compare("borehole") == 0) { - GeoLib::StationBorehole* s = GeoLib::StationBorehole::createStation( + GeoLib::StationBorehole* s = + GeoLib::StationBorehole::createStation( stationName, station.attribute("x").toDouble(), station.attribute("y").toDouble(), @@ -177,16 +181,19 @@ void XmlStnInterface::readStations( const QDomNode &stationsRoot, } else { - WARN("XmlStnInterface::readStations(): Attribute missing in <station> tag."); + WARN( + "XmlStnInterface::readStations(): Attribute missing in " + "<station> tag."); } station = station.nextSiblingElement(); } } -void XmlStnInterface::readStratigraphy( const QDomNode &stratRoot, - GeoLib::StationBorehole* borehole ) +void XmlStnInterface::readStratigraphy(const QDomNode& stratRoot, + GeoLib::StationBorehole* borehole) { - //borehole->addSoilLayer((*borehole)[0], (*borehole)[1], (*borehole)[2], ""); + // borehole->addSoilLayer((*borehole)[0], (*borehole)[1], (*borehole)[2], + // ""); double depth_check((*borehole)[2]); QDomElement horizon = stratRoot.firstChildElement(); while (!horizon.isNull()) @@ -207,7 +214,7 @@ void XmlStnInterface::readStratigraphy( const QDomNode &stratRoot, } /* add other horizon features here */ - double depth (horizon.attribute("z").toDouble()); + double depth(horizon.attribute("z").toDouble()); if (std::abs(depth - depth_check) > std::numeric_limits<double>:: epsilon()) // skip soil-layer if its thickness is zero @@ -228,7 +235,9 @@ void XmlStnInterface::readStratigraphy( const QDomNode &stratRoot, } else { - WARN("XmlStnInterface::readStratigraphy(): Attribute missing in <horizon> tag."); + WARN( + "XmlStnInterface::readStratigraphy(): Attribute missing in " + "<horizon> tag."); } horizon = horizon.nextSiblingElement(); } @@ -249,13 +258,14 @@ bool XmlStnInterface::write() QDomDocument doc("OGS-STN-DOM"); QDomElement root = doc.createElement("OpenGeoSysSTN"); - root.setAttribute( "xmlns:ogs", "http://www.opengeosys.org" ); - root.setAttribute( "xmlns:xsi", "http://www.w3.org/2001/XMLSchema-instance" ); + root.setAttribute("xmlns:ogs", "http://www.opengeosys.org"); + root.setAttribute("xmlns:xsi", "http://www.w3.org/2001/XMLSchema-instance"); const std::vector<GeoLib::Point*>* stations( _geo_objs.getStationVec(export_name)); - bool const is_borehole = static_cast<GeoLib::Station*>((*stations)[0])->type() == - GeoLib::Station::StationType::BOREHOLE; + bool const is_borehole = + static_cast<GeoLib::Station*>((*stations)[0])->type() == + GeoLib::Station::StationType::BOREHOLE; doc.appendChild(root); QDomElement stationListTag = doc.createElement("stationlist"); @@ -271,12 +281,13 @@ bool XmlStnInterface::write() stationListTag.appendChild(stationsTag); bool useStationValue(false); - double sValue = static_cast<GeoLib::Station*>((*stations)[0])->getStationValue(); + double sValue = + static_cast<GeoLib::Station*>((*stations)[0])->getStationValue(); std::size_t nStations(stations->size()); for (std::size_t i = 1; i < nStations; i++) { - if ((static_cast<GeoLib::Station*>((*stations)[i])->getStationValue() - sValue) < - std::numeric_limits<double>::epsilon()) + if ((static_cast<GeoLib::Station*>((*stations)[i])->getStationValue() - + sValue) < std::numeric_limits<double>::epsilon()) { useStationValue = true; break; @@ -285,29 +296,33 @@ bool XmlStnInterface::write() for (std::size_t i = 0; i < nStations; i++) { - QString stationType = is_borehole ? "borehole" : "station"; + QString stationType = is_borehole ? "borehole" : "station"; QDomElement stationTag = doc.createElement(stationType); - stationTag.setAttribute( "id", QString::number(i) ); - stationTag.setAttribute( "x", QString::number((*(*stations)[i])[0], 'f', - std::numeric_limits<double>::digits10)); - stationTag.setAttribute( "y", QString::number((*(*stations)[i])[1], 'f', - std::numeric_limits<double>::digits10)); - stationTag.setAttribute( "z", QString::number((*(*stations)[i])[2], 'f', - std::numeric_limits<double>::digits10)); + stationTag.setAttribute("id", QString::number(i)); + stationTag.setAttribute( + "x", QString::number((*(*stations)[i])[0], 'f', + std::numeric_limits<double>::digits10)); + stationTag.setAttribute( + "y", QString::number((*(*stations)[i])[1], 'f', + std::numeric_limits<double>::digits10)); + stationTag.setAttribute( + "z", QString::number((*(*stations)[i])[2], 'f', + std::numeric_limits<double>::digits10)); stationsTag.appendChild(stationTag); QDomElement stationNameTag = doc.createElement("name"); stationTag.appendChild(stationNameTag); - QDomText stationNameText = - doc.createTextNode(QString::fromStdString(static_cast<GeoLib::Station*>((*stations)[i])->getName())); + QDomText stationNameText = doc.createTextNode(QString::fromStdString( + static_cast<GeoLib::Station*>((*stations)[i])->getName())); stationNameTag.appendChild(stationNameText); if (useStationValue) { QDomElement stationValueTag = doc.createElement("value"); stationTag.appendChild(stationValueTag); - QDomText stationValueText = - doc.createTextNode(QString::number(static_cast<GeoLib::Station*>((*stations)[i])->getStationValue())); + QDomText stationValueText = doc.createTextNode( + QString::number(static_cast<GeoLib::Station*>((*stations)[i]) + ->getStationValue())); stationValueTag.appendChild(stationValueText); } @@ -324,21 +339,21 @@ bool XmlStnInterface::write() return true; } -void XmlStnInterface::writeBoreholeData(QDomDocument &doc, - QDomElement &boreholeTag, +void XmlStnInterface::writeBoreholeData(QDomDocument& doc, + QDomElement& boreholeTag, GeoLib::StationBorehole* borehole) const { QDomElement stationDepthTag = doc.createElement("bdepth"); boreholeTag.appendChild(stationDepthTag); - QDomText stationDepthText = doc.createTextNode(QString::number(borehole->getDepth(), 'f')); + QDomText stationDepthText = + doc.createTextNode(QString::number(borehole->getDepth(), 'f')); stationDepthTag.appendChild(stationDepthText); if (std::abs(borehole->getDate()) > 0) { QDomElement stationDateTag = doc.createElement("bdate"); boreholeTag.appendChild(stationDateTag); - QDomText stationDateText = - doc.createTextNode(QString::fromStdString(BaseLib::date2string(borehole-> - getDate()))); + QDomText stationDateText = doc.createTextNode( + QString::fromStdString(BaseLib::date2string(borehole->getDate()))); stationDateTag.appendChild(stationDateText); } @@ -351,18 +366,23 @@ void XmlStnInterface::writeBoreholeData(QDomDocument &doc, QDomElement stratTag = doc.createElement("strat"); boreholeTag.appendChild(stratTag); - for (std::size_t j = 1; j < nHorizons; j++) /// the first entry in the profile vector is just the position of the borehole + for (std::size_t j = 1; j < nHorizons; + j++) /// the first entry in the profile vector is just the + /// position of the borehole { QDomElement horizonTag = doc.createElement("horizon"); - horizonTag.setAttribute( "id", QString::number(j) ); - horizonTag.setAttribute( "x", QString::number((*profile[j])[0], 'f') ); - horizonTag.setAttribute( "y", QString::number((*profile[j])[1], 'f') ); - horizonTag.setAttribute( "z", QString::number((*profile[j])[2], 'f') ); + horizonTag.setAttribute("id", QString::number(j)); + horizonTag.setAttribute("x", + QString::number((*profile[j])[0], 'f')); + horizonTag.setAttribute("y", + QString::number((*profile[j])[1], 'f')); + horizonTag.setAttribute("z", + QString::number((*profile[j])[2], 'f')); stratTag.appendChild(horizonTag); QDomElement horizonNameTag = doc.createElement("name"); horizonTag.appendChild(horizonNameTag); QDomText horizonNameText = - doc.createTextNode(QString::fromStdString(soilNames[j])); + doc.createTextNode(QString::fromStdString(soilNames[j])); horizonNameTag.appendChild(horizonNameText); } } diff --git a/GeoLib/LineSegment.cpp b/GeoLib/LineSegment.cpp index 20962f9de6451091dc72a393317e0ffdb37b32e5..0978b23fd16fabf8408341504212caac24bf0d2e 100644 --- a/GeoLib/LineSegment.cpp +++ b/GeoLib/LineSegment.cpp @@ -17,13 +17,15 @@ LineSegment::LineSegment(Point* const a, Point* const b, _b(b), _point_mem_management_by_line_segment( point_mem_management_by_line_segment) -{} +{ +} LineSegment::LineSegment(LineSegment const& line_segment) : _a(new Point(line_segment.getBeginPoint())), _b(new Point(line_segment.getEndPoint())), _point_mem_management_by_line_segment(true) -{} +{ +} LineSegment::LineSegment(LineSegment&& line_segment) : _a(line_segment._a), @@ -38,7 +40,8 @@ LineSegment::LineSegment(LineSegment&& line_segment) LineSegment::~LineSegment() { - if (_point_mem_management_by_line_segment) { + if (_point_mem_management_by_line_segment) + { delete _b; delete _a; } @@ -65,7 +68,7 @@ Point const& LineSegment::getBeginPoint() const return *_a; } -Point & LineSegment::getBeginPoint() +Point& LineSegment::getBeginPoint() { return *_a; } @@ -75,12 +78,12 @@ Point const& LineSegment::getEndPoint() const return *_b; } -Point & LineSegment::getEndPoint() +Point& LineSegment::getEndPoint() { return *_b; } -std::ostream& operator<< (std::ostream& os, LineSegment const& s) +std::ostream& operator<<(std::ostream& os, LineSegment const& s) { os << "{(" << s.getBeginPoint() << "), (" << s.getEndPoint() << ")}"; return os; @@ -98,8 +101,8 @@ bool operator==(LineSegment const& s0, LineSegment const& s1) { double const tol(std::numeric_limits<double>::epsilon()); return (MathLib::sqrDist(s0.getBeginPoint(), s1.getBeginPoint()) < tol && - MathLib::sqrDist(s0.getEndPoint(), s1.getEndPoint()) < tol) || - (MathLib::sqrDist(s0.getBeginPoint(), s1.getEndPoint()) < tol && - MathLib::sqrDist(s0.getEndPoint(), s1.getBeginPoint()) < tol); + MathLib::sqrDist(s0.getEndPoint(), s1.getEndPoint()) < tol) || + (MathLib::sqrDist(s0.getBeginPoint(), s1.getEndPoint()) < tol && + MathLib::sqrDist(s0.getEndPoint(), s1.getBeginPoint()) < tol); } } // namespace GeoLib diff --git a/GeoLib/MinimalBoundingSphere.cpp b/GeoLib/MinimalBoundingSphere.cpp index 3e989cfab3a66843ce766b6864828b9022e3db3e..f1424620adada5e24ad1a14c049de864e8a1ff06 100644 --- a/GeoLib/MinimalBoundingSphere.cpp +++ b/GeoLib/MinimalBoundingSphere.cpp @@ -16,22 +16,23 @@ #include <ctime> -#include "MathLib/Point3d.h" #include "MathLib/GeometricBasics.h" #include "MathLib/MathTools.h" +#include "MathLib/Point3d.h" -namespace GeoLib { +namespace GeoLib +{ MinimalBoundingSphere::MinimalBoundingSphere() = default; -MinimalBoundingSphere::MinimalBoundingSphere( - MathLib::Point3d const& p, double radius) -: _radius(radius), _center(p) +MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p, + double radius) + : _radius(radius), _center(p) { } -MinimalBoundingSphere::MinimalBoundingSphere( - MathLib::Point3d const& p, MathLib::Point3d const& q) -: _radius(std::numeric_limits<double>::epsilon()), _center(p) +MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p, + MathLib::Point3d const& q) + : _radius(std::numeric_limits<double>::epsilon()), _center(p) { auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords()); auto const vq = Eigen::Map<Eigen::Vector3d const>(q.getCoords()); @@ -44,7 +45,8 @@ MinimalBoundingSphere::MinimalBoundingSphere( } MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p, - MathLib::Point3d const& q, MathLib::Point3d const& r) + MathLib::Point3d const& q, + MathLib::Point3d const& r) { auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords()); auto const vq = Eigen::Map<Eigen::Vector3d const>(q.getCoords()); @@ -67,7 +69,7 @@ MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p, MinimalBoundingSphere two_pnts_sphere; if (a.squaredNorm() > b.squaredNorm()) { - two_pnts_sphere = MinimalBoundingSphere(p,r); + two_pnts_sphere = MinimalBoundingSphere(p, r); } else { @@ -79,9 +81,9 @@ MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p, } MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p, - MathLib::Point3d const& q, - MathLib::Point3d const& r, - MathLib::Point3d const& s) + MathLib::Point3d const& q, + MathLib::Point3d const& r, + MathLib::Point3d const& s) { auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords()); auto const vq = Eigen::Map<Eigen::Vector3d const>(q.getCoords()); @@ -106,10 +108,10 @@ MinimalBoundingSphere::MinimalBoundingSphere(MathLib::Point3d const& p, } else { - MinimalBoundingSphere const pqr(p, q , r); - MinimalBoundingSphere const pqs(p, q , s); - MinimalBoundingSphere const prs(p, r , s); - MinimalBoundingSphere const qrs(q, r , s); + MinimalBoundingSphere const pqr(p, q, r); + MinimalBoundingSphere const pqs(p, q, s); + MinimalBoundingSphere const prs(p, r, s); + MinimalBoundingSphere const qrs(q, r, s); _radius = pqr.getRadius(); _center = pqr.getCenter(); if (_radius < pqs.getRadius()) @@ -135,62 +137,72 @@ MinimalBoundingSphere::MinimalBoundingSphere( : _radius(-1), _center({0, 0, 0}) { const std::vector<MathLib::Point3d*>& sphere_points(points); - MinimalBoundingSphere const bounding_sphere = recurseCalculation(sphere_points, 0, sphere_points.size(), 0); + MinimalBoundingSphere const bounding_sphere = + recurseCalculation(sphere_points, 0, sphere_points.size(), 0); _center = bounding_sphere.getCenter(); _radius = bounding_sphere.getRadius(); } -MinimalBoundingSphere -MinimalBoundingSphere::recurseCalculation( +MinimalBoundingSphere MinimalBoundingSphere::recurseCalculation( std::vector<MathLib::Point3d*> sphere_points, std::size_t start_idx, std::size_t length, std::size_t n_boundary_points) { MinimalBoundingSphere sphere; - switch(n_boundary_points) + switch (n_boundary_points) { - case 0: - sphere = MinimalBoundingSphere(); - break; - case 1: - sphere = MinimalBoundingSphere(*sphere_points[start_idx-1]); - break; - case 2: - sphere = MinimalBoundingSphere(*sphere_points[start_idx-1], *sphere_points[start_idx-2]); - break; - case 3: - sphere = MinimalBoundingSphere(*sphere_points[start_idx-1], *sphere_points[start_idx-2], *sphere_points[start_idx-3]); - break; - case 4: - sphere = MinimalBoundingSphere(*sphere_points[start_idx-1], *sphere_points[start_idx-2], *sphere_points[start_idx-3], *sphere_points[start_idx-4]); - return sphere; + case 0: + sphere = MinimalBoundingSphere(); + break; + case 1: + sphere = MinimalBoundingSphere(*sphere_points[start_idx - 1]); + break; + case 2: + sphere = MinimalBoundingSphere(*sphere_points[start_idx - 1], + *sphere_points[start_idx - 2]); + break; + case 3: + sphere = MinimalBoundingSphere(*sphere_points[start_idx - 1], + *sphere_points[start_idx - 2], + *sphere_points[start_idx - 3]); + break; + case 4: + sphere = MinimalBoundingSphere( + *sphere_points[start_idx - 1], *sphere_points[start_idx - 2], + *sphere_points[start_idx - 3], *sphere_points[start_idx - 4]); + return sphere; } - for(std::size_t i=0; i<length; ++i) + for (std::size_t i = 0; i < length; ++i) { // current point is located outside of sphere - if (sphere.pointDistanceSquared(*sphere_points[start_idx+i]) > 0) + if (sphere.pointDistanceSquared(*sphere_points[start_idx + i]) > 0) { - if (i>start_idx) + if (i > start_idx) { - using DiffType = std::vector<MathLib::Point3d*>::iterator::difference_type; + using DiffType = + std::vector<MathLib::Point3d*>::iterator::difference_type; std::vector<MathLib::Point3d*> const tmp_ps( sphere_points.cbegin() + static_cast<DiffType>(start_idx), - sphere_points.cbegin() + static_cast<DiffType>(start_idx + i + 1)); + sphere_points.cbegin() + + static_cast<DiffType>(start_idx + i + 1)); std::copy(tmp_ps.cbegin(), --tmp_ps.cend(), - sphere_points.begin() + static_cast<DiffType>(start_idx + 1)); + sphere_points.begin() + + static_cast<DiffType>(start_idx + 1)); sphere_points[start_idx] = tmp_ps.back(); } - sphere = recurseCalculation(sphere_points, start_idx+1, i, n_boundary_points+1); + sphere = recurseCalculation(sphere_points, start_idx + 1, i, + n_boundary_points + 1); } } return sphere; } -double MinimalBoundingSphere::pointDistanceSquared(MathLib::Point3d const& pnt) const +double MinimalBoundingSphere::pointDistanceSquared( + MathLib::Point3d const& pnt) const { - return MathLib::sqrDist(_center, pnt)-(_radius*_radius); + return MathLib::sqrDist(_center, pnt) - (_radius * _radius); } } // namespace GeoLib diff --git a/GeoLib/PointVec.cpp b/GeoLib/PointVec.cpp index 5889627f03075d33e3d09b8814cfc4990cc583ec..6822171f434cd1f1a421e3b177deed1a24649ff9 100644 --- a/GeoLib/PointVec.cpp +++ b/GeoLib/PointVec.cpp @@ -12,12 +12,11 @@ * */ +#include "PointVec.h" + #include <numeric> #include "BaseLib/Logging.h" - -#include "PointVec.h" - #include "MathLib/MathTools.h" namespace GeoLib diff --git a/GeoLib/Polygon.cpp b/GeoLib/Polygon.cpp index 5c55877ec65acaf42e070521a15c6ed3279096d9..07dba6090bb16f5ba1a0e9e83e15ffeac8f53559 100644 --- a/GeoLib/Polygon.cpp +++ b/GeoLib/Polygon.cpp @@ -16,14 +16,13 @@ #include <boost/math/constants/constants.hpp> -#include "BaseLib/quicksort.h" - #include "AnalyticalGeometry.h" +#include "BaseLib/quicksort.h" namespace GeoLib { -Polygon::Polygon(const Polyline &ply, bool init) : - Polyline(ply), _aabb(ply.getPointsVec(), ply._ply_pnt_ids) +Polygon::Polygon(const Polyline& ply, bool init) + : Polyline(ply), _aabb(ply.getPointsVec(), ply._ply_pnt_ids) { if (init) { @@ -32,8 +31,7 @@ Polygon::Polygon(const Polyline &ply, bool init) : _simple_polygon_list.push_back(this); } -Polygon::Polygon(Polygon const& other) - : Polyline(other), _aabb(other._aabb) +Polygon::Polygon(Polygon const& other) : Polyline(other), _aabb(other._aabb) { _simple_polygon_list.push_back(this); auto sub_polygon_it(other._simple_polygon_list.begin()); @@ -59,9 +57,10 @@ Polygon::~Polygon() } } -bool Polygon::initialise () +bool Polygon::initialise() { - if (this->isClosed()) { + if (this->isClosed()) + { ensureCCWOrientation(); return true; } @@ -128,18 +127,19 @@ bool Polygon::isPntInPolygon(GeoLib::Point const& pnt) const bool Polygon::isPntInPolygon(double x, double y, double z) const { - const GeoLib::Point pnt(x,y,z); + const GeoLib::Point pnt(x, y, z); return isPntInPolygon(pnt); } std::vector<GeoLib::Point> Polygon::getAllIntersectionPoints( - GeoLib::LineSegment const& segment) const + GeoLib::LineSegment const& segment) const { std::vector<GeoLib::Point> intersections; GeoLib::Point s; for (auto&& seg_it : *this) { - if (GeoLib::lineSegmentIntersect(seg_it, segment, s)) { + if (GeoLib::lineSegmentIntersect(seg_it, segment, s)) + { intersections.push_back(s); } } @@ -154,21 +154,25 @@ bool Polygon::containsSegment(GeoLib::LineSegment const& segment) const GeoLib::Point const& a{segment.getBeginPoint()}; GeoLib::Point const& b{segment.getEndPoint()}; // no intersections -> check if at least one point of segment is in polygon - if (s.empty()) { + if (s.empty()) + { return (isPntInPolygon(a)); } const double tol(std::numeric_limits<float>::epsilon()); // one intersection, intersection in line segment end point - if (s.size() == 1) { - const double sqr_dist_as(MathLib::sqrDist(a,s[0])); - if (sqr_dist_as < tol) { + if (s.size() == 1) + { + const double sqr_dist_as(MathLib::sqrDist(a, s[0])); + if (sqr_dist_as < tol) + { return (isPntInPolygon(b)); } - const double sqr_dist_bs(MathLib::sqrDist(b,s[0])); - if (sqr_dist_bs < tol) { + const double sqr_dist_bs(MathLib::sqrDist(b, s[0])); + if (sqr_dist_bs < tol) + { return (isPntInPolygon(a)); } } @@ -176,16 +180,19 @@ bool Polygon::containsSegment(GeoLib::LineSegment const& segment) const // Sorting the intersection with respect to the distance to the point a. // This induces a partition of the line segment into sub segments. std::sort(s.begin(), s.end(), - [&a] (GeoLib::Point const& p0, GeoLib::Point const& p1) { - return MathLib::sqrDist(a, p0) < MathLib::sqrDist(a, p1); - } - ); + [&a](GeoLib::Point const& p0, GeoLib::Point const& p1) { + return MathLib::sqrDist(a, p0) < MathLib::sqrDist(a, p1); + }); // remove sub segments with almost zero length - for (std::size_t k(0); k<s.size()-1; ) { - if (MathLib::sqrDist(s[k], s[k+1]) < tol) { - s.erase(s.begin()+k+1); - } else { + for (std::size_t k(0); k < s.size() - 1;) + { + if (MathLib::sqrDist(s[k], s[k + 1]) < tol) + { + s.erase(s.begin() + k + 1); + } + else + { k++; } } @@ -197,8 +204,9 @@ bool Polygon::containsSegment(GeoLib::LineSegment const& segment) const { return false; } - const std::size_t n_sub_segs(s.size()-1); - for (std::size_t k(0); k<n_sub_segs; k++) { + const std::size_t n_sub_segs(s.size() - 1); + for (std::size_t k(0); k < n_sub_segs; k++) + { if (!isPntInPolygon(GeoLib::Point(0.5 * (s[k][0] + s[k + 1][0]), 0.5 * (s[k][1] + s[k + 1][1]), 0.5 * (s[k][2] + s[k + 1][2])))) @@ -220,7 +228,7 @@ bool Polygon::isPolylineInPolygon(const Polyline& ply) const bool Polygon::isPartOfPolylineInPolygon(const Polyline& ply) const { - const std::size_t ply_size (ply.getNumberOfPoints()); + const std::size_t ply_size(ply.getNumberOfPoints()); // check points for (std::size_t k(0); k < ply_size; k++) { @@ -247,15 +255,19 @@ bool Polygon::getNextIntersectionPointPolygonLine( GeoLib::LineSegment const& seg, GeoLib::Point& intersection_pnt, std::size_t& seg_num) const { - if (_simple_polygon_list.size() == 1) { - for (auto seg_it(begin()+seg_num); seg_it != end(); ++seg_it) { + if (_simple_polygon_list.size() == 1) + { + for (auto seg_it(begin() + seg_num); seg_it != end(); ++seg_it) + { if (GeoLib::lineSegmentIntersect(*seg_it, seg, intersection_pnt)) { seg_num = seg_it.getSegmentNumber(); return true; } } - } else { + } + else + { for (auto polygon : _simple_polygon_list) { for (auto seg_it(polygon->begin()); seg_it != polygon->end(); @@ -277,40 +289,42 @@ EdgeType Polygon::getEdgeType(std::size_t k, GeoLib::Point const& pnt) const { switch (getLocationOfPoint(k, pnt)) { - case Location::LEFT: { - const GeoLib::Point & v (*(getPoint(k))); - const GeoLib::Point & w (*(getPoint(k + 1))); - if (v[1] < pnt[1] && pnt[1] <= w[1]) + case Location::LEFT: { - return EdgeType::CROSSING; - } + const GeoLib::Point& v(*(getPoint(k))); + const GeoLib::Point& w(*(getPoint(k + 1))); + if (v[1] < pnt[1] && pnt[1] <= w[1]) + { + return EdgeType::CROSSING; + } - return EdgeType::INESSENTIAL; - } - case Location::RIGHT: { - const GeoLib::Point & v (*(getPoint(k))); - const GeoLib::Point & w (*(getPoint(k + 1))); - if (w[1] < pnt[1] && pnt[1] <= v[1]) - { - return EdgeType::CROSSING; + return EdgeType::INESSENTIAL; } + case Location::RIGHT: + { + const GeoLib::Point& v(*(getPoint(k))); + const GeoLib::Point& w(*(getPoint(k + 1))); + if (w[1] < pnt[1] && pnt[1] <= v[1]) + { + return EdgeType::CROSSING; + } - return EdgeType::INESSENTIAL; - } - case Location::BETWEEN: - case Location::SOURCE: - case Location::DESTINATION: - return EdgeType::TOUCHING; - default: - return EdgeType::INESSENTIAL; + return EdgeType::INESSENTIAL; + } + case Location::BETWEEN: + case Location::SOURCE: + case Location::DESTINATION: + return EdgeType::TOUCHING; + default: + return EdgeType::INESSENTIAL; } } -void Polygon::ensureCCWOrientation () +void Polygon::ensureCCWOrientation() { // *** pre processing: rotate points to xy-plan // *** copy points to vector - last point is identical to the first - std::size_t n_pnts (this->getNumberOfPoints() - 1); + std::size_t n_pnts(this->getNumberOfPoints() - 1); std::vector<GeoLib::Point*> tmp_polygon_pnts; for (std::size_t k(0); k < n_pnts; k++) { @@ -327,10 +341,11 @@ void Polygon::ensureCCWOrientation () } // *** get the left most upper point - std::size_t min_x_max_y_idx (0); // for orientation check + std::size_t min_x_max_y_idx(0); // for orientation check for (std::size_t k(0); k < n_pnts; k++) { - if ((*(tmp_polygon_pnts[k]))[0] <= (*(tmp_polygon_pnts[min_x_max_y_idx]))[0]) + if ((*(tmp_polygon_pnts[k]))[0] <= + (*(tmp_polygon_pnts[min_x_max_y_idx]))[0]) { if ((*(tmp_polygon_pnts[k]))[0] < (*(tmp_polygon_pnts[min_x_max_y_idx]))[0]) @@ -371,8 +386,9 @@ void Polygon::ensureCCWOrientation () if (orient != GeoLib::CCW) { // switch orientation - std::size_t tmp_n_pnts (n_pnts); - tmp_n_pnts++; // include last point of polygon (which is identical to the first) + std::size_t tmp_n_pnts(n_pnts); + tmp_n_pnts++; // include last point of polygon (which is identical to + // the first) for (std::size_t k(0); k < tmp_n_pnts / 2; k++) { std::swap(_ply_pnt_ids[k], _ply_pnt_ids[tmp_n_pnts - 1 - k]); @@ -405,9 +421,9 @@ void Polygon::splitPolygonAtIntersection( std::size_t idx0(seg_it0.getSegmentNumber()); std::size_t idx1(seg_it1.getSegmentNumber()); // adding intersection point to pnt_vec - std::size_t const intersection_pnt_id (_ply_pnts.size()); - const_cast<std::vector<Point*>&>(_ply_pnts) - .push_back(new GeoLib::Point(intersection_pnt)); + std::size_t const intersection_pnt_id(_ply_pnts.size()); + const_cast<std::vector<Point*>&>(_ply_pnts).push_back( + new GeoLib::Point(intersection_pnt)); // split Polygon if (idx0 > idx1) @@ -449,18 +465,19 @@ void Polygon::splitPolygonAtIntersection( splitPolygonAtIntersection(polygon1_it); } -void Polygon::splitPolygonAtPoint (const std::list<GeoLib::Polygon*>::iterator& polygon_it) +void Polygon::splitPolygonAtPoint( + const std::list<GeoLib::Polygon*>::iterator& polygon_it) { std::size_t const n((*polygon_it)->getNumberOfPoints() - 1); std::vector<std::size_t> id_vec(n); std::vector<std::size_t> perm(n); for (std::size_t k(0); k < n; k++) { - id_vec[k] = (*polygon_it)->getPointID (k); + id_vec[k] = (*polygon_it)->getPointID(k); perm[k] = k; } - BaseLib::quicksort (id_vec, 0, n, perm); + BaseLib::quicksort(id_vec, 0, n, perm); for (std::size_t k(0); k < n - 1; k++) { @@ -525,8 +542,10 @@ bool operator==(Polygon const& lhs, Polygon const& rhs) // search start point of first polygon in second polygon bool nfound(true); std::size_t k(0); - for (; k < n-1 && nfound; k++) { - if (start_pnt == rhs.getPointID(k)) { + for (; k < n - 1 && nfound; k++) + { + if (start_pnt == rhs.getPointID(k)) + { nfound = false; break; } @@ -540,26 +559,35 @@ bool operator==(Polygon const& lhs, Polygon const& rhs) // *** determine direction // opposite direction - if (k == n-2) { - for (k=1; k<n-1; k++) { - if (lhs.getPointID(k) != rhs.getPointID(n-1-k)) { + if (k == n - 2) + { + for (k = 1; k < n - 1; k++) + { + if (lhs.getPointID(k) != rhs.getPointID(n - 1 - k)) + { return false; } } return true; } - // same direction - start point of first polygon at arbitrary position in second polygon - if (lhs.getPointID(1) == rhs.getPointID(k+1)) { - std::size_t j(k+2); - for (; j<n-1; j++) { - if (lhs.getPointID(j-k) != rhs.getPointID(j)) { + // same direction - start point of first polygon at arbitrary position in + // second polygon + if (lhs.getPointID(1) == rhs.getPointID(k + 1)) + { + std::size_t j(k + 2); + for (; j < n - 1; j++) + { + if (lhs.getPointID(j - k) != rhs.getPointID(j)) + { return false; } } - j=0; // new start point at second polygon - for (; j<k+1; j++) { - if (lhs.getPointID(n-(k+2)+j+1) != rhs.getPointID(j)) { + j = 0; // new start point at second polygon + for (; j < k + 1; j++) + { + if (lhs.getPointID(n - (k + 2) + j + 1) != rhs.getPointID(j)) + { return false; } } @@ -598,4 +626,4 @@ bool operator==(Polygon const& lhs, Polygon const& rhs) return false; } -} // end namespace GeoLib +} // end namespace GeoLib diff --git a/GeoLib/PolygonWithSegmentMarker.cpp b/GeoLib/PolygonWithSegmentMarker.cpp index fd7d85bb6a0fe0bdd8370af49abe25980a10921f..707cfad4b140563783e78d7e6a1494af991d70ee 100644 --- a/GeoLib/PolygonWithSegmentMarker.cpp +++ b/GeoLib/PolygonWithSegmentMarker.cpp @@ -9,7 +9,8 @@ #include "PolygonWithSegmentMarker.h" -namespace GeoLib { +namespace GeoLib +{ PolygonWithSegmentMarker::PolygonWithSegmentMarker( GeoLib::Polyline const& polyline) : GeoLib::Polygon(polyline, true), @@ -29,7 +30,8 @@ bool PolygonWithSegmentMarker::isSegmentMarked(std::size_t seg_num) const bool PolygonWithSegmentMarker::addPoint(std::size_t pnt_id) { - if (Polyline::addPoint(pnt_id)) { + if (Polyline::addPoint(pnt_id)) + { _marker.push_back(false); return true; } @@ -38,8 +40,9 @@ bool PolygonWithSegmentMarker::addPoint(std::size_t pnt_id) bool PolygonWithSegmentMarker::insertPoint(std::size_t pos, std::size_t pnt_id) { - if (Polyline::insertPoint(pos, pnt_id)) { - _marker.insert(_marker.begin()+pos, _marker[pos]); + if (Polyline::insertPoint(pos, pnt_id)) + { + _marker.insert(_marker.begin() + pos, _marker[pos]); return true; } return false; diff --git a/GeoLib/Polyline.cpp b/GeoLib/Polyline.cpp index 5f9ecd222a98a45ff8e78b70aa063432b7f1a810..c6feec89399eed3fbc1b38cb42507de0cbc62824 100644 --- a/GeoLib/Polyline.cpp +++ b/GeoLib/Polyline.cpp @@ -15,6 +15,7 @@ #include "Polyline.h" #include <algorithm> + #include "AnalyticalGeometry.h" #include "BaseLib/Error.h" #include "BaseLib/Logging.h" @@ -24,16 +25,17 @@ namespace GeoLib { Polyline::Polyline(const std::vector<Point*>& pnt_vec) : _ply_pnts(pnt_vec) { - _length.push_back (0.0); + _length.push_back(0.0); } Polyline::Polyline(const Polyline& ply) : _ply_pnts(ply._ply_pnts), _ply_pnt_ids(ply._ply_pnt_ids), _length(ply._length) -{} +{ +} -void Polyline::write(std::ostream &os) const +void Polyline::write(std::ostream& os) const { std::size_t size(_ply_pnt_ids.size()); for (std::size_t k(0); k < size; k++) @@ -47,7 +49,8 @@ bool Polyline::addPoint(std::size_t pnt_id) assert(pnt_id < _ply_pnts.size()); std::size_t const n_pnts(_ply_pnt_ids.size()); - // don't insert point if this would result in identical IDs for two adjacent points + // don't insert point if this would result in identical IDs for two adjacent + // points if (n_pnts > 0 && _ply_pnt_ids.back() == pnt_id) { return false; @@ -58,7 +61,7 @@ bool Polyline::addPoint(std::size_t pnt_id) if (n_pnts > 0) { double const act_dist(std::sqrt(MathLib::sqrDist( - *_ply_pnts[_ply_pnt_ids[n_pnts-1]], *_ply_pnts[pnt_id]))); + *_ply_pnts[_ply_pnt_ids[n_pnts - 1]], *_ply_pnts[pnt_id]))); double dist_until_now(0.0); if (n_pnts > 1) { @@ -75,12 +78,15 @@ bool Polyline::insertPoint(std::size_t pos, std::size_t pnt_id) assert(pnt_id < _ply_pnts.size()); assert(pos <= _ply_pnt_ids.size()); - if (pos == _ply_pnt_ids.size()) { + if (pos == _ply_pnt_ids.size()) + { return addPoint(pnt_id); } - // check if inserting pnt_id would result in two identical IDs for adjacent points - if (pos == 0 && pnt_id == _ply_pnt_ids[0]) { + // check if inserting pnt_id would result in two identical IDs for adjacent + // points + if (pos == 0 && pnt_id == _ply_pnt_ids[0]) + { return false; } if (pos != 0) @@ -93,16 +99,18 @@ bool Polyline::insertPoint(std::size_t pos, std::size_t pnt_id) { return false; } - } + } auto const pos_dt( static_cast<std::vector<std::size_t>::difference_type>(pos)); auto it(_ply_pnt_ids.begin() + pos_dt); _ply_pnt_ids.insert(it, pnt_id); - if (_ply_pnt_ids.size() > 1) { + if (_ply_pnt_ids.size() > 1) + { // update the _length vector - if (pos == 0) { + if (pos == 0) + { // insert at first position double const act_dist(std::sqrt(MathLib::sqrDist( *_ply_pnts[_ply_pnt_ids[1]], *_ply_pnts[pnt_id]))); @@ -112,13 +120,16 @@ bool Polyline::insertPoint(std::size_t pos, std::size_t pnt_id) { _length[k] += _length[1]; } - } else { - if (pos == _ply_pnt_ids.size() - 1) { + } + else + { + if (pos == _ply_pnt_ids.size() - 1) + { // insert at last position double const act_dist(std::sqrt(MathLib::sqrDist( *_ply_pnts[_ply_pnt_ids[_ply_pnt_ids.size() - 2]], *_ply_pnts[pnt_id]))); - double dist_until_now (0.0); + double dist_until_now(0.0); if (_ply_pnt_ids.size() > 2) { dist_until_now = _length[_ply_pnt_ids.size() - 2]; @@ -126,21 +137,21 @@ bool Polyline::insertPoint(std::size_t pos, std::size_t pnt_id) _length.insert(_length.begin() + pos_dt, dist_until_now + act_dist); - } else { + } + else + { // insert at arbitrary position within the vector - double dist_until_now (0.0); + double dist_until_now(0.0); if (pos > 1) { dist_until_now = _length[pos - 1]; } double len_seg0(std::sqrt(MathLib::sqrDist( - *_ply_pnts[_ply_pnt_ids[pos - 1]], - *_ply_pnts[pnt_id]))); + *_ply_pnts[_ply_pnt_ids[pos - 1]], *_ply_pnts[pnt_id]))); double len_seg1(std::sqrt(MathLib::sqrDist( - *_ply_pnts[_ply_pnt_ids[pos + 1]], - *_ply_pnts[pnt_id]))); - double update_dist( - len_seg0 + len_seg1 - (_length[pos] - dist_until_now)); + *_ply_pnts[_ply_pnt_ids[pos + 1]], *_ply_pnts[pnt_id]))); + double update_dist(len_seg0 + len_seg1 - + (_length[pos] - dist_until_now)); _length[pos] = dist_until_now + len_seg0; auto it1(_length.begin() + pos_dt + 1); _length.insert(it1, _length[pos] + len_seg1); @@ -173,19 +184,22 @@ void Polyline::removePoint(std::size_t pos) } const std::size_t n_ply_pnt_ids(_ply_pnt_ids.size()); - if (pos == 0) { + if (pos == 0) + { double seg_length(_length[0]); for (std::size_t k(0); k < n_ply_pnt_ids; k++) { _length[k] = _length[k + 1] - seg_length; } _length.pop_back(); - } else { + } + else + { const double len_seg0(_length[pos] - _length[pos - 1]); const double len_seg1(_length[pos + 1] - _length[pos]); _length.erase(_length.begin() + pos_dt); - const double len_new_seg(std::sqrt(MathLib::sqrDist(*_ply_pnts[_ply_pnt_ids[pos - 1]], - *_ply_pnts[_ply_pnt_ids[pos]]))); + const double len_new_seg(std::sqrt(MathLib::sqrDist( + *_ply_pnts[_ply_pnt_ids[pos - 1]], *_ply_pnts[_ply_pnt_ids[pos]]))); double seg_length_diff(len_new_seg - len_seg0 - len_seg1); for (std::size_t k(pos); k < n_ply_pnt_ids; k++) @@ -202,7 +216,7 @@ std::size_t Polyline::getNumberOfPoints() const std::size_t Polyline::getNumberOfSegments() const { - return _ply_pnt_ids.empty() ? 0 : _ply_pnt_ids.size()-1; + return _ply_pnt_ids.empty() ? 0 : _ply_pnt_ids.size() - 1; } bool Polyline::isClosed() const @@ -217,16 +231,16 @@ bool Polyline::isClosed() const bool Polyline::isCoplanar() const { - std::size_t const n_points (_ply_pnt_ids.size()); + std::size_t const n_points(_ply_pnt_ids.size()); if (n_points < 4) { return true; } - GeoLib::Point const& p0 (*this->getPoint(0)); - GeoLib::Point const& p1 (*this->getPoint(1)); - GeoLib::Point const& p2 (*this->getPoint(2)); - for (std::size_t i=3; i<n_points; ++i) + GeoLib::Point const& p0(*this->getPoint(0)); + GeoLib::Point const& p1(*this->getPoint(1)); + GeoLib::Point const& p2(*this->getPoint(2)); + for (std::size_t i = 3; i < n_points; ++i) { if (!MathLib::isCoplanar(p0, p1, p2, *this->getPoint(i))) { @@ -242,7 +256,8 @@ bool Polyline::isCoplanar() const bool Polyline::isPointIDInPolyline(std::size_t pnt_id) const { - return std::find(_ply_pnt_ids.begin(), _ply_pnt_ids.end(), pnt_id) != _ply_pnt_ids.end(); + return std::find(_ply_pnt_ids.begin(), _ply_pnt_ids.end(), pnt_id) != + _ply_pnt_ids.end(); } std::size_t Polyline::getPointID(std::size_t i) const @@ -277,19 +292,19 @@ const Point* Polyline::getPoint(std::size_t i) const return _ply_pnts[_ply_pnt_ids[i]]; } -std::vector<Point*> const& Polyline::getPointsVec () const +std::vector<Point*> const& Polyline::getPointsVec() const { return _ply_pnts; } -double Polyline::getLength (std::size_t k) const +double Polyline::getLength(std::size_t k) const { assert(k < _length.size()); return _length[k]; } -Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> &ply_vec, - double prox) +Polyline* Polyline::constructPolylineFromSegments( + const std::vector<Polyline*>& ply_vec, double prox) { std::size_t nLines = ply_vec.size(); @@ -305,14 +320,14 @@ Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> & while (!local_ply_vec.empty()) { bool ply_found(false); - prox *= prox; // square distance once to save time later + prox *= prox; // square distance once to save time later for (auto it = local_ply_vec.begin(); it != local_ply_vec.end(); ++it) { if (pnt_vec == (*it)->getPointsVec()) { std::size_t nPoints((*it)->getNumberOfPoints()); - //if (new_ply->getPointID(0) == (*it)->getPointID(0)) + // if (new_ply->getPointID(0) == (*it)->getPointID(0)) if (pointsAreIdentical(pnt_vec, new_ply->getPointID(0), (*it)->getPointID(0), prox)) { @@ -331,9 +346,11 @@ Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> & new_ply = tmp; ply_found = true; } - //else if (new_ply->getPointID(0) == (*it)->getPointID(nPoints-1)) + // else if (new_ply->getPointID(0) == + // (*it)->getPointID(nPoints-1)) else if (pointsAreIdentical(pnt_vec, new_ply->getPointID(0), - (*it)->getPointID(nPoints - 1), prox)) + (*it)->getPointID(nPoints - 1), + prox)) { auto* tmp = new Polyline(**it); std::size_t new_ply_size(new_ply->getNumberOfPoints()); @@ -345,12 +362,13 @@ Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> & new_ply = tmp; ply_found = true; } - //else if (new_ply->getPointID(new_ply->getNumberOfPoints()-1) == (*it)->getPointID(0)) - else if (pointsAreIdentical(pnt_vec, - new_ply->getPointID(new_ply-> - getNumberOfPoints() - - 1), - (*it)->getPointID(0), prox)) + // else if (new_ply->getPointID(new_ply->getNumberOfPoints()-1) + // == (*it)->getPointID(0)) + else if (pointsAreIdentical( + pnt_vec, + new_ply->getPointID(new_ply->getNumberOfPoints() - + 1), + (*it)->getPointID(0), prox)) { for (std::size_t k = 1; k < nPoints; k++) { @@ -358,12 +376,13 @@ Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> & } ply_found = true; } - //else if (new_ply->getPointID(new_ply->getNumberOfPoints()-1) == (*it)->getPointID(nPoints-1)) - else if (pointsAreIdentical(pnt_vec, - new_ply->getPointID(new_ply-> - getNumberOfPoints() - - 1), - (*it)->getPointID(nPoints - 1), prox)) + // else if (new_ply->getPointID(new_ply->getNumberOfPoints()-1) + // == (*it)->getPointID(nPoints-1)) + else if (pointsAreIdentical( + pnt_vec, + new_ply->getPointID(new_ply->getNumberOfPoints() - + 1), + (*it)->getPointID(nPoints - 1), prox)) { for (std::size_t k = 1; k < nPoints; k++) { @@ -379,13 +398,15 @@ Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> & } else { - ERR("Error in Polyline::contructPolylineFromSegments() - Line segments use different point vectors."); + ERR("Error in Polyline::contructPolylineFromSegments() - Line " + "segments use different point vectors."); } } if (!ply_found) { - ERR("Error in Polyline::contructPolylineFromSegments() - Not all segments are connected."); + ERR("Error in Polyline::contructPolylineFromSegments() - Not all " + "segments are connected."); delete new_ply; new_ply = nullptr; break; @@ -396,24 +417,28 @@ Polyline* Polyline::constructPolylineFromSegments(const std::vector<Polyline*> & void Polyline::closePolyline() { - if (getNumberOfPoints() < 2) { - ERR("Polyline::closePolyline(): Input polyline needs to be composed of at least three points."); + if (getNumberOfPoints() < 2) + { + ERR("Polyline::closePolyline(): Input polyline needs to be composed of " + "at least three points."); } - if (! isClosed()) { + if (!isClosed()) + { addPoint(getPointID(0)); } } -Location Polyline::getLocationOfPoint (std::size_t k, GeoLib::Point const & pnt) const +Location Polyline::getLocationOfPoint(std::size_t k, + GeoLib::Point const& pnt) const { - assert (k < _ply_pnt_ids.size() - 1); + assert(k < _ply_pnt_ids.size() - 1); - GeoLib::Point const& source (*(_ply_pnts[_ply_pnt_ids[k]])); - GeoLib::Point const& dest (*(_ply_pnts[_ply_pnt_ids[k + 1]])); - long double a[2] = {dest[0] - source[0], dest[1] - source[1]}; // vector - long double b[2] = {pnt[0] - source[0], pnt[1] - source[1]}; // vector + GeoLib::Point const& source(*(_ply_pnts[_ply_pnt_ids[k]])); + GeoLib::Point const& dest(*(_ply_pnts[_ply_pnt_ids[k + 1]])); + long double a[2] = {dest[0] - source[0], dest[1] - source[1]}; // vector + long double b[2] = {pnt[0] - source[0], pnt[1] - source[1]}; // vector - long double det_2x2 (a[0] * b[1] - a[1] * b[0]); + long double det_2x2(a[0] * b[1] - a[1] * b[0]); if (det_2x2 > std::numeric_limits<double>::epsilon()) { @@ -445,16 +470,18 @@ Location Polyline::getLocationOfPoint (std::size_t k, GeoLib::Point const & pnt) } double Polyline::getDistanceAlongPolyline(const MathLib::Point3d& pnt, - const double epsilon_radius) const + const double epsilon_radius) const { double dist(-1.0); double lambda; bool found = false; // loop over all line segments of the polyline - for (std::size_t k = 0; k < getNumberOfSegments(); k++) { + for (std::size_t k = 0; k < getNumberOfSegments(); k++) + { // is the orthogonal projection of the j-th node to the - // line g(lambda) = _ply->getPoint(k) + lambda * (_ply->getPoint(k+1) - _ply->getPoint(k)) - // at the k-th line segment of the polyline, i.e. 0 <= lambda <= 1? + // line g(lambda) = _ply->getPoint(k) + lambda * (_ply->getPoint(k+1) - + // _ply->getPoint(k)) at the k-th line segment of the polyline, i.e. 0 + // <= lambda <= 1? if (MathLib::calcProjPntToLineAndDists(pnt, *getPoint(k), *getPoint(k + 1), lambda, dist) <= epsilon_radius) @@ -471,7 +498,7 @@ double Polyline::getDistanceAlongPolyline(const MathLib::Point3d& pnt, break; } // end if lambda } - } // end line segment loop + } // end line segment loop if (!found) { @@ -485,11 +512,13 @@ Polyline::SegmentIterator::SegmentIterator(Polyline const& polyline, : _polyline(&polyline), _segment_number( static_cast<std::vector<GeoLib::Point*>::size_type>(segment_number)) -{} +{ +} Polyline::SegmentIterator::SegmentIterator(SegmentIterator const& src) : _polyline(src._polyline), _segment_number(src._segment_number) -{} +{ +} Polyline::SegmentIterator& Polyline::SegmentIterator::operator=( SegmentIterator const& rhs) @@ -539,10 +568,13 @@ bool Polyline::SegmentIterator::operator!=(SegmentIterator const& other) const Polyline::SegmentIterator& Polyline::SegmentIterator::operator+=( std::vector<GeoLib::Point>::difference_type n) { - if (n < 0) { + if (n < 0) + { _segment_number -= static_cast<std::vector<GeoLib::Point>::size_type>(-n); - } else { + } + else + { _segment_number += static_cast<std::vector<GeoLib::Point>::size_type>(n); } @@ -564,10 +596,13 @@ Polyline::SegmentIterator Polyline::SegmentIterator::operator+( Polyline::SegmentIterator& Polyline::SegmentIterator::operator-=( std::vector<GeoLib::Point>::difference_type n) { - if (n >= 0) { + if (n >= 0) + { _segment_number -= static_cast<std::vector<GeoLib::Point>::size_type>(n); - } else { + } + else + { _segment_number += static_cast<std::vector<GeoLib::Point>::size_type>(-n); } @@ -586,13 +621,13 @@ Polyline::SegmentIterator Polyline::SegmentIterator::operator-( return t; } -std::ostream& operator<< (std::ostream &os, const Polyline &pl) +std::ostream& operator<<(std::ostream& os, const Polyline& pl) { - pl.write (os); + pl.write(os); return os; } -bool containsEdge (const Polyline& ply, std::size_t id0, std::size_t id1) +bool containsEdge(const Polyline& ply, std::size_t id0, std::size_t id1) { if (id0 == id1) { @@ -603,11 +638,11 @@ bool containsEdge (const Polyline& ply, std::size_t id0, std::size_t id1) { std::swap(id0, id1); } - const std::size_t n (ply.getNumberOfPoints() - 1); + const std::size_t n(ply.getNumberOfPoints() - 1); for (std::size_t k(0); k < n; k++) { - std::size_t ply_pnt0 (ply.getPointID (k)); - std::size_t ply_pnt1 (ply.getPointID (k + 1)); + std::size_t ply_pnt0(ply.getPointID(k)); + std::size_t ply_pnt1(ply.getPointID(k + 1)); if (ply_pnt0 > ply_pnt1) { std::swap(ply_pnt0, ply_pnt1); @@ -639,7 +674,7 @@ bool operator==(Polyline const& lhs, Polyline const& rhs) return true; } -bool pointsAreIdentical(const std::vector<Point*> &pnt_vec, +bool pointsAreIdentical(const std::vector<Point*>& pnt_vec, std::size_t i, std::size_t j, double prox) @@ -650,4 +685,4 @@ bool pointsAreIdentical(const std::vector<Point*> &pnt_vec, } return MathLib::sqrDist(*pnt_vec[i], *pnt_vec[j]) < prox; } -} // end namespace GeoLib +} // end namespace GeoLib diff --git a/GeoLib/PolylineWithSegmentMarker.cpp b/GeoLib/PolylineWithSegmentMarker.cpp index a1debe18a18c081d0e7dd009f597c34b284bc531..67d22cfd00bfe48ddec88b9b57cf657f213dcd9c 100644 --- a/GeoLib/PolylineWithSegmentMarker.cpp +++ b/GeoLib/PolylineWithSegmentMarker.cpp @@ -14,7 +14,8 @@ #include "PolylineWithSegmentMarker.h" -namespace GeoLib { +namespace GeoLib +{ PolylineWithSegmentMarker::PolylineWithSegmentMarker( GeoLib::Polyline const& polyline) : GeoLib::Polyline(polyline), _marker(polyline.getNumberOfSegments(), false) @@ -33,7 +34,8 @@ bool PolylineWithSegmentMarker::isSegmentMarked(std::size_t seg_num) const bool PolylineWithSegmentMarker::addPoint(std::size_t pnt_id) { - if (Polyline::addPoint(pnt_id)) { + if (Polyline::addPoint(pnt_id)) + { _marker.push_back(false); return true; } @@ -42,8 +44,9 @@ bool PolylineWithSegmentMarker::addPoint(std::size_t pnt_id) bool PolylineWithSegmentMarker::insertPoint(std::size_t pos, std::size_t pnt_id) { - if (Polyline::insertPoint(pos, pnt_id)) { - _marker.insert(_marker.begin()+pos, _marker[pos]); + if (Polyline::insertPoint(pos, pnt_id)) + { + _marker.insert(_marker.begin() + pos, _marker[pos]); return true; } return false; diff --git a/GeoLib/Raster.cpp b/GeoLib/Raster.cpp index 456a935a8562a358d7d3cf043d4eb6e6f5a2eeaa..979f414e9c3627ae379e1be78e365a135606de2b 100644 --- a/GeoLib/Raster.cpp +++ b/GeoLib/Raster.cpp @@ -11,30 +11,37 @@ * http://www.opengeosys.org/project/license */ -#include <fstream> - #include "Raster.h" +#include <fstream> + // BaseLib #include "BaseLib/FileTools.h" #include "BaseLib/StringTools.h" - #include "Triangle.h" -namespace GeoLib { - +namespace GeoLib +{ void Raster::refineRaster(std::size_t scaling) { auto* new_raster_data( new double[_header.n_rows * _header.n_cols * scaling * scaling]); - for (std::size_t row(0); row<_header.n_rows; row++) { - for (std::size_t col(0); col<_header.n_cols; col++) { - const std::size_t idx(row*_header.n_cols+col); - for (std::size_t new_row(row*scaling); new_row<(row+1)*scaling; new_row++) { - const std::size_t idx0(new_row*_header.n_cols*scaling); - for (std::size_t new_col(col*scaling); new_col<(col+1)*scaling; new_col++) { - new_raster_data[idx0+new_col] = _raster_data[idx]; + for (std::size_t row(0); row < _header.n_rows; row++) + { + for (std::size_t col(0); col < _header.n_cols; col++) + { + const std::size_t idx(row * _header.n_cols + col); + for (std::size_t new_row(row * scaling); + new_row < (row + 1) * scaling; + new_row++) + { + const std::size_t idx0(new_row * _header.n_cols * scaling); + for (std::size_t new_col(col * scaling); + new_col < (col + 1) * scaling; + new_col++) + { + new_raster_data[idx0 + new_col] = _raster_data[idx]; } } } @@ -45,18 +52,20 @@ void Raster::refineRaster(std::size_t scaling) _header.n_cols *= scaling; _header.n_rows *= scaling; - delete [] new_raster_data; + delete[] new_raster_data; } Raster::~Raster() { - delete [] _raster_data; + delete[] _raster_data; } -double Raster::getValueAtPoint(const MathLib::Point3d &pnt) const +double Raster::getValueAtPoint(const MathLib::Point3d& pnt) const { - if (pnt[0]>=_header.origin[0] && pnt[0]<(_header.origin[0]+(_header.cell_size*_header.n_cols)) && - pnt[1]>=_header.origin[1] && pnt[1]<(_header.origin[1]+(_header.cell_size*_header.n_rows))) + if (pnt[0] >= _header.origin[0] && + pnt[0] < (_header.origin[0] + (_header.cell_size * _header.n_cols)) && + pnt[1] >= _header.origin[1] && + pnt[1] < (_header.origin[1] + (_header.cell_size * _header.n_rows))) { auto cell_x = static_cast<int>( std::floor((pnt[0] - _header.origin[0]) / _header.cell_size)); @@ -65,12 +74,14 @@ double Raster::getValueAtPoint(const MathLib::Point3d &pnt) const // use raster boundary values if node is outside raster due to rounding // errors or floating point arithmetic - cell_x = (cell_x < 0) ? 0 : ((cell_x > static_cast<int>(_header.n_cols)) - ? static_cast<int>(_header.n_cols - 1) - : cell_x); - cell_y = (cell_y < 0) ? 0 : ((cell_y > static_cast<int>(_header.n_rows)) - ? static_cast<int>(_header.n_rows - 1) - : cell_y); + cell_x = (cell_x < 0) ? 0 + : ((cell_x > static_cast<int>(_header.n_cols)) + ? static_cast<int>(_header.n_cols - 1) + : cell_x); + cell_y = (cell_y < 0) ? 0 + : ((cell_y > static_cast<int>(_header.n_rows)) + ? static_cast<int>(_header.n_rows - 1) + : cell_y); const std::size_t index = cell_y * _header.n_cols + cell_x; return _raster_data[index]; @@ -81,11 +92,11 @@ double Raster::getValueAtPoint(const MathLib::Point3d &pnt) const double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const { // position in raster - double const xPos ((pnt[0] - _header.origin[0]) / _header.cell_size); - double const yPos ((pnt[1] - _header.origin[1]) / _header.cell_size); + double const xPos((pnt[0] - _header.origin[0]) / _header.cell_size); + double const yPos((pnt[1] - _header.origin[1]) / _header.cell_size); // raster cell index - double const xIdx (std::floor(xPos)); //carry out computions in double - double const yIdx (std::floor(yPos)); // so not to over- or underflow. + double const xIdx(std::floor(xPos)); // carry out computions in double + double const yIdx(std::floor(yPos)); // so not to over- or underflow. // weights for bilinear interpolation double const xShift = std::abs((xPos - xIdx) - 0.5); @@ -97,13 +108,13 @@ double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const // neighbors to include in interpolation int const xShiftIdx = (xPos - xIdx >= 0.5) ? 1 : -1; int const yShiftIdx = (yPos - yIdx >= 0.5) ? 1 : -1; - std::array<int,4> const x_nb = {{ 0, xShiftIdx, xShiftIdx, 0 }}; - std::array<int,4> const y_nb = {{ 0, 0, yShiftIdx, yShiftIdx }}; + std::array<int, 4> const x_nb = {{0, xShiftIdx, xShiftIdx, 0}}; + std::array<int, 4> const y_nb = {{0, 0, yShiftIdx, yShiftIdx}}; // get pixel values - Eigen::Vector4d pix_val{}; - unsigned no_data_count (0); - for (unsigned j=0; j<4; ++j) + Eigen::Vector4d pix_val{}; + unsigned no_data_count(0); + for (unsigned j = 0; j < 4; ++j) { // check if neighbour pixel is still on the raster, otherwise substitute // a no data value. This also allows the cast to unsigned type. @@ -142,7 +153,7 @@ double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const // new value return weight.dot(pix_val); - } +} bool Raster::isPntOnRaster(MathLib::Point3d const& pnt) const { @@ -153,4 +164,4 @@ bool Raster::isPntOnRaster(MathLib::Point3d const& pnt) const (pnt[1] > _header.origin[1] + (_header.n_rows * _header.cell_size))); } -} // end namespace GeoLib +} // end namespace GeoLib diff --git a/GeoLib/SensorData.cpp b/GeoLib/SensorData.cpp index 728f75c5bdcd1973fecc7ce8b0e1be882d05f9a8..f2464ccc39f91874b95ae21ebbf00ed75cd021c4 100644 --- a/GeoLib/SensorData.cpp +++ b/GeoLib/SensorData.cpp @@ -17,21 +17,24 @@ #include <cstdlib> #include <fstream> +#include "BaseLib/DateTools.h" #include "BaseLib/Logging.h" - #include "BaseLib/StringTools.h" -#include "BaseLib/DateTools.h" -SensorData::SensorData(const std::string &file_name) -: _start(0), _end(0), _step_size(0), _time_unit(TimeStepType::NONE) +SensorData::SensorData(const std::string& file_name) + : _start(0), _end(0), _step_size(0), _time_unit(TimeStepType::NONE) { this->readDataFromFile(file_name); } SensorData::SensorData(std::vector<std::size_t> time_steps) -: _start(time_steps[0]), _end(time_steps[time_steps.size()-1]), _step_size(0), _time_unit(TimeStepType::NONE), _time_steps(time_steps) + : _start(time_steps[0]), + _end(time_steps[time_steps.size() - 1]), + _step_size(0), + _time_unit(TimeStepType::NONE), + _time_steps(time_steps) { - for (std::size_t i=1; i<time_steps.size(); i++) + for (std::size_t i = 1; i < time_steps.size(); i++) { if (time_steps[i - 1] >= time_steps[i]) { @@ -40,8 +43,13 @@ SensorData::SensorData(std::vector<std::size_t> time_steps) } } -SensorData::SensorData(std::size_t first_timestep, std::size_t last_timestep, std::size_t step_size) -: _start(first_timestep), _end(last_timestep), _step_size(step_size), _time_unit(TimeStepType::NONE) +SensorData::SensorData(std::size_t first_timestep, + std::size_t last_timestep, + std::size_t step_size) + : _start(first_timestep), + _end(last_timestep), + _step_size(step_size), + _time_unit(TimeStepType::NONE) { } @@ -53,22 +61,36 @@ SensorData::~SensorData() } } - -void SensorData::addTimeSeries( const std::string &data_name, std::vector<float> *data, const std::string &data_unit_string ) +void SensorData::addTimeSeries(const std::string& data_name, + std::vector<float>* data, + const std::string& data_unit_string) { - this->addTimeSeries(SensorData::convertString2SensorDataType(data_name), data, data_unit_string); + this->addTimeSeries(SensorData::convertString2SensorDataType(data_name), + data, + data_unit_string); } -void SensorData::addTimeSeries(SensorDataType data_name, std::vector<float> *data, const std::string &data_unit_string) +void SensorData::addTimeSeries(SensorDataType data_name, + std::vector<float>* data, + const std::string& data_unit_string) { - if (_step_size>0) { - if (((_end-_start)/_step_size) != data->size()) { - WARN("Warning in SensorData::addTimeSeries() - Lengths of time series does not match number of time steps."); + if (_step_size > 0) + { + if (((_end - _start) / _step_size) != data->size()) + { + WARN( + "Warning in SensorData::addTimeSeries() - Lengths of time " + "series does not match number of time steps."); return; } - } else { - if (data->size() != _time_steps.size()) { - WARN("Warning in SensorData::addTimeSeries() - Lengths of time series does not match number of time steps."); + } + else + { + if (data->size() != _time_steps.size()) + { + WARN( + "Warning in SensorData::addTimeSeries() - Lengths of time " + "series does not match number of time steps."); return; } } @@ -78,9 +100,10 @@ void SensorData::addTimeSeries(SensorDataType data_name, std::vector<float> *dat _data_unit_string.push_back(data_unit_string); } -const std::vector<float>* SensorData::getTimeSeries(SensorDataType time_series_name) const +const std::vector<float>* SensorData::getTimeSeries( + SensorDataType time_series_name) const { - for (std::size_t i=0; i<_vec_names.size(); i++) + for (std::size_t i = 0; i < _vec_names.size(); i++) { if (time_series_name == _vec_names[i]) { @@ -92,9 +115,9 @@ const std::vector<float>* SensorData::getTimeSeries(SensorDataType time_series_n return nullptr; } -int SensorData::readDataFromFile(const std::string &file_name) +int SensorData::readDataFromFile(const std::string& file_name) { - std::ifstream in( file_name.c_str() ); + std::ifstream in(file_name.c_str()); if (!in.is_open()) { @@ -108,7 +131,7 @@ int SensorData::readDataFromFile(const std::string &file_name) /* first line contains field names */ std::getline(in, line); std::list<std::string> fields = BaseLib::splitString(line, '\t'); - std::list<std::string>::const_iterator it (fields.begin()); + std::list<std::string>::const_iterator it(fields.begin()); std::size_t nFields = fields.size(); if (nFields < 2) @@ -116,12 +139,13 @@ int SensorData::readDataFromFile(const std::string &file_name) return 0; } - std::size_t nDataArrays(nFields-1); + std::size_t nDataArrays(nFields - 1); - //create vectors necessary to hold the data - for (std::size_t i=0; i<nDataArrays; i++) + // create vectors necessary to hold the data + for (std::size_t i = 0; i < nDataArrays; i++) { - this->_vec_names.push_back(SensorData::convertString2SensorDataType(*++it)); + this->_vec_names.push_back( + SensorData::convertString2SensorDataType(*++it)); this->_data_unit_string.emplace_back(""); auto* data = new std::vector<float>; this->_data_vecs.push_back(data); @@ -135,7 +159,9 @@ int SensorData::readDataFromFile(const std::string &file_name) { it = fields.begin(); std::size_t pos(it->rfind(".")); - std::size_t current_time_step = (pos == std::string::npos) ? atoi((it++)->c_str()) : BaseLib::strDate2int(*it++); + std::size_t current_time_step = (pos == std::string::npos) + ? atoi((it++)->c_str()) + : BaseLib::strDate2int(*it++); this->_time_steps.push_back(current_time_step); for (std::size_t i = 0; i < nDataArrays; i++) @@ -153,7 +179,7 @@ int SensorData::readDataFromFile(const std::string &file_name) in.close(); this->_start = this->_time_steps[0]; - this->_end = this->_time_steps[this->_time_steps.size()-1]; + this->_end = this->_time_steps[this->_time_steps.size() - 1]; return 1; } @@ -176,7 +202,7 @@ std::string SensorData::convertSensorDataType2String(SensorDataType t) return "Unknown"; } -SensorDataType SensorData::convertString2SensorDataType(const std::string &s) +SensorDataType SensorData::convertString2SensorDataType(const std::string& s) { if (s == "Evaporation" || s == "EVAPORATION") { @@ -192,4 +218,3 @@ SensorDataType SensorData::convertString2SensorDataType(const std::string &s) } return SensorDataType::OTHER; } - diff --git a/GeoLib/SimplePolygonTree.cpp b/GeoLib/SimplePolygonTree.cpp index ff4ef96db30734ae28a6877b2d7e4130aa8a60d5..5cb1fdf382fb804b7d5fc178448a47d15287758a 100644 --- a/GeoLib/SimplePolygonTree.cpp +++ b/GeoLib/SimplePolygonTree.cpp @@ -16,42 +16,52 @@ namespace GeoLib { -SimplePolygonTree::SimplePolygonTree(Polygon * polygon, SimplePolygonTree * parent) : - _node_polygon (polygon), _parent (parent) -{} +SimplePolygonTree::SimplePolygonTree(Polygon* polygon, + SimplePolygonTree* parent) + : _node_polygon(polygon), _parent(parent) +{ +} SimplePolygonTree::~SimplePolygonTree() { - for (auto * child : _children) { + for (auto* child : _children) + { delete child; } } -bool SimplePolygonTree::isPolygonInside (const SimplePolygonTree* polygon_hierarchy) const +bool SimplePolygonTree::isPolygonInside( + const SimplePolygonTree* polygon_hierarchy) const { - return _node_polygon->isPolylineInPolygon(*(polygon_hierarchy->getPolygon())); + return _node_polygon->isPolylineInPolygon( + *(polygon_hierarchy->getPolygon())); } -void SimplePolygonTree::insertSimplePolygonTree (SimplePolygonTree* polygon_hierarchy) +void SimplePolygonTree::insertSimplePolygonTree( + SimplePolygonTree* polygon_hierarchy) { - const Polygon* polygon (polygon_hierarchy->getPolygon()); - bool nfound (true); - for (std::list<SimplePolygonTree*>::const_iterator it (_children.begin()); - it != _children.end() && nfound; ++it) { - if (((*it)->getPolygon())->isPolylineInPolygon (*(polygon))) { - (*it)->insertSimplePolygonTree (polygon_hierarchy); + const Polygon* polygon(polygon_hierarchy->getPolygon()); + bool nfound(true); + for (std::list<SimplePolygonTree*>::const_iterator it(_children.begin()); + it != _children.end() && nfound; + ++it) + { + if (((*it)->getPolygon())->isPolylineInPolygon(*(polygon))) + { + (*it)->insertSimplePolygonTree(polygon_hierarchy); nfound = false; } } - if (nfound) { - _children.push_back (polygon_hierarchy); + if (nfound) + { + _children.push_back(polygon_hierarchy); polygon_hierarchy->setParent(this); } } -const Polygon* SimplePolygonTree::getPolygon () const +const Polygon* SimplePolygonTree::getPolygon() const { return _node_polygon; } -} // end namespace GeoLib +} // end namespace GeoLib diff --git a/GeoLib/Station.cpp b/GeoLib/Station.cpp index 0be09c81eb7945dbcc275fe165da615fb019c0aa..b32c38e8a497d631f3d08bc53402783bc4a6e7a5 100644 --- a/GeoLib/Station.cpp +++ b/GeoLib/Station.cpp @@ -18,7 +18,6 @@ #include <utility> #include "BaseLib/Logging.h" - #include "BaseLib/StringTools.h" namespace GeoLib @@ -26,25 +25,28 @@ namespace GeoLib Station::Station(double x, double y, double z, std::string name) : Point(x, y, z), _name(std::move(name)) -{} +{ +} Station::Station(Point* coords, std::string name) : Point(*coords), _name(std::move(name)) -{} +{ +} Station::Station(Station const& src) : Point(src), _name(src._name), _type(src._type), _station_value(src._station_value) -{} +{ +} Station::~Station() { delete this->_sensor_data; } -Station* Station::createStation(const std::string & line) +Station* Station::createStation(const std::string& line) { Station* station = new Station(); std::list<std::string> fields = BaseLib::splitString(line, '\t'); @@ -53,8 +55,10 @@ Station* Station::createStation(const std::string & line) { auto it = fields.begin(); station->_name = *it; - (*station)[0] = std::strtod((BaseLib::replaceString(",", ".", *(++it))).c_str(), nullptr); - (*station)[1] = std::strtod((BaseLib::replaceString(",", ".", *(++it))).c_str(), nullptr); + (*station)[0] = std::strtod( + (BaseLib::replaceString(",", ".", *(++it))).c_str(), nullptr); + (*station)[1] = std::strtod( + (BaseLib::replaceString(",", ".", *(++it))).c_str(), nullptr); if (++it != fields.end()) { (*station)[2] = std::strtod( diff --git a/GeoLib/StationBorehole.cpp b/GeoLib/StationBorehole.cpp index 43761f69c8b7e118576c46cb2c6a60cef0123bb7..ab69083709b35c286f9d31aeb50976bd5d71d5fb 100644 --- a/GeoLib/StationBorehole.cpp +++ b/GeoLib/StationBorehole.cpp @@ -19,14 +19,12 @@ #include <cstdlib> #include <fstream> +#include "BaseLib/DateTools.h" #include "BaseLib/Logging.h" - #include "BaseLib/StringTools.h" -#include "BaseLib/DateTools.h" namespace GeoLib { - //////////////////////// // The Borehole class // //////////////////////// @@ -54,7 +52,7 @@ StationBorehole::~StationBorehole() } } -StationBorehole* StationBorehole::createStation(const std::string &line) +StationBorehole* StationBorehole::createStation(const std::string& line) { StationBorehole* borehole = new StationBorehole(); std::list<std::string> fields = BaseLib::splitString(line, '\t'); @@ -63,13 +61,17 @@ StationBorehole* StationBorehole::createStation(const std::string &line) { borehole->_name = fields.front(); fields.pop_front(); - (*borehole)[0] = strtod(BaseLib::replaceString(",", ".", fields.front()).c_str(), nullptr); + (*borehole)[0] = strtod( + BaseLib::replaceString(",", ".", fields.front()).c_str(), nullptr); fields.pop_front(); - (*borehole)[1] = strtod(BaseLib::replaceString(",", ".", fields.front()).c_str(), nullptr); + (*borehole)[1] = strtod( + BaseLib::replaceString(",", ".", fields.front()).c_str(), nullptr); fields.pop_front(); - (*borehole)[2] = strtod(BaseLib::replaceString(",", ".", fields.front()).c_str(), nullptr); + (*borehole)[2] = strtod( + BaseLib::replaceString(",", ".", fields.front()).c_str(), nullptr); fields.pop_front(); - borehole->_depth = strtod(BaseLib::replaceString(",", ".", fields.front()).c_str(), nullptr); + borehole->_depth = strtod( + BaseLib::replaceString(",", ".", fields.front()).c_str(), nullptr); fields.pop_front(); if (fields.empty()) { @@ -90,18 +92,18 @@ StationBorehole* StationBorehole::createStation(const std::string &line) return borehole; } -StationBorehole* StationBorehole::createStation(const std::string &name, +StationBorehole* StationBorehole::createStation(const std::string& name, double x, double y, double z, double depth, - const std::string &date) + const std::string& date) { StationBorehole* station = new StationBorehole(); - station->_name = name; - (*station)[0] = x; - (*station)[1] = y; - (*station)[2] = z; + station->_name = name; + (*station)[0] = x; + (*station)[1] = y; + (*station)[2] = z; station->_depth = depth; if (date != "0000-00-00") { @@ -110,7 +112,8 @@ StationBorehole* StationBorehole::createStation(const std::string &name, return station; } -void StationBorehole::addSoilLayer ( double thickness, const std::string &soil_name) +void StationBorehole::addSoilLayer(double thickness, + const std::string& soil_name) { /* // TF - Altmark @@ -132,16 +135,19 @@ void StationBorehole::addSoilLayer ( double thickness, const std::string &soil_n addSoilLayer((*this)[0], (*this)[1], (*this)[2], ""); } - std::size_t idx (_profilePntVec.size()); + std::size_t idx(_profilePntVec.size()); double x((*_profilePntVec[idx - 1])[0]); double y((*_profilePntVec[idx - 1])[1]); double z((*_profilePntVec[0])[2] - thickness); - addSoilLayer (x, y, z, soil_name); + addSoilLayer(x, y, z, soil_name); } -void StationBorehole::addSoilLayer ( double x, double y, double z, const std::string &soil_name) +void StationBorehole::addSoilLayer(double x, + double y, + double z, + const std::string& soil_name) { - _profilePntVec.push_back (new Point (x, y, z)); + _profilePntVec.push_back(new Point(x, y, z)); _soilName.push_back(soil_name); } diff --git a/GeoLib/Surface.cpp b/GeoLib/Surface.cpp index 7fd2b5084342c63614521c276bd5575e5a5e1c60..92dd2b558d56f0eb6ddd5c2eef2b2278a79a477b 100644 --- a/GeoLib/Surface.cpp +++ b/GeoLib/Surface.cpp @@ -7,18 +7,17 @@ * http://www.opengeosys.org/project/license */ -#include <list> - #include "Surface.h" +#include <list> + // GeoLib #include "AABB.h" +#include "AnalyticalGeometry.h" #include "Polygon.h" +#include "Polyline.h" #include "SurfaceGrid.h" -#include "AnalyticalGeometry.h" - #include "Triangle.h" -#include "Polyline.h" namespace GeoLib { diff --git a/GeoLib/SurfaceGrid.cpp b/GeoLib/SurfaceGrid.cpp index c48462e66161f1258a235d516150e07f3b881919..9b05727a9ee22c2ed2e888c296ed80ac2b5b1087 100644 --- a/GeoLib/SurfaceGrid.cpp +++ b/GeoLib/SurfaceGrid.cpp @@ -14,47 +14,50 @@ #include <algorithm> #include <bitset> -#include "BaseLib/Logging.h" - #include "BaseLib/Error.h" - +#include "BaseLib/Logging.h" #include "MathLib/Point3d.h" - #include "Surface.h" #include "Triangle.h" -namespace GeoLib { - -SurfaceGrid::SurfaceGrid(Surface const*const sfc) : - AABB(sfc->getAABB()), _n_steps({{1,1,1}}) +namespace GeoLib +{ +SurfaceGrid::SurfaceGrid(Surface const* const sfc) + : AABB(sfc->getAABB()), _n_steps({{1, 1, 1}}) { // enlarge the bounding, such that the points with maximal coordinates // fits into the grid - for (std::size_t k(0); k<3; ++k) { + for (std::size_t k(0); k < 3; ++k) + { _max_pnt[k] += std::abs(_max_pnt[k]) * 1e-6; - if (std::abs(_max_pnt[k]) < std::numeric_limits<double>::epsilon()) { + if (std::abs(_max_pnt[k]) < std::numeric_limits<double>::epsilon()) + { _max_pnt[k] = (_max_pnt[k] - _min_pnt[k]) * (1.0 + 1e-6); } } - std::array<double, 3> delta{{ _max_pnt[0] - _min_pnt[0], - _max_pnt[1] - _min_pnt[1], _max_pnt[2] - _min_pnt[2] }}; + std::array<double, 3> delta{{_max_pnt[0] - _min_pnt[0], + _max_pnt[1] - _min_pnt[1], + _max_pnt[2] - _min_pnt[2]}}; - if (delta[0] < std::numeric_limits<double>::epsilon()) { + if (delta[0] < std::numeric_limits<double>::epsilon()) + { const double max_delta(std::max(delta[1], delta[2])); _min_pnt[0] -= max_delta * 0.5e-3; _max_pnt[0] += max_delta * 0.5e-3; delta[0] = _max_pnt[0] - _min_pnt[0]; } - if (delta[1] < std::numeric_limits<double>::epsilon()) { + if (delta[1] < std::numeric_limits<double>::epsilon()) + { const double max_delta(std::max(delta[0], delta[2])); _min_pnt[1] -= max_delta * 0.5e-3; _max_pnt[1] += max_delta * 0.5e-3; delta[1] = _max_pnt[1] - _min_pnt[1]; } - if (delta[2] < std::numeric_limits<double>::epsilon()) { + if (delta[2] < std::numeric_limits<double>::epsilon()) + { const double max_delta(std::max(delta[0], delta[1])); _min_pnt[2] -= max_delta * 0.5e-3; _max_pnt[2] += max_delta * 0.5e-3; @@ -64,7 +67,7 @@ SurfaceGrid::SurfaceGrid(Surface const*const sfc) : const std::size_t n_tris(sfc->getNumberOfTriangles()); const std::size_t n_tris_per_cell(5); - std::bitset<3> dim; // all bits are set to zero. + std::bitset<3> dim; // all bits are set to zero. for (std::size_t k(0); k < 3; ++k) { if (std::abs(delta[k]) >= std::numeric_limits<double>::epsilon()) @@ -75,55 +78,74 @@ SurfaceGrid::SurfaceGrid(Surface const*const sfc) : // *** condition: n_tris / n_cells < n_tris_per_cell // where n_cells = _n_steps[0] * _n_steps[1] * _n_steps[2] - // *** with _n_steps[0] = ceil(pow(n_tris*delta[0]*delta[0]/(n_tris_per_cell*delta[1]*delta[2]), 1/3.))); + // *** with _n_steps[0] = + // ceil(pow(n_tris*delta[0]*delta[0]/(n_tris_per_cell*delta[1]*delta[2]), + // 1/3.))); // _n_steps[1] = _n_steps[0] * delta[1]/delta[0], // _n_steps[2] = _n_steps[0] * delta[2]/delta[0] - auto sc_ceil = [](double v){ + auto sc_ceil = [](double v) { return static_cast<std::size_t>(std::ceil(v)); }; - switch (dim.count()) { - case 3: // 3d case - _n_steps[0] = sc_ceil(std::cbrt( - n_tris*delta[0]*delta[0]/(n_tris_per_cell*delta[1]*delta[2]))); - _n_steps[1] = sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0)); - _n_steps[2] = sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0)); - break; - case 2: // 2d cases - if (dim[0] && dim[2]) { // 2d case: xz plane, y = const - _n_steps[0] = sc_ceil(std::sqrt(n_tris*delta[0]/(n_tris_per_cell*delta[2]))); - _n_steps[2] = sc_ceil(_n_steps[0]*delta[2]/delta[0]); - } - else if (dim[0] && dim[1]) { // 2d case: xy plane, z = const - _n_steps[0] = sc_ceil(std::sqrt(n_tris*delta[0]/(n_tris_per_cell*delta[1]))); - _n_steps[1] = sc_ceil(_n_steps[0] * delta[1]/delta[0]); - } - else if (dim[1] && dim[2]) { // 2d case: yz plane, x = const - _n_steps[1] = sc_ceil(std::sqrt(n_tris*delta[1]/(n_tris_per_cell*delta[2]))); - _n_steps[2] = sc_ceil(n_tris * delta[2] / (n_tris_per_cell*delta[1])); - } - break; - case 1: // 1d cases - for (std::size_t k(0); k<3; ++k) { - if (dim[k]) { - _n_steps[k] = sc_ceil(static_cast<double>(n_tris)/n_tris_per_cell); + switch (dim.count()) + { + case 3: // 3d case + _n_steps[0] = + sc_ceil(std::cbrt(n_tris * delta[0] * delta[0] / + (n_tris_per_cell * delta[1] * delta[2]))); + _n_steps[1] = + sc_ceil(_n_steps[0] * std::min(delta[1] / delta[0], 100.0)); + _n_steps[2] = + sc_ceil(_n_steps[0] * std::min(delta[2] / delta[0], 100.0)); + break; + case 2: // 2d cases + if (dim[0] && dim[2]) + { // 2d case: xz plane, y = const + _n_steps[0] = sc_ceil(std::sqrt(n_tris * delta[0] / + (n_tris_per_cell * delta[2]))); + _n_steps[2] = sc_ceil(_n_steps[0] * delta[2] / delta[0]); + } + else if (dim[0] && dim[1]) + { // 2d case: xy plane, z = const + _n_steps[0] = sc_ceil(std::sqrt(n_tris * delta[0] / + (n_tris_per_cell * delta[1]))); + _n_steps[1] = sc_ceil(_n_steps[0] * delta[1] / delta[0]); + } + else if (dim[1] && dim[2]) + { // 2d case: yz plane, x = const + _n_steps[1] = sc_ceil(std::sqrt(n_tris * delta[1] / + (n_tris_per_cell * delta[2]))); + _n_steps[2] = + sc_ceil(n_tris * delta[2] / (n_tris_per_cell * delta[1])); + } + break; + case 1: // 1d cases + for (std::size_t k(0); k < 3; ++k) + { + if (dim[k]) + { + _n_steps[k] = + sc_ceil(static_cast<double>(n_tris) / n_tris_per_cell); + } } - } } // some frequently used expressions to fill the grid vectors - for (std::size_t k(0); k<3; k++) { + for (std::size_t k(0); k < 3; k++) + { _step_sizes[k] = delta[k] / _n_steps[k]; _inverse_step_sizes[k] = 1.0 / _step_sizes[k]; } - _triangles_in_grid_box.resize(_n_steps[0]*_n_steps[1]*_n_steps[2]); + _triangles_in_grid_box.resize(_n_steps[0] * _n_steps[1] * _n_steps[2]); sortTrianglesInGridCells(sfc); } -void SurfaceGrid::sortTrianglesInGridCells(Surface const*const sfc) +void SurfaceGrid::sortTrianglesInGridCells(Surface const* const sfc) { - for (std::size_t l(0); l<sfc->getNumberOfTriangles(); l++) { - if (! sortTriangleInGridCells((*sfc)[l])) { + for (std::size_t l(0); l < sfc->getNumberOfTriangles(); l++) + { + if (!sortTriangleInGridCells((*sfc)[l])) + { Point const& p0(*((*sfc)[l]->getPoint(0))); Point const& p1(*((*sfc)[l]->getPoint(1))); Point const& p2(*((*sfc)[l]->getPoint(2))); @@ -137,7 +159,7 @@ void SurfaceGrid::sortTrianglesInGridCells(Surface const*const sfc) } } -bool SurfaceGrid::sortTriangleInGridCells(Triangle const*const triangle) +bool SurfaceGrid::sortTriangleInGridCells(Triangle const* const triangle) { // compute grid coordinates for each triangle point std::optional<std::array<std::size_t, 3> const> c_p0( @@ -160,20 +182,30 @@ bool SurfaceGrid::sortTriangleInGridCells(Triangle const*const triangle) } // determine interval in grid (grid cells the triangle will be inserted) - std::size_t const i_min(std::min(std::min((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0])); - std::size_t const i_max(std::max(std::max((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0])); - std::size_t const j_min(std::min(std::min((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1])); - std::size_t const j_max(std::max(std::max((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1])); - std::size_t const k_min(std::min(std::min((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2])); - std::size_t const k_max(std::max(std::max((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2])); + std::size_t const i_min( + std::min(std::min((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0])); + std::size_t const i_max( + std::max(std::max((*c_p0)[0], (*c_p1)[0]), (*c_p2)[0])); + std::size_t const j_min( + std::min(std::min((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1])); + std::size_t const j_max( + std::max(std::max((*c_p0)[1], (*c_p1)[1]), (*c_p2)[1])); + std::size_t const k_min( + std::min(std::min((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2])); + std::size_t const k_max( + std::max(std::max((*c_p0)[2], (*c_p1)[2]), (*c_p2)[2])); - const std::size_t n_plane(_n_steps[0]*_n_steps[1]); + const std::size_t n_plane(_n_steps[0] * _n_steps[1]); // insert the triangle into the grid cells - for (std::size_t i(i_min); i<=i_max; i++) { - for (std::size_t j(j_min); j<=j_max; j++) { - for (std::size_t k(k_min); k<=k_max; k++) { - _triangles_in_grid_box[i+j*_n_steps[0]+k*n_plane].push_back(triangle); + for (std::size_t i(i_min); i <= i_max; i++) + { + for (std::size_t j(j_min); j <= j_max; j++) + { + for (std::size_t k(k_min); k <= k_max; k++) + { + _triangles_in_grid_box[i + j * _n_steps[0] + k * n_plane] + .push_back(triangle); } } } @@ -184,13 +216,17 @@ bool SurfaceGrid::sortTriangleInGridCells(Triangle const*const triangle) std::optional<std::array<std::size_t, 3>> SurfaceGrid::getGridCellCoordinates( MathLib::Point3d const& p) const { - std::array<std::size_t, 3> coords{{ - static_cast<std::size_t>((p[0]-_min_pnt[0]) * _inverse_step_sizes[0]), - static_cast<std::size_t>((p[1]-_min_pnt[1]) * _inverse_step_sizes[1]), - static_cast<std::size_t>((p[2]-_min_pnt[2]) * _inverse_step_sizes[2]) - }}; + std::array<std::size_t, 3> coords{ + {static_cast<std::size_t>((p[0] - _min_pnt[0]) * + _inverse_step_sizes[0]), + static_cast<std::size_t>((p[1] - _min_pnt[1]) * + _inverse_step_sizes[1]), + static_cast<std::size_t>((p[2] - _min_pnt[2]) * + _inverse_step_sizes[2])}}; - if (coords[0]>=_n_steps[0] || coords[1]>=_n_steps[1] || coords[2]>=_n_steps[2]) { + if (coords[0] >= _n_steps[0] || coords[1] >= _n_steps[1] || + coords[2] >= _n_steps[2]) + { DBUG( "Computed indices ({:d},{:d},{:d}), max grid cell indices " "({:d},{:d},{:d})", @@ -206,13 +242,16 @@ bool SurfaceGrid::isPointInSurface(MathLib::Point3d const& pnt, { std::optional<std::array<std::size_t, 3>> optional_c( getGridCellCoordinates(pnt)); - if (!optional_c) { + if (!optional_c) + { return false; } std::array<std::size_t, 3> c(optional_c.value()); - std::size_t const grid_cell_idx(c[0]+c[1]*_n_steps[0]+c[2]*_n_steps[0]*_n_steps[1]); - std::vector<Triangle const*> const& triangles(_triangles_in_grid_box[grid_cell_idx]); + std::size_t const grid_cell_idx(c[0] + c[1] * _n_steps[0] + + c[2] * _n_steps[0] * _n_steps[1]); + std::vector<Triangle const*> const& triangles( + _triangles_in_grid_box[grid_cell_idx]); auto const it = std::find_if(triangles.begin(), triangles.end(), [eps, pnt](auto const* triangle) { return triangle->containsPoint(pnt, eps); @@ -220,4 +259,4 @@ bool SurfaceGrid::isPointInSurface(MathLib::Point3d const& pnt, return it != triangles.end(); } -} // end namespace GeoLib +} // end namespace GeoLib diff --git a/GeoLib/Triangle.cpp b/GeoLib/Triangle.cpp index 3de2c15894eb15efbaafc3a828ea67f4212f9de9..e6e37a1c5f91f4527d2f304515ec31d5b70e0aca 100644 --- a/GeoLib/Triangle.cpp +++ b/GeoLib/Triangle.cpp @@ -11,19 +11,19 @@ #include <Eigen/Dense> -#include "Point.h" #include "AnalyticalGeometry.h" - #include "MathLib/GeometricBasics.h" +#include "Point.h" -namespace GeoLib { - -Triangle::Triangle (std::vector<Point *> const &pnt_vec, - std::size_t pnt_a, std::size_t pnt_b, std::size_t pnt_c) : - _pnts(pnt_vec), _pnt_ids( {{pnt_a, pnt_b, pnt_c}} ) +namespace GeoLib +{ +Triangle::Triangle(std::vector<Point*> const& pnt_vec, std::size_t pnt_a, + std::size_t pnt_b, std::size_t pnt_c) + : _pnts(pnt_vec), _pnt_ids({{pnt_a, pnt_b, pnt_c}}) { assert(!_pnts.empty()); - assert (pnt_a < _pnts.size() && pnt_b < _pnts.size() && pnt_c < _pnts.size()); + assert(pnt_a < _pnts.size() && pnt_b < _pnts.size() && + pnt_c < _pnts.size()); } bool Triangle::containsPoint(MathLib::Point3d const& q, double eps) const