fast accurate terrain linesegment intersection

Started by
17 comments, last by zedz 16 years ago
the terrain is an evenly spaced triangle grid (pretty standard) ignoring the scenegraph phase for a sec. now whats the best method to check for a intersection? y == up im thinking of taking the XZ positions of the start + end points of the linesegment + creating a 2d line from this + check the triangles along this line, now if i use bresnhams line algorithm it will miss some quads, eg 1,1 hit, 2,1 hit but 2,2 the line goes through it but its not registered whats need is some sort of antialiased line. whats the best method? Xiaolin Wu's line algorithm or DDA or something else? cheers zed
Advertisement
Is this really a 2D problem? Or are you projecting the ray and terrain into 2D so you can do a DDA algorithm? If it is truly a 3D situation, I'd recommend a quadtree where you store each cell's min/max heights, and so can actually step "over" entire cells by doing a simple ray/box calc (rejecting many tris at once). If you get rid of the the 3rd dimension, you'll have to literally walk every triangle from the start of the projected ray to the end. Granted, you'll be walking along the ray, but that could still be a number of tris. I have tried the way you're talking about (a bresenham-style cell-walking alg), it wasn't amazingly elegant and I don't have the source code in front of me, but I recall just checking for whether my current x or y value had changed and adding in the skipped quad for consideration. It worked okay, but I believe the quadtree method is faster.
yes see the second line of my post :)
'ignoring the scenegraph phase for a sec.'
this concerns the testing of the tris in the final leaf eg say the final 16x16x2 triangles

the bresenham line will miss triangles, does the DDA method etc also miss them?

Quote:but I recall just checking for whether my current x or y value had changed and adding in the skipped quad for consideration
so u werent doing a true bresenham method, hmm i might have a look at this closer.
ta
Quote:yes see the second line of my post :)
'ignoring the scenegraph phase for a sec.'
this concerns the testing of the tris in the final leaf eg say the final 16x16x2 triangles
Apologies, wasn't sure if that's what you meant :).
Quote:so u werent doing a true bresenham method, hmm i might have a look at this closer.
No, it was not a true bresenham, but I guess I lump all those algorithms together (DDA, others). DDA will also miss those quads I believe. The part about checking for a change in x or y was where I diverged from those algorithms (i.e., if you stepped in both x and y, then you know you've possibly missed a quad by moving along a diagonal; I just handled this case manually within the line drawing algorithm, checking for which direction I just stepped, calculating the potentially missed quad; this is the part that might have a more elegant/natural solution, but it seemed to work fine and wasn't the most performance sensitive part of the loop).
Here's a variation of bresenham that seems like it would do what you want.

First, as you know, bresenham steps 1 unit at a time over the "long axis", and then updates the "short axis" by a fractional delta value. In your diagram, x is the long axis, y is the short. I'll just refer to them as that for the explanation, and the delta as dy/dx.

So my idea is to first define the x coordinate as the left edge of the "pixel", and the y coordinate as the top edge. So when iterating along the line, x is directly on the left edge (has no fractional portion), and y is the intersection of the line with that edge. Draw the pixel using that left edge coordinate (call this point curX,curY). Then step y (before stepping x), and if the integer portion of it changes, plot another pixel at (curX,nextY). Then step x and repeat.

That way you get a more "solid" line, where when the y coordinate moves down, the line steps down a pixel before stepping over.

One more thing to be aware of incase you haven't ever written a totally acurate line algorithm before: In order to have the left coordinate refer to the left edge of the pixel, you may need to "back up" along the line until you get to that edge, before you start iterating.
In your diagram, the line starts near the center of the pixel, so you'll need to step about -0.5 in the x direction, and -0.5*dy/dx in the y direction. That will get you to the intersection point of the line at the left edge.


...and during all that typing emeyex said basically the same thing. But hopefully this will help you with the specifics.
That's what I use.

