Skip to content

038 Ray‐triangle intersection

badcoiq edited this page Sep 16, 2023 · 2 revisions

Ray-triangle intersection

Определяем пересекает ли луч треугольник.

Добавлю класс в котором нужно указать три координаты, и в нём будут нужные методы.

class bqTriangle
{
public:
	bqTriangle() {}
	bqTriangle(const bqVec4_t<bqReal>& _v1, const bqVec4_t<bqReal>& _v2, const bqVec4_t<bqReal>& _v3)
		:
		v1(_v1),
		v2(_v2),
		v3(_v3)
	{
	}

	bqTriangle(const bqVec3f& _v1, const bqVec3f& _v2, const bqVec3f& _v3)
		:
		v1(_v1),
		v2(_v2),
		v3(_v3)
	{
	}

	bqVec4_t<bqReal> v1;
	bqVec4_t<bqReal> v2;
	bqVec4_t<bqReal> v3;
	bqVec4_t<bqReal> e1;
	bqVec4_t<bqReal> e2;

	void Update()
	{
		e1 = bqVec4_t<bqReal>(v2.x - v1.x,
			v2.y - v1.y,
			v2.z - v1.z,
			0.f);
		e2 = bqVec4_t<bqReal>(v3.x - v1.x,
			v3.y - v1.y,
			v3.z - v1.z,
			0.f);
		//	e1.cross(e2, faceNormal);
	}

	void Center(bqVec4_t<bqReal>& out)
	{
		out = (v1 + v2 + v3) * 0.3333333;
	}

	bool RayTest_MT(const bqRay& ray, bool withBackFace, bqReal& T, bqReal& U, bqReal& V, bqReal& W)
	{
		bqVec4_t<bqReal>  pvec;
		bqMath::Cross(ray.m_direction, e2, pvec);
		bqReal det = e1.Dot(pvec);

		if (withBackFace)
		{
			if (std::fabs(det) < bqEpsilon)
				return false;
		}
		else
		{
			if (det < bqEpsilon && det > -bqEpsilon)
				return false;
		}

		bqVec4_t<bqReal> tvec(
			ray.m_origin.x - v1.x,
			ray.m_origin.y - v1.y,
			ray.m_origin.z - v1.z,
			0.f);

		bqReal inv_det = 1.f / det;
		U = tvec.Dot(pvec) * inv_det;

		if (U < 0.f || U > 1.f)
			return false;

		bqVec4_t<bqReal>  qvec;
		bqMath::Cross(tvec, e1, qvec);

		V = ray.m_direction.Dot(qvec) * inv_det;

		if (V < 0.f || U + V > 1.f)
			return false;

		T = e2.Dot(qvec) * inv_det;

		if (T < bqEpsilon) return false;

		W = 1.f - U - V;
		return true;
	}

	bool RayTest_Watertight(const bqRay& ray, bool withBackFace, bqReal& T, bqReal& U, bqReal& V, bqReal& W)
	{
		v1.w = 1.f;
		v2.w = 1.f;
		v3.w = 1.f;
		const auto A = v2 - ray.m_origin;
		const auto B = v3 - ray.m_origin;
		const auto C = v1 - ray.m_origin;

		const bqReal Ax = A[ray.m_kx] - (ray.m_Sx * A[ray.m_kz]);
		const bqReal Ay = A[ray.m_ky] - (ray.m_Sy * A[ray.m_kz]);
		const bqReal Bx = B[ray.m_kx] - (ray.m_Sx * B[ray.m_kz]);
		const bqReal By = B[ray.m_ky] - (ray.m_Sy * B[ray.m_kz]);
		const bqReal Cx = C[ray.m_kx] - (ray.m_Sx * C[ray.m_kz]);
		const bqReal Cy = C[ray.m_ky] - (ray.m_Sy * C[ray.m_kz]);

		U = (Cx * By) - (Cy * Bx);
		V = (Ax * Cy) - (Ay * Cx);
		W = (Bx * Ay) - (By * Ax);

		if (U == 0.f || V == 0.f || W == 0.f)
		{
			double CxBy = (double)Cx * (double)By;
			double CyBx = (double)Cy * (double)Bx;
			U = (float)(CxBy - CyBx);

			double AxCy = (double)Ax * (double)Cy;
			double AyCx = (double)Ay * (double)Cx;
			V = (float)(AxCy - AyCx);

			double BxAy = (double)Bx * (double)Ay;
			double ByAx = (double)By * (double)Ax;
			W = (float)(BxAy - ByAx);
		}

		if (withBackFace)
		{
			if ((U < 0.f || V < 0.f || W < 0.f) &&
				(U > 0.f || V > 0.f || W > 0.f))
				return false;
		}
		else
		{
			if (U < 0.f || V < 0.f || W < 0.f)
				return false;
		}

		bqReal det = U + V + W;

		if (det == 0.f)
			return false;

		const bqReal Az = ray.m_Sz * A[ray.m_kz];
		const bqReal Bz = ray.m_Sz * B[ray.m_kz];
		const bqReal Cz = ray.m_Sz * C[ray.m_kz];
		const bqReal Ts = (U * Az) + (V * Bz) + (W * Cz);

		if (!withBackFace) // CULL
		{
			if (Ts < 0.f || Ts > bqInfinity * det)
				return false;
		}
		else
		{
			if (det < 0.f && (Ts >= 0.f || Ts < bqInfinity * det))
				return false;
			else if (det > 0.f && (Ts <= 0.f || Ts > bqInfinity * det))
				return false;
		}

		const bqReal invDet = 1.f / det;
		U = U * invDet;
		V = V * invDet;
		W = W * invDet;
		T = Ts * invDet;
		if (T < bqEpsilon)
			return false;
		return true;
	}
};

