Bézier curves/Intersections: Difference between revisions
Added FreeBASIC
(Added FreeBASIC) |
|||
| (12 intermediate revisions by 6 users not shown) | |||
Line 594:
pt2 : opoint;
begin
pt0 := eval (gcurve, t0);
pt1 := ((t1 - t0 - t0) * t1 + (t0 * t0)) * gcurve.pt1;
Line 1,794 ⟶ 1,789:
(0.654983, 2.854982)</pre>
=={{header|
{{trans|Java}}
<syntaxhighlight lang="C#">
using System;
using System.Collections.Generic;
public class BezierCurveIntersection
{
public static void Main(string[] args)
{
QuadCurve vertical = new QuadCurve(new QuadSpline(-1.0, 0.0, 1.0), new QuadSpline(0.0, 10.0, 0.0));
// QuadCurve vertical represents the Bezier curve having control points (-1, 0), (0, 10) and (1, 0)
QuadCurve horizontal = new QuadCurve(new QuadSpline(2.0, -8.0, 2.0), new QuadSpline(1.0, 2.0, 3.0));
// QuadCurve horizontal represents the Bezier curve having control points (2, 1), (-8, 2) and (2, 3)
Console.WriteLine("The points of intersection are:");
List<Point> intersects = FindIntersects(vertical, horizontal);
foreach (Point intersect in intersects)
{
Console.WriteLine($"( {intersect.X,9:0.000000}, {intersect.Y,9:0.000000} )");
}
}
private static List<Point> FindIntersects(QuadCurve p, QuadCurve q)
{
List<Point> result = new List<Point>();
Stack<QuadCurve> stack = new Stack<QuadCurve>();
stack.Push(p);
stack.Push(q);
while (stack.Count > 0)
{
QuadCurve pp = stack.Pop();
QuadCurve qq = stack.Pop();
List<object> objects = TestIntersection(pp, qq);
bool accepted = (bool)objects[0];
bool excluded = (bool)objects[1];
Point intersect = (Point)objects[2];
if (accepted)
{
if (!SeemsToBeDuplicate(result, intersect))
{
result.Add(intersect);
}
}
else if (!excluded)
{
QuadCurve p0 = new QuadCurve();
QuadCurve q0 = new QuadCurve();
QuadCurve p1 = new QuadCurve();
QuadCurve q1 = new QuadCurve();
SubdivideQuadCurve(pp, 0.5, p0, p1);
SubdivideQuadCurve(qq, 0.5, q0, q1);
stack.Push(p0);
stack.Push(q0);
stack.Push(p0);
stack.Push(q1);
stack.Push(p1);
stack.Push(q0);
stack.Push(p1);
stack.Push(q1);
}
}
return result;
}
private static bool SeemsToBeDuplicate(List<Point> intersects, Point point)
{
foreach (Point intersect in intersects)
{
if (Math.Abs(intersect.X - point.X) < Spacing && Math.Abs(intersect.Y - point.Y) < Spacing)
{
return true;
}
}
return false;
}
private static List<object> TestIntersection(QuadCurve p, QuadCurve q)
{
double pxMin = Math.Min(Math.Min(p.X.C0, p.X.C1), p.X.C2);
double pyMin = Math.Min(Math.Min(p.Y.C0, p.Y.C1), p.Y.C2);
double pxMax = Math.Max(Math.Max(p.X.C0, p.X.C1), p.X.C2);
double pyMax = Math.Max(Math.Max(p.Y.C0, p.Y.C1), p.Y.C2);
double qxMin = Math.Min(Math.Min(q.X.C0, q.X.C1), q.X.C2);
double qyMin = Math.Min(Math.Min(q.Y.C0, q.Y.C1), q.Y.C2);
double qxMax = Math.Max(Math.Max(q.X.C0, q.X.C1), q.X.C2);
double qyMax = Math.Max(Math.Max(q.Y.C0, q.Y.C1), q.Y.C2);
bool accepted = false;
bool excluded = true;
Point intersect = new Point(0.0, 0.0);
if (RectanglesOverlap(pxMin, pyMin, pxMax, pyMax, qxMin, qyMin, qxMax, qyMax))
{
excluded = false;
double xMin = Math.Max(pxMin, qxMin);
double xMax = Math.Min(pxMax, pxMax);
if (xMax - xMin <= Tolerance)
{
double yMin = Math.Max(pyMin, qyMin);
double yMax = Math.Min(pyMax, qyMax);
if (yMax - yMin <= Tolerance)
{
accepted = true;
intersect = new Point(0.5 * (xMin + xMax), 0.5 * (yMin + yMax));
}
}
}
return new List<object> { accepted, excluded, intersect };
}
private static bool RectanglesOverlap(double xa0, double ya0, double xa1, double ya1,
double xb0, double yb0, double xb1, double yb1)
{
return xb0 <= xa1 && xa0 <= xb1 && yb0 <= ya1 && ya0 <= yb1;
}
private static void SubdivideQuadCurve(QuadCurve q, double t, QuadCurve u, QuadCurve v)
{
SubdivideQuadSpline(q.X, t, u.X, v.X);
SubdivideQuadSpline(q.Y, t, u.Y, v.Y);
}
// de Casteljau's algorithm
private static void SubdivideQuadSpline(QuadSpline q, double t, QuadSpline u, QuadSpline v)
{
double s = 1.0 - t;
u.C0 = q.C0;
v.C2 = q.C2;
u.C1 = s * q.C0 + t * q.C1;
v.C1 = s * q.C1 + t * q.C2;
u.C2 = s * u.C1 + t * v.C1;
v.C0 = u.C2;
}
public struct Point
{
public double X { get; }
public double Y { get; }
public Point(double x, double y)
{
X = x;
Y = y;
}
}
public class QuadSpline
{
public double C0 { get; set; }
public double C1 { get; set; }
public double C2 { get; set; }
public QuadSpline(double c0, double c1, double c2)
{
C0 = c0;
C1 = c1;
C2 = c2;
}
public QuadSpline() : this(0.0, 0.0, 0.0) { }
}
public class QuadCurve
{
public QuadSpline X { get; set; }
public QuadSpline Y { get; set; }
public QuadCurve(QuadSpline x, QuadSpline y)
{
X = x;
Y = y;
}
public QuadCurve() : this(new QuadSpline(), new QuadSpline()) { }
}
private const double Tolerance = 0.000_000_1;
private const double Spacing = 10 * Tolerance;
}
</syntaxhighlight>
{{out}}
<pre>
The points of intersection are:
( 0.654983, 2.854983 )
( -0.681025, 2.681025 )
( 0.881025, 1.118975 )
( -0.854983, 1.345017 )
</pre>
=={{header|C++}}==
<syntaxhighlight lang="c++">
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stack>
#include <tuple>
#include <vector>
const double TOLERANCE = 0.000'000'1;
const double SPACING = 10 * TOLERANCE;
typedef std::pair<double, double> point;
class quad_spline {
public:
quad_spline(double c0, double c1, double c2) : c0(c0), c1(c1), c2(c2) {};
quad_spline() : c0(0.0), c1(0.0), c2(0.0) {};
double c0, c1, c2;
};
class quad_curve {
public:
quad_curve(quad_spline x, quad_spline y) : x(x), y(y) {};
quad_curve() : x(quad_spline()), y(quad_spline()) {};
quad_spline x, y;
};
// de Casteljau's algorithm
void subdivide_quad_spline(const quad_spline& q, const double& t, quad_spline& u, quad_spline& v) {
const double s = 1.0 - t;
u.c0 = q.c0;
v.c2 = q.c2;
u.c1 = s * q.c0 + t * q.c1;
v.c1 = s * q.c1 + t * q.c2;
u.c2 = s * u.c1 + t * v.c1;
v.c0 = u.c2;
}
void subdivide_quad_curve(const quad_curve& q, const double t, quad_curve& u, quad_curve& v) {
subdivide_quad_spline(q.x, t, u.x, v.x);
subdivide_quad_spline(q.y, t, u.y, v.y);
}
bool rectangles_overlap(const double& xa0, const double& ya0, const double& xa1, const double& ya1,
const double& xb0, const double& yb0, const double& xb1, const double& yb1) {
return xb0 <= xa1 && xa0 <= xb1 && yb0 <= ya1 && ya0 <= yb1;
}
std::tuple<bool, bool, point> test_intersection(const quad_curve& p, const quad_curve& q) {
const double px_min = std::min(std::min(p.x.c0, p.x.c1), p.x.c2);
const double py_min = std::min(std::min(p.y.c0, p.y.c1), p.y.c2);
const double px_max = std::max(std::max(p.x.c0, p.x.c1), p.x.c2);
const double py_max = std::max(std::max(p.y.c0, p.y.c1), p.y.c2);
const double qx_min = std::min(std::min(q.x.c0, q.x.c1), q.x.c2);
const double qy_min = std::min(std::min(q.y.c0, q.y.c1), q.y.c2);
const double qx_max = std::max(std::max(q.x.c0, q.x.c1), q.x.c2);
const double qy_max = std::max(std::max(q.y.c0, q.y.c1), q.y.c2);
bool accepted = false;
bool excluded = true;
point intersect = std::make_pair(0.0, 0.0);
if ( rectangles_overlap(px_min, py_min, px_max, py_max, qx_min, qy_min, qx_max, qy_max) ) {
excluded = false;
const double x_min = std::max(px_min, qx_min);
const double x_max = std::min(px_max, px_max);
if ( x_max - x_min <= TOLERANCE ) {
const double y_min = std::max(py_min, qy_min);
const double y_max = std::min(py_max, qy_max);
if ( y_max - y_min <= TOLERANCE ) {
accepted = true;
intersect = std::make_pair(0.5 * ( x_min + x_max ), 0.5 * ( y_min + y_max));
}
}
}
return std::make_tuple(accepted, excluded, intersect);
}
bool seems_to_be_duplicate(const std::vector<point>& intersects, const point& pt) {
for ( point intersect : intersects ) {
if ( fabs(intersect.first - pt.first) < SPACING && fabs(intersect.second - pt.second) < SPACING ) {
return true;
}
}
return false;
}
std::vector<point> find_intersects(const quad_curve& p, const quad_curve& q) {
std::vector<point> result;
std::stack<quad_curve> stack;
stack.push(p);
stack.push(q);
while ( ! stack.empty() ) {
quad_curve pp = stack.top();
stack.pop();
quad_curve qq = stack.top();
stack.pop();
std::tuple<bool, bool, point> objects = test_intersection(pp, qq);
bool accepted, excluded;
point intersect;
std::tie(accepted, excluded, intersect) = objects;
if ( accepted ) {
if ( ! seems_to_be_duplicate(result, intersect) ) {
result.push_back(intersect);
}
} else if ( ! excluded ) {
quad_curve p0{}, q0{}, p1{}, q1{};
subdivide_quad_curve(pp, 0.5, p0, p1);
subdivide_quad_curve(qq, 0.5, q0, q1);
for ( quad_curve item : { p0, q0, p0, q1, p1, q0, p1, q1 } ) {
stack.push(item);
}
}
}
return result;
}
int main() {
quad_curve vertical(quad_spline(-1.0, 0.0, 1.0), quad_spline(0.0, 10.0, 0.0));
// QuadCurve vertical represents the Bezier curve having control points (-1, 0), (0, 10) and (1, 0)
quad_curve horizontal(quad_spline(2.0, -8.0, 2.0), quad_spline(1.0, 2.0, 3.0));
// QuadCurve horizontal represents the Bezier curve having control points (2, 1), (-8, 2) and (2, 3)
std::cout << "The points of intersection are:" << std::endl;
std::vector<point> intersects = find_intersects(vertical, horizontal);
for ( const point& intersect : intersects ) {
std::cout << "( " << std::setw(9) << std::setprecision(6) << intersect.first
<< ", " << intersect.second << " )" << std::endl;
}
}
</syntaxhighlight>
{{ out }}
<pre>
The points of intersection are:
( 0.654983, 2.85498 )
( -0.681025, 2.68102 )
( 0.881025, 1.11898 )
( -0.854983, 1.34502 )
</pre>
=={{header|D}}==
This program subdivides both curves by de Casteljau's algorithm, until only very small subdivisions with overlapping control polygons remain. (''You could use recursion instead of the workload container. With the container it is easier to terminate early, and also the program then uses only constant stack space.'')
Line 2,010 ⟶ 2,345:
(-0.681025, 2.681025)
(-0.854983, 1.345017)</pre>
=={{header|FreeBASIC}}==
===FreeBASIC Implementation 1===
{{trans|C}}
<syntaxhighlight lang="vbnet">' The control points of a planar quadratic Bézier curve form a
' triangle--called the "control polygon"--that completely contains
' the curve. Furthermore, the rectangle formed by the minimum and
' maximum x and y values of the control polygon completely contain
' the polygon, and therefore also the curve.
'
' Thus a simple method for narrowing down where intersections might
' be is: subdivide both curves until you find "small enough" regions
' where these rectangles overlap.
#define Min(a, b) iif((a) < (b), (a), (b))
#define Max(a, b) iif((a) > (b), (a), (b))
' Note that these are all mutable as we need to pass by reference.
Type Punto
x As Double
y As Double
End Type
Type QuadSpline
c0 As Double
c1 As Double
c2 As Double ' Non-parametric spline
End Type
Type QuadCurve
x As QuadSpline
y As QuadSpline ' Planar parametric spline
End Type
' Subdivision by de Casteljau's algorithm
Sub subdivideQuadSpline(q As QuadSpline, t As Double, u As QuadSpline, v As QuadSpline)
Dim As Double s = 1.0 - t
Dim As Double c0 = q.c0
Dim As Double c1 = q.c1
Dim As Double c2 = q.c2
u.c0 = c0
v.c2 = c2
u.c1 = s * c0 + t * c1
v.c1 = s * c1 + t * c2
u.c2 = s * u.c1 + t * v.c1
v.c0 = u.c2
End Sub
Sub subdivideQuadCurve(q As QuadCurve, t As Double, u As QuadCurve, v As QuadCurve)
subdivideQuadSpline(q.x, t, u.x, v.x)
subdivideQuadSpline(q.y, t, u.y, v.y)
End Sub
' It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1.
Function rectsOverlap(xa0 As Double, ya0 As Double, xa1 As Double, ya1 As Double, _
xb0 As Double, yb0 As Double, xb1 As Double, yb1 As Double) As Boolean
Return (xb0 <= xa1 And xa0 <= xb1 And yb0 <= ya1 And ya0 <= yb1)
End Function
Function max3(x As Double, y As Double, z As Double) As Double
Return Max(Max(x, y), z)
End Function
Function min3(x As Double, y As Double, z As Double) As Double
Return Min(Min(x, y), z)
End Function
' This accepts the point as an intersection if the boxes are small enough.
Sub testIntersect(p As QuadCurve, q As QuadCurve, tol As Double, _
Byref exclude As Boolean, Byref accept As Boolean, Byref intersect As Punto)
Dim As Double pxmin = min3(p.x.c0, p.x.c1, p.x.c2)
Dim As Double pymin = min3(p.y.c0, p.y.c1, p.y.c2)
Dim As Double pxmax = max3(p.x.c0, p.x.c1, p.x.c2)
Dim As Double pymax = max3(p.y.c0, p.y.c1, p.y.c2)
Dim As Double qxmin = min3(q.x.c0, q.x.c1, q.x.c2)
Dim As Double qymin = min3(q.y.c0, q.y.c1, q.y.c2)
Dim As Double qxmax = max3(q.x.c0, q.x.c1, q.x.c2)
Dim As Double qymax = max3(q.y.c0, q.y.c1, q.y.c2)
exclude = True
accept = False
If rectsOverlap(pxmin, pymin, pxmax, pymax, qxmin, qymin, qxmax, qymax) Then
exclude = False
Dim As Double xmin = Max(pxmin, qxmin)
Dim As Double xmax = Min(pxmax, qxmax)
Assert(xmax >= xmin)
If xmax - xmin <= tol Then
Dim As Double ymin = Max(pymin, qymin)
Dim As Double ymax = Min(pymax, qymax)
Assert(ymax >= ymin)
If ymax - ymin <= tol Then
accept = True
intersect.x = 0.5 * xmin + 0.5 * xmax
intersect.y = 0.5 * ymin + 0.5 * ymax
End If
End If
End If
End Sub
Function seemsToBeDuplicate(intersects() As Punto, icount As Integer, _
xy As Punto, spacing As Double) As Boolean
Dim As Boolean seemsToBeDup = False
Dim As Integer i = 0
While Not seemsToBeDup And i <> icount
Dim As Punto pt = intersects(i)
seemsToBeDup = Abs(pt.x - xy.x) < spacing And Abs(pt.y - xy.y) < spacing
i += 1
Wend
Return seemsToBeDup
End Function
Sub findIntersects(p As QuadCurve, q As QuadCurve, tol As Double, _
spacing As Double, intersects() As Punto)
Dim As Integer numIntersects = 0
Type workset
p As QuadCurve
q As QuadCurve
End Type
Dim As workset workload(64)
Dim As Integer numWorksets = 1
workload(0) = Type<workset>(p, q)
' Quit looking after having emptied the workload.
While numWorksets <> 0
Dim As workset work = workload(numWorksets-1)
numWorksets -= 1
Dim As Boolean exclude, accept
Dim As Punto intersect
testIntersect(work.p, work.q, tol, exclude, accept, intersect)
If accept Then
' To avoid detecting the same intersection twice, require some
' space between intersections.
If Not seemsToBeDuplicate(intersects(), numIntersects, intersect, spacing) Then
intersects(numIntersects) = intersect
numIntersects += 1
End If
Elseif Not exclude Then
Dim As QuadCurve p0, p1, q0, q1
subdivideQuadCurve(work.p, 0.5, p0, p1)
subdivideQuadCurve(work.q, 0.5, q0, q1)
workload(numWorksets) = Type<workset>(p0, q0)
numWorksets += 1
workload(numWorksets) = Type<workset>(p0, q1)
numWorksets += 1
workload(numWorksets) = Type<workset>(p1, q0)
numWorksets += 1
workload(numWorksets) = Type<workset>(p1, q1)
numWorksets += 1
End If
Wend
End Sub
Dim As QuadCurve p, q
p.x = Type<QuadSpline>(-1.0, 0.0, 1.0)
p.y = Type<QuadSpline>( 0.0, 10.0, 0.0)
q.x = Type<QuadSpline>( 2.0, -8.0, 2.0)
q.y = Type<QuadSpline>( 1.0, 2.0, 3.0)
Dim As Double tol = 0.0000001
Dim As Double spacing = tol * 10.0
Dim As Punto intersects(4)
findIntersects(p, q, tol, spacing, intersects())
For i As Integer = 0 To 3
Print "("; intersects(i).x; ", "; intersects(i).y; ")"
Next i
Sleep</syntaxhighlight>
{{out}}
<pre>( 0.6549834311008453, 2.854983419179916)
( 0.8810249865055084, 1.118975013494492)
(-0.6810249347860431, 2.681024998426437)
(-0.8549834191799164, 1.345016568899155)</pre>
===FreeBASIC Implementation 2===
{{trans|C}}
<syntaxhighlight lang="vbnet">' In this program, both of the curves are adaptively "flattened":
' that is, converted to a piecewise linear approximation. Then the
' problem is reduced to finding intersections of line segments.
'
' How efficient or inefficient the method is I will not try to answer.
' (And I do sometimes compute things "too often", although a really good
' optimizer might fix that.)
'
' I will use the symmetric power basis that was introduced by J. Sánchez-Reyes:
'
' J. Sánchez-Reyes, ‘The symmetric analogue of the polynomial power
' basis’, ACM Transactions on Graphics, vol 16 no 3, July 1997, page 319.
'
' J. Sánchez-Reyes, ‘Applications of the polynomial s-power basis
' in geometry processing’, ACM Transactions on Graphics, vol 19
' no 1, January 2000, page 35.
'
' Flattening a quadratic that is represented in this basis has a few
' advantages, which I will not go into here. */
Type bernstein_spline
b0 As Double
b1 As Double
b2 As Double
End Type
Type spower_spline
c0 As Double
c1 As Double
c2 As Double
End Type
Type spower_curve
x As spower_spline
y As spower_spline
End Type
' Convert a non-parametric spline from Bernstein basis to s-power.
Function bernstein_spline_to_spower(S As bernstein_spline) As spower_spline
Dim As spower_spline T
T.c0 = S.b0
T.c1 = (2 * S.b1) - S.b0 - S.b2
T.c2 = S.b2
Return T
End Function
' Compose (c0, c1, c2) with (t0, t1). This will map the portion [t0,t1] onto
' [0,1]. (To get these expressions, I did not use the general-degree methods
' described by Sánchez-Reyes, but instead used Maxima, some while ago.)
'
' This method is an alternative to de Casteljau subdivision, and can be done
' with the coefficients in any basis. Instead of breaking the spline into two
' pieces at a parameter value t, it gives you the portion lying between two
' parameter values. In general that requires two applications of de Casteljau
' subdivision. On the other hand, subdivision requires two applications of the
' following.
Function spower_spline_portion(S As spower_spline, t0 As Double, t1 As Double) As spower_spline
Dim As Double t0_t0 = t0 * t0
Dim As Double t0_t1 = t0 * t1
Dim As Double t1_t1 = t1 * t1
Dim As Double c2p1m0 = S.c2 + S.c1 - S.c0
Dim As spower_spline T
T.c0 = S.c0 + (c2p1m0 * t0) - (S.c1 * t0_t0)
T.c1 = (S.c1 * t1_t1) - (2 * S.c1 * t0_t1) + (S.c1 * t0_t0)
T.c2 = S.c0 + (c2p1m0 * t1) - (S.c1 * t1_t1)
Return T
End Function
Function spower_curve_portion(C As spower_curve, t0 As Double, t1 As Double) As spower_curve
Dim As spower_curve D
D.x = spower_spline_portion(C.x, t0, t1)
D.y = spower_spline_portion(C.y, t0, t1)
Return D
End Function
' Given a parametric curve, is it "flat enough" to have its quadratic
' terms removed?
Function flat_enough(C As spower_curve, tol As Double) As Boolean
' The degree-2 s-power polynomials are 1-t, t(1-t), t. We want to
' remove the terms in t(1-t). The maximum of t(1-t) is 1/4, reached
' at t=1/2. That accounts for the 1/8=0.125 in the following:
Dim As Double cx0 = C.x.c0
Dim As Double cx1 = C.x.c1
Dim As Double cx2 = C.x.c2
Dim As Double cy0 = C.y.c0
Dim As Double cy1 = C.y.c1
Dim As Double cy2 = C.y.c2
Dim As Double dx = cx2 - cx0
Dim As Double dy = cy2 - cy0
Dim As Double error_squared = 0.125 * ((cx1 * cx1) + (cy1 * cy1))
Dim As Double length_squared = (dx * dx) + (dy * dy)
Return (error_squared <= length_squared * tol * tol)
End Function
' Given two line segments, do they intersect? One solution to this problem is
' to use the implicitization method employed in the Maxima example, except to
' do it with linear instead of quadratic curves. That is what I do here, with
' the the roles of who gets implicitized alternated. If both ways you get as
' answer a parameter in [0,1], then the segments intersect.
Sub test_line_segment_intersection(ax0 As Double, ax1 As Double, _
ay0 As Double, ay1 As Double, bx0 As Double, bx1 As Double, _
by0 As Double, by1 As Double, Byref they_intersect As Boolean, _
Byref x As Double, Byref y As Double)
Dim As Double anumer = ((bx1 - bx0) * ay0 - (by1 - by0) * ax0 _
+ bx0 * by1 - bx1 * by0)
Dim As Double bnumer = -((ax1 - ax0) * by0 - (ay1 - ay0) * bx0 _
+ ax0 * ay1 - ax1 * ay0)
Dim As Double denom = ((ax1 - ax0) * (by1 - by0) _
- (ay1 - ay0) * (bx1 - bx0))
Dim As Double ta = anumer / denom ' Parameter of segment a.
Dim As Double tb = bnumer / denom ' Parameter of segment b.
they_intersect = (0 <= ta And ta <= 1 And 0 <= tb And tb <= 1)
If they_intersect Then
x = ((1 - ta) * ax0) + (ta * ax1)
y = ((1 - ta) * ay0) + (ta * ay1)
End If
End Sub
Function too_close(x As Double, y As Double, xs() As Double, ys() As Double, _
num_points As Integer, spacing As Double) As Boolean
Dim As Boolean tooclose = False
Dim As Integer i = 0
While Not tooclose And i <> num_points
tooclose = (Abs(x - xs(i)) < spacing And Abs(y - ys(i)) < spacing)
i += 1
Wend
Return tooclose
End Function
Sub recursion(tp0 As Double, tp1 As Double, tq0 As Double, tq1 As Double, _
P As spower_curve, Q As spower_curve, tol As Double, spacing As Double, _
max_points As Integer, xs() As Double, ys() As Double, Byref num_points As Integer)
If num_points = max_points Then
' do nothing
Elseif Not flat_enough(spower_curve_portion(P, tp0, tp1), tol) Then
Dim As Double tp_half = (0.5 * tp0) + (0.5 * tp1)
If Not flat_enough(spower_curve_portion(Q, tq0, tq1), tol) Then
Dim As Double tq_half = (0.5 * tq0) + (0.5 * tq1)
recursion(tp0, tp_half, tq0, tq_half, P, Q, tol, _
spacing, max_points, xs(), ys(), num_points)
recursion(tp0, tp_half, tq_half, tq1, P, Q, tol, _
spacing, max_points, xs(), ys(), num_points)
recursion(tp_half, tp1, tq0, tq_half, P, Q, tol, _
spacing, max_points, xs(), ys(), num_points)
recursion(tp_half, tp1, tq_half, tq1, P, Q, tol, _
spacing, max_points, xs(), ys(), num_points)
Else
recursion(tp0, tp_half, tq0, tq1, P, Q, tol, _
spacing, max_points, xs(), ys(), num_points)
recursion(tp_half, tp1, tq0, tq1, P, Q, tol, _
spacing, max_points, xs(), ys(), num_points)
End If
Elseif Not flat_enough(spower_curve_portion(Q, tq0, tq1), tol) Then
Dim As Double tq_half = (0.5 * tq0) + (0.5 * tq1)
recursion(tp0, tp1, tq0, tq_half, P, Q, tol, _
spacing, max_points, xs(), ys(), num_points)
recursion(tp0, tp1, tq_half, tq1, P, Q, tol, _
spacing, max_points, xs(), ys(), num_points)
Else
Dim As spower_curve P1 = spower_curve_portion(P, tp0, tp1)
Dim As spower_curve Q1 = spower_curve_portion(Q, tq0, tq1)
Dim As Boolean they_intersect
Dim As Double x, y
test_line_segment_intersection(P1.x.c0, P1.x.c2, _
P1.y.c0, P1.y.c2, _
Q1.x.c0, Q1.x.c2, _
Q1.y.c0, Q1.y.c2, _
they_intersect, x, y)
If they_intersect And Not too_close(x, y, xs(), ys(), num_points, spacing) Then
xs(num_points) = x
ys(num_points) = y
num_points += 1
End If
End If
End Sub
Sub find_intersections(P As spower_curve, Q As spower_curve, _
flatness_tolerance As Double, point_spacing As Double, _
max_points As Integer, xs() As Double, ys() As Double, _
Byref num_points As Integer)
num_points = 0
recursion(0, 1, 0, 1, P, Q, flatness_tolerance, point_spacing, _
max_points, xs(), ys(), num_points)
End Sub
Dim As bernstein_spline bPx = Type<bernstein_spline>(-1, 0, 1)
Dim As bernstein_spline bPy = Type<bernstein_spline>(0, 10, 0)
Dim As bernstein_spline bQx = Type<bernstein_spline>(2, -8, 2)
Dim As bernstein_spline bQy = Type<bernstein_spline>(1, 2, 3)
Dim As spower_spline Px = bernstein_spline_to_spower(bPx)
Dim As spower_spline Py = bernstein_spline_to_spower(bPy)
Dim As spower_spline Qx = bernstein_spline_to_spower(bQx)
Dim As spower_spline Qy = bernstein_spline_to_spower(bQy)
Dim As spower_curve P = Type<spower_curve>(Px, Py)
Dim As spower_curve Q = Type<spower_curve>(Qx, Qy)
Dim As Double flatness_tolerance = 0.001
Dim As Double point_spacing = 0.000001 ' Max norm minimum spacing.
Const max_points As Integer = 10
Dim As Double xs(max_points)
Dim As Double ys(max_points)
Dim As Integer num_points
find_intersections(P, Q, flatness_tolerance, point_spacing, _
max_points, xs(), ys(), num_points)
For i As Integer = 0 To num_points - 1
Print "("; xs(i); ", "; ys(i); ")"
Next i
Sleep</syntaxhighlight>
=={{header|Go}}==
Line 2,172 ⟶ 2,896:
(-0.681025, 2.681025)
(-0.854983, 1.345017)
</pre>
=={{header|Java}}==
<syntaxhighlight lang="java">
import java.util.ArrayList;
import java.util.List;
import java.util.Stack;
public final class BezierCurveIntersection {
public static void main(String[] aArgs) {
QuadCurve vertical = new QuadCurve( new QuadSpline(-1.0, 0.0, 1.0), new QuadSpline(0.0, 10.0, 0.0) );
// QuadCurve vertical represents the Bezier curve having control points (-1, 0), (0, 10) and (1, 0)
QuadCurve horizontal = new QuadCurve( new QuadSpline(2.0, -8.0, 2.0), new QuadSpline(1.0, 2.0, 3.0) );
// QuadCurve horizontal represents the Bezier curve having control points (2, 1), (-8, 2) and (2, 3)
System.out.println("The points of intersection are:");
List<Point> intersects = findIntersects(vertical, horizontal);
for ( Point intersect : intersects ) {
System.out.println(String.format("%s%9.6f%s%9.6f%s", "( ", intersect.aX, ", ", intersect.aY, " )"));
}
}
private static List<Point> findIntersects(QuadCurve aP, QuadCurve aQ) {
List<Point> result = new ArrayList<Point>();
Stack<QuadCurve> stack = new Stack<QuadCurve>();
stack.push(aP);
stack.push(aQ);
while ( ! stack.isEmpty() ) {
QuadCurve pp = stack.pop();
QuadCurve qq = stack.pop();
List<Object> objects = testIntersection(pp, qq);
final boolean accepted = (boolean) objects.get(0);
final boolean excluded = (boolean) objects.get(1);
Point intersect = (Point) objects.get(2);
if ( accepted ) {
if ( ! seemsToBeDuplicate(result, intersect) ) {
result.add(intersect);
}
} else if ( ! excluded ) {
QuadCurve p0 = new QuadCurve();
QuadCurve q0 = new QuadCurve();
QuadCurve p1 = new QuadCurve();
QuadCurve q1 = new QuadCurve();
subdivideQuadCurve(pp, 0.5, p0, p1);
subdivideQuadCurve(qq, 0.5, q0, q1);
stack.addAll(List.of( p0, q0, p0, q1, p1, q0, p1, q1 ));
}
}
return result;
}
private static boolean seemsToBeDuplicate(List<Point> aIntersects, Point aPoint) {
for ( Point intersect : aIntersects ) {
if ( Math.abs(intersect.aX - aPoint.aX) < SPACING && Math.abs(intersect.aY - aPoint.aY) < SPACING ) {
return true;
}
}
return false;
}
private static List<Object> testIntersection(QuadCurve aP, QuadCurve aQ) {
final double pxMin = Math.min(Math.min(aP.x.c0, aP.x.c1), aP.x.c2);
final double pyMin = Math.min(Math.min(aP.y.c0, aP.y.c1), aP.y.c2);
final double pxMax = Math.max(Math.max(aP.x.c0, aP.x.c1), aP.x.c2);
final double pyMax = Math.max(Math.max(aP.y.c0, aP.y.c1), aP.y.c2);
final double qxMin = Math.min(Math.min(aQ.x.c0, aQ.x.c1), aQ.x.c2);
final double qyMin = Math.min(Math.min(aQ.y.c0, aQ.y.c1), aQ.y.c2);
final double qxMax = Math.max(Math.max(aQ.x.c0, aQ.x.c1), aQ.x.c2);
final double qyMax = Math.max(Math.max(aQ.y.c0, aQ.y.c1), aQ.y.c2);
boolean accepted = false;
boolean excluded = true;
Point intersect = new Point(0.0, 0.0);
if ( rectanglesOverlap(pxMin, pyMin, pxMax, pyMax, qxMin, qyMin, qxMax, qyMax) ) {
excluded = false;
final double xMin = Math.max(pxMin, qxMin);
final double xMax = Math.min(pxMax, pxMax);
if ( xMax - xMin <= TOLERANCE ) {
final double yMin = Math.max(pyMin, qyMin);
final double yMax = Math.min(pyMax, qyMax);
if ( yMax - yMin <= TOLERANCE ) {
accepted = true;
intersect = new Point(0.5 * ( xMin + xMax), 0.5 * ( yMin + yMax));
}
}
}
return List.of( accepted, excluded, intersect );
}
private static boolean rectanglesOverlap(double aXa0, double aYa0, double aXa1, double aYa1,
double aXb0, double aYb0, double aXb1, double aYb1) {
return aXb0 <= aXa1 && aXa0 <= aXb1 && aYb0 <= aYa1 && aYa0 <= aYb1;
}
private static void subdivideQuadCurve(QuadCurve aQ, double aT, QuadCurve aU, QuadCurve aV) {
subdivideQuadSpline(aQ.x, aT, aU.x, aV.x);
subdivideQuadSpline(aQ.y, aT, aU.y, aV.y);
}
// de Casteljau's algorithm
private static void subdivideQuadSpline(QuadSpline aQ, double aT, QuadSpline aU, QuadSpline aV) {
final double s = 1.0 - aT;
aU.c0 = aQ.c0;
aV.c2 = aQ.c2;
aU.c1 = s * aQ.c0 + aT * aQ.c1;
aV.c1 = s * aQ.c1 + aT * aQ.c2;
aU.c2 = s * aU.c1 + aT * aV.c1;
aV.c0 = aU.c2;
}
private static record Point(double aX, double aY) {}
private static class QuadSpline {
public QuadSpline(double aC0, double aC1, double aC2) {
c0 = aC0; c1 = aC1; c2 = aC2;
}
public QuadSpline() {
this(0.0, 0.0, 0.0);
}
private double c0, c1, c2;
}
private static class QuadCurve {
public QuadCurve(QuadSpline aX, QuadSpline aY) {
x = aX; y = aY;
}
public QuadCurve() {
this( new QuadSpline(), new QuadSpline() );
}
private QuadSpline x, y;
}
private static final double TOLERANCE = 0.000_000_1;
private static final double SPACING = 10 * TOLERANCE;
}
</syntaxhighlight>
{{ out }}
<pre>
The points of intersection are:
( 0.654983, 2.854983 )
( -0.681025, 2.681025 )
( 0.881025, 1.118975 )
( -0.854983, 1.345017 )
</pre>
Line 2,538 ⟶ 3,420:
0.9405124667 (0.8810249334, 1.118975334) 0.05948753331 (0.8810246662, 1.118975067)
</pre>
=={{header|Julia}}==
{{trans|Go}}
<syntaxhighlight lang="julia">#= The control points of a planar quadratic Bézier curve form a
triangle--called the "control polygon"--that completely contains
the curve. Furthermore, the rectangle formed by the minimum and
maximum x and y values of the control polygon completely contain
the polygon, and therefore also the curve.
Thus a simple method for narrowing down where intersections might
be is: subdivide both curves until you find "small enough" regions
where these rectangles overlap.
=#
mutable struct Point
x::Float64
y::Float64
Point() = new(0, 0)
end
mutable struct QuadSpline # Non-parametric spline.
c0::Float64
c1::Float64
c2::Float64
QuadSpline(a = 0.0, b = 0.0, c = 0.0) = new(a, b, c)
end
mutable struct QuadCurve # Planar parametric spline.
x::QuadSpline
y::QuadSpline
QuadCurve(x = QuadSpline(), y = QuadSpline()) = new(x, y)
end
const Workset = Tuple{QuadCurve, QuadCurve}
""" Subdivision by de Casteljau's algorithm. """
function subdivide!(q::QuadSpline, t::Float64, u::QuadSpline, v::QuadSpline)
s = 1.0 - t
c0 = q.c0
c1 = q.c1
c2 = q.c2
u.c0 = c0
v.c2 = c2
u.c1 = s*c0 + t*c1
v.c1 = s*c1 + t*c2
u.c2 = s*u.c1 + t*v.c1
v.c0 = u.c2
end
function subdivide!(q::QuadCurve, t::Float64, u::QuadCurve, v::QuadCurve)
subdivide!(q.x, t, u.x, v.x)
subdivide!(q.y, t, u.y, v.y)
end
""" It is assumed that xa0 <= xa1, ya0 <= ya1, xb0 <= xb1, and yb0 <= yb1. """
function rectsoverlap(xa0, ya0, xa1, ya1, xb0, yb0, xb1, yb1)
return xb0 <= xa1 && xa0 <= xb1 && yb0 <= ya1 && ya0 <= yb1
end
""" This accepts the point as an intersection if the boxes are small enough. """
function testintersect(p::QuadCurve, q::QuadCurve, tol::Float64, intersect::Point)
pxmin = min(p.x.c0, p.x.c1, p.x.c2)
pymin = min(p.y.c0, p.y.c1, p.y.c2)
pxmax = max(p.x.c0, p.x.c1, p.x.c2)
pymax = max(p.y.c0, p.y.c1, p.y.c2)
qxmin = min(q.x.c0, q.x.c1, q.x.c2)
qymin = min(q.y.c0, q.y.c1, q.y.c2)
qxmax = max(q.x.c0, q.x.c1, q.x.c2)
qymax = max(q.y.c0, q.y.c1, q.y.c2)
exclude = true
accept = false
if rectsoverlap(pxmin, pymin, pxmax, pymax, qxmin, qymin, qxmax, qymax)
exclude = false
xmin = max(pxmin, qxmin)
xmax = min(pxmax, pxmax)
@assert xmin < xmax "Assertion failure: $xmax < $xmin"
if xmax-xmin <= tol
ymin = max(pymin, qymin)
ymax = min(pymax, qymax)
@assert ymin < ymax "Assertion failure: $ymax < $ymin"
if ymax-ymin <= tol
accept = true
intersect.x = 0.5*xmin + 0.5*xmax
intersect.y = 0.5*ymin + 0.5*ymax
end
end
end
return exclude, accept
end
""" return true if there seems to be a duplicate of xy in the intersects `isects` """
seemsduplicate(isects, xy, sp) = any(p -> abs(p.x-xy.x) < sp && abs(p.y-xy.y) < sp, isects)
""" find intersections of p and q """
function findintersects(p, q, tol, spacing)
intersects = Point[]
workload = Workset[(p, q)]
while !isempty(workload)
work = popfirst!(workload)
intersect = Point()
exclude, accept = testintersect(first(work), last(work), tol, intersect)
if accept
# To avoid detecting the same intersection twice, require some
# space between intersections.
if !seemsduplicate(intersects, intersect, spacing)
pushfirst!(intersects, intersect)
end
elseif !exclude
p0, p1, q0, q1 = QuadCurve(), QuadCurve(), QuadCurve(), QuadCurve()
subdivide!(first(work), 0.5, p0, p1)
subdivide!(last(work), 0.5, q0, q1)
pushfirst!(workload, (p0, q0), (p0, q1), (p1, q0), (p1, q1))
end
end
return intersects
end
""" test the methods """
function testintersections()
p, q = QuadCurve(), QuadCurve()
p.x = QuadSpline(-1.0, 0.0, 1.0)
p.y = QuadSpline(0.0, 10.0, 0.0)
q.x = QuadSpline(2.0, -8.0, 2.0)
q.y = QuadSpline(1.0, 2.0, 3.0)
tol = 0.0000001
spacing = tol * 10
for intersect in findintersects(p, q, tol, spacing)
println("($(intersect.x) $(intersect.y))")
end
end
testintersections()
</syntaxhighlight>{{out}}
<pre>
(0.6549834311008453 2.8549834191799164)
(0.8810249865055084 1.1189750134944916)
(-0.6810249984264374 2.681024934786043)
(-0.8549834191799164 1.3450165688991547)
</pre>
Line 4,776 ⟶ 5,799:
</pre>
=={{header|Raku}}==
{{trans|Go}}
<syntaxhighlight lang="raku" line># 20231025 Raku programming solution
class Point { has ($.x, $.y) is rw is default(0) }
class QuadSpline { has ($.c0, $.c1, $.c2) is rw is default(0) }
class QuadCurve { has QuadSpline ($.x, $.y) is rw }
class Workset { has QuadCurve ($.p, $.q) }
sub subdivideQuadSpline($q, $t) {
my $s = 1.0 - $t;
my ($c0,$c1,$c2) = do given $q { .c0, .c1, .c2 };
my $u_c1 = $s*$c0 + $t*$c1;
my $v_c1 = $s*$c1 + $t*$c2;
my $u_c2 = $s*$u_c1 + $t*$v_c1;
return ($c0, $u_c1, $u_c2), ($u_c2, $v_c1, $c2)
}
sub subdivideQuadCurve($q, $t, $u is rw, $v is rw) {
with (subdivideQuadSpline($q.x,$t),subdivideQuadSpline($q.y,$t))».List.flat {
$u=QuadCurve.new(x => QuadSpline.new(c0 => .[0],c1 => .[1],c2 => .[2]),
y => QuadSpline.new(c0 => .[6],c1 => .[7],c2 => .[8]));
$v=QuadCurve.new(x => QuadSpline.new(c0 => .[3],c1 => .[4],c2 => .[5]),
y => QuadSpline.new(c0 => .[9],c1 => .[10],c2 => .[11]))
}
}
sub seemsToBeDuplicate(@intersects, $xy, $spacing) {
my $seemsToBeDup = False;
for @intersects {
$seemsToBeDup = abs(.x - $xy.x) < $spacing && abs(.y - $xy.y) < $spacing;
last if $seemsToBeDup;
}
return $seemsToBeDup;
}
sub rectsOverlap($xa0, $ya0, $xa1, $ya1, $xb0, $yb0, $xb1, $yb1) {
return $xb0 <= $xa1 && $xa0 <= $xb1 && $yb0 <= $ya1 && $ya0 <= $yb1
}
sub testIntersect($p,$q,$tol,$exclude is rw,$accept is rw,$intersect is rw) {
my $pxmin = min($p.x.c0, $p.x.c1, $p.x.c2);
my $pymin = min($p.y.c0, $p.y.c1, $p.y.c2);
my $pxmax = max($p.x.c0, $p.x.c1, $p.x.c2);
my $pymax = max($p.y.c0, $p.y.c1, $p.y.c2);
my $qxmin = min($q.x.c0, $q.x.c1, $q.x.c2);
my $qymin = min($q.y.c0, $q.y.c1, $q.y.c2);
my $qxmax = max($q.x.c0, $q.x.c1, $q.x.c2);
my $qymax = max($q.y.c0, $q.y.c1, $q.y.c2);
($exclude, $accept) = True, False;
if rectsOverlap($pxmin,$pymin,$pxmax,$pymax,$qxmin,$qymin,$qxmax,$qymax) {
$exclude = False;
my ($xmin,$xmax) = max($pxmin, $qxmin), min($pxmax, $pxmax);
if ($xmax < $xmin) { die "Assertion failure: $xmax < $xmin\n" }
my ($ymin,$ymax) = max($pymin, $qymin), min($pymax, $qymax);
if ($ymax < $ymin) { die "Assertion failure: $ymax < $ymin\n" }
if $xmax - $xmin <= $tol and $ymax - $ymin <= $tol {
$accept = True;
given $intersect { (.x, .y) = 0.5 X* $xmin+$xmax, $ymin+$ymax }
}
}
}
sub find-intersects($p, $q, $tol, $spacing) {
my Point @intersects;
my @workload = Workset.new(p => $p, q => $q);
while my $work = @workload.pop {
my ($intersect,$exclude,$accept) = Point.new, False, False;
testIntersect($work.p, $work.q, $tol, $exclude, $accept, $intersect);
if $accept {
unless seemsToBeDuplicate(@intersects, $intersect, $spacing) {
@intersects.push: $intersect;
}
} elsif not $exclude {
my QuadCurve ($p0, $p1, $q0, $q1);
subdivideQuadCurve($work.p, 0.5, $p0, $p1);
subdivideQuadCurve($work.q, 0.5, $q0, $q1);
for $p0, $p1 X $q0, $q1 {
@workload.push: Workset.new(p => .[0], q => .[1])
}
}
}
return @intersects;
}
my $p = QuadCurve.new( x => QuadSpline.new(c0 => -1.0,c1 => 0.0,c2 => 1.0),
y => QuadSpline.new(c0 => 0.0,c1 => 10.0,c2 => 0.0));
my $q = QuadCurve.new( x => QuadSpline.new(c0 => 2.0,c1 => -8.0,c2 => 2.0),
y => QuadSpline.new(c0 => 1.0,c1 => 2.0,c2 => 3.0));
my $spacing = ( my $tol = 0.0000001 ) * 10;
.say for find-intersects($p, $q, $tol, $spacing);</syntaxhighlight>
You may [https://ato.pxeger.com/run?1=nVdfi9tGEH8_6HcYggjSRRaW07TXOC5pUwqFQloaaKA9imyvcyI6WdJKtpbjPkle8tB-ir71C_S9n6azM7valc_JXU5wutXOzG_n_47f_dlkb7v37__q2s3k7L_P_l0VmZTw0zYvW7iCi0xCGCR9DEGiIsglNHv9XotN1hVtOI3g-uSEZX7usvUvVZGXwgmuplpyldJ7dgeAF12zs_Ie4A0dBqFft81bKVpPhCFQotISNR0guyXg3zrf5WvhcMOgRp42gqsTALhUEEhYQJpMYYLbc7MZBmhGgEYE2oQFrLfwJt-JEoIajyUbyUS0EK7nYKG6P1YpcgfyFOXhEQLiIrWgwc4jp5Y8m3vSM0MmIGbQQsTSiLZrSlaNj-J_syjGTb2I-Qj8h5snx3xAjjIu0NLsWi3HK-OWfd5eQHjcexgVdF_8AaLSxOifv5Mfc9kmmyJrGRGfoFsMKiSl2Ic9LL72Ik576DbcTH6bnsfaV3qZ4nLGy9l5FBu0G4_6CNoXDu1Lh3Z2HkVzq9zuE5R77OA-d3BP7qncV56pUweXpqiexrseIinEpXy1_VZ81yHKKmtF-BxrVjRSrFqJQewVvmSVrfLyjZ_gnhzm1_dZIQXZvdk24CG4SB1IZEsZJr0ukF4lfQTPhlPg4UOmKkNVPtU6F6u2hXwzxp2zbS6xD6jG6EZr9nInmiKrwqDPdO4revdZSmv97pe0T-9-SfvL1LjAwiMPPFuQnFZbY_H3kr-VoStDV4aOSFaZVsj2B-uvMKhirKSg3RZxIPpV0a2FqacgW61E1dqvwcWjItOhqfrLvEQP4xvhkp67Jy1Su5hFQ4eo1IhdWXZl2dWYvb_Mes2e9XdE99lvQ6993WuLXlv0-gC9ViN2ZdmVZb-B7pS5E7rP_jH00AYL9zlOusO_ajrcMLWh2TBfx7lHsYo5BjH7NmafxeyLmG2MWfeYdYpcVdkc8UvQXjcs35OAjQDtGT9jj-egE7QJ7dC9UNeQhHXxETveUutcwINvpBRNm29L2GR50TXiKYwYfy8fcBFaRdgCNVJEGUWUr4hiRdjIkSLK4KvbFPEZfUV0qyAlJ6wklSGWGWTl2khNWGqgXLnOa4uPYzp3BHODu2q8glDPGLppLWCaPIHXp3zgo8C4WdEHnWiVux715E1erieuheqeAHy9Ylc40ox5yvKars3g53ucaopttkZVzIBDt0SlbwONWtOijjg79xd5ISjztRzKDPJJta0Gd1BIh8OGNhV7iU8a6aNM9sfj7DxoefoUmrJo4Qw9LKnY83PEA5IJrImOF7CuLASOdbdeb86OQ8cOjyeQVJ28eOpJeZkwxBJEIVGpctu68vQQ0X3-eFlRP6R-Qt0ljTzMY3OW9RbmlpZk8TsJ1VboyEn62rZg8HpgwXQe-8JlBHniRlbRlMV5paesCI446PCSHmUu1gDdG5hF4-EJPjw9TXDUNvMOWjg14w5u3md4YgQCSx0YrvRgR9fCp6gGswFtcjagze6pmmfnbAB77FSzE9QCQqpj3cV0F5rSk0IEp2jU_CSRmaKQ37HVzPl3nfl5Z3_m_Q8 Attempt This Online!]
=={{header|Scheme}}==
Line 6,017 ⟶ 7,138:
{{libheader|Wren-seq}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="
triangle--called the "control polygon"--that completely contains
the curve. Furthermore, the rectangle formed by the minimum and
| |||