After extensive searches i havn't been able to find a solution.

Working in 2d, i want to find the intersection point of a line with a catmull rom spline curve.

At the moment, i'm dividing the curve segment into numerous straight lines and performing line to line intersection tests with each.

It works, but isn't very efficient and not very accurate, depending on how many sub divisions of the curve i use.

Does anybody know a much more efficient and accurate solution to this problem.

Any help would be much appreciated.

**0**

# Line To Catmull Rom Intersection

Started by K_J_M, Aug 24 2012 05:27 PM

12 replies to this topic

Sponsor:

###
#2
Crossbones+ - Reputation: **17464**

Posted 25 August 2012 - 12:08 AM

To find the intersection between a straight line and a cubic piece of a spline, write the line as Ax+By+C=0, then substitute x and y by their values as a function of the parameter t. You should get a cubic equation on t, which you can either solve in terms of radicals (tricky) or using Newton's method (easy).

If you have trouble with that plan, explain how far you got and what part you can't get to work.

If you have trouble with that plan, explain how far you got and what part you can't get to work.

###
#3
Members - Reputation: **401**

Posted 25 August 2012 - 01:33 AM

Hello Alvaro.

Thanks for the response.

I've been able to google a lot of information i wasn't aware of before thanks to your explanation.

To clarify.

Using Newton's method to find the root of (Ax+By+C=0), where X and Y is a point on the curve segment.

Taking an initial guess for X ,Y (i presume t = 0.5)

Subsequent iterations narrow down the accuracy of the curve point X , Y until the root is found.

In code terms, a quick google search provided (Ax+By+C=0) to be

(y1 – y2)x + (x2 – x1)y + (x1y2 – x2y1) = 0

If my above assumptions are correct, all thats needed now is to write Newtons Method to find the root.

Another google search provided this wikipedia page explaining Newton's Method.

http://en.wikipedia.org/wiki/Newton%27s_method

Thanks for the response.

I've been able to google a lot of information i wasn't aware of before thanks to your explanation.

To clarify.

Using Newton's method to find the root of (Ax+By+C=0), where X and Y is a point on the curve segment.

Taking an initial guess for X ,Y (i presume t = 0.5)

Subsequent iterations narrow down the accuracy of the curve point X , Y until the root is found.

In code terms, a quick google search provided (Ax+By+C=0) to be

(y1 – y2)x + (x2 – x1)y + (x1y2 – x2y1) = 0

If my above assumptions are correct, all thats needed now is to write Newtons Method to find the root.

Another google search provided this wikipedia page explaining Newton's Method.

http://en.wikipedia.org/wiki/Newton%27s_method

###
#4
Members - Reputation: **401**

Posted 25 August 2012 - 05:35 AM

I've had some success after a bit of experimentation.

I initially choose a T value of 0.5 and calculate the spline X,Y position of T.

Then test for the of root of the line (Ax+By+C=0) using the spline X,Y values.

This returns either a positive or negative number of arbitry value.

Then i iterate N number of times (depending on the accuracy required ) and either add or subtract to T a value that gets halved each time depending on the Sign of the root value.

Since T initially has a value of 0.5, the next update of T is either +- 0.25 , +-0.125 , +- 0.0625 ect

Then recalculate the spline X,Y position for the new value of T which gets plugged back into the (Ax+By+C=0) equation.

Each iteration gets closer to the solution.

I don't know if this is the correct method to use, but so far it seems to work well, except for when the line is reversed.

So, I'm assuming this special case will need to be checked for first.

I also need to run further tests to see how robust this technique is.

I initially choose a T value of 0.5 and calculate the spline X,Y position of T.

Then test for the of root of the line (Ax+By+C=0) using the spline X,Y values.

This returns either a positive or negative number of arbitry value.

Then i iterate N number of times (depending on the accuracy required ) and either add or subtract to T a value that gets halved each time depending on the Sign of the root value.

Since T initially has a value of 0.5, the next update of T is either +- 0.25 , +-0.125 , +- 0.0625 ect

Then recalculate the spline X,Y position for the new value of T which gets plugged back into the (Ax+By+C=0) equation.

Each iteration gets closer to the solution.

I don't know if this is the correct method to use, but so far it seems to work well, except for when the line is reversed.

So, I'm assuming this special case will need to be checked for first.

I also need to run further tests to see how robust this technique is.

###
#5
Crossbones+ - Reputation: **17464**

Posted 25 August 2012 - 06:35 AM