#include <stdio.h>#include <stdlib.h>#include <math.h>#include <gl/glut.h>inline float frand(float x=1.0f){	return (rand() / (float) RAND_MAX) * x;}//===========================================================================// VECTORS//===========================================================================class Vector{public:	float x,y;public:	inline Vector(void)	{}	inline Vector(float Ix,float Iy)	: x(Ix)	, y(Iy)	{}	inline Vector &operator /=(const float Scalar)	{ x /= Scalar; y /= Scalar;		return *this; }	inline Vector &operator *=(const float Scalar)	{ x *= Scalar; y *= Scalar;		return *this; }		inline Vector &operator +=(const Vector &Other) { x += Other.x;	y += Other.y;	return *this; }	inline Vector &operator -=(const Vector &Other) { x -= Other.x;	y -= Other.y;	return *this;	}	inline float operator ^ (const Vector &V)	const	{	return (x * V.y) - (y * V.x); } // cross product	inline float operator * (const Vector &V)	const	{	return (x*V.x) + (y*V.y); } // dot product	inline Vector operator * (float  s)			const	{	return Vector(x*s, y*s); }		inline Vector operator / (float  s)			const	{	return Vector(x/s, y/s); }		inline Vector operator + (const Vector &V)	const	{	return Vector(x+V.x, y+V.y); }			inline Vector operator - (const Vector &V)	const	{	return Vector(x-V.x, y-V.y); }	friend Vector operator * (float k, const Vector& V) {	return Vector(V.x*k, V.y*k); }		inline Vector operator -(void) const { return Vector(-x, -y); }		inline float length(void) const { return (float) sqrt(x*x + y*y); }	void randomise(const Vector& min, const Vector& max)	{		x = frand(max.x - min.x) + min.x;		y = frand(max.y - min.y) + min.y;	}};//--------------------------------------------------------------------------// window size//--------------------------------------------------------------------------float width  = 640;float height = 480;float zoom = 30.0f;Vector start;Vector end;void Reset(){	start.randomise(Vector(0, 0), Vector(width / zoom, height / zoom));	end.randomise(Vector(0, 0), Vector(width / zoom, height / zoom));}void renderPixel(int x, int y){	glColor4f(1, 1, 1, 1);	glBegin(GL_LINE_LOOP);	glVertex2f(x, y);	glVertex2f(x+1, y);	glVertex2f(x+1, y+1);	glVertex2f(x, y+1);	glEnd();}void dda(const Vector& start, const Vector& end){	const float tolerance = 1.0E-8f;		// render segment	glColor4f(1, 1, 1, 1);	glBegin(GL_LINES);	glVertex2f(start.x, start.y);	glVertex2f(end.x, end.y);	glEnd();	// increment	Vector delta = end - start;	Vector inc;	inc.x = (fabs(delta.x) < tolerance)? 1.0f / tolerance : 1.0f / fabs(delta.x);	inc.y = (fabs(delta.y) < tolerance)? 1.0f / tolerance : 1.0f / fabs(delta.y);	// pixel coords	int x = (int)(start.x);	int y = (int)(start.y);	int dx = (delta.x < 0.0f)? -1 : (delta.x > 0.0f)? 1 : 0;	int dy = (delta.y < 0.0f)? -1 : (delta.y > 0.0f)? 1 : 0;	Vector accum;	accum.x = (delta.x < 0.0f)? (start.x - (x)) * inc.x : ((x+1) - start.x) * inc.x;	accum.y = (delta.y < 0.0f)? (start.y - (y)) * inc.y : ((y+1) - start.y) * inc.y;	float t = 0.0f;	// dda	while (t <= 1.0f)	{		renderPixel(x, y);		if(accum.x < accum.y)		{			t		 = accum.x;			accum.x += inc.x;			x		+= dx;		}		else		{			t		 = accum.y;			accum.y += inc.y;			y		+= dy;		}	}}//-----------------------------------------------------// displays the objects//-----------------------------------------------------void Display(){	//--------------------------------------------------------------------------	// render stuff	//--------------------------------------------------------------------------	glClearColor(0, 0, 0, 0);	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);	glMatrixMode(GL_PROJECTION);	glLoadIdentity();	glOrtho(0, width / zoom, 0, height / zoom, -100, 100);		//-----------------------------------------------------------------	// Setup the model view matrix	//-----------------------------------------------------------------	glMatrixMode(GL_MODELVIEW);	glLoadIdentity();		dda(start, end);	glutSwapBuffers();}void Timer(int t){	Display();	glutTimerFunc(t, Timer, (int) 500.0f / 60.0f);}	void Reshape(int w, int h){	width  = w;	height = h;	glViewport(	0, 0, w, h);}void Keyboard(unsigned char key, int x, int y){	if (key == 27)		exit(0);	if(key == ' ')		Reset();}void main(int argc, char** argv){	//--------------------------------------------------------------------------	// OpenGL / GLUT init	//--------------------------------------------------------------------------    glutInit( &argc, argv );	glutInitDisplayMode		(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);		glutInitWindowSize		(width, height);	glutInitWindowPosition	(0, 0);	glutCreateWindow		("dda");		glEnable				(GL_BLEND);	glBlendFunc				(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);	glDisable				(GL_DEPTH_TEST);	glDisable				(GL_LIGHTING);		glutDisplayFunc			(Display);	glutReshapeFunc			(Reshape);	glutTimerFunc			(0, Timer, (int) 100.0f / 60.0f);	glutKeyboardFunc		(Keyboard);	Reset					();	glutMainLoop			();}

Everything is better with Metal.

3D version (I was bored)

