Revision 798f79be491b04a7cbc8912f7676dd96f11a7aab authored by Dmitri Naumov on 24 September 2021, 14:05:12 UTC, committed by Dmitry Yu. Naumov on 30 September 2021, 19:47:14 UTC
1 parent 75dcfa9
AnalyticalGeometry.cpp
/**
* \file
* \author Thomas Fischer
* \date 2010-03-17
* \brief Implementation of analytical geometry functions.
*
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#include "AnalyticalGeometry.h"
#include <Eigen/Dense>
#include <algorithm>
#include <cmath>
#include <limits>
#include "BaseLib/StringTools.h"
#include "MathLib/GeometricBasics.h"
#include "PointVec.h"
#include "Polyline.h"
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)
{
return orient2d(const_cast<double*>(a.getCoords()),
const_cast<double*>(b.getCoords()),
const_cast<double*>(c.getCoords()));
}
double getOrientation2dFast(MathLib::Point3d const& a,
MathLib::Point3d const& b,
MathLib::Point3d const& c)
{
return orient2dfast(const_cast<double*>(a.getCoords()),
const_cast<double*>(b.getCoords()),
const_cast<double*>(c.getCoords()));
}
} // namespace ExactPredicates
namespace GeoLib
{
Orientation getOrientation(MathLib::Point3d const& p0,
MathLib::Point3d const& p1,
MathLib::Point3d const& p2)
{
double const orientation = ExactPredicates::getOrientation2d(p0, p1, p2);
if (orientation > 0)
{
return CCW;
}
if (orientation < 0)
{
return CW;
}
return COLLINEAR;
}
Orientation getOrientationFast(MathLib::Point3d const& p0,
MathLib::Point3d const& p1,
MathLib::Point3d const& p2)
{
double const orientation =
ExactPredicates::getOrientation2dFast(p0, p1, p2);
if (orientation > 0)
{
return CCW;
}
if (orientation < 0)
{
return CW;
}
return COLLINEAR;
}
bool parallel(Eigen::Vector3d v, Eigen::Vector3d w)
{
const double eps(std::numeric_limits<double>::epsilon());
double const eps_squared = eps * eps;
// check degenerated cases
if (v.squaredNorm() < eps_squared)
{
return false;
}
if (w.squaredNorm() < eps_squared)
{
return false;
}
v.normalize();
w.normalize();
bool parallel(true);
if (std::abs(v[0] - w[0]) > eps)
{
parallel = false;
}
if (std::abs(v[1] - w[1]) > eps)
{
parallel = false;
}
if (std::abs(v[2] - w[2]) > eps)
{
parallel = false;
}
if (!parallel)
{
parallel = true;
// change sense of direction of v_normalised
v *= -1.0;
// check again
if (std::abs(v[0] - w[0]) > eps)
{
parallel = false;
}
if (std::abs(v[1] - w[1]) > eps)
{
parallel = false;
}
if (std::abs(v[2] - w[2]) > eps)
{
parallel = false;
}
}
return parallel;
}
bool lineSegmentIntersect(GeoLib::LineSegment const& s0,
GeoLib::LineSegment const& s1,
GeoLib::Point& s)
{
GeoLib::Point const& pa{s0.getBeginPoint()};
GeoLib::Point const& pb{s0.getEndPoint()};
GeoLib::Point const& pc{s1.getBeginPoint()};
GeoLib::Point const& pd{s1.getEndPoint()};
if (!isCoplanar(pa, pb, pc, pd))
{
return false;
}
auto const a =
Eigen::Map<Eigen::Vector3d const>(s0.getBeginPoint().getCoords());
auto const b =
Eigen::Map<Eigen::Vector3d const>(s0.getEndPoint().getCoords());
auto const c =
Eigen::Map<Eigen::Vector3d const>(s1.getBeginPoint().getCoords());
auto const d =
Eigen::Map<Eigen::Vector3d const>(s1.getEndPoint().getCoords());
Eigen::Vector3d const v = b - a;
Eigen::Vector3d const w = d - c;
Eigen::Vector3d const qp = c - a;
Eigen::Vector3d const pq = a - c;
double const eps = std::numeric_limits<double>::epsilon();
double const squared_eps = eps * eps;
// handle special cases here to avoid computing intersection numerical
if (qp.squaredNorm() < squared_eps || (d - a).squaredNorm() < squared_eps)
{
s = pa;
return true;
}
if ((c - b).squaredNorm() < squared_eps ||
(d - b).squaredNorm() < squared_eps)
{
s = pb;
return true;
}
auto isLineSegmentIntersectingAB = [&v](Eigen::Vector3d const& ap,
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
// 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
// possible to express c as c(t) = a + t * (b-a) (analog for the
// point d). Since the evaluation of all three coordinate equations
// (x,y,z) have to lead to the same solution for the parameter t it
// is sufficient to evaluate t only once.
// Search id of coordinate with largest absolute value which is will
// be used in the subsequent computations. This prevents division by
// zero in case the line segments are parallel to one of the
// 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))
{
s = pc;
return true;
}
Eigen::Vector3d const ad = d - a;
if (isLineSegmentIntersectingAB(ad, i_max))
{
s = pd;
return true;
}
return false;
}
return false;
}
// general case
const double sqr_len_v(v.squaredNorm());
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);
Eigen::Vector2d rhs{v.dot(qp), w.dot(pq)};
rhs = mat.partialPivLu().solve(rhs);
// 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());
// 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])
{
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]);
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)
{
s[0] = 0.5 * (p0[0] + p1[0]);
s[1] = 0.5 * (p0[1] + p1[1]);
s[2] = 0.5 * (p0[2] + p1[2]);
return true;
}
return false;
}
bool lineSegmentsIntersect(const GeoLib::Polyline* ply,
GeoLib::Polyline::SegmentIterator& seg_it0,
GeoLib::Polyline::SegmentIterator& seg_it1,
GeoLib::Point& intersection_pnt)
{
std::size_t const n_segs(ply->getNumberOfSegments());
// Neighbouring segments always intersects at a common vertex. The algorithm
// checks for intersections of non-neighbouring segments.
for (seg_it0 = ply->begin(); seg_it0 != ply->end() - 2; ++seg_it0)
{
seg_it1 = seg_it0 + 2;
std::size_t const seg_num_0 = seg_it0.getSegmentNumber();
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))
{
return true;
}
}
}
}
return false;
}
void rotatePoints(Eigen::Matrix3d const& rot_mat,
std::vector<GeoLib::Point*>& pnts)
{
rotatePoints(rot_mat, pnts.begin(), pnts.end());
}
Eigen::Matrix3d computeRotationMatrixToXY(Eigen::Vector3d const& n)
{
Eigen::Matrix3d rot_mat = Eigen::Matrix3d::Zero();
// check if normal points already in the right direction
if (n[0] == 0 && n[1] == 0)
{
rot_mat(1, 1) = 1.0;
if (n[2] > 0)
{
// identity matrix
rot_mat(0, 0) = 1.0;
rot_mat(2, 2) = 1.0;
}
else
{
// rotate by pi about the y-axis
rot_mat(0, 0) = -1.0;
rot_mat(2, 2) = -1.0;
}
return rot_mat;
}
// sqrt (n_1^2 + n_3^2)
double const h0(std::sqrt(n[0] * n[0] + n[2] * n[2]));
// In case the x and z components of the normal are both zero the rotation
// to the x-z-plane is not required, i.e. only the rotation in the z-axis is
// required. The angle is either pi/2 or 3/2*pi. Thus the components of
// rot_mat are as follows.
if (h0 < std::numeric_limits<double>::epsilon())
{
rot_mat(0, 0) = 1.0;
if (n[1] > 0)
{
rot_mat(1, 2) = -1.0;
rot_mat(2, 1) = 1.0;
}
else
{
rot_mat(1, 2) = 1.0;
rot_mat(2, 1) = -1.0;
}
return rot_mat;
}
double const h1(1 / n.norm());
// general case: calculate entries of rotation matrix
rot_mat(0, 0) = n[2] / h0;
rot_mat(0, 1) = 0;
rot_mat(0, 2) = -n[0] / h0;
rot_mat(1, 0) = -n[1] * n[0] / h0 * h1;
rot_mat(1, 1) = h0 * h1;
rot_mat(1, 2) = -n[1] * n[2] / h0 * h1;
rot_mat(2, 0) = n[0] * h1;
rot_mat(2, 1) = n[1] * h1;
rot_mat(2, 2) = n[2] * h1;
return rot_mat;
}
Eigen::Matrix3d rotatePointsToXY(std::vector<GeoLib::Point*>& pnts)
{
return rotatePointsToXY(pnts.begin(), pnts.end(), pnts.begin(), pnts.end());
}
std::unique_ptr<GeoLib::Point> triangleLineIntersection(
MathLib::Point3d const& a, MathLib::Point3d const& b,
MathLib::Point3d const& c, MathLib::Point3d const& p,
MathLib::Point3d const& q)
{
auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
auto const vq = Eigen::Map<Eigen::Vector3d const>(q.getCoords());
Eigen::Vector3d const pq = vq - vp;
Eigen::Vector3d const pa = va - vp;
Eigen::Vector3d const pb = vb - vp;
Eigen::Vector3d const pc = vc - vp;
double u = pq.cross(pc).dot(pb);
if (u < 0)
{
return nullptr;
}
double v = pq.cross(pa).dot(pc);
if (v < 0)
{
return nullptr;
}
double w = pq.cross(pb).dot(pa);
if (w < 0)
{
return nullptr;
}
const double denom(1.0 / (u + v + w));
u *= denom;
v *= denom;
w *= denom;
return std::make_unique<GeoLib::Point>(u * a[0] + v * b[0] + w * c[0],
u * a[1] + v * b[1] + w * c[1],
u * a[2] + v * b[2] + w * c[2]);
}
void computeAndInsertAllIntersectionPoints(GeoLib::PointVec& pnt_vec,
std::vector<GeoLib::Polyline*>& plys)
{
auto computeSegmentIntersections = [&pnt_vec](GeoLib::Polyline& poly0,
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)
{
GeoLib::Point s(0.0, 0.0, 0.0, pnt_vec.size());
if (lineSegmentIntersect(*seg0_it, *seg1_it, s))
{
std::size_t const id(
pnt_vec.push_back(new GeoLib::Point(s)));
poly0.insertPoint(seg0_it.getSegmentNumber() + 1, id);
poly1.insertPoint(seg1_it.getSegmentNumber() + 1, id);
}
}
}
};
for (auto it0(plys.begin()); it0 != plys.end(); ++it0)
{
auto it1(it0);
++it1;
for (; it1 != plys.end(); ++it1)
{
computeSegmentIntersections(*(*it0), *(*it1));
}
}
}
GeoLib::Polygon rotatePolygonToXY(GeoLib::Polygon const& polygon_in,
Eigen::Vector3d& plane_normal)
{
// 1 copy all points
auto* polygon_pnts(new std::vector<GeoLib::Point*>);
for (std::size_t k(0); k < polygon_in.getNumberOfPoints(); k++)
{
polygon_pnts->push_back(new GeoLib::Point(*(polygon_in.getPoint(k))));
}
// 2 rotate points
double d_polygon;
std::tie(plane_normal, d_polygon) = GeoLib::getNewellPlane(*polygon_pnts);
Eigen::Matrix3d const rot_mat =
GeoLib::computeRotationMatrixToXY(plane_normal);
GeoLib::rotatePoints(rot_mat, *polygon_pnts);
// 3 set z coord to zero
std::for_each(polygon_pnts->begin(), polygon_pnts->end(),
[](GeoLib::Point* p) { (*p)[2] = 0.0; });
// 4 create new polygon
GeoLib::Polyline rot_polyline(*polygon_pnts);
for (std::size_t k(0); k < polygon_in.getNumberOfPoints(); k++)
{
rot_polyline.addPoint(k);
}
rot_polyline.addPoint(0);
return GeoLib::Polygon(rot_polyline);
}
std::vector<MathLib::Point3d> lineSegmentIntersect2d(
GeoLib::LineSegment const& ab, GeoLib::LineSegment const& cd)
{
GeoLib::Point const& a{ab.getBeginPoint()};
GeoLib::Point const& b{ab.getEndPoint()};
GeoLib::Point const& c{cd.getBeginPoint()};
GeoLib::Point const& d{cd.getEndPoint()};
double const orient_abc(getOrientation(a, b, c));
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))
{
return std::vector<MathLib::Point3d>();
}
// check: (cd) and (ab) are on the same line
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)
{
return {{a, b}};
}
if (MathLib::sqrDist2d(a, d) < eps && MathLib::sqrDist2d(b, c) < eps)
{
return {{a, b}};
}
// 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) {
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]))
{
// check if a in (cd)
if (isPointOnSegment(a[0], c[0], d[0]))
{
return {{a, c}};
}
// check b == c
if (MathLib::sqrDist2d(b, c) < eps)
{
return {{b}};
}
// check if b in (cd)
if (isPointOnSegment(b[0], c[0], d[0]))
{
return {{b, c}};
}
// check d in (ab)
if (isPointOnSegment(d[0], a[0], b[0]))
{
return {{c, d}};
}
std::stringstream err;
err.precision(std::numeric_limits<double>::digits10);
err << ab << " x " << cd;
OGS_FATAL(
"The case of parallel line segments ({:s}) is not handled yet. "
"Aborting.",
err.str());
}
// check if d in (ab)
if (isPointOnSegment(d[0], a[0], b[0]))
{
// check if a in (cd)
if (isPointOnSegment(a[0], c[0], d[0]))
{
return {{a, d}};
}
// check if b==d
if (MathLib::sqrDist2d(b, d) < eps)
{
return {{b}};
}
// check if b in (cd)
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]))
{
return {{c, d}};
}
std::stringstream err;
err.precision(std::numeric_limits<double>::digits10);
err << ab << " x " << cd;
OGS_FATAL(
"The case of parallel line segments ({:s}) is not handled yet. "
"Aborting.",
err.str());
}
return std::vector<MathLib::Point3d>();
}
// precondition: points a, b, c are collinear
// the function checks if the point c is onto the line segment (a,b)
auto isCollinearPointOntoLineSegment = [](MathLib::Point3d const& a,
MathLib::Point3d const& b,
MathLib::Point3d const& c) {
if (b[0] - a[0] != 0)
{
double const t = (c[0] - a[0]) / (b[0] - a[0]);
return 0.0 <= t && t <= 1.0;
}
if (b[1] - a[1] != 0)
{
double const t = (c[1] - a[1]) / (b[1] - a[1]);
return 0.0 <= t && t <= 1.0;
}
if (b[2] - a[2] != 0)
{
double const t = (c[2] - a[2]) / (b[2] - a[2]);
return 0.0 <= t && t <= 1.0;
}
return false;
};
if (orient_abc == 0.0)
{
if (isCollinearPointOntoLineSegment(a, b, c))
{
return {{c}};
}
return std::vector<MathLib::Point3d>();
}
if (orient_abd == 0.0)
{
if (isCollinearPointOntoLineSegment(a, b, d))
{
return {{d}};
}
return std::vector<MathLib::Point3d>();
}
// 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))
{
return std::vector<MathLib::Point3d>();
}
// at this point it is sure that there is an intersection and the system of
// 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];
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])}}}};
}
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)
{
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) {
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();
findNextSegment(seg_beg_pnt, sub_segments, seg_it);
while (seg_it != sub_segments.end())
{
MathLib::Point3d& new_seg_beg_pnt(seg_it->getEndPoint());
seg_it++;
if (seg_it != sub_segments.end())
{
findNextSegment(new_seg_beg_pnt, sub_segments, seg_it);
}
}
}
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;
return rot_mat;
}
Eigen::Matrix3d compute3DRotationMatrixToX(Eigen::Vector3d const& v)
{
// a vector on the plane
Eigen::Vector3d yy = Eigen::Vector3d::Zero();
auto const eps = std::numeric_limits<double>::epsilon();
if (std::abs(v[0]) > 0.0 && std::abs(v[1]) + std::abs(v[2]) < eps)
{
yy[2] = 1.0;
}
else if (std::abs(v[1]) > 0.0 && std::abs(v[0]) + std::abs(v[2]) < eps)
{
yy[0] = 1.0;
}
else if (std::abs(v[2]) > 0.0 && std::abs(v[0]) + std::abs(v[1]) < eps)
{
yy[1] = 1.0;
}
else
{
for (unsigned i = 0; i < 3; i++)
{
if (std::abs(v[i]) > 0.0)
{
yy[i] = -v[i];
break;
}
}
}
// z"_vec
Eigen::Vector3d const zz = v.cross(yy).normalized();
// y"_vec
yy = zz.cross(v).normalized();
Eigen::Matrix3d rot_mat;
rot_mat.row(0) = v;
rot_mat.row(1) = yy;
rot_mat.row(2) = zz;
return rot_mat;
}
} // end namespace GeoLib
![swh spinner](/static/img/swh-spinner.gif)
Computing file changes ...