I think you are not doing it the way I described, so maybe I'll describe it in more detail.

A cubic spline is piece-wise cubic, which means that it consists of several pieces of parametric curve (each for a range of values of t) that look like this:

x = a3*t^3+a2*t^2+a1*t+a0

y = b3*t^3+b2*t^2+b1*t+b0

If you now have an equation Ax+By+C=0, you can substitute the values of x and y from the parametrization above and you'll get

(A*a3+B*b3)*t^3 + (A*a2+B*b2)*t^2 + (A*a1+B*b1)*t + A*a0+B*b0 = 0

This is a cubic equation in t. Compute its coefficients:

C3 = A*a3+B*b3

C2 = A*a2+B*b2

C1 = A*a1+B*b1

C0 = A*a0+B*b0

Now the equation is C3*t^3 + C2*t^2 + C1*t + C0 = 0. Plug this into a solver that can get you all the real roots and look for roots in the range of values of t that this piece is for. The place where I was suggesting to use Newton's method was in the solver.

Is that more clear?

A cubic spline is piece-wise cubic, which means that it consists of several pieces of parametric curve (each for a range of values of t) that look like this:

x = a3*t^3+a2*t^2+a1*t+a0

y = b3*t^3+b2*t^2+b1*t+b0

If you now have an equation Ax+By+C=0, you can substitute the values of x and y from the parametrization above and you'll get

(A*a3+B*b3)*t^3 + (A*a2+B*b2)*t^2 + (A*a1+B*b1)*t + A*a0+B*b0 = 0

This is a cubic equation in t. Compute its coefficients:

C3 = A*a3+B*b3

C2 = A*a2+B*b2

C1 = A*a1+B*b1

C0 = A*a0+B*b0

Now the equation is C3*t^3 + C2*t^2 + C1*t + C0 = 0. Plug this into a solver that can get you all the real roots and look for roots in the range of values of t that this piece is for. The place where I was suggesting to use Newton's method was in the solver.

Is that more clear?

###
#6
Members - Reputation: **401**

Posted 25 August 2012 - 12:35 PM

Hello Alvaro.

Thanks for the more detailed explanation.

That makes a lot more sense to me.

Unfortunately, i'm running into a few problems.

I need clarification whether

x = a3*t^3+a2*t^2+a1*t+a0

y = b3*t^3+b2*t^2+b1*t+b0

returns a point on a Catmull Rom Spline.

I'm currently using the equation below to return x,y points on the Catmull Rom spline for a given T value.

x = 0.5 * ((2 * p1x) + (p2x - p0x) * t + (2 * p0x - 5 * p1x + 4 * p2x - p3x) * t * t + (3 * p1x + p3x - p0x - 3 * p2x) * t * t * t)

y = 0.5 * ((2 * p1y) + (p2y - p0y) * t + (2 * p0y - 5 * p1y + 4 * p2y - p3y) * t * t + (3 * p1y + p3y - p0y - 3 * p2y) * t * t * t)

After coding your more detailed explanation i found i was getting wild results when trying to obtain the root.

Further exploration seems to reveal different spline formula's.

To verify this, i simply replaced my above spline formula's with yours to see if the same curve was drawn.

And the results were not the same.

Although both of our spline equations require 4 control points each, my spline gets drawn from point 1 to 2 ( with no drawing between points 0 and 3), and your spline was being drawn from point 0.

Which obviously leads me to ask if we are dealing with 2 different sets of spline equations ?

Thanks for the more detailed explanation.

That makes a lot more sense to me.

Unfortunately, i'm running into a few problems.

I need clarification whether

x = a3*t^3+a2*t^2+a1*t+a0

y = b3*t^3+b2*t^2+b1*t+b0

returns a point on a Catmull Rom Spline.

I'm currently using the equation below to return x,y points on the Catmull Rom spline for a given T value.

x = 0.5 * ((2 * p1x) + (p2x - p0x) * t + (2 * p0x - 5 * p1x + 4 * p2x - p3x) * t * t + (3 * p1x + p3x - p0x - 3 * p2x) * t * t * t)

y = 0.5 * ((2 * p1y) + (p2y - p0y) * t + (2 * p0y - 5 * p1y + 4 * p2y - p3y) * t * t + (3 * p1y + p3y - p0y - 3 * p2y) * t * t * t)

After coding your more detailed explanation i found i was getting wild results when trying to obtain the root.

Further exploration seems to reveal different spline formula's.

