Commit 002021cc authored by Dennis Gläser's avatar Dennis Gläser
Browse files

outsource algorithms

parent 0216dee5
add_subdirectory(algorithms)
install(FILES
intersect.hh
intersectiontraits.hh
......
install(FILES
intersection_disk_disk.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/frackit/intersection/algorithms)
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Contains the intersection algorithm
* between a lateral cylinder surface and a disk.
*/
#ifndef FRACKIT_CYLINDERSURFACE_DISK_INTERSECTION_HH
#define FRACKIT_CYLINDERSURFACE_DISK_INTERSECTION_HH
#include <vector>
#include <algorithm>
#include <frackit/common/utilities.hh>
#include <frackit/intersection/intersectiontraits.hh>
#include <frackit/intersection/emptyintersection.hh>
namespace Frackit {
namespace IntersectionAlgorithms {
//! Intersect a lateral cylinder surface and a disk
//! The result can be:
//! - an ellipse
//! - ellipse arc(s)
//! - segment(s)
//! - touching points
template<class ctype>
Intersection< CylindricalSurface<ctype>, Disk<ctype> >
intersect_cylinderSurface_disk(const CylindricalSurface<ctype>& cylSurface,
const Disk<ctype>& disk,
ctype eps)
{
// possible result geometries
using ResultType = Intersection< CylindricalSurface<ctype>, Disk<ctype> >;
using Point = Frackit::Point<ctype, 3>;
using Segment = Frackit::Segment<ctype, 3>;
using EllipseArc = Frackit::EllipseArc<ctype, 3>;
using Ellipse = Frackit::Ellipse<ctype, 3>;
using Direction = typename Ellipse::Direction;
using Vector = typename Direction::Vector;
// Intersect the disk with the cylinder
const auto diskFace = OCCUtilities::getShape(disk);
const auto cylinder = OCCUtilities::getShape(cylSurface.cylinder());
const auto cylLateralFace = OCCUtilities::getShape(cylSurface);
const auto containedFaceShape = OCCUtilities::intersect(diskFace, cylinder, 0.1*eps);
auto containedFaces = OCCUtilities::getFaces(containedFaceShape);
assert(containedFaces.size() <= 1);
// intersect the disk boundary with the cylinder surface (detect touching points)
const auto ellipseShape = OCCUtilities::getShape(disk.boundingEllipse());
const auto ellipseCut = OCCUtilities::cut(ellipseShape, cylLateralFace, 0.1*eps);
const auto ellipseCutVertices = OCCUtilities::getVertices(ellipseCut);
// only touching points are possible
if (containedFaces.size() == 0)
{
// find vertices on surface (avoid duplicates)
std::vector<Point> touchingPoints;
for (const auto& v : ellipseCutVertices)
{
auto p = OCCUtilities::point(v);
if (cylSurface.contains(p, eps)
&& std::none_of(touchingPoints.begin(),
touchingPoints.end(),
[&p, eps] (const auto& pnt) { return pnt.isEqual(p, eps); }))
touchingPoints.emplace_back(std::move(p));
}
return ResultType({touchingPoints.begin(), touchingPoints.end()});
}
// there are intersection edges, get orientation of the geometries
gp_Vec dn(OCCUtilities::direction(disk.normal()));
gp_Vec ca(OCCUtilities::direction(cylSurface.direction()));
const bool diskIsParallel = dn.IsNormal(ca, Precision<ctype>::angular());
const bool diskIsOrthogonal = dn.IsParallel(ca, Precision<ctype>::angular());
// ... and (maybe) the ellipse on which the results live
Ellipse infEllipse;
if (diskIsOrthogonal)
infEllipse = Ellipse(disk.supportingPlane().projection(cylSurface.centerSegment().source()),
cylSurface.base1(), cylSurface.base2(),
cylSurface.radius(), cylSurface.radius());
else if (!diskIsParallel)
{
Vector cylDir(cylSurface.direction());
cylDir *= disk.majorAxisLength();
const auto& diskPlane = disk.supportingPlane();
Direction majAxis(Vector(disk.center(), diskPlane.projection(disk.center() + cylDir)));
Direction minAxis(crossProduct(Vector(disk.normal()), Vector(majAxis)));
if (!isRightHandSystem(Vector(majAxis), Vector(minAxis), Vector(disk.normal())))
majAxis.invert();
using std::cos;
const ctype majAxisLength = cylSurface.radius()/cos(dn.Angle(ca));
const auto& cylAxisLine = cylSurface.centerSegment().supportingLine();
const auto center = std::get<Point>(intersect(diskPlane, cylAxisLine, eps));
infEllipse = Ellipse(center, majAxis, minAxis, majAxisLength, cylSurface.radius());
}
// intersect the wire of the contained disk with cylinder surface
const auto containedFaceWires = OCCUtilities::getWires(containedFaceShape);
assert(containedFaceWires.size() == 1);
const auto wireIntersection = OCCUtilities::intersect(containedFaceWires[0], cylLateralFace, 0.1*eps);
const auto wiresOnCylSurface = OCCUtilities::getWires(wireIntersection);
// parse each wire into its basic geometry
std::vector<EllipseArc> isArcs;
std::vector<Segment> isSegments;
for (const auto& wire : wiresOnCylSurface)
{
const auto wireEdges = OCCUtilities::getEdges(wire);
const auto corners = OCCUtilities::getBoundaryVertices(wireEdges);
const auto& c1 = corners.first[1] == 0 ? TopExp::FirstVertex(wireEdges[corners.first[0]])
: TopExp::LastVertex(wireEdges[corners.first[0]]);
const auto& c2 = corners.second[1] == 0 ? TopExp::FirstVertex(wireEdges[corners.second[0]])
: TopExp::LastVertex(wireEdges[corners.second[0]]);
const auto p1 = OCCUtilities::point(c1);
const auto p2 = OCCUtilities::point(c2);
// if both corners are the same the full ellipse is the intersection
if (p1.isEqual(p2, eps))
return ResultType({ infEllipse });
if (!diskIsParallel)
{
// get corners of the wire
EllipseArc arc1(infEllipse, p1, p2);
EllipseArc arc2(infEllipse, p2, p1);
const auto center1 = OCCUtilities::point(arc1.getPoint(0.5));
const auto center2 = OCCUtilities::point(arc2.getPoint(0.5));
// select the arc whose center is on the set of given edges
unsigned int resultArcIndex = 0;
for (const auto& edge : wireEdges)
{
// get unbounded curve and parameter bounds (umin, umax), then trim
Standard_Real uMin, uMax;
Handle(Geom_Curve) curve = BRep_Tool::Curve(edge, uMin, uMax);
curve = new Geom_TrimmedCurve(curve, uMin, uMax);
GeomAPI_ProjectPointOnCurve c1OnCurve(center1, curve);
GeomAPI_ProjectPointOnCurve c2OnCurve(center2, curve);
if (c1OnCurve.LowerDistance() < eps) { resultArcIndex = 1; break; }
if (c2OnCurve.LowerDistance() < eps) { resultArcIndex = 2; break; }
}
assert(resultArcIndex != 0);
if (resultArcIndex == 1) isArcs.emplace_back( std::move(arc1) );
else isArcs.emplace_back( std::move(arc2) );
}
else
isSegments.emplace_back(p1, p2);
}
// There might still be intersection points (avoid duplicates!)
std::vector<Point> isPoints;
if (wiresOnCylSurface.size() < 2)
{
for (const auto& v : ellipseCutVertices)
{
auto p = OCCUtilities::point(v);
if (cylSurface.contains(p, eps))
{
const auto onSeg = std::any_of(isSegments.begin(),
isSegments.end(),
[&p, eps] (const auto& seg) { return seg.contains(p, eps); });
const auto onArc = std::any_of(isArcs.begin(),
isArcs.end(),
[&p, eps] (const auto& arc) { return arc.contains(p, eps); });
if (!onSeg && !onArc && std::none_of(isPoints.begin(),
isPoints.end(),
[&p, eps] (const auto isP) { return isP.isEqual(p, eps); }))
isPoints.emplace_back(std::move(p));
}
}
}
ResultType result;
for (auto&& p : isPoints)
result.emplace_back(std::move(p));
for (auto&& s : isSegments)
result.emplace_back(std::move(s));
for (auto&& arc : isArcs)
result.emplace_back(std::move(arc));
return result.empty() ? ResultType({ EmptyIntersection<3>() }) : result;
}
} // end namespace IntersectionAlgorithms
} // end namespace Frackit
#endif // FRACKIT_CYLINDERSURFACE_DISK_INTERSECTION_HH
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Contains the intersection algorithm
* between two disks.
*/
#ifndef FRACKIT_DISK_DISK_INTERSECTION_HH
#define FRACKIT_DISK_DISK_INTERSECTION_HH
#include <variant>
// Geometries of intersections
#include <frackit/geometry/point.hh>
#include <frackit/geometry/line.hh>
#include <frackit/geometry/segment.hh>
#include <frackit/geometry/plane.hh>
#include <frackit/geometry/disk.hh>
#include <frackit/intersection/intersectiontraits.hh>
#include <frackit/intersection/emptyintersection.hh>
#include "intersection_segment_segment.hh"
#include "intersection_plane_plane.hh"
#include "intersection_disk_line.hh"
namespace Frackit {
namespace IntersectionAlgorithms {
//! Intersect two disks
//! The result can be:
//! - a surface bounded by segments and/or elliptical arcs
//! - a segment
//! - a point
//! - no intersection
template<class ctype>
Intersection< Disk<ctype>, Disk<ctype> >
intersect_disk_disk(const Disk<ctype>& disk1,
const Disk<ctype>& disk2,
ctype eps)
{
using ResultType = Intersection< Disk<ctype>, Disk<ctype> >;
static constexpr int worldDim = 3;
// first intersect the supporting planes
std::string isGeomName;
const auto planeIS = intersect_plane_plane(disk1.supportingPlane(), disk2.supportingPlane(), eps);
// if the planes don't intersect, the disks don't either
if (std::holds_alternative<EmptyIntersection<worldDim>>(planeIS))
return ResultType( EmptyIntersection<worldDim>() );
// if the result is a line, find the possible segment or point intersection
else if (std::holds_alternative<Line<ctype, worldDim>>(planeIS))
{
const auto& isLine = std::get<Line<ctype, worldDim>>(planeIS);
const auto is1 = intersect_disk_line(disk1, isLine, eps);
const auto is2 = intersect_disk_line(disk2, isLine, eps);
// both disks intersect the support plane of the other disks
if (std::holds_alternative<Segment<ctype, worldDim>>(is1)
&& std::holds_alternative<Segment<ctype, worldDim>>(is2))
{
const auto& seg1 = std::get<Segment<ctype, worldDim>>(is1);
const auto& seg2 = std::get<Segment<ctype, worldDim>>(is2);
const auto segmentIS = intersect_segment_segment(seg1, seg2, eps);
if (std::holds_alternative<Segment<ctype, worldDim>>(segmentIS))
return ResultType( std::get<Segment<ctype, worldDim>>(segmentIS) );
if (std::holds_alternative<Point<ctype, worldDim>>(segmentIS))
return ResultType( std::get<Point<ctype, worldDim>>(segmentIS) );
else
return ResultType( EmptyIntersection<worldDim>() );
}
// one of the disks might still touch the other disk
if (std::holds_alternative<Segment<ctype, worldDim>>(is1)
&& std::holds_alternative<Point<ctype, worldDim>>(is2))
{
const auto& seg1 = std::get<Segment<ctype, worldDim>>(is1);
const auto& p2 = std::get<Point<ctype, worldDim>>(is2);
if (seg1.contains(p2, eps))
return ResultType( p2 );
}
else if (std::holds_alternative<Point<ctype, worldDim>>(is1)
&& std::holds_alternative<Segment<ctype, worldDim>>(is2))
{
const auto& p1 = std::get<Point<ctype, worldDim>>(is1);
const auto& seg2 = std::get<Segment<ctype, worldDim>>(is2);
if (seg2.contains(p1, eps))
return ResultType( p1 );
}
// intersection is empty
return ResultType( EmptyIntersection<worldDim>() );
}
// if the result is a plane, find the intersection surface
else if (std::holds_alternative<Plane<ctype, worldDim>>(planeIS))
throw std::runtime_error(std::string("NotImplemented: planar disk-disk intersections"));
throw std::runtime_error(std::string("Unexpected plane-plane intersection result"));
}
} // end namespace IntersectionAlgorithms
} // end namespace Frackit
#endif // FRACKIT_DISK_DISK_INTERSECTION_HH
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Contains the intersection algorithm
* between a disk and a line.
*/
#ifndef FRACKIT_DISK_LINE_INTERSECTION_HH
#define FRACKIT_DISK_LINE_INTERSECTION_HH
#include <variant>
#include <frackit/geometry/disk.hh>
#include <frackit/geometry/line.hh>
#include <frackit/common/utilities.hh>
#include <frackit/intersection/intersectiontraits.hh>
#include <frackit/intersection/emptyintersection.hh>
#include "intersection_plane_line.hh"
namespace Frackit {
namespace IntersectionAlgorithms {
//! Intersect a disk and a line
//! The result can be:
//! - a segment
//! - a point
//! - no intersection
template<class ctype>
Intersection< Disk<ctype>, Line<ctype, 3> >
intersect_disk_line(const Disk<ctype>& disk,
const Line<ctype, 3>& line,
ctype eps)
{
using ResultType = Intersection< Disk<ctype>, Line<ctype, 3> >;
// intersect the line with the plane
const auto planeLineIS = intersect_plane_line(disk.supportingPlane(), line, eps);
// empty result
if (std::holds_alternative<EmptyIntersection<3>>(planeLineIS))
return ResultType( EmptyIntersection<3>() );
// potential point intersection
else if (std::holds_alternative<Point<ctype, 3>>(planeLineIS))
{
// check if point is on the disk
const auto p = std::get<Point<ctype, 3>>(planeLineIS);
if (disk.contains(p, Precision<ctype>::confusion(), false))
return ResultType( p );
else
return ResultType( EmptyIntersection<3>() );
}
// potential segment intersection
else if (std::holds_alternative<Line<ctype, 3>>(planeLineIS))
{
// Intersection might be a segment, a touching point on the rim or empty
// Use the BRep algorithm package to determine the part of the line on the ellipse.
// For this, we create a segment of the line in the neighborhood of ellipse which is
// large enough to cover the entire disk in case they intersect.
auto source = line.projection(disk.center());
auto target = source;
auto dir = Vector<ctype, 3>(line.direction());
dir *= disk.majorAxisLength()*20.0;
source += dir;
target -= dir;
// find the part of the segment on the disk
const auto segmentShape = OCCUtilities::makeEdge(source, target);
const auto diskShape = OCCUtilities::getShape(disk);
const auto isShape = OCCUtilities::intersect(segmentShape, diskShape, 0.1*eps);
const auto isEdges = OCCUtilities::getEdges(isShape);
assert(isEdges.size() <= 1);
if (isEdges.size() == 1)
{
const auto vertices = OCCUtilities::getVertices(isEdges[0]);
assert(vertices.size() == 2);
return ResultType( Segment<ctype, 3>(OCCUtilities::point(vertices[0]),
OCCUtilities::point(vertices[1])) );
}
// the line can still touch the ellipse on the rim
const auto cutShape = OCCUtilities::cut(segmentShape, diskShape, 0.1*eps);
const auto vertices = OCCUtilities::getVertices(cutShape);
assert(vertices.size() == 2 || vertices.size() == 4);
if (vertices.size() == 2)
return ResultType( EmptyIntersection<3>() );
else
{
// find the new point
for (const auto& v : vertices)
{
const auto curPoint = OCCUtilities::point(v);
if (!curPoint.isEqual(source, eps) && !curPoint.isEqual(target, eps))
return ResultType( curPoint );
}
}
throw std::runtime_error(std::string("Unexpected code behaviour"));
}
else
throw std::runtime_error(std::string("Unexpected Plane-Line intersection result"));
}
} // end namespace IntersectionAlgorithms
} // end namespace Frackit
#endif // FRACKIT_INTERSECT_HH
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Contains the intersection algorithm
* between a plane and a line.
*/
#ifndef FRACKIT_PLANE_LINE_INTERSECTION_HH
#define FRACKIT_PLANE_LINE_INTERSECTION_HH
#include <Geom_Line.hxx>
#include <Geom_Surface.hxx>
#include <GeomAPI_IntCS.hxx>
#include <Standard_Handle.hxx>
#include <frackit/geometry/plane.hh>
#include <frackit/geometry/line.hh>
#include <frackit/common/utilities.hh>
#include <frackit/intersection/intersectiontraits.hh>
#include <frackit/intersection/emptyintersection.hh>
namespace Frackit {
namespace IntersectionAlgorithms {
//! Intersect a plane and a line
//! The result can be:
//! - a line
//! - a point
//! - no intersection
template<class ctype>
Intersection< Plane<ctype, 3>, Line<ctype, 3> >
intersect_plane_line(const Plane<ctype, 3>& plane,
const Line<ctype, 3>& line,
ctype eps)
{
using ResultType = Intersection< Plane<ctype, 3>, Line<ctype, 3> >;
// if the line is not parallel to the plane, the result can only be a point
const auto lineDir = OCCUtilities::direction(line.direction());
const auto planeNormal = OCCUtilities::direction(plane.normal());
if (!lineDir.IsNormal(planeNormal, Precision<ctype>::angular()))
{
const auto linePoint = OCCUtilities::point(line.supportingPoint());
const auto planePoint = OCCUtilities::point(plane.supportingPoint());
// let the the geometry package compute the intersection
Handle(Geom_Surface) planeHandle = new Geom_Plane(planePoint, planeNormal);
Handle(Geom_Line) lineHandle = new Geom_Line(linePoint, lineDir);
GeomAPI_IntCS interSection(lineHandle, planeHandle);
if (!interSection.IsDone())
throw std::runtime_error(std::string("Could not perform disk-line intersection"));
assert(interSection.NbSegments() == 0);
assert(interSection.NbPoints() == 1);
// check if point is on the disk
const auto p = OCCUtilities::point(interSection.Point(1));
if (plane.contains(p, eps))
return ResultType( p );
else
return ResultType( EmptyIntersection<3>() );
}
// The line is parallel. If the distance is > eps, there is no intersection
const auto d = Vector<ctype, 3>(line.supportingPoint(),
plane.projection(line.supportingPoint()));
if (d.length() > eps)
return ResultType( EmptyIntersection<3>() );
// the intersection is the line
return ResultType( line );
}
} // end namespace IntersectionAlgorithms
} // end namespace Frackit
#endif // FRACKIT_PLANE_LINE_INTERSECTION_HH
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *