Line segment intersecting curve

Started by
37 comments, last by apatriarca 16 years, 4 months ago
The use of #defines to declare constants was a common practice in C, but it’s discouraged in C++ (and also in some standard for critical code in C).

The best way to declare a constant in C++ is using a const variable. The advantages are:
1. It’s local to the block
2. It’s type-safe
3. It’s easier to debug
4. There isn’t any risk to introduce errors very hard to find
...

Another acceptable way to declare constants is using enums but the constants must have integer types and the type is chosen by the compiler. So in C++ is better to use the const modifier.

Why you verifies if the roots are equals. You only add a case that you have to verifies in the code that use the function. It would be more useful to order the roots based on the x (I don’t know if your code gives the roots in order).
In the last part of your code it seem that you calculate the complex conjugate solutions. But you are storing that values in a single double (for a complex number you need 2). There is no reason to do that calculation.
I haven’t test the code. It gives the right intersections?
Advertisement
Well, new SolveCubic implementation correct solve cubical equations...
About intersections I'm not sure, because now we need 4 coefficients instead of 3 as it was before.

I try to fill four coefficient, but it isn't correct :

bool IsLineSegmentIntersectCurve( BezierCurvePoints bCurve, Vector vSrc, Vector vEnd, Vector &M, Vector &N, float &t ){        Vector tangent; // TODO: Implement tangent equation and compute normal later!        // correct part of code :	Vector N = (vEnd-vSrc).Normal();	float c = vSrc.y*N.x - vSrc.x*N.y;	float D0 = dist(N, c, bCurve.P0);	float D1 = dist(N, c, bCurve.P1);	float D2 = dist(N, c, bCurve.P2);	float D3 = dist(N, c, bCurve.P3);        // incorrect part of code : 	float a = 3.0f*(D2 - 2.0f*D1 + D0)/(D3 - 3.0f*D2 + 3.0f*D1 - D0);	float b = (D3 - 3.0f*D2 + 3.0f*D1 - D0);	float c = D0/(D3 - 3.0f*D2 + 3.0f*D1 - D0);	float d = 3.0f*(D1 - D0);        // ax^3 + bx^2 + cx + d = 0	double results[3]; // array with taken results	int roots = SolveCubic( a,b,c,d, results );                ... // Root's check code still same    }


So with this four coefficients (wich can be incorrect), I take intersections, and some time it corret, but some time interection occur when it should'n, or I take intersection point on curve but they not in segment range.

[Edited by - Lolmen on December 7, 2007 12:31:47 PM]
The correct coefficients are:
a = (D3 - 3D2 + 3D1 - D0)
b = 3(D2 - 2D1 + D0)
c = 3(D1 - D0)
d = D0
Damn! Now it absolutely correct working for vertical lines!!!
But if line is horisontal or non vertical and non horisontal, so then no intersections was returns...

Can you give me analitical conclusion why this happen?
Is that because d coefficient equal zero and we should check that inside SolveCubic( ... )?

Also, can you give me you ICQ or say me time when I can cath you for couple of hours???
Because wait for 5-10 hours to meet you again is too long :) I need to talk closer, There is few thing I have to understand, but I need to talk about in real time...
I’m sorry but I haven’t ICQ.
I’ll send to you a PM with the IM contacts that I have.

You verify the roots in order (before results[0], then results[1] and finally results[2]). But your solutions aren’t ordered. What you want is the solution nearer to A. So you have to test the distance between every intersection point inside the bézier curve and the segment and A.

You have to test if your SolveCubic works correctly. You can substitute your function with a function that you know works correctly and see if the intersections are the correct ones (try to show all the intersections between the segment and the Bézier curve). If your function will have the same behaviors that the problem it’s elsewhere.

