I am trying to implement one shot contact manifold generation for convex hulls using SAT, based on Dirk Gregorius' solution presented at GDC 2013 and 2015. I implemented the Edge-Edge detection case using the Gauss Map, as per the presentation. I also implemented the brute force method for this Edge-Edge case, where I project the hulls on the cross product of the edges, mostly to cross reference that my Gauss map version is correct. For most positions, they produce the same output, but for some configurations, the brute force method produces values that are greater - meaning smaller penetration. In order to debug this, I looked at the edges that the brute force found to see what the Gauss map did. The Gauss map method rejected the pair as not possible. Then I drew a debug mock gauss unit sphere with the normals involved, and they seemed to confirm the gauss map as not possible. I don't understand why that is, I have repeatedly double checked both methods and they seem correct.
This is the code for the Gauss map test (the AB_normal and CD_normal are simply the edge directions, in the hull construction I make sure to be in the same direction as the cross product of the face normals for that edge)
static SIMDVectorMask EdgeGaussMapTest(const Vector3& A, const Vector3& B, const Vector3& C, const Vector3& D, const Vector3& AB_normal, const Vector3& CD_normal) {
Vec8f zero_tolerance_vector = 0.0f;
// We still perform the normalization here since
// We can get some small values for the normals
// That would result in faulty values
Vector3 N_ab = Normalize(AB_normal);
Vector3 N_cd = Normalize(CD_normal);
Vec8f t1_1 = Dot(C, N_ab);
Vec8f t1_2 = Dot(D, N_ab);
Vec8f t2_1 = Dot(A, N_cd);
Vec8f t2_2 = Dot(B, N_cd);
Vec8f t1_val = t1_1 * t1_2;
Vec8f t2_val = t2_1 * t2_2;
Vec8f hemisphere_val = t2_2 * t1_1;
SIMDVectorMask t1_mask = t1_val < zero_tolerance_vector;
SIMDVectorMask t2_mask = t2_val < zero_tolerance_vector;
SIMDVectorMask hemisphere_mask = hemisphere_val > -zero_tolerance_vector;
SIMDVectorMask result = t1_mask && t2_mask && hemisphere_mask;
return result;
}This is the brute force method
float global_max_distance = -FLT_MAX;
float global_squared_edge_distance = FLT_MAX;
unsigned int global_first_index = -1;
unsigned int global_second_index = -1;
for (size_t index = 0; index < first->edges.size; index++) {
Line3D first_line = first->GetEdgePoints(index);
float3 edge = first_line.B - first_line.A;
for (size_t subindex = 0; subindex < second->edges.size; subindex++) {
Line3D second_line = second->GetEdgePoints(subindex);
float3 second_edge = second_line.B - second_line.A;
float3 cross = Cross(Normalize(edge), Normalize(second_edge));
if (!CompareMask(cross, float3::Splat(0.0f), float3::Splat(0.000001f))) {
float3 normalized_cross = Normalize(cross);
float2 first_projection_range = { FLT_MAX, -FLT_MAX };
for (size_t first_index = 0; first_index < first->vertex_size; first_index++) {
float dot = Dot(normalized_cross, first->GetPoint(first_index));
first_projection_range.x = min(dot, first_projection_range.x);
first_projection_range.y = max(dot, first_projection_range.y);
}
float2 second_projection_range = { FLT_MAX, -FLT_MAX };
for (size_t second_index = 0; second_index < second->vertex_size; second_index++) {
float dot = Dot(normalized_cross, second->GetPoint(second_index));
second_projection_range.x = min(dot, second_projection_range.x);
second_projection_range.y = max(dot, second_projection_range.y);
}
float separation = -FLT_MAX;
if (IsInRange(first_projection_range.x, second_projection_range.x, second_projection_range.y)) {
if (IsInRange(first_projection_range.y, second_projection_range.y, second_projection_range.y)) {
separation = -(first_projection_range.y - first_projection_range.x);
}
else {
separation = -(second_projection_range.y - first_projection_range.x);
}
}
else {
if (IsInRange(second_projection_range.x, first_projection_range.x, first_projection_range.y)) {
if (IsInRange(second_projection_range.y, first_projection_range.y, first_projection_range.y)) {
separation = -(second_projection_range.y - second_projection_range.x);
}
else {
separation = -(first_projection_range.y - second_projection_range.x);
}
}
else {
if (first_projection_range.y < second_projection_range.x) {
separation = first_projection_range.y - second_projection_range.x;
}
else {
separation = second_projection_range.y - first_projection_range.x;
}
}
}
if (separation > global_max_distance) {
global_max_distance = separation;
global_first_index = index;
global_second_index = subindex;
float squared_distance = SquaredDistanceBetweenSegmentPoints(first_line.A, first_line.B, second_line.A, second_line.B);
global_squared_edge_distance = squared_distance;
}
else if (FloatCompare(separation, global_max_distance, 0.000001f)) {
// It can happen that multiple edges have a similar or equal value
// Just to be sure, use a small epsilon to test for this equality
// In that case, choose the pair that has the least distance between
// The 2 edges
float squared_distance = SquaredDistanceBetweenSegmentPoints(first_line.A, first_line.B, second_line.A, second_line.B);
if (squared_distance < global_squared_edge_distance) {
global_first_index = index;
global_second_index = subindex;
global_squared_edge_distance = squared_distance;
}
}
}
}
}
return { global_max_distance, global_first_index, global_second_index };This is a screenshot with the brute force method, to see what pairs it finds
This is what the Gauss map finds (the edge-edge case is rejected in favor of a face-face contact)
The brute force one finds one of the edges that in the other case is from a face contact. To me the face manifold makes more sense, which would agree with the gauss map. Does anyone have any idea?
