Commit 9fd7c9f0 authored by Libor Moravčík's avatar Libor Moravčík
Browse files

Fixed barycentric function in geometry/utils.cpp

parent f25428c3
Loading
Loading
Loading
Loading
+6 −1
Original line number Diff line number Diff line
@@ -63,9 +63,14 @@ ClosestPoints closest_points_face_face(const std::vector<vec3> &vertices_a,
bool convex_test(const std::vector<vec3> &points,
		 const std::vector<std::array<size_t, 3>> &indices);

inline float triangle_area(scalar x1, scalar y1, scalar x2, scalar y2, scalar x3, scalar y3)
{
	return (x1 - x2) * (y2 - y3) - (x2 - x3) * (y1 - y2);
}

// Compute barycentric coordinates for
// point p with respect to triangle (a, b, c)
// code taken from (page 2): https://ceng2.ktu.edu.tr/~cakir/files/grafikler/Texture_Mapping.pdf
// code taken from (page 52): https://ceng2.ktu.edu.tr/~cakir/files/grafikler/Texture_Mapping.pdf
vec3 barycentric(const vec3 &a, const vec3 &b, const vec3 &c, const vec3 &p = {0, 0, 0});

scalar distance_point_triangle(const vec3 &p, const vec3 &v1, const vec3 &v2, const vec3 &v3);
+33 −14
Original line number Diff line number Diff line
@@ -325,20 +325,39 @@ bool convex_test(const std::vector<vec3> &points, const std::vector<std::array<s

vec3 barycentric(const vec3 &a, const vec3 &b, const vec3 &c, const vec3 &p)
{
	vec3 v0 = b - a, v1 = c - a, v2 = p - a;
	scalar d00 = dot(v0, v0);
	scalar d01 = dot(v0, v1);
	scalar d11 = dot(v1, v1);
	scalar d20 = dot(v2, v0);

	scalar d21 = dot(v2, v1);
	scalar denom = d00 * d11 - d01 * d01;

	scalar v, w, u;
	v = (d11 * d20 - d01 * d21) / denom;
	w = (d00 * d21 - d01 * d20) / denom;
	u = 1.0f - v - w;
	return {u, v, w};
	// Unnormalized triangle normal
	vec3 m = cross(b - a, c - a);
	// Nominators and one-over-denominator for u and v ratios
	float nu, nv, ood;
	// Absolute components for determining projection plane
	float x = std::abs(m.x), y = std::abs(m.y), z = std::abs(m.z);
	// Compute areas in plane of largest projection
	if (x >= y && x >= z) 
	{
		// x is largest, project to the yz plane
		nu = triangle_area(p.y, p.z, b.y, b.z, c.y, c.z); // Area of PBC in yz plane
		nv = triangle_area(p.y, p.z, c.y, c.z, a.y, a.z); // Area of PCA in yz plane
		ood = 1.0f / m.x; // 1/(2*area of ABC in yz plane
	} 
	else if (y >= x && y >= z)
	{
		// y is largest, project to the xz plane
		nu = triangle_area(p.x, p.z, b.x, b.z, c.x, c.z);
		nv = triangle_area(p.x, p.z, c.x, c.z, a.x, a.z);
		ood = 1.0f / -m.y;
	} else 
	{
		// z is largest, project to the xy plane
		nu = triangle_area(p.x, p.y, b.x, b.y, c.x, c.y);
		nv = triangle_area(p.x, p.y, c.x, c.y, a.x, a.y);
		ood = 1.0f / m.z;
	}

	vec3 result;
	result.x = nu * ood;
	result.y = nv * ood;
	result.z = 1.0f - result.x - result.y;
	return result;
}

scalar distance_point_triangle(const vec3 &p, const vec3 &v1, const vec3 &v2, const vec3 &v3)