#include <stdio.h>#include <stdlib.h>#include <math.h>#include <gl/glut.h>inline float frand(float x=1.0f){	return (rand() / (float) RAND_MAX) * x;}//===========================================================================// VECTORS//===========================================================================template<typename TYPE>class Vector3{public:		union	{		struct		{			TYPE x,y,z;		};		TYPE coords[3];	};public:	inline Vector3(void)	{}	inline Vector3(TYPE Ix,TYPE Iy, TYPE Iz)	: x(Ix)	, y(Iy), z(Iz)	{}	inline TYPE  operator[] (const int i) const { return coords; }	inline TYPE& operator[] (const int i)       { return coords; }		inline Vector3 &operator /=(const TYPE k)		{ x /= k; y /= k; z /= k;		return *this; }	inline Vector3 &operator *=(const TYPE k)		{ x *= k; y *= k; z /= k; 		return *this; }	inline Vector3 &operator +=(const Vector3 &v)	{ x += v.x;	y += v.y; z += v.z; return *this; }	inline Vector3 &operator -=(const Vector3 &v)	{ x -= v.x;	y -= v.y; z -= v.z;	return *this; }	inline Vector3 operator * (TYPE  s)				const	{	return Vector3(x*s, y*s, z*s); }	inline Vector3 operator / (TYPE  s)				const	{	return Vector3(x/s, y/s, z/s); }	inline Vector3 operator + (const Vector3 &v)	const	{	return Vector3(x+v.x, y+v.y, z+v.z); }	inline Vector3 operator - (const Vector3 &v)	const	{	return Vector3(x-v.x, y-v.y, z-v.z); }	friend Vector3 operator * (TYPE k, const Vector3& v)	{	return Vector3(v.x*k, v.y*k, v.z*k); }	inline Vector3 operator ^ (const Vector3 &v)	const	{ return Vector3((y*v.z)-(z*v.y), (z*v.x)-(x*v.z), (x*v.y)-(y*v.x)); } // cross product	inline TYPE    operator * (const Vector3 &v)	const	{ return (x*v.x) + (y*v.y) + (z*v.z); } // dot product	inline Vector3 operator - (void)				const	{ return Vector3(-x, -y, -z); }	inline TYPE length(void) const { return (TYPE) sqrt(x*x + y*y + z*z); }};typedef Vector3<float> Vector;typedef Vector3<int>   Voxel;Vector randomise(const Vector& min, const Vector& max){	Vector v;	v.x = (frand(max.x - min.x) + min.x);	v.y = (frand(max.y - min.y) + min.y);	v.z = (frand(max.z - min.z) + min.z);	return v;}//--------------------------------------------------------------------------// window size//--------------------------------------------------------------------------float width    = 640;float height   = 480;float size	   = 30.0f;float camdist  = size * 0.5f;float camtheta = 0.0f;float camphi   = 0.0f;Vector start;Vector end;void Reset(){	start = randomise(Vector(0, 0, 0), Vector(size, size, size));	end   = randomise(Vector(0, 0, 0), Vector(size, size, size));}void renderSegment(const Vector& start, const Vector& end){	// render segment	glLineWidth(3.0f);	glColor4f(1, 0.5f, 0.5f, 1);	glBegin(GL_LINES);	glVertex3f(start.x, start.y, start.z);	glVertex3f(end.x, end.y, end.z);	glEnd();}void renderVoxel(const Voxel& voxelPos){	glPushMatrix();	glLineWidth(1.0f);	glTranslatef(voxelPos.x + 0.5f, voxelPos.y + 0.5f, voxelPos.z + 0.5f);		glColor4f(1, 1, 1, 1);	glutWireCube(1.0f);	glColor4f(0.2f, 1.0f, 0.5f, 0.1f);	glutSolidCube(1.0f);		glPopMatrix();}void renderDDA(const Vector& start, const Vector& end){	const float tolerance = 1.0E-8f;		// dda increment	Vector delta = end - start;	Vector inc;	inc.x = (fabs(delta.x) < tolerance)? 1.0f / tolerance : 1.0f / fabs(delta.x);	inc.y = (fabs(delta.y) < tolerance)? 1.0f / tolerance : 1.0f / fabs(delta.y);	inc.z = (fabs(delta.z) < tolerance)? 1.0f / tolerance : 1.0f / fabs(delta.z);	// voxel coords	Voxel  voxelpos((int)(start.x), (int)(start.y), (int)(start.z));	Voxel  voxelinc((delta.x < 0.0f)? -1 : (delta.x > 0.0f)? 1 : 0,					(delta.y < 0.0f)? -1 : (delta.y > 0.0f)? 1 : 0,					(delta.z < 0.0f)? -1 : (delta.z > 0.0f)? 1 : 0);	// dda accumulator	Vector accum;	accum.x = (delta.x < 0.0f)? (start.x - (voxelpos.x)) * inc.x : ((voxelpos.x+1) - start.x) * inc.x;	accum.y = (delta.y < 0.0f)? (start.y - (voxelpos.y)) * inc.y : ((voxelpos.y+1) - start.y) * inc.y;	accum.z = (delta.z < 0.0f)? (start.z - (voxelpos.z)) * inc.z : ((voxelpos.z+1) - start.z) * inc.z;	float t = 0.0f; // terminating condition	// dda	while (t <= 1.0f)	{		// callback when the ray hits a new voxel.		renderVoxel(voxelpos);		// find axis with minimum accumulated value		int min = 0;		if(accum[1] < accum[min]) min = 1;		if(accum[2] < accum[min]) min = 2;				// add the increments		t				 =    accum[min];		accum[min]		+=      inc[min];		voxelpos[min]   += voxelinc[min];	}}//-----------------------------------------------------// displays the objects//-----------------------------------------------------void Display(){	//--------------------------------------------------------------------------	// render stuff	//--------------------------------------------------------------------------	glClearColor(0.1, 0.1, 0.1, 1);	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);	glMatrixMode(GL_PROJECTION);	glLoadIdentity();	gluPerspective(80, (float)width/(float)(height+1), 0.1f, 1000.0f);		//-----------------------------------------------------------------	// Setup the model view matrix	//-----------------------------------------------------------------	glMatrixMode(GL_MODELVIEW);	glLoadIdentity();		Vector camDir	= Vector(	cos(camphi) * cos(camtheta), 								sin(camphi), 								cos(camphi) * sin(camtheta));	Vector camRight = Vector(	-sin(camtheta), 								0.0f,								cos(camtheta));	Vector camUp	= (camDir ^ camRight);	Vector camTarget = (start + end) * 0.5f;	Vector camPos	= camTarget + camDir * -camdist;	gluLookAt(camPos.x, camPos.y, camPos.z, camTarget.x, camTarget.y, camTarget.z, camUp.x, camUp.y, camUp.z);	// show stuff	renderSegment(start, end);	renderDDA(start, end);	glutSwapBuffers();}void Timer(int t){	Display();	glutTimerFunc(t, Timer, (int) 500.0f / 60.0f);}	void Reshape(int w, int h){	width  = w;	height = h;	glViewport(	0, 0, w, h);}void ExtKeyboard(int key, int x, int y){	if(key == GLUT_KEY_DOWN)		camdist += size / 100.0f;	if(key == GLUT_KEY_UP)		camdist -= size / 100.0f;}void Keyboard(unsigned char key, int x, int y){	if (key == 27)		exit(0);	if(key == ' ')		Reset();}void Mouse(int x, int y){	static int mousex = -1;	static int mousey = -1;	if(mousex == -1 || mousey == -1)	{		mousex = x;		mousey = y;	}	else	{		const float speed = 0.005f;		int deltax = x - mousex;		int deltay = y - mousey;		mousex = x;		mousey = y;		camtheta += deltax * speed;		camphi   -= deltay * speed;	}}void main(int argc, char** argv){	//--------------------------------------------------------------------------	// OpenGL / GLUT init	//--------------------------------------------------------------------------    glutInit( &argc, argv );	glutInitDisplayMode		(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);		glutInitWindowSize		(width, height);	glutInitWindowPosition	(0, 0);	glutCreateWindow		("dda 3D");		glEnable				(GL_BLEND);	glBlendFunc				(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);	glDisable				(GL_DEPTH_TEST);	glDisable				(GL_LIGHTING);		glutDisplayFunc			(Display);	glutReshapeFunc			(Reshape);	glutTimerFunc			(0, Timer, (int) 100.0f / 60.0f);	glutSpecialFunc			(ExtKeyboard);	glutKeyboardFunc		(Keyboard);	glutPassiveMotionFunc	(Mouse);	Reset					();	glutMainLoop			();}

Everything is better with Metal.

thanks very much for that guys

i see how that fits in with what im doing. actually im thinking, perhaps ignoring the scenegraph is quickest just cast a line over the whole terrain, in theory it is
I dont really see how you could do a faster trace. At the end of the day, it's only a few additions and comparaisons, and it will scan every single cell along the path. However, if you subdivided your heightmap with a quad tree and store the height value of each quadtree cell, there is potential for fast rejects if for example, the ray fires high above the terrain.

Note that for a ray/segment, it's usually useful to cache the inverse of the direction vector (or the delta in my case), as it's costly to compute, and it's used in every ray interection routines.

Everything is better with Metal.

Are you doing this for generating shadow maps? If so this might be a good read.

http://www.gamedev.net/reference/articles/article1817.asp

This topic is closed to new replies.

Advertisement