To verify this, i simply replaced my above spline formula's with yours to see if the same curve was drawn.

And the results were not the same.

Although both of our spline equations require 4 control points each, my spline gets drawn from point 1 to 2 ( with no drawing between points 0 and 3), and your spline was being drawn from point 0.

Which obviously leads me to ask if we are dealing with 2 different sets of spline equations ?

###
#8
Members - Reputation: **401**

Posted 25 August 2012 - 03:32 PM

Ok, i just assumed they were the control points for the spline.

In that case i'm at a loss as to where to go from here.

Can you elaborate on what the a0,a1,a2,a3,b0,b1,b2,b3 terms are relating to my specific problem, as i have no idea.

Thanks in advance

In that case i'm at a loss as to where to go from here.

Can you elaborate on what the a0,a1,a2,a3,b0,b1,b2,b3 terms are relating to my specific problem, as i have no idea.

Thanks in advance

###
#9
Members - Reputation: **401**

Posted 26 August 2012 - 02:47 AM

Heres what i have so far.

The code might serve to help others.

But it's not working correctly and the root returned isn't correct.

The terms a0,a1,a2,a3, b0,b1,b2,b3 are what i've found from google searches, i think they are correct for a Catmull Rom spline, but i'm not sure.

\\ Spline Control Points

\\ p0x#

\\ p0y#

\\ p1x#

\\ p1y#

\\ p2x#

\\ p2y#

\\ p3x#

\\ p3y#

t# = 0.5

spline_x_pos# = 0.5 * ((2 * p1x#) + (p2x# - p0x#) * t# + (2 * p0x# - 5 * p1x# + 4 * p2x# - p3x#) * t# * t# + (3 * p1x# + p3x# - p0x# - 3 * p2x#) * t# * t# * t#)

spline_y_pos# = 0.5 * ((2 * p1y#) + (p2y# - p0y#) * t# + (2 * p0y# - 5 * p1y# + 4 * p2y# - p3y#) * t# * t# + (3 * p1y# + p3y# - p0y# - 3 * p2y#) * t# * t# * t#)

A# = (line_point_1_y# - line_point_2_y#) * spline_x_pos#

B# = (line_point_2_x# - line_point_1_x#) * spline_y_pos#

a0# = -0.5 * p0x# + 1.5 * p1x# - 1.5 * p2x# + 0.5 * p3x#

a1# = p0x# - 2.5 * p1x# + 2 * p2x# - 0.5 * p3x#

a2# = -0.5 * p0x# + 0.5 * p2x#

a3# = p1x#

b0# = -0.5 * p0y# + 1.5 * p1y# - 1.5 * p2y# + 0.5 * p3y#

b1# = p0y# - 2.5 * p1y# + 2 * p2y# - 0.5 * p3y#

b2# = -0.5 * p0y# + 0.5 * p2y#

b3# = p1y#

c0# = (A# * a0#) + (B# * b0#)

c1# = (A# * a1#) + (B# * b1#)

c2# = (A# * a2#) + (B# * b2#)

c3# = (A# * a3#) + (B# * b3#)

root# = ((c3# * t#) ^ 3) + ((c2# * t#) ^ 2) + (c1# * t#) + c0#

And heres my previous code which does return a correct root.

t# = 0.5

spline_x_pos# = 0.5 * ((2 * p1x#) + (p2x# - p0x#) * t# + (2 * p0x# - 5 * p1x# + 4 * p2x# - p3x#) * t# * t# + (3 * p1x# + p3x# - p0x# - 3 * p2x#) * t# * t# * t#)

spline_y_pos# = 0.5 * ((2 * p1y#) + (p2y# - p0y#) * t# + (2 * p0y# - 5 * p1y# + 4 * p2y# - p3y#) * t# * t# + (3 * p1y# + p3y# - p0y# - 3 * p2y#) * t# * t# * t#)

temp_1# = (line_point_1_y# - line_point_2_y#) * spline_x_pos#

temp_2# = (line_point_2_x# - line_point_1_x#) * spline_y_pos#

temp_3# = ((line_point_1_x# * line_point_2_y#) - (line_point_2_x# * line_point_1_y#))

root# = temp_1# + temp_2# + temp_3#

The code might serve to help others.

But it's not working correctly and the root returned isn't correct.

The terms a0,a1,a2,a3, b0,b1,b2,b3 are what i've found from google searches, i think they are correct for a Catmull Rom spline, but i'm not sure.

\\ Spline Control Points