Well, Seems like I'm check code as you say... Because It works for vertical lines, I debug it a little and see that for t < 0.5 SolveCubic works correctly on vertical lines, and for t > 0.5 SolveCubic call SolveQuadratic, why other cases not work I have no clue.

	roots = SolveCubic( a,b,c,d, results );		// local definition	#define IS_BAD_ROOT( x )	( (x) < 0.0f || (x) > 1.0f )	if ( roots < 0 ){ return false; } // no real roots out there	float u=0;	bool badRoot1 = IS_BAD_ROOT(results[0]);	bool badRoot2 = IS_BAD_ROOT(results[1]);	bool badRoot3 = IS_BAD_ROOT(results[2]);		if ( roots == 3 && badRoot1 && badRoot2 && badRoot3 )		return false;	if ( roots == 2 && badRoot1 && badRoot2 )		return false;	if ( roots == 1 && badRoot1 )		return false;	if ( roots < 1 )		return false;	if (!badRoot1 && roots == 1 )	{		// Bezier: (1-t)^3 * P0 + 3t(1-t)^2 * P1 + 3t^2 * (1-t)P2 + t^3 * P3, t -> [0, 1]		M = B( results[0], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );		if ( IsLineSegmentIntersectPoint(vSrc, vEnd, M, u) )			return true;		return false;	}	if ( roots == 2 )	{		if ( !badRoot1 && badRoot2 )		{			M = B( results[0], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			if ( IsLineSegmentIntersectPoint(vSrc, vEnd, M, u) )				return true;			return false;		}		if ( badRoot1 && !badRoot2 )		{			// Bezier: (1-t)^3 * P0 + 3t(1-t)^2 * P1 + 3t^2 * (1-t)P2 + t^3 * P3, t -> [0, 1]			M = B( results[1], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			if ( IsLineSegmentIntersectPoint(vSrc, vEnd, M, u) )				return true;			return false;		}		if ( !badRoot1 && !badRoot2 )		{			Vector M1 = B( results[0], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			Vector M2 = B( results[1], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			int root = 0; // This variable for further implementation			Vector d1 = vSrc - M1;			Vector d2 = vSrc - M2;			if ( d1.Length() >= d2.Length() ){ M = M2; root = 1; }			else { M = M1; }			if ( IsLineSegmentIntersectPoint(vSrc, vEnd, M, u) )				return true;			return false;		}	}	if ( roots == 3 )	{		// root 1 is good		if ( !badRoot1 && badRoot2 && badRoot3 )		{			M = B( results[0], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			if ( IsLineSegmentIntersectPoint(vSrc, vEnd, M, u) )				return true;			return false;		}		// root 2 is good		if ( badRoot1 && !badRoot2 && badRoot3 )		{			M = B( results[1], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			if ( IsLineSegmentIntersectPoint(vSrc, vEnd, M, u) )				return true;			return false;		}		// root 3 is good		if ( badRoot1 && badRoot2 && !badRoot3 )		{			M = B( results[2], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			if ( IsLineSegmentIntersectPoint(vSrc, vEnd, M, u) )				return true;			return false;		}		// root 1 & 2 is good		if ( badRoot1 && badRoot2 && !badRoot3 )		{			Vector M1 = B( results[0], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			Vector M2 = B( results[1], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			int root = 0;			Vector d1 = vSrc - M1;			Vector d2 = vSrc - M2;			if ( d1.Length() >= d2.Length() ){ M = M2; root = 1; }			else { M = M1; }			if ( IsLineSegmentIntersectPoint(vSrc, vEnd, M, u) )				return true;			return false;		}		// root 2 & 3 is good		if ( badRoot1 && !badRoot2 && !badRoot3 )		{			// Bezier: (1-t)^3 * P0 + 3t(1-t)^2 * P1 + 3t^2 * (1-t)P2 + t^3 * P3, t -> [0, 1]			Vector M1 = B( results[1], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			Vector M2 = B( results[2], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			int root = 1;			Vector d1 = vSrc - M1;			Vector d2 = vSrc - M2;			if ( d1.Length() >= d2.Length() ){ M = M2; root = 2; }			else { M = M1; }			if ( IsLineSegmentIntersectPoint(vSrc, vEnd, M, u) )				return true;			return false;		}		// root 1 & 3 is good		if ( !badRoot1 && badRoot2 && !badRoot3 )		{			// Bezier: (1-t)^3 * P0 + 3t(1-t)^2 * P1 + 3t^2 * (1-t)P2 + t^3 * P3, t -> [0, 1]			Vector M1 = B( results[0], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			Vector M2 = B( results[2], bCurve.P0, bCurve.P1, bCurve.P2, bCurve.P3 );			int root = 0;			Vector d1 = vSrc - M1;			Vector d2 = vSrc - M2;			if ( d1.Length() >= d2.Length() ){ M = M2; root = 2; }			else { M = M1; }			if ( IsLineSegmentIntersectPoint(vSrc, vEnd, M, u) )				return true;			return false;		}	}
I haven't understanded what you mean. SolveCubic should work the same for every value of t.
Can you show me images of the intersections found by the code?

As I've said the normal is equals to the tangent vector rotated by 90°. Given the tangent vector (x,y) the normals are (-y,x) and (y,-x). You have to choose the one that is faced to A.
I have done some math (that I hope is correct) and I'm arrived to the following formula for the tangent vector:
T(t) = 3((P1 - P0)(1-t)2 + 2(P2 - P1)(1-t)t + (P3 - P2)t2)
I controll line segment by mouse in real time, can move end point by cursor or translate whole segment by press LMB.

For vertical lines all ok, not matter 1 root, 2 or 3 I take...
For horisontal everytime is -2 roots, wich mean it's complex.
For random oriented and faced line is noithing at all, it returns 1 or 3 roots, but this roots is incorrect...

So here you go :

Vertical line:

Verical line faced up

Vertical line other direction and position :

Vertical line faced down

Horisontal line example :

Horisontal line

Random oriented line example :

Random oriented line

So, if cubical solver and quadratic solver can solve equation for vertical line, so I'm assure that something goes wrong with next things :
( look at comments )

                // this ( line can have vEnd and vSrc at same point, so delta then will be 0,0 and it wreck all.	Vector delta = (vEnd-vSrc).Normal();  	float C = vSrc.y*delta.x - vSrc.x*delta.y;                 // this : Maybe we should add some checking and         // reconsider Dx depends  on case of line and curve points....	float D0 = dist(delta, C, bCurve.P0);	float D1 = dist(delta, C, bCurve.P1);	float D2 = dist(delta, C, bCurve.P2);	float D3 = dist(delta, C, bCurve.P3);        // this : if horisontal line throw two complex roots because        // coefficient "a" equal zero and Discriminant of quadratic equation        // is less than zero, then we should retransform this stuff,         // also maybe something wrong there and also maybe this code should        // be whole retransformed, and also maybe there sould be some        // checkings for specific cases and tricks like do some :        // swap( a, c ) swap( b, d ) because playing with chenge structure of        // this formulas can give us sometimes an correct results for random oriented lines...	double a = (D3 - 3.0*D2 + 3.0*D1 - D0);	double b = 3.0*(D2 - 2.0*D1 + D0);	double c = 3.0*(D1 - D0);	double d = D0;


Roots checking work correct, and equation solver as we see too...

[Edited by - Lolmen on December 9, 2007 2:11:41 PM]
The equation of the line is correct. If B-A is equal to 0 than it isn’t a segment but a point. In this case there are infinite lines that contains A=B.
This case must be handled in a different way.

If there are two complex roots that the cubic curve never intersect the line.

In your code I am unable to find where you set the parameter t that I think is used to return the intersection point.

This topic is closed to new replies.

Advertisement