Curves
DelaunayTriangulation.AbstractParametricCurve
— Typeabstract type AbstractParametricCurve <: Function end
Abstract type for representing a parametric curve parametrised over 0 ≤ t ≤ 1
. The curves represented by this abstract type should not be self-intersecting, with the exception of allowing for closed curves.
The structs that subtype this abstract type must implement are:
differentiate
.twice_differentiate
.thrice_differentiate
(only if you have not manually definedtotal_variation
).- The struct must be callable so that
c(t)
, wherec
an instance of the struct, returns the associated value of the curve att
. - If the struct does not implement
point_position_relative_to_curve
, then the struct must implementget_closest_point
. Alternatively, rather than implementingget_closest_point
, the struct should have alookup_table
field as aVector{NTuple{2,Float64}}
, which returns values on the curve at a set of points, wherelookup_table[i]
is the value of the curve att = (i - 1) / (length(lookup_table) - 1)
.
Functions that are defined for all AbstractParametricCurve
subtypes are:
The curves in this package evaluate the total variation not by evaluating the integral itself, but by taking care of the changes in orientation in the curve to efficiently compute it. This is done by using the orientation markers of the curves, obtained using orientation_markers
, that stored in the field orientation_markers
of these curves. The function marked_total_variation
is then used to evaluate it. You may like to consider using these functions for any curve you wish to implement yourself, using e.g. the BezierCurve
struct's implementation as a reference.
DelaunayTriangulation.LineSegment
— TypeLineSegment <: AbstractParametricCurve
Curve for representing a line segment, parametrised over 0 ≤ t ≤ 1
. This curve can be using line_segment(t)
and returns a tuple (x, y)
of the coordinates of the point on the curve at t
.
Fields
first::NTuple{2,Float64}
: The first point of the line segment.last::NTuple{2,Float64}
: The last point of the line segment.length::Float64
: The length of the line segment.
Constructor
You can construct a LineSegment
using
LineSegment(first, last)
DelaunayTriangulation.CircularArc
— TypeCircularArc <: AbstractParametricCurve
Curve for representing a circular arc, parametrised over 0 ≤ t ≤ 1
. This curve can be evaluated using circular_arc(t)
and returns a tuple (x, y)
of the coordinates of the point on the curve at t
.
Fields
center::NTuple{2,Float64}
: The center of the arc.radius::Float64
: The radius of the arc.start_angle::Float64
: The angle of the initial point of the arc, in radians.sector_angle::Float64
: The angle of the sector of the arc, in radians. This is given byend_angle - start_angle
, whereend_angle
is the angle atlast
, and so might be negative for negatively oriented arcs.first::NTuple{2,Float64}
: The first point of the arc.last::NTuple{2,Float64}
: The last point of the arc.pqr::NTuple{3, NTuple{2, Float64}}
: Three points on the circle through the arc. This is needed forpoint_position_relative_to_curve
.
The angles start_angle
and end_angle
should be setup such that start_angle > end_angle
implies a positively oriented arc, and start_angle < end_angle
implies a negatively oriented arc. Moreover, they must be in [0°, 2π°)
.
Constructor
You can construct a CircularArc
using
CircularArc(first, last, center; positive=true)
It is up to you to ensure that first
and last
are equidistant from center
- the radius used will be the distance between center
and first
. The positive
keyword argument is used to determine if the arc is positively oriented or negatively oriented.
DelaunayTriangulation.EllipticalArc
— TypeEllipticalArc <: AbstractParametricCurve
Curve for representing an elliptical arc, parametrised over 0 ≤ t ≤ 1
. This curve can be evaluated using elliptical_arc(t)
and returns a tuple (x, y)
of the coordinates of the point on the curve at t
.
Fields
center::NTuple{2,Float64}
: The center of the ellipse.horz_radius::Float64
: The horizontal radius of the ellipse.vert_radius::Float64
: The vertical radius of the ellipse.rotation_scales::NTuple{2,Float64}
: Ifθ
is the angle of rotation of the ellipse, then this is(sin(θ), cos(θ))
.start_angle::Float64
: The angle of the initial point of the arc measured fromcenter
, in radians. This angle is measured from the center prior to rotating the ellipse.sector_angle::Float64
: The angle of the sector of the arc, in radians. This is given byend_angle - start_angle
, whereend_angle
is the angle atlast
, and so might be negative for negatively oriented arcs.first::NTuple{2,Float64}
: The first point of the arc.last::NTuple{2,Float64}
: The last point of the arc.
Constructor
You can construct an EllipticalArc
using
EllipticalArc(first, last, center, major_radius, minor_radius, rotation; positive=true)
where rotation
is the angle of rotation of the ellipse, in degrees. The positive
keyword argument is used to determine if the arc is positively oriented or negatively oriented.
DelaunayTriangulation.BezierCurve
— TypeBezierCurve <: AbstractParametricCurve
Curve for representing a Bezier curve, parametrised over 0 ≤ t ≤ 1
. This curve can be evaluated using bezier_curve(t)
and returns a tuple (x, y)
of the coordinates of the point on the curve at t
.
A good reference on Bezier curves is this.
See also BSpline
and CatmullRomSpline
.
This curve is only tested on loop-free curves (and closed curves that otherwise have no self-intersections). It is not guaranteed to work on curves with loops, especially for finding the nearest point on the curve to a given point.
Remember that Bezier curves are not interpolation curves. They only go through the first and last control points, but not the intermediate ones. If you want an interpolation curve, use CatmullRomSpline
.
Fields
control_points::Vector{NTuple{2,Float64}}
: The control points of the Bezier curve. The curve goes through the first and last control points, but not the intermediate ones.cache::Vector{NTuple{2,Float64}}
: A cache of the points on the curve. This is used to speed up evaluation of the curve using de Casteljau's algorithm.lookup_table::Vector{NTuple{2,Float64}}
: A lookup table for the Bezier curve, used for finding the point on the curve closest to a given point. Thei
th entry of the lookup table corresponds to thet
-valuei / (length(lookup_table) - 1)
.orientation_markers::Vector{Float64}
: The orientation markers of the curve. These are defined so that the orientation of the curve is monotone between any two consecutive markers. The first and last markers are always0
and1
, respectively. Seeorientation_markers
.
The cache is not thread-safe, and so you should not evaluate this curve in parallel.
Constructor
You can construct a BezierCurve
using
BezierCurve(control_points::Vector{NTuple{2,Float64}}; lookup_steps=5000, kwargs...)
The keyword argument lookup_steps=100
controls how many time points in [0, 1]
are used for the lookup table. The kwargs...
are keyword arguments passed to orientation_markers
.
DelaunayTriangulation.BSpline
— TypeBSpline <: AbstractParametricCurve
Curve for representing a BSpline, parametrised over 0 ≤ t ≤ 1
. This curve can be evaluated using b_spline(t)
and returns a tuple (x, y)
of the coordinates of the point on the curve at t
.
See also BezierCurve
and CatmullRomSpline
.
Our implementation of a BSpline is based on https://github.com/thibauts/b-spline.
This curve is only tested on loop-free curves (and closed curves that otherwise have no self-intersections). It is not guaranteed to work on curves with loops, especially for finding the nearest point on the curve to a given point.
Remember that B-spline curves are not interpolation curves. They only go through the first and last control points, but not the intermediate ones. For an interpolating spline, see CatmullRomSpline
.
Fields
control_points::Vector{NTuple{2,Float64}}
: The control points of the BSpline. The curve goes through the first and last control points, but not the intermediate ones.knots::Vector{Int}
: The knots of the BSpline. You should not modify or set this field directly (in particular, do not expect any support for non-uniform B-splines).cache::Vector{NTuple{2,Float64}}
: A cache of the points on the curve. This is used to speed up evaluation of the curve using de Boor's algorithm.lookup_table::Vector{NTuple{2,Float64}}
: A lookup table for the B-spline curve, used for finding the point on the curve closest to a given point. Thei
th entry of the lookup table corresponds to thet
-valuei / (length(lookup_table) - 1)
.orientation_markers::Vector{Float64}
: The orientation markers of the curve. These are defined so that the orientation of the curve is monotone between any two consecutive markers. The first and last markers are always0
and1
, respectively. Seeorientation_markers
.
Constructor
You can construct a BSpline
using
BSpline(control_points::Vector{NTuple{2,Float64}}; degree=3, lookup_steps=5000, kwargs...)
The keyword argument lookup_steps
is used to build the lookup table for the curve. Note that the default degree=3
corresponds to a cubic B-spline curve. The kwargs...
are keyword arguments passed to orientation_markers
.
DelaunayTriangulation.CatmullRomSpline
— TypeCatmullRomSpline <: AbstractParametricCurve
Curve for representing a Catmull-Rom spline, parametrised over 0 ≤ t ≤ 1
. This curve can be evaluated using catmull_rom_spline(t)
and returns a tuple (x, y)
of the coordinates of the point on the curve at t
.
For information on these splines, see e.g. this article and this article. Additionally, this article lists some nice properties of these splines.
This curve is only tested on loop-free curves (and closed curves that otherwise have no self-intersections). It is not guaranteed to work on curves with loops, especially for finding the nearest point on the curve to a given point.
Typically, Catmull-Rom splines are defined on segments of four control points, and drawn between the two interior control points. This creates an issue in that the first and last control points will not be joined to the spline. To overcome this, we extend the spline to the left and right during the evaluation of a spline, using the fields left
and right
defined below. The rules used for extending these points come from CatmullRom.jl, which extrapolates based on a Thiele-like cubic polynomial.
Fields
control_points::Vector{NTuple{2,Float64}}
: The control points of the Catmull-Rom spline. The curve goes through each point.knots::Vector{Float64}
: The parameter values of the Catmull-Rom spline. Thei
th entry of this vector corresponds to thet
-value associated with thei
th control point. With an alpha parameterα
, these values are given byknots[i+1] = knots[i] + dist(control_points[i], control_points[i+1])^α
, whereknots[1] = 0
, and the vector is the normalised by dividing byknots[end]
.lookup_table::Vector{NTuple{2,Float64}}
: A lookup table for the Catmull-Rom spline, used for finding the point on the curve closest to a given point. Thei
th entry of the lookup table corresponds to thet
-valuei / (length(lookup_table) - 1)
.alpha::Float64
: The alpha parameter of the Catmull-Rom spline. This controls the type of the parametrisation, wherealpha = 0
corresponds to uniform parametrisation,alpha = 1/2
corresponds to centripetal parametrisation, andalpha = 1
corresponds to chordal parametrisation. Must be in[0, 1]
. For reasons similar to what we describe fortension
below, we only supportalpha = 1/2
for now. (If you do really want to change it, use the_alpha
keyword argument in the constructor.)tension::Float64
: The tension parameter of the Catmull-Rom spline. This controls the tightness of the spline, withtension = 0
being the least tight, andtension = 1
leading to straight lines between the control points. Must be in[0, 1]
. You can not currently set this to anything except0.0
due to numerical issues with boundary refinement. (For example, equivariation splits are not possible iftension=1
since the curve is piecewise linear in that case, and fortension
very close to1
, the equivariation split is not always between the provided times. If you really want to change it, then you can use the_tension
keyword argument in the constructor - but be warned that this may lead to numerical issues and potentially infinite loops.)left::NTuple{2,Float64}
: The left extension of the spline. This is used to evaluate the spline on the first segment.right::NTuple{2,Float64}
: The right extension of the spline. This is used to evaluate the spline on the last segment.lengths::Vector{Float64}
: The lengths of the individual segments of the spline.segments::Vector{CatmullRomSplineSegment}
: The individual segments of the spline.orientation_markers::Vector{Float64}
: The orientation markers of the curve. These are defined so that the orientation of the curve is monotone between any two consecutive markers. The first and last markers are always0
and1
, respectively. Seeorientation_markers
.
Constructor
To construct a CatmullRomSpline
, use
CatmullRomSpline(control_points::Vector{NTuple{2,Float64}}; lookup_steps=5000, kwargs...)
The keyword argument lookup_steps
is used to build the lookup table for the curve, with lookup_steps
giving the number of time points in [0, 1]
used for the lookup table. The kwargs...
are keyword arguments passed to orientation_markers
.
DelaunayTriangulation.arc_length
— Functionarc_length(c::AbstractParametricCurve) -> Float64
arc_length(c::AbstractParametricCurve, t₁, t₂) -> Float64
Returns the arc length of the [AbstractParametricCurve
] c
. The second method returns the arc length in the interval [t₁, t₂]
, where 0 ≤ t₁ ≤ t₂ ≤ 1
.
DelaunayTriangulation.differentiate
— Functiondifferentiate(c::AbstractParametricCurve, t) -> NTuple{2, Float64}
Evaluates the derivative of c
at t
.
DelaunayTriangulation.twice_differentiate
— Functiontwice_differentiate(c::AbstractParametricCurve, t) -> NTuple{2, Float64}
Evaluates the second derivative of c
at t
.
DelaunayTriangulation.thrice_differentiate
— Functionthrice_differentiate(c::AbstractParametricCurve, t) -> NTuple{2, Float64}
Evaluates the third derivative of c
at t
.
DelaunayTriangulation.curvature
— Functioncurvature(c::AbstractParametricCurve, t) -> Float64
Returns the curvature of the [AbstractParametricCurve
] c
at t
.
DelaunayTriangulation.total_variation
— Functiontotal_variation(c::AbstractParametricCurve) -> Float64
total_variation(c::AbstractParametricCurve, t₁, t₂) -> Float64
Returns the total variation of a curve c
, or the subcurve over [t₁, t₂]
with 0 ≤ t₁ ≤ t₂ ≤ 1
, defined as the integral of the absolute curvature over this interval. (This is also known as the total absolute curvature.)
DelaunayTriangulation.point_position_relative_to_curve
— Functionpoint_position_relative_to_curve([kernel::AbstractPredicateKernel=AdaptiveKernel(),] e::AbstractParametricCurve, p) -> Certificate
Returns the position of the point p
relative to the curve c
. This function returns a [Certificate
]:
Left
:p
is to the left ofc
.Right
:p
is to the right ofc
.On
:p
is onc
.
The kernel
argument determines how this result is computed, and should be one of ExactKernel
, FastKernel
, and AdaptiveKernel
(the default). See the documentation for more information about these choices.
point_position_relative_to_curve([kernel::AbstractPredicateKernel=AdaptiveKernel(),] L::LineSegment, p) -> Certificate
Returns the position of p
relative to L
, returning a Certificate
:
Left
:p
is to the left ofL
.Right
:p
is to the right ofL
.On
:p
is onL
.
See also point_position_relative_to_line
.
The kernel
argument determines how this result is computed, and should be one of ExactKernel
, FastKernel
, and AdaptiveKernel
(the default). See the documentation for more information about these choices.
point_position_relative_to_curve([kernel::AbstractPredicateKernel=AdaptiveKernel(),] enricher::BoundaryEnricher, curve_index, p) -> Certificate
Returns a Certificate
which is
Left
: Ifp
is to the left of thecurve_index
th curve.Right
: Ifp
is to the right of thecurve_index
th curve.On
: Ifp
is on thecurve_index
th curve.
The kernel
argument determines how this result is computed, and should be one of ExactKernel
, FastKernel
, and AdaptiveKernel
(the default). See the documentation for more information about these choices.
DelaunayTriangulation.get_closest_point
— Functionget_closest_point(b::AbstractParametricCurve p) -> (Float64, NTuple{2,Float64})
Returns the t
-value and the associated point q
on the curve b
that is nearest to p
using a binary search. The search is done until the binary search interval is smaller than 1e-12
. This function will only work if the curve b
has a lookup table.
This function is only tested on loop-free curves. It is not guaranteed to work on curves with loops. Moreover, for this function to be accurate, you want the lookup table in b
to be sufficiently dense.
DelaunayTriangulation.get_equidistant_split
— Functionget_equidistant_split(c::AbstractParametricCurve, t₁, t₂) -> Float64
Returns a value of t
such that the arc length along c
from t₁
to t
is equal to the arc length along c
from t
to t₂
. Uses the bisection method to compute the t
-value.
get_equidistant_split(enricher::BoundaryEnricher, curve_index, t₁, t₂) -> Float64
Returns the equidistant split of the curve_index
th curve between t₁
and t₂
.
DelaunayTriangulation.get_equivariation_split
— Functionget_equivariation_split(c::AbstractParametricCurve, t₁, t₂) -> Float64, Float64
Returns a value of t
such that the total variation of c
from t₁
to t
is equal to the total variation of c
from t
to t₂
. Uses the bisection method to compute the t
-value. Also returns the new total variation of the two pieces.
get_equivariation_split(enricher::BoundaryEnricher, curve_index, t₁, t₂) -> Float64, Float64
Returns the equivariation split of the curve_index
th curve between t₁
and t₂
. Also returns the total variation of the two pieces.
DelaunayTriangulation.get_inverse
— Functionget_inverse(c::AbstractParametricCurve, p) -> Float64
Given a point p
on c
, returns the t
-value such that c(t) ≈ p
.
get_inverse(enricher::BoundaryEnricher, curve_index, q) -> Float64
Returns the inverse of the curve_index
th curve at q
.
DelaunayTriangulation.angle_between
— Functionangle_between(c₁::AbstractParametricCurve, c₂::AbstractParametricCurve) -> Float64
Given two curves c₁
and c₂
such that c₁(1) == c₂(0)
, returns the angle between the two curves, treating the interior of the curves as being left of both.
angle_between(L₁::LineSegment, L₂::LineSegment) -> Float64
Returns the angle between L₁
and L₂
, assuming that L₁.last == L₂.first
(this is not checked). For consistency with If the segments are part of some domain, then the line segments should be oriented so that the interior is to the left of both segments.
angle_between(enricher::BoundaryEnricher, curve_index1, curve_index2) -> Float64
Evaluates angle_between
on the curves with indices curve_index1
and curve_index2
in enricher
.
angle_between(p, q) -> Number
Returns the angle between the vectors p
and q
in radians, treating q
as the base. See this article. The returned angle is in [0, 2π)
.
DelaunayTriangulation.get_circle_intersection
— Functionget_circle_intersection(c::AbstractParametricCurve, t₁, t₂, r) -> (Float64, NTuple{2,Float64})
Given a circle centered at c(t₁)
with radius r
, finds the first intersection of the circle with the curve after t₁
and less than t₂
. It is assumed that such an intersection exists. The returned value is (t, q)
, where t
is the parameter value of the intersection and q
is the point of intersection.
get_circle_intersection(enricher::BoundaryEnricher, curve_index, t₁, t₂, r) -> (Float64, NTuple{2,Float64})
Finds the intersection of the curve_index
th curve with the circle centered at the curve evaluated at t₁
with radius r
. The argument t₂
defines the end of the subcurve to consider. The returned tuple is (t, p)
where t
is the parameter value of the intersection and p
is the point of intersection.
DelaunayTriangulation.orientation_markers
— Functionorientation_markers(c::AbstractParametricCurve; steps=200, iters=50, tol=1e-5) -> Vector{Float64}
Finds all orientation markers of the AbstractParametricCurve
c
. These are points t
where any of the following conditions hold (not necessarily simultaneously), letting c(t) = (x(t), y(t))
:
x'(t) = 0
y'(t) = 0
κ(t; x) = 0
, whereκ(t; x)
is the curvature of the component functionx(t)
κ(t; y) = 0
, whereκ(t; y)
is the curvature of the component functiony(t)
κ(t) = 0
, whereκ
is the curvature ofc(t)
Note that the third and fourth conditions give all the inflection points of the component functions, and similarly for the fifth condition.
See also horizontal_turning_points
, vertical_turning_points
, horizontal_inflection_points
, vertical_inflection_points
, and inflection_points
.
For curves of very high degree, such as Bezier curves with steps
control points or greater, this function might fail to return all inflection points.
Arguments
c::AbstractParametricCurve
: TheAbstractParametricCurve
.
Keyword Arguments
steps=200
: The number of equally spaced points to use for initialising Newton's method.iters=50
: How many iterations to use for Newton's method.tol=1e-5
: The tolerance used for determining if twot
-values are the same.
Output
markers::Vector{Float64}
: Thet
-values of the orientation markers ofb
. The returned vector is sorted, and also includes the endpoints0
and1
; anyt
-values outside of[0, 1]
are discarded, and multiplicity of anyt
is not considered (so thet
-values in the returned vector are unique). These values can be used to split the curve into monotone pieces, meaning the orientation is monotone. These markers also guarantee that, over any monotone piece, the orientation changes by an angle of at mostπ/2
.
DelaunayTriangulation.marked_total_variation
— Functionmarked_total_variation(b::AbstractParametricCurve, t₁, t₂)
Returns the total variation of the curve b
over the interval [t₁, t₂]
using the orientation markers of b
.