Есть два алгоритма, MT самый быстрый, Watertight самый точный.

Примеры из демо. Проходимся по bqPolygonMesh

bool ExampleBasicsRayTri::_getTriangle(bqTriangle* out, bqRay* ray, bqVec4& ip)
{
	bool ret = false;
	bqReal len = 9999.0;
	if (m_polygonMesh->m_polygons.m_head)
	{
		for (auto p : m_polygonMesh->m_polygons)
		{
			auto V1 = p->m_vertices.m_head;
			auto last = V1->m_left;
			while (true)
			{
				auto V2 = V1->m_right;
				auto V3 = V2->m_right;

				bqTriangle tri;
				tri.v1 = V1->m_data->m_data.BaseData.Position;
				tri.v2 = V2->m_data->m_data.BaseData.Position;
				tri.v3 = V3->m_data->m_data.BaseData.Position;
				tri.Update();

				bqReal T = 0.f;
				bqReal U = 0.f;
				bqReal V = 0.f;
				bqReal W = 0.f;
				ray->Update();
				if (tri.RayTest_MT(*ray, true, T, U, V, W))
				{
					if (T < len)
					{
						ip = ray->m_origin + (T * ray->m_direction);
						len = T;
						*out = tri;
						ret = true;
					}
				}

				if (V1 == last)
					break;
				V1 = V1->m_right;
			}
		}
	}
	return ret;
}

Проходимся по bqMesh

  • оптимизация, сначала проверяем пересекает ли луч AABB
bool ExampleBasicsRayTri2::_getRayHit(bqAabb& aabb, bqTriangle& outTri, bqRay& ray, bqVec4& ip)
{
	bool ret = false;
	bqReal len = 9999.0;

	for (size_t i = 0; i < m_model->m_meshBuffers.m_size; ++i)
	{
		auto inf = m_model->m_meshBuffers.m_data[i]->m_CPUMmesh->GetInfo();
		if (inf.m_aabb.RayTest(ray))
		{
			auto VB = m_model->m_meshBuffers.m_data[i]->m_CPUMmesh->GetVBuffer();
			auto IB = m_model->m_meshBuffers.m_data[i]->m_CPUMmesh->GetIBuffer();
			
			bqVertexTriangle* vt = (bqVertexTriangle*)VB;
			bqVertexTriangleSkinned* vts = (bqVertexTriangleSkinned*)VB;
			uint32_t* i32 = (uint32_t*)IB;
			uint16_t* i16 = (uint16_t*)IB;

			for (uint32_t ii = 0; ii < inf.m_iCount; )
			{
				bqTriangle tri;

				uint32_t ind1 = 0;
				uint32_t ind2 = 0;
				uint32_t ind3 = 0;

				if (inf.m_indexType == bqMeshIndexType::u16)
				{
					ind1 = i16[ii];
					ind2 = i16[ii+1];
					ind3 = i16[ii+2];
				}
				else
				{
					ind1 = i32[ii];
					ind2 = i32[ii + 1];
					ind3 = i32[ii + 2];
				}

				if (inf.m_skinned)
				{
					tri.v1 = vts[ind1].BaseData.Position;
					tri.v2 = vts[ind2].BaseData.Position;
					tri.v3 = vts[ind3].BaseData.Position;
				}
				else
				{
					tri.v1 = vt[ind1].Position;
					tri.v2 = vt[ind2].Position;
					tri.v3 = vt[ind3].Position;
				}

				ii += 3;

				tri.Update();

				bqReal T = 0.f;
				bqReal U = 0.f;
				bqReal V = 0.f;
				bqReal W = 0.f;
				ray.Update();
				if (tri.RayTest_Watertight(ray, true, T, U, V, W))
				{
					if (T < len)
					{
						ip = ray.m_origin + (T * ray.m_direction);
						len = T;
						outTri = tri;
						ret = true;
						aabb = inf.m_aabb;
					}
				}
			}
		}
	}

	return ret;
}

Значение T это расстояние до точки пересечения. Чтобы получить саму точку нужно взять начало луча и прибавить (T * направление)

intersectionPoint = ray.m_origin + (T * ray.m_direction);

https://github.com/badcoiq/badcoiq/tree/1b8a223b4d0cad35e9e1b31573ac2bdb471c5c7c