🎉 Celebrating 25 Years of GameDev.net! 🎉

Not many can claim 25 years on the Internet! Join us in celebrating this milestone. Learn more about our history, and thank you for being a part of our community!

Ray to Octree Intersection for boolean geometry

Started by
48 comments, last by JoeJ 3 years, 4 months ago

Hi,

I'm trying to do boolean geometry that I have two meshes and would like to compute the intersection between them.

So I construct an octree from mesh A, and I check the vertices from mesh B against the octants, if there is an intersection, check the octant triangles for intersection, then I add the triangles, construct a mesh.

auto intersected_octant_faces = mn::buf_new<musa::Triangle>();
	std::stack < cg::Octant *> stack;
	stack.push(octree_sphere.root);
	for (size_t i = 0; i < tri_mesh_cube.vertices.count; i++)
	{
		while (!stack.empty())
		{
			cg::Octant* node = stack.top();
			if (!node)
				break;
			stack.pop();

			musa::Ray ray;
			ray.origin = { tri_mesh_cube.vertices[i] };
			ray.dir = { 1,0,0 };

			musa::Intersection_Points pts = {};
			pts = musa::ray_box_intersection(ray, node->region);
			if (pts.count >= 1)
			{
				musa::Intersection_Points t = {};
				auto vertices = node->faces;
				for (size_t j = 0; j < vertices.count; j += 3)
				{
					musa::Triangle tri{ vertices[j], vertices[j + 1], vertices[j + 2] };
					t = musa::ray_triangle_intersection(ray, tri);
					if (t.count == 1)
					{
						mn::buf_push(intersected_octant_faces, tri);
					}
				}

			}
			for (auto& n : node->octants)
			{
				stack.push(n);
			}
		}
	}

Right now I get just two faces as shown in the figure and not sure where is the problem.

Game Programming is the process of converting dead pictures to live ones .
Advertisement

The use of raytracing makes no sense to me here, looks like a logic mistake.

This would make sense for me: For each triangle of mesh B, make bounding box for this triangle, find all octree nodes intersecting this box, then do the triangle intersections.

Though, doing boolean operations on meshes is very difficult. It's hard to know what is inside / outside, and there are precision issues limiting robustness. Risk you invest some time but then have to give up is high. What do you try to achieve?

@joej

Thanks for your answer. Would you tell me after finding all the octree nodes interesecting the box which test to do ? Ray to triangle ? or Triangle to Triangle ?

Here is what I have done so far, of course not compiling..

and for splitting:

