🔒 Locked

union of two spheres

Started by w_poons Mar 24, 2006 at 5:30 AM 10 replies 2.3k views

Original Post

Mar 24, 2006 at 5:30 AM

hi, i trying to get the union of two polygonial spheres. i know that it could be solved by using CSG or metaballs and then construct the polygonal representation. but maybe there is a simpler way for this problem. the spheres have the same size and are all positioned at the same z-coordinate. does anybody know about a solution or can point me to some resources? thanks!

haegarr

Mar 24, 2006 at 5:54 AM

The easiest way is IMHO to construct the union's mesh from an implicit description (you know all needed data for an explicit description AFAIK). This would yield in a principal solution, of course, but not in an exact solution w.r.t. the before given original spheres.

Working with the polygonal data directly is more general, of course, but requires more effort, too.
(1) Create a new mesh.
(2) Iterate all polygons of sphere A.
(a) Copy the polygon into the new mesh if it is fully contained by sphere B.
(b) Ignore the polygon if it is fully outside of sphere B.
(c) Split the polygon and copy the contained part if it lies partly inside sphere B.
(3) Perform all of (2) by exchanging spheres A and B.
(4) Iterate the vertices of the new mesh and merge duplicate vertices.

Checking whether a polygon is inside/outside could be done by first bbox test, followed by plane tests if needed. Splitting could also be done by plane tests (see e.g. Sutherland-Hodgeman algorithm).

w_poons

Mar 24, 2006 at 6:49 AM

thanks for your suggestions. i think i will rather stick to the polygonal approach. but could it be that your solution is creating the intersection of the two spheres instead of the union?

do you think the polygonal approach is possible in a realtime application?

haegarr

Mar 24, 2006 at 7:08 AM

Quote:
Original post by w_poons
but could it be that your solution is creating the intersection of the two spheres instead of the union?

Ups, right you are. I've confused union w/ intersection, sorry.

Well, for the union algorithm one needs to save the polygons that are dropped in case of intersection, and vice-versa.

