Notifications
Clear all
Jul 07, 2022 9:43 pm
plane bezier seg intersection code some may find useful…
enjoy
--*************************************************
fn build_spline =
(
spl = splineShape();
crv = addNewSpline spl;
addKnot spl crv #bezier #curve [-40.9457,-47.6849,28.5503] [-36.5278,-55.1807,114.356] [-45.3636,-40.1891,-57.2554];
addKnot spl crv #bezier #curve [22.5917,-37.6313,-13.6633] [13.1029,-40.9434,60.2389] [32.0804,-34.3193,-87.5654];
spl.steps = 32;
updateShape spl;
spl;
)
--**********************************************
fn build_plane =
(
p = Plane length:100 width:100 pos:[-8.80521,-37.87,2.34646] dir:[0.306637,-0.0996364,0.946597] isSelected:off lengthsegs:1 widthsegs:1
)
--**********************************************
fn cuberoot x =
(
y = pow (abs x) (1.0/3.0);
if x < 0 then -y else y;
)
--*************************************************
fn solve_linear a b =
(
roots = #();
if abs (a as float) > 0.000001 then -- non Degenerate case
roots = #(-1.0 * b/a);
roots;
)
--****************************************************
fn solve_quad a b c =
(
roots = #();
if abs (a as float) < 0.000001 then -- Linear case, ax+b=0
roots = solve_linear b c;
else
(
v = 1.0 * b * b - 4.0 * a * c;
if abs v < 0.0000001 then
roots = #(-b/(2*a));
else if v > 0 then
(
temp = sqrt v;
roots = #((-b + temp)/(2 * a), (-b - temp)/(2 * a));
)
)
roots;
)
--****************************************************
fn solve_cubic a b c d =
(
roots = #();
if abs (a as float) < 0.000001 then -- quadratic case
(
roots = solve_quad b c d;
)
else -- Convert to depressed cubic t^3+pt+q = 0 (subst x = t - b/3a)
(
p = (3.0 * a * c - 1.0 * b * b)/(3.0 * a * a);
q = (2.0 * b * b * b - 9.0 * a * b * c + 27.0 * a * a * d)/(27.0 * a * a * a);
if abs p < 0.0000001 then -- p = 0 -> t^3 = -q -> t = -q^1/3
roots = #(cuberoot -q);
else if abs q < 0.0000001 then -- q = 0 -> t^3 + pt = 0 -> t(t^2+p)=0
roots = if p < 0 then #(sqrt -p, -sqrt -p) else #();
else
(
v = q*q/4 + p*p*p/27;
if abs v < 0.0000001 then -- val = 0 -> two roots
roots = #(-1.5*q/p, 3*q/p);
else if v > 0 then -- Only one real root
(
u = cuberoot (-q/2.0 - sqrt v);
roots = #(u - p/(3*u));
)
else -- val < 0, three roots, but needs to use complex numbers/trigonometric solution
(
u = 2.0 * sqrt (-p/3.0);
t = degToRad ((acos (3.0 * q/p/u))/3.0); -- val < 0 implies p < 0 and acos argument in [-1..1]
k = 2.0 * PI / 3.0;
roots = #(u * cos (radToDeg t), u * cos (radToDeg(t-k)), u * cos (radToDeg (t-2*k)));
)
)
for i = 1 to roots.count do roots[i] -= b/(3*a); -- Convert back from depressed cubic
)
roots;
)
--****************************************************
fn comp_plane_bez_seg_intersection pln spl crv seg =
(
-- get plane info
n = pln.dir
q = pln.pos
-- extract bezier info (no range checking)
p0 = getknotpoint spl crv seg;
p1 = getOutVec spl crv seg;
p2 = getinVec spl crv (seg + 1);
p3 = getknotpoint spl crv (seg + 1);
-- create the polynomial coefficients where possible intersections occur n.b(t) = n.q
a = dot n (3 * p1 - 3 * p2 + p3 - p0);
b = dot n (3 * p0 - 6 * p1 + 3 * p2);
c = dot n (3 * p1 - 3 * p0);
d = dot n (p0 - q);
-- solve for t
roots = solve_cubic a b c d;
for i in roots where i >= 0 and i <= 1 do point pos:(pathInterp spl crv i)
)
--****************************************************
delete objects
comp_plane_bez_seg_intersection (build_plane()) (build_spline()) 1 1 10 Replies
Jul 07, 2022 9:43 pm
Great job!
as soon as I have time, I will make a c++ version and add it to my MXS extension.
thanks!
Jul 07, 2022 9:43 pm
cheers, I was quite chuffed with it though trying to expand (1-x)^3 after 30 algebra free years was quite comical, thank god for the internet
also I haven’t checked what happens when the all the bezier curve is entirely in the plane
Jul 07, 2022 9:43 pm
#define is_valid_param(t) (t > 0.0f && t < 1.0f)
class CubicRootsMath
{
public:
static double cube_root(double x)
{
double y = pow(abs(x), 1.0/3) * sign(x);
return y;
}
static void linear_solver(double a, double b, double_vector& roots)
{
if (abs(a) > FLT_EPSILON)
{
auto t = -b / a;
roots = { t };
}
}
static void quadratic_solver(double a, double b, double c, double_vector& roots)
{
if (abs(a) < FLT_EPSILON)
{
linear_solver(b, c, roots);
}
else
{
auto v = 1 * b * b - 4 * a * c;
if (abs(v) < FLT_EPSILON)
{
auto t = -b / (2 * a);
roots = { t };
}
else if (v > 0)
{
auto t0 = (-b + sqrt(v)) / (2 * a);
auto t1 = (-b - sqrt(v)) / (2 * a);
roots = { t0, t1 };
}
}
}
static void cubic_solver(double a, double b, double c, double d, double_vector& roots)
{
roots = { };
if (abs(a) < FLT_EPSILON)
{
quadratic_solver(b, c, d, roots);
}
else
{
auto p = (3 * a * c - 1 * b * b) / (3 * a * a);
auto q = (2 * b * b * b - 9 * a * b * c + 27 * a * a * d) / (27 * a * a * a);
if (abs(p) < FLT_EPSILON)
{
auto t = cube_root(-q);
roots = { t };
}
else if (abs(q) < FLT_EPSILON)
{
if (p < 0)
{
auto t0 = sqrt(-p);
auto t1 = -t0;
roots = { t0, t1 };
}
}
else
{
auto v = q * q / 4 + p * p * p / 27;
if (abs(v) < FLT_EPSILON)
{
auto t0 = -1.5 * q / p;
auto t1 = 3 * q / p;
roots = { t0, t1 };
}
else if (v > 0)
{
auto u = cube_root(-q / 2 - sqrt(v));
auto t = u - p / (3 * u);
roots = { t };
}
else
{
auto u = 2 * sqrt(-p / 3);
auto t = acos(3 * q / p / u) / 3;
auto k = (double)TWOPI / 3;
auto t0 = u * cos(t);
auto t1 = u * cos(t - k);
auto t2 = u * cos(t - 2 * k);
roots = { t0, t1, t2 };
}
}
}
for(auto& it : roots) it -= b / (3 * a);
}
static void bezie_spline_plane_intersection(Point3Tab pp, Ray plane, double_vector& roots)
{
Point3 n = plane.dir;
Point3 q = plane.p;
// create the polynomial coefficients where possible intersections occur n.b(t) = n.q
double a = DotProd(n, 3 * pp[1] - 3 * pp[2] + pp[3] - pp[0]);
double b = DotProd(n, 3 * pp[0] - 6 * pp[1] + 3 * pp[2]);
double c = DotProd(n, 3 * pp[1] - 3 * pp[0]);
double d = DotProd(n, pp[0] - q);
// solve for t
cubic_solver(a, b, c, d, roots);
}
};I hope this is clear enough… thanks again. works really well, so as soon as I have time I will make a plugin “cut spline by plane”
Jul 07, 2022 9:43 pm
found some issues with this approach, “straight” beziers (lines) at acute angles to the plane don’t give very good results so perhaps testing that p1 & p2 are colinear with p0 & p3 first then it’s a trivial line plane intersection instead.
or just reduce the “tolerance” in script
Jul 07, 2022 9:43 pm
c++ version works fine with both “curve” and “linear” segments, only needs to change type of affected knots to BEZIER_CORNER
Jul 07, 2022 9:43 pm
finally:
template <class T> double double_dotproduct(T vector_a, T vector_b, int size)
{
double product(0);
for (int i = 0; i < size; i++)
product += (double)vector_a[i] * (double)vector_b[i];
return product;
}
#define dot_double3(a, b) double_dotproduct(a, b, 3)
int cubic_equation(double a, double b, double c, double* x)
{
double q, r, r2, q3;
q = (a * a - 3 * b) / 9;
r = (a * (2 * a * a - 9 * b) + 27 * c) / 54;
r2 = r * r;
q3 = q * q * q;
if (r2 < q3)
{
double t = acos(r / sqrt(q3));
a /= 3;
q = -2 * sqrt(q);
x[0] = q * cos(t / 3) - a;
x[1] = q * cos((t + TWOPI) / 3) - a;
x[2] = q * cos((t - TWOPI) / 3) - a;
return 3;
}
else
{
double p, aa, bb;
if (r <= 0) r = -r;
p = r + sqrt(r2 - q3);
aa = pow(abs(p), 1.0 / 3) * sign(p);
bb = (aa == 0) ? 0 : q / aa;
a /= 3;
q = aa + bb;
r = aa - bb;
x[0] = q - a;
x[1] = (-0.5) * q - a;
x[2] = (sqrt(3.0) * 0.5) * fabs(r);
return (x[2] == 0) ? 2 : 1;
}
}
int bezie_spline_plane_intersection(Point3Tab pp, Ray plane, double* roots)
{
Point3 n = plane.dir;
Point3 q = plane.p;
// create the polynomial coefficients where possible intersections occur n.b(t) = n.q
double a = dot_double3(n, 3 * pp[1] - 3 * pp[2] + pp[3] - pp[0]);
double b = dot_double3(n, 3 * pp[0] - 6 * pp[1] + 3 * pp[2]);
double c = dot_double3(n, 3 * pp[1] - 3 * pp[0]);
double d = dot_double3(n, pp[0] - q);
// solve for t
return cubic_equation(b / a, c / a, d / a, roots);
}
… and, of course, it is necessary to choose solutions in the range [0,1]
if it is supposed to refine (or cut) spline segments, then the solutions must also be sorted
Jul 07, 2022 9:43 pm
you’ll need to be careful as a can equal 0.0
Jul 07, 2022 9:43 pm
int quadratic_equation(double a, double b, double c, double* x)
{
double dis = b * b - 4 * a * c;
if (dis < 0) return 0;
double sqrt_d = sqrt(dis);
x[0] = (-b + sqrt_d) / (2 * a);
x[1] = (-b - sqrt_d) / (2 * a);
return 2;
}
...
if (a == 0)
return quadratic_equation(b, c, d, roots);
else
return cubic_equation(b / a, c / a, d / a, roots);how rigorous you are!
Jul 07, 2022 9:43 pm
I’ve had this bug in my library for ten years… so far, none of the clients have complained.
Page 1 / 2
Next