by Paul Murrell http://orcid.org/0000-0002-3224-8858
Wednesday 14 June 2017
This document
is licensed under a Creative
Commons Attribution 4.0 International License.
This document describes
an algorithm for generating offset curves
for X-splines, where the width of the X-spline
is allowed to vary along its length.
The implementation is provided by the
grid.offsetXspline
function in the
'vwline' package for R.
X-splines (Blanc and Schlick, 1995) are a family of curves defined by control points (like a Bezier curve), with a shape parameter at each control point that allows the curve to vary between interpolation and approximation at each control point. In the diagram below, there are three X-splines, each drawn relative to five control points. The top X-spline has a shape parameter of 1 (except at the end points) so it approximates the control points. The bottom X-spline has a shape parameter of -1 (except at the end points) so it interpolates the control points. The middle X-spline has a shape parameter of 0, so it produces straight lines and sharp corners.
The R graphics system provides functions for drawing X-splines to produce smooth lines or closed shapes with smooth boundaries. The following code shows an example (used for the top X-spline in the diagram above).
library(grid)
x <- c(.1, .3, .5, .7, .9) y <- c(.4, .6, .4, .6, .4) grid.xspline(x, y, shape=1)
The R package 'vwline' provides several functions that draw curves based on X-splines, but with the width of the X-spline allowed to vary along the length of the curve. The following R code shows an example, where the width of the X-spline increases smoothly along the length of the line.
library(vwline)
w <- c(0, .1) grid.offsetXspline(x, y, w, shape=1)
In order to render a variable-width X-spline, we must be able to
calculate locations along the boundary of the variable-width X-spline.
This document focuses on the function grid.offsetXspline
from the 'vwline' package, which uses a generalised offset curve approach
to calculate the boundary
of a variable-width X-spline.
In the Offset curves Section, we establish a general definition of offset curves and the X-spline curves Section describes the equations for X-spline curves. The X-spline tangent functions Section generates expressions for the tangent to an X-spline curve and the Rendering X-spline offset curves Section describes how the tangent is used to draw an offset curve.
Following Lin and Chen, 2014, if we have a planar (two-dimensional) parametric curve \(\mathbf{r}=\mathbf{r}(t)\), \(t \in [0, 1]\), then an offset curve with variable offset is given by $$\mathbf{r}_0(t)=\mathbf{r}(t) + d(t)\mathbf{n}(t)$$ where \(\mathbf{n}(t)\) is the unit normal vector at each point on the original curve and \(d(t)\) is a function that defines the offset at any point on the original curve.
In order to find the unit normal function, we must first obtain the unit tangent function. This requires differentiating the original curve and then dividing by the magnitude of the derivative: $$\mathbf{e}(t) = \frac{\mathbf{r}'(t)}{\left\Vert\mathbf{r}'(t)\right\Vert}$$
Because the unit tangent function has fixed length, its derivative is perpendicular to it, so the unit normal function is then found by (Apostol, 2007): $$\mathbf{n}(t) = \frac{\mathbf{e}'(t)}{\left\Vert\mathbf{e}'(t)\right\Vert}$$
We can also observe that the unit normal is just a 90 degree rotation of the unit tangent, so we can obtain the unit normal function simply as: $$\begin{array}{l} n_x(t) = e_y(t) \\ n_y(t) = -e_x(t) \\ \end{array}$$
Following Blanc and Schlick, 1995, a (two-dimensional) spline is a planar parametric curve that is based on a discrete set of control points, \(P_k \in \mathbb{R}^2\), and a discrete set of blending functions, \(F_k : [0, 1] \rightarrow \mathbb{R}\), as follows: $$C(t) = \sum_{k=0}^n F_k(t)P_k, t \in [0, 1]$$
The definition of an X-spline reparameterises the curve and adds a set of shape parameters, \(s_k \in [-1, 1]\), which are used to both select and influence the blending functions. The X-spline curve is calculated piecewise for each (overlapping) set of four control points. The section of the curve between control points \(k+1\) and \(k+2\) is calculated as follows: $$C(t) = \frac{A_0(t)P_k + A_1(t)P_{k+1} + A_2(t)P_{k+2} + A_3(t)P_{k+3}}{A_0(t) + A_1(t) + A_2(t) + A_3(t)}, t \in [0, 1]$$ where the value and choice of blending functons used for \(A_0\) to \(A_3\) depend on the shape parameters \(s_{k+1}\) and \(s_{k+2}\): $$\begin{array}{l} A_0 = & \begin{cases} h(-t, -s_{k+1}), \\ f(t - s_{k+1}, -1 - s_{k+1}), \\ 0, \end{cases} & \begin{array}{l} \text{if \(s_{k+1} \lt 0\)} \\ \text{if \(s_{k+1} \ge 0\) and \(t \lt s_{k+1}\)} \\ \text{otherwise} \end{array} \\[1em] A_1 = & \begin{cases} g(1 - t, -s_{k+2}), \\ f(t - 1 - s_{k+2}, -1 - s_{k+2}), \end{cases} & \begin{array}{l} \text{if \(s_{k+2} \lt 0\)} \\ \text{otherwise} \end{array} \\[1em] A_2 = & \begin{cases} g(t, -s_{k+1}), \\ f(t + s_{k+1}, 1 + s_{k+1}), \end{cases} & \begin{array}{l} \text{if \(s_{k+1} \lt 0\)} \\ \text{otherwise} \end{array} \\[1em] A_3 = & \begin{cases} h(t - 1, -s_{k+2}), \\ f(t - 1 + s_{k+2}, 1 + s_{k+2}), \\ 0, \end{cases} & \begin{array}{l} \text{if \(s_{k+2} \lt 0\)} \\ \text{if \(s_{k+2} \ge 0\) and \(t \gt 1 - s_{k+2}\)} \\ \text{otherwise} \end{array} \end{array}$$ and the blending functions \(f\), \(g\), and \(h\) are defined as follows: $$\begin{array}{l} F(u, p) = u^3(10 - p + (2p - 15)u + (6 - p)u^2) \\ f(n, d) = F(n/d, 2d^2) \\ g(u, q) = u(q + u(2q + u(8 - 12q + u(14q - 11 + u(4 - 5q))))) \\ h(u, q) = u(q + u(2q + u^2(-2q - uq))) \\ \end{array}$$
The notation above, particularly the definition of the blending function \(f\), is also based on the C code implementation of X-splines from the XFig drawing program (Smith, 2001).
In order to generate an offset curve for an X-spline, we need to differentiate the function \(C(t)\) to find its tangent function, \(C'(t)\).
In the absence of the requisite discipline and mathematical wit, one approach to calculating the X-spline tangent function is to brute force the result via computational methods. This approach also has the benefit of directly producing output in the form of R code that can be evaluated and, because the X-spline equation is piecewise, it efficiently replicates across multiple equations.
The first step is to generate an expression for \(C(t)\) that is just in terms of \(t\) (and \(P_k\) and \(s_k\), both of which are constants for a particular X-spline). The idea is to expand the X-spline equations so that, for example ... $$A0 = h(-t, -s_{k+1})$$ ... and ... $$h(u, q) = u(q + u(2q + u^2(-2q - uq)))$$ ... becomes ... $$A0 = (-t)((-s_{k+1}) + (-t)(2(-s_{k+1}) + (-t)^2(-2(-s_{k+1}) - (-t)(-s_{k+1}))))$$
The following R code performs the above transformation using
text substitution and demonstrates that the result can be
executed (for specific values of t
and s1
).
hblend <- "u*(q + u*(2*q + u*u*(-2*q - u*q)))" A0 <- gsub("u", "(-t)", gsub("q", "(-s1)", hblend)) A0
[1] "(-t)*((-s1) + (-t)*(2*(-s1) + (-t)*(-t)*(-2*(-s1) - (-t)*(-s1))))"
eval(parse(text=A0), list(t=0.5, s1=1))
[1] 0.09375
To get a complete expression for \(C(t)\), we must repeat this sort of expansion for \(A1\), \(A2\), and \(A3\) as well, and there are 16 different scenarios to create equations for: \(s_{k+1}\) can be negative or non-negative, \(s_{k+2}\) can be negative or non-negative, \(A0\) can be zero or non-zero, and \(A3\) can be zero or non-zero. Furthermore, we need to generate both \(x\) and \(y\) versions to provide the two components of the vector function \(C(t)\). Not all of those scenarios can actually occur, but because we are programmatically generating the equations, it is easier, and costs no more, to generate all 16 scenarios.
The xsplineFunGenerator
function was written to
generate an expression for a specific scenario. The following
code demonstrates its use for the scenario where
\(s_{k+1}\) is negative,
\(s_{k+2}\) is negative,
\(A0\) is zero, and
\(A3\) is zero. The expression that is generated is the
equation for \(C_x(t)\) in that scenario (the px
terms in the expression are the x-values from the four
control points that are controlling the curve, e.g.,
px0
corresponds to the x-component of \(P_k\);
s1
and s2
correspond to
\(s_{k+1}\) and \(s_{k+2}\)).
xsplineFunGenerator("s1neg", "s2neg", "noA0", "noA3")("x")
expression(((-t) * ((-s1) + (-t) * (2 * (-s1) + (-t) * (-t) * (-2 * (-s1) - (-t) * (-s1)))) * px0 + (1 - t) * ((-s2) + (1 - t) * (2 * (-s2) + (1 - t) * (8 - 12 * (-s2) + (1 - t) * (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2)))))) * px1 + t * ((-s1) + t * (2 * (-s1) + t * (8 - 12 * (-s1) + t * (14 * (-s1) - 11 + t * (4 - 5 * (-s1)))))) * px2 + (t - 1) * ((-s2) + (t - 1) * (2 * (-s2) + (t - 1) * (t - 1) * (-2 * (-s2) - (t - 1) * (-s2)))) * px3)/((-t) * ((-s1) + (-t) * (2 * (-s1) + (-t) * (-t) * (-2 * (-s1) - (-t) * (-s1)))) + (1 - t) * ((-s2) + (1 - t) * (2 * (-s2) + (1 - t) * (8 - 12 * (-s2) + (1 - t) * (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2)))))) + t * ((-s1) + t * (2 * (-s1) + t * (8 - 12 * (-s1) + t * (14 * (-s1) - 11 + t * (4 - 5 * (-s1)))))) + (t - 1) * ((-s2) + (t - 1) * (2 * (-s2) + (t - 1) * (t - 1) * (-2 * (-s2) - (t - 1) * (-s2))))))
A full set of equations for \(C(t)\) for all scenarios can be generated with
a simple loop. A link to the complete definition of
xsplineFunGenerator
is given in the
Resources.
Now that we have an equation for \(C(t)\) in terms of \(t\),
we need to find its derivative. Again, we can take a computational
approach using the D
(deriv
) function from the
'Ryacas' package (Goedman et al., 2016), which provides an
interface to the Yacas
computer algebra system (Pinkus and Winitzki, 2002).
The xsplineTangentExpr
function was written
to take the output of xsplineFunGenerator
and differentiate
it. The following code demonstrates (the x-component of) the result from
xsplineTangentExpr
for the scenario above.
xsplineTangentExpr(xsplineFunGenerator("s1neg", "s2neg", "noA0", "noA3"))$x
(((-t) * ((-t) * ((-t) * (-t) * (-s1) - ((-t) + (-t)) * (-2 * (-s1) - (-t) * (-s1))) - (2 * (-s1) + (-t) * (-t) * (-2 * (-s1) - (-t) * (-s1)))) - ((-s1) + (-t) * (2 * (-s1) + (-t) * (-t) * (-2 * (-s1) - (-t) * (-s1))))) * px0 - ((1 - t) * ((1 - t) * ((1 - t) * ((1 - t) * (4 - 5 * (-s2)) + (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2)))) + (8 - 12 * (-s2) + (1 - t) * (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2))))) + (2 * (-s2) + (1 - t) * (8 - 12 * (-s2) + (1 - t) * (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2)))))) + ((-s2) + (1 - t) * (2 * (-s2) + (1 - t) * (8 - 12 * (-s2) + (1 - t) * (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2))))))) * px1 + (((-s1) + t * (2 * (-s1) + t * (8 - 12 * (-s1) + t * (14 * (-s1) - 11 + t * (4 - 5 * (-s1)))))) + t * ((2 * (-s1) + t * (8 - 12 * (-s1) + t * (14 * (-s1) - 11 + t * (4 - 5 * (-s1))))) + t * ((8 - 12 * (-s1) + t * (14 * (-s1) - 11 + t * (4 - 5 * (-s1)))) + t * ((14 * (-s1) - 11 + t * (4 - 5 * (-s1))) + t * (4 - 5 * (-s1)))))) * px2 + (((-s2) + (t - 1) * (2 * (-s2) + (t - 1) * (t - 1) * (-2 * (-s2) - (t - 1) * (-s2)))) + (t - 1) * ((2 * (-s2) + (t - 1) * (t - 1) * (-2 * (-s2) - (t - 1) * (-s2))) + (t - 1) * (((t - 1) + (t - 1)) * (-2 * (-s2) - (t - 1) * (-s2)) - (t - 1) * (t - 1) * (-s2)))) * px3)/((-t) * ((-s1) + (-t) * (2 * (-s1) + (-t) * (-t) * (-2 * (-s1) - (-t) * (-s1)))) + (1 - t) * ((-s2) + (1 - t) * (2 * (-s2) + (1 - t) * (8 - 12 * (-s2) + (1 - t) * (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2)))))) + t * ((-s1) + t * (2 * (-s1) + t * (8 - 12 * (-s1) + t * (14 * (-s1) - 11 + t * (4 - 5 * (-s1)))))) + (t - 1) * ((-s2) + (t - 1) * (2 * (-s2) + (t - 1) * (t - 1) * (-2 * (-s2) - (t - 1) * (-s2))))) - ((-t) * ((-s1) + (-t) * (2 * (-s1) + (-t) * (-t) * (-2 * (-s1) - (-t) * (-s1)))) * px0 + (1 - t) * ((-s2) + (1 - t) * (2 * (-s2) + (1 - t) * (8 - 12 * (-s2) + (1 - t) * (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2)))))) * px1 + t * ((-s1) + t * (2 * (-s1) + t * (8 - 12 * (-s1) + t * (14 * (-s1) - 11 + t * (4 - 5 * (-s1)))))) * px2 + (t - 1) * ((-s2) + (t - 1) * (2 * (-s2) + (t - 1) * (t - 1) * (-2 * (-s2) - (t - 1) * (-s2)))) * px3) * ((-t) * ((-t) * ((-t) * (-t) * (-s1) - ((-t) + (-t)) * (-2 * (-s1) - (-t) * (-s1))) - (2 * (-s1) + (-t) * (-t) * (-2 * (-s1) - (-t) * (-s1)))) - ((-s1) + (-t) * (2 * (-s1) + (-t) * (-t) * (-2 * (-s1) - (-t) * (-s1)))) - ((1 - t) * ((1 - t) * ((1 - t) * ((1 - t) * (4 - 5 * (-s2)) + (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2)))) + (8 - 12 * (-s2) + (1 - t) * (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2))))) + (2 * (-s2) + (1 - t) * (8 - 12 * (-s2) + (1 - t) * (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2)))))) + ((-s2) + (1 - t) * (2 * (-s2) + (1 - t) * (8 - 12 * (-s2) + (1 - t) * (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2))))))) + (((-s1) + t * (2 * (-s1) + t * (8 - 12 * (-s1) + t * (14 * (-s1) - 11 + t * (4 - 5 * (-s1)))))) + t * ((2 * (-s1) + t * (8 - 12 * (-s1) + t * (14 * (-s1) - 11 + t * (4 - 5 * (-s1))))) + t * ((8 - 12 * (-s1) + t * (14 * (-s1) - 11 + t * (4 - 5 * (-s1)))) + t * ((14 * (-s1) - 11 + t * (4 - 5 * (-s1))) + t * (4 - 5 * (-s1)))))) + (((-s2) + (t - 1) * (2 * (-s2) + (t - 1) * (t - 1) * (-2 * (-s2) - (t - 1) * (-s2)))) + (t - 1) * ((2 * (-s2) + (t - 1) * (t - 1) * (-2 * (-s2) - (t - 1) * (-s2))) + (t - 1) * (((t - 1) + (t - 1)) * (-2 * (-s2) - (t - 1) * (-s2)) - (t - 1) * (t - 1) * (-s2)))))/((-t) * ((-s1) + (-t) * (2 * (-s1) + (-t) * (-t) * (-2 * (-s1) - (-t) * (-s1)))) + (1 - t) * ((-s2) + (1 - t) * (2 * (-s2) + (1 - t) * (8 - 12 * (-s2) + (1 - t) * (14 * (-s2) - 11 + (1 - t) * (4 - 5 * (-s2)))))) + t * ((-s1) + t * (2 * (-s1) + t * (8 - 12 * (-s1) + t * (14 * (-s1) - 11 + t * (4 - 5 * (-s1)))))) + (t - 1) * ((-s2) + (t - 1) * (2 * (-s2) + (t - 1) * (t - 1) * (-2 * (-s2) - (t - 1) * (-s2)))))^2
Again, a simple loop can be used to generate tangent expressions for all X-spline scenarios.
It is easy to generate an expression for the unit tangent from
the X-spline
tangent expression by dividing by the square-root of the sum of
squares of the x-component and y-component;
the xsplineUnitTangentExpr
function can do this. We can then repeat the process,
differentiating the unit tangent to get an expression for the normal
function and then
generating an expression for the unit normal. This is provided by the
xsplineUnitNormal
function. However, these expressions
become uncomfortably large. For the scenario demonstrated above,
where the tangent expression consists of 65 lines of R code,
the unit tangent expression is 193 lines of R code and the unit normal
expression is 4,458 lines of R code.
Consequently, the implementation
in the 'vwline' package only uses the R expressions that were generated for
X-spline tangents and then divides by the tangent lengths and performs
a 90 degree rotation to obtain unit normals.
An X-spline offset curve is a continuous mathematical function. Rendering such a curve is most easily achieved by drawing a series of many short straight line segments that give the appearance of a smooth curve. More segments are used for longer curves and for curves with sharper bends. For example, the image below shows an X-spline drawn in R with alternating straight line segments coloured black and white; the segments are shorter on the sharper bends.
The offset curve is drawn by evaluating the tangent at the end of each of these straight line segments, dividing by the length of the tangent, rotating 90 degrees to get the unit normal, then multiplying by the appropriate offset. The diagram below shows the calculation of the left offset for a fixed offset (the unit normals are green lines and the left offset curve is red).
The offset curve equation consists of the original curve plus the unit normal function mulitplied by a function, \(d(t)\), that specifies that amount of offset at any point on the original curve: $$\mathbf{r}_0(t)=\mathbf{r}(t) + d(t)\mathbf{n}(t)$$
The function \(d(t)\) can in theory be any function that takes a
single numeric argument in the range 0 to 1.
In the 'vwline' package, the offset can be specified
using the widthSpline
function.
This describes the offset as an X-spline, where the x-values for
the control points are distances along the main curve and the
y-values for the control points are amounts of offset.
The following code gives an example where the offset rises smoothly
from 0 to 1cm at half-way along the main curve, then back to 0
(the main curve is shown as a white line).
x <- 1:4/5 y <- c(.2, .7, .2, .7) w <- widthSpline(unit(c(0, 1, 0), "cm")) grid.offsetXspline(x, y, w, shape=-1)
The offset curve only produces left and right offsets for the main curve, as shown below (the red lines).
If we want to produce a closed shape, we need to connect the ends of the offset curves. In the 'vwline' package, there are options for "butt", "square", "round", and "mitre" endings. The example below shows "round" ends.
The calculation of these line ends in described in detail in Murrell, 2017a.
The offset curve on the inside of a sharp bend in the main curve can form loops, as shown in the example below, where the left offset curve (the red line) forms a loop when the main curve takes a sharp left turn.
The 'vwline' package eliminates those loops by calling
polysimplify
from
the 'polyclip' package (Johnson and Baddeley, 2017) on the
final closed shape (after adding line ends). In the example below,
the left image shows offset curves that contain loops and the
right image shows the final result with round ends added and the loops
eliminated.
The offset curve is obtained by breaking the main curve into straight line segments, calculating offset points at the ends of those segments, and then joining up those offset points to create segments on the offset curve. This means that, if there is a sharp bend in the main curve and the offset is large, the segments in the offset curve may become long. This may result in a visibly non-smooth offset curve, as shown below.
This problem could be avoided if the offset curve function was evaluated directly, but this approach has not even been attempted because of the cost of evaluating the enormous offset curve expressions.
The expressions for the X-spline tangent functions, even though they are "only" 65 lines of R code long, can be slow to evaluate. For this reason, the 'vwline' package has the 'ByteCompile' option set so that the package code is byte-compiled at installation to improve performance.
This document has described
an algorithm for generating offset curves for X-splines,
where the amount of offset can vary along the length
of the X-spline. The algorithm relies on having expressions
for the X-spline tangent function and this was generated
computationally. The algorithm involves flattening the
X-spline to a series of straight line segments, evaluating
the tangent function for an X-spline at the end of each segment,
calculating the unit normal at each point,
and finally multiplying the unit normals by an offset.
The algorithm is implemented as the
grid.offsetXspline
function in the
'vwline' package for R.
The examples and discussion in this document relate to 'vwline' version 0.1-1.
This document was generated within a Docker container (see Resources below).
xsplineFunGenerator
and
xsplineTangentExpr
.
Murrell, P. (2017). Offset Curves for Variable-Width X-splines. Technical Report 2017-03, University of Auckland. [ bib ]
This document
is licensed under a Creative
Commons Attribution 4.0 International License.