\\ p0x#

\\ p0y#

\\ p1x#

\\ p1y#

\\ p2x#

\\ p2y#

\\ p3x#

\\ p3y#

t# = 0.5

spline_x_pos# = 0.5 * ((2 * p1x#) + (p2x# - p0x#) * t# + (2 * p0x# - 5 * p1x# + 4 * p2x# - p3x#) * t# * t# + (3 * p1x# + p3x# - p0x# - 3 * p2x#) * t# * t# * t#)

spline_y_pos# = 0.5 * ((2 * p1y#) + (p2y# - p0y#) * t# + (2 * p0y# - 5 * p1y# + 4 * p2y# - p3y#) * t# * t# + (3 * p1y# + p3y# - p0y# - 3 * p2y#) * t# * t# * t#)

A# = (line_point_1_y# - line_point_2_y#) * spline_x_pos#

B# = (line_point_2_x# - line_point_1_x#) * spline_y_pos#

a0# = -0.5 * p0x# + 1.5 * p1x# - 1.5 * p2x# + 0.5 * p3x#

a1# = p0x# - 2.5 * p1x# + 2 * p2x# - 0.5 * p3x#

a2# = -0.5 * p0x# + 0.5 * p2x#

a3# = p1x#

b0# = -0.5 * p0y# + 1.5 * p1y# - 1.5 * p2y# + 0.5 * p3y#

b1# = p0y# - 2.5 * p1y# + 2 * p2y# - 0.5 * p3y#

b2# = -0.5 * p0y# + 0.5 * p2y#

b3# = p1y#

c0# = (A# * a0#) + (B# * b0#)

c1# = (A# * a1#) + (B# * b1#)

c2# = (A# * a2#) + (B# * b2#)

c3# = (A# * a3#) + (B# * b3#)

root# = ((c3# * t#) ^ 3) + ((c2# * t#) ^ 2) + (c1# * t#) + c0#

And heres my previous code which does return a correct root.

t# = 0.5

spline_x_pos# = 0.5 * ((2 * p1x#) + (p2x# - p0x#) * t# + (2 * p0x# - 5 * p1x# + 4 * p2x# - p3x#) * t# * t# + (3 * p1x# + p3x# - p0x# - 3 * p2x#) * t# * t# * t#)

spline_y_pos# = 0.5 * ((2 * p1y#) + (p2y# - p0y#) * t# + (2 * p0y# - 5 * p1y# + 4 * p2y# - p3y#) * t# * t# + (3 * p1y# + p3y# - p0y# - 3 * p2y#) * t# * t# * t#)

temp_1# = (line_point_1_y# - line_point_2_y#) * spline_x_pos#

temp_2# = (line_point_2_x# - line_point_1_x#) * spline_y_pos#

temp_3# = ((line_point_1_x# * line_point_2_y#) - (line_point_2_x# * line_point_1_y#))

root# = temp_1# + temp_2# + temp_3#

###
#10
Members - Reputation: **401**

Posted 26 August 2012 - 05:34 PM

For anyone thats interested i've created a small demo with source code demonstrating the intersection method described in my third post here.

655k. Windows only.

http://www.trinosis.com/files/line_to_catmul_rom_spline_intersection_demo.zip

655k. Windows only.

http://www.trinosis.com/files/line_to_catmul_rom_spline_intersection_demo.zip

###
#11
Members - Reputation: **136**

Posted 27 August 2012 - 12:59 PM

Thanks K_J_M.

It seems to me that basic, building-block code like this should be available in a standard form we could all share.

Instead it seems we all implement the same algorithms and (never having common test data in most cases) we often don't know about the bugs we have in our version.

The alternative is to buy some API that has more than what you need.

It seems to me that basic, building-block code like this should be available in a standard form we could all share.

Instead it seems we all implement the same algorithms and (never having common test data in most cases) we often don't know about the bugs we have in our version.

The alternative is to buy some API that has more than what you need.

###
#12
Members - Reputation: **401**

Posted 27 August 2012 - 05:07 PM

Hi Catmull Dog.

Yes, i totally agree with you.

Theres nothing worse than having to re invent the wheel.

I'm really suprised i havn't been able to find an already established algorithm for solving this problem.

But if i find one, i'll definately share it with the community.

Yes, i totally agree with you.

Theres nothing worse than having to re invent the wheel.

I'm really suprised i havn't been able to find an already established algorithm for solving this problem.

But if i find one, i'll definately share it with the community.