New Upstream Snapshot - lib2geom
Ready changes
Summary
Merged new upstream version: 1.2.2+git20230107.1.6f7dfdc (was: 1.2.2).
Resulting package
Built on 2023-01-19T11:30 (took 6m31s)
The resulting binary packages can be installed (if you have the apt repository enabled) by running one of:
apt install -t fresh-snapshots lib2geom-devapt install -t fresh-snapshots lib2geom1.2.0-dbgsymapt install -t fresh-snapshots lib2geom1.2.0
Lintian Result
- lib2geom-dev_1.2.2+git20230107.1.6f7dfdc-1~jan+nus1_amd64.deb
- lib2geom1.2.0-dbgsym_1.2.2+git20230107.1.6f7dfdc-1~jan+nus1_amd64.deb
- lib2geom1.2.0_1.2.2+git20230107.1.6f7dfdc-1~jan+nus1_amd64.deb
- lib2geom_1.2.2+git20230107.1.6f7dfdc-1~jan+nus1.dsc
- lib2geom_1.2.2+git20230107.1.6f7dfdc-1~jan+nus1_amd64.buildinfo
- lib2geom_1.2.2+git20230107.1.6f7dfdc-1~jan+nus1_amd64.changes
Diff
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
deleted file mode 100644
index 9ef9005..0000000
--- a/.gitlab-ci.yml
+++ /dev/null
@@ -1,32 +0,0 @@
-# Dependencies are managed in the Dockerfile in the inkscape-ci-docker
-# Git repository. Recycle the one for the master Inkscape branch for now
-image: registry.gitlab.com/inkscape/inkscape-ci-docker/master
-
-variables:
- GIT_DEPTH: "10"
- GIT_SUBMODULE_STRATEGY: recursive
-
-
-2geom:
- stage: build
- before_script:
- - mkdir -p ccache
- - export CCACHE_BASEDIR=${PWD}
- - export CCACHE_DIR=${PWD}/ccache
- script:
- - mkdir -p opt
- - cmake . -DCMAKE_C_COMPILER_LAUNCHER=ccache
- -DCMAKE_CXX_COMPILER_LAUNCHER=ccache
- -DCMAKE_CXX_FLAGS="-fsanitize=address -fno-omit-frame-pointer"
- -DCMAKE_BUILD_TYPE=Debug
- -DCMAKE_INSTALL_PREFIX=./opt
- -D2GEOM_BOOST_PYTHON=ON
- - make -j3
- - make -j3 test
- - make -j3 perf
- - make py2geom
- - make install
- - ( cd tests/dependent-project ; bash build-and-run.sh )
- cache:
- paths:
- - ccache/
diff --git a/AUTHORS.md b/AUTHORS.md
index 6297026..2a54918 100644
--- a/AUTHORS.md
+++ b/AUTHORS.md
@@ -21,7 +21,6 @@
- Jelle R. Moulder
- Alvin Penner
- Jan Pulmann
- - Rafał M. Siejakowski
- Michael G. Sloan
- Aaron Spike
- Michael Wybrow
diff --git a/CMakeLists.txt b/CMakeLists.txt
index d2b06b9..f2de358 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -3,8 +3,8 @@ cmake_policy(SET CMP0054 NEW)
set(2GEOM_MAJOR_VERSION 1)
-set(2GEOM_MINOR_VERSION 2)
-set(2GEOM_PATCH_VERSION 2)
+set(2GEOM_MINOR_VERSION 1)
+set(2GEOM_PATCH_VERSION 0)
set(2GEOM_VERSION ${2GEOM_MAJOR_VERSION}.${2GEOM_MINOR_VERSION}.${2GEOM_PATCH_VERSION})
set(2GEOM_ABI_VERSION ${2GEOM_MAJOR_VERSION}.${2GEOM_MINOR_VERSION}.0)
@@ -17,6 +17,7 @@ project(lib2geom
)
set(2GEOM_INCLUDE_DIR "${CMAKE_CURRENT_LIST_DIR}/include" CACHE INTERNAL "")
+include_directories("${CMAKE_CURRENT_LIST_DIR}/src/2geom") # for private headers/template support.
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
@@ -127,7 +128,6 @@ add_custom_target(uninstall_${PROJECT_NAME}
add_subdirectory(src)
-
install(EXPORT 2geom_targets
FILE 2GeomTargets.cmake
NAMESPACE 2Geom::
diff --git a/NEWS.md b/NEWS.md
index e142f90..4d3ccb9 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -1,14 +1,3 @@
-lib2geom v1.2.2
-===============
-
-2Geom v1.2.2 is a point release providing minor bugfixes.
-Its ABI is compatible with v1.2.
-
-Changes:
-
-- Do not raise assertions when intersecting `Geom::EllipticalArc`s.
-
-
lib2geom v1.1.0
===============
diff --git a/debian/changelog b/debian/changelog
index 4eb0ef4..6df565d 100644
--- a/debian/changelog
+++ b/debian/changelog
@@ -1,3 +1,9 @@
+lib2geom (1.2.2+git20230107.1.6f7dfdc-1) UNRELEASED; urgency=low
+
+ * New upstream snapshot.
+
+ -- Debian Janitor <janitor@jelmer.uk> Thu, 19 Jan 2023 11:25:32 -0000
+
lib2geom (1.2.2-3) unstable; urgency=medium
* Upload to unstable.
diff --git a/include/2geom/bezier-curve.h b/include/2geom/bezier-curve.h
index 2fe3d9c..754c9cc 100644
--- a/include/2geom/bezier-curve.h
+++ b/include/2geom/bezier-curve.h
@@ -105,11 +105,12 @@ public:
Point initialPoint() const override { return inner.at0(); }
Point finalPoint() const override { return inner.at1(); }
bool isDegenerate() const override;
- bool isLineSegment() const override { return size() == 2; }
+ bool isLineSegment() const override;
void setInitial(Point const &v) override { setPoint(0, v); }
void setFinal(Point const &v) override { setPoint(order(), v); }
Rect boundsFast() const override { return *bounds_fast(inner); }
Rect boundsExact() const override { return *bounds_exact(inner); }
+ void expandToTransformed(Rect &bbox, Affine const &transform) const override;
OptRect boundsLocal(OptInterval const &i, unsigned deg) const override {
if (!i) return OptRect();
if(i->min() == 0 && i->max() == 1) return boundsFast();
@@ -122,9 +123,9 @@ public:
Curve *duplicate() const override {
return new BezierCurve(*this);
}
- Curve *portion(Coord f, Coord t) const override {
- return new BezierCurve(Geom::portion(inner, f, t));
- }
+
+ Curve *portion(Coord f, Coord t) const override;
+
Curve *reverse() const override {
return new BezierCurve(Geom::reverse(inner));
}
@@ -250,7 +251,11 @@ public:
}
bool isLineSegment() const override {
- return size() == 2;
+ if constexpr (degree == 1) {
+ return true;
+ } else {
+ return BezierCurve::isLineSegment();
+ }
}
Curve *duplicate() const override {
@@ -286,6 +291,10 @@ public:
// call super. this is implemented only to allow specializations
BezierCurve::feed(sink, moveto_initial);
}
+ void expandToTransformed(Rect &bbox, Affine const &transform) const override {
+ // call super. this is implemented only to allow specializations
+ BezierCurve::expandToTransformed(bbox, transform);
+ }
};
// BezierCurveN<0> is meaningless; specialize it out
@@ -319,10 +328,15 @@ template <> inline bool BezierCurveN<1>::isLineSegment() const { return true; }
template <> Curve *BezierCurveN<1>::derivative() const;
template <> Coord BezierCurveN<1>::nearestTime(Point const &, Coord, Coord) const;
template <> std::vector<CurveIntersection> BezierCurveN<1>::intersect(Curve const &, Coord) const;
+template <> std::vector<CurveIntersection> BezierCurveN<2>::intersect(Curve const &, Coord) const;
+template <> std::vector<CurveIntersection> BezierCurveN<3>::intersect(Curve const &, Coord) const;
template <> int BezierCurveN<1>::winding(Point const &) const;
template <> void BezierCurveN<1>::feed(PathSink &sink, bool moveto_initial) const;
template <> void BezierCurveN<2>::feed(PathSink &sink, bool moveto_initial) const;
template <> void BezierCurveN<3>::feed(PathSink &sink, bool moveto_initial) const;
+template <> void BezierCurveN<1>::expandToTransformed(Rect &bbox, Affine const &transform) const;
+template <> void BezierCurveN<2>::expandToTransformed(Rect &bbox, Affine const &transform) const;
+template <> void BezierCurveN<3>::expandToTransformed(Rect &bbox, Affine const &transform) const;
inline Point middle_point(LineSegment const& _segment) {
return ( _segment.initialPoint() + _segment.finalPoint() ) / 2;
diff --git a/include/2geom/bezier.h b/include/2geom/bezier.h
index 07f9e9d..746a27b 100644
--- a/include/2geom/bezier.h
+++ b/include/2geom/bezier.h
@@ -157,6 +157,25 @@ public:
return *this;
}
+ bool operator==(Bezier const &other) const
+ {
+ if (degree() != other.degree()) {
+ return false;
+ }
+
+ for (size_t i = 0; i < c_.size(); i++) {
+ if (c_[i] != other.c_[i]) {
+ return false;
+ }
+ }
+ return true;
+ }
+
+ bool operator!=(Bezier const &other) const
+ {
+ return !(*this == other);
+ }
+
struct Order {
unsigned order;
explicit Order(Bezier const &b) : order(b.order()) {}
@@ -302,6 +321,14 @@ public:
}
Bezier &operator+=(Bezier const &other);
Bezier &operator-=(Bezier const &other);
+
+ /// Unary minus
+ Bezier operator-() const
+ {
+ Bezier result;
+ result.c_ = -c_;
+ return result;
+ }
};
@@ -339,6 +366,10 @@ OptInterval bounds_fast(Bezier const &b);
OptInterval bounds_exact(Bezier const &b);
OptInterval bounds_local(Bezier const &b, OptInterval const &i);
+/// Expand an interval to the image of a Bézier-Bernstein polynomial, assuming it already contains the initial point x0.
+void bezier_expand_to_image(Interval &range, Coord x0, Coord x1, Coord x2);
+void bezier_expand_to_image(Interval &range, Coord x0, Coord x1, Coord x2, Coord x3);
+
inline std::ostream &operator<< (std::ostream &os, const Bezier & b) {
os << "Bezier(";
for(unsigned i = 0; i < b.order(); i++) {
@@ -348,7 +379,8 @@ inline std::ostream &operator<< (std::ostream &os, const Bezier & b) {
return os;
}
-}
+} // namespace Geom
+
#endif // LIB2GEOM_SEEN_BEZIER_H
/*
diff --git a/include/2geom/conicsec.h b/include/2geom/conicsec.h
index 7ea2262..bfd5f36 100644
--- a/include/2geom/conicsec.h
+++ b/include/2geom/conicsec.h
@@ -411,6 +411,25 @@ public:
bool decompose (Line& l1, Line& l2) const;
+ /**
+ * @brief Division-free decomposition of a degenerate conic section, without
+ * degeneration test.
+ *
+ * When the conic is degenerate, it consists of 0, 1 or 2 lines in the xy-plane.
+ * This function returns these lines. But it does not check whether the conic
+ * is really degenerate, so calling it on a non-degenerate conic produces a
+ * meaningless result.
+ *
+ * If the number of lines is less than two, the trailing lines in the returned
+ * array will be degenerate. Use Line::isDegenerate() to test for that.
+ *
+ * This version of the decomposition is division-free, which improves numerical
+ * stability compared to decompose().
+ *
+ * @param epsilon The numerical threshold for floating point comparison of discriminants.
+ */
+ std::array<Line, 2> decompose_df(Coord epsilon = EPSILON) const;
+
/*
* Generate a RatQuad object from a conic arc.
*
diff --git a/include/2geom/coord.h b/include/2geom/coord.h
index 9cc220d..40db84e 100644
--- a/include/2geom/coord.h
+++ b/include/2geom/coord.h
@@ -49,14 +49,14 @@ enum Dim2 { X=0, Y=1 };
/** @brief Get the other (perpendicular) dimension.
* @ingroup Primitives */
-inline Dim2 other_dimension(Dim2 d) { return d == Y ? X : Y; }
+inline constexpr Dim2 other_dimension(Dim2 d) { return Dim2(int(d) ^ 1); }
// TODO: make a smarter implementation with C++11
template <typename T>
struct D2Traits {
- typedef typename T::D1Value D1Value;
- typedef typename T::D1Reference D1Reference;
- typedef typename T::D1ConstReference D1ConstReference;
+ using D1Value = typename T::D1Value;
+ using D1Reference = typename T::D1Reference;
+ using D1ConstReference = typename T::D1ConstReference;
};
/** @brief Axis extraction functor.
@@ -64,8 +64,8 @@ struct D2Traits {
* @ingroup Utilities */
template <Dim2 D, typename T>
struct GetAxis {
- typedef typename D2Traits<T>::D1Value result_type;
- typedef T argument_type;
+ using result_type = typename D2Traits<T>::D1Value;
+ using argument_type = T;
typename D2Traits<T>::D1Value operator()(T const &a) const {
return a[D];
}
@@ -73,28 +73,28 @@ struct GetAxis {
/** @brief Floating point type used to store coordinates.
* @ingroup Primitives */
-typedef double Coord;
+using Coord = double;
/** @brief Type used for integral coordinates.
* @ingroup Primitives */
-typedef int IntCoord;
+using IntCoord = int;
/** @brief Default "acceptably small" value.
* @ingroup Primitives */
-const Coord EPSILON = 1e-6; //1e-18;
+constexpr Coord EPSILON = 1e-6;
/** @brief Get a value representing infinity.
* @ingroup Primitives */
-inline Coord infinity() { return std::numeric_limits<Coord>::infinity(); }
+inline constexpr Coord infinity() { return std::numeric_limits<Coord>::infinity(); }
/** @brief Nearness predicate for values.
* @ingroup Primitives */
-inline bool are_near(Coord a, Coord b, double eps=EPSILON) { return a-b <= eps && a-b >= -eps; }
-inline bool rel_error_bound(Coord a, Coord b, double eps=EPSILON) { return a <= eps*b && a >= -eps*b; }
+inline constexpr bool are_near(Coord a, Coord b, double eps=EPSILON) { return std::abs(a-b) <= eps; }
+inline constexpr bool rel_error_bound(Coord a, Coord b, double eps=EPSILON) { return std::abs(a) <= eps*b; }
/** @brief Numerically stable linear interpolation.
* @ingroup Primitives */
-inline Coord lerp(Coord t, Coord a, Coord b) {
+inline constexpr Coord lerp(Coord t, Coord a, Coord b) {
return (1 - t) * a + t * b;
}
@@ -103,24 +103,22 @@ inline Coord lerp(Coord t, Coord a, Coord b) {
* @ingroup Utilities */
template <typename C>
struct CoordTraits {
- typedef D2<C> PointType;
- typedef GenericInterval<C> IntervalType;
- typedef GenericOptInterval<C> OptIntervalType;
- typedef GenericRect<C> RectType;
- typedef GenericOptRect<C> OptRectType;
+ using PointType = D2<C>;
+ using IntervalType = GenericInterval<C>;
+ using OptIntervalType = GenericOptInterval<C>;
+ using RectType = GenericRect<C>;
+ using OptRectType = GenericOptRect<C>;
- typedef
+ using IntervalOps =
boost::equality_comparable< IntervalType
, boost::orable< IntervalType
- > >
- IntervalOps;
+ >>;
- typedef
+ using RectOps =
boost::equality_comparable< RectType
, boost::orable< RectType
, boost::orable< RectType, OptRectType
- > > >
- RectOps;
+ >>>;
};
// NOTE: operator helpers for Rect and Interval are defined here.
@@ -128,56 +126,52 @@ struct CoordTraits {
template<>
struct CoordTraits<IntCoord> {
- typedef IntPoint PointType;
- typedef IntInterval IntervalType;
- typedef OptIntInterval OptIntervalType;
- typedef IntRect RectType;
- typedef OptIntRect OptRectType;
+ using PointType = IntPoint;
+ using IntervalType = IntInterval;
+ using OptIntervalType = OptIntInterval;
+ using RectType = IntRect;
+ using OptRectType = OptIntRect;
- typedef
+ using IntervalOps =
boost::equality_comparable< IntInterval
, boost::additive< IntInterval
, boost::additive< IntInterval, IntCoord
, boost::orable< IntInterval
- > > > >
- IntervalOps;
+ >>>>;
- typedef
+ using RectOps =
boost::equality_comparable< IntRect
, boost::orable< IntRect
, boost::orable< IntRect, OptIntRect
, boost::additive< IntRect, IntPoint
- > > > >
- RectOps;
+ >>>>;
};
template<>
struct CoordTraits<Coord> {
- typedef Point PointType;
- typedef Interval IntervalType;
- typedef OptInterval OptIntervalType;
- typedef Rect RectType;
- typedef OptRect OptRectType;
+ using PointType = Point;
+ using IntervalType = Interval;
+ using OptIntervalType = OptInterval;
+ using RectType = Rect;
+ using OptRectType = OptRect;
- typedef
+ using IntervalOps =
boost::equality_comparable< Interval
, boost::equality_comparable< Interval, IntInterval
, boost::additive< Interval
, boost::multipliable< Interval
, boost::orable< Interval
, boost::arithmetic< Interval, Coord
- > > > > > >
- IntervalOps;
+ >>>>>>;
- typedef
+ using RectOps =
boost::equality_comparable< Rect
, boost::equality_comparable< Rect, IntRect
, boost::orable< Rect
, boost::orable< Rect, OptRect
, boost::additive< Rect, Point
, boost::multipliable< Rect, Affine
- > > > > > >
- RectOps;
+ >>>>>>;
};
/** @brief Convert coordinate to shortest possible string.
@@ -198,7 +192,7 @@ std::string format_coord_nice(Coord x);
* @relates Coord */
Coord parse_coord(std::string const &s);
-} // end namespace Geom
+} // namespace Geom
#endif // LIB2GEOM_SEEN_COORD_H
diff --git a/include/2geom/curve.h b/include/2geom/curve.h
index 4470366..9519144 100644
--- a/include/2geom/curve.h
+++ b/include/2geom/curve.h
@@ -156,10 +156,16 @@ public:
virtual Rect boundsFast() const = 0;
/** @brief Compute the curve's exact bounding box.
- * This method can be dramatically slower than boundsExact() depending on the curve type.
+ * This method can be dramatically slower than boundsFast() depending on the curve type.
* @return The smallest possible rectangle containing all of the curve's points. */
virtual Rect boundsExact() const = 0;
+ /** @brief Expand the given rectangle to include the transformed curve,
+ * assuming it already contains its initial point.
+ * @param bbox[in,out] bbox The rectangle to expand; it is assumed to already contain (initialPoint() * transform).
+ * @param transform The transform to apply to the curve before taking its bounding box. */
+ virtual void expandToTransformed(Rect &bbox, Affine const &transform) const = 0;
+
// I have no idea what the 'deg' parameter is for, so this is undocumented for now.
virtual OptRect boundsLocal(OptInterval const &i, unsigned deg) const = 0;
diff --git a/include/2geom/ellipse.h b/include/2geom/ellipse.h
index f6089c9..0d1567a 100644
--- a/include/2geom/ellipse.h
+++ b/include/2geom/ellipse.h
@@ -174,6 +174,13 @@ public:
/// Get the tight-fitting bounding box of the ellipse.
Rect boundsExact() const;
+ /** @brief Get a fast to compute bounding box which contains the ellipse.
+ *
+ * The returned rectangle engulfs the ellipse but it may not be the smallest
+ * axis-aligned rectangle with this property.
+ */
+ Rect boundsFast() const;
+
/// Get the coefficients of the ellipse's implicit equation.
std::vector<double> coefficients() const;
void coefficients(Coord &A, Coord &B, Coord &C, Coord &D, Coord &E, Coord &F) const;
diff --git a/include/2geom/elliptical-arc.h b/include/2geom/elliptical-arc.h
index b95cfb7..567e207 100644
--- a/include/2geom/elliptical-arc.h
+++ b/include/2geom/elliptical-arc.h
@@ -270,6 +270,7 @@ public:
return boundsExact();
}
Rect boundsExact() const override;
+ void expandToTransformed(Rect &bbox, Affine const &transform) const override;
// TODO: native implementation of the following methods
OptRect boundsLocal(OptInterval const &i, unsigned int deg) const override {
return SBasisCurve(toSBasis()).boundsLocal(i, deg);
@@ -308,6 +309,7 @@ private:
void _updateCenterAndAngles();
std::vector<ShapeIntersection> _filterIntersections(std::vector<ShapeIntersection> &&xs, bool is_first) const;
bool _validateIntersection(ShapeIntersection &xing, bool is_first) const;
+ std::vector<ShapeIntersection> _intersectSameEllipse(EllipticalArc const *other) const;
Point _initial_point, _final_point;
Ellipse _ellipse;
diff --git a/include/2geom/exception.h b/include/2geom/exception.h
index 88ba517..b472aae 100644
--- a/include/2geom/exception.h
+++ b/include/2geom/exception.h
@@ -116,6 +116,18 @@ public:
};
#define THROW_INFINITESOLUTIONS(i) throw(InfiniteSolutions(__FILE__, __LINE__))
+class InfinitelyManySolutions : public RangeError {
+private:
+ char const *const _message;
+public:
+ InfinitelyManySolutions(const char *file, const int line, char const *message)
+ : RangeError("There are infinitely many solutions", file, line)
+ , _message{message}
+ {}
+ char const *what() const noexcept override { return _message; }
+};
+#define THROW_INFINITELY_MANY_SOLUTIONS(msg) throw(InfinitelyManySolutions(__FILE__, __LINE__, msg))
+
class ContinuityError : public RangeError {
public:
ContinuityError(const char *file, const int line)
diff --git a/include/2geom/generic-rect.h b/include/2geom/generic-rect.h
index be10019..4524d43 100644
--- a/include/2geom/generic-rect.h
+++ b/include/2geom/generic-rect.h
@@ -192,6 +192,12 @@ public:
/** @brief Get the smaller extent (width or height) of the rectangle. */
C minExtent() const { return std::min(f[X].extent(), f[Y].extent()); }
+ /** @brief Get rectangle's distance SQUARED away from the given point **/
+ C distanceSq(const CPoint pt) const {
+ auto v = clamp(pt) - pt;
+ return v.x() * v.x() + v.y() * v.y();
+ }
+
/** @brief Clamp point to the rectangle. */
CPoint clamp(CPoint const &p) const {
CPoint result(f[X].clamp(p[X]), f[Y].clamp(p[Y]));
@@ -416,6 +422,11 @@ public:
bool contains(CPoint const &p) const { return *this && (*this)->contains(p); }
/// @}
+ /** @brief Returns an empty optional (testing false) if the rectangle has zero area. */
+ OptCRect regularized() const {
+ return *this && !(*this)->hasZeroArea() ? *this : OptCRect();
+ }
+
/// @name Modify the potentially empty rectangle.
/// @{
/** @brief Enlarge the rectangle to contain the argument.
diff --git a/include/2geom/int-point.h b/include/2geom/int-point.h
index 2a951a3..6dbed11 100644
--- a/include/2geom/int-point.h
+++ b/include/2geom/int-point.h
@@ -49,17 +49,20 @@ namespace Geom {
class IntPoint
: boost::additive< IntPoint
, boost::totally_ordered< IntPoint
- > >
+ , boost::multiplicative< IntPoint, IntCoord
+ , boost::multiplicative< IntPoint
+ > > > >
{
- IntCoord _pt[2];
+ IntCoord _pt[2] = { 0, 0 };
public:
- /// @name Creating integer points
+ /// @name Create integer points
/// @{
- IntPoint() { }
- IntPoint(IntCoord x, IntCoord y) {
- _pt[X] = x;
- _pt[Y] = y;
- }
+ /** Construct a point at the origin. */
+ IntPoint() = default;
+ /** Construct a point from its coordinates. */
+ IntPoint(IntCoord x, IntCoord y)
+ : _pt{ x, y }
+ {}
/// @}
/// @name Access the coordinates of a point
@@ -84,8 +87,7 @@ public:
/// @name Vector-like arithmetic operations
/// @{
IntPoint operator-() const {
- IntPoint ret(-_pt[X], -_pt[Y]);
- return ret;
+ return IntPoint(-_pt[X], -_pt[Y]);
}
IntPoint &operator+=(IntPoint const &o) {
_pt[X] += o._pt[X];
@@ -97,6 +99,26 @@ public:
_pt[Y] -= o._pt[Y];
return *this;
}
+ IntPoint &operator*=(IntPoint const &o) {
+ _pt[X] *= o._pt[X];
+ _pt[Y] *= o._pt[Y];
+ return *this;
+ }
+ IntPoint &operator*=(IntCoord o) {
+ _pt[X] *= o;
+ _pt[Y] *= o;
+ return *this;
+ }
+ IntPoint &operator/=(IntPoint const &o) {
+ _pt[X] /= o._pt[X];
+ _pt[Y] /= o._pt[Y];
+ return *this;
+ }
+ IntPoint &operator/=(IntCoord o) {
+ _pt[X] /= o;
+ _pt[Y] /= o;
+ return *this;
+ }
/// @}
/// @name Various utilities
diff --git a/include/2geom/intersection-graph.h b/include/2geom/intersection-graph.h
index 5b187c1..940c43c 100644
--- a/include/2geom/intersection-graph.h
+++ b/include/2geom/intersection-graph.h
@@ -43,26 +43,114 @@
namespace Geom {
+/** @class PathIntersectionGraph
+ * @brief Intermediate data for computing Boolean operations on paths.
+ *
+ * This class implements the Greiner-Hormann clipping algorithm,
+ * with improvements inspired by Foster and Overfelt as well as some
+ * original contributions.
+ *
+ * For the purposes of boolean operations, a shape is defined as a PathVector
+ * using the "even-odd" rule, i.e., regions with odd winding are considered part
+ * of the shape, whereas regions with even winding are not.
+ *
+ * For this reason, the two path-vectors are sometimes called "shapes" or "operands" of
+ * the boolean operation. Each path-vector may contain several paths, which are called
+ * either "paths" or "components" in the documentation.
+ *
+ * @ingroup Paths
+ */
class PathIntersectionGraph
{
// this is called PathIntersectionGraph so that we can also have a class for polygons,
// e.g. PolygonIntersectionGraph, which is going to be significantly faster
public:
+ /** @brief Construct a path intersection graph for two shapes described via their boundaries.
+ * The boundaries are passed as path-vectors.
+ *
+ * @param a – the first operand, also referred to as operand A.
+ * @param b – the second operand, also referred to as operand B.
+ * @param precision – precision setting used for intersection calculations.
+ */
PathIntersectionGraph(PathVector const &a, PathVector const &b, Coord precision = EPSILON);
+ /**
+ * @brief Get the union of the shapes, A ∪ B.
+ *
+ * A point belongs to the union if and only if it belongs to at least one of the operands.
+ *
+ * @return A path-vector describing the union of the operands A and B.
+ */
PathVector getUnion();
+
+ /**
+ * @brief Get the intersection of the shapes, A ∩ B.
+ *
+ * A point belongs to the intersection if and only if it belongs to both shapes.
+ *
+ * @return A path-vector describing the intersection of the operands A and B.
+ */
PathVector getIntersection();
+
+ /**
+ * @brief Get the difference of the shapes, A ∖ B.
+ *
+ * A point belongs to the difference if and only if it belongs to A but not to B.
+ *
+ * @return A path-vector describing the difference of the operands A and B.
+ */
PathVector getAminusB();
+
+ /**
+ * @brief Get the opposite difference of the shapes, B ∖ A.
+ *
+ * A point belongs to the difference if and only if it belongs to B but not to A.
+ *
+ * @return A path-vector describing the difference of the operands B and A.
+ */
PathVector getBminusA();
+
+ /**
+ * @brief Get the symmetric difference of the shapes, A ∆ B.
+ *
+ * A point belongs to the symmetric difference if and only if it belongs to one of the two
+ * shapes A or B, but not both. This is equivalent to the logical XOR operation: the elements
+ * of A ∆ B are points which are in A XOR in B.
+ *
+ * @return A path-vector describing the symmetric difference of the operands A and B.
+ */
PathVector getXOR();
/// Returns the number of intersections used when computing Boolean operations.
std::size_t size() const;
+
+ /**
+ * @brief Get the geometric points where the two path-vectors intersect.
+ *
+ * Degenerate intersection points, where the shapes merely "kiss", are not retured.
+ *
+ * @param defective – whether to return only the defective crossings or only the true crossings.
+ * @return If defective is true, returns a vector containing all defective intersection points,
+ * i.e., points that are neither true transverse intersections nor degenerate intersections.
+ * If defective is false, returns all true transverse intersections.
+ */
std::vector<Point> intersectionPoints(bool defective = false) const;
+
+ /**
+ * @brief Get the geometric points located on path portions between consecutive intersections.
+ *
+ * These points were used for the winding number calculations which determined which path portions
+ * lie inside the other shape and which lie outside.
+ *
+ * @return A vector containing all sample points used for winding calculations.
+ */
std::vector<Point> windingPoints() const {
return _winding_points;
}
+
void fragments(PathVector &in, PathVector &out) const;
+
+
bool valid() const { return _graph_valid; }
private:
@@ -75,11 +163,18 @@ private:
struct IntersectionVertex {
boost::intrusive::list_member_hook<> _hook;
boost::intrusive::list_member_hook<> _proc_hook;
- PathVectorTime pos;
- Point p; // guarantees that endpoints are exact
- IntersectionVertex *neighbor;
+ PathVectorTime pos; ///< Intersection time.
+ Point p; ///< Geometric position of the intersection point; guarantees that endpoints are exact.
+ IntersectionVertex *neighbor; ///< A pointer to the corresponding vertex on the other shape.
+ /** Tells us whether the edge originating at this intersection lies inside or outside of
+ * the shape given by the other path-vector. The "edge originating" at this intersection is
+ * the portion of the path between this intersection and the next intersection, in the
+ * direction of increasing path time. */
InOutFlag next_edge;
- unsigned which;
+ unsigned which; ///< Index of the operand path-vector that this intersection vertex lies on.
+ /** Whether the intersection is defective, which means that for some reason the paths
+ * neither cross transversally through each other nor "kiss" at a common tangency point.
+ */
bool defective;
};
@@ -101,10 +196,14 @@ private:
>
> UnprocessedList;
+ /// Stores processed intersection information for a single path in an operand path-vector.
struct PathData {
- IntersectionList xlist;
- std::size_t path_index;
- int which;
+ IntersectionList xlist; ///< List of crossings on this particular path.
+ std::size_t path_index; ///< Index of the path in its path-vector.
+ int which; ///< Index of the path-vector (in PathIntersectionGraph::_pv) that the path belongs to.
+ /** Whether this path as a whole is contained INSIDE or OUTSIDE relative to the other path-vector.
+ * The value BOTH means that some portions of the path are inside while others are outside.
+ */
InOutFlag status;
PathData(int w, std::size_t pi)
@@ -130,11 +229,14 @@ private:
ILIter _getNeighbor(ILIter iter);
PathData &_getPathData(ILIter iter);
- PathVector _pv[2];
- boost::ptr_vector<IntersectionVertex> _xs;
- boost::ptr_vector<PathData> _components[2];
- UnprocessedList _ulist;
- bool _graph_valid;
+ PathVector _pv[2]; ///< Stores the two operand path-vectors, A at _pv[0] and B at _pv[1].
+ boost::ptr_vector<IntersectionVertex> _xs; ///< Stores all crossings between the two shapes.
+ boost::ptr_vector<PathData> _components[2]; ///< Stores the crossing information for the operands.
+ UnprocessedList _ulist; ///< Temporarily holds all unprocessed during a boolean operation.
+ bool _graph_valid; ///< Whether all intersections are regular.
+ /** Stores sample points located on paths of the operand path-vectors,
+ * between consecutive intersections.
+ */
std::vector<Point> _winding_points;
friend std::ostream &operator<<(std::ostream &, PathIntersectionGraph const &);
diff --git a/include/2geom/interval.h b/include/2geom/interval.h
index 3401d6a..6731e5b 100644
--- a/include/2geom/interval.h
+++ b/include/2geom/interval.h
@@ -98,16 +98,16 @@ public:
}
/** @brief Map the interval [0,1] onto this one.
* This method simply performs 1D linear interpolation between endpoints. */
- Coord valueAt(Coord t) {
+ Coord valueAt(Coord t) const {
return lerp(t, min(), max());
}
/** @brief Compute a time value that maps to the given value.
* The supplied value does not need to be in the interval for this method to work. */
- Coord timeAt(Coord v) {
+ Coord timeAt(Coord v) const {
return (v - min()) / extent();
}
/// Find closest time in [0,1] that maps to the given value. */
- Coord nearestTime(Coord v) {
+ Coord nearestTime(Coord v) const {
if (v <= min()) return 0;
if (v >= max()) return 1;
return timeAt(v);
diff --git a/include/2geom/line.h b/include/2geom/line.h
index e05660e..9a56602 100644
--- a/include/2geom/line.h
+++ b/include/2geom/line.h
@@ -127,7 +127,8 @@ public:
/// Get the line's origin point.
Point origin() const { return _initial; }
/** @brief Get the line's raw direction vector.
- * The retrieved vector is normalized to unit length. */
+ * The length of the retrieved vector is equal to the length of a segment parametrized by
+ * a time interval of length 1. */
Point vector() const { return _final - _initial; }
/** @brief Get the line's normalized direction vector.
* The retrieved vector is normalized to unit length. */
diff --git a/include/2geom/math-utils.h b/include/2geom/math-utils.h
index ee923f5..4c35a80 100644
--- a/include/2geom/math-utils.h
+++ b/include/2geom/math-utils.h
@@ -37,6 +37,7 @@
#define LIB2GEOM_SEEN_MATH_UTILS_H
#include <math.h> // sincos is usually only available in math.h
+#include <array>
#include <cmath>
#include <utility> // for std::pair
#include <boost/math/special_functions/fpclassify.hpp>
@@ -66,7 +67,7 @@ template <class T> inline const T& between (const T& min, const T& max, const T&
Implemented in terms of round, i.e. we make no guarantees as to what happens if x is
half way between two rounded numbers.
-
+
Note: places is the number of decimal places without using scientific (e) notation, not the
number of significant figures. This function may not be suitable for values of x whose
magnitude is so far from 1 that one would want to use scientific (e) notation.
@@ -94,6 +95,35 @@ inline void sincos(double angle, double &sin_, double &cos_) {
#endif
}
+/** @brief Scale the doubles in the passed array to make them "reasonably large".
+ *
+ * All doubles in the passed array will get scaled by the same power of 2 (which is
+ * a lossless operation) in such a way that their geometric average gets closer to 1.
+ *
+ * @tparam N The size of the passed array.
+ * @param[in,out] values The doubles to be rescaled in place.
+ * @return The exponent in the power of two by which the doubles got scaled.
+ */
+template <size_t N>
+inline int rescale_homogenous(std::array<double, N> &values)
+{
+ if constexpr (N == 0) {
+ return 0;
+ }
+ std::array<int, N> exponents;
+ std::array<double, N> mantissas;
+ int average = 0;
+ for (size_t i = 0; i < N; i++) {
+ mantissas[i] = std::frexp(values[i], &exponents[i]);
+ average += exponents[i];
+ }
+ average /= (int)N;
+ for (size_t i = 0; i < N; i++) {
+ values[i] = std::ldexp(mantissas[i], exponents[i] - average);
+ }
+ return -average;
+}
+
} // end namespace Geom
#endif // LIB2GEOM_SEEN_MATH_UTILS_H
diff --git a/include/2geom/path.h b/include/2geom/path.h
index 9a71ee0..f3042d7 100644
--- a/include/2geom/path.h
+++ b/include/2geom/path.h
@@ -35,14 +35,18 @@
#ifndef LIB2GEOM_SEEN_PATH_H
#define LIB2GEOM_SEEN_PATH_H
+#include <cstddef>
#include <iterator>
#include <algorithm>
#include <iostream>
#include <memory>
#include <optional>
#include <utility>
+#include <vector>
+
#include <boost/operators.hpp>
#include <boost/ptr_container/ptr_vector.hpp>
+
#include <2geom/intersection.h>
#include <2geom/curve.h>
#include <2geom/bezier-curve.h>
@@ -106,6 +110,10 @@ class BaseIterator
index -= d;
return *this;
}
+ std::ptrdiff_t operator-(Self const &other) const {
+ assert(path == other.path);
+ return (std::ptrdiff_t)index - (std::ptrdiff_t)other.index;
+ }
private:
P *path;
@@ -124,7 +132,7 @@ class BaseIterator
* in long paths, either use this class and related methods instead of the standard methods
* pointAt(), nearestTime() and so on, or use curveAt() to first obtain the curve, then
* call the method again to obtain a high precision result.
- *
+ *
* @ingroup Paths */
struct PathTime
: boost::totally_ordered<PathTime>
@@ -277,6 +285,26 @@ struct ShapeTraits<Path> {
typedef PathIntersection IntersectionType;
};
+/** @brief Stores information about the extremum points on a path, with respect
+ * to one of the coordinate axes.
+ * @relates Path::extrema()
+ */
+struct PathExtrema {
+ /** Points with the minimum and maximum value of a coordinate. */
+ Point min_point, max_point;
+
+ /** Directions in which the OTHER coordinates change at the extremum points.
+ *
+ * - equals +1.0 if the other coordinate increases,
+ * - equals 0.0 if the other coordinate is constant (e.g., for an axis-aligned segment),
+ * - equals -1.0 if the other coordinate decreases.
+ */
+ float glance_direction_at_min, glance_direction_at_max;
+
+ /** Path times corresponding to minimum and maximum points. */
+ PathTime min_time, max_time;
+};
+
/** @brief Sequence of contiguous curves, aka spline.
*
* Path represents a sequence of contiguous curves, also known as a spline.
@@ -536,7 +564,7 @@ public:
* Note that this method has reduced precision with respect to calling pointAt()
* directly on the curve. If you want high precision results, use the version
* that takes a PathTime parameter.
- *
+ *
* Allowed time values range from zero to the number of curves; you can retrieve
* the allowed range of values with timeRange(). */
Point pointAt(Coord t) const;
@@ -553,14 +581,24 @@ public:
Point operator()(Coord t) const { return pointAt(t); }
+ /** @brief Find the extrema of the specified coordinate.
+ *
+ * Returns a PathExtrema struct describing "witness" points on the path
+ * where the specified coordinate attains its minimum and maximum values.
+ */
+ PathExtrema extrema(Dim2 d) const;
+
/// Compute intersections with axis-aligned line.
std::vector<PathTime> roots(Coord v, Dim2 d) const;
/// Compute intersections with another path.
std::vector<PathIntersection> intersect(Path const &other, Coord precision = EPSILON) const;
+ /// Compute intersections of the path with itself.
+ std::vector<PathIntersection> intersectSelf(Coord precision = EPSILON) const;
+
/** @brief Determine the winding number at the specified point.
- *
+ *
* The winding number is the number of full turns made by a ray that connects the passed
* point and the path's value (i.e. the result of the pointAt() method) as the time increases
* from 0 to the maximum valid value. Positive numbers indicate turns in the direction
@@ -643,7 +681,7 @@ public:
Path reversed() const;
void insert(iterator pos, Curve const &curve);
-
+
template <typename Iter>
void insert(iterator pos, Iter first, Iter last) {
_unshare();
@@ -670,6 +708,29 @@ public:
* If the path is closed, this is always the same as the initial point. */
Point finalPoint() const { return (*_closing_seg)[_closed ? 1 : 0]; }
+ /** @brief Get the unit tangent vector at the start of the path,
+ * or the zero vector if undefined. */
+ Point initialUnitTangent() const {
+ for (auto const &curve : *this) {
+ if (!curve.isDegenerate()) {
+ return curve.unitTangentAt(0.0);
+ }
+ }
+ return Point();
+ }
+
+ /** @brief Get the unit tangent vector at the end of the path,
+ * or the zero vector if undefined. */
+ Point finalUnitTangent() const {
+ for (auto it = end(); it != begin();) {
+ --it;
+ if (!it->isDegenerate()) {
+ return it->unitTangentAt(1.0);
+ }
+ }
+ return Point();
+ }
+
void setInitial(Point const &p) {
_unshare();
_closed = false;
@@ -752,6 +813,10 @@ public:
/// Append a stitching segment ending at the specified point.
void stitchTo(Point const &p);
+ /** @brief Return a copy of the path without degenerate curves, except possibly for a
+ * degenerate closing segment. */
+ Path withoutDegenerateCurves() const;
+
/** @brief Verify the continuity invariant.
* If the path is not contiguous, this will throw a CountinuityError. */
void checkContinuity() const;
@@ -807,6 +872,32 @@ inline Coord nearest_time(Point const &p, Path const &c) {
bool are_near(Path const &a, Path const &b, Coord precision = EPSILON);
+/**
+ * @brief Find the first point where two paths diverge away from one another.
+ *
+ * If the two paths have a common starting point, the algorithm follows them for as long as the
+ * images of the paths coincide and finds the first point where they stop coinciding. Note that
+ * only the images of paths in the plane are compared, and not their parametrizations, so this
+ * is not a functional (parametric) coincidence. If you want to test parametric coincidence, use
+ * bool are_near(Path const&, Path const&, Coord) instead.
+ *
+ * The function returns the point where the traces of the two paths finally diverge up to the
+ * specified precision. If the traces (images) of the paths are nearly identical until the end,
+ * the returned point is their (almost) common endpoint. If however the image of one of the paths
+ * is completely contained in the image of the other path, the returned point is the endpoint of
+ * the shorter path.
+ *
+ * If the paths have different starting points, then the returned intersection has the special
+ * time values of -1.0 on both paths and the returned intersection point is the midpoint of the
+ * line segment connecting the two starting points.
+ *
+ * @param first The first path to follow; corresponds to .first in the return value.
+ * @param second The second path to follow; corresponds to .second in the return value.
+ * @param precision How close the paths' images need to be in order to be considered as overlapping.
+ * @return A path intersection specifying the point and path times where the two paths part ways.
+ */
+PathIntersection parting_point(Path const &first, Path const &second, Coord precision = EPSILON);
+
std::ostream &operator<<(std::ostream &out, Path const &path);
} // end namespace Geom
diff --git a/include/2geom/pathvector.h b/include/2geom/pathvector.h
index 5781663..3fd7d36 100644
--- a/include/2geom/pathvector.h
+++ b/include/2geom/pathvector.h
@@ -222,6 +222,9 @@ public:
Point finalPoint() const {
return _data.back().finalPoint();
}
+ /** @brief Get all intersections of the path-vector with itself. This includes both
+ * self-intersections of constituent paths and intersections between different paths. */
+ std::vector<PathVectorIntersection> intersectSelf(Coord precision = EPSILON) const;
Path &pathAt(Coord t, Coord *rest = NULL);
Path const &pathAt(Coord t, Coord *rest = NULL) const;
Curve const &curveAt(Coord t, Coord *rest = NULL) const;
diff --git a/include/2geom/point.h b/include/2geom/point.h
index 1e0b1eb..3a29066 100644
--- a/include/2geom/point.h
+++ b/include/2geom/point.h
@@ -50,6 +50,8 @@ class Point
: boost::additive< Point
, boost::totally_ordered< Point
, boost::multiplicative< Point, Coord
+ , boost::multiplicative< Point
+ , boost::multiplicative< Point, IntPoint
, MultipliableNoncommutative< Point, Affine
, MultipliableNoncommutative< Point, Translate
, MultipliableNoncommutative< Point, Rotate
@@ -57,29 +59,26 @@ class Point
, MultipliableNoncommutative< Point, HShear
, MultipliableNoncommutative< Point, VShear
, MultipliableNoncommutative< Point, Zoom
- > > > > > > > > > > // base class chaining, see documentation for Boost.Operator
+ > > > > > > > > > > > > // base class chaining, see documentation for Boost.Operator
{
- Coord _pt[2];
+ Coord _pt[2] = { 0, 0 };
public:
- typedef Coord D1Value;
- typedef Coord &D1Reference;
- typedef Coord const &D1ConstReference;
+ using D1Value = Coord;
+ using D1Reference = Coord &;
+ using D1ConstReference = Coord const &;
/// @name Create points
/// @{
- /** Construct a point on the origin. */
- Point()
- { _pt[X] = _pt[Y] = 0; }
-
+ /** Construct a point at the origin. */
+ Point() = default;
/** Construct a point from its coordinates. */
- Point(Coord x, Coord y) {
- _pt[X] = x; _pt[Y] = y;
- }
+ Point(Coord x, Coord y)
+ : _pt{ x, y }
+ {}
/** Construct from integer point. */
- Point(IntPoint const &p) {
- _pt[X] = p[X];
- _pt[Y] = p[Y];
- }
+ Point(IntPoint const &p)
+ : Point(p[X], p[Y])
+ {}
/** @brief Construct a point from its polar coordinates.
* The angle is specified in radians, in the mathematical convention (increasing
* counter-clockwise from +X). */
@@ -112,7 +111,7 @@ public:
/// @{
/** @brief Compute the distance from origin.
* @return Length of the vector from origin to this point */
- Coord length() const { return hypot(_pt[0], _pt[1]); }
+ Coord length() const { return std::hypot(_pt[0], _pt[1]); }
void normalize();
Point normalized() const {
Point ret(*this);
@@ -141,26 +140,44 @@ public:
return Point(-_pt[X], -_pt[Y]);
}
Point &operator+=(Point const &o) {
- for ( unsigned i = 0 ; i < 2 ; ++i ) {
- _pt[i] += o._pt[i];
- }
+ _pt[X] += o._pt[X];
+ _pt[Y] += o._pt[Y];
return *this;
}
Point &operator-=(Point const &o) {
- for ( unsigned i = 0 ; i < 2 ; ++i ) {
- _pt[i] -= o._pt[i];
- }
+ _pt[X] -= o._pt[X];
+ _pt[Y] -= o._pt[Y];
return *this;
}
Point &operator*=(Coord s) {
for (double & i : _pt) i *= s;
return *this;
}
+ Point &operator*=(Point const &o) {
+ _pt[X] *= o._pt[X];
+ _pt[Y] *= o._pt[Y];
+ return *this;
+ }
+ Point &operator*=(IntPoint const &o) {
+ _pt[X] *= o.x();
+ _pt[Y] *= o.y();
+ return *this;
+ }
Point &operator/=(Coord s) {
//TODO: s == 0?
for (double & i : _pt) i /= s;
return *this;
}
+ Point &operator/=(Point const &o) {
+ _pt[X] /= o._pt[X];
+ _pt[Y] /= o._pt[Y];
+ return *this;
+ }
+ Point &operator/=(IntPoint const &o) {
+ _pt[X] /= o.x();
+ _pt[Y] /= o.y();
+ return *this;
+ }
/// @}
/// @name Affine transformations
diff --git a/include/2geom/sbasis-curve.h b/include/2geom/sbasis-curve.h
index 4c1bd5b..93d6772 100644
--- a/include/2geom/sbasis-curve.h
+++ b/include/2geom/sbasis-curve.h
@@ -103,6 +103,9 @@ public:
}
Rect boundsFast() const override { return *bounds_fast(inner); }
Rect boundsExact() const override { return *bounds_exact(inner); }
+ void expandToTransformed(Rect &bbox, Affine const &transform) const override {
+ bbox |= bounds_exact(inner * transform);
+ }
OptRect boundsLocal(OptInterval const &i, unsigned deg) const override {
return bounds_local(inner, i, deg);
}
diff --git a/include/2geom/utils.h b/include/2geom/utils.h
index 2efcddd..87f49bb 100644
--- a/include/2geom/utils.h
+++ b/include/2geom/utils.h
@@ -45,9 +45,6 @@ enum Errors : ErrorCode {
GEOM_ERR_INTERSECGRAPH,
};
-// proper logical xor
-inline bool logical_xor (bool a, bool b) { return (a || b) && !(a && b); }
-
void binomial_coefficients(std::vector<size_t>& bc, std::size_t n);
struct EmptyClass {};
diff --git a/src/2geom/CMakeLists.txt b/src/2geom/CMakeLists.txt
index df03fc1..734ee96 100755
--- a/src/2geom/CMakeLists.txt
+++ b/src/2geom/CMakeLists.txt
@@ -55,6 +55,8 @@ add_library(2geom ${LIB_TYPE}
numeric/matrix.cpp
parallelogram.cpp
+ parting-point.cpp
+ path-extrema.cpp
path-intersection.cpp
path-sink.cpp
path.cpp
@@ -73,6 +75,7 @@ add_library(2geom ${LIB_TYPE}
sbasis-roots.cpp
sbasis-to-bezier.cpp
sbasis.cpp
+ self-intersect.cpp
solve-bezier.cpp
solve-bezier-one-d.cpp
solve-bezier-parametric.cpp
@@ -85,6 +88,10 @@ add_library(2geom ${LIB_TYPE}
utils.cpp
# headers (for IDE support only)
+ # private:
+ planar-graph.h
+
+ # public:
${2GEOM_INCLUDE_DIR}/2geom/affine.h
${2GEOM_INCLUDE_DIR}/2geom/angle.h
diff --git a/src/2geom/bezier-curve.cpp b/src/2geom/bezier-curve.cpp
index 2ea0dbe..ca0f787 100644
--- a/src/2geom/bezier-curve.cpp
+++ b/src/2geom/bezier-curve.cpp
@@ -35,6 +35,7 @@
#include <2geom/path-sink.h>
#include <2geom/basic-intersection.h>
#include <2geom/nearest-time.h>
+#include <2geom/polynomial.h>
namespace Geom
{
@@ -126,6 +127,29 @@ bool BezierCurve::isDegenerate() const
return true;
}
+/** Return false if there are at least 3 distinct control points, true otherwise. */
+bool BezierCurve::isLineSegment() const
+{
+ auto const last_idx = size() - 1;
+ if (last_idx == 1) {
+ return true;
+ }
+ auto const start = controlPoint(0);
+ auto const end = controlPoint(last_idx);
+ for (unsigned i = 1; i < last_idx; ++i) {
+ auto const pi = controlPoint(i);
+ if (pi != start && pi != end) {
+ return false;
+ }
+ }
+ return true;
+}
+
+void BezierCurve::expandToTransformed(Rect &bbox, Affine const &transform) const
+{
+ bbox |= bounds_exact(inner * transform);
+}
+
Coord BezierCurve::length(Coord tolerance) const
{
switch (order())
@@ -198,9 +222,40 @@ bool BezierCurve::isNear(Curve const &c, Coord precision) const
}
return true;
} else {
- // TODO: comparison after degree elevation
- return false;
+ // Must equalize the degrees before comparing
+ BezierCurve elevated_this, elevated_other;
+ for (size_t dim = 0; dim < 2; dim++) {
+ unsigned const our_degree = inner[dim].degree();
+ unsigned const other_degree = other->inner[dim].degree();
+
+ if (our_degree < other_degree) {
+ // Elevate our degree
+ elevated_this.inner[dim] = inner[dim].elevate_to_degree(other_degree);
+ elevated_other.inner[dim] = other->inner[dim];
+ } else if (our_degree > other_degree) {
+ // Elevate the other's degree
+ elevated_this.inner[dim] = inner[dim];
+ elevated_other.inner[dim] = other->inner[dim].elevate_to_degree(our_degree);
+ } else {
+ // Equal degrees: just copy
+ elevated_this.inner[dim] = inner[dim];
+ elevated_other.inner[dim] = other->inner[dim];
+ }
+ }
+ assert(elevated_other.size() == elevated_this.size());
+ return elevated_this.isNear(elevated_other, precision);
+ }
+}
+
+Curve *BezierCurve::portion(Coord f, Coord t) const
+{
+ if (f == 0.0 && t == 1.0) {
+ return duplicate();
}
+ if (f == 1.0 && t == 0.0) {
+ return reverse();
+ }
+ return new BezierCurve(Geom::portion(inner, f, t));
}
bool BezierCurve::operator==(Curve const &c) const
@@ -294,6 +349,7 @@ Coord BezierCurveN<1>::nearestTime(Point const& p, Coord from, Coord to) const
else return from + t*(to-from);
}
+/* Specialized intersection routine for line segments. */
template <>
std::vector<CurveIntersection> BezierCurveN<1>::intersect(Curve const &other, Coord eps) const
{
@@ -314,6 +370,93 @@ std::vector<CurveIntersection> BezierCurveN<1>::intersect(Curve const &other, Co
return result;
}
+/** @brief Find intersections of a low-degree Bézier curve with a line segment.
+ *
+ * Uses algebraic solutions to low-degree polynomial equations which may be faster
+ * and more precise than iterative methods.
+ *
+ * @tparam degree The degree of the Bézier curve; must be 2 or 3.
+ * @param curve A Bézier curve of the given degree.
+ * @param line A line (but really a segment).
+ * @return Intersections between the passed curve and the fundamental segment of the line
+ * (the segment where the time parameter lies in the unit interval).
+ */
+template <unsigned degree>
+static std::vector<CurveIntersection> bezier_line_intersections(BezierCurveN<degree> const &curve, Line const &line)
+{
+ static_assert(degree == 2 || degree == 3, "bezier_line_intersections<degree>() error: degree must be 2 or 3.");
+
+ auto const length = distance(line.initialPoint(), line.finalPoint());
+ if (length == 0) {
+ return {};
+ }
+ std::vector<CurveIntersection> result;
+
+ // Find the isometry mapping the line to the x-axis, taking the initial point to the origin
+ // and the final point to (length, 0). Apply this transformation to the Bézier curve and
+ // extract the y-coordinate polynomial.
+ auto const transform = line.rotationToZero(Y);
+ Bezier const y = (curve.fragment() * transform)[Y];
+ std::vector<double> roots;
+
+ // Find roots of the polynomial y.
+ {
+ double const c2 = y[0] + y[2] - 2.0 * y[1];
+ double const c1 = y[1] - y[0];
+ double const c0 = y[0];
+
+ if constexpr (degree == 2) {
+ roots = solve_quadratic(c2, 2.0 * c1, c0);
+ } else if constexpr (degree == 3) {
+ double const c3 = y[3] - y[0] + 3.0 * (y[1] - y[2]);
+ roots = solve_cubic(c3, 3.0 * c2, 3.0 * c1 , c0);
+ }
+ }
+
+ // Filter the roots and assemble intersections.
+ for (double root : roots) {
+ if (root < 0.0 || root > 1.0) {
+ continue;
+ }
+ Coord x = (curve.pointAt(root) * transform)[X];
+ if (x < 0.0 || x > length) {
+ continue;
+ }
+ result.emplace_back(curve, line, root, x / length);
+ }
+ return result;
+}
+
+/* Specialized intersection routine for quadratic Bézier curves. */
+template <>
+std::vector<CurveIntersection> BezierCurveN<2>::intersect(Curve const &other, Coord eps) const
+{
+ if (auto other_bezier = dynamic_cast<BezierCurve const *>(&other)) {
+ auto const other_degree = other_bezier->order();
+ if (other_degree == 1) {
+ // Use the exact method to intersect a quadratic Bézier with a line segment.
+ auto line = Line(other_bezier->initialPoint(), other_bezier->finalPoint());
+ return bezier_line_intersections<2>(*this, line);
+ }
+ // TODO: implement exact intersection of two quadratic Béziers using the method of resultants.
+ }
+ return BezierCurve::intersect(other, eps);
+}
+
+/* Specialized intersection routine for cubic Bézier curves. */
+template <>
+std::vector<CurveIntersection> BezierCurveN<3>::intersect(Curve const &other, Coord eps) const
+{
+ if (auto other_bezier = dynamic_cast<BezierCurve const *>(&other)) {
+ if (other_bezier->order() == 1) {
+ // Use the exact method to intersect a cubic Bézier with a line segment.
+ auto line = Line(other_bezier->initialPoint(), other_bezier->finalPoint());
+ return bezier_line_intersections<3>(*this, line);
+ }
+ }
+ return BezierCurve::intersect(other, eps);
+}
+
template <>
int BezierCurveN<1>::winding(Point const &p) const
{
@@ -358,6 +501,42 @@ void BezierCurveN<3>::feed(PathSink &sink, bool moveto_initial) const
sink.curveTo(controlPoint(1), controlPoint(2), controlPoint(3));
}
+static void bezier_expand_to_image(Rect &range, Point const &x0, Point const &x1, Point const &x2)
+{
+ for (auto i : { X, Y }) {
+ bezier_expand_to_image(range[i], x0[i], x1[i], x2[i]);
+ }
+}
+
+static void bezier_expand_to_image(Rect &range, Point const &x0, Point const &x1, Point const &x2, Point const &x3)
+{
+ for (auto i : { X, Y }) {
+ bezier_expand_to_image(range[i], x0[i], x1[i], x2[i], x3[i]);
+ }
+}
+
+template <>
+void BezierCurveN<1>::expandToTransformed(Rect &bbox, Affine const &transform) const
+{
+ bbox.expandTo(finalPoint() * transform);
+}
+
+template <>
+void BezierCurveN<2>::expandToTransformed(Rect &bbox, Affine const &transform) const
+{
+ bezier_expand_to_image(bbox, controlPoint(0) * transform,
+ controlPoint(1) * transform,
+ controlPoint(2) * transform);
+}
+
+template <>
+void BezierCurveN<3>::expandToTransformed(Rect &bbox, Affine const &transform) const
+{
+ bezier_expand_to_image(bbox, controlPoint(0) * transform,
+ controlPoint(1) * transform,
+ controlPoint(2) * transform,
+ controlPoint(3) * transform);
+}
static Coord bezier_length_internal(std::vector<Point> &v1, Coord tolerance, int level)
{
diff --git a/src/2geom/bezier.cpp b/src/2geom/bezier.cpp
index f9f2beb..aa5a331 100644
--- a/src/2geom/bezier.cpp
+++ b/src/2geom/bezier.cpp
@@ -208,8 +208,6 @@ Bezier &Bezier::operator-=(Bezier const &other)
return *this;
}
-
-
Bezier operator*(Bezier const &f, Bezier const &g)
{
unsigned m = f.order();
@@ -310,7 +308,86 @@ OptInterval bounds_local(Bezier const &b, OptInterval const &i)
}
}
-} // end namespace Geom
+/*
+ * The general bézier of degree n is
+ *
+ * p(t) = sum_{i = 0...n} binomial(n, i) t^i (1 - t)^(n - i) x[i]
+ *
+ * It can be written explicitly as a polynomial in t as
+ *
+ * p(t) = sum_{i = 0...n} binomial(n, i) t^i [ sum_{j = 0...i} binomial(i, j) (-1)^(i - j) x[j] ]
+ *
+ * Its derivative is
+ *
+ * p'(t) = n sum_{i = 1...n} binomial(n - 1, i - 1) t^(i - 1) [ sum_{j = 0...i} binomial(i, j) (-1)^(i - j) x[j] ]
+ *
+ * This is used by the various specialisations below as an optimisation for low degree n <= 3.
+ * In the remaining cases, the generic implementation is used which resorts to iteration.
+ */
+
+void bezier_expand_to_image(Interval &range, Coord x0, Coord x1, Coord x2)
+{
+ range.expandTo(x2);
+
+ if (range.contains(x1)) {
+ // The interval contains all control points, and therefore the entire curve.
+ return;
+ }
+
+ // p'(t) / 2 = at + b
+ auto a = (x2 - x1) - (x1 - x0);
+ auto b = x1 - x0;
+
+ // t = -b / a
+ if (std::abs(a) > EPSILON) {
+ auto t = -b / a;
+ if (t > 0.0 && t < 1.0) {
+ auto s = 1.0 - t;
+ auto x = s * s * x0 + 2 * s * t * x1 + t * t * x2;
+ range.expandTo(x);
+ }
+ }
+}
+
+void bezier_expand_to_image(Interval &range, Coord x0, Coord x1, Coord x2, Coord x3)
+{
+ range.expandTo(x3);
+
+ if (range.contains(x1) && range.contains(x2)) {
+ // The interval contains all control points, and therefore the entire curve.
+ return;
+ }
+
+ // p'(t) / 3 = at^2 + 2bt + c
+ auto a = (x3 - x0) - 3 * (x2 - x1);
+ auto b = (x2 - x1) - (x1 - x0);
+ auto c = x1 - x0;
+
+ auto expand = [&] (Coord t) {
+ if (t > 0.0 && t < 1.0) {
+ auto s = 1.0 - t;
+ auto x = s * s * s * x0 + 3 * s * s * t * x1 + 3 * t * t * s * x2 + t * t * t * x3;
+ range.expandTo(x);
+ }
+ };
+
+ // t = (-b ± sqrt(b^2 - ac)) / a
+ if (std::abs(a) < EPSILON) {
+ if (std::abs(b) > EPSILON) {
+ expand(-c / (2 * b));
+ }
+ } else {
+ auto d2 = b * b - a * c;
+ if (d2 >= 0.0) {
+ auto bsign = b >= 0.0 ? 1 : -1;
+ auto tmp = -(b + bsign * std::sqrt(d2));
+ expand(tmp / a);
+ expand(c / tmp); // Using Vieta's formula: product of roots == c/a
+ }
+ }
+}
+
+} // namespace Geom
/*
Local Variables:
diff --git a/src/2geom/conicsec.cpp b/src/2geom/conicsec.cpp
index d519f4a..0865c0e 100644
--- a/src/2geom/conicsec.cpp
+++ b/src/2geom/conicsec.cpp
@@ -1325,6 +1325,149 @@ bool xAx::decompose (Line& l1, Line& l2) const
return true;
}
+std::array<Line, 2> xAx::decompose_df(Coord epsilon) const
+{
+ // For the classification of degenerate conics, see https://mathworld.wolfram.com/QuadraticCurve.html
+ using std::sqrt, std::abs;
+
+ // Create 2 degenerate lines
+ auto const origin = Point(0, 0);
+ std::array<Line, 2> result = {Line(origin, origin), Line(origin, origin)};
+
+ double A = c[0];
+ double B = c[1];
+ double C = c[2];
+ double D = c[3];
+ double E = c[4];
+ double F = c[5];
+ Coord discriminant = sqr(B) - 4 * A * C;
+ if (discriminant < -epsilon) {
+ return result;
+ }
+
+ bool single_line = false; // In the generic case, there will be 2 lines.
+ bool parallel_lines = false;
+ if (discriminant < epsilon) {
+ discriminant = 0;
+ parallel_lines = true;
+ // Check the secondary discriminant
+ Coord const secondary = sqr(D) + sqr(E) - 4 * F * (A + C);
+ if (secondary < -epsilon) {
+ return result;
+ }
+ single_line = (secondary < epsilon);
+ }
+
+ if (abs(A) > epsilon || abs(C) > epsilon) {
+ // This is the typical case: either x² or y² come with a nonzero coefficient.
+ // To guard against numerical errors, we check which of the coefficients A, C has larger absolute value.
+
+ bool const swap_xy = abs(C) > abs(A);
+ if (swap_xy) {
+ std::swap(A, C);
+ std::swap(D, E);
+ }
+
+ // From now on, we may assume that A is "reasonably large".
+ if (parallel_lines) {
+ if (single_line) {
+ // Special case: a single line.
+ std::array<double, 3> coeffs = {sqrt(abs(A)), sqrt(abs(C)), sqrt(abs(F))};
+ if (swap_xy) {
+ std::swap(coeffs[0], coeffs[1]);
+ }
+ rescale_homogenous(coeffs);
+ result[0].setCoefficients(coeffs[0], coeffs[1], coeffs[2]);
+ return result;
+ }
+
+ // Two parallel lines.
+ Coord const quotient_discriminant = sqr(D) - 4 * A * F;
+ if (quotient_discriminant < 0) {
+ return result;
+ }
+ Coord const sqrt_disc = sqrt(quotient_discriminant);
+ double const c1 = 0.5 * (D - sqrt_disc);
+ double const c2 = c1 + sqrt_disc;
+ std::array<double, 3> coeffs = {A, 0.5 * B, c1};
+ if (swap_xy) {
+ std::swap(coeffs[0], coeffs[1]);
+ }
+ rescale_homogenous(coeffs);
+ result[0].setCoefficients(coeffs[0], coeffs[1], coeffs[2]);
+
+ coeffs = {A, 0.5 * B, c2};
+ if (swap_xy) {
+ std::swap(coeffs[0], coeffs[1]);
+ }
+ rescale_homogenous(coeffs);
+ result[1].setCoefficients(coeffs[0], coeffs[1], coeffs[2]);
+ return result;
+ }
+
+ // Now for the typical case of 2 non-parallel lines.
+
+ // We know that A is further away from 0 than C is.
+ // The mathematical derivation of the solution is as follows:
+ // let Δ = B² - 4AC (the discriminant); we know Δ > 0.
+ // Write δ = sqrt(Δ); we know that this is also positive.
+ // Then the product AΔ is nonzero, so the equation
+ // Ax² + Bxy + Cy² + Dx + Ey + F = 0
+ // is equivalent to
+ // AΔ (Ax² + Bxy + Cy² + Dx + Ey + F) = 0.
+ // Consider the two factors
+ // L_1 = Aδx + 0.5 (Bδ-Δ)y + EA - 0.5 D(B-δ)
+ // L_2 = Aδx + 0.5 (Bδ+Δ)y - EA + 0.5 D(B+δ)
+ // With a bit of algebra, you can show that L_1 * L_2 expands
+ // to AΔ (Ax² + Bxy + Cy² + Dx + Ey + F) (in order to get the
+ // correct value of F, you have to use the fact that the conic
+ // is degenerate). Therefore, the factors L_1 and L_2 are in
+ // fact equations of the two lines to be found.
+ Coord const delta = sqrt(discriminant);
+ std::array<double, 3> coeffs1 = {A * delta, 0.5 * (B * delta - discriminant), E * A - 0.5 * D * (B - delta)};
+ std::array<double, 3> coeffs2 = {coeffs1[0], coeffs1[1] + discriminant, D * delta - coeffs1[2]};
+ if (swap_xy) { // We must unswap the coefficients of x and y
+ std::swap(coeffs1[0], coeffs1[1]);
+ std::swap(coeffs2[0], coeffs2[1]);
+ }
+
+ unsigned index = 0;
+ if (coeffs1[0] != 0 || coeffs1[1] != 0) {
+ rescale_homogenous(coeffs1);
+ result[index++].setCoefficients(coeffs1[0], coeffs1[1], coeffs1[2]);
+ }
+ if (coeffs2[0] != 0 || coeffs2[1] != 0) {
+ rescale_homogenous(coeffs2);
+ result[index].setCoefficients(coeffs2[0], coeffs2[1], coeffs2[2]);
+ }
+ return result;
+ }
+
+ // If we're here, then A==0 and C==0.
+ if (abs(B) < epsilon) { // A == B == C == 0, so the conic reduces to Dx + Ey + F.
+ if (D == 0 && E == 0) {
+ return result;
+ }
+ std::array<double, 3> coeffs = {D, E, F};
+ rescale_homogenous(coeffs);
+ result[0].setCoefficients(coeffs[0], coeffs[1], coeffs[2]);
+ return result;
+ }
+
+ // OK, so A == C == 0 but B != 0. In other words, the conic has the form
+ // Bxy + Dx + Ey + F. Since B != 0, the zero set stays the same if we multiply the
+ // equation by B, which gives us this equation:
+ // B²xy + BDx + BEy + BF = 0.
+ // The above factors as (Bx + E)(By + D) = 0.
+ std::array<double, 2> nonzero_coeffs = {B, E};
+ rescale_homogenous(nonzero_coeffs);
+ result[0].setCoefficients(nonzero_coeffs[0], 0, nonzero_coeffs[1]);
+
+ nonzero_coeffs = {B, D};
+ rescale_homogenous(nonzero_coeffs);
+ result[1].setCoefficients(0, nonzero_coeffs[0], nonzero_coeffs[1]);
+ return result;
+}
/*
* Return the rectangle that bound the conic section arc characterized by
diff --git a/src/2geom/coord.cpp b/src/2geom/coord.cpp
index aea6187..205a82f 100644
--- a/src/2geom/coord.cpp
+++ b/src/2geom/coord.cpp
@@ -76,7 +76,7 @@ namespace Geom {
std::string format_coord_shortest(Coord x)
{
- static double_conversion::DoubleToStringConverter conv(
+ static const double_conversion::DoubleToStringConverter conv(
double_conversion::DoubleToStringConverter::UNIQUE_ZERO,
"inf", "NaN", 'e', -3, 6, 0, 0);
std::string ret(' ', 32);
@@ -88,7 +88,7 @@ std::string format_coord_shortest(Coord x)
std::string format_coord_nice(Coord x)
{
- static double_conversion::DoubleToStringConverter conv(
+ static const double_conversion::DoubleToStringConverter conv(
double_conversion::DoubleToStringConverter::UNIQUE_ZERO,
"inf", "NaN", 'e', -6, 21, 0, 0);
std::string ret(' ', 32);
@@ -100,7 +100,7 @@ std::string format_coord_nice(Coord x)
Coord parse_coord(std::string const &s)
{
- static double_conversion::StringToDoubleConverter conv(
+ static const double_conversion::StringToDoubleConverter conv(
double_conversion::StringToDoubleConverter::ALLOW_LEADING_SPACES |
double_conversion::StringToDoubleConverter::ALLOW_TRAILING_SPACES |
double_conversion::StringToDoubleConverter::ALLOW_SPACES_AFTER_SIGN,
diff --git a/src/2geom/curve.cpp b/src/2geom/curve.cpp
index d0c2849..f79edb3 100644
--- a/src/2geom/curve.cpp
+++ b/src/2geom/curve.cpp
@@ -4,6 +4,7 @@
* MenTaLguY <mental@rydia.net>
* Marco Cecchetti <mrcekets at gmail.com>
* Krzysztof Kosiński <tweenk.pl@gmail.com>
+ * Rafał Siejakowski <rs@rs-math.net>
*
* Copyright 2007-2009 Authors
*
@@ -39,8 +40,6 @@
#include <2geom/ord.h>
#include <2geom/path-sink.h>
-//#include <iostream>
-
namespace Geom
{
@@ -103,46 +102,96 @@ std::vector<CurveIntersection> Curve::intersect(Curve const &/*other*/, Coord /*
std::vector<CurveIntersection> Curve::intersectSelf(Coord eps) const
{
- std::vector<CurveIntersection> result;
- // Monotonic segments cannot have self-intersections.
- // Thus, we can split the curve at roots and intersect the portions.
- std::vector<Coord> splits;
- std::unique_ptr<Curve> deriv(derivative());
- splits = deriv->roots(0, X);
- if (splits.empty()) {
+ /// Represents a sub-arc of the curve.
+ struct Subcurve
+ {
+ std::unique_ptr<Curve> curve;
+ Interval parameter_range;
+
+ Subcurve(Curve *piece, Coord from, Coord to)
+ : curve{piece}
+ , parameter_range{from, to}
+ {}
+ };
+
+ /// A closure to split the curve into portions at the prescribed split points.
+ auto const split_into_subcurves = [=](std::vector<Coord> const &splits) {
+ std::vector<Subcurve> result;
+ result.reserve(splits.size() + 1);
+ Coord previous = 0;
+ for (Coord split : splits) {
+ // Use global EPSILON since we're operating on normalized curve times.
+ if (split < EPSILON || split > 1.0 - EPSILON) {
+ continue;
+ }
+ result.emplace_back(portion(previous, split), previous, split);
+ previous = split;
+ }
+ result.emplace_back(portion(previous, 1.0), previous, 1.0);
return result;
+ };
+
+ /// A closure to find pairwise intersections between the passed subcurves.
+ auto const pairwise_intersect = [=](std::vector<Subcurve> const &subcurves) {
+ std::vector<CurveIntersection> result;
+ for (unsigned i = 0; i < subcurves.size(); i++) {
+ for (unsigned j = i + 1; j < subcurves.size(); j++) {
+ auto const xings = subcurves[i].curve->intersect(*subcurves[j].curve, eps);
+ for (auto const &xing : xings) {
+ // To avoid duplicate intersections, skip values at exactly 1.
+ if (xing.first == 1. || xing.second == 1.) {
+ continue;
+ }
+ Coord const ti = subcurves[i].parameter_range.valueAt(xing.first);
+ Coord const tj = subcurves[j].parameter_range.valueAt(xing.second);
+ result.emplace_back(ti, tj, xing.point());
+ }
+ }
+ }
+ std::sort(result.begin(), result.end());
+ return result;
+ };
+
+ // Monotonic segments cannot have self-intersections. Thus, we can split
+ // the curve at critical points of the X or Y coordinate and intersect
+ // the portions. However, there's the risk that a juncture between two
+ // adjacent portions is mistaken for an intersection due to numerical errors.
+ // Hence, we run the algorithm for both the X and Y coordinates and only
+ // keep the intersections that show up in both intersection lists.
+
+ // Find the critical points of both coordinates.
+ std::unique_ptr<Curve> deriv{derivative()};
+ auto const crits_x = deriv->roots(0, X);
+ auto const crits_y = deriv->roots(0, Y);
+ if (crits_x.empty() || crits_y.empty()) {
+ return {};
}
- deriv.reset();
- splits.push_back(1.);
-
- boost::ptr_vector<Curve> parts;
- Coord previous = 0;
- for (double split : splits) {
- if (split == 0.) continue;
- parts.push_back(portion(previous, split));
- previous = split;
- }
-
- Coord prev_i = 0;
- for (unsigned i = 0; i < parts.size()-1; ++i) {
- Interval dom_i(prev_i, splits[i]);
- prev_i = splits[i];
-
- Coord prev_j = 0;
- for (unsigned j = i+1; j < parts.size(); ++j) {
- Interval dom_j(prev_j, splits[j]);
- prev_j = splits[j];
- std::vector<CurveIntersection> xs = parts[i].intersect(parts[j], eps);
- for (auto & x : xs) {
- // to avoid duplicated intersections, skip values at exactly 1
- if (x.first == 1. || x.second == 1.) continue;
-
- Coord ti = dom_i.valueAt(x.first);
- Coord tj = dom_j.valueAt(x.second);
+ // Split into pieces in two ways and find self-intersections.
+ auto const pieces_x = split_into_subcurves(crits_x);
+ auto const pieces_y = split_into_subcurves(crits_y);
+ auto const crossings_from_x = pairwise_intersect(pieces_x);
+ auto const crossings_from_y = pairwise_intersect(pieces_y);
+ if (crossings_from_x.empty() || crossings_from_y.empty()) {
+ return {};
+ }
- CurveIntersection real(ti, tj, x.point());
- result.push_back(real);
+ // Filter the results, only keeping self-intersections found by both approaches.
+ std::vector<CurveIntersection> result;
+ unsigned index_y = 0;
+ for (auto &&candidate_x : crossings_from_x) {
+ // Find a crossing corresponding to this one in the y-method collection.
+ while (index_y != crossings_from_y.size()) {
+ auto const gap = crossings_from_y[index_y].first - candidate_x.first;
+ if (std::abs(gap) < EPSILON) {
+ // We found the matching intersection!
+ result.emplace_back(candidate_x);
+ index_y++;
+ break;
+ } else if (gap < 0.0) {
+ index_y++;
+ } else {
+ break;
}
}
}
diff --git a/src/2geom/ellipse.cpp b/src/2geom/ellipse.cpp
index f076cfa..5ec7c38 100644
--- a/src/2geom/ellipse.cpp
+++ b/src/2geom/ellipse.cpp
@@ -31,6 +31,7 @@
* the specific language governing rights and limitations.
*/
+#include <2geom/conicsec.h>
#include <2geom/ellipse.h>
#include <2geom/elliptical-arc.h>
#include <2geom/numeric/fitting-tool.h>
@@ -161,6 +162,17 @@ Rect Ellipse::boundsExact() const
return result;
}
+Rect Ellipse::boundsFast() const
+{
+ // Every ellipse is contained in the circle with the same center and radius
+ // equal to the larger of the two rays. We return the bounding square
+ // of this circle (this is really fast but only exact for circles).
+ auto const larger_ray = (ray(X) > ray(Y) ? ray(X) : ray(Y));
+ assert(larger_ray >= 0.0);
+ auto const rr = Point(larger_ray, larger_ray);
+ return Rect(_center - rr, _center + rr);
+}
+
std::vector<double> Ellipse::coefficients() const
{
std::vector<double> c(6);
@@ -242,9 +254,6 @@ Ellipse::arc(Point const &ip, Point const &inner, Point const &fp)
large_arc_flag = true;
}
- //cross(-iv, fv) && large_arc_flag
-
-
// Determination of sweep flag:
// For clarity, let's assume that Y grows up. Then the cross product
// is positive for points on the left side of a vector and negative
@@ -403,25 +412,77 @@ bool Ellipse::contains(Point const &p) const
return tp.length() <= 1;
}
-std::vector<ShapeIntersection> Ellipse::intersect(Line const &line) const
+/** @brief Convert curve time on the major axis to the corresponding angle
+ * parameters on a degenerate ellipse collapsed onto that axis.
+ * @param t The curve time on the major axis of an ellipse.
+ * @param vertical If true, the major axis goes from angle -π/2 to +π/2;
+ * otherwise, the major axis connects angles π and 0.
+ * @return The two angles at which the collapsed ellipse passes through the
+ * major axis point corresponding to the given time \f$t \in [0, 1]\f$.
+ */
+static std::array<Coord, 2> axis_time_to_angles(Coord t, bool vertical)
{
+ Coord const to_unit = std::clamp(2.0 * t - 1.0, -1.0, 1.0);
+ if (vertical) {
+ double const arcsin = std::asin(to_unit);
+ return {arcsin, M_PI - arcsin};
+ } else {
+ double const arccos = std::acos(to_unit);
+ return {arccos, -arccos};
+ }
+}
+/** @brief For each intersection of some shape with the major axis of an ellipse, produce one or two
+ * intersections of a degenerate ellipse (collapsed onto that axis) with the same shape.
+ *
+ * @param axis_intersections The intersections of some shape with the major axis.
+ * @param vertical Whether this is the vertical major axis (in the ellipse's natural coordinates).
+ * @return A vector with doubled intersections (corresponding to the two passages of the squashed
+ * ellipse through that point) and swapped order of the intersected shapes.
+*/
+static std::vector<ShapeIntersection> double_axis_intersections(std::vector<ShapeIntersection> &&axis_intersections,
+ bool vertical)
+{
+ if (axis_intersections.empty()) {
+ return {};
+ }
std::vector<ShapeIntersection> result;
+ result.reserve(2 * axis_intersections.size());
- if (line.isDegenerate()) return result;
- if (ray(X) == 0 || ray(Y) == 0) {
- // TODO intersect with line segment.
+ for (auto const &x : axis_intersections) {
+ for (auto a : axis_time_to_angles(x.second, vertical)) {
+ result.emplace_back(a, x.first, x.point()); // Swap first <-> converted second.
+ if (x.second == 0.0 || x.second == 1.0) {
+ break; // Do not double up endpoint intersections.
+ }
+ }
+ }
+ return result;
+}
+
+std::vector<ShapeIntersection> Ellipse::intersect(Line const &line) const
+{
+ std::vector<ShapeIntersection> result;
+
+ if (line.isDegenerate()) {
return result;
}
+ if (ray(X) == 0 || ray(Y) == 0) {
+ return double_axis_intersections(line.intersect(majorAxis()), ray(X) == 0);
+ }
// Ax^2 + Bxy + Cy^2 + Dx + Ey + F
- Coord A, B, C, D, E, F;
- coefficients(A, B, C, D, E, F);
+ std::array<Coord, 6> coeffs;
+ coefficients(coeffs[0], coeffs[1], coeffs[2], coeffs[3], coeffs[4], coeffs[5]);
+ rescale_homogenous(coeffs);
+ auto [A, B, C, D, E, F] = coeffs;
Affine iuct = inverseUnitCircleTransform();
// generic case
- Coord a, b, c;
- line.coefficients(a, b, c);
+ std::array<Coord, 3> line_coeffs;
+ line.coefficients(line_coeffs[0], line_coeffs[1], line_coeffs[2]);
+ rescale_homogenous(line_coeffs);
+ auto [a, b, c] = line_coeffs;
Point lv = line.vector();
if (fabs(lv[X]) > fabs(lv[Y])) {
@@ -458,6 +519,10 @@ std::vector<ShapeIntersection> Ellipse::intersect(Line const &line) const
std::vector<ShapeIntersection> Ellipse::intersect(LineSegment const &seg) const
{
+ if (!boundsFast().intersects(seg.boundsFast())) {
+ return {};
+ }
+
// we simply re-use the procedure for lines and filter out
// results where the line time value is outside of the unit interval.
std::vector<ShapeIntersection> result = intersect(Line(seg));
@@ -467,49 +532,80 @@ std::vector<ShapeIntersection> Ellipse::intersect(LineSegment const &seg) const
std::vector<ShapeIntersection> Ellipse::intersect(Ellipse const &other) const
{
- // handle degenerate cases first
- if (ray(X) == 0 || ray(Y) == 0) {
-
+ // Handle degenerate cases first.
+ if (ray(X) == 0 || ray(Y) == 0) { // Degenerate ellipse, collapsed to the major axis.
+ return double_axis_intersections(other.intersect(majorAxis()), ray(X) == 0);
+ }
+ if (*this == other) { // Two identical ellipses.
+ THROW_INFINITELY_MANY_SOLUTIONS("The two ellipses are identical.");
+ }
+ if (!boundsFast().intersects(other.boundsFast())) {
+ return {};
}
- // intersection of two ellipses can be solved analytically.
- // http://maptools.home.comcast.net/~maptools/BivariateQuadratics.pdf
-
- Coord A, B, C, D, E, F;
- Coord a, b, c, d, e, f;
-
- // NOTE: the order of coefficients is different to match the convention in the PDF above
- // Ax^2 + Bx^2 + Cx + Dy + Exy + F
- this->coefficients(A, E, B, C, D, F);
- other.coefficients(a, e, b, c, d, f);
- // Assume that Q is the ellipse equation given by uppercase letters
- // and R is the equation given by lowercase ones. An intersection exists when
- // there is a coefficient mu such that
- // mu Q + R = 0
+ // Find coefficients of the implicit equations of the two ellipses and rescale
+ // them (losslessly) for better numerical conditioning.
+ std::array<double, 6> coeffs;
+ coefficients(coeffs[0], coeffs[1], coeffs[2], coeffs[3], coeffs[4], coeffs[5]);
+ rescale_homogenous(coeffs);
+ auto [A, B, C, D, E, F] = coeffs;
+
+ std::array<double, 6> otheffs;
+ other.coefficients(otheffs[0], otheffs[1], otheffs[2], otheffs[3], otheffs[4], otheffs[5]);
+ rescale_homogenous(otheffs);
+ auto [a, b, c, d, e, f] = otheffs;
+
+ // Assume that Q(x, y) = 0 is the ellipse equation given by uppercase letters
+ // and R(x, y) = 0 is the equation given by lowercase ones.
+ // In other words, Q is the quadratic function describing this ellipse and
+ // R is the quadratic function for the other ellipse.
//
- // This can be written in the following way:
+ // A point (x, y) is common to both ellipses if and only if it solves the system
+ // { Q(x, y) = 0,
+ // { R(x, y) = 0.
//
- // | ff cc/2 dd/2 | |1|
- // mu Q + R = [1 x y] | cc/2 aa ee/2 | |x| = 0
- // | dd/2 ee/2 bb | |y|
+ // If µ is any real number, we can multiply the first equation by µ and add that
+ // to the first equation, obtaining the new system of equations:
+ // { Q(x, y) = 0,
+ // { µQ(x, y) + R(x, y) = 0.
//
- // where aa = mu A + a and so on. The determinant can be explicitly written out,
- // giving an equation which is cubic in mu and can be solved analytically.
-
- Coord I, J, K, L;
- I = (-E*E*F + 4*A*B*F + C*D*E - A*D*D - B*C*C) / 4;
- J = -((E*E - 4*A*B) * f + (2*E*F - C*D) * e + (2*A*D - C*E) * d +
- (2*B*C - D*E) * c + (C*C - 4*A*F) * b + (D*D - 4*B*F) * a) / 4;
- K = -((e*e - 4*a*b) * F + (2*e*f - c*d) * E + (2*a*d - c*e) * D +
- (2*b*c - d*e) * C + (c*c - 4*a*f) * B + (d*d - 4*b*f) * A) / 4;
- L = (-e*e*f + 4*a*b*f + c*d*e - a*d*d - b*c*c) / 4;
+ // The first equation still says that (x, y) is a point on this ellipse, but the
+ // second equation uses the new expression (µQ + R) instead of the original R.
+ //
+ // Why do we do this? The reason is that the set of functions {µQ + R : µ real}
+ // is a "real system of conics" and there's a theorem which guarantees that such a system
+ // always contains a "degenerate conic" [proof below].
+ // In practice, the degenerate conic will describe a line or a pair of lines, and intersecting
+ // a line with an ellipse is much easier than intersecting two ellipses directly.
+ //
+ // But in order to be able to do this, we must find a value of µ for which µQ + R is degenerate.
+ // We can write the expression (µQ + R)(x, y) in the following way:
+ //
+ // | aa bb/2 dd/2 | |x|
+ // (µQ + R)(x, y) = [x y 1] | bb/2 cc ee/2 | |y|
+ // | dd/2 ee/2 ff | |1|
+ //
+ // where aa = µA + a and so on. The determinant can be explicitly written out,
+ // giving an equation which is cubic in µ and can be solved analytically.
+ // The conic µQ + R is degenerate if and only if this determinant is 0.
+ //
+ // Proof that there's always a degenerate conic: a cubic real polynomial always has a root,
+ // and if the polynomial in µ isn't cubic (coefficient of µ^3 is zero), then the starting
+ // conic is already degenerate.
+
+ Coord I, J, K, L; // Coefficients of µ in the expression for the determinant.
+ I = (-B*B*F + 4*A*C*F + D*E*B - A*E*E - C*D*D) / 4;
+ J = -((B*B - 4*A*C) * f + (2*B*F - D*E) * b + (2*A*E - D*B) * e +
+ (2*C*D - E*B) * d + (D*D - 4*A*F) * c + (E*E - 4*C*F) * a) / 4;
+ K = -((b*b - 4*a*c) * F + (2*b*f - d*e) * B + (2*a*e - d*b) * E +
+ (2*c*d - e*b) * D + (d*d - 4*a*f) * C + (e*e - 4*c*f) * A) / 4;
+ L = (-b*b*f + 4*a*c*f + d*e*b - a*e*e - c*d*d) / 4;
std::vector<Coord> mus = solve_cubic(I, J, K, L);
Coord mu = infinity();
- std::vector<ShapeIntersection> result;
- // Now that we have solved for mu, we need to check whether the conic
- // determined by mu Q + R is reducible to a product of two lines. If it's not,
+ // Now that we have solved for µ, we need to check whether the conic
+ // determined by µQ + R is reducible to a product of two lines. If it's not,
// it means that there are no intersections. If it is, the intersections of these
// lines with the original ellipses (if there are any) give the coordinates
// of intersections.
@@ -521,68 +617,72 @@ std::vector<ShapeIntersection> Ellipse::intersect(Ellipse const &other) const
if (mus.size() == 3) {
std::swap(mus[1], mus[0]);
}
- for (double i : mus) {
- Coord aa = i * A + a;
- Coord bb = i * B + b;
- Coord ee = i * E + e;
- Coord delta = ee*ee - 4*aa*bb;
- if (delta < 0) continue;
- mu = i;
+ /// Discriminant within this radius of 0 will be considered zero.
+ static Coord const discriminant_precision = 1e-10;
+
+ for (Coord candidate_mu : mus) {
+ Coord const aa = std::fma(candidate_mu, A, a);
+ Coord const bb = std::fma(candidate_mu, B, b);
+ Coord const cc = std::fma(candidate_mu, C, c);
+ Coord const delta = sqr(bb) - 4*aa*cc;
+ if (delta < -discriminant_precision) {
+ continue;
+ }
+ mu = candidate_mu;
break;
}
// if no suitable mu was found, there are no intersections
- if (mu == infinity()) return result;
-
- Coord aa = mu * A + a;
- Coord bb = mu * B + b;
- Coord cc = mu * C + c;
- Coord dd = mu * D + d;
- Coord ee = mu * E + e;
- Coord ff = mu * F + f;
-
- unsigned line_num = 0;
- Line lines[2];
-
- if (aa != 0) {
- bb /= aa; cc /= aa; dd /= aa; ee /= aa; /*ff /= aa;*/
- Coord s = (ee + std::sqrt(ee*ee - 4*bb)) / 2;
- Coord q = ee - s;
- Coord alpha = (dd - cc*q) / (s - q);
- Coord beta = cc - alpha;
-
- line_num = 2;
- lines[0] = Line(1, q, alpha);
- lines[1] = Line(1, s, beta);
- } else if (bb != 0) {
- cc /= bb; /*dd /= bb;*/ ee /= bb; ff /= bb;
- Coord s = ee;
- Coord q = 0;
- Coord alpha = cc / ee;
- Coord beta = ff * ee / cc;
-
- line_num = 2;
- lines[0] = Line(q, 1, alpha);
- lines[1] = Line(s, 1, beta);
- } else if (ee != 0) {
- line_num = 2;
- lines[0] = Line(ee, 0, dd);
- lines[1] = Line(0, 1, cc/ee);
- } else if (cc != 0 || dd != 0) {
- line_num = 1;
- lines[0] = Line(cc, dd, ff);
+ if (mu == infinity()) {
+ return {};
}
+ // Create the degenerate conic and decompose it into lines.
+ std::array<double, 6> degen = {std::fma(mu, A, a), std::fma(mu, B, b), std::fma(mu, C, c),
+ std::fma(mu, D, d), std::fma(mu, E, e), std::fma(mu, F, f)};
+ rescale_homogenous(degen);
+ auto const lines = xAx(degen[0], degen[1], degen[2],
+ degen[3], degen[4], degen[5]).decompose_df(discriminant_precision);
+
// intersect with the obtained lines and report intersections
- for (unsigned li = 0; li < line_num; ++li) {
- std::vector<ShapeIntersection> as = intersect(lines[li]);
- std::vector<ShapeIntersection> bs = other.intersect(lines[li]);
-
- if (!as.empty() && as.size() == bs.size()) {
- for (unsigned i = 0; i < as.size(); ++i) {
- ShapeIntersection ix(as[i].first, bs[i].first,
- middle_point(as[i].point(), bs[i].point()));
- result.push_back(ix);
+ std::vector<ShapeIntersection> result;
+ for (auto const &line : lines) {
+ if (line.isDegenerate()) {
+ continue;
+ }
+ auto as = intersect(line);
+ // NOTE: If we only cared about the intersection points, we could simply
+ // intersect this ellipse with the lines and ignore the other ellipse.
+ // But we need the time coordinates on the other ellipse as well.
+ auto bs = other.intersect(line);
+ if (as.empty() || bs.empty()) {
+ continue;
+ }
+ // Due to numerical errors, a tangency may sometimes be found as 1 intersection
+ // on one ellipse and 2 intersections on the other. If this happens, we average
+ // the points of the two intersections.
+ auto const intersection_average = [](ShapeIntersection const &i,
+ ShapeIntersection const &j) -> ShapeIntersection
+ {
+ return ShapeIntersection(i.first, j.first, middle_point(i.point(), j.point()));
+ };
+ auto const synthesize_intersection = [&](ShapeIntersection const &i,
+ ShapeIntersection const &j) -> void
+ {
+ result.emplace_back(i.first, j.first, middle_point(i.point(), j.point()));
+ };
+ if (as.size() == 2) {
+ if (bs.size() == 2) {
+ synthesize_intersection(as[0], bs[0]);
+ synthesize_intersection(as[1], bs[1]);
+ } else if (bs.size() == 1) {
+ synthesize_intersection(intersection_average(as[0], as[1]), bs[0]);
+ }
+ } else if (as.size() == 1) {
+ if (bs.size() == 2) {
+ synthesize_intersection(as[0], intersection_average(bs[0], bs[1]));
+ } else if (bs.size() == 1) {
+ synthesize_intersection(as[0], bs[0]);
}
}
}
@@ -594,6 +694,7 @@ std::vector<ShapeIntersection> Ellipse::intersect(D2<Bezier> const &b) const
Coord A, B, C, D, E, F;
coefficients(A, B, C, D, E, F);
+ // We plug the X and Y curves into the implicit equation and solve for t.
Bezier x = A*b[X]*b[X] + B*b[X]*b[Y] + C*b[Y]*b[Y] + D*b[X] + E*b[Y] + F;
std::vector<Coord> r = x.roots();
diff --git a/src/2geom/elliptical-arc.cpp b/src/2geom/elliptical-arc.cpp
index 96837bc..f388480 100644
--- a/src/2geom/elliptical-arc.cpp
+++ b/src/2geom/elliptical-arc.cpp
@@ -54,7 +54,7 @@ namespace Geom
* @brief Elliptical arc curve
*
* Elliptical arc is a curve taking the shape of a section of an ellipse.
- *
+ *
* The arc function has two forms: the regular one, mapping the unit interval to points
* on 2D plane (the linear domain), and a second form that maps some interval
* \f$A \subseteq [0,2\pi)\f$ to the same points (the angular domain). The interval \f$A\f$
@@ -100,9 +100,9 @@ namespace Geom
* by computing a partial derivative with respect to the angle
* and equating that to zero:
* \f{align*}{
- x &= r_x \cos \varphi \cos \theta - r_y \sin \varphi \sin \theta + c_x \\
- \frac{\partial x}{\partial \theta} &= -r_x \cos \varphi \sin \theta - r_y \sin \varphi \cos \theta = 0 \\
- \frac{\sin \theta}{\cos \theta} &= \tan\theta = -\frac{r_y \sin \varphi}{r_x \cos \varphi} \\
+ x &= r_x \cos \varphi \cos \theta - r_y \sin \varphi \sin \theta + c_x \\
+ \frac{\partial x}{\partial \theta} &= -r_x \cos \varphi \sin \theta - r_y \sin \varphi \cos \theta = 0 \\
+ \frac{\sin \theta}{\cos \theta} &= \tan\theta = -\frac{r_y \sin \varphi}{r_x \cos \varphi} \\
\theta &= \tan^{-1} \frac{-r_y \sin \varphi}{r_x \cos \varphi}
\f}
* The local extremes correspond to two angles separated by \f$\pi\f$.
@@ -146,6 +146,12 @@ Rect EllipticalArc::boundsExact() const
return result;
}
+void EllipticalArc::expandToTransformed(Rect &bbox, Affine const &transform) const
+{
+ auto c = *this;
+ c *= transform;
+ bbox |= c.boundsExact();
+}
Point EllipticalArc::pointAtAngle(Coord t) const
{
@@ -302,16 +308,20 @@ Coord EllipticalArc::valueAt(Coord t, Dim2 d) const
Curve* EllipticalArc::portion(double f, double t) const
{
// fix input arguments
- if (f < 0) f = 0;
- if (f > 1) f = 1;
- if (t < 0) t = 0;
- if (t > 1) t = 1;
+ f = std::clamp(f, 0.0, 1.0);
+ t = std::clamp(t, 0.0, 1.0);
if (f == t) {
EllipticalArc *arc = new EllipticalArc();
arc->_initial_point = arc->_final_point = pointAt(f);
return arc;
}
+ if (f == 0.0 && t == 1.0) {
+ return duplicate();
+ }
+ if (f == 1.0 && t == 0.0) {
+ return reverse();
+ }
EllipticalArc *arc = static_cast<EllipticalArc*>(duplicate());
arc->_initial_point = pointAt(f);
@@ -625,8 +635,14 @@ std::vector<CurveIntersection> EllipticalArc::intersect(Curve const &other, Coor
}
if (auto arc = dynamic_cast<EllipticalArc const *>(&other)) {
- auto result = _filterIntersections(_ellipse.intersect(arc->_ellipse), true);
- return arc->_filterIntersections(std::move(result), false);
+ std::vector<CurveIntersection> crossings;
+ try {
+ crossings = _ellipse.intersect(arc->_ellipse);
+ } catch (InfinitelyManySolutions &) {
+ // This could happen if the two arcs come from the same ellipse.
+ return _intersectSameEllipse(arc);
+ }
+ return arc->_filterIntersections(_filterIntersections(std::move(crossings), true), false);
}
// in case someone wants to make a custom curve type
@@ -635,6 +651,43 @@ std::vector<CurveIntersection> EllipticalArc::intersect(Curve const &other, Coor
return result;
}
+/** @brief Check if two arcs on the same ellipse intersect/overlap.
+ *
+ * @param other Another elliptical arc on the same ellipse as this one.
+ * @return If the arcs overlap, the returned vector contains synthesized intersections
+ * at the start and end of the overlap.
+ * If the arcs do not overlap, an empty vector is returned.
+ */
+std::vector<ShapeIntersection> EllipticalArc::_intersectSameEllipse(EllipticalArc const *other) const
+{
+ assert(_ellipse == other->_ellipse);
+ auto const &other_angles = other->angularInterval();
+ std::vector<ShapeIntersection> result;
+
+ /// A closure to create an "intersection" at the prescribed angle.
+ auto const synthesize_intersection = [&](Angle angle) {
+ auto const time = timeAtAngle(angle);
+ if (result.end() == std::find_if(result.begin(), result.end(),
+ [=](ShapeIntersection const &xing) -> bool {
+ return xing.first == time;
+ }))
+ {
+ result.emplace_back(time, other->timeAtAngle(angle), _ellipse.pointAt(angle));
+ }
+ };
+
+ for (auto a : {_angles.initialAngle(), _angles.finalAngle()}) {
+ if (other_angles.contains(a)) {
+ synthesize_intersection(a);
+ }
+ }
+ for (auto a : {other_angles.initialAngle(), other_angles.finalAngle()}) {
+ if (_angles.contains(a)) {
+ synthesize_intersection(a);
+ }
+ }
+ return result;
+}
void EllipticalArc::_updateCenterAndAngles()
{
diff --git a/src/2geom/intersection-graph.cpp b/src/2geom/intersection-graph.cpp
index a294c72..524267e 100644
--- a/src/2geom/intersection-graph.cpp
+++ b/src/2geom/intersection-graph.cpp
@@ -40,35 +40,30 @@
namespace Geom {
+/// Function object for comparing intersection vertices based on the intersection time.
struct PathIntersectionGraph::IntersectionVertexLess {
bool operator()(IntersectionVertex const &a, IntersectionVertex const &b) const {
return a.pos < b.pos;
}
};
-/** @class PathIntersectionGraph
- * @brief Intermediate data for computing Boolean operations on paths.
- *
- * This class implements the Greiner-Hormann clipping algorithm,
- * with improvements inspired by Foster and Overfelt as well as some
- * original contributions.
- *
- * @ingroup Paths
- */
-
PathIntersectionGraph::PathIntersectionGraph(PathVector const &a, PathVector const &b, Coord precision)
: _graph_valid(true)
{
- if (a.empty() || b.empty()) return;
-
_pv[0] = a;
_pv[1] = b;
+ if (a.empty() || b.empty()) return;
+
_prepareArguments();
bool has_intersections = _prepareIntersectionLists(precision);
if (!has_intersections) return;
_assignEdgeWindingParities(precision);
+
+ // If a path has only degenerate intersections, assign its status now.
+ // This protects against later accidentally picking a point for winding
+ // determination that is exactly at a removed intersection.
_assignComponentStatusFromDegenerateIntersections();
_removeDegenerateIntersections();
if (_graph_valid) {
@@ -76,6 +71,9 @@ PathIntersectionGraph::PathIntersectionGraph(PathVector const &a, PathVector con
}
}
+/** Prepare the operands stored in PathIntersectionGraph::_pv by closing all of their constituent
+ * paths and removing degenerate segments from them.
+ */
void PathIntersectionGraph::_prepareArguments()
{
// all paths must be closed, otherwise we will miss some intersections
@@ -100,6 +98,10 @@ void PathIntersectionGraph::_prepareArguments()
}
}
+/** @brief Compute the lists of intersections between the constituent paths of both operands.
+ * @param precision – the precision setting for the sweepline algorithm.
+ * @return Whether any intersections were found.
+ */
bool PathIntersectionGraph::_prepareIntersectionLists(Coord precision)
{
std::vector<PVIntersection> pxs = _pv[0].intersect(_pv[1], precision);
@@ -134,7 +136,7 @@ bool PathIntersectionGraph::_prepareIntersectionLists(Coord precision)
_components[1][xb->pos.path_index].xlist.push_back(*xb);
}
- // sort components according to time value of intersections
+ // sort intersections in each component according to time value
for (auto & _component : _components) {
for (std::size_t i = 0; i < _component.size(); ++i) {
_component[i].xlist.sort(IntersectionVertexLess());
@@ -144,18 +146,21 @@ bool PathIntersectionGraph::_prepareIntersectionLists(Coord precision)
return true;
}
+/** Determine whether path portions between consecutive intersections lie inside or outside
+ * of the other path-vector.
+ */
void PathIntersectionGraph::_assignEdgeWindingParities(Coord precision)
{
- // determine the winding numbers of path portions between intersections
for (unsigned w = 0; w < 2; ++w) {
- unsigned ow = (w+1) % 2;
+ unsigned ow = (w+1) % 2; ///< The index of the other operand
- for (unsigned li = 0; li < _components[w].size(); ++li) {
+ for (unsigned li = 0; li < _components[w].size(); ++li) { // Traverse all paths in the component
IntersectionList &xl = _components[w][li].xlist;
- for (ILIter i = xl.begin(); i != xl.end(); ++i) {
+ for (ILIter i = xl.begin(); i != xl.end(); ++i) { // Traverse all intersections in the path
ILIter n = cyclic_next(i, xl);
std::size_t pi = i->pos.path_index;
+ /// Path time interval from the current crossing to the next one
PathInterval ival = forward_interval(i->pos, n->pos, _pv[w][pi].size());
PathTime mid = ival.inside(precision);
@@ -172,11 +177,11 @@ void PathIntersectionGraph::_assignEdgeWindingParities(Coord precision)
}
}
+/** Detect the situation where a path is either entirely inside or entirely outside of the other
+ * path-vector and set the status flag accordingly.
+ */
void PathIntersectionGraph::_assignComponentStatusFromDegenerateIntersections()
{
- // If a path has only degenerate intersections, assign its status now.
- // This protects against later accidentally picking a point for winding
- // determination that is exactly at a removed intersection.
for (auto & _component : _components) {
for (unsigned li = 0; li < _component.size(); ++li) {
IntersectionList &xl = _component[li].xlist;
@@ -196,32 +201,36 @@ void PathIntersectionGraph::_assignComponentStatusFromDegenerateIntersections()
}
}
+/** Remove intersections that don't change between in/out.
+ *
+ * In general, a degenerate intersection can happen at a point where
+ * two shapes "kiss" (are tangent) but do not cross into each other.
+ */
void PathIntersectionGraph::_removeDegenerateIntersections()
{
- // remove intersections that don't change between in/out
for (auto & _component : _components) {
for (unsigned li = 0; li < _component.size(); ++li) {
IntersectionList &xl = _component[li].xlist;
for (ILIter i = xl.begin(); i != xl.end();) {
ILIter n = cyclic_next(i, xl);
- if (i->next_edge == n->next_edge) {
- bool last_node = (i == n);
+ if (i->next_edge == n->next_edge) { // Both edges inside or both outside
+ bool last_node = (i == n); ///< Whether this is the last remaining crossing.
ILIter nn = _getNeighbor(n);
IntersectionList &oxl = _getPathData(nn).xlist;
// When exactly 3 out of 4 edges adjacent to an intersection
// have the same winding, we have a defective intersection,
// which is neither degenerate nor normal. Those can occur in paths
- // that contain overlapping segments. We cannot handle that case
- // for now, so throw an exception.
+ // that contain overlapping segments.
if (cyclic_prior(nn, oxl)->next_edge != nn->next_edge) {
+ // Not a backtrack - set the defective flag.
_graph_valid = false;
n->defective = true;
nn->defective = true;
++i;
continue;
}
-
+ // Erase the degenerate or defective crossings
oxl.erase(nn);
xl.erase(n);
if (last_node) break;
@@ -233,6 +242,9 @@ void PathIntersectionGraph::_removeDegenerateIntersections()
}
}
+/** Verify that all paths contain an even number of intersections and that
+ * the intersection graph does not contain leaves (degree one vertices).
+ */
void PathIntersectionGraph::_verify()
{
#ifndef NDEBUG
@@ -335,12 +347,25 @@ void PathIntersectionGraph::fragments(PathVector &in, PathVector &out) const
}
}
+/** @brief Compute the partial result of a boolean operation by looking at components containing
+ * intersections and stitching the correct path portions between them, depending on the truth
+ * table of the operation.
+ *
+ * @param enter_a – whether the path portions contained inside operand A should be part of the boundary
+ * of the boolean operation's result.
+ * @param enter_b – whether the path portions contained inside operand B should be part of the boundary
+ * of the boolean operation's result.
+ *
+ * These two flags completely determine how to resolve the crossings when building the result
+ * and therefore encode which boolean operation we are performing. For example, the boolean intersection
+ * corresponds to enter_a == true and enter_b == true, as can be seen by looking at a Venn diagram.
+ */
PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
{
PathVector result;
if (_xs.empty()) return result;
- // reset processed status
+ // Create the list of intersections to process
_ulist.clear();
for (auto & _component : _components) {
for (auto & li : _component) {
@@ -361,19 +386,18 @@ PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
result.push_back(Path(i->p));
result.back().setStitching(true);
- bool reverse = false;
+ bool reverse = false; ///< Whether to traverse the current component in the backwards direction.
while (i->_proc_hook.is_linked()) {
ILIter prev = i;
- std::size_t pi = i->pos.path_index;
+ std::size_t pi = i->pos.path_index; ///< Index of the path in its PathVector
// determine which direction to go
// union: always go outside
// intersection: always go inside
// a minus b: go inside in b, outside in a
// b minus a: go inside in a, outside in b
- reverse = false;
- if (w == 0) {
+ if (w == 0) { // The path we're on is a part of A
reverse = (i->next_edge == INSIDE) ^ enter_a;
- } else {
+ } else { // The path we're on is a part of B
reverse = (i->next_edge == INSIDE) ^ enter_b;
}
@@ -384,16 +408,14 @@ PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
i = cyclic_next(i, _components[w][pi].xlist);
}
- // append portion of path
+ // append portion of path to the result
PathInterval ival = PathInterval::from_direction(
prev->pos.asPathTime(), i->pos.asPathTime(),
reverse, _pv[i->which][pi].size());
_pv[i->which][pi].appendPortionTo(result.back(), ival, prev->p, i->p);
- // mark both vertices as processed
- //prev->processed = true;
- //i->processed = true;
+ // count both vertices as processed
n_processed += 2;
if (prev->_proc_hook.is_linked()) {
_ulist.erase(_ulist.iterator_to(*prev));
@@ -424,14 +446,20 @@ PathVector PathIntersectionGraph::_getResult(bool enter_a, bool enter_b)
return result;
}
+/** @brief Select intersection-free path components ahead of a boolean operation based on whether
+ * they should be a part of that operation's result.
+ *
+ * Every component that has intersections will be processed by _getResult().
+ * Here we take care of paths that don't have any intersections. They are either
+ * completely inside or completely outside the other path-vector.
+ *
+ * @param result – output parameter to store the selected components.
+ * @param which – which of the two operands to search for intersection-free paths.
+ * @param inside – If set to true, add paths entirely contained inside the other path-vector to
+ * the result. If set to false, add paths entirely outside of the other path-vector instead.
+ */
void PathIntersectionGraph::_handleNonintersectingPaths(PathVector &result, unsigned which, bool inside)
{
- /* Every component that has any intersections will be processed by _getResult.
- * Here we take care of paths that don't have any intersections. They are either
- * completely inside or completely outside the other pathvector. We test this by
- * evaluating the winding rule at the initial point. If inside is true and
- * the path is inside, we add it to the result.
- */
unsigned w = which;
unsigned ow = (w+1) % 2;
@@ -442,12 +470,13 @@ void PathIntersectionGraph::_handleNonintersectingPaths(PathVector &result, unsi
if (has_path_data && !_components[w][i].xlist.empty()) continue;
bool path_inside = false;
- // Use the in/out determination from constructor, if available
+ // Use the status flag set in the constructor if available.
if (has_path_data && _components[w][i].status == INSIDE) {
path_inside = true;
} else if (has_path_data && _components[w][i].status == OUTSIDE) {
path_inside = false;
} else {
+ // The status flag is ambiguous: we evaluate the winding number of the initial point.
int wdg = _pv[ow].winding(_pv[w][i].initialPoint());
path_inside = wdg % 2 != 0;
}
@@ -458,18 +487,25 @@ void PathIntersectionGraph::_handleNonintersectingPaths(PathVector &result, unsi
}
}
+/** @brief Get an iterator to the corresponding crossing on the other path-vector.
+ *
+ * @param ILIter – an iterator to a list of intersections in one of the path-vectors.
+ * @return An iterator to the corresponding intersection in the other path-vector.
+ */
PathIntersectionGraph::ILIter PathIntersectionGraph::_getNeighbor(ILIter iter)
{
unsigned ow = (iter->which + 1) % 2;
return _components[ow][iter->neighbor->pos.path_index].xlist.iterator_to(*iter->neighbor);
}
+/** Get the path data for the path containing the intersection given an iterator to the intersection */
PathIntersectionGraph::PathData &
PathIntersectionGraph::_getPathData(ILIter iter)
{
return _components[iter->which][iter->pos.path_index];
}
+/** Format the PathIntersectionGraph for output. */
std::ostream &operator<<(std::ostream &os, PathIntersectionGraph const &pig)
{
os << "Intersection graph:\n"
diff --git a/src/2geom/parting-point.cpp b/src/2geom/parting-point.cpp
new file mode 100644
index 0000000..3e3e803
--- /dev/null
+++ b/src/2geom/parting-point.cpp
@@ -0,0 +1,280 @@
+/** @file Implementation of parting_point(Path const&, Path const&, Coord)
+ */
+/* An algorithm to find the first parting point of two paths.
+ *
+ * Authors:
+ * Rafał Siejakowski <rs@rs-math.net>
+ *
+ * Copyright 2022 the Authors.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#include <2geom/path.h>
+#include <2geom/point.h>
+
+namespace Geom
+{
+
+PathIntersection parting_point(Path const &first, Path const &second, Coord precision)
+{
+ Path const *paths[2] = { &first, &second };
+ Point const starts[2] = { first.initialPoint(), second.initialPoint() };
+
+ if (!are_near(starts[0], starts[1], precision)) {
+ auto const invalid = PathTime(0, -1.0);
+ return PathIntersection(invalid, invalid, middle_point(starts[0], starts[1]));
+ }
+
+ if (first.empty() || second.empty()) {
+ auto const start_time = PathTime(0, 0.0);
+ return PathIntersection(start_time, start_time, middle_point(starts[0], starts[1]));
+ }
+
+ size_t const curve_count[2] = { first.size(), second.size() };
+ Coord const max_time[2] = { first.timeRange().max(), second.timeRange().max() };
+
+ /// Curve indices up until which the paths are known to overlap
+ unsigned pos[2] = { 0, 0 };
+ /// Curve times on the curves with indices pos[] up until which the
+ /// curves are known to overlap ahead of the nodes.
+ Coord curve_times[2] = { 0.0, 0.0 };
+
+ bool leg = 0; ///< Flag indicating which leg is stepping on the ladder
+ bool just_changed_legs = false;
+
+ /* The ladder algorithm takes steps along the two paths, as if they the stiles of
+ * an imaginary ladder. Note that the nodes (X) on boths paths may not coincide:
+ *
+ * paths[0] START--------X-----------X-----------------------X---------X----> ...
+ * paths[1] START--------------X-----------------X-----------X--------------> ...
+ *
+ * The variables pos[0], pos[1] are the indices of the nodes we've cleared so far;
+ * i.e., we know that the portions before pos[] overlap.
+ *
+ * In each iteration of the loop, we move to the next node along one of the paths;
+ * the variable `leg` tells us which path. We find the point nearest to that node
+ * on the first unprocessed curve of the other path and check the curve time.
+ *
+ * Suppose the current node positions are denoted by P; one possible location of
+ * the nearest point (N) to the next node is:
+ *
+ * ----P------------------N--X---- paths[!leg]
+ * --------P--------------X------- paths[leg] (we've stepped forward from P to X)
+ *
+ * We detect this situation when we find that the curve time of N is < 1.0.
+ * We then create a trimmed version of the top curve so that it corresponds to
+ * the current bottom curve:
+ *
+ * ----P----------------------N--X---------- paths[!leg]
+ * [------------------] trimmed curve
+ * --------P------------------X------------- paths[leg]
+ *
+ * Using isNear(), we can compare the trimmed curve to the front curve (P--X) on
+ * paths[leg]; if they are indeed near, then pos[leg] can be incremented.
+ *
+ * Another possibility is that we overstep the end of the other curve:
+ *
+ * ----P-----------------X------------------ paths[!leg]
+ * N
+ * --------P------------------X------------- paths[leg]
+ *
+ * so the nearest point N now coincides with a node on the top path. We detect
+ * this situation by observing that the curve time of N is close to 1. In case
+ * of such overstep, we change legs by flipping the `leg` variable:
+ *
+ * ----P-----------------X------------------ paths[leg]
+ * --------P------------------X------------- paths[!leg]
+ *
+ * We can now continue the stepping procedure, but the next step will be taken on
+ * the path `paths[leg]`, so it should be a shorter step (if it isn't, the paths
+ * must have diverged and we're done):
+ *
+ * ----P-----------------X------------------ paths[leg]
+ * --------P-------------N----X------------- paths[!leg]
+ *
+ * Another piece of data we hold on to are the curve times on the current curves
+ * up until which the paths have been found to coincide. In other words, at every
+ * step of the algorithm we know that the curves agree up to the path-times
+ * PathTime(pos[i], curve_times[i]).
+ *
+ * In the situation mentioned just above, the times (T) will be as follows:
+ *
+ * ----P---T-------------X------------------ paths[leg]
+ *
+ * --------P-------------N----X------------- paths[!leg]
+ * T
+ *
+ * In this example, the time on top path is > 0.0, since the T mark is further
+ * ahead than P on that path. This value of the curve time is needed to correctly
+ * crop the top curve for the purpose of the isNear() comparison:
+ *
+ * ----P---T-------------X---------- paths[leg]
+ * [-------------] comparison curve (cropped from paths[leg])
+ * [-------------] comparison curve (cropped from paths[!leg])
+ * --------P-------------N----X----- paths[!leg]
+ * T
+ *
+ * In fact, the lower end of the curve time range for cropping is always
+ * given by curve_times[i].
+ *
+ * The iteration ends when we find that the two paths have diverged or when we
+ * reach the end. When that happens, the positions and curve times will be
+ * the PathTime components of the actual point of divergence on both paths.
+ */
+
+ /// A closure to crop and compare the curve pieces ([----] in the diagrams above).
+ auto const pieces_agree = [&](Coord time_on_other) -> bool {
+ Curve *pieces[2];
+ // The leg-side curve is always cropped to the end:
+ pieces[ leg] = paths[ leg]->at(pos[ leg]).portion(curve_times[ leg], 1.0);
+ // The other one is cropped to a variable curve time:
+ pieces[!leg] = paths[!leg]->at(pos[!leg]).portion(curve_times[!leg], time_on_other);
+ bool ret = pieces[0]->isNear(*pieces[1], precision);
+ delete pieces[0];
+ delete pieces[1];
+ return ret;
+ };
+
+ /// A closure to skip degenerate curves; returns true if we reached the end.
+ auto const skip_degenerates = [&](size_t which) -> bool {
+ while (paths[which]->at(pos[which]).isDegenerate()) {
+ ++pos[which];
+ curve_times[which] = 0.0;
+ if (pos[which] == curve_count[which]) {
+ return true; // We've reached the end
+ }
+ }
+ return false;
+ };
+
+ // Main loop of the ladder algorithm.
+ while (pos[0] < curve_count[0] && pos[1] < curve_count[1]) {
+ // Skip degenerate curves if any.
+ if (skip_degenerates(0)) {
+ break;
+ }
+ if (skip_degenerates(1)) {
+ break;
+ }
+
+ // Try to step to the next node with the current leg and see what happens.
+ Coord forward_coord = (Coord)(pos[leg] + 1);
+ if (forward_coord > max_time[leg]) {
+ forward_coord = max_time[leg];
+ }
+ auto const step_point = paths[leg]->pointAt(forward_coord);
+ auto const time_on_other = paths[!leg]->at(pos[!leg]).nearestTime(step_point);
+
+ if (are_near(time_on_other, 1.0, precision) &&
+ are_near(step_point, paths[!leg]->pointAt(pos[!leg] + 1), precision))
+ { // The step took us very near to the first uncertified node on the other path.
+ just_changed_legs = false;
+ //
+ // -------PT-----------------X---------- paths[!leg]
+ // --P-----T-----------------X---------- paths[leg]
+ // ^
+ // endpoints (almost) coincide
+ //
+ // We should compare the curves cropped to the end:
+ //
+ // --------T-----------------X---------- paths[!leg]
+ // [-----------------]
+ // [-----------------]
+ // --------T-----------------X---------- paths[leg]
+ if (pieces_agree(1.0)) {
+ // The curves are nearly identical, so we advance both positions
+ // and zero out the forward curve times.
+ for (size_t i = 0; i < 2; i++) {
+ pos[i]++;
+ curve_times[i] = 0.0;
+ }
+ } else { // We've diverged.
+ break;
+ }
+ } else if (time_on_other < 1.0 - precision) {
+ just_changed_legs = false;
+
+ // The other curve is longer than our step! We trim the other curve to the point
+ // nearest to the step point and compare the resulting pieces.
+ //
+ // --------T-----------------N--------X---- paths[!leg]
+ // [-----------------]
+ // [-----------------]
+ // --------T-----------------X------------- paths[leg]
+ //
+ if (pieces_agree(time_on_other)) { // The curve pieces are near to one another!
+ // We can advance our position and zero out the curve time:
+ pos[leg]++;
+ curve_times[leg] = 0.0;
+ // But on the other path, we can only advance the time, not the curve index:
+ curve_times[!leg] = time_on_other;
+ } else { // We've diverged.
+ break;
+ }
+ } else {
+ // The other curve is shorter than ours, which means that we've overstepped.
+ // We change legs and try to take a shorter step in the next iteration.
+ if (just_changed_legs) {
+ // We already changed legs before and it didn't help, i.e., we've diverged.
+ break;
+ } else {
+ leg = !leg;
+ just_changed_legs = true;
+ }
+ }
+ }
+
+ // Compute the parting time on both paths
+ PathTime path_times[2];
+ for (size_t i = 0; i < 2; i++) {
+ path_times[i] = (pos[i] == curve_count[i]) ? PathTime(curve_count[i] - 1, 1.0)
+ : PathTime(pos[i], curve_times[i]);
+ }
+
+ // Get the parting point from the numerically nicest source
+ Point parting_pt;
+ if (curve_times[0] == 0.0) {
+ parting_pt = paths[0]->pointAt(path_times[0]);
+ } else if (curve_times[1] == 0.0) {
+ parting_pt = paths[1]->pointAt(path_times[1]);
+ } else {
+ parting_pt = middle_point(first.pointAt(path_times[0]), second.pointAt(path_times[1]));
+ }
+
+ return PathIntersection(path_times[0], path_times[1], std::move(parting_pt));
+}
+
+} // namespace Geom
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/path-extrema.cpp b/src/2geom/path-extrema.cpp
new file mode 100644
index 0000000..319ccec
--- /dev/null
+++ b/src/2geom/path-extrema.cpp
@@ -0,0 +1,156 @@
+/** @file Implementation of Path::extrema()
+ */
+/* An algorithm to find the points on a path where a given coordinate
+ * attains its minimum and maximum values.
+ *
+ * Authors:
+ * Rafał Siejakowski <rs@rs-math.net>
+ *
+ * Copyright 2022 the Authors.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#include <2geom/curve.h>
+#include <2geom/path.h>
+#include <2geom/point.h>
+
+namespace Geom {
+
+/** Returns +1 for positive numbers, -1 for negative numbers, and 0 otherwise. */
+inline static float sign(double number)
+{
+ if (number > 0.0) {
+ return 1.0;
+ } else if (number < 0.0) {
+ return -1.0;
+ }
+ return 0.0;
+}
+
+/** @brief Determine whether the d-coordinate increases or decreases at the given path time.
+ *
+ * @param path A path.
+ * @param time A forward-normalized time on the given path.
+ * @param d The coordinate about which we want to know whether it increases.
+ * @return +1.0 if the coordinate increases, -1.0 if it decreases, 0.0 if it remains constant.
+*/
+static float find_direction_of_travel(Path const &path, PathTime const &time, Dim2 d)
+{
+ if (time.t == 0.0) { // We're at a node point
+ if (time.curve_index == 0) { // Starting point of the path.
+ if (path.closed()) {
+ return sign(path.initialUnitTangent()[d] + path.finalUnitTangent()[d]);
+ } else {
+ return sign(path.initialUnitTangent()[d]);
+ }
+ } else if (time.curve_index == path.size()) { // End point of the path.
+ if (path.closed()) {
+ return sign(path.initialUnitTangent()[d] + path.finalUnitTangent()[d]);
+ } else {
+ return sign(path.finalUnitTangent()[d]);
+ }
+ }
+
+ // Otherwise, check the average of the two unit tangent vectors.
+ auto const outgoing_tangent = path.curveAt(time.curve_index).unitTangentAt(0.0);
+ auto const incoming_tangent = path.curveAt(time.curve_index - 1).unitTangentAt(1.0);
+ return sign(outgoing_tangent[d] + incoming_tangent[d]);
+ }
+ // We're in the middle of a curve
+ return sign(path.curveAt(time.curve_index).unitTangentAt(time.t)[d]);
+}
+
+/* Find information about the points on the path where the specified
+ * coordinate attains its minimum and maximum values.
+ */
+PathExtrema Path::extrema(Dim2 d) const
+{
+ auto const ZERO_TIME = PathTime(0, 0);
+
+ // Handle the trivial case of an empty path.
+ if (empty()) {
+ auto const origin = initialPoint();
+ return PathExtrema{
+ .min_point = origin,
+ .max_point = origin,
+ .glance_direction_at_min = 0.0,
+ .glance_direction_at_max = 0.0,
+ .min_time = ZERO_TIME,
+ .max_time = ZERO_TIME
+ };
+ }
+
+ // Set up the simultaneous min-max search
+ Point min_point = initialPoint(), max_point = min_point;
+ auto min_time = ZERO_TIME, max_time = ZERO_TIME;
+ unsigned curve_counter = 0;
+
+ /// A closure to update the current minimum and maximum values.
+ auto const update_minmax = [&](Point const &new_point, Coord t) {
+ if (new_point[d] < min_point[d]) {
+ min_point = new_point;
+ min_time = PathTime(curve_counter, t);
+ } else if (new_point[d] > max_point[d]) {
+ max_point = new_point;
+ max_time = PathTime(curve_counter, t);
+ }
+ };
+
+ // Iterate through the curves, searching for min and max.
+ for (auto const &curve : _data->curves) {
+ // Check the starting node first
+ update_minmax(curve.initialPoint(), 0.0);
+
+ // Check the critical points (zeroes of the derivative).
+ std::unique_ptr<Curve> const derivative{curve.derivative()};
+ for (auto root : derivative->roots(0.0, d)) {
+ update_minmax(curve.pointAt(root), root);
+ }
+ curve_counter++;
+ }
+
+ auto const other = other_dimension(d);
+ return PathExtrema{
+ .min_point = min_point,
+ .max_point = max_point,
+ .glance_direction_at_min = find_direction_of_travel(*this, min_time, other),
+ .glance_direction_at_max = find_direction_of_travel(*this, max_time, other),
+ .min_time = min_time,
+ .max_time = max_time
+ };
+}
+
+} // namespace Geom
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/path.cpp b/src/2geom/path.cpp
index 0bffb01..aeff503 100644
--- a/src/2geom/path.cpp
+++ b/src/2geom/path.cpp
@@ -63,17 +63,33 @@ PathInterval::PathInterval(PathTime const &from, PathTime const &to, bool cross_
, _to(to)
, _path_size(path_size)
, _cross_start(cross_start)
- , _reverse(cross_start ? to >= from : to < from)
+ , _reverse((to < from) ^ cross_start)
{
if (_reverse) {
_to.normalizeForward(_path_size);
+ if (cross_start && _to < to) {
+ // Normalization made us cross start (closed path),
+ // so we don't need to cross the start anymore.
+ _cross_start = false;
+ }
if (_from != _to) {
_from.normalizeBackward(_path_size);
+ if (cross_start && _from > from) {
+ // Normalization backwards made us logically cross
+ // the start – we shouldn't cross the start again.
+ _cross_start = false;
+ }
}
} else {
_from.normalizeForward(_path_size);
+ if (cross_start && _from < from) {
+ _cross_start = false;
+ }
if (_from != _to) {
_to.normalizeBackward(_path_size);
+ if (cross_start && _to > to) {
+ _cross_start = false;
+ }
}
}
@@ -960,8 +976,24 @@ void Path::snapEnds(Coord precision)
}
}
-// replace curves between first and last with contents of source,
-//
+Path Path::withoutDegenerateCurves() const
+{
+ Sequence cleaned;
+ cleaned.reserve(size());
+
+ for (auto it = begin(); it != end_open(); ++it) {
+ if (!it->isDegenerate()) {
+ cleaned.push_back(it->duplicate());
+ }
+ }
+
+ Path result;
+ result._closed = _closed;
+ result.do_update(result._data->curves.begin(), result._data->curves.end(), cleaned);
+ return result;
+}
+
+// Replace curves between first and last with the contents of source.
void Path::do_update(Sequence::iterator first, Sequence::iterator last, Sequence &source)
{
// TODO: handle cases where first > last in closed paths?
diff --git a/src/2geom/planar-graph.h b/src/2geom/planar-graph.h
new file mode 100644
index 0000000..a9c11de
--- /dev/null
+++ b/src/2geom/planar-graph.h
@@ -0,0 +1,1251 @@
+/** @file PlanarGraph – a graph geometrically embedded in the plane.
+ */
+/*
+ * Authors:
+ * Rafał Siejakowski <rs@rs-math.net>
+ *
+ * Copyright 2022 the Authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+// WARNING: This is a private header. Do not include it directly.
+
+#ifndef LIB2GEOM_SEEN_PLANAR_GRAPH_H
+#define LIB2GEOM_SEEN_PLANAR_GRAPH_H
+
+#include <algorithm>
+#include <iterator>
+#include <list>
+
+#include <2geom/angle.h>
+#include <2geom/coord.h>
+#include <2geom/line.h>
+#include <2geom/point.h>
+#include <2geom/path.h>
+#include <2geom/path-intersection.h>
+#include <2geom/utils.h>
+
+namespace Geom {
+
+/**
+ * \class PlanarGraph
+ * \brief Planar graph - a graph geometrically embedded in the plane.
+ *
+ * A planar graph is composed of vertices with assigned locations (as points in the plane)
+ * and of edges (arcs), which are imagined as non-intersecting paths in the plane connecting
+ * the vertices. The edges can hold user-supplied labels (e.g., weights) which support event
+ * callbacks for when the graph is reconfigured, allowing the labels to be updated accordingly.
+ *
+ * \tparam EdgeLabel A user-supplied type; an object of this type will be attached to each
+ * edge of the planar graph (e.g., a "weight" of the edge). The type must
+ * satisfy requirements described further below.
+ *
+ * In order to construct a planar graph, you should specify the clumping precision (passed as
+ * a constructor argument) and then use the method insertEdge() to add edges to the graph, as
+ * many times as necessary. The graph will automatically figure out the locations of the
+ * vertices based on the endpoints of the inserted edges. Vertices will be combined into one
+ * when they are positioned within the distance specified as the clumping threshold, and the
+ * inserted edges will be attached to them accordingly. It is also possible to insert paths
+ * (typically, closed) not attached to any vertices, using the method insertDetached().
+ *
+ * After the edges are inserted, the graph is in a potentially degenerate state, where several
+ * edges may exactly coincide in part or in full. If this is not desired, you can regularize
+ * the graph by calling regularize(). During the regularization process, any overlapping edges
+ * are combined into one. Partially overlapping edges are first split into overlapping and
+ * non-overlapping portions, after which the overlapping portions are combined. If the edges
+ * or their parts overlap but run in opposite directions, one of them will be reversed before
+ * being merged with the other one. The overlaps are detected using the precision setting passed
+ * as the clumping precision in the constructor argument.
+ *
+ * Note however that the regularization procedure does NOT detect transverse intersections
+ * between the edge paths: if such intersections are not desired, the user must pass non-\
+ * intersecting paths to the insertEdge() method (the paths may still have common endpoints,
+ * which is fine: that's how common vertices are created).
+ *
+ * The insertion of new edges invalidates the regularized status, which you can check at any
+ * time by calling isRegularized().
+ *
+ * The vertices stored by the graph are sorted by increasing X-coordinate, and if they have
+ * equal X-coordinates, by increasing Y-coordinate. Even before regularization, incidences of
+ * edges to each vertex are sorted by increasing azimuth of the outgoing tangent (departure
+ * heading, but in radians, in the interval \f$(-\pi, \pi]\f$). After regularization, the edges
+ * around each vertex are guaranteed to be sorted counterclockwise (when the Y-axis points up)
+ * by where they end up going eventually, even if they're tangent at the vertex and therefore
+ * have equal or nearly equal departure azimuths.
+ *
+ * \note
+ * Requirements on the \c EdgeLabel template parameter type.
+ * In order for the template to instantiate correctly, the following must be satisfied:
+ * \li The \c EdgeLabel type provides a method \c onReverse() which gets called whenever
+ * the orientation of the labeled edge is reversed. This is useful when implementing
+ * a directed graph, since the label can keep track of the logical direction.
+ * \li The \c EdgeLabel type provides a method \c onMergeWith(EdgeLabel const&) which gets
+ * called when the labeled edge is combined with a geometrically identical (coinciding)
+ * edge (both combined edges having the same orientations). The label of the edge merged
+ * with the current one is provided as an argument to the method. This is useful when
+ * implementing a graph with weights: for example, when two edges are merged, you may
+ * want to combine their weights in some way.
+ * \li There is a method \c onDetach() called when the edge is removed from the graph. The
+ * edge objects are never destroyed but may be disconnected from the graph when they're no
+ * longer needed; this allows the user to put the labels of such edges in a "dead" state.
+ * \li The \c EdgeLabel objects must be copy-constructible and copy-assignable. This is
+ * because when an edge is subdivided into two, the new edges replacing it get decorated
+ * with copies of the original edge's label.
+ */
+template<typename EdgeLabel>
+#if __cplusplus >= 202002L
+requires requires(EdgeLabel el, EdgeLabel const &other) {
+ el.onReverse();
+ el.onMergeWith(other);
+ el.onDetach();
+ el = other;
+}
+#endif
+class PlanarGraph
+{
+public:
+
+ /** Represents the joint between an edge and a vertex. */
+ struct Incidence
+ {
+ using Sign = bool;
+ inline static Sign const START = false;
+ inline static Sign const END = true;
+
+ Angle azimuth; ///< Angle of the edge's departure.
+ unsigned index; ///< Index of the edge in the parent graph.
+ Sign sign; ///< Whether this is the start or end of the edge.
+ bool invalid = false; ///< Whether this incidence has been marked for deletion.
+
+ Incidence(unsigned edge_index, double departure_azimuth, Sign which_end)
+ : azimuth{departure_azimuth}
+ , index{edge_index}
+ , sign{which_end}
+ {
+ }
+ ~Incidence() = default;
+
+ /// Compare incidences based on their azimuths in radians.
+ inline bool operator<(Incidence const &other) const { return azimuth < other.azimuth; }
+
+ /// Compare the azimuth of an incidence with the given angle.
+ inline bool operator<(double angle) const { return azimuth < angle; }
+
+ /// Check equality (only tests edges and their ends).
+ inline bool operator==(Incidence const &other) const
+ {
+ return index == other.index && sign == other.sign;
+ }
+ };
+ using IncIt = typename std::list<Incidence>::iterator;
+ using IncConstIt = typename std::list<Incidence>::const_iterator;
+
+ /** Represents the vertex of a planar graph. */
+ class Vertex
+ {
+ private:
+ Point const _position; ///< Geometric position of the vertex.
+ std::list<Incidence> _incidences; ///< List of incidences of edges to this vertex.
+ unsigned mutable _flags = 0; ///< User-settable flags.
+
+ inline static Point::LexLess<X> const _cmp; ///< Point sorting function object.
+
+ public:
+ Vertex(Point const &pos)
+ : _position{pos}
+ {
+ }
+
+ /** Get the geometric position of the vertex. */
+ Point const &point() const { return _position; }
+
+ /** Get the list of incidences to the vertex. */
+ auto const &getIncidences() const { return _incidences; }
+
+ /** Compare vertices based on their coordinates (lexicographically). */
+ bool operator<(Vertex const &other) const { return _cmp(_position, other._position); }
+
+ unsigned flags() const { return _flags; } ///< Get the user flags.
+ void setFlags(unsigned flags) const { _flags = flags; } ///< Set the user flags.
+
+ /** Get the cyclic-next incidence after the passed one, in the CCW direction. */
+ IncConstIt cyclicNextIncidence(IncConstIt it) const { return cyclic_next(it, _incidences); }
+
+ /** Get the cyclic-next incidence after the passed one, in the CW direction. */
+ IncConstIt cyclicPrevIncidence(IncConstIt it) const { return cyclic_prior(it, _incidences); }
+
+ /** The same but with pointers. */
+ Incidence *cyclicNextIncidence(Incidence *from)
+ {
+ return &(*cyclicNextIncidence(_incidencePtr2It(from)));
+ }
+ Incidence *cyclicPrevIncidence(Incidence *from)
+ {
+ return &(*cyclicPrevIncidence(_incidencePtr2It(from)));
+ }
+
+ private:
+ /** Same as above, but not const (so only for private use). */
+ IncIt cyclicNextIncidence(IncIt it) { return cyclic_next(it, _incidences); }
+ IncIt cyclicPrevIncidence(IncIt it) { return cyclic_prior(it, _incidences); }
+
+ /** Insert an incidence; for internal use by the PlanarGraph class. */
+ Incidence &_addIncidence(unsigned edge_index, double azimuth, typename Incidence::Sign sign)
+ {
+ auto where = std::find_if(_incidences.begin(), _incidences.end(), [=](auto &inc) -> bool {
+ return inc.azimuth >= azimuth;
+ });
+ return *(_incidences.emplace(where, edge_index, azimuth, sign));
+ }
+
+ /** Return a valid iterator to an incidence passed by pointer;
+ * if the pointer is invalid, return a start iterator. */
+ IncIt _incidencePtr2It(Incidence *pointer)
+ {
+ auto it = std::find_if(_incidences.begin(), _incidences.end(),
+ [=](Incidence const &i) -> bool { return &i == pointer; });
+ return (it == _incidences.end()) ? _incidences.begin() : it;
+ }
+
+ friend class PlanarGraph<EdgeLabel>;
+ };
+ using VertexIterator = typename std::list<Vertex>::iterator;
+
+ /** Represents an edge of the planar graph. */
+ struct Edge
+ {
+ Vertex *start = nullptr, *end = nullptr; ///< Start and end vertices.
+ Path path; ///< The path associated to the edge.
+ bool detached = false; ///< Whether the edge is detached from the graph.
+ bool inserted_as_detached = false; ///< Whether the edge was inserted as detached.
+ EdgeLabel mutable label; ///< The user-supplied label of the edge.
+
+ /** Construct an edge with a given label. */
+ Edge(Path &&movein_path, EdgeLabel &&movein_label)
+ : path{movein_path}
+ , label{movein_label}
+ {
+ }
+
+ /** Detach the edge from the graph. */
+ void detach()
+ {
+ detached = true;
+ label.onDetach();
+ }
+ };
+ using EdgeIterator = typename std::vector<Edge>::iterator;
+ using EdgeConstIterator = typename std::vector<Edge>::const_iterator;
+
+private:
+ double const _precision; ///< Numerical epsilon for vertex clumping.
+ std::list<Vertex> _vertices; ///< Vertices of the graph.
+ std::vector<Edge> _edges; ///< Edges of the graph.
+ std::vector< std::pair<Vertex *, Incidence *> > _junk; ///< Incidences that should be purged.
+ bool _regularized = true; // An empty graph is (trivially) regularized.
+
+public:
+ PlanarGraph(Coord precision = EPSILON)
+ : _precision{precision}
+ {
+ }
+
+ std::list<Vertex> const &getVertices() const { return _vertices; }
+ std::vector<Edge> const &getEdges() const { return _edges; }
+ Edge const &getEdge(size_t index) const { return _edges.at(index); }
+ size_t getEdgeIndex(Edge const &edge) const { return &edge - _edges.data(); }
+ double getPrecision() const { return _precision; }
+ size_t numVertices() const { return _vertices.size(); }
+ size_t numEdges(bool include_detached = true) const
+ {
+ if (include_detached) {
+ return _edges.size();
+ }
+ return std::count_if(_edges.begin(), _edges.end(),
+ [](Edge const &e) -> bool { return !e.detached; });
+ }
+
+ /** Check if the graph has been regularized. */
+ bool isRegularized() const { return _regularized; }
+
+ // 0x1p-50 is about twice the distance between M_PI and the next representable double.
+ void regularize(double angle_precision = 0x1p-50, bool remove_collapsed_loops = true);
+
+ /** Allocate memory to store the specified number of edges. */
+ void reserveEdgeCapacity(size_t capacity) { _edges.reserve(capacity); }
+
+ unsigned insertEdge(Path &&path, EdgeLabel &&edge = EdgeLabel());
+ unsigned insertDetached(Path &&path, EdgeLabel &&edge = EdgeLabel());
+
+ /** Edge insertion with a copy of the path. */
+ unsigned insertEdge(Path const &path, EdgeLabel &&edge = EdgeLabel())
+ {
+ return insertEdge(Path(path), std::forward<EdgeLabel>(edge));
+ }
+ unsigned insertDetached(Path const &path, EdgeLabel &&edge = EdgeLabel())
+ {
+ return insertDetached(Path(path), std::forward<EdgeLabel>(edge));
+ }
+
+ /** \brief Find the incidence at the specified endpoint of the edge.
+ *
+ * \param edge_index The index of the edge whose incidence we wish to query.
+ * \param sign Which end of the edge do we want an incidence of.
+ * \return A pair consisting of pointers to the vertex and the incidence.
+ * If not found, both pointers will be null.
+ */
+ std::pair<Vertex *, Incidence *>
+ getIncidence(unsigned edge_index, typename Incidence::Sign sign) const
+ {
+ if (edge_index >= _edges.size() || _edges[edge_index].detached) {
+ return {nullptr, nullptr};
+ }
+ Vertex *vertex = (sign == Incidence::START) ? _edges[edge_index].start
+ : _edges[edge_index].end;
+ if (!vertex) {
+ return {nullptr, nullptr};
+ }
+ auto it = std::find(vertex->_incidences.begin(), vertex->_incidences.end(),
+ Incidence(edge_index, 42, sign)); // azimuth doesn't matter.
+ if (it == vertex->_incidences.end()) {
+ return {nullptr, nullptr};
+ }
+ return {vertex, &(*it)};
+ }
+
+ /**
+ * \brief Go clockwise or counterclockwise around the vertex and find the next incidence.
+ * The notions of "clockwise"/"counterclockwise" correspond to the y-axis pointing up.
+ *
+ * \param vertex The vertex around which to orbit.
+ * \param incidence The incidence from which to start traversal.
+ * \param clockwise Whether to go clockwise instead of (default) counterclockwise.
+ * \return The next incidence encountered going in the specified direction.
+ */
+ inline Incidence const &nextIncidence(VertexIterator const &vertex, IncConstIt const &incidence,
+ bool clockwise = false) const
+ {
+ return clockwise ? *(vertex->_cyclicPrevIncidence(incidence))
+ : *(vertex->_cyclicNextIncidence(incidence));
+ }
+
+ /** As above, but taking references instead of iterators. */
+ inline Incidence const &nextIncidence(Vertex const &vertex, Incidence const &incidence,
+ bool clockwise = false) const
+ {
+ IncConstIt it = std::find(vertex._incidences.begin(), vertex._incidences.end(), incidence);
+ if (it == vertex._incidences.end()) {
+ return incidence;
+ }
+ return clockwise ? *(vertex.cyclicPrevIncidence(it))
+ : *(vertex.cyclicNextIncidence(it));
+ }
+
+ /** As above, but return an iterator to a const incidence. */
+ inline IncConstIt nextIncidenceIt(Vertex const &vertex, Incidence const &incidence,
+ bool clockwise = false) const
+ {
+ IncConstIt it = std::find(vertex._incidences.begin(), vertex._incidences.end(), incidence);
+ if (it == vertex._incidences.end()) {
+ return vertex._incidences.begin();
+ }
+ return clockwise ? vertex.cyclicPrevIncidence(it)
+ : vertex.cyclicNextIncidence(it);
+ }
+ inline IncConstIt nextIncidenceIt(Vertex const &vertex, IncConstIt const &incidence,
+ bool clockwise = false) const
+ {
+ return clockwise ? vertex.cyclicPrevIncidence(incidence)
+ : vertex.cyclicNextIncidence(incidence);
+ }
+
+ /** As above, but start at the prescribed departure azimuth (in radians).
+ *
+ * \return A pointer to the incidence emanating from the vertex at or immediately after
+ * the specified azimuth, when going around the vertex in the specified direction.
+ * If the vertex has no incidences, return value is nullptr.
+ */
+ Incidence *nextIncidence(VertexIterator const &vertex, double azimuth,
+ bool clockwise = false) const;
+
+ /** Get the incident path, always oriented away from the vertex. */
+ Path getOutgoingPath(Incidence const *incidence) const
+ {
+ return incidence ? _getPathImpl(incidence, Incidence::START) : Path();
+ }
+
+ /** Get the incident path, always oriented towards the vertex. */
+ Path getIncomingPath(Incidence const *incidence) const
+ {
+ return incidence ? _getPathImpl(incidence, Incidence::END) : Path();
+ }
+
+private:
+ inline Path _getPathImpl(Incidence const *incidence, typename Incidence::Sign origin) const
+ {
+ return (incidence->sign == origin) ? _edges[incidence->index].path
+ : _edges[incidence->index].path.reversed();
+ }
+
+ /** Earmark an incidence for future deletion. */
+ inline void _throwAway(Vertex *vertex, Incidence *incidence)
+ {
+ if (!vertex || !incidence) {
+ return;
+ }
+ incidence->invalid = true;
+ _junk.emplace_back(vertex, incidence);
+ }
+
+ // Topological reconfiguration functions; see their definitions for documentation.
+ bool _compareAndReglue(Vertex &vertex, Incidence *first, Incidence *second, bool deloop);
+ Vertex *_ensureVertexAt(Point const &position);
+ void _mergeCoincidingEdges(Incidence *first, Incidence *second);
+ void _mergeShorterLonger(Vertex &vertex, Incidence *shorter, Incidence *longer,
+ PathTime const &time_on_longer);
+ void _mergeWyeConfiguration(Vertex &vertex, Incidence *first, Incidence *second,
+ PathIntersection const &split);
+ void _purgeJunkIncidences();
+ void _reglueLasso(Vertex &vertex, Incidence *first, Incidence *second,
+ PathIntersection const &split);
+ bool _reglueTeardrop(Vertex &vertex, Incidence *first, Incidence *second, bool deloop);
+ void _reglueTangentFan(Vertex &vertex, IncIt const &first, IncIt const &last, bool deloop);
+ void _regularizeVertex(Vertex &vertex, double angle_precision, bool deloop);
+
+ // === Static stuff ===
+
+ /** Return the angle between the vector and the positive X axis or 0 if undefined. */
+ inline static double _getAzimuth(Point const &vec) { return vec.isZero() ? 0.0 : atan2(vec); }
+
+ /** Return path time corresponding to the same point on the reversed path. */
+ inline static PathTime _reversePathTime(PathTime const &time, Path const &path)
+ {
+ int new_index = path.size() - time.curve_index - 1;
+ Coord new_time = 1.0 - time.t;
+ if (new_index < 0) {
+ new_index = 0;
+ new_time = 0;
+ }
+ return PathTime(new_index, new_time);
+ }
+
+ /** Return path time at the end of the path. */
+ inline static PathTime _pathEnd(Path const &path) { return PathTime(path.size() - 1, 1.0); }
+ inline static auto const PATH_START = PathTime(0, 0);
+
+public:
+ static double closedPathArea(Path const &path);
+ static bool deviatesLeft(Path const &first, Path const &second);
+};
+
+/**
+ * \brief Insert a new vertex or reuse an existing one.
+ *
+ * Ensures that there is a vertex at or near the specified position
+ * (within the distance of _precision).
+ *
+ * \param pos The desired geometric position of the new vertex.
+ * \return A pointer to the inserted vertex or a pre-existing vertex near the
+ * desired position.
+ */
+template<typename EL>
+typename PlanarGraph<EL>::Vertex *PlanarGraph<EL>::_ensureVertexAt(Point const &pos)
+{
+ auto const insert_at_front = [&, this]() -> Vertex* {
+ _vertices.emplace_front(pos);
+ return &(_vertices.front());
+ };
+
+ if (_vertices.empty()) {
+ return insert_at_front();
+ }
+
+ // TODO: Use a heap?
+ auto it = std::find_if(_vertices.begin(), _vertices.end(), [&](Vertex const &v) -> bool {
+ return Vertex::_cmp(pos, v._position); // existing vertex position > pos.
+ });
+
+ if (it != _vertices.end()) {
+ if (are_near(pos, it->_position, _precision)) {
+ return &(*it); // Reuse existing vertex.
+ }
+ if (it == _vertices.begin()) {
+ return insert_at_front();
+ }
+ }
+ // Look at the previous element, reuse if near, insert before `it` otherwise.
+ return &(*(are_near(pos, std::prev(it)->_position, _precision) ? std::prev(it)
+ : _vertices.emplace(it, pos)));
+}
+
+/**
+ * \brief Move-insert a new labeled edge into the planar graph.
+ *
+ * \param path The geometric path of the edge.
+ * \param label Optionally, the label (extra user data) associated to this edge.
+ * If absent, a default-constructed label will be used.
+ * \return The index of the inserted edge.
+ */
+template<typename EdgeLabel>
+unsigned PlanarGraph<EdgeLabel>::insertEdge(Path &&path, EdgeLabel &&label)
+{
+ unsigned edge_index = _edges.size();
+ auto &inserted = _edges.emplace_back(std::forward<Path>(path),
+ std::forward<EdgeLabel>(label));
+
+ // Calculate the outgoing azimuths at both endpoints.
+ double const start_azimuth = _getAzimuth(inserted.path.initialUnitTangent());
+ double const end_azimuth = _getAzimuth(-inserted.path.finalUnitTangent());
+
+ // Get the endpoints into the graph.
+ auto start = _ensureVertexAt(inserted.path.initialPoint());
+ auto end = _ensureVertexAt(inserted.path.finalPoint());
+
+ // Inform the edge about its endpoints.
+ inserted.start = start;
+ inserted.end = end;
+
+ // Add incidences at the endpoints.
+ start->_addIncidence(edge_index, start_azimuth, Incidence::START);
+ end->_addIncidence(edge_index, end_azimuth, Incidence::END);
+
+ _regularized = false;
+ return edge_index;
+}
+
+/**
+ * \brief Move-insert a new labeled edge but do not connect it to the graph.
+ *
+ * Although the graph will hold the path data of an edge inserted in this way, the edge
+ * will not be connected to any vertex. This can be used to store information about closed
+ * paths (loops) in the instance, without having to specify starting points for them.
+ *
+ * \param path The geometric path of the edge.
+ * \param label Optionally, the label (extra user data) associated to this edge; if absent,
+ * the label will be default-constructed.
+ * \return The index of the inserted edge.
+ */
+template<typename EdgeLabel>
+unsigned PlanarGraph<EdgeLabel>::insertDetached(Path &&path, EdgeLabel &&label)
+{
+ unsigned edge_index = _edges.size();
+ auto &inserted = _edges.emplace_back(std::forward<Path>(path),
+ std::forward<EdgeLabel>(label));
+ inserted.detached = true;
+ inserted.inserted_as_detached = true;
+ return edge_index;
+}
+
+/** Remove incidences previously marked as junk. */
+template<typename EdgeLabel>
+void PlanarGraph<EdgeLabel>::_purgeJunkIncidences()
+{
+ for (auto &[vertex, incidence] : _junk) {
+ Incidence *to_remove = incidence;
+ auto it = std::find_if(vertex->_incidences.begin(), vertex->_incidences.end(),
+ [=](Incidence const &inc) -> bool { return &inc == to_remove; });
+ if (it != vertex->_incidences.end()) {
+ vertex->_incidences.erase(it);
+ }
+ }
+ _junk.clear();
+}
+
+/**
+ * \brief Merge overlapping edges or their portions, adding vertices if necessary.
+ *
+ * \param angle_precision The numerical epsilon for radian angle comparisons.
+ * \param remove_collapsed_loops Whether to detach edges with both ends incident to the same
+ * vertex (loops) when these loops don't enclose any area.
+ *
+ * This function performs the following operations:
+ * \li Edges that are tangent at a vertex but don't otherwise overlap are sorted correctly
+ * in the counterclockwise cyclic order around the vertex.
+ * \li Degenerate loops which don't enclose any area are removed if the argument is true.
+ * \li Edges that coincide completely are reversed if needed and merged into one.
+ * \li Edges that coincide partially are split into overlapping and non-overlapping portions.
+ * Any overlapping portions are oriented consistently and then merged.
+ * \li As a sub-case of the above, any non-degenerate loop with an initial self-everlap
+ * (a "lasso") is replaced with a shorter non-overlapping loop and a simple path leading
+ * to it.
+ */
+template<typename EdgeLabel>
+void PlanarGraph<EdgeLabel>::regularize(double angle_precision, bool remove_collapsed_loops)
+{
+ for (auto it = _vertices.begin(); it != _vertices.end(); ++it) {
+ // Note: the list of vertices may grow during the execution of this loop,
+ // so don't replace it with a range-for (which stores the end iterator).
+ // We want the loop to carry on going over the elements it inserted.
+ if (it->_incidences.size() < 2) {
+ continue;
+ }
+ _regularizeVertex(*it, angle_precision, remove_collapsed_loops);
+ }
+ _purgeJunkIncidences();
+ _regularized = true;
+}
+
+/**
+ * \brief Analyze and regularize all edges emanating from a given vertex.
+ *
+ * This function goes through the list of incidences at the vertex (roughly sorted by
+ * azimuth, i.e., departure heading in radians), picking out runs of mutually tangent
+ * edges and calling _reglueTangentFan() on each run. The algorithm is quite complicated
+ * because the incidences have to be treated as a cyclic container and a run of mutually
+ * tangent edges may straddle the "end" of the list, including the possibility that the
+ * entire list is a single such run.
+ *
+ * \param vertex The vertex whose incidences should be analyzed.
+ * \param angle_precision The numerical epsilon for radian angle comparisons.
+ * \param deloop Whether loops that don't enclose any area should be detached.
+ */
+template<typename EdgeLabel>
+void PlanarGraph<EdgeLabel>::_regularizeVertex(typename PlanarGraph<EdgeLabel>::Vertex &vertex,
+ double angle_precision, bool deloop)
+{
+ auto &incidences = vertex._incidences;
+
+ /// Compare two polar angles in the interval [-π, π] modulo 2π to within angle_precision:
+ auto const angles_equal = [=](double az1, double az2) -> bool {
+ static double const twopi = 2.0 * M_PI;
+ return are_near(std::fmod(az1 + twopi, twopi), std::fmod(az2 + twopi, twopi),
+ angle_precision);
+ };
+
+ IncIt run_begin; // First element in the last discovered run of equal azimuths.
+
+ /// Find and reglue runs of nearly identical azimuths in the specified range.
+ auto const process_runs = [&](IncIt begin, IncIt end) -> bool
+ {
+ double current_azimuth = 42; // Invalid radian angle.
+ bool in_a_run = false;
+
+ for (auto it = begin; it != end; ++it) {
+ bool const equal = angles_equal(it->azimuth, current_azimuth);
+ if (equal && !in_a_run) {
+ run_begin = std::prev(it); // Save to enclosing scope.
+ in_a_run = true;
+ } else if (!equal && in_a_run) {
+ _reglueTangentFan(vertex, run_begin, std::prev(it), deloop);
+ in_a_run = false;
+ }
+ current_azimuth = it->azimuth;
+ }
+ return in_a_run;
+ };
+
+ double const last_azimuth = incidences.back().azimuth;
+
+ if (angles_equal(incidences.front().azimuth, last_azimuth)) {
+ // The cyclic list contains a run of equal azimuths which straddles the "end".
+ // This means that we must skip the part of this run on the "begin" side on the
+ // first pass and handle it once we've traversed the remainder of the list.
+
+ bool processed = false; ///< Whether we've cleared the straddling run.
+ double previous_azimuth = last_azimuth;
+ IncIt straddling_run_last;
+
+ for (auto it = incidences.begin(); it != incidences.end(); ++it) {
+ if (!angles_equal(it->azimuth, previous_azimuth)) {
+ straddling_run_last = std::prev(it);
+ process_runs(it, incidences.end());
+ processed = true;
+ break;
+ }
+ previous_azimuth = it->azimuth;
+ }
+ if (processed) {
+ // Find the first element of the straddling run.
+ auto it = std::prev(incidences.end());
+ while (angles_equal(it->azimuth, last_azimuth)) {
+ --it;
+ }
+ ++it; // Now we're at the start of the straddling run.
+ _reglueTangentFan(vertex, it, straddling_run_last, deloop);
+ } else {
+ // We never encountered anything outside of the straddling run: reglue everything.
+ _reglueTangentFan(vertex, incidences.begin(), std::prev(incidences.end()), deloop);
+ }
+ } else if (process_runs(incidences.begin(), incidences.end())) {
+ // Our run got rudely interrupted by the end of the container; reglue till the end.
+ _reglueTangentFan(vertex, run_begin, std::prev(incidences.end()), deloop);
+ }
+}
+
+/**
+ * \brief Regularize a fan of mutually tangent edges emanating from a vertex.
+ *
+ * This function compares the tangent edges pairwise and ensures that the sequence of their
+ * incidences to the vertex ends up being sorted by the ultimate direction in which the
+ * emanating edges fan out, in the counterclockwise order.
+ *
+ * If a partial or complete overlap between edges is detected, these edges are reglued.
+ *
+ * \param vertex The vertex from which the fan emanates.
+ * \param first An iterator pointing to the first incidence in the fan.
+ * \param last An iterator pointing to the last incidence in the fan.
+ * NOTE: This iterator must point to the actual last incidence, not "past" it.
+ * The reason is that we're iterating over a cyclic collection, so there
+ * isn't really a meaningful end.
+ * \param deloop Whether loops that don't enclose any area should be detached.
+ */
+template<typename EL>
+void PlanarGraph<EL>::_reglueTangentFan(typename PlanarGraph<EL>::Vertex &vertex,
+ typename PlanarGraph<EL>::IncIt const &first,
+ typename PlanarGraph<EL>::IncIt const &last, bool deloop)
+{
+ // Search all pairs (triangular pattern), skipping invalid incidences.
+ for (auto it = first; it != last; it = vertex.cyclicNextIncidence(it)) {
+ if (it->invalid) {
+ continue;
+ }
+ for (auto is = vertex.cyclicNextIncidence(it); true; is = vertex.cyclicNextIncidence(is)) {
+ if (!is->invalid && _compareAndReglue(vertex, &(*it), &(*is), deloop)) {
+ // Swap the incidences, effectively implementing "bubble sort".
+ std::swap(*it, *is);
+ }
+ if (is == last) {
+ break;
+ }
+ }
+ }
+}
+
+/**
+ * \brief Compare a pair of edges emanating from the same vertex in the same direction.
+ *
+ * If the edges overlap in part or in full, they get reglued, which means that the topology
+ * of the graph may get modified. Otherwise, if the detailed comparison shows that the edges
+ * aren't correctly ordered around the vertex (because the second edge deviates to the right
+ * instead of to the left of the first, when looking away from the vertex), then the function
+ * will return true, signalling that the incidences should be swapped.
+ *
+ * \param vertex The vertex where the mutually tangent paths meet.
+ * \param first The incidence appearing as the first one in the provisional cyclic order.
+ * \param second The incidence appearing as the second one in the provisional cyclic order.
+ * \param deloop Whether to detach collapsed loops (backtracks) which don't enclose any area.
+ * \return Whether the incidences should be swapped.
+ */
+template<typename EL>
+bool PlanarGraph<EL>::_compareAndReglue(typename PlanarGraph<EL>::Vertex &vertex,
+ typename PlanarGraph<EL>::Incidence *first,
+ typename PlanarGraph<EL>::Incidence *second, bool deloop)
+{
+ if (first->index == second->index) {
+ return _reglueTeardrop(vertex, first, second, deloop);
+ }
+
+ // Get paths corresponding to the edges but travelling away from the vertex.
+ auto first_path_out = getOutgoingPath(first);
+ auto second_path_out = getOutgoingPath(second);
+ auto split = parting_point(first_path_out, second_path_out, _precision);
+
+ if (are_near(split.point(), vertex.point(), _precision)) {
+ // Paths deviate immediately, so no gluing is needed. The incidences should
+ // be swapped if the first edge path departs to the left of the second one.
+ return deviatesLeft(first_path_out, second_path_out);
+ }
+
+ // Determine the nature of the initial overlap between the paths.
+ bool const till_end_of_1st = are_near(split.point(), first_path_out.finalPoint(), _precision);
+ bool const till_end_of_2nd = are_near(split.point(), second_path_out.finalPoint(), _precision);
+
+ if (till_end_of_1st && till_end_of_2nd) { // Paths coincide completely.
+ _mergeCoincidingEdges(first, second);
+ } else if (till_end_of_1st) {
+ // Paths coincide until the end of the 1st one, which however isn't the end of the
+ // 2nd one; for example, the first one could be the vertical riser of the letter L
+ // whereas the second one – the entire letter stroke.
+ _mergeShorterLonger(vertex, first, second, split.second);
+ } else if (till_end_of_2nd) {
+ // The same but with with the second edge shorter than the first one.
+ _mergeShorterLonger(vertex, second, first, split.first);
+ } else { // A Y-shaped split.
+ _mergeWyeConfiguration(vertex, first, second, split);
+ }
+ return false; // We've glued so no need to swap anything.
+}
+
+/**
+ * \brief Analyze a loop path a with self-tangency at the attachment point (a teardrop).
+ *
+ * The following steps are taken:
+ * \li If the loop encloses zero area and \c deloop is true, the loop is detached.
+ * \li If the two arms of the loop split out immediately, the loop is left alone and we
+ * only check whether the incidences should be swapped.
+ * \li If the loop overlaps itself near the vertex, resembling a lasso, we split it into
+ * a shorter simple path and a smaller loop attached to the end of the shorter path.
+ *
+ * \param vertex The vertex at which the teardrop originates.
+ * \param first The first incidence of the loop to the vertex.
+ * \param second The second incidence of the loop to the vertex.
+ * \param deloop Whether the loop should be removed if it doesn't enclose any area
+ * (i.e., the path exactly backtracks on itself).
+ * \return Whether the two incidences of the loop to the vertex should be swapped.
+ */
+template<typename EL>
+bool PlanarGraph<EL>::_reglueTeardrop(typename PlanarGraph<EL>::Vertex &vertex,
+ typename PlanarGraph<EL>::Incidence *first,
+ typename PlanarGraph<EL>::Incidence *second, bool deloop)
+{
+ // Calculate the area enclosed by the teardrop.
+ // The convention is that the unit circle (cos(t), sint(t)), t from 0 to 2pi,
+ // encloses an area of +pi.
+ auto &edge = _edges[first->index];
+ Path loop = edge.path; loop.close();
+ double signed_area = closedPathArea(loop);
+
+ if (deloop && are_near(signed_area, 0.0, _precision)) {
+ edge.detach();
+ _throwAway(&vertex, first);
+ _throwAway(&vertex, second);
+ return false;
+ }
+
+ auto split = parting_point(loop, loop.reversed(), _precision);
+ if (are_near(split.point(), vertex.point(), _precision)) {
+ // The loop spreads out immediately. We simply check if the incidences should be swapped.
+ // We want them to be ordered such that the signed area encircled by the path going out
+ // at the first incidence and coming back at the second (with this orientation) is > 0.
+ return (first->sign == Incidence::START) ^ (signed_area > 0.0);
+ }
+
+ // The loop encloses a nonzero area, but the two branches don't separate at the starting
+ // point. Instead, they travel together for a while before they split like a lasso.
+ _reglueLasso(vertex, first, second, split);
+ return false;
+}
+
+/**
+ * \brief Reglue a lasso-shaped loop, separating it into the "free rope" and the "hoop".
+ *
+ * The lasso is an edge looping back to the same vertex, where the closed path encloses
+ * a non-zero area, but its two branches don't separate at the starting point. Instead,
+ * they travel together for a while (forming the doubled-up "free rope") before they
+ * split like a lasso. This function cuts the lasso at the split point:
+ * \code{.unparsed}
+ * ____ ____
+ * / \ / \
+ * VERTEX =====< | ==> VERTEX ------ NEW + NEW < |
+ * \____/ (lasso) (rope) \____/ (hoop)
+ *
+ * \endcode
+ *
+ * \param vertex A reference to the vertex where the lasso is attached.
+ * \param first The first incidence of the lasso to the vertex.
+ * \param second The second incidence of the lasso to the vertex.
+ * \param split The point where the free rope of the lasso ends and the hoop begins.
+ */
+template<typename EL>
+void PlanarGraph<EL>::_reglueLasso(typename PlanarGraph<EL>::Vertex &vertex,
+ typename PlanarGraph<EL>::Incidence *first,
+ typename PlanarGraph<EL>::Incidence *second,
+ PathIntersection const &split)
+{
+ unsigned lasso = first->index;
+
+ // Create the "free rope" path.
+ auto rope = _edges[lasso].path.portion(PATH_START, split.first);
+ rope.setInitial(vertex.point());
+ rope.setFinal(split.point());
+ double const rope_final_backward_azimuth = _getAzimuth(-rope.finalUnitTangent());
+
+ // Compute the new label of the rope edge.
+ auto oriented_as_loop = _edges[lasso].label;
+ auto reversed = oriented_as_loop; reversed.onReverse();
+ oriented_as_loop.onMergeWith(reversed);
+
+ // Insert the rope and its endpoint.
+ unsigned const rope_index = _edges.size();
+ auto &rope_edge = _edges.emplace_back(std::move(rope), std::move(oriented_as_loop));
+ auto const new_split_vertex = _ensureVertexAt(split.point());
+
+ // Reuse lasso's first incidence as the incidence to the rope (azimuth can stay).
+ first->index = rope_index;
+ first->sign = Incidence::START;
+
+ // Connect the rope to the newly created split vertex.
+ new_split_vertex->_addIncidence(rope_index, rope_final_backward_azimuth, Incidence::END);
+ rope_edge.start = &vertex;
+ rope_edge.end = new_split_vertex;
+
+ // Insert the hoop
+ auto hoop = _edges[lasso].path.portion(split.first,
+ _reversePathTime(split.second, _edges[lasso].path));
+ hoop.setInitial(split.point());
+ hoop.setFinal(split.point());
+ insertEdge(std::move(hoop), EL(_edges[lasso].label));
+
+ // Detach the original lasso edge and mark the second incidence for cleanup.
+ _edges[lasso].detach();
+ _throwAway(&vertex, second);
+}
+
+/**
+ * \brief Completely coallesce two fully overlapping edges.
+ *
+ * In practice, the first edge stays and the second one gets detached from the graph.
+ *
+ * \param first An iterator to the first edge's incidence to a common vertex.
+ * \param second An iterator to the second edge's incidence to a common vertex.
+ */
+template<typename EL>
+void PlanarGraph<EL>::_mergeCoincidingEdges(typename PlanarGraph<EL>::Incidence *first,
+ typename PlanarGraph<EL>::Incidence *second)
+{
+ auto &surviver = _edges[first->index];
+ auto &casualty = _edges[second->index];
+
+ auto other_label = casualty.label;
+ if (first->sign != second->sign) { // Logically reverse the label before merging.
+ other_label.onReverse();
+ }
+ surviver.label.onMergeWith(other_label);
+
+ // Mark both incidences of the second edge as junk and detach it.
+ auto [start_vertex, start_inc] = getIncidence(second->index, Incidence::START);
+ _throwAway(start_vertex, start_inc);
+ auto [end_vertex, end_inc] = getIncidence(second->index, Incidence::END);
+ _throwAway(end_vertex, end_inc);
+ casualty.detach();
+}
+
+/**
+ * \brief Merge a longer edge with a shorter edge that overlaps it.
+ *
+ * In practice, the shorter edge remains unchanged and the longer one is trimmed to
+ * become just the part extending past the shorter one.
+ *
+ * \param vertex The vertex where the overlap starts.
+ * \param shorter The incidence of the shorter edge to the common vertex.
+ * \param longer The incidence of the longer edge to the common vertex.
+ * \param time_on_longer The PathTime on the longer edge at which it passes through
+ * the endpoint of the shorter edge.
+ */
+template<typename EL>
+void PlanarGraph<EL>::_mergeShorterLonger(typename PlanarGraph<EL>::Vertex &vertex,
+ typename PlanarGraph<EL>::Incidence *shorter,
+ typename PlanarGraph<EL>::Incidence *longer,
+ PathTime const &time_on_longer)
+{
+ auto &shorter_edge = _edges[shorter->index];
+ auto &longer_edge = _edges[longer->index];
+
+ // Get the vertices at the far ends of both edges.
+ auto shorter_far_end = (shorter->sign == Incidence::START) ? shorter_edge.end
+ : shorter_edge.start;
+ /// Whether the longer edge heads out of the vertex.
+ bool const longer_out = (longer->sign == Incidence::START);
+ auto longer_far_end = longer_out ? longer_edge.end : longer_edge.start;
+
+ // Copy the longer edge's label and merge with that of the shorter.
+ auto longer_label = longer_edge.label;
+ if (shorter->sign != longer->sign) {
+ longer_label.onReverse();
+ }
+ shorter_edge.label.onMergeWith(longer_label);
+
+ // Create the trimmed path (longer minus shorter).
+ Path trimmed;
+ double trimmed_departure_azimuth;
+ if (longer_out) {
+ trimmed = longer_edge.path.portion(time_on_longer, _pathEnd(longer_edge.path));
+ longer_edge.start = shorter_far_end;
+ trimmed.setInitial(shorter_far_end->point());
+ trimmed.setFinal(longer_far_end->point());
+ trimmed_departure_azimuth = _getAzimuth(trimmed.initialUnitTangent());
+ } else {
+ trimmed = longer_edge.path.portion(PATH_START, _reversePathTime(time_on_longer,
+ longer_edge.path));
+ longer_edge.end = shorter_far_end;
+ trimmed.setInitial(longer_far_end->point());
+ trimmed.setFinal(shorter_far_end->point());
+ trimmed_departure_azimuth = _getAzimuth(-trimmed.finalUnitTangent());
+ }
+
+ // Set the trimmed path as the new path of the longer edge and set up the incidences:
+ longer_edge.path = std::move(trimmed);
+ shorter_far_end->_addIncidence(longer->index, trimmed_departure_azimuth, longer->sign);
+
+ // Throw away the old start incidence of the longer edge.
+ _throwAway(&vertex, longer);
+}
+
+/**
+ * \brief Merge a pair of partially overlapping edges, producing a Y-split at a new vertex.
+ *
+ * This topological modification is performed by inserting a new vertex at the three-way
+ * point (where the two paths separate) and clipping the original edges to that point.
+ * In this way, the original edges become the "arms" of the Y-shape. In addition, a new
+ * edge is inserted, forming the "stem" of the Y.
+ *
+ * \param vertex The vertex from which the partially overlapping edges originate (bottom of Y).
+ * \param first The incidence to the first edge (whose path is the stem and one arm of the Y).
+ * \param second The incidence to the second edge (stem and the other arm of the Y).
+ * \param fork The splitting point of the two paths.
+ */
+template<typename EL>
+void PlanarGraph<EL>::_mergeWyeConfiguration(typename PlanarGraph<EL>::Vertex &vertex,
+ typename PlanarGraph<EL>::Incidence *first,
+ typename PlanarGraph<EL>::Incidence *second,
+ PathIntersection const &fork)
+{
+ bool const first_is_out = (first->sign == Incidence::START);
+ bool const second_is_out = (second->sign == Incidence::START);
+
+ auto &first_edge = _edges[first->index];
+ auto &second_edge = _edges[second->index];
+
+ // Calculate the path forming the stem of the Y:
+ auto stem_path = getOutgoingPath(first).portion(PATH_START, fork.first);
+ stem_path.setInitial(vertex.point());
+ stem_path.setFinal(fork.point());
+
+ /// A closure to clip the path of an original edge to the fork point.
+ auto const clip_to_fork = [&](PathTime const &t, Edge &e, bool out) {
+ if (out) { // Trim from time to end
+ e.path = e.path.portion(t, _pathEnd(e.path));
+ e.path.setInitial(fork.point());
+ } else { // Trim from reverse-end to reverse-time
+ e.path = e.path.portion(PATH_START, _reversePathTime(t, e.path));
+ e.path.setFinal(fork.point());
+ }
+ };
+
+ /// A closure to find the departing azimuth of an edge at the fork point.
+ auto const departing_azimuth = [&](Edge const &e, bool out) -> double {
+ return _getAzimuth((out) ? e.path.initialUnitTangent()
+ : -e.path.finalUnitTangent());
+ };
+
+ // Clip the paths obtaining the arms of the Y.
+ clip_to_fork(fork.first, first_edge, first_is_out);
+ clip_to_fork(fork.second, second_edge, second_is_out);
+
+ // Create the fork vertex and set up its incidences.
+ auto const fork_vertex = _ensureVertexAt(fork.point());
+ fork_vertex->_addIncidence(first->index, departing_azimuth(first_edge, first_is_out),
+ first->sign);
+ fork_vertex->_addIncidence(second->index, departing_azimuth(second_edge, second_is_out),
+ second->sign);
+
+ // Repoint the ends of the edges that were clipped
+ (first_is_out ? first_edge.start : first_edge.end) = fork_vertex;
+ (second_is_out ? second_edge.start : second_edge.end) = fork_vertex;
+
+ /// A closure to get a consistently oriented label of an edge.
+ auto upwards_oriented_label = [&](Edge const &e, bool out) -> EL {
+ auto label = e.label;
+ if (!out) {
+ label.onReverse();
+ }
+ return label;
+ };
+
+ auto stem_label = upwards_oriented_label(first_edge, first_is_out);
+ stem_label.onMergeWith(upwards_oriented_label(second_edge, second_is_out));
+ auto stem_departure_from_fork = _getAzimuth(-stem_path.finalUnitTangent());
+
+ // Insert the stem of the Y-configuration.
+ unsigned const stem_index = _edges.size();
+ auto stem_edge = _edges.emplace_back(std::move(stem_path), std::move(stem_label));
+ stem_edge.start = &vertex;
+ stem_edge.end = fork_vertex;
+
+ // Set up the incidences.
+ fork_vertex->_addIncidence(stem_index, stem_departure_from_fork, Incidence::END);
+ first->index = stem_index;
+ first->sign = Incidence::START;
+ _throwAway(&vertex, second);
+}
+
+template<typename EL>
+typename PlanarGraph<EL>::Incidence*
+PlanarGraph<EL>::nextIncidence(typename PlanarGraph<EL>::VertexIterator const &vertex,
+ double azimuth, bool clockwise) const
+{
+ auto &incidences = vertex._incidences;
+ Incidence *result = nullptr;
+
+ if (incidences.empty()) {
+ return result;
+ }
+ // Normalize azimuth to the interval [-pi; pi].
+ auto angle = Angle(azimuth);
+
+ if (clockwise) { // Go backwards and find a lower bound
+ auto it = std::find_if(incidences.rbegin(), incidences.rend(), [=](auto inc) -> bool {
+ return inc.azimuth <= angle;
+ });
+ if (it == incidences.rend()) {
+ // azimuth is lower than the azimuths of all incidences;
+ // going clockwise we wrap back to the highest azimuth (last incidence).
+ return &incidences.back();
+ }
+ result = &(*it);
+ } else {
+ auto it = std::find_if(incidences.begin(), incidences.end, [=](auto inc) -> bool {
+ return inc.azimuth >= angle;
+ });
+ if (it == incidences.end()) {
+ // azimuth is higher than the azimuths of all incidences;
+ // going counterclockwise we wrap back to the lowest azimuth.
+ return &incidences.front();
+ }
+ result = &(*it);
+ }
+ return result;
+}
+
+/** Return the signed area enclosed by a closed path. */
+template<typename EL>
+double PlanarGraph<EL>::closedPathArea(Path const &path)
+{
+ double area;
+ Point _;
+ centroid(path.toPwSb(), _, area);
+ return -area; // Our convention is that the Y-axis points up
+}
+
+/** \brief Determine whether the first path deviates to the left of the second.
+ *
+ * The two paths are assumed to have identical or nearly identical starting points
+ * but not an overlapping initial portion. The concept of "left" is based on the
+ * y-axis pointing up.
+ *
+ * \param first The first path.
+ * \param second The second path.
+ *
+ * \return True if the first path deviates towards the left of the second;
+ * False if the first path deviates towards the right of the second.
+ */
+template<typename EL>
+bool PlanarGraph<EL>::deviatesLeft(Path const &first, Path const &second)
+{
+ auto start = middle_point(first.initialPoint(), second.initialPoint());
+ auto tangent_between = middle_point(first.initialUnitTangent(), second.initialUnitTangent());
+ if (tangent_between.isZero()) {
+ return false;
+ }
+ auto tangent_line = Line::from_origin_and_vector(start, tangent_between);
+
+ // Find the first non-degenerate curves on both paths
+ std::unique_ptr<Curve> c[2];
+ auto const find_first_nondegen = [](std::unique_ptr<Curve> &pointer, Path const &path) {
+ for (auto const &c : path) {
+ if (!c.isDegenerate()) {
+ pointer.reset(c.duplicate());
+ return;
+ }
+ }
+ };
+
+ find_first_nondegen(c[0], first);
+ find_first_nondegen(c[1], second);
+ if (!c[0] || !c[1]) {
+ return false;
+ }
+
+ // Find the bounding boxes
+ Rect const bounding_boxes[] {
+ c[0]->boundsExact(),
+ c[1]->boundsExact()
+ };
+
+ // For a bounding box, find the corner that goes the furthest in the direction of the
+ // tangent vector.
+ auto const furthest_corner = [&](Rect const &r) -> unsigned {
+ Coord max_dot = dot(r.corner(0) - start, tangent_between);
+ unsigned result = 0;
+ for (unsigned i = 1; i < 4; i++) {
+ auto current_dot = dot(r.corner(i), tangent_between);
+ if (current_dot > max_dot) {
+ max_dot = current_dot;
+ result = i;
+ } else if (current_dot == max_dot) {
+ // Disambiguate based on proximity to the tangent line.
+ auto const offset = start + tangent_between;
+ if (distance(offset, r.corner(i)) < distance(offset, r.corner(result))) {
+ result = i;
+ }
+ }
+ }
+ return result;
+ };
+
+ // Calculate the corner points overlooking the "rift" between the paths.
+ Point corner_points[2];
+ for (size_t i : {0, 1}) {
+ corner_points[i] = bounding_boxes[i].corner(furthest_corner(bounding_boxes[i]));
+ }
+
+ // Find a vantage point from which we can best observe the splitting paths.
+ Point vantage_point;
+ bool found = false;
+ if (corner_points[0] != corner_points[1]) {
+ auto line_connecting_corners = Line(corner_points[0], corner_points[1]);
+ auto xing = line_connecting_corners.intersect(tangent_line);
+ if (!xing.empty()) {
+ vantage_point = xing[0].point();
+ found = true;
+ }
+ }
+ if (!found) {
+ vantage_point = tangent_line.pointAt(tangent_line.timeAtProjection(corner_points[0]));
+ }
+
+ // Move to twice as far in the direction of the vantage point.
+ vantage_point += vantage_point - start;
+
+ // Find the points on both curves that are nearest to the vantage point.
+ Coord nearest[2];
+ for (size_t i : {0, 1}) {
+ nearest[i] = c[i]->nearestTime(vantage_point);
+ }
+
+ // Clip to the nearest points and examine the closed contour.
+ Path closed_contour(start);
+ closed_contour.setStitching(true);
+ closed_contour.append(c[0]->portion(0, nearest[0]));
+ closed_contour = closed_contour.reversed();
+ closed_contour.append(c[1]->portion(0, nearest[1]));
+ closed_contour.close();
+ return !path_direction(closed_contour); // Reverse to match the convention that y-axis is up.
+}
+
+} // namespace Geom
+
+#endif // LIB2GEOM_SEEN_PLANAR_GRAPH_H
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/2geom/self-intersect.cpp b/src/2geom/self-intersect.cpp
new file mode 100644
index 0000000..4fe4d9e
--- /dev/null
+++ b/src/2geom/self-intersect.cpp
@@ -0,0 +1,313 @@
+/**
+ * @file Implementation of Path::intersectSelf() and PathVector::intersectSelf().
+ */
+/* An algorithm for finding self-intersections of paths and path-vectors.
+ *
+ * Authors:
+ * Rafał Siejakowski <rs@rs-math.net>
+ *
+ * (C) Copyright 2022 Authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#include <list>
+
+#include <2geom/coord.h>
+#include <2geom/curve.h>
+#include <2geom/path.h>
+#include <2geom/pathvector.h>
+#include <2geom/point.h>
+#include <2geom/sweeper.h>
+
+namespace Geom {
+
+/** @brief The PathSelfIntersector class is a sweepset class used for intersecting curves in the
+ * same path with one another. It is intended to be used as the template parameter of Sweeper.
+ */
+class PathSelfIntersector
+{
+public:
+ using ItemIterator = Path::iterator;
+
+private:
+ Path _path; ///< The path searched for self-crossings, cleaned of degenerate curves.
+ std::list<ItemIterator> _active; ///< List of active curves during the sweepline passage.
+ std::vector<PathIntersection> _crossings; ///< Stores the crossings found.
+ std::vector<size_t> _original_indices; ///< Curve indices before removal of degenerate curves.
+ double const _precision; ///< Numerical epsilon.
+
+public:
+ PathSelfIntersector(Path const &path, double precision)
+ : _path{path.initialPoint()}
+ , _precision{precision}
+ {
+ _original_indices.reserve(path.size());
+ for (size_t i = 0; i < path.size(); i++) {
+ if (!path[i].isDegenerate()) {
+ _path.append(path[i]);
+ _original_indices.push_back(i);
+ }
+ }
+ _path.close(path.closed());
+ }
+
+ // === SweepSet API ===
+ auto &items() { return _path; }
+ Interval itemBounds(ItemIterator curve) const { return curve->boundsFast()[X]; }
+ /// Callback for when the sweepline starts intersecting a new item.
+ void addActiveItem(ItemIterator incoming)
+ {
+ _intersectWithActive(incoming);
+ _intersectWithSelf(incoming);
+ _active.push_back(incoming);
+ }
+ /// Callback for when the sweepline stops intersecting an item.
+ void removeActiveItem(ItemIterator to_remove)
+ {
+ auto it = std::find(_active.begin(), _active.end(), to_remove);
+ _active.erase(it);
+ }
+ // ===
+
+ std::vector<PathIntersection> &&moveOutCrossings() { return std::move(_crossings); }
+
+private:
+ /** Find and store all intersections of a curve with itself. */
+ void _intersectWithSelf(ItemIterator curve)
+ {
+ size_t const index = std::distance(_path.begin(), curve);
+ for (auto &&self_x : curve->intersectSelf(_precision)) {
+ _appendCurveCrossing(std::move(self_x), index, index);
+ }
+ }
+
+ /** Find and store all intersections of a curve with the active curves. */
+ void _intersectWithActive(ItemIterator curve)
+ {
+ size_t const index = std::distance(_path.begin(), curve);
+ for (auto const &other : _active) {
+ if (!curve->boundsFast().intersects(other->boundsFast())) {
+ continue;
+ }
+
+ size_t const other_index = std::distance(_path.begin(), other);
+ auto const &[smaller, larger] = std::minmax(index, other_index);
+ /// Whether the curves meet at a common node in the path.
+ bool consecutive = smaller + 1 == larger;
+ /// Whether the curves meet at the closure point of the path.
+ bool wraparound = _path.closed() && smaller == 0 && larger + 1 == _path.size();
+ for (auto &&xing : curve->intersect(*other, _precision)) {
+ _appendCurveCrossing(std::move(xing), index, other_index, consecutive, wraparound);
+ }
+ }
+ }
+
+ /** Append a curve crossing to the store as long as it satisfies nondegeneracy criteria. */
+ void _appendCurveCrossing(CurveIntersection &&xing, size_t first_index, size_t second_index,
+ bool consecutive = false, bool wraparound = false)
+ {
+ // Filter out crossings that aren't real but rather represent the agreement of final
+ // and initial points of consecutive curves – a consequence of the path's continuity.
+ auto const should_exclude = [&](bool flipped) -> bool {
+ // Filter out spurious self-intersections by using squared geometric average.
+ bool const first_is_first = (first_index < second_index) ^ flipped;
+ double const geom2 = first_is_first ? (1.0 - xing.first) * xing.second
+ : (1.0 - xing.second) * xing.first;
+ return geom2 < EPSILON;
+ };
+
+ if ((consecutive && should_exclude(false)) || (wraparound && should_exclude(true))) {
+ return;
+ }
+
+ // Convert curve indices to the original ones (before the removal of degenerate curves).
+ _crossings.emplace_back(PathTime(_original_indices[first_index], xing.first),
+ PathTime(_original_indices[second_index], xing.second),
+ xing.point());
+ }
+};
+
+// Compute all crossings of a path with itself.
+std::vector<PathIntersection> Path::intersectSelf(Coord precision) const
+{
+ auto intersector = PathSelfIntersector(*this, precision);
+ Sweeper(intersector).process();
+ auto result = intersector.moveOutCrossings();
+ std::sort(result.begin(), result.end());
+ return result;
+}
+
+/**
+ * @brief The PathVectorSelfIntersector class is an implementation of a SweepSet whose intended
+ * use is the search for self-intersections in a single PathVector. It's designed to be used as
+ * the template parameter for the Sweeper class template.
+ */
+class PathVectorSelfIntersector
+{
+public:
+ using ItemIterator = PathVector::const_iterator;
+
+private:
+ PathVector const &_pathvector; ///< A reference to the path-vector searched for self-crossings.
+ std::list<ItemIterator> _active; ///< A list of active paths during sweepline passage.
+ std::vector<PathVectorIntersection> _crossings; ///< Stores the crossings found.
+ double const _precision; ///< Numerical epsilon.
+
+public:
+ PathVectorSelfIntersector(PathVector const &subject, double precision)
+ : _pathvector{subject}
+ , _precision{precision}
+ {
+ }
+
+ // == SweepSet API ===
+ auto const &items() { return _pathvector; }
+ Interval itemBounds(ItemIterator path)
+ {
+ auto const r = path->boundsFast();
+ return r ? (*r)[X] : Interval(); // Sweeplines are vertical
+ }
+
+ /// Callback for when the sweepline starts intersecting a new item.
+ void addActiveItem(ItemIterator incoming)
+ {
+ _intersectWithActive(incoming);
+ _intersectWithSelf(incoming);
+ _active.push_back(incoming);
+ }
+
+ /// Callback for when the sweepline stops intersecting an item.
+ void removeActiveItem(ItemIterator to_remove)
+ {
+ auto it = std::find(_active.begin(), _active.end(), to_remove);
+ _active.erase(it);
+ }
+ // ===
+
+ std::vector<PathVectorIntersection> &&moveOutCrossings() { return std::move(_crossings); }
+
+private:
+ /**
+ * @brief Find all intersections of the path pointed to by the given
+ * iterator with all currently active paths and store results
+ * in the instance of the class.
+ *
+ * @param it An iterator to a path to be intersected with the active ones.
+ */
+ void _intersectWithActive(ItemIterator &it);
+
+ /**
+ * @brief Find all intersections of the path pointed to by the given
+ * iterator with itself and store the results in the class instance.
+ *
+ * @param it An iterator to a path which will be intersected with itself.
+ */
+ void _intersectWithSelf(ItemIterator &it);
+
+ /// Append a path crossing to the store.
+ void _appendPathCrossing(PathIntersection const &xing, size_t first_index, size_t second_index)
+ {
+ auto const first_time = PathVectorTime(first_index, xing.first);
+ auto const second_time = PathVectorTime(second_index, xing.second);
+ _crossings.emplace_back(first_time, second_time, xing.point());
+ }
+
+public:
+
+ std::vector<PathVectorIntersection>
+ filterDeduplicate(std::vector<PathVectorIntersection> &&xings) const;
+};
+
+/** Remove duplicate intersections (artifacts of the path/curve crossing algorithms). */
+std::vector<PathVectorIntersection>
+PathVectorSelfIntersector::filterDeduplicate(std::vector<PathVectorIntersection> &&xings) const
+{
+ std::vector<PathVectorIntersection> result;
+ result.reserve(xings.size());
+
+ auto const are_same_times = [&](Coord a1, Coord a2, Coord b1, Coord b2) -> bool {
+ return (are_near(a1, b1) && are_near(a2, b2)) ||
+ (are_near(a1, b2) && are_near(a2, b1));
+ };
+
+ Coord last_time_1 = -1.0, last_time_2 = -1.0; // Invalid path times
+ for (auto &&x : xings) {
+ auto const current_1 = x.first.asFlatTime(), current_2 = x.second.asFlatTime();
+ if (!are_same_times(current_1, current_2, last_time_1, last_time_2)) {
+ result.push_back(std::move(x));
+ }
+ last_time_1 = current_1;
+ last_time_2 = current_2;
+ }
+
+ return result;
+}
+
+/** Compute and store intersections of a path with all active paths. */
+void PathVectorSelfIntersector::_intersectWithActive(ItemIterator &it)
+{
+ auto const start = _pathvector.begin();
+ for (auto &path : _active) {
+ if (!path->boundsFast().intersects(it->boundsFast())) {
+ continue;
+ }
+ for (auto &&xing : path->intersect(*it, _precision)) {
+ _appendPathCrossing(std::move(xing), std::distance(start, path),
+ std::distance(start, it));
+ }
+ }
+}
+
+/** Compute and store intersections of a constituent path with itself. */
+void PathVectorSelfIntersector::_intersectWithSelf(ItemIterator &it)
+{
+ size_t const path_index = std::distance(_pathvector.begin(), it);
+ for (auto &&xing : it->intersectSelf(_precision)) {
+ _appendPathCrossing(std::move(xing), path_index, path_index);
+ }
+}
+
+// Compute self-intersections in a path-vector.
+std::vector<PathVectorIntersection> PathVector::intersectSelf(Coord precision) const
+{
+ auto intersector = PathVectorSelfIntersector(*this, precision);
+ Sweeper(intersector).process();
+ auto result = intersector.moveOutCrossings();
+ std::sort(result.begin(), result.end());
+ return (result.size() > 1) ? intersector.filterDeduplicate(std::move(result)) : result;
+}
+
+} // namespace Geom
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/src/toys/CMakeLists.txt b/src/toys/CMakeLists.txt
index 9499739..c4929b7 100644
--- a/src/toys/CMakeLists.txt
+++ b/src/toys/CMakeLists.txt
@@ -15,6 +15,7 @@ SET(2GEOM_TOYS-2_SRC
aa
arc-bez
arc-length-param
+auto-cross
boolops-toy
bound-path
bounds-test
@@ -53,7 +54,7 @@ evolute
filet-minion
find-derivative
gear
-hatches
+#hatches
implicit-toy
ineaa
inner-product-clip
@@ -90,9 +91,9 @@ rect_02
rect_03
rect-toy
root-finder-comparer
-rtree-toy
+#rtree-toy
sanitize
-sb1d
+#sb1d
sb2d
sb2d-solver
sbasisdim
@@ -110,7 +111,7 @@ squiggles
sweep
sweeper-toy
# these ones have only had a trivial rewrite to toy-2
-uncross
+#uncross
winding-test
worms
)
diff --git a/src/toys/auto-cross.cpp b/src/toys/auto-cross.cpp
new file mode 100644
index 0000000..00b5b25
--- /dev/null
+++ b/src/toys/auto-cross.cpp
@@ -0,0 +1,321 @@
+/* @brief
+ * A toy for playing around with Path::intersectSelf().
+ *
+ * Authors:
+ * Rafał Siejakowski <rs@rs-math.net>
+ *
+ * Copyright 2022 the Authors.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#include <toys/toy-framework-2.h>
+#include <2geom/path.h>
+#include <2geom/elliptical-arc.h>
+#include <2geom/svg-path-parser.h>
+
+using namespace Geom;
+using Color = uint32_t;
+
+Color const RED = 0x80000000;
+Color const GREEN = 0x00800000;
+Color const BROWN = 0x90500000;
+Color const BLUE = 0x0000ff00;
+Color const BLACK = 0x00000000;
+
+static void set_cairo_rgb(cairo_t *c, Color rgb)
+{
+ cairo_set_source_rgba(c, (double)((rgb & 0xFF000000) >> 24) / 255.0,
+ (double)((rgb & 0x00FF0000) >> 16) / 255.0,
+ (double)((rgb & 0x0000FF00) >> 8) / 255.0,
+ 1.0);
+}
+
+static void write_text(cairo_t *c, const char *text, Point const &position, Color color)
+{
+ cairo_move_to(c, position);
+ cairo_set_font_size(c, 12);
+ set_cairo_rgb(c, color);
+ cairo_show_text(c, text);
+}
+
+static std::string format_point(Point const &pt)
+{
+ std::ostringstream ss;
+ ss.precision(4);
+ ss << pt;
+ return ss.str();
+}
+
+static EllipticalArc random_arc(Point from, Point to)
+{
+ double const dist = distance(from, to);
+ auto angle = atan2(to - from);
+ bool sweep = std::abs(angle) > M_PI_2;
+ angle *= 2;
+ angle = std::fmod(angle, 2.0 * M_PI);
+ return EllipticalArc(from, Point(0.5 * dist, 2.0 * dist), angle, false, sweep, to);
+}
+
+class Item
+{
+private:
+ Path _path;
+ Color _color;
+ std::string _d;
+
+public:
+ Item(Color color)
+ : _color{color}
+ {}
+
+ void setPath(Path &&new_path)
+ {
+ _path = std::forward<Path>(new_path);
+ std::ostringstream oss;
+ oss << _path;
+ _d = oss.str();
+ }
+
+ void draw(cairo_t *cr) const
+ {
+ cairo_set_line_width(cr, 2);
+ set_cairo_rgb(cr, _color);
+ cairo_path(cr, _path);
+ cairo_stroke(cr);
+ _drawBezierTangents(cr);
+ _drawSelfIntersections(cr);
+ }
+
+ void write(cairo_t *cr, Point const &pos) const
+ {
+ write_text(cr, _d.c_str(), pos, _color);
+ }
+
+ std::string const& getSVGD() const { return _d; }
+
+private:
+ void _drawBezierTangents(cairo_t *c) const
+ {
+ cairo_set_line_width(c, 1);
+ set_cairo_rgb(c, 0x0000b000);
+ // Draw tangents for Beziers:
+ for (auto const &curve : _path) {
+ if (auto const *bezier = dynamic_cast<BezierCurve const *>(&curve)) {
+ if (bezier->order() > 1) {
+ auto points = bezier->controlPoints();
+ cairo_move_to(c, points[0]);
+ cairo_line_to(c, points[1]);
+ cairo_stroke(c);
+ cairo_move_to(c, points.back());
+ cairo_line_to(c, points[points.size() - 2]);
+ cairo_stroke(c);
+ }
+ }
+ }
+ }
+
+ void _drawSelfIntersections(cairo_t *cr) const
+ {
+ set_cairo_rgb(cr, BLACK);
+ for (auto const &xing : _path.intersectSelf()) {
+ draw_cross(cr, xing.point());
+ auto const coords = format_point(xing.point());
+ write_text(cr, coords.c_str(), xing.point() + Point(8, -8), BLACK);
+ }
+ }
+};
+
+class AutoCross : public Toy
+{
+public:
+ AutoCross()
+ : items{Item(RED), Item(GREEN), Item(BROWN)}
+ {
+ bezier_handles.pts = { {200, 400}, {100, 300}, {300, 400}, {300, 300}, {450, 300}, {500, 500}, {400, 400} };
+ elliptical_handles.pts = { {500, 200}, {700, 400}, {600, 500} };
+ mixed_handles.pts = { {100, 600}, {120, 690}, {300, 650}, {330, 600}, {500, 800} };
+ handles.push_back(&bezier_handles);
+ handles.push_back(&elliptical_handles);
+ handles.push_back(&mixed_handles);
+ }
+
+ void draw(cairo_t *cr, std::ostringstream *notify, int width, int height, bool save,
+ std::ostringstream *timer_stream) override
+ {
+ if (crashed) {
+ draw_error(cr, width, height);
+ } else {
+ try {
+ draw_impl(cr, width, height);
+ } catch (Exception &e) {
+ error = e.what();
+ handles.clear();
+ crashed = true;
+ }
+ }
+ Toy::draw(cr, notify, width, height, save, timer_stream);
+ }
+
+ void key_hit(GdkEventKey *ev) override
+ {
+ if (ev->keyval == GDK_KEY_space) {
+ print_path_d();
+ } else if ((ev->keyval == GDK_KEY_V || ev->keyval == GDK_KEY_v) && (ev->state & GDK_CONTROL_MASK)) {
+ paste_d();
+ }
+
+ }
+
+private:
+ std::string error;
+ std::vector<Item> items;
+ PointSetHandle bezier_handles, elliptical_handles, mixed_handles;
+ bool crashed = false;
+
+ void paste_d()
+ {
+ auto *clipboard = gtk_clipboard_get(GDK_SELECTION_CLIPBOARD);
+ auto *text = gtk_clipboard_wait_for_text(clipboard);
+ if (!text) {
+ return;
+ }
+ PathVector pv;
+ try {
+ pv = parse_svg_path(text);
+ } catch (SVGPathParseError &error) {
+ std::cerr << "Error pasting path d: " << error.what() << std::endl;
+ return;
+ }
+ if (pv.empty()) {
+ return;
+ }
+ Item paste_item{RED}; // TODO: cycle through a color palette.
+ paste_item.setPath(std::move(pv[0]));
+ items.push_back(paste_item);
+ redraw();
+ }
+
+ void print_path_d()
+ {
+ std::cout << "Path snapshots:\n";
+ for (auto it = items.rbegin(); it != items.rend(); ++it) {
+ std::cout << it->getSVGD() << '\n';
+ }
+ }
+
+ void refresh_geometry()
+ {
+ // Construct the 2-segment Bézier path
+ auto const &cp = bezier_handles.pts;
+ Path bezier;
+ bezier.append(BezierCurveN<3>(cp[0].round(), cp[1].round(), cp[2].round(), cp[3].round()));
+ bezier.append(BezierCurveN<3>(cp[3].round(), cp[4].round(), cp[5].round(), cp[6].round()));
+ items[0].setPath(std::move(bezier));
+
+ // Construct the elliptical arcs
+ auto const &ae = elliptical_handles.pts;
+ Path elliptical;
+ elliptical.append(random_arc(ae[0], ae[1]));
+ elliptical.append(random_arc(ae[1], ae[2]));
+ items[1].setPath(std::move(elliptical));
+
+ // Construct a mixed path
+ auto const &mh = mixed_handles.pts;
+ Path mixed;
+ mixed.append(BezierCurveN<3>(mh[0], mh[1], mh[2], mh[3]));
+ mixed.append(random_arc(mh[3], mh[4]));
+ mixed.close();
+ items[2].setPath(std::move(mixed));
+ }
+
+ void draw_impl(cairo_t *cr, int width, int height)
+ {
+ refresh_geometry();
+ write_title(cr);
+
+ auto text_pos = Point(20, height - 20);
+ for (auto const &item : items) {
+ item.draw(cr);
+ item.write(cr, text_pos);
+ text_pos -= Point(0, 20);
+ }
+ }
+
+ void write_title(cairo_t *c)
+ {
+ cairo_move_to(c, 10, 40);
+ cairo_set_font_size(c, 30);
+ set_cairo_rgb(c, 0x0);
+ cairo_show_text(c, "Self-intersection of paths in lib2geom!");
+ cairo_set_font_size(c, 14);
+ cairo_move_to(c, 10, 60);
+ cairo_show_text(c, "[Space]: Print SVG 'd' attributes to stdout");
+ cairo_move_to(c, 10, 80);
+ cairo_show_text(c, "[Ctrl-V]: Paste a 'd' attribute from clipboard");
+ }
+
+ void draw_error(cairo_t *cr, int width, int height)
+ {
+ auto center = Point(0.5 * (double)width, 0.5 * (double)height);
+ cairo_move_to(cr, center + Point(-90, -100));
+ cairo_line_to(cr, center + Point(100, 90));
+ cairo_line_to(cr, center + Point(90, 100));
+ cairo_line_to(cr, center + Point(-100, -90));
+ cairo_close_path(cr);
+ cairo_set_source_rgb(cr, 1, 0, 0);
+ cairo_fill(cr);
+
+ cairo_move_to(cr, center + Point(90, -100));
+ cairo_line_to(cr, center + Point(100, -90));
+ cairo_line_to(cr, center + Point(-90, 100));
+ cairo_line_to(cr, center + Point(-100, 90));
+ cairo_close_path(cr);
+ cairo_set_source_rgb(cr, 1, 0, 0);
+ cairo_fill(cr);
+
+ cairo_move_to(cr, center + Point(-90, 120));
+ cairo_show_text(cr, "Sorry, your toy has just broken :-/");
+ cairo_move_to(cr, Point(10, center[Y] + 150));
+ cairo_show_text(cr, error.c_str());
+ }
+};
+
+int main(int argc, char **argv)
+{
+ auto toy = AutoCross();
+ init(argc, argv, &toy, 800, 800);
+ return 0;
+}
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+ */
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
\ No newline at end of file
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index e96f768..95e10d4 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -28,10 +28,12 @@ line-test
nl-vector-test
parallelogram-test
path-test
+planar-graph-test
point-test
polynomial-test
rect-test
sbasis-test
+self-intersections-test
)
foreach(source ${2GEOM_GTESTS_SRC})
diff --git a/tests/bezier-test.cpp b/tests/bezier-test.cpp
index 513352c..507a4d1 100644
--- a/tests/bezier-test.cpp
+++ b/tests/bezier-test.cpp
@@ -6,7 +6,7 @@
* Nathan Hurst <njh@njhurst.com>
* Krzysztof Kosiński <tweenk.pl@gmail.com>
* Johan Engelen <j.b.c.engelen@alumnus.utwente.nl>
- *
+ *
* Copyright 2010 Authors
*
* This library is free software; you can redistribute it and/or
@@ -44,7 +44,7 @@
#include <iterator>
#include <glib.h>
-using namespace std;
+using std::vector, std::min, std::max;
using namespace Geom;
Poly lin_poly(double a, double b) { // ax + b
@@ -57,7 +57,7 @@ Poly lin_poly(double a, double b) { // ax + b
bool are_equal(Bezier A, Bezier B) {
int maxSize = max(A.size(), B.size());
double t = 0., dt = 1./maxSize;
-
+
for(int i = 0; i <= maxSize; i++) {
EXPECT_FLOAT_EQ(A.valueAt(t), B.valueAt(t));// return false;
t += dt;
@@ -85,10 +85,10 @@ protected:
};
TEST_F(BezierTest, Basics) {
-
+
//std::cout << unit <<std::endl;
//std::cout << hump <<std::endl;
-
+
EXPECT_TRUE(Bezier(0,0,0,0).isZero());
EXPECT_TRUE(Bezier(0,1,2,3).isFinite());
@@ -96,13 +96,13 @@ TEST_F(BezierTest, Basics) {
///cout << " Bezier::Bezier(const Bezier& b);\n";
//cout << Bezier(wiggle) << " == " << wiggle << endl;
-
+
//cout << "explicit Bezier(unsigned ord);\n";
//cout << Bezier(10) << endl;
-
+
//cout << "Bezier(Coord c0, Coord c1);\n";
//cout << Bezier(0.0,1.0) << endl;
-
+
//cout << "Bezier(Coord c0, Coord c1, Coord c2);\n";
//cout << Bezier(0,1, 2) << endl;
@@ -121,9 +121,9 @@ TEST_F(BezierTest, ValueAt) {
EXPECT_EQ(3.0, wiggle.at1());
EXPECT_EQ(0.0, wiggle.valueAt(0.5));
-
+
EXPECT_EQ(0.0, wiggle(0.5));
-
+
//cout << "SBasis toSBasis();\n";
//cout << unit.toSBasis() << endl;
//cout << hump.toSBasis() << endl;
@@ -280,7 +280,7 @@ TEST_F(BezierTest, Deflate) {
TEST_F(BezierTest, Roots) {
expect_array((const double[]){0, 0.5, 0.5}, wiggle.roots());
-
+
/*Bezier bigun(Bezier::Order(30));
for(unsigned i = 0; i < bigun.size(); i++) {
bigun.setCoeff(i,rand()-0.5);
@@ -303,7 +303,7 @@ TEST_F(BezierTest, Roots) {
tests.push_back(vector_from_array((const double[]){0, 0.2, 0.6, 0.6, 1}));
tests.push_back(vector_from_array((const double[]){.1,.2,.3,.4,.5,.6}));
tests.push_back(vector_from_array((const double[]){0.25,0.25,0.25,0.75,0.75,0.75}));
-
+
for(auto & test : tests) {
Bezier b = array_roots(test);
//std::cout << tests[test_i] << ": " << b << std::endl;
@@ -327,11 +327,24 @@ TEST_F(BezierTest, BoundsExact) {
}
TEST_F(BezierTest, Operators) {
- /*cout << "scalar operators\n";
- cout << hump + 3 << endl;
- cout << hump - 3 << endl;
- cout << hump*3 << endl;
- cout << hump/3 << endl;*/
+ // Test equality operators
+ EXPECT_EQ(zero, zero);
+ EXPECT_EQ(hump, hump);
+ EXPECT_EQ(wiggle, wiggle);
+ EXPECT_EQ(unit, unit);
+
+ EXPECT_NE(zero, hump);
+ EXPECT_NE(hump, zero);
+ EXPECT_NE(wiggle, hump);
+ EXPECT_NE(zero, wiggle);
+ EXPECT_NE(wiggle, unit);
+
+ // Recall that hump == Bezier(0,1,0);
+ EXPECT_EQ(hump + 3, Bezier(3, 4, 3));
+ EXPECT_EQ(hump - 3, Bezier(-3, -2, -3));
+ EXPECT_EQ(hump * 3, Bezier(0, 3, 0));
+ EXPECT_EQ(hump / 3, Bezier(0, 1.0/3.0, 0));
+ EXPECT_EQ(-hump, Bezier(0, -1, 0));
Bezier reverse_wiggle = reverse(wiggle);
EXPECT_EQ(reverse_wiggle.at0(), wiggle.at1());
@@ -502,6 +515,152 @@ TEST_F(BezierTest, Intersection) {
#endif
}
+/** Basic test for intersecting a quadratic Bézier with a line segment. */
+TEST_F(BezierTest, QuadraticIntersectLineSeg)
+{
+ double const EPS = 1e-12;
+ auto const bow = QuadraticBezier({0, 0}, {1, 1}, {2, 0});
+ auto const highhoriz = LineSegment(Point(0, 0), Point(2, 0));
+ auto const midhoriz = LineSegment(Point(0, 0.25), Point(2, 0.25));
+ auto const lowhoriz = LineSegment(Point(0, 0.5), Point(2, 0.5));
+ auto const noninters = LineSegment(Point(0, 0.5 + EPS), Point(2, 0.5 + EPS));
+ auto const noninters2 = LineSegment(Point(1, 0), Point(1, 0.5 - EPS));
+
+ auto const endpoint_intersections = bow.intersect(highhoriz, EPS);
+ EXPECT_EQ(endpoint_intersections.size(), 2);
+ EXPECT_intersections_valid(bow, highhoriz, endpoint_intersections, EPS);
+ for (auto const &ex : endpoint_intersections) {
+ EXPECT_DOUBLE_EQ(ex.point()[Y], 0.0);
+ }
+
+ auto const mid_intersections = bow.intersect(midhoriz, EPS);
+ EXPECT_EQ(mid_intersections.size(), 2);
+ EXPECT_intersections_valid(bow, midhoriz, mid_intersections, EPS);
+ for (auto const &mx : mid_intersections) {
+ EXPECT_DOUBLE_EQ(mx.point()[Y], 0.25);
+ }
+
+ auto const tangent_intersection = bow.intersect(lowhoriz, EPS);
+ EXPECT_EQ(tangent_intersection.size(), 1);
+ EXPECT_intersections_valid(bow, lowhoriz, tangent_intersection, EPS);
+ for (auto const &tx : tangent_intersection) {
+ EXPECT_DOUBLE_EQ(tx.point()[Y], 0.5);
+ }
+
+ auto no_intersections = bow.intersect(noninters, EPS);
+ EXPECT_TRUE(no_intersections.empty());
+
+ no_intersections = bow.intersect(noninters2, EPS);
+ EXPECT_TRUE(no_intersections.empty());
+}
+
+TEST_F(BezierTest, QuadraticIntersectLineRandom)
+{
+ g_random_set_seed(0xB747A380);
+ auto const diagonal = LineSegment(Point(0, 0), Point(1, 1));
+ double const EPS = 1e-12;
+
+ for (unsigned i = 0; i < 10'000; i++) {
+ auto q = QuadraticBezier({0, 1}, {g_random_double_range(0.0, 1.0), g_random_double_range(0.0, 1.0)}, {1, 0});
+ auto xings = q.intersect(diagonal, EPS);
+ ASSERT_EQ(xings.size(), 1);
+ auto pt = xings[0].point();
+ EXPECT_TRUE(are_near(pt[X], pt[Y], EPS));
+ EXPECT_intersections_valid(q, diagonal, xings, EPS);
+ }
+}
+
+/** Basic test for intersecting a cubic Bézier with a line segment. */
+TEST_F(BezierTest, CubicIntersectLine)
+{
+ double const EPS = 1e-12;
+ auto const wavelet = CubicBezier({0, 0}, {1, 2}, {0, -2}, {1, 0});
+
+ auto const unit_seg = LineSegment(Point(0, 0), Point(1, 0));
+ auto const expect3 = wavelet.intersect(unit_seg, EPS);
+ EXPECT_EQ(expect3.size(), 3);
+ EXPECT_intersections_valid(wavelet, unit_seg, expect3, EPS);
+
+ auto const half_seg = LineSegment(Point(0, 0), Point(0.5, 0));
+ auto const expect2 = wavelet.intersect(half_seg, EPS);
+ EXPECT_EQ(expect2.size(), 2);
+ EXPECT_intersections_valid(wavelet, half_seg, expect2, EPS);
+
+ auto const less_than_half = LineSegment(Point(0, 0), Point(0.5 - EPS, 0));
+ auto const expect1 = wavelet.intersect(less_than_half, EPS);
+ EXPECT_EQ(expect1.size(), 1);
+ EXPECT_intersections_valid(wavelet, less_than_half, expect1, EPS);
+
+ auto const dollar_stroke = LineSegment(Point(0, 0.5), Point(1, -0.5));
+ auto const dollar_xings = wavelet.intersect(dollar_stroke, EPS);
+ EXPECT_EQ(dollar_xings.size(), 3);
+ EXPECT_intersections_valid(wavelet, dollar_stroke, dollar_xings, EPS);
+}
+
+TEST_F(BezierTest, CubicIntersectLineRandom)
+{
+ g_random_set_seed(0xCAFECAFE);
+ auto const diagonal = LineSegment(Point(0, 0), Point(1, 1));
+ double const EPS = 1e-8;
+
+ for (unsigned i = 0; i < 10'000; i++) {
+ double a1 = g_random_double_range(0.0, 1.0);
+ double a2 = g_random_double_range(a1, 1.0);
+ double b1 = g_random_double_range(0.0, 1.0);
+ double b2 = g_random_double_range(0.0, b1);
+
+ auto c = CubicBezier({0, 1}, {a1, a2}, {b1, b2}, {1, 0});
+ auto xings = c.intersect(diagonal, EPS);
+ ASSERT_EQ(xings.size(), 1);
+ auto pt = xings[0].point();
+ EXPECT_TRUE(are_near(pt[X], pt[Y], EPS));
+ EXPECT_intersections_valid(c, diagonal, xings, EPS);
+ }
+}
+
+/** Regression test for issue https://gitlab.com/inkscape/lib2geom/-/issues/47 . */
+TEST_F(BezierTest, Balloon)
+{
+ auto const loop = CubicBezier({0, 0}, {4, -2}, {4, 2}, {0, 0});
+ auto const seghoriz = LineSegment(Point(-1, 0), Point(0, 0));
+
+ for (double EPS : {1e-6, 1e-9, 1e-12}) {
+ // We expect that 2 intersections are found: one at each end of the loop,
+ // both at the coordinates (0, 0).
+ auto xings_horiz = loop.intersect(seghoriz, EPS);
+ EXPECT_EQ(xings_horiz.size(), 2);
+ EXPECT_intersections_valid(loop, seghoriz, xings_horiz, EPS);
+ }
+}
+
+TEST_F(BezierTest, ExpandToTransformedTest)
+{
+ auto test_curve = [] (Curve const &c) {
+ constexpr int N = 50;
+ for (int i = 0; i < N; i++) {
+ auto angle = 2 * M_PI * i / N;
+ auto transform = Affine(Rotate(angle));
+
+ auto copy = std::unique_ptr<Curve>(c.duplicate());
+ *copy *= transform;
+ auto box1 = copy->boundsExact();
+
+ auto pt = c.initialPoint() * transform;
+ auto box2 = Rect(pt, pt);
+ c.expandToTransformed(box2, transform);
+
+ for (auto i : { X, Y }) {
+ EXPECT_DOUBLE_EQ(box1[i].min(), box2[i].min());
+ EXPECT_DOUBLE_EQ(box1[i].max(), box2[i].max());
+ }
+ }
+ };
+
+ test_curve(LineSegment(Point(-1, 0), Point(1, 2)));
+ test_curve(QuadraticBezier(Point(-1, 0), Point(1, 1), Point(3, 0)));
+ test_curve(CubicBezier(Point(-1, 0), Point(1, 1), Point(2, -2), Point(3, 0)));
+}
+
/*
Local Variables:
mode:c++
diff --git a/tests/choose-test.cpp b/tests/choose-test.cpp
index 766fa0c..57719fe 100644
--- a/tests/choose-test.cpp
+++ b/tests/choose-test.cpp
@@ -47,7 +47,7 @@ TEST(ChooseTest, PascalsTriangle) {
double b = choose<double>(n-1, k);
double c = choose<double>(n-1, k-1);
- EXPECT_EQ(a, b + c);
+ EXPECT_DOUBLE_EQ(a, b + c);
}
}
diff --git a/tests/dependent-project/.gitignore b/tests/dependent-project/.gitignore
new file mode 100644
index 0000000..0bb214c
--- /dev/null
+++ b/tests/dependent-project/.gitignore
@@ -0,0 +1,2 @@
+build-as-subproject
+build-with-find-package
\ No newline at end of file
diff --git a/tests/dependent-project/build-and-run.sh b/tests/dependent-project/build-and-run.sh
deleted file mode 100644
index b7bb2db..0000000
--- a/tests/dependent-project/build-and-run.sh
+++ /dev/null
@@ -1,25 +0,0 @@
-#!/bin/bash
-
-set -x
-mkdir -p ../../opt2
-mkdir -p build-as-subproject
-(
- cd build-as-subproject;
- cmake .. -D2GEOM_AS_SUBPROJECT=ON \
- -DCMAKE_INSTALL_PREFIX=../../../opt2
- make main -j 2
- ./main
- make install
-)
-
-mkdir -p build-with-find-package
-(
- cd build-with-find-package;
- cmake .. -D2GEOM_AS_SUBPROJECT=OFF \
- -DCMAKE_EXE_LINKER_FLAGS="-fsanitize=address" \
- -DCMAKE_INSTALL_PREFIX=../../../opt2 \
- -D2Geom_DIR="$PWD/../../../opt/lib/cmake/2Geom"
- make main -j 2
- ./main
- make install
-)
diff --git a/tests/ellipse-test.cpp b/tests/ellipse-test.cpp
index 2720398..38eca0e 100644
--- a/tests/ellipse-test.cpp
+++ b/tests/ellipse-test.cpp
@@ -4,7 +4,7 @@
*//*
* Authors:
* Krzysztof Kosiński <tweenk.pl@gmail.com>
- *
+ *
* Copyright 2015 Authors
*
* This library is free software; you can redistribute it and/or
@@ -42,7 +42,7 @@
#include "testing.h"
#ifndef M_SQRT2
-# define M_SQRT2 1.41421356237309504880
+# define M_SQRT2 1.41421356237309504880
#endif
using namespace Geom;
@@ -79,7 +79,7 @@ TEST(EllipseTest, Arcs) {
// exactly half arc
std::unique_ptr<EllipticalArc> arc3(e.arc(Point(5,0), Point(0,10), Point(5,20)));
-
+
EXPECT_EQ(arc3->boundsExact(), Rect::from_xywh(0,0,5,20));
EXPECT_EQ(arc3->largeArc(), false);
EXPECT_EQ(arc3->sweep(), false);
@@ -167,6 +167,37 @@ TEST(EllipseTest, LineIntersection) {
EXPECT_NEAR(xs[1].point()[Y], 8./5, 1e-15);
EXPECT_intersections_valid(e, l, xs, 1e-15);
+
+ // Test with a degenerate ellipse
+ auto degen = Ellipse({0, 0}, {3, 2}, 0);
+ degen *= Scale(1.0, 0.0); // Squash to the X-axis interval [-3, 3].
+
+ g_random_set_seed(0xCAFECAFE);
+ // Intersect with a line
+ for (size_t _ = 0; _ < 10'000; _++) {
+ auto line = Line(Point(g_random_double_range(-3.0, 3.0), g_random_double_range(-3.0, -1.0)),
+ Point(g_random_double_range(-3.0, 3.0), g_random_double_range(1.0, 3.0)));
+ auto xings = degen.intersect(line);
+ EXPECT_EQ(xings.size(), 2u);
+ EXPECT_intersections_valid(degen, line, xings, 1e-14);
+ }
+ // Intersect with another, non-degenerate ellipse
+ for (size_t _ = 0; _ < 10'000; _++) {
+ auto other = Ellipse(Point(g_random_double_range(-1.0, 1.0), g_random_double_range(-1.0, 1.0)),
+ Point(g_random_double_range(1.0, 2.0), g_random_double_range(1.0, 3.0)), 0);
+ auto xings = degen.intersect(other);
+ EXPECT_intersections_valid(degen, other, xings, 1e-14);
+ }
+ // Intersect with another ellipse which is also degenerate
+ for (size_t _ = 0; _ < 10'000; _++) {
+ auto other = Ellipse({0, 0}, {1, 1}, 0); // Unit circle
+ other *= Scale(0.0, g_random_double_range(0.5, 4.0)); // Squash to Y axis
+ other *= Rotate(g_random_double_range(-1.5, 1.5)); // Rotate a little (still passes through the origin)
+ other *= Translate(g_random_double_range(-2.9, 2.9), 0.0);
+ auto xings = degen.intersect(other);
+ EXPECT_EQ(xings.size(), 4u);
+ EXPECT_intersections_valid(degen, other, xings, 1e-14);
+ }
}
TEST(EllipseTest, EllipseIntersection) {
@@ -178,7 +209,7 @@ TEST(EllipseTest, EllipseIntersection) {
e2.set(Point(250, 300), Point(230, 90), 1.321);
xs = e1.intersect(e2);
EXPECT_EQ(xs.size(), 4ul);
- EXPECT_intersections_valid(e1, e2, xs, 1e-10);
+ EXPECT_intersections_valid(e1, e2, xs, 4e-10);
e1.set(Point(0, 0), Point(1, 1), 0);
e2.set(Point(0, 1), Point(1, 1), 0);
@@ -191,6 +222,57 @@ TEST(EllipseTest, EllipseIntersection) {
xs = e1.intersect(e2);
EXPECT_EQ(xs.size(), 2ul);
EXPECT_intersections_valid(e1, e2, xs, 1e-10);
+
+ // === Test detection of external tangency between ellipses ===
+ // Perpendicular major axes
+ e1.set({0, 0}, {5, 3}, 0); // rightmost point (5, 0)
+ e2.set({6, 0}, {1, 2}, 0); // leftmost point (5, 0)
+ xs = e1.intersect(e2);
+ ASSERT_GT(xs.size(), 0);
+ EXPECT_intersections_valid(e1, e2, xs, 1e-10);
+ EXPECT_TRUE(are_near(xs[0].point(), Point(5, 0)));
+
+ // Collinear major axes
+ e1.set({30, 0}, {9, 1}, 0); // leftmost point (21, 0)
+ e2.set({18, 0}, {3, 2}, 0); // rightmost point (21, 0)
+ xs = e1.intersect(e2);
+ ASSERT_GT(xs.size(), 0);
+ EXPECT_intersections_valid(e1, e2, xs, 1e-10);
+ EXPECT_TRUE(are_near(xs[0].point(), Point(21, 0)));
+
+ // Circles not aligned to an axis (Pythagorean triple: 3^2 + 4^2 == 5^2)
+ e1.set({0, 0}, {3, 3}, 0); // radius 3
+ e2.set({3, 4}, {2, 2}, 0); // radius 2
+ // We know 2 + 3 == 5 == distance((0, 0), (3, 4)) so there's an external tangency
+ // between these circles, at a point at distance 3 from the origin, on the line x = 0.75 y.
+ xs = e1.intersect(e2);
+ ASSERT_GT(xs.size(), 0);
+ EXPECT_intersections_valid(e1, e2, xs, 1e-6);
+
+ // === Test the detection of internal tangency between ellipses ===
+ // Perpendicular major axes
+ e1.set({0, 0}, {8, 17}, 0); // rightmost point (8, 0)
+ e2.set({6, 0}, {2, 1}, 0); // rightmost point (8, 0)
+ xs = e1.intersect(e2);
+ ASSERT_GT(xs.size(), 0);
+ EXPECT_intersections_valid(e1, e2, xs, 1e-10);
+ EXPECT_TRUE(are_near(xs[0].point(), Point(8, 0)));
+
+ // Collinear major axes
+ e1.set({30, 0}, {9, 5}, 0); // rightmost point (39, 0)
+ e2.set({36, 0}, {3, 1}, 0); // rightmost point (39, 0)
+ xs = e1.intersect(e2);
+ ASSERT_GT(xs.size(), 0);
+ EXPECT_intersections_valid(e1, e2, xs, 1e-6);
+ EXPECT_TRUE(are_near(xs[0].point(), Point(39, 0)));
+
+ // Circles not aligned to an axis (Pythagorean triple: 3^2 + 4^2 == 5^2)
+ e1.set({4, 3}, {5, 5}, 0); // Passes through (0, 0), center on the line y = 0.75 x
+ e2.set({8, 6}, {10, 10}, 0); // Also passes through (0, 0), center on the same line.
+ xs = e1.intersect(e2);
+ ASSERT_GT(xs.size(), 0);
+ EXPECT_intersections_valid(e1, e2, xs, 1e-6);
+ EXPECT_TRUE(are_near(xs[0].point(), Point(0, 0)));
}
TEST(EllipseTest, BezierIntersection) {
@@ -273,7 +355,9 @@ TEST(EllipseTest, UnitTangentAt) {
EXPECT_near(b.unitTangentAt(3*M_PI/2), Point(M_SQRT2/2, M_SQRT2/2), 1e-12);
}
-TEST(EllipseTest, BoundsExact) {
+TEST(EllipseTest, Bounds)
+{
+ // Create example ellipses
std::vector<Ellipse> es;
es.emplace_back(Point(-15,25), Point(10,15), Angle::from_degrees(45));
es.emplace_back(Point(-10,33), Point(40,20), M_PI);
@@ -285,27 +369,42 @@ TEST(EllipseTest, BoundsExact) {
for (auto & e : es) {
Rect r = e.boundsExact();
+ Rect f = e.boundsFast();
for (unsigned j = 0; j < 10000; ++j) {
Coord t = g_random_double_range(-M_PI, M_PI);
- EXPECT_TRUE(r.contains(e.pointAt(t)));
+ auto const p = e.pointAt(t);
+ EXPECT_TRUE(r.contains(p));
+ EXPECT_TRUE(f.contains(p));
}
}
Ellipse e(Point(0,0), Point(10, 10), M_PI);
Rect bounds = e.boundsExact();
+ Rect coarse = e.boundsFast();
EXPECT_EQ(bounds, Rect(Point(-10,-10), Point(10,10)));
EXPECT_TRUE(bounds.contains(e.pointAt(0)));
EXPECT_TRUE(bounds.contains(e.pointAt(M_PI/2)));
EXPECT_TRUE(bounds.contains(e.pointAt(M_PI)));
EXPECT_TRUE(bounds.contains(e.pointAt(3*M_PI/2)));
EXPECT_TRUE(bounds.contains(e.pointAt(2*M_PI)));
+ EXPECT_TRUE(coarse.contains(e.pointAt(0)));
+ EXPECT_TRUE(coarse.contains(e.pointAt(M_PI/2)));
+ EXPECT_TRUE(coarse.contains(e.pointAt(M_PI)));
+ EXPECT_TRUE(coarse.contains(e.pointAt(3*M_PI/2)));
+ EXPECT_TRUE(coarse.contains(e.pointAt(2*M_PI)));
e = Ellipse(Point(0,0), Point(10, 10), M_PI/2);
bounds = e.boundsExact();
+ coarse = e.boundsFast();
EXPECT_EQ(bounds, Rect(Point(-10,-10), Point(10,10)));
EXPECT_TRUE(bounds.contains(e.pointAt(0)));
EXPECT_TRUE(bounds.contains(e.pointAt(M_PI/2)));
EXPECT_TRUE(bounds.contains(e.pointAt(M_PI)));
EXPECT_TRUE(bounds.contains(e.pointAt(3*M_PI/2)));
EXPECT_TRUE(bounds.contains(e.pointAt(2*M_PI)));
+ EXPECT_TRUE(coarse.contains(e.pointAt(0)));
+ EXPECT_TRUE(coarse.contains(e.pointAt(M_PI/2)));
+ EXPECT_TRUE(coarse.contains(e.pointAt(M_PI)));
+ EXPECT_TRUE(coarse.contains(e.pointAt(3*M_PI/2)));
+ EXPECT_TRUE(coarse.contains(e.pointAt(2*M_PI)));
}
diff --git a/tests/elliptical-arc-test.cpp b/tests/elliptical-arc-test.cpp
index 8d1ac47..ff0923d 100644
--- a/tests/elliptical-arc-test.cpp
+++ b/tests/elliptical-arc-test.cpp
@@ -4,7 +4,7 @@
*//*
* Authors:
* Krzysztof Kosiński <tweenk.pl@gmail.com>
- *
+ *
* Copyright 2015 Authors
*
* This library is free software; you can redistribute it and/or
@@ -82,6 +82,79 @@ TEST(EllipticalArcTest, LineSegmentIntersection) {
r1 = a3.intersect(ls);
EXPECT_EQ(r1.size(), 2u);
EXPECT_intersections_valid(a3, ls, r1, 1e-10);
+
+ g_random_set_seed(0xB747A380);
+ // Test with randomized arcs and segments.
+ for (size_t _ = 0; _ < 10'000; _++) {
+ auto arc = EllipticalArc({g_random_double_range(1.0, 5.0), 0.0},
+ {g_random_double_range(6.0, 8.0), g_random_double_range(2.0, 7.0)},
+ g_random_double_range(-0.5, 0.5), true, g_random_boolean(),
+ {g_random_double_range(-5.0, -1.0), 0.0});
+ Coord x = g_random_double_range(15, 30);
+ Coord y = g_random_double_range(10, 20);
+ auto seg = LineSegment(Point(-x, y), Point(x, -y));
+ auto xings = arc.intersect(seg);
+ EXPECT_EQ(xings.size(), 1u);
+ EXPECT_intersections_valid(arc, seg, xings, 1e-12);
+ }
+
+ // Test with degenerate arcs
+ EllipticalArc x_squash_pos{{3.0, 0.0}, {3.0, 2.0}, 0, true, true, {-3.0, 0.0}};
+ EllipticalArc x_squash_neg{{3.0, 0.0}, {3.0, 2.0}, 0, true, false, {-3.0, 0.0}};
+ auto const squash_to_x = Scale(1.0, 0.0);
+ x_squash_pos *= squash_to_x; // squash to X axis interval [-3, 3].
+ x_squash_neg *= squash_to_x;
+
+ for (size_t _ = 0; _ < 10'000; _++) {
+ auto seg = LineSegment(Point(g_random_double_range(-3.0, 3.0), g_random_double_range(-3.0, -1.0)),
+ Point(g_random_double_range(-3.0, 3.0), g_random_double_range(1.0, 3.0)));
+ auto xings = x_squash_pos.intersect(seg);
+ EXPECT_EQ(xings.size(), 1u);
+ EXPECT_intersections_valid(x_squash_pos, seg, xings, 1e-12);
+
+ std::unique_ptr<Curve> rev{x_squash_pos.reverse()};
+ xings = rev->intersect(seg);
+ EXPECT_EQ(xings.size(), 1u);
+ EXPECT_intersections_valid(*rev, seg, xings, 1e-12);
+
+ xings = x_squash_neg.intersect(seg);
+ EXPECT_EQ(xings.size(), 1u);
+ EXPECT_intersections_valid(x_squash_neg, seg, xings, 1e-12);
+
+ rev.reset(x_squash_neg.reverse());
+ xings = rev->intersect(seg);
+ EXPECT_EQ(xings.size(), 1u);
+ EXPECT_intersections_valid(*rev, seg, xings, 1e-12);
+ }
+
+ // Now test with an arc squashed to the Y-axis.
+ EllipticalArc y_squash_pos{{0.0, -2.0}, {3.0, 2.0}, 0, true, true, {0.0, 2.0}};
+ EllipticalArc y_squash_neg{{0.0, -2.0}, {3.0, 2.0}, 0, true, false, {0.0, 2.0}};
+ auto const squash_to_y = Scale(0.0, 1.0);
+ y_squash_pos *= squash_to_y; // Y-axis interval [-2, 2].
+ y_squash_neg *= squash_to_y;
+
+ for (size_t _ = 0; _ < 10'000; _++) {
+ auto seg = LineSegment(Point(g_random_double_range(-3.0, -1.0), g_random_double_range(-2.0, 2.0)),
+ Point(g_random_double_range(1.0, 3.0), g_random_double_range(-2.0, 2.0)));
+ auto xings = y_squash_pos.intersect(seg, 1e-10);
+ EXPECT_EQ(xings.size(), 1u);
+ EXPECT_intersections_valid(y_squash_pos, seg, xings, 1e-12);
+
+ std::unique_ptr<Curve> rev{y_squash_pos.reverse()};
+ xings = rev->intersect(seg, 1e-12);
+ EXPECT_EQ(xings.size(), 1u);
+ EXPECT_intersections_valid(*rev, seg, xings, 1e-12);
+
+ xings = y_squash_neg.intersect(seg, 1e-12);
+ EXPECT_EQ(xings.size(), 1u);
+ EXPECT_intersections_valid(y_squash_neg, seg, xings, 1e-12);
+
+ rev.reset(y_squash_neg.reverse());
+ xings = rev->intersect(seg, 1e-12);
+ EXPECT_EQ(xings.size(), 1u);
+ EXPECT_intersections_valid(*rev, seg, xings, 1e-12);
+ }
}
TEST(EllipticalArcTest, ArcIntersection) {
@@ -98,6 +171,39 @@ TEST(EllipticalArcTest, ArcIntersection) {
r2 = a3.intersect(a4);
EXPECT_EQ(r2.size(), 3u);
EXPECT_intersections_valid(a3, a4, r2, 1e-10);
+
+ // Make sure intersections are found between two identical arcs on the unit circle.
+ EllipticalArc const upper(Point(1, 0), Point(1, 1), 0, true, true, Point(-1, 0));
+ auto self_intersect = upper.intersect(upper);
+ EXPECT_EQ(self_intersect.size(), 2u);
+
+ // Make sure intersections are found between overlapping arcs.
+ EllipticalArc const right(Point(0, -1), Point(1, 1), 0, true, true, Point(0, 1));
+ auto quartering_overlap_xings = right.intersect(upper);
+ EXPECT_EQ(quartering_overlap_xings.size(), 2u);
+
+ // Make sure intersecections are found between an arc and its sub-arc.
+ EllipticalArc const middle(upper.pointAtAngle(0.25 * M_PI), Point(1, 1), 0, true, true, upper.pointAtAngle(-0.25 * M_PI));
+ EXPECT_EQ(middle.intersect(upper).size(), 2u);
+
+ // Make sure intersections are NOT found between non-overlapping sub-arcs of the same circle.
+ EllipticalArc const arc1{Point(1, 0), Point(1, 1), 0, true, true, Point(0, 1)};
+ EllipticalArc const arc2{Point(-1, 0), Point(1, 1), 0, true, true, Point(0, -1)};
+ EXPECT_EQ(arc1.intersect(arc2).size(), 0u);
+
+ // Overlapping sub-arcs but on an Ellipse with different rays.
+ EllipticalArc const eccentric{Point(2, 0), Point(2, 1), 0, true, true, Point(-2, 0)};
+ EllipticalArc const subarc{eccentric.pointAtAngle(0.8), Point(2, 1), 0, true, true, eccentric.pointAtAngle(2)};
+ EXPECT_EQ(eccentric.intersect(subarc).size(), 2u);
+
+ // Check intersection times for two touching arcs.
+ EllipticalArc const lower{Point(-1, 0), Point(1, 1), 0, false, true, Point(0, -1)};
+ auto expected_neg_x = upper.intersect(lower);
+ ASSERT_EQ(expected_neg_x.size(), 1);
+ auto const &left_pt = expected_neg_x[0];
+ EXPECT_EQ(left_pt.point(), Point(-1, 0));
+ EXPECT_DOUBLE_EQ(left_pt.first, 1.0); // Expect (-1, 0) reached at the end of upper
+ EXPECT_DOUBLE_EQ(left_pt.second, 0.0); // Expect (-1, 0) passed at the start of lower
}
TEST(EllipticalArcTest, BezierIntersection) {
diff --git a/tests/intersection-graph-test.cpp b/tests/intersection-graph-test.cpp
index 5335ad1..19fb25c 100644
--- a/tests/intersection-graph-test.cpp
+++ b/tests/intersection-graph-test.cpp
@@ -212,6 +212,30 @@ TEST_F(IntersectionGraphTest, RhombusInSquare) {
checkRandomPoints(square, rhombus, s2, B_MINUS_A);
}
+TEST_F(IntersectionGraphTest, EmptyOperand) {
+ PathVector square = string_to_path("M 0,0 L 20, 0 L 20, 20 L 0, 20 Z");
+ PathVector empty;
+
+ auto graph = PathIntersectionGraph(square, empty);
+ // Taking union with the empty set should be a no-op: A ∪ ∅ = A
+ PathVector u = graph.getUnion();
+ EXPECT_EQ(u.size(), 1u);
+ EXPECT_EQ(u.curveCount(), 4u);
+
+ // Intersection with empty should produce empty: A ∩ ∅ = ∅
+ PathVector i = graph.getIntersection();
+ EXPECT_EQ(i.size(), 0u);
+
+ // Subtracting empty set should be a no-op: A ∖ ∅ = A
+ PathVector rd = graph.getAminusB();
+ EXPECT_EQ(rd.size(), 1u);
+ EXPECT_EQ(rd.curveCount(), 4u);
+
+ // Subtracting FROM the empty set should produce the empty set: ∅ ∖ A = ∅
+ PathVector ld = graph.getBminusA();
+ EXPECT_EQ(ld.size(), 0u);
+}
+
// this test is disabled, since we cannot handle overlapping segments for now.
#if 0
TEST_F(IntersectionGraphTest, EqualUnionAndIntersection) {
diff --git a/tests/path-test.cpp b/tests/path-test.cpp
index a4a4243..cd5db23 100644
--- a/tests/path-test.cpp
+++ b/tests/path-test.cpp
@@ -1,4 +1,6 @@
-#include "testing.h"
+#include <cmath>
+#include <vector>
+#include <iterator>
#include <iostream>
#include <2geom/bezier.h>
@@ -7,8 +9,8 @@
#include <2geom/path-intersection.h>
#include <2geom/svg-path-parser.h>
#include <2geom/svg-path-writer.h>
-#include <vector>
-#include <iterator>
+
+#include "testing.h"
using namespace std;
using namespace Geom;
@@ -408,6 +410,45 @@ TEST_F(PathTest, AppendPath) {
EXPECT_NO_THROW(p_closed.checkContinuity());
}
+TEST_F(PathTest, AppendPortion) {
+ // A closed path with two curves:
+ Path bigon = string_to_path("M 0,0 Q 1,1 2,0 Q 1,-1 0,0 Z");
+ Path target{Point(0, 0)};
+
+ PathTime end_time{1, 1.0}; // End of the closed path
+ PathTime mid_time{1, 0.0}; // Middle of the closed path (juncture between the two curves)
+ bigon.appendPortionTo(target, end_time, mid_time, true /* do cross start */);
+
+ // We expect that the target path now contains the entire first curve "M 0,0 Q 1,1 2,0",
+ // since we started at the end of a closed path and requested to cross its start.
+ EXPECT_EQ(target.size(), 1);
+ EXPECT_EQ(target, string_to_path("M 0,0 Q 1,1 2,0"));
+
+ // Similar test but with reversal (swapped times)
+ Path target_reverse{Point(2, 0)};
+ bigon.appendPortionTo(target_reverse, mid_time, end_time, true /* do cross start please */);
+ // What do we expect? To cross start going from the midpoint to the endpoint requires
+ // not taking the obvious route (bigon[1]) but rather taking bigon[0] in reverse.
+ EXPECT_EQ(target_reverse.size(), 1);
+ EXPECT_EQ(target_reverse, string_to_path("M 2,0 Q 1,1 0,0"));
+
+ // Similar test but using start time
+ PathTime start_time{0, 0.0};
+ Path mid_target{Point(2, 0)};
+ bigon.appendPortionTo(mid_target, mid_time, start_time, true /* cross start to 0:0 */);
+ // We expect to go forward from mid_time and cross over the start to start_time.
+ EXPECT_EQ(mid_target.size(), 1);
+ EXPECT_EQ(mid_target, string_to_path("M 2,0 Q 1,-1 0,0"));
+
+ // Use start time with reversal
+ Path mid_reverse{Point(0, 0)};
+ bigon.appendPortionTo(mid_reverse, start_time, mid_time, true /* Cross start, going backwards. */);
+ // We expect that we don't go forwards from start_time to mid_time, but rather cross over the starting
+ // point and backtrack over bigon[1] to the midpoint.
+ EXPECT_EQ(mid_reverse.size(), 1);
+ EXPECT_EQ(mid_reverse, string_to_path("M 0,0 Q 1,-1 2,0"));
+}
+
TEST_F(PathTest, ReplaceMiddle) {
p_open.replace(p_open.begin() + 1, p_open.begin() + 2, p_add);
EXPECT_EQ(p_open.size(), 5u);
@@ -495,6 +536,378 @@ TEST_F(PathTest, Roots) {
EXPECT_EQ(roots.size(), 2u);
}
+TEST_F(PathTest, PartingPoint)
+{
+ // === Test complete overlaps between identical curves ===
+ // Line segment
+ auto line = string_to_path("M 0,0 L 3.33, 7.77");
+ auto pt = parting_point(line, line);
+ EXPECT_TRUE(are_near(pt.point(), line.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.t, 1.0));
+
+ // Cubic Bézier
+ auto bezier = string_to_path("M 0,0 C 1,1 14,1 15,0");
+ pt = parting_point(bezier, bezier);
+ EXPECT_TRUE(are_near(pt.point(), bezier.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.t, 1.0));
+
+ // Eliptical arc
+ auto const arc = string_to_path("M 0,0 A 100,20 0,0,0 200,0");
+ pt = parting_point(arc, arc);
+ EXPECT_TRUE(are_near(pt.point(), arc.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.t, 1.0));
+
+ // === Test complete overlap between degree-elevated and degree-shrunk Béziers ===
+ auto artificially_cubic = string_to_path("M 0,0 C 10,10 20,10 30,0");
+ auto really_quadratic = string_to_path("M 0,0 Q 15,15 30,0");
+ pt = parting_point(artificially_cubic, really_quadratic);
+ EXPECT_TRUE(are_near(pt.point(), artificially_cubic.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.asFlatTime(), 1.0));
+ EXPECT_TRUE(are_near(pt.second.asFlatTime(), 1.0));
+
+ // === Test complete overlaps between a curve and its subdivision ===
+ // Straight line
+ line = string_to_path("M 0,0 L 15,15");
+ auto subdivided_line = string_to_path("M 0,0 L 3,3 L 4,4 L 9,9 L 15,15");
+ pt = parting_point(line, subdivided_line);
+ EXPECT_TRUE(are_near(pt.point(), line.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.t, 1.0));
+
+ // Cubic Bézier
+ bezier = string_to_path("M 0,0 C 0,40 50,40 50,0");
+ auto de_casteljau = string_to_path("M 0,0 C 0,10 3.125,17.5 7.8125,22.5 12.5,27.5 18.75,30 25,30"
+ " 31.25,30 37.5,27.5 42.1875,22.5 46.875,17.5 50,10 50,0");
+ pt = parting_point(bezier, de_casteljau);
+ EXPECT_TRUE(are_near(pt.point(), bezier.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.t, 1.0));
+
+ // Eliptical arc
+ auto subdivided_arc = string_to_path("M 0,0 A 100,20, 0,0,0 100,20 A 100,20 0,0,0 200,0");
+ pt = parting_point(arc, subdivided_arc);
+ EXPECT_TRUE(are_near(pt.point(), arc.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.t, 1.0));
+
+ // === Test complete overlap between different subdivisions ===
+ auto line1 = string_to_path("M 0,0 L 3,3 L 5,5 L 10,10");
+ auto line2 = string_to_path("M 0,0 L 2,2 L 4.2,4.2 L 4.5,4.5 L 6,6 L 10,10");
+ pt = parting_point(line1, line2);
+ EXPECT_TRUE(are_near(pt.point(), line1.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.asFlatTime(), line1.timeRange().max()));
+ EXPECT_TRUE(are_near(pt.second.asFlatTime(), line2.timeRange().max()));
+
+ // === Test complete overlaps in the presence of degenerate segments ===
+ // Straight line
+ line = string_to_path("M 0,0 L 15,15");
+ subdivided_line = string_to_path("M 0,0 L 3,3 H 3 V 3 L 3,3 L 4,4 H 4 V 4 L 4,4 L 9,9 H 9 L 15,15");
+ pt = parting_point(line, subdivided_line);
+ EXPECT_TRUE(are_near(pt.point(), line.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.asFlatTime(), 1.0));
+
+ // Eliptical arc
+ auto arc_degen = string_to_path("M 0,0 A 100,20, 0,0,0 100,20 H 100 V 20 L 100,20 A 100,20 0,0,0 200,0");
+ pt = parting_point(arc, arc_degen);
+ EXPECT_TRUE(are_near(pt.point(), arc.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.asFlatTime(), 1.0));
+
+ // === Paths that overlap but one is shorter than the other ===
+ // Straight lines
+ auto long_line = string_to_path("M 0,0 L 20,10");
+ auto short_line = string_to_path("M 0,0 L 4,2");
+ pt = parting_point(long_line, short_line);
+ EXPECT_TRUE(are_near(pt.point(), short_line.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.t, 0.2));
+ EXPECT_TRUE(are_near(pt.second.t, 1.0));
+
+ // Cubic Bézier
+ auto const s_shape = string_to_path("M 0,0 C 10, 0 0,10 10,10");
+ auto half_s = string_to_path("M 0,0 C 5,0 5,2.5 5,5");
+ pt = parting_point(s_shape, half_s);
+ EXPECT_TRUE(are_near(pt.first.t, 0.5));
+ EXPECT_TRUE(are_near(pt.second.t, 1.0));
+
+ // Elliptical arc
+ auto quarter_ellipse = string_to_path("M 0,0 A 100,20, 0,0,0 100,20");
+ pt = parting_point(arc, quarter_ellipse);
+ EXPECT_TRUE(are_near(pt.point(), quarter_ellipse.finalPoint()));
+ EXPECT_TRUE(are_near(pt.first.t, 0.5));
+ EXPECT_TRUE(are_near(pt.second.t, 1.0));
+
+ // === Paths that overlap initially but then they split ===
+ // Straight lines
+ auto boring_line = string_to_path("M 0,0 L 50,10");
+ auto line_then_arc = string_to_path("M 0,0 L 5,1 A 1,1 0,0,0 7,1");
+ pt = parting_point(boring_line, line_then_arc);
+ EXPECT_TRUE(are_near(pt.point(), Point(5, 1)));
+ EXPECT_TRUE(are_near(pt.first.t, 0.1));
+ EXPECT_TRUE(are_near(pt.second.asFlatTime(), 1.0));
+
+ // Cubic Bézier
+ auto half_s_then_line = string_to_path("M 0,0 C 5,0 5,2.5 5,5 L 10,10");
+ pt = parting_point(s_shape, half_s_then_line);
+ EXPECT_TRUE(are_near(pt.point(), Point(5, 5)));
+ EXPECT_TRUE(are_near(pt.first.t, 0.5));
+ EXPECT_TRUE(are_near(pt.second.asFlatTime(), 1.0));
+
+ // Elliptical arc
+ auto quarter_ellipse_then_quadratic = string_to_path("M 0,0 A 100,20, 0,0,0 100,20 Q 120,40 140,60");
+ pt = parting_point(arc, quarter_ellipse_then_quadratic);
+ EXPECT_TRUE(are_near(pt.point(), Point(100, 20)));
+ EXPECT_TRUE(are_near(pt.first.t, 0.5));
+ EXPECT_TRUE(are_near(pt.second.asFlatTime(), 1.0));
+
+ // === Paths that split at a common node ===
+ // Polylines
+ auto branch_90 = string_to_path("M 0,0 H 3 H 6 V 7");
+ auto branch_45 = string_to_path("M 0,0 H 2 H 6 L 7,7");
+ pt = parting_point(branch_90, branch_45);
+ EXPECT_TRUE(are_near(pt.point(), Point(6, 0)));
+ EXPECT_TRUE(are_near(pt.first.asFlatTime(), 2.0));
+ EXPECT_TRUE(are_near(pt.second.asFlatTime(), 2.0));
+
+ // Arcs
+ auto quarter_circle_then_horiz = string_to_path("M 0,0 A 1,1 0,0,0 1,1 H 10");
+ auto quarter_circle_then_slant = string_to_path("M 0,0 A 1,1 0,0,0 1,1 L 10, 1.1");
+ pt = parting_point(quarter_circle_then_horiz, quarter_circle_then_slant);
+ EXPECT_TRUE(are_near(pt.point(), Point(1, 1)));
+ EXPECT_TRUE(are_near(pt.first.asFlatTime(), 1.0));
+ EXPECT_TRUE(are_near(pt.second.asFlatTime(), 1.0));
+
+ // Last common nodes followed by degenerates
+ auto degen_horiz = string_to_path("M 0,0 A 1,1 0,0,0 1,1 V 1 H 1 L 1,1 H 10");
+ auto degen_slant = string_to_path("M 0,0 A 1,1 0,0,0 1,1 V 1 H 1 L 1,1 L 10, 1.1");
+ pt = parting_point(quarter_circle_then_horiz, quarter_circle_then_slant);
+ EXPECT_TRUE(are_near(pt.point(), Point(1, 1)));
+
+ // === Paths that split at the starting point ===
+ auto vertical = string_to_path("M 0,0 V 1");
+ auto quarter = string_to_path("M 0,0 A 1,1 0,0,0, 1,1");
+ pt = parting_point(vertical, quarter);
+ EXPECT_TRUE(are_near(pt.point(), Point(0, 0)));
+ EXPECT_TRUE(are_near(pt.first.asFlatTime(), 0.0));
+ EXPECT_TRUE(are_near(pt.second.asFlatTime(), 0.0));
+
+ // === Symmetric split (both legs of the same length) ===
+ auto left_leg = string_to_path("M 1,0 L 0,10");
+ auto right_leg = string_to_path("M 1,0 L 2,10");
+ pt = parting_point(left_leg, right_leg);
+ EXPECT_TRUE(are_near(pt.point(), Point(1, 0)));
+ EXPECT_TRUE(are_near(pt.first.asFlatTime(), 0.0));
+ EXPECT_TRUE(are_near(pt.second.asFlatTime(), 0.0));
+
+ // === Different starting points ===
+ auto start_at_0_0 = string_to_path("M 0,0 C 1,0 0,1 1,1");
+ auto start_at_10_10 = string_to_path("M 10,10 L 50,50");
+ pt = parting_point(start_at_0_0, start_at_10_10);
+ EXPECT_TRUE(are_near(pt.point(), Point (5,5)));
+ EXPECT_DOUBLE_EQ(pt.first.t, -1.0);
+ EXPECT_DOUBLE_EQ(pt.second.t, -1.0);
+ EXPECT_EQ(pt.first.curve_index, 0);
+ EXPECT_EQ(pt.second.curve_index, 0);
+}
+
+TEST_F(PathTest, InitialFinalTangents) {
+ // Test tangents for an open path
+ auto L_shape = string_to_path("M 1,1 H 0 V 0");
+ EXPECT_EQ(L_shape.initialUnitTangent(), Point(-1.0, 0.0));
+ EXPECT_EQ(L_shape.finalUnitTangent(), Point(0.0, -1.0));
+
+ // Closed path with non-degenerate closing segment
+ auto triangle = string_to_path("M 0,0 H 2 L 0,3 Z");
+ EXPECT_EQ(triangle.initialUnitTangent(), Point(1.0, 0.0));
+ EXPECT_EQ(triangle.finalUnitTangent(), Point(0.0, -1.0));
+
+ // Closed path with a degenerate closing segment
+ auto full360 = string_to_path("M 0,0 A 1,1, 0,1,1, 0,2 A 1,1 0,1,1 0,0 Z");
+ EXPECT_EQ(full360.initialUnitTangent(), Point(1.0, 0.0));
+ EXPECT_EQ(full360.finalUnitTangent(), Point(1.0, 0.0));
+
+ // Test multiple degenerate segments at the start
+ auto start_degen = string_to_path("M 0,0 L 0,0 H 0 V 0 Q 1,0 1,1");
+ EXPECT_EQ(start_degen.initialUnitTangent(), Point(1.0, 0.0));
+
+ // Test multiple degenerate segments at the end
+ auto end_degen = string_to_path("M 0,0 L 1,1 H 1 V 1 L 1,1");
+ double comp = 1.0 / sqrt(2.0);
+ EXPECT_EQ(end_degen.finalUnitTangent(), Point(comp, comp));
+
+ // Test a long and complicated path with both tangents along the positive x-axis.
+ auto complicated = string_to_path("M 0,0 H 0 L 1,0 C 2,1 3,2 1,0 L 1,0 H 1 Q 2,3 0,5 H 2");
+ EXPECT_EQ(complicated.initialUnitTangent(), Point(1.0, 0.0));
+ EXPECT_EQ(complicated.finalUnitTangent(), Point(1.0, 0.0));
+}
+
+TEST_F(PathTest, WithoutDegenerates) {
+ // Ensure nothing changes when there are no degenerate segments to remove.
+ auto plain_open = string_to_path("M 0,0 Q 5,5 10,10");
+ EXPECT_EQ(plain_open, plain_open.withoutDegenerateCurves());
+
+ auto closed_nondegen_closing = string_to_path("M 0,0 L 5,5 H 0 Z");
+ EXPECT_EQ(closed_nondegen_closing,closed_nondegen_closing.withoutDegenerateCurves());
+
+ // Ensure that a degenerate closing segment is left alone.
+ auto closed_degen_closing = string_to_path("M 0,0 L 2,4 H 0 L 0,0 Z");
+ EXPECT_EQ(closed_degen_closing, closed_degen_closing.withoutDegenerateCurves());
+
+ // Ensure that a trivial path is left alone (both open and closed).
+ auto trivial_open = string_to_path("M 0,0");
+ EXPECT_EQ(trivial_open, trivial_open.withoutDegenerateCurves());
+
+ auto trivial_closed = string_to_path("M 0,0 Z");
+ EXPECT_EQ(trivial_closed, trivial_closed.withoutDegenerateCurves());
+
+ // Ensure that initial degenerate segments are removed
+ auto degen_start = string_to_path("M 0,0 L 0,0 H 0 V 0 Q 5,5 10,10");
+ auto degen_start_cleaned = degen_start.withoutDegenerateCurves();
+ EXPECT_EQ(degen_start_cleaned, string_to_path("M 0,0 Q 5,5 10,10"));
+ EXPECT_NE(degen_start.size(), degen_start_cleaned.size());
+
+ // Ensure that degenerate segments are removed from the middle
+ auto degen_middle = string_to_path("M 0,0 L 1,1 H 1 V 1 L 1,1 Q 6,6 10,10");
+ auto degen_middle_cleaned = degen_middle.withoutDegenerateCurves();
+ EXPECT_EQ(degen_middle_cleaned, string_to_path("M 0,0 L 1,1 Q 6,6 10,10"));
+ EXPECT_NE(degen_middle.size(), degen_middle_cleaned.size());
+
+ // Ensure that degenerate segment are removed from the end of an open path
+ auto end_open = string_to_path("M 0,0 L 1,1 H 1 V 1 L 1,1");
+ auto end_open_cleaned = end_open.withoutDegenerateCurves();
+ EXPECT_EQ(end_open_cleaned, string_to_path("M 0,0 L 1,1"));
+ EXPECT_NE(end_open.size(), end_open_cleaned.size());
+
+ // Ensure removal of degenerates just before the closing segment
+ auto end_nondegen = string_to_path("M 0,0 L 1,1 L 0,1 H 0 V 1 Z");
+ auto end_nondegen_cleaned = end_nondegen.withoutDegenerateCurves();
+ EXPECT_EQ(end_nondegen_cleaned, string_to_path("M 0,0 L 1,1 L 0,1 Z"));
+ EXPECT_NE(end_nondegen.size(), end_nondegen_cleaned.size());
+}
+
+/** Test Path::extrema() */
+TEST_F(PathTest, GetExtrema) {
+
+ // Circle of radius 4.5 centered at (-4.5, 0).
+ auto extrema_x = circle.extrema(X);
+ EXPECT_EQ(extrema_x.min_point, Point(-9, 0));
+ EXPECT_EQ(extrema_x.max_point, Point( 0, 0));
+ EXPECT_DOUBLE_EQ(extrema_x.min_time.asFlatTime(), 1.0);
+ EXPECT_DOUBLE_EQ(extrema_x.max_time.asFlatTime(), 0.0);
+ EXPECT_EQ(extrema_x.glance_direction_at_min, -1.0);
+ EXPECT_EQ(extrema_x.glance_direction_at_max, 1.0);
+
+ auto extrema_y = circle.extrema(Y);
+ EXPECT_EQ(extrema_y.min_point, Point(-4.5, -4.5));
+ EXPECT_EQ(extrema_y.max_point, Point(-4.5, 4.5));
+ EXPECT_DOUBLE_EQ(extrema_y.min_time.asFlatTime(), 1.5);
+ EXPECT_DOUBLE_EQ(extrema_y.max_time.asFlatTime(), 0.5);
+ EXPECT_FLOAT_EQ(extrema_y.glance_direction_at_min, 1.0);
+ EXPECT_FLOAT_EQ(extrema_y.glance_direction_at_max, -1.0);
+
+ // Positively oriented unit square
+ extrema_x = square.extrema(X);
+ EXPECT_DOUBLE_EQ(extrema_x.min_point[X], 0.0);
+ EXPECT_DOUBLE_EQ(extrema_x.max_point[X], 1.0);
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_min, -1.0);
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_max, 1.0);
+
+ extrema_y = square.extrema(Y);
+ EXPECT_DOUBLE_EQ(extrema_y.min_point[Y], 0.0);
+ EXPECT_DOUBLE_EQ(extrema_y.max_point[Y], 1.0);
+ EXPECT_FLOAT_EQ(extrema_y.glance_direction_at_min, 1.0);
+ EXPECT_FLOAT_EQ(extrema_y.glance_direction_at_max, -1.0);
+
+ // Path glancing its min X line while going towards negative Y
+ auto down_glance = string_to_path("M 1,18 L 0,0 1,-20");
+ extrema_x = down_glance.extrema(X);
+ EXPECT_EQ(extrema_x.min_point, Point(0, 0));
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_min, -1.0);
+ EXPECT_DOUBLE_EQ(extrema_x.min_time.asFlatTime(), 1.0);
+
+ // Similar but not at a node
+ auto down_glance_smooth = string_to_path("M 1,20 C 0,20 0,-20 1,-20");
+ extrema_x = down_glance_smooth.extrema(X);
+ EXPECT_TRUE(are_near(extrema_x.min_point[Y], 0.0));
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_min, -1.0);
+ EXPECT_DOUBLE_EQ(extrema_x.min_time.asFlatTime(), 0.5);
+
+ // Path coming down to the min X and then retreating horizontally
+ auto retreat = string_to_path("M 1,20 L 0,0 H 5 L 4,-20");
+ extrema_x = retreat.extrema(X);
+ EXPECT_EQ(extrema_x.min_point, Point(0, 0));
+ EXPECT_EQ(extrema_x.max_point, Point(5, 0));
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_min, -1.0);
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_max, -1.0);
+ EXPECT_DOUBLE_EQ(extrema_x.min_time.asFlatTime(), 1.0);
+ EXPECT_DOUBLE_EQ(extrema_x.max_time.asFlatTime(), 2.0);
+
+ // Perfectly horizontal path
+ auto horizontal = string_to_path("M 0,0 H 12");
+ extrema_x = horizontal.extrema(X);
+ extrema_y = horizontal.extrema(Y);
+ EXPECT_EQ(extrema_x.min_point, Point(0, 0));
+ EXPECT_EQ(extrema_x.max_point, Point(12, 0));
+ EXPECT_DOUBLE_EQ(extrema_y.min_point[Y], 0.0);
+ EXPECT_DOUBLE_EQ(extrema_y.max_point[Y], 0.0);
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_min, 0.0);
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_max, 0.0);
+ EXPECT_FLOAT_EQ(extrema_y.glance_direction_at_min, 1.0);
+ EXPECT_FLOAT_EQ(extrema_y.glance_direction_at_max, 1.0);
+ EXPECT_DOUBLE_EQ(extrema_x.min_time.asFlatTime(), 0.0);
+ EXPECT_DOUBLE_EQ(extrema_x.max_time.asFlatTime(), 1.0);
+
+ // Perfectly vertical path
+ auto vertical = string_to_path("M 0,0 V 42");
+ extrema_y = vertical.extrema(Y);
+ extrema_x = vertical.extrema(X);
+ EXPECT_DOUBLE_EQ(extrema_x.min_point[Y], 0.0);
+ EXPECT_DOUBLE_EQ(extrema_x.max_point[Y], 0.0);
+ EXPECT_EQ(extrema_y.min_point, Point(0, 0));
+ EXPECT_EQ(extrema_y.max_point, Point(0, 42));
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_min, 1.0);
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_max, 1.0);
+ EXPECT_FLOAT_EQ(extrema_y.glance_direction_at_min, 0.0);
+ EXPECT_FLOAT_EQ(extrema_y.glance_direction_at_max, 0.0);
+ EXPECT_DOUBLE_EQ(extrema_y.min_time.asFlatTime(), 0.0);
+ EXPECT_DOUBLE_EQ(extrema_y.max_time.asFlatTime(), 1.0);
+
+ // Detect downward glance at the closing point (degenerate closing segment)
+ auto closed = string_to_path("M 0,0 L 1,-2 H 3 V 5 H 1 L 0,0 Z");
+ extrema_x = closed.extrema(X);
+ EXPECT_EQ(extrema_x.min_point, Point(0, 0));
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_min, -1.0);
+
+ // Same but with a non-degenerate closing segment
+ auto closed_nondegen = string_to_path("M 0,0 L 1,-2 H 3 V 5 H 1 Z");
+ extrema_x = closed_nondegen.extrema(X);
+ EXPECT_EQ(extrema_x.min_point, Point(0, 0));
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_min, -1.0);
+
+ // Collapsed Bezier not glancing up nor down
+ auto collapsed = string_to_path("M 10, 0 Q -10 0 10, 0");
+ extrema_x = collapsed.extrema(X);
+ EXPECT_EQ(extrema_x.min_point, Point(0, 0));
+ EXPECT_EQ(extrema_x.max_point, Point(10, 0));
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_min, 0.0);
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_max, 0.0);
+
+ // Degenerate segments at min X
+ auto degen = string_to_path("M 0.01,20 L 0, 0 H 0 V 0 L 0,0 V 0 L 0.02 -30");
+ extrema_x = degen.extrema(X);
+ EXPECT_EQ(extrema_x.min_point, Point(0, 0));
+ EXPECT_FLOAT_EQ(extrema_x.glance_direction_at_min, -1.0);
+}
+
+/** Regression test for issue https://gitlab.com/inkscape/lib2geom/-/issues/50 */
+TEST_F(PathTest, PizzaSlice)
+{
+ auto pv = parse_svg_path("M 0 0 L 0.30901699437494745 0.9510565162951535 "
+ "A 1 1 0 0 1 -0.8090169943749473 0.5877852522924732 z");
+ auto §or = pv[0];
+ Path piece;
+ EXPECT_NO_THROW(piece = sector.portion(PathTime(0, 0.0), PathTime(2, 0.0), false));
+ EXPECT_FALSE(piece.closed());
+ EXPECT_TRUE(piece.size() == 2 ||
+ (piece.size() == 3 && piece[2].isDegenerate()));
+ EXPECT_EQ(piece.finalPoint(), Point(-0.8090169943749473, 0.5877852522924732));
+}
+
/*
Local Variables:
mode:c++
diff --git a/tests/planar-graph-test.cpp b/tests/planar-graph-test.cpp
new file mode 100644
index 0000000..f19e2eb
--- /dev/null
+++ b/tests/planar-graph-test.cpp
@@ -0,0 +1,457 @@
+/** @file
+ * @brief Unit tests for PlanarGraph class template
+ */
+/*
+ * Authors:
+ * Rafał Siejakowski <rs@rs-math.net>
+ *
+ * Copyright 2022 the Authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#include <gtest/gtest.h>
+#include <iostream>
+
+#include <2geom/point.h>
+#include <2geom/pathvector.h>
+#include <2geom/svg-path-parser.h>
+#include <2geom/svg-path-writer.h>
+
+#include "planar-graph.h"
+#include "testing.h"
+
+using namespace Geom;
+
+#define PV(d) (parse_svg_path(d))
+#define PTH(d) (std::move(PV(d)[0]))
+#define REV(d) ((PV(d)[0]).reversed())
+
+/** An edge label for the purpose of tests. */
+struct TestLabel
+{
+ unsigned reversal_count = 0, merge_count = 0, detachment_count = 0;
+ void onReverse() { reversal_count++; }
+ void onMergeWith(TestLabel const &) { merge_count++; }
+ void onDetach() { detachment_count++; }
+};
+
+using TestGraph = PlanarGraph<TestLabel>;
+
+static std::vector<TestLabel> extract_labels(TestGraph const &graph)
+{
+ // Find labels of edges remaining in the graph.
+ std::vector<TestLabel> result;
+ for (auto &e : graph.getEdges()) {
+ if (!e.detached) {
+ result.push_back(e.label);
+ }
+ }
+ return result;
+}
+
+class PlanarGraphTest : public ::testing::Test
+{
+};
+
+/** Test edge insertion and vertex clumping to within the tolerance. */
+TEST(PlanarGraphTest, EdgeInsertion)
+{
+ double const precision = 1e-3;
+ auto graph = TestGraph(precision);
+ graph.insertEdge(PTH("M 0, 0 L 1, 0"));
+ graph.insertEdge(PTH("M 0, 1 L 1, 1")); // } Endpoints near
+ graph.insertEdge(PTH("M 1, 0 L 1, 1.0009")); // } but not exact.
+
+ auto vertices = graph.getVertices();
+
+ // Test vertex clumping within the given precision
+ EXPECT_EQ(vertices.size(), 4);
+ EXPECT_EQ(graph.numEdges(), 3);
+
+ // Test lexicographic vertex position sorting by X and then Y
+ EXPECT_EQ(vertices.front().point(), Point(0, 0));
+ auto after = std::next(vertices.begin());
+ EXPECT_EQ(after->point(), Point(0, 1));
+ ++after;
+ EXPECT_EQ(after->point(), Point(1, 0));
+ EXPECT_TRUE(are_near(vertices.back().point(), Point(1, 1), precision));
+
+ EXPECT_FALSE(graph.isRegularized());
+}
+
+/** Test PlanarGraph<T>::insertDetached(). */
+TEST(PlanarGraphTest, InsertDetached)
+{
+ TestGraph graph;
+ auto detached = graph.insertDetached(PTH("M 0,0 A 1,1 0,0,1 2,0 V -2 H 0 Z"));
+
+ auto const &edges = graph.getEdges();
+ EXPECT_EQ(edges.size(), 1);
+ EXPECT_TRUE(edges.at(detached).detached);
+ EXPECT_TRUE(edges.at(detached).inserted_as_detached);
+
+ EXPECT_EQ(graph.numVertices(), 0);
+ EXPECT_EQ(graph.numEdges(false), 0);
+ EXPECT_TRUE(graph.isRegularized());
+}
+
+/** Test signed area calculation. */
+TEST(PlanarGraphTest, ClosedPathArea)
+{
+ // Square with counter-clockwise oriented boundary, when imagining that the y-axis
+ // points up – expect the area to be +1.
+ auto square_positive = PTH("M 0,0 H 1 V 1 H 0 Z");
+ EXPECT_DOUBLE_EQ(TestGraph::closedPathArea(square_positive), 1.0);
+
+ // Expect negative area for a negatively oriented path.
+ auto triangle_negative = PTH("M 0,0 V 1 L 1,1 Z");
+ EXPECT_DOUBLE_EQ(TestGraph::closedPathArea(triangle_negative), -0.5);
+}
+
+/** Test the detection of direction of deviation of initially tangent paths. */
+TEST(PlanarGraphTest, Deviation)
+{
+ auto vertical_up = PTH("M 0,0 V 1");
+ auto arc_right1 = PTH("M 0,0 A 1,1 0,1,0 2,0");
+ auto arc_left1 = PTH("M 0,0 A 1,1 0,1,1 -2,0");
+ auto arc_right2 = PTH("M 0,0 A 2,2 0,1,0, 4,0");
+ auto arc_left2 = PTH("M 0,0 A 2,2 0,1,1 -4,0");
+ // A very "flat" Bézier curve deviating to the right but slower than the large arc
+ auto bezier_right = PTH("M 0,0 C 0,50 1,20 2,10");
+
+ EXPECT_TRUE(TestGraph::deviatesLeft(arc_left1, arc_left2));
+ EXPECT_TRUE(TestGraph::deviatesLeft(arc_left2, vertical_up));
+ EXPECT_TRUE(TestGraph::deviatesLeft(vertical_up, arc_right2));
+ EXPECT_TRUE(TestGraph::deviatesLeft(vertical_up, bezier_right));
+ EXPECT_TRUE(TestGraph::deviatesLeft(bezier_right, arc_right2));
+ EXPECT_TRUE(TestGraph::deviatesLeft(arc_right2, arc_right1));
+ EXPECT_TRUE(TestGraph::deviatesLeft(arc_left1, arc_right1));
+ EXPECT_TRUE(TestGraph::deviatesLeft(arc_left2, arc_right1));
+
+ EXPECT_FALSE(TestGraph::deviatesLeft(arc_right1, vertical_up));
+ EXPECT_FALSE(TestGraph::deviatesLeft(arc_right1, arc_right2));
+ EXPECT_FALSE(TestGraph::deviatesLeft(vertical_up, arc_left2));
+ EXPECT_FALSE(TestGraph::deviatesLeft(arc_left2, arc_left1));
+ EXPECT_FALSE(TestGraph::deviatesLeft(arc_right1, arc_left1));
+ EXPECT_FALSE(TestGraph::deviatesLeft(arc_right1, arc_left2));
+}
+
+/** Test sorting of incidences at a vertex by the outgoing heading. */
+TEST(PlanarGraphTest, BasicAzimuthalSort)
+{
+ TestGraph graph;
+
+ // Imagine the Y-axis pointing up (as in mathematics)!
+ bool const clockwise = true;
+ unsigned const num_rays = 9;
+ unsigned edges[num_rays];
+
+ // Insert the edges randomly but store them in what we know to be the
+ // clockwise order of outgoing azimuths from the vertex at the origin.
+ edges[7] = graph.insertEdge(PTH("M -0.2, -1 L 0, 0"));
+ edges[1] = graph.insertEdge(PTH("M -1, 0.2 L 0, 0"));
+ edges[4] = graph.insertEdge(PTH("M 0, 0 L 1, 0.2"));
+ edges[6] = graph.insertEdge(PTH("M 0.1, -1 L 0, 0"));
+ edges[2] = graph.insertEdge(PTH("M 0, 0 L -0.3, 1"));
+ edges[0] = graph.insertEdge(PTH("M -1, 0 H 0"));
+ edges[5] = graph.insertEdge(PTH("M 0, 0 L 1, -0.2"));
+ edges[3] = graph.insertEdge(PTH("M 0.2, 1 L 0, 0"));
+ edges[8] = graph.insertEdge(PTH("M -1, -0.1 L 0, 0"));
+
+ // We expect the incidence to edges[0] to be the last one
+ // in the sort order so it should appear first when going clockwise.
+ auto [origin, incidence] = graph.getIncidence(edges[0], TestGraph::Incidence::END);
+ ASSERT_TRUE(origin);
+ ASSERT_TRUE(incidence);
+
+ // Expect ±pi as the azimuth
+ EXPECT_DOUBLE_EQ(std::abs(incidence->azimuth), M_PI);
+
+ // Test sort order
+ for (unsigned i = 0; i < num_rays; i++) {
+ EXPECT_EQ(incidence->index, edges[i]);
+ incidence = (TestGraph::Incidence *)&graph.nextIncidence(*origin, *incidence, clockwise);
+ }
+}
+
+/** Test retrieval of a path inserted as an edge in both orientations. */
+TEST(PlanarGraphTest, PathRetrieval)
+{
+ TestGraph graph;
+
+ Path const path = PTH("M 0,0 L 1,1 C 2,2 4,2 5,1");
+ Path const htap = path.reversed();
+
+ auto edge = graph.insertEdge(path);
+
+ ASSERT_EQ(graph.numEdges(), 1);
+
+ auto [start_point, start_incidence] = graph.getIncidence(edge, TestGraph::Incidence::START);
+ ASSERT_TRUE(start_point);
+ ASSERT_TRUE(start_incidence);
+ EXPECT_EQ(graph.getOutgoingPath(start_incidence), path);
+ EXPECT_EQ(graph.getIncomingPath(start_incidence), htap);
+
+ auto [end_point, end_incidence] = graph.getIncidence(edge, TestGraph::Incidence::END);
+ ASSERT_TRUE(end_point);
+ ASSERT_TRUE(end_incidence);
+ EXPECT_EQ(graph.getIncomingPath(end_incidence), path);
+ EXPECT_EQ(graph.getOutgoingPath(end_incidence), htap);
+}
+
+/** Make sure the edge labels are correctly stored. */
+TEST(PlanarGraphTest, LabelRetrieval)
+{
+ TestGraph graph;
+ TestLabel label;
+
+ label.reversal_count = 420;
+ label.merge_count = 69;
+ label.detachment_count = 111;
+
+ auto edge = graph.insertEdge(PTH("M 0,0 L 1,1"), std::move(label));
+
+ auto retrieved = graph.getEdge(edge).label;
+ EXPECT_EQ(retrieved.reversal_count, 420);
+ EXPECT_EQ(retrieved.merge_count, 69);
+ EXPECT_EQ(retrieved.detachment_count, 111);
+}
+
+/** Regularization of duplicate edges. */
+TEST(PlanarGraphTest, MergeDuplicate)
+{
+ char const *const d = "M 2, 3 H 0 C 1,4 1,5 0,6 H 10 L 8, 0";
+ char const *const near_d = "M 2.0009,3 H 0 C 1,4 1,5 0,6 H 10.0009 L 8, 0.0005";
+
+ // Test removal of perfect overlap:
+ TestGraph graph;
+ graph.insertEdge(PTH(d));
+ graph.insertEdge(PTH(d)); // exact duplicate
+ graph.regularize();
+
+ EXPECT_TRUE(graph.isRegularized());
+
+ auto remaining = extract_labels(graph);
+
+ // Expect there to be only 1 edge after regularization.
+ ASSERT_EQ(remaining.size(), 1);
+
+ EXPECT_EQ(remaining[0].merge_count, 1); // expect one merge,
+ EXPECT_EQ(remaining[0].reversal_count, 0); // no reversals,
+ EXPECT_EQ(remaining[0].detachment_count, 0); // no detachments.
+
+ // Test removal of imperfect overlaps within numerical precision
+ TestGraph fuzzy{1e-3};
+ fuzzy.insertEdge(PTH(d));
+ fuzzy.insertEdge(PTH(near_d));
+ fuzzy.regularize();
+
+ EXPECT_TRUE(fuzzy.isRegularized());
+
+ auto fuzmaining = extract_labels(fuzzy);
+ ASSERT_EQ(fuzmaining.size(), 1);
+
+ EXPECT_EQ(fuzmaining[0].merge_count, 1); // expect one merge,
+ EXPECT_EQ(fuzmaining[0].reversal_count, 0); // no reversals,
+ EXPECT_EQ(fuzmaining[0].detachment_count, 0); // no detachments.
+
+ // Test overlap of edges with oppositie orientations.
+ TestGraph twoway;
+ twoway.insertEdge(PTH(d));
+ twoway.insertEdge(REV(d));
+ twoway.regularize();
+
+ EXPECT_TRUE(twoway.isRegularized());
+
+ auto left = extract_labels(twoway);
+ ASSERT_EQ(left.size(), 1);
+
+ EXPECT_EQ(left[0].merge_count, 1); // expect one merge,
+ EXPECT_TRUE(left[0].reversal_count == 0 || left[0].reversal_count == 1); // 0 or 1 reversals
+ EXPECT_EQ(left[0].detachment_count, 0); // no detachments.
+}
+
+/** Regularization of a shorter edge overlapping a longer one. */
+TEST(PlanarGraphTest, MergePartial)
+{
+ TestGraph graph;
+ auto longer = graph.insertEdge(PTH("M 0, 0 L 10, 10"));
+ auto shorter = graph.insertEdge(PTH("M 0, 0 L 6, 6"));
+
+ EXPECT_EQ(graph.numVertices(), 3);
+
+ graph.regularize();
+
+ EXPECT_EQ(graph.numVertices(), 3);
+ EXPECT_TRUE(graph.isRegularized());
+
+ auto labels = extract_labels(graph);
+ ASSERT_EQ(labels.size(), 2);
+
+ EXPECT_EQ(labels[longer].merge_count, 0);
+ EXPECT_EQ(labels[longer].reversal_count, 0);
+ EXPECT_EQ(labels[longer].detachment_count, 0);
+
+ EXPECT_EQ(labels[shorter].merge_count, 1);
+ EXPECT_EQ(labels[shorter].reversal_count, 0);
+ EXPECT_EQ(labels[shorter].detachment_count, 0);
+
+ // Now the same thing but with edges of opposite orientations.
+ TestGraph graphopp;
+ longer = graphopp.insertEdge(PTH("M 0, 0 L 10, 0"));
+ shorter = graphopp.insertEdge(PTH("M 10, 0 L 5, 0"));
+
+ EXPECT_EQ(graphopp.numVertices(), 3);
+
+ graphopp.regularize();
+
+ EXPECT_EQ(graphopp.numVertices(), 3);
+ EXPECT_TRUE(graphopp.isRegularized());
+
+ labels = extract_labels(graphopp);
+ ASSERT_EQ(labels.size(), 2);
+
+ EXPECT_EQ(labels[longer].merge_count, 0);
+ EXPECT_EQ(labels[longer].reversal_count, 0);
+ EXPECT_EQ(labels[longer].detachment_count, 0);
+
+ EXPECT_EQ(labels[shorter].merge_count, 1);
+ EXPECT_EQ(labels[shorter].reversal_count, 0);
+ EXPECT_EQ(labels[shorter].detachment_count, 0);
+}
+
+/** Regularization of a Y-split. */
+TEST(PlanarGraphTest, MergeY)
+{
+ TestGraph graph;
+ auto left = graph.insertEdge(PTH("M 1 0 V 1 L 0, 2"));
+ auto right = graph.insertEdge(PTH("M 1,0 V 1 L 2, 2"));
+
+ EXPECT_EQ(graph.numVertices(), 3);
+ graph.regularize();
+ EXPECT_EQ(graph.numVertices(), 4);
+
+ auto edges = graph.getEdges();
+ EXPECT_EQ(edges.size(), 3);
+
+ EXPECT_TRUE(are_near(edges[right].start->point(), Point(1, 1)));
+}
+
+/** Test reversal of a wrongly oriented teardrop */
+TEST(PlanarGraphTest, Teardrop)
+{
+ TestGraph graph;
+ auto loop = graph.insertEdge(PTH("M 1,0 A 1,1, 0,0,1 0,1 L 2,2 V 1 H 1 V 0"));
+ // Insert a few unrelated edges
+ auto before = graph.insertEdge(PTH("M 1,0 H 10"));
+ auto after = graph.insertEdge(PTH("M 1,0 H -10"));
+
+ EXPECT_EQ(graph.numVertices(), 3);
+
+ graph.regularize();
+
+ EXPECT_EQ(graph.numVertices(), 3);
+ auto [start_vertex, start_incidence] = graph.getIncidence(loop, TestGraph::Incidence::START);
+ auto [end_vertex, end_incidence] = graph.getIncidence(loop, TestGraph::Incidence::END);
+
+ EXPECT_EQ(start_vertex, end_vertex);
+ ASSERT_NE(start_vertex, nullptr);
+
+ // Check that the incidences have been swapped
+ EXPECT_EQ(start_vertex->cyclicNextIncidence(end_incidence), start_incidence);
+ EXPECT_EQ(start_vertex->cyclicPrevIncidence(start_incidence), end_incidence);
+ auto [b, before_incidence] = graph.getIncidence(before, TestGraph::Incidence::START);
+ EXPECT_EQ(start_vertex->cyclicNextIncidence(before_incidence), end_incidence);
+ auto [a, after_incidence] = graph.getIncidence(after, TestGraph::Incidence::START);
+ EXPECT_EQ(start_vertex->cyclicPrevIncidence(after_incidence), start_incidence);
+}
+
+/** Test the regularization of a lasso-shaped path. */
+TEST(PlanarGraphTest, ReglueLasso)
+{
+ TestGraph graph;
+ // Insert a lasso-shaped path (a teardrop with initial self-overlap).
+ auto original_lasso = graph.insertEdge(PTH("M 0,0 V 1 C 0,2 1,3 1,4 "
+ "A 1,1 0,1,1 -1,4 C -1,3 0,2 0,1 V 0"));
+ EXPECT_EQ(graph.numVertices(), 1);
+
+ graph.regularize();
+ EXPECT_EQ(graph.numVertices(), 2);
+ EXPECT_EQ(graph.numEdges(false), 2);
+ EXPECT_TRUE(graph.getEdge(original_lasso).detached);
+
+ auto const &edges = graph.getEdges();
+ // Find the edge from origin and ensure it got glued.
+ auto from_origin = std::find_if(edges.begin(), edges.end(), [](auto const &edge) -> bool {
+ return !edge.detached && (edge.start->point() == Point(0, 0) ||
+ edge.end->point() == Point(0, 0));
+ });
+ ASSERT_NE(from_origin, edges.end());
+ ASSERT_EQ(from_origin->label.merge_count, 1);
+}
+
+/** Test the removal of a collapsed loop. */
+TEST(PlanarGraphTest, RemoveCollapsed)
+{
+ TestGraph graph;
+ // Insert a collapsed loop
+ auto collapsed = graph.insertEdge(PTH("M 0,0 L 1,1 L 0,0"));
+ ASSERT_EQ(graph.numEdges(), 1);
+ graph.regularize();
+ ASSERT_EQ(graph.numEdges(false), 0);
+ ASSERT_TRUE(graph.getEdge(collapsed).detached);
+
+ TestGraph fuzzy(1e-3);
+ // Insert a nearly collapsed loop
+ auto nearly = fuzzy.insertEdge(PTH("M 0,0 H 2 V 0.001 L 1,0 H 0"));
+ ASSERT_EQ(fuzzy.numEdges(), 1);
+ fuzzy.regularize();
+ ASSERT_EQ(fuzzy.numEdges(false), 0);
+ ASSERT_TRUE(fuzzy.getEdge(nearly).detached);
+}
+
+/** Test regularization of straddling runs. */
+TEST(PlanarGraphTest, RemoveWisp)
+{
+ TestGraph graph;
+ // Insert a horizontal segment at the origin towards positive X:
+ graph.insertEdge(PTH("M 0 0 H 1"));
+ // Insert a path with a collapsed Bézier curve towards negative X:
+ graph.insertEdge(PTH("M 0 0 C -1 0 -1 0 0 0"));
+ graph.regularize();
+
+ // Ensure that the folded Bézier is removed (and no segfault occurs).
+ EXPECT_EQ(graph.numEdges(false), 1);
+}
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
diff --git a/tests/point-test.cpp b/tests/point-test.cpp
index a996d13..16596a5 100644
--- a/tests/point-test.cpp
+++ b/tests/point-test.cpp
@@ -74,6 +74,37 @@ TEST(PointTest, Near) {
EXPECT_FALSE(are_near_rel(Point(100, 0), Point(100, 1e-2)));
}
+TEST(PointTest, Multiplicative) {
+ EXPECT_EQ(Point(2, 3) * Point(4, 5), Point(8, 15));
+ EXPECT_EQ(IntPoint(2, 3) * IntPoint(4, 5), IntPoint(8, 15));
+ EXPECT_EQ(Point(10, 11) / Point(2, 3), Point(5, 11.0 / 3.0));
+ EXPECT_EQ(IntPoint(10, 11) / IntPoint(2, 3), IntPoint(5, 11 / 3));
+}
+
+TEST(PointTest, PointCtors) {
+ Point a(2, 3);
+ EXPECT_EQ(a[X], 2);
+ EXPECT_EQ(a[Y], 3);
+
+ a.~Point();
+ new (&a) Point;
+ EXPECT_EQ(a, Point(0, 0));
+
+ a = Point(IntPoint(4, 5));
+ EXPECT_EQ(a[X], 4);
+ EXPECT_EQ(a[Y], 5);
+}
+
+TEST(PointTest, IntPointCtors) {
+ IntPoint a(2, 3);
+ EXPECT_EQ(a[X], 2);
+ EXPECT_EQ(a[Y], 3);
+
+ a.~IntPoint();
+ new (&a) IntPoint;
+ EXPECT_EQ(a, IntPoint(0, 0));
+}
+
} // end namespace Geom
/*
diff --git a/tests/self-intersections-test.cpp b/tests/self-intersections-test.cpp
new file mode 100644
index 0000000..268273f
--- /dev/null
+++ b/tests/self-intersections-test.cpp
@@ -0,0 +1,219 @@
+/** @file
+ * @brief Unit tests for PathVector::intersectSelf()
+ */
+/*
+ * Authors:
+ * Rafał Siejakowski <rs@rs-math.net>
+ *
+ * Copyright 2022 Authors
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it either under the terms of the GNU Lesser General Public
+ * License version 2.1 as published by the Free Software Foundation
+ * (the "LGPL") or, at your option, under the terms of the Mozilla
+ * Public License Version 1.1 (the "MPL"). If you do not alter this
+ * notice, a recipient may use your version of this file under either
+ * the MPL or the LGPL.
+ *
+ * You should have received a copy of the LGPL along with this library
+ * in the file COPYING-LGPL-2.1; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * You should have received a copy of the MPL along with this library
+ * in the file COPYING-MPL-1.1
+ *
+ * The contents of this file are subject to the Mozilla Public License
+ * Version 1.1 (the "License"); you may not use this file except in
+ * compliance with the License. You may obtain a copy of the License at
+ * http://www.mozilla.org/MPL/
+ *
+ * This software is distributed on an "AS IS" basis, WITHOUT WARRANTY
+ * OF ANY KIND, either express or implied. See the LGPL or the MPL for
+ * the specific language governing rights and limitations.
+ */
+
+#include <gtest/gtest.h>
+#include <2geom/pathvector.h>
+#include <2geom/svg-path-parser.h>
+
+using namespace Geom;
+
+#define PV(d) (parse_svg_path(d))
+#define PTH(d) (PV(d)[0])
+
+class PVSelfIntersections : public testing::Test
+{
+protected:
+ PathVector const _rectangle, _bowtie, _bowtie_curved, _bowtie_node, _openpath,
+ _open_closed_nonintersecting, _open_closed_intersecting, _tangential, _degenerate_segments,
+ _degenerate_closing, _degenerate_multiple;
+
+ PVSelfIntersections()
+ // A simple rectangle.
+ : _rectangle{PV("M 0,0 L 5,0 5,8 0,8 Z")}
+ // A polyline path with a self-intersection @(2,1).
+ , _bowtie{PV("M 0,0 L 4,2 V 0 L 0,2 Z")}
+ // A curved bow-tie path with a self-intersection @(10,5) between cubic Béziers.
+ , _bowtie_curved{PV("M 0,0 V 10 C 10,10 10,0 20,0 V 10 C 10,10 10,0 0,0 Z")}
+ // As above, but twice as large and the self-intersection @(20,10) happens at a node.
+ , _bowtie_node{PV("M 0,0 V 20 C 0,20 10,20 20,10 25,5 30,0 40,0 V 20 "
+ "C 30,20 25,15 20,10 10,0 0,0 0,0 Z")}
+ // An open path with no self-intersections ◠―◡
+ , _openpath{PV("M 0,0 A 10,10 0,0,1 20,0 L 40,0 Q 50,10 60,0")}
+ // A line and a square with no intersections | □
+ , _open_closed_nonintersecting{PV("M 0,0 V 20 M 10,0 V 20 H 30 V 0 Z")}
+ // A line slicing through a square; two self-intersections ⎅
+ , _open_closed_intersecting{PV("M 10,0 V 40 M 0,10 V 30 H 20 V 10 Z")}
+ // A circle whose diameter precisely coincides with the top side of a rectangle.
+ , _tangential{PV("M 0,0 A 10,10 0,0,1 20,0 A 10,10, 0,0,1 0,0 Z M 0,0 H 20 V 30 H 0 Z")}
+ // A rectangle containing degenerate segments.
+ , _degenerate_segments{PV("M 0,0 H 5 V 4 L 5,4 V 8 H 5 L 5,8 H 0 Z")}
+ // A rectangle with a degenerate closing segment.
+ , _degenerate_closing{PV("M 0,0 H 5 V 8 H 0 L 0,0 Z")}
+ // Multiple consecutive degenerate segments, with a degenerate closing segment in the middle.
+ , _degenerate_multiple{PV("M 0,0 L 0,0 V 0 H 0 L 5,0 V 8 H 0 L 0,0 V 0 H 0 Z")}
+ {
+ }
+};
+
+/* Ensure that no spurious intersections are returned. */
+TEST_F(PVSelfIntersections, NoSpurious)
+{
+ auto empty = PathVector();
+ EXPECT_EQ(empty.intersectSelf().size(), 0u);
+
+ auto r = _rectangle.intersectSelf();
+ EXPECT_EQ(r.size(), 0u);
+
+ auto o = _openpath.intersectSelf();
+ EXPECT_EQ(o.size(), 0u);
+
+ auto n = _open_closed_nonintersecting.intersectSelf();
+ EXPECT_EQ(n.size(), 0u);
+
+ auto d = _degenerate_segments.intersectSelf();
+ EXPECT_EQ(d.size(), 0u);
+
+ auto dc = _degenerate_closing.intersectSelf();
+ EXPECT_EQ(dc.size(), 0u);
+
+ auto dm = _degenerate_multiple.intersectSelf();
+ EXPECT_EQ(dm.size(), 0u);
+
+ auto cusp_node = PTH("M 1 3 C 12 8 42 101 86 133 C 78 168 136 83 80 64");
+ EXPECT_EQ(cusp_node.intersectSelf().size(), 0u);
+}
+
+/* Test figure-eight shaped paths */
+TEST_F(PVSelfIntersections, Bowties)
+{
+ // Simple triangular bowtie: intersection between straight lines
+ auto triangular = _bowtie.intersectSelf();
+ EXPECT_EQ(triangular.size(), 1u);
+ ASSERT_GT(triangular.size(), 0u); // To ensure access to [0]
+ EXPECT_TRUE(are_near(triangular[0].point(), Point(2, 1)));
+
+ // Curved bowtie: intersection between cubic Bézier curves
+ auto curved_intersections = _bowtie_curved.intersectSelf();
+ EXPECT_EQ(curved_intersections.size(), 1u);
+ ASSERT_GT(curved_intersections.size(), 0u);
+ EXPECT_TRUE(are_near(curved_intersections[0].point(), Point(10, 5)));
+
+ // Curved bowtie but the intersection point is a node on both paths
+ auto node_case_intersections = _bowtie_node.intersectSelf();
+ EXPECT_EQ(node_case_intersections.size(), 1u);
+ ASSERT_GT(node_case_intersections.size(), 0u);
+ EXPECT_TRUE(are_near(node_case_intersections[0].point(), Point(20, 10)));
+}
+
+/* Test intersecting an open path with a closed one */
+TEST_F(PVSelfIntersections, OpenClosed)
+{
+ // Square cut by a vertical line
+ auto open_closed = _open_closed_intersecting.intersectSelf();
+ auto const P1 = Point(10, 10);
+ auto const P2 = Point(10, 30);
+
+ ASSERT_EQ(open_closed.size(), 2u); // Prevent crash on out-of-bounds access
+ // This test doesn't care about how the intersections are ordered.
+ bool points_as_expected = (are_near(open_closed[0].point(), P1) && are_near(open_closed[1].point(), P2))
+ || (are_near(open_closed[0].point(), P2) && are_near(open_closed[1].point(), P1));
+ EXPECT_TRUE(points_as_expected);
+}
+
+/* Test some nasty, tangential crossings: a circle with a rectangle built on its diameter. */
+TEST_F(PVSelfIntersections, Tangential)
+{
+ auto circle_x_rect = _tangential.intersectSelf();
+ auto const P1 = Point(0, 0);
+ auto const P2 = Point(20, 0);
+
+ ASSERT_EQ(circle_x_rect.size(), 2u); // Prevent crash on out-of-bounds access
+ // This test doesn't care how the intersections are ordered.
+ bool points_as_expected = (are_near(circle_x_rect[0].point(), P1) && are_near(circle_x_rect[1].point(), P2))
+ || (are_near(circle_x_rect[0].point(), P2) && are_near(circle_x_rect[1].point(), P1));
+ EXPECT_TRUE(points_as_expected);
+}
+
+/* Regression test for issue https://gitlab.com/inkscape/lib2geom/-/issues/33 */
+TEST_F(PVSelfIntersections, Regression33)
+{
+ // Test case provided by Pascal Bies in the issue description.
+ auto const line = LineSegment(Point(486, 597), Point(313, 285));
+ Point const c{580.1377046525328, 325.5830744834947};
+ Point const d{289.35338528516013, 450.62476639303753};
+ auto const curve = CubicBezier(c, c, d, d);
+
+ EXPECT_EQ(curve.intersect(line).size(), 1);
+}
+
+/* Regression test for issue https://gitlab.com/inkscape/lib2geom/-/issues/46 */
+TEST_F(PVSelfIntersections, NumericalInstability)
+{
+ // Test examples provided by M.B. Fraga in the issue report.
+ auto missing_intersection = PTH("M 138 237 C 293 207 129 12 167 106 Q 205 200 309 198 z");
+ auto missing_xings = missing_intersection.intersectSelf();
+ EXPECT_EQ(missing_xings.size(), 2);
+
+ auto duplicate_intersection = PTH("M 60 280 C 60 213 236 227 158 178 S 174 306 127 310 Q 80 314 60 280 z");
+ auto const only_expected = Point(130.9693916417836, 224.587385497877);
+ auto duplicate_xings = duplicate_intersection.intersectSelf();
+ ASSERT_EQ(duplicate_xings.size(), 1);
+ EXPECT_TRUE(are_near(duplicate_xings[0].point(), only_expected));
+}
+
+/* Check various numerically challenging paths consisting of 2 cubic Béziers. */
+TEST_F(PVSelfIntersections, NumericallyChallenging)
+{
+ auto two_kinks = PTH("M 85 88 C 4 425 19 6 72 426 C 128 6 122 456 68 96");
+ EXPECT_EQ(two_kinks.intersectSelf().size(), 3);
+
+ auto omega = PTH("M 47 132 C 179 343 0 78 106 74 C 187 74 0 358 174 106");
+ EXPECT_EQ(omega.intersectSelf().size(), 0);
+
+ auto spider = PTH("M 47 132 C 203 339 0 78 106 74 C 187 74 0 358 174 106");
+ EXPECT_EQ(spider.intersectSelf().size(), 4);
+
+ auto egret = PTH("M 38 340 C 183 141 16 76 255 311 C 10 79 116 228 261 398");
+ EXPECT_EQ(egret.intersectSelf().size(), 0);
+}
+
+/* Test a regression from 88040ea2aeab8ccec2b0e96c7bda2fc7d500d5ec */
+TEST_F(PVSelfIntersections, BigonFiltering)
+{
+ auto const lens = PTH("M 0,0 C 2,1 3,1 5,0 A 2.5,1 0 1 0 0,0 Z");
+ auto const xings = lens.intersectSelf();
+ // This is a simple closed path, so we expect that no self-intersections are reported.
+ EXPECT_EQ(xings.size(), 0);
+}
+
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:fileencoding=utf-8:textwidth=99 :
Debdiff
Debdiff is too long (more than 200 lines). Download the raw debdiff.