(However, are your spheres transparent? Often a union is simulated by simply letting the parts overlap and doing nothing special. But I assume you've already rejected that possibility ...)

Quote:
Original post by w_poons
do you think the polygonal approach is possible in a realtime application?

In your special case you have exactly 2 shapes only, both being convex. This a-priori knowledge should allow some optimizations. IMHO, if the resolution of the polygonal approximation isn't too high, it should be possible to be done in realtime on a halfway modern computer.

w_poons

Mar 24, 2006 at 8:00 AM

Quote:
Original post by haegarr
(However, are your spheres transparent? Often a union is simulated by simply letting the parts overlap and doing nothing special. But I assume you've already rejected that possibility ...)

yes they are and another issue is that a shader applies some wobbling effect to the vertices. so i'm not sure if this is visible if the two spheres "melt together".

Quote:
Original post by haegarr
In your special case you have exactly 2 shapes only, both being convex. This a-priori knowledge should allow some optimizations. IMHO, if the resolution of the polygonal approximation isn't too high, it should be possible to be done in realtime on a halfway modern computer.

there are more spheres in the scene. but only between pairs of them the union operation would be performed!

ApochPiQ

Mar 25, 2006 at 8:00 PM

I'm assuming that these are truly spheres and not deformed spheroids - i.e. they are definitely a surface of points a constant radius away from a center point. If this is the case... why not just re-tesselate the spheres on the fly as needed? You can do some simple cheats here:

- Tesselate both spheres separately and cache the polygons that result

- If two spheres intersect, retesselate them

- When tesselating one sphere, check if the current point is inside the other
sphere. If so, do not form a polygon at that point; instead, mark the other vertices of that polygon as "edges" and move on

- This step should result in the union, plus a jagged area in the overlapping region where the two spheres touch; in that region there will be no polygons

- Find the plane of intersection between the two spheres (i.e. the plane where the two surfaces intersect) and compute the radius of the circle that lies on that plane. This should be pretty trivial geometry.

- For each sphere, iterate over the marked edge points. When you find a pair of edge points, generate a point on the intersection circle that can form a triangle with the two edge points, and add that to the polygon list.

The final result should be an optimal tesselation of both spheres to within the precise degree of detail you require, with no requirement to split or merge polygons across the two spheres. By caching the tesselations of non-intersecting spheres, you can avoid any performance hit at all for the general case; you only need to retesselate if two intersecting spheres move relative to each other. Sphere tesselation is a pretty trivial process as well, and easy to optimize, so you should end up with good performance overall.

w_poons

Mar 27, 2006 at 2:55 AM

thanks for your answer!

Quote:
Original post by ApochPiQ
I'm assuming that these are truly spheres and not deformed spheroids - i.e. they are definitely a surface of points a constant radius away from a center point.

yes they are true spheres, but the vertices get deformed in the vertex shader to apply some wobbling effect to them.

Quote:
Original post by ApochPiQ
- Tesselate both spheres separately and cache the polygons that result

since tessellation is a new topic to me i did some google search on sphere tesselation and found out that several algorithms exist. is it worth to implement a more sophisticated algorithm for small spheres (say about 300 facets) or will it do a simple lattitude-longitude tessellation?

SKREAMZ

Mar 27, 2006 at 3:15 AM

i did this about 2 years back using csg, i know your want something simpler but i cant help u there except by posting my code

BSPCSG.h

struct TRI{	sklib::XVec3 p1;	sklib::XVec3 p2;	sklib::XVec3 p3;	void operator =(const TRI &tri)	{		p1 = tri.p1;		p2 = tri.p2;		p3 = tri.p3;	}	sklib::XVec3 &operator [](int i)	{		switch(i)		{		case 0:			return p1;		case 1:			return p2;		case 2:			return p3;		default:			return p3;		};	}};class BspCsg{public:	BspCsg();	~BspCsg();	void Build(std::vector<TRI> &polys);	int Clip(TRI tri, std::vector<TRI> &fpolys, std::vector<TRI> &bpolys);	int ClassifyTriangle(TRI tri);	int SplitTriangle(TRI tri, TRI &front, TRI &back, TRI &extra);	void FindNormal(TRI tri, sklib::XVec3 &normal);public:	BspCsg *front;	BspCsg *back;	sklib::XPlane plane;};

and BSPCSG.cpp

BspCsg::BspCsg(){	front = NULL;	back = NULL;}BspCsg::~BspCsg(){	if(front != NULL)	{		delete front;		front = NULL;	}	if(back != NULL)	{		delete back;		back = NULL;	}}void BspCsg::Build(std::vector<TRI> &polys){	std::vector<TRI> fpolys;	std::vector<TRI> bpolys;	for(int i = 0; i < (int)polys.size(); i++)	{		TRI tri = polys;		int re = ClassifyTriangle(tri);		if(re == FRONT)		{			if(front == NULL)			{				front = new BspCsg;				front->plane.FromPoints(tri.p1, tri.p2, tri.p3);			}			fpolys.push_back(tri);		}		else if(re == BACK)		{			if(back == NULL)			{				back = new BspCsg;				back->plane.FromPoints(tri.p1, tri.p2, tri.p3);			}			bpolys.push_back(tri);		}		else if(re == SPANNING)		{			if(front == NULL)			{				front = new BspCsg;				front->plane.FromPoints(tri.p1, tri.p2, tri.p3);			}			fpolys.push_back(tri);			if(back == NULL)			{				back = new BspCsg;				back->plane.FromPoints(tri.p1, tri.p2, tri.p3);			}			bpolys.push_back(tri);		}	}	if((int)fpolys.size() > 0)        		front->Build(fpolys);	if((int)bpolys.size() > 0)		back->Build(bpolys);}int BspCsg::Clip(TRI tri, std::vector<TRI> &fpolys, std::vector<TRI> &bpolys){	int re = ClassifyTriangle(tri);	if(re == FRONT)	{		if(front == NULL)		{			fpolys.push_back(tri);			return FRONT;		}		else			return front->Clip(tri, fpolys, bpolys);	}	else if(re == BACK)	{		if(back == NULL)		{			bpolys.push_back(tri);			return BACK;		}		else			return back->Clip(tri, fpolys, bpolys);	}	else if(re == SPANNING)	{		TRI pfront, pback, pextra;		int fback = -1;		int ffront = -1;		int fextra = -1;		int r = SplitTriangle(tri, pfront, pback, pextra);		if(front == NULL)		{			fpolys.push_back(pfront);			ffront = FRONT;		}		else			ffront = front->Clip(pfront, fpolys, bpolys);		if(back == NULL)		{			bpolys.push_back(pback);			fback = BACK;		}		else			fback = back->Clip(pback, fpolys, bpolys);		if(r == FRONT)		{			if(front == NULL)			{				fpolys.push_back(pextra);				fextra = FRONT;			}			else				fextra = front->Clip(pextra, fpolys, bpolys);		}		else if(r == BACK)		{			if(back == NULL)			{				bpolys.push_back(pextra);				fextra = BACK;			}			else				fextra = back->Clip(pextra, fpolys, bpolys);		}		if(ffront == FRONT && fback == FRONT && fextra == FRONT)		{			fpolys.pop_back();			fpolys.pop_back();			fpolys.pop_back();			fpolys.push_back(tri);			return FRONT;		}		if(ffront == BACK && fback == BACK && fextra == BACK)		{			bpolys.pop_back();			bpolys.pop_back();			bpolys.pop_back();			bpolys.push_back(tri);			return BACK;		}		return SPANNING;	}	else	{		sklib::XVec3 norm(plane.a, plane.b, plane.c);		sklib::XVec3 norm2;		FindNormal(tri, norm2);		float dot = norm.Dot(norm2);		if(dot > 0 + EPSILON) // facing in the same direction		{			if(front == NULL)			{				fpolys.push_back(tri);				return FRONT;			}			else				return front->Clip(tri, fpolys, bpolys);		}		else if(dot < 0 - EPSILON) //facing in oppisite directions		{			if(back == NULL)			{				bpolys.push_back(tri);				return BACK;			}			else				return back->Clip(tri, fpolys, bpolys);		}	}	return -1;}int BspCsg::ClassifyTriangle(TRI tri){	bool front = false;	bool back = false;	for(int i = 0; i < 3; i++)	{		float re = plane.Distance(tri);		if(re > 0 + EPSILON)			front = true;				if(re < 0 - EPSILON)			back = true;	}	if(front == true && back == false)		return FRONT;	if(front == false && back == true)		return BACK;	if(front == true && back == true)		return SPANNING;	return ONPLANE;}int BspCsg::SplitTriangle(TRI tri, TRI &front, TRI &back, TRI &extra){	sklib::XVec3 point;	sklib::XVec3 fpoints[4];	sklib::XVec3 bpoints[4];	int f = 0;	int b = 0;	for(int i = 0; i < 3; i++)	{		float re = plane.Distance(tri);		if(re > 0 + EPSILON) //FRONT		{			fpoints[f++] = tri;			if(plane.Distance(tri[i == 2 ? 0 : i + 1]) < 0 - EPSILON)			{				if(D3DXPlaneIntersectLine(point, plane, tri, tri[i == 2 ? 0 : i + 1]) != NULL)				{					fpoints[f++] = point;					bpoints[b++] = point;									}			}		}		else if(re < 0 - EPSILON) //BACK		{			bpoints[b++] = tri;			if(plane.Distance(tri[i == 2 ? 0 : i + 1]) > 0 + EPSILON)			{				if(D3DXPlaneIntersectLine(point, plane, tri, tri[i == 2 ? 0 : i + 1]) != NULL)				{					bpoints[b++] = point;					fpoints[f++] = point;				}			}		}		else		{			bpoints[b++] = tri;			fpoints[f++] = tri;		}	}	front.p1 = fpoints[0];	front.p2 = fpoints[1];	front.p3 = fpoints[2];	back.p1 = bpoints[0];	back.p2 = bpoints[1];	back.p3 = bpoints[2];	if(f == 4 && b == 3)	{		extra.p1 = fpoints[2];		extra.p2 = fpoints[3];		extra.p3 = fpoints[0];		return FRONT; //has been split into 3 triangles (2 front & 1 back)	}	if(f == 3 && b == 4)	{		extra.p1 = bpoints[2];		extra.p2 = bpoints[3];		extra.p3 = bpoints[0];		return BACK; //has been split into 3 triangles (1 front & 1 back)	}	return SPANNING; //has been split into 2 triangles (1 front & 1 back)}void BspCsg::FindNormal(TRI tri, sklib::XVec3 &normal){	sklib::XVec3 vert = tri.p3 - tri.p2;	normal = tri.p2 - tri.p1;	normal = normal.Cross(vert);	normal.Normalize();}

and code to use the class

void CMain::CSGSubtract(){	if(rootnode == NULL)	{		rootnode = new BspCsg;		TRI tri = brush1.front();		rootnode->plane.FromPoints(tri.p1, tri.p2, tri.p3);		rootnode->Build(brush1);		std::vector<TRI> fpolys;		std::vector<TRI> bpolys;		for(int i = 0; i < (int)brush2.size(); i++)            rootnode->Clip(brush2, fpolys, bpolys);		std::vector<TRI> tlist = fpolys;		delete rootnode;		rootnode = NULL;		rootnode = new BspCsg;		tri = brush2.front();		rootnode->plane.FromPoints(tri.p1, tri.p2, tri.p3);		rootnode->Build(brush2);		fpolys.clear();		bpolys.clear();		for(int i = 0; i < (int)brush1.size(); i++)            rootnode->Clip(brush1, fpolys, bpolys);		delete rootnode;		rootnode = NULL;		brush2.clear();		brush2 = tlist;		for(int i = 0; i < (int)bpolys.size(); i++)		{			TRI ftri = bpolys;			Flip(ftri);			brush2.push_back(ftri);		}	}}void CMain::CSGIntersection(){	if(rootnode == NULL)	{		rootnode = new BspCsg;		TRI tri = brush1.front();		rootnode->plane.FromPoints(tri.p1, tri.p2, tri.p3);		rootnode->Build(brush1);		std::vector<TRI> fpolys;		std::vector<TRI> bpolys;		for(int i = 0; i < (int)brush2.size(); i++)            rootnode->Clip(brush2, fpolys, bpolys);		std::vector<TRI> tlist = bpolys;		delete rootnode;		rootnode = NULL;		rootnode = new BspCsg;		tri = brush2.front();		rootnode->plane.FromPoints(tri.p1, tri.p2, tri.p3);		rootnode->Build(brush2);		fpolys.clear();		bpolys.clear();		for(int i = 0; i < (int)brush1.size(); i++)            rootnode->Clip(brush1, fpolys, bpolys);		delete rootnode;		rootnode = NULL;		brush2.clear();		brush2 = tlist;		for(int i = 0; i < (int)bpolys.size(); i++)			brush2.push_back(bpolys);	}}void CMain::CSGUnion(){	if(rootnode == NULL)	{		rootnode = new BspCsg;		TRI tri = brush1.front();		rootnode->plane.FromPoints(tri.p1, tri.p2, tri.p3);		rootnode->Build(brush1);		std::vector<TRI> fpolys;		std::vector<TRI> bpolys;		for(int i = 0; i < (int)brush2.size(); i++)            rootnode->Clip(brush2, fpolys, bpolys);		std::vector<TRI> tlist = fpolys;		delete rootnode;		rootnode = NULL;		rootnode = new BspCsg;		tri = brush2.front();		rootnode->plane.FromPoints(tri.p1, tri.p2, tri.p3);		rootnode->Build(brush2);		fpolys.clear();		bpolys.clear();		for(int i = 0; i < (int)brush1.size(); i++)            rootnode->Clip(brush1, fpolys, bpolys);		delete rootnode;		rootnode = NULL;		brush2.clear();		brush2 = tlist;		for(int i = 0; i < (int)fpolys.size(); i++)			brush2.push_back(fpolys);	}}void CMain::Flip(TRI &tri){	sklib::XVec3 v = tri.p1;	tri.p1 = tri.p3;	tri.p3 = v;}

now i know the code is poorly written i wasnt expecting anyone else to read it

and u going have to change a few types like XVec3 to (D3DVECTOR3) or to ur own libray and so on

i hope this helps

w_poons

Mar 27, 2006 at 3:31 AM

Quote:
Original post by SKREAMZ
i did this about 2 years back using csg, i know your want something simpler but i cant help u there except by posting my code

thanks! i will have a deep look at your code after work.

w_poons

Mar 31, 2006 at 7:37 AM

after some thinking i came to the conclusion that i could solve the problem not in a geometrical way, because i'm only interested whats shown in the framebuffer. following solution came to my mind:

render the triangles and add 1 to the stencil buffer if it is occupied by a sphere. if the value is than greater as 1 this must be an area where spheres overlap. i than enable writing to the frame buffer where the value in the stencil buffer is less than 2.

what do you think? would this work?

edit:
ok, this approach has a problem. it leaves a hole where the spheres overlap ...

taby

Mar 31, 2006 at 8:08 AM

Quote:
Original post by SKREAMZ
i did this about 2 years back using csg, i know your want something simpler but i cant help u there except by posting my code
...

i hope this helps

Amazing. Thank you so much for sharing this!

taby

Mar 31, 2006 at 8:13 AM

Quote:
Original post by w_poons
after some thinking i came to the conclusion ...

what do you think? would this work?

edit:
ok, this approach has a problem. it leaves a hole where the spheres overlap ...

Advanced Graphics Programming Using OpenGL has an entire section (16.10) regarding CSG via stencil buffer. This book to me is an indespensible resource for just about every fundamental graphics algorithm I can think of.

I just wanted to point out that you may very well be going in the right direction, so to speak.

Topic Locked

This topic has been locked by a moderator. New replies are not allowed.