inline static void
	_cg_plane_split_polygon(CG_Plane& self, CG_Polygon& polygon, mn::Buf<CG_Polygon>& coplanarFront, mn::Buf<CG_Polygon>& coplanarBack, mn::Buf<CG_Polygon>& front, mn::Buf<CG_Polygon>& back)
	{
		enum
		{
			COPLANAR = 0,
			FRONT = 1,
			BACK = 2,
			SPANNING = 3
		};

		// Classify each point as well as the entire polygon into one of the above
		// four classes.
		int polygonType = 0;
		mn::Buf<int> types = mn::buf_new<int>();

		for (size_t i = 0; i < polygon.vertices.count; i++)
		{
			float t = dot(self.normal, polygon.vertices[i].pos) - self.w;
			int type = (t < -CG_EPSILON) ? BACK : ((t > CG_EPSILON) ? FRONT : COPLANAR);
			polygonType |= type;
			mn::buf_push(types, type);

		}

		// Put the polygon in the correct list, splitting it when necessary.
		switch (polygonType)
		{
			case COPLANAR:
			{
				if (dot(self.normal, polygon.plane.normal) > 0)
					mn::buf_push(coplanarFront, polygon);
				else
					mn::buf_push(coplanarBack, polygon);
				break;
			}
			case FRONT:
			{
				mn::buf_push(front, polygon);
				break;
			}
			case BACK:
			{
				mn::buf_push(back, polygon);
				break;
			}
			case SPANNING:
			{
				mn::Buf<CG_Vertex> f = mn::buf_with_allocator<CG_Vertex>(mn::memory::tmp());
				mn::Buf<CG_Vertex> b = mn::buf_with_allocator<CG_Vertex>(mn::memory::tmp());

				for (size_t i = 0; i < polygon.vertices.count; i++)
				{
					size_t j = (i + 1) % polygon.vertices.count;
					int ti = types[i], tj = types[j];
					CG_Vertex vi = polygon.vertices[i], vj = polygon.vertices[j];
					if (ti != BACK) mn::buf_push(f, vi);
					if (ti != FRONT) mn::buf_push(b, vi);
					if ((ti | tj) == SPANNING)
					{
						float t = (self.w - dot(self.normal, vi.pos)) / dot(self.normal, vj.pos - vi.pos);
						CG_Vertex v = _cg_vertex_interpolate(vi, vj, t);
						mn::buf_push(f, v);
						mn::buf_push(b, v);
					}
				}
				if (f.count >= 3) mn::buf_push(front, _cg_polygon_from_vertices(f));
				if (b.count >= 3) mn::buf_push(back, _cg_polygon_from_vertices(b));
				break;
			}
		}
		buf_free(types);
	}
	void split_octree_mesh(Octant* root, mn::Buf<CG_Polygon> ilist, mn::Buf<CG_Polygon>& front)
	{

		std::stack < cg::Octant*> stack;
		auto intersected_octant_faces = mn::buf_new<musa::Triangle>();

		for (size_t i = 0; i < ilist.count; i ++)
		{
		
			stack.push(root);
			while (!stack.empty())
			{
				cg::Octant* node = stack.top();
				stack.pop();


				musa::Triangle tri{ ilist[i].vertices[0].pos, ilist[i].vertices[1].pos,
					ilist[i].vertices[2].pos};
			
				if(box_box_intersect(musa::triangle_bounding_box(tri), node->region))
				{
					auto vertices = node->faces;
					musa::Intersection_Points t = {};
					musa::Ray ray;

					for (int a = 0; a < vertices.count; a += 3)
					{
						
						ray.origin = {vertices[a]};
						ray.dir = { 1,0,0 };

						t = musa::ray_triangle_intersection(ray, tri);
						if (t.count == 1)
						{
							mn::buf_push(intersected_octant_faces, tri);
						}
					}

				}
				 
				for (auto& n : node->octants)
				{
					if (n == nullptr)
						continue;
					stack.push(n);
				}
			}
		}
		 
	}
Game Programming is the process of converting dead pictures to live ones .

AhmedSaleh said:
Would you tell me after finding all the octree nodes interesecting the box which test to do ? Ray to triangle ? or Triangle to Triangle ?

Triangle to triangle. Rays are of no use here.

Looks you work on poly by poly clipping. That's surely a necessary detail, but maybe i can help with some thoughts from a more general point of view.

Starting with some potential requirements, issues, worries:

  1. For visual geometry, games usually use indexed meshes. They have issues like duplicated vertices on UV seams, non manifold, lacking definition of inside / outside, and most importantly lacking adjacency information.
    For geometry processing you typically want well defined meshes with adjacency so you can do things like e.g. iterating all edges around a vertex in clockwise order, find shortest paths, flip edges, insert new vertices, etc. The most widely used data structure seems ‘half-edge’, see here: https://www.graphics.rwth-aachen.de/media/openmesh_static/Documentations/OpenMesh-5.2-Documentation/a00016.html

(editor crashed, new post…)

I use this data structure as well, and i like it. But it's a lot of work to implement, getting used to it, and then implement mesh editing operations like edge split, flip, etc. Multiple weeks of work, depending on what you need. The linked library is free IIRC, and it should have all this functionality already. Maybe worth to consider. There are more such libraries, e.g. CGAL which also has boolean mesh operations. https://www.cgal.org/​ No idea about license. To deal with floating point issues CGAL uses some extended precision representation of numbers. (Your use of epsilon will still leave you with some very close vertices which should be merged, and this is where the rabbit hole opens up.)

