# brush-strokes
This library provides a toolkit for dealing with Bézier curves and Bézier splines,
with a particular focus on computing the outline of a brush stroke in the spirit
of Knuth's [METAFONT](https://en.wikipedia.org/wiki/Metafont).
<p align="center"><img src="img/example_stroke.svg" width="400" alt="Brush stroke example"/></p>
---
## Bézier curves
This library mainly deals with quadratic and cubic Bézier curves,
with `Math.Bezier.Quadratic` and `Math.Bezier.Cubic` providing basic operations
such as differentiation, curvature, subdivision, and closest-point queries.
### Bézier splines
`Math.Bezier.Spline` contains functionality for working with Bézier splines
(made up of straight lines and quadratic/cubic Bézier curve). For example:
- subdividing a spline,
- deleting points from a spline,
- distinguishing between open and closed splines.
### Arc-length parametrisation
`Math.Bezier.ArcLength` implements several numerical integration strategies
for computing the arc length of a Bézier curve, and inverts the arc-length
map via Newton–Raphson for arc-length re-parametrisation:
| Integrator | Description |
|-------------------|-------------|
| `GaussLegendre` | Gauss–Legendre quadrature |
| `ClenshawCurtis` | Clenshaw–Curtis quadrature |
| `TanhSinh` | Tanh–sinh (double-exponential) quadrature |
| `Gravesen` | Gravesen's adaptive polygon/chord method |
These methods are also extended to Bézier splines.
## Brush stroking
`Math.Bezier.Stroke` implements brush stroking, in which:
- the base path is a Bézier spline,
- the brush shape is a closed Bézier spline whose parameters can vary as the
brush moves along the path; the parameters are themselves Bézier functions
of the base path Bézier parameter.
Brush shapes are defined in `Calligraphy.Brushes`; predefined shapes include
circular, elliptical, and teardrop brushes.
The stroking algorithm:
- Computes outline points by numerically solving the envelope equation
using (forward-mode) automatic differentiation.
(See `Math.Bezier.Stroke.EnvelopeEquation` for details.)
- Detects cusps via real root isolation algorithms implemented using
interval arithmetic.
- Splits up the base path at cuspidal points.
- Fits a cubic Bézier spline to the computed outline, using Hoschek's
least-squares curve fitting algorithm (see `Math.Bezier.Cubic.Fit`).
- Joins up adjacent stroke segments at corners.
### Cusps
The most complex (and computationally intensive) step in the above algorithm is
the computation of cusps, using real root isolation with interval arithmetic.
The cusp equation is formulated in `Math.Bezier.Stroke.EnvelopeEquation` using
implicit differentiation of the envelope equation. (There is also a secondary
cusp equation corresponding to corners of the brush, if the brush has corners.)
The `cusps` benchmark suite compares the performance of the root isolation
algorithms for cusp finding on a selection of brush strokes that arose in
testing.
## Root finding
`Math.Roots` implements a few basic 1D root-finding algorithms:
- A quadratic equation solver,
- A Newton–Raphson implementation with Armijo line search,
- Laguerre's method,
- A modified Halley M2 method by [Cordero, Ramos, Torregrosa (2020)](https://doi.org/10.1007/s10910-020-01108-3).
A very basic `n`-dimensional Newton's method is also implemented,
calling out to `eigen` for inversion of the Jacobian.
### Real root isolation
`Math.Root.Isolation` provides general-purpose real root isolation algorithms
for general multi-dimensional non-linear systems.
Available algorithms include bisection, the interval Newton method,
and narrowing methods (e.g. [Kubica (2017)](https://doi.org/10.1016/j.jpdc.2017.03.009)
and [Goldsztejn & Goualard (2010)](https://dl.acm.org/doi/10.1145/1774088.1774519)).
For interval Newton, two methods are available for inverting the interval
Jacobian:
- Gauss–Seidel step (with option to use the complete Gauss–Seidel method,
[Montanher, Domes, Schichl, Neumaier (2017)](https://link.springer.com/article/10.1007/s10543-017-0657-x)),
- (2D only) a highly-performant linear-programming approach, kindly contributed by
Walter Figueiredo Mascarenhas; see [`cbits/lp_2d.cpp`](cbits/lp_2d.cpp).
It is possible to customise the pipeline to choose in which order and when to apply the different algorithms when performing branch-and-bound style search.
### Interval arithmetic
`Math.Interval` provides rigorously rounded intervals and *n*-dimensional
interval boxes.
Three back-ends for computing correctly-rounded interval arithmetic operations
are selectable at build time via Cabal flags:
| Flag | Description |
|-------------|-------------|
| *(default)* | Use `rounded-hw` |
| `use-fma` | Use Kahan TwoSum and FMA TwoProd |
| `use-simd` | Use 128-bit SIMD registers |
### Automatic differentiation
`Math.Algebra.Dual` implements *k*-th order dual numbers for forward-mode
automatic differentiation in *n* variables, used throughout the library for
derivatives, curvature, and critical-point detection.
---
## Other utilities
### Arc-length test suite
The `test-arc-length` test executable can output relative error data to a CSV
file, useful for visualising which Béziers each integrator struggles with:
```
cabal run test-arc-length -- csv --output errors.csv INTEGRATOR [--grid-size STEP]
```
<p align="center"><img src="img/clenshaw_error_landscape.png" alt="Clenshaw–Curtis arc-length integrator error landscape"/></p>
### Arc-length benchmark suite
The `bench-arc-length` executable compares the speed and accuracy of various
arc-length integrators across a selection of quadratic and cubic Béziers.
On top of command-line output, it also writes the results to `bench_results.csv`.
These results can be visualised using the included helper Python script
`plot_bench.py` (requires `numpy` and `matplotlib`):
```
python src/arc-length/bench/plot_bench.py bench_results.csv
```
<p align="center"><img src="img/bench_quadratics.png" alt="Integrator time & error histograms – Quadratic Béziers"/></p>
<p align="center"><img src="img/bench_cubics.png" alt="Benchmark time & error histograms – Cubic Béziers"/></p>