2. Let's assume you want to keep things simple: ‘Mesh adjacency should not be necessary - just cut those polygons along the intersections. Should not be that hard, and if there are some errors in the result it does not matter.’
Not sure if this can work because i've never done it. But i'm optimistic and would do it like this:

Find intersection polygons, split them so the new edges are added to both meshes. (As you currently do, i guess)
After that, determinate which polygons should be deleted.

The second step is maybe harder than you think. But i know of a very robust method, which is also easy to implement. I use it to voxelize non manifold meshes with holes, or confusing ‘inside / outside’ situations.
(E.g. think of a column on a planar floor. The modeler guy will not cut out a hole into the floor where the column is. He will just place a cylinder (column) on a big quad (floor). This gives us a problem, because inside the column we see front faces of the floor but back faces of the column. So we do not know if we are in solid or empty space - only human intelligence creates this information for us, but writing an algorithm becomes hard. )

So, the method i propose does something that i call ‘solid angle integration’. It basically gives you the best guess we can do without requiring a form of intelligence. (A very similar method has been described in papers of Alec Jacobson https://www.dgp.toronto.edu/projects/fast-winding-numbers/​ but they somehow did it more complicated than necessary)

Here is the code where i figured it out:

static bool visInsideOutsideTest = 0; ImGui::Checkbox("visInsideOutsideTest", &visInsideOutsideTest);
		if (visInsideOutsideTest)
		{
			static int init = 1;
			static HEMesh mesh;
			if (init) 
			{
				//MeshConverter::TestTrianglesToHE (mesh, mesh);
				//((HEMesh_Serialization&)mesh).AutoLoadModel ("C:\\dev\\pengII\\mod\\bunny closed.lxo"); // hangs on tessellation (delauney)
				((HEMesh_Serialization&)mesh).AutoLoadModel ("C:\\dev\\pengII\\mod\\Armadillo_input.obj");
				//((HEMesh_Serialization&)mesh).AutoLoadModel ("C:\\dev\\pengII\\mod\\Gear.lxo");
				//((HEMesh_Serialization&)mesh).AutoLoadModel ("C:\\dev\\pengII\\mod\\MaxPlanck_input.obj");
				//((HEMesh_Serialization&)mesh).AutoLoadModel ("C:\\dev\\pengII\\mod\\Beetle_input.obj"); // hangs on tessellation (delauney)
				//((HEMesh_Serialization&)mesh).AutoLoadModel ("C:\\dev\\pengII\\mod\\Buddha_input.obj"); // hangs on tessellation (delauney)
			}

			static vec query(0.46f, 0, 0);
			ImGui::DragFloat3("query", (float*)&query, 0.001f);
			static bool surfels = 0; ImGui::Checkbox("surfels", &surfels);

			float integral = 0;

			if (surfels) for (int i=0; i<mesh.GetPolyCount(); i++) // surfels
			{
				float area;
				vec center = mesh.CalcPolyCenterAndArea(area, i);
				vec norm = mesh.CalcPolyNormal(i);

				vec diff = query - center;
				float sql = diff.SqL(); // SqL is just dot(diff,diff)
				if (sql > FP_EPSILON2)	
				{
					float solidAngle = diff.Dot(norm * area) / (sqrt(sql) * sql);
					integral += solidAngle;
				}
			}

			if (!surfels) for (int i=0; i<mesh.GetPolyCount(); i++) // triangles
			{
				vec triVD[3];
				bool onSurface = false;

				int eI = mesh.GetPolyEdge(i);
				for (int j=0; j<3; j++)
				{
					vec diff = mesh.GetEdgeVertexPos(eI) - query;
					float dist = diff.Length();
					if (dist < FP_EPSILON2)
					{
						onSurface = true;
						break;
					}
					triVD[j] = diff / dist;
					eI = mesh.GetEdgeNext(eI);
				}

				if (!onSurface)
				{
					float numer = triVD[0].Dot(
						vec(triVD[1]-triVD[0]).Cross(triVD[2]-triVD[0]));
					float denom = 1 + 
						triVD[0].Dot(triVD[1]) + 
						triVD[0].Dot(triVD[2]) + 
						triVD[1].Dot(triVD[2]);

					float solidAngle = 2 * atan2(numer, denom);
					integral += solidAngle;
				}
			}

			float winding = 1 / float(4*PI) * integral;
			bool inside = winding > 0.5f;
			RenderCircle(0.01f, query, vec(0), !inside, inside, 0);
			ImGui::Text("winding %f", winding);

			((HEMesh_Vis&)mesh).RenderMesh(0,0,0);


			init = 0;
			
		}

We can use a query point, e.g. the center of a triangle from mesh A, and see if it is inside or outside of mesh B. Then delete accordingly.

Then i have two options. First option is fast by approximating each triangle with a disc ('surfel'). So i precompute those surfels to be on the triangle center, having equal normal and area. This gives approximate results, but good enough if the triangles are not close to each other.

Second option is exact by using triangles instead surfels.

This works pretty well, but i still need some further heuristices to get a solid voxelization of the Sponza model, for example. For some meshes which are even worse, user interaction is necessary to decide which parts should be removed.

@JoeJ Thanks a lot for your help and for your ideas.

The problem is that I have to stick with octree implementation. Would you tell me the correct algorithm for doing it ? What are your comments about the implementation detail that I wrote ? I still can't get it working correctly.

For example triangle to triangle and the result of it how would I use it ?

Game Programming is the process of converting dead pictures to live ones .

Here is a rebuild of your idea

void split_octree_mesh(Octant* root, mn::Buf<CG_Polygon> ilist, mn::Buf<CG_Polygon>& front)
	{

		std::stack < cg::Octant*> stack;
		auto intersected_octant_faces = mn::buf_new<musa::Triangle>();

		for (size_t i = 0; i < ilist.count; i ++)
		{
		
			stack.push(root);
			while (!stack.empty())
			{
				cg::Octant* node = stack.top();
				stack.pop();


				musa::Triangle tri{ ilist[i].vertices[0].pos, ilist[i].vertices[1].pos,
					ilist[i].vertices[2].pos};
			
				if(box_box_intersect(musa::triangle_bounding_box(tri), node->region))
				{
					auto vertices = node->faces;
				
					for (size_t a = 0; a < vertices.count; a += 3)
					{
						musa::Triangle verTri;
						verTri.p0 = vertices[a];
						verTri.p1 = vertices[a + 1];
						verTri.p2 = vertices[a + 1];
						
						if (triangle_triangle_intersection_check(tri, verTri))
						{
							mn::buf_push(intersected_octant_faces, tri);
						}
					}

				}
				 
				for (auto& n : node->octants)
				{
					if (n == nullptr)
						continue;
					stack.push(n);
				}
			}
		}
		 
	}
Game Programming is the process of converting dead pictures to live ones .

Here is full code:

	        auto tri_mesh_a = cg::trimesh_from_indexed_mesh(a);
			auto tri_mesh_b = cg::trimesh_from_indexed_mesh(b);
			auto octree_a = octree_from_mesh(tri_mesh_a);
			auto octree_b = octree_from_mesh(tri_mesh_b);
			return _cg_operation(a, b, octree_b.root, octree_a.root, _cg_intersect);

	
	inline static CG_Node*
	_cg_node_from_polygons(mn::Buf<CG_Polygon> list, Octant * root)
	{

		CG_Node* self = _cg_node_new();
		_cg_node_build(*self, list, root);
		return self;
	}


	// Algorithm.
	typedef CG_Node* _cg_function(CG_Node* a1, CG_Node* b1, Octant* rootA, Octant* rootB);

	inline static cg::Indexed_Mesh
	_cg_operation(const cg::Indexed_Mesh& a, const cg::Indexed_Mesh& b, Octant *rootA, Octant *rootB, _cg_function fun)
	{

		CG_Node* A = _cg_node_from_polygons(_cg_model_to_polygons(a), rootA);
		CG_Node* B = _cg_node_from_polygons(_cg_model_to_polygons(b), rootB);
		CG_Node* AB = fun(A, B, rootA, rootB);
		mn::Buf<CG_Polygon> polygons = _cg_node_all_polygons(*AB);
	/*	_cg_node_free(*A);
		buf_free(A->polygons);

		A = 0;
		_cg_node_free(*B);
		buf_free(B->polygons);
		B = 0;
		_cg_node_free(*AB);
		buf_free(AB->polygons);

		AB = 0;
		*/
		return _cg_model_from_polygons(polygons);
	}

	inline static CG_Node*
	_cg_intersect(CG_Node* a1, CG_Node* b1, Octant *rootA, Octant *rootB)
	{

	
		CG_Node* a = _cg_node_clone(*a1);
		CG_Node* b = _cg_node_clone(*b1);
		_cg_node_invert(*a);
		_cg_node_clip_to(*b, a);
		_cg_node_invert(*b);
		_cg_node_clip_to(*a, b);
		_cg_node_clip_to(*b, a);
		_cg_node_build(*a, _cg_node_all_polygons(*b), rootA);
		_cg_node_invert(*a);
		CG_Node* ret = _cg_node_from_polygons(_cg_node_all_polygons(*a), rootB);
		mn::buf_free(a->polygons);
		_cg_node_free(*a);
		mn::buf_free(b->polygons);
		_cg_node_free(*b);
	
		
		return ret;
	}
inline static void
	_cg_node_invert(CG_Node& self)
	{
		mn::Deque<CG_Node*> nodes = mn::deque_new<CG_Node*>();
		mn::deque_push_back(nodes, &self);

		while (nodes.count)
		{
			CG_Node* me = mn::deque_front(nodes);
			mn::deque_pop_front(nodes);

			for (size_t i = 0; i < me->polygons.count; i++)
				_cg_polygon_flip(me->polygons[i]);
			_cg_plane_flip(me->plane);
			std::swap(me->front, me->back);
			if (me->front)
				mn::deque_push_back(nodes, me->front);
			if (me->back)
				mn::deque_push_back(nodes, me->back);
		}
		deque_free(nodes);
	}

	inline static mn::Buf<CG_Polygon>
	_cg_node_clip_polygons(CG_Node& self, mn::Buf<CG_Polygon> ilist)
	{
		mn::Buf<CG_Polygon> result = mn::buf_with_allocator<CG_Polygon>(mn::memory::tmp());
		auto clips = mn::deque_new<std::pair<CG_Node*, mn::Buf<CG_Polygon>>>();
		mn::deque_push_back(clips, std::make_pair(&self, ilist));

		while (clips.count)
		{
			CG_Node* me = mn::deque_front(clips).first;
			mn::Buf<CG_Polygon> list = mn::deque_front(clips).second;
			deque_pop_back(clips);

			if (!_cg_plane_ok(me->plane))
			{
				for (int i = 0; i < list.count; i++)
				{
					mn::buf_push(result, list[i]);
				}
				continue;
			}

			mn::Buf<CG_Polygon> list_front = mn::buf_with_allocator<CG_Polygon>(mn::memory::tmp());
			mn::Buf<CG_Polygon> list_back  = mn::buf_with_allocator<CG_Polygon>(mn::memory::tmp());

			for (size_t i = 0; i < list.count; i++)
				_cg_plane_split_polygon(me->plane, list[i], list_front, list_back, list_front, list_back);

			if (me->front)
				deque_push_back(clips,std::make_pair(me->front, list_front));
			else
			{
				for (int i = 0; i < list_front.count; i++)
				{
					mn::buf_push(result, list_front[i]);
				}
			}

			if (me->back)
				deque_push_back(clips, std::make_pair(me->back, list_back));

		}
		mn::deque_free(clips);
		mn::buf_free(ilist);
		return result;
	}

	// Remove all polygons in this BSP tree that are inside the other BSP tree
	inline static void
	_cg_node_clip_to(CG_Node& self, CG_Node* other)
	{
		mn::Deque<CG_Node*> nodes = mn::deque_new<CG_Node*>();
		mn::deque_push_back(nodes, &self);

		while (nodes.count)
		{
			CG_Node* me = mn::deque_front(nodes);
			mn::deque_pop_front(nodes);

			me->polygons = _cg_node_clip_polygons(*other, me->polygons);

			if (me->front)
				mn::deque_push_back(nodes, me->front);

			if (me->back)
				mn::deque_push_back(nodes, me->back);
		}
		mn::deque_free(nodes);
	}

	inline static mn::Buf<CG_Polygon>
	_cg_node_all_polygons(CG_Node& self)
	{
		mn::Buf<CG_Polygon> result = mn::buf_new<CG_Polygon>();

		auto nodes = mn::deque_new<CG_Node*>();
		mn::deque_push_back(nodes, &self);

		while (nodes.count)
		{
			auto me = mn::deque_front(nodes);
			mn::deque_pop_front(nodes);

			for (int i = 0; i < me->polygons.count; i++)
			{
				mn::buf_push(result, me->polygons[i]);

			}
			if (me->front)
				mn::deque_push_back(nodes, me->front);
			if (me->back)
				mn::deque_push_back(nodes, me->back);
		}

		mn::deque_free(nodes);
		return result;
	}

The current work out and problems I think here:

void split_octree_mesh(Octant* root, mn::Buf<CG_Polygon> ilist, mn::Buf<CG_Polygon> &list_front)
	{

		std::stack < cg::Octant*> stack;

		for (size_t i = 0; i < ilist.count; i ++)
		{
		
			stack.push(root);
			while (!stack.empty())
			{
				cg::Octant* node = stack.top();
				stack.pop();


				musa::Triangle tri{ ilist[i].vertices[0].pos, ilist[i].vertices[1].pos,
					ilist[i].vertices[2].pos};
				
				if(box_box_intersect(musa::triangle_bounding_box(tri), node->region))
				{
					auto vertices = node->faces;

					cg::CG_Polygon poly;
					poly.vertices = mn::buf_new<CG_Vertex>();
					for (size_t a = 0; a < vertices.count; a += 3)
					{
						musa::Triangle verTri;
						verTri.p0 = vertices[a];
						verTri.p1 = vertices[a + 1];
						verTri.p2 = vertices[a + 1];

						CG_Vertex vertex1 = { vertices[a] };
						CG_Vertex vertex2 = { vertices[a + 1] };
						CG_Vertex vertex3 = { vertices[a + 2] };
					
						if (triangle_triangle_intersection_check(tri, verTri))
						{
							mn::buf_push(poly.vertices, vertex1);
							mn::buf_push(poly.vertices, vertex2);
							mn::buf_push(poly.vertices, vertex3);

							buf_push(list_front, poly);
						}
					
					}

				}
				 
				for (auto& n : node->octants)
				{
					if (n == nullptr)
						continue;
					stack.push(n);
				}
			}
		}
		 
	}
Game Programming is the process of converting dead pictures to live ones .

AhmedSaleh said:
The current work out and problems I think here:

hmm… it's way too much complexity i could help much just by looking at the code. The only thing catches my eye is about efficiency: You should visit the children of a node only if the parent intersects the query box (bounding box of a triangle). Currently using an octree does not give you any benefit, as you always process the whole tree for each triangle.

It's not clear to me what you try to do, and what problem you have exactly. You would need to elaborate in more detail. I see ‘BSP’ in the comments. Do you work on some low poly level editor like Quake? If so, BSP is also commonly used to do CSG, so why don't you convert the models to BSP and do the boolean stuff with BSP trees? No need for octree then.

@joej

I did it with bsp, now I wanna do it using octree for practices. Would you elaborate more about the efficacy code ?

Game Programming is the process of converting dead pictures to live ones .

This topic is closed to new replies.

Advertisement