Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Testing Various Cases, and Potential fix #24

Draft
wants to merge 2 commits into
base: master
Choose a base branch
from
Draft

Testing Various Cases, and Potential fix #24

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
25385e2
Testing Various Cases, and potential fix
Kushal-Shah-03 Jul 10, 2024
08452dc
Merge branch 'master' into master
Kushal-Shah-03 Jul 10, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 5 additions & 2 deletions MathUtils.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -21,7 +21,8 @@ namespace quickhull {
// (divide by N.getNormalized() to get actual Euclidian distance).
template <typename T>
inline T getSignedDistanceToPlane(const Vector3<T>& v, const Plane<T>& p) {
return p.m_N.dotProduct(v) + p.m_D;
return p.m_N.dotProduct(v) + p.m_D;
// return (p.m_N.getNormalized().dotProduct(v) + p.m_D);
}

template <typename T>
Expand All @@ -36,7 +37,9 @@ namespace quickhull {
T px = y * rhsz - z * rhsy ;
T py = z * rhsx - x * rhsz ;
T pz = x * rhsy - y * rhsx ;
return Vector3<T>(px,py,pz);
// errorNormal = Vector3<T>(std::abs(y*rhsz)+std::abs(z*rhsy),std::abs(z*rhsx)+std::abs(x*rhsz),std::abs(x*rhsy)+std::abs(y*rhsx));
// return Vector3<T>(px,py,pz);
return Vector3<T>(px,py,pz).getNormalized();
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think this is your proposed fix?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, so basically normalizing the output of getTriangleNormal function, since we were comparing the distance (unnormalized) directly to epsilon in

quickhull/QuickHull.cpp

Line 462 in 4ef66c6

if (d>maxD) {
[maxD was initially set to m_epsilon]
and based on my understanding of the code that wasn't intended.

Also, the above line and

quickhull/QuickHull.hpp

Lines 196 to 206 in 4ef66c6

const T D = mathutils::getSignedDistanceToPlane(m_vertexData[ pointIndex ],f.m_P);
if (D>0 && D*D > m_epsilonSquared*f.m_P.m_sqrNLength) {
if (!f.m_pointsOnPositiveSide) {
f.m_pointsOnPositiveSide = std::move(getIndexVectorFromPool());
}
f.m_pointsOnPositiveSide->push_back( pointIndex );
if (D > f.m_mostDistantPointDist) {
f.m_mostDistantPointDist = D;
f.m_mostDistantPoint = pointIndex;
}
return true;

Both compare distances to find the maximum distance, in such a situation when the distance is not the Euclidean distance the comparison could select an "incorrect point" (Now I know that we had discussed that choosing the absolute maximum point doesn't make a difference in floating point arithmetic, but it's possible that due to this non-Euclidean distance we pick a point which is very close, which could lead to the errors [this is just speculation on my end] leading to those errors)

PS: I'm sorry if I wasn't able to explain this properly in the original comment. I hope it's clearer now

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good find! Since normalization is somewhat expensive, you might want to check all the places this function is called. If they all need to be normalized, then this is perfect, but if not, it would be ideal to only add normalization where needed. Have you noticed any difference in speed from this change?

}


Expand Down
85 changes: 82 additions & 3 deletions QuickHull.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,31 @@ namespace quickhull {
return ConvexHull<T>(m_mesh,m_vertexData, CCW, useOriginalIndices);
}

// Check if the mesh remains convex after adding new faces
template<typename T>
bool QuickHull<T>::isMeshConvex(const std::vector<size_t>& newFaces) {
for (const size_t& faceIndex : newFaces) {
const auto& face = m_mesh.m_faces[faceIndex];
bool positive = false;
bool negative = false;
for (auto i : tempPoints) {
Vector3 vertex = m_vertexData[i];
double distance = face.m_P.m_N.dotProduct(vertex) + face.m_P.m_D;
if (distance > m_epsilon) {
positive = true;
} else if (distance < -m_epsilon) {
negative = true;
}

// If we find both positive and negative distances, the polyhedron is not convex
if (positive && negative) {
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should there be any positive distances for a valid convex hulll?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yeah, so the way I'm checking it is if all the other points that are a part of the convex hull, don't lie on the same side of the face it's not a convex hull. I wrote this to account for incorrect orientation of the normal, since I am testing it after each iteration, and I wasn't sure whether the normals were oriented correctly after each step or just at the last step.

Although, now that you've pointed it out, I think I should have tried checking only for the negatives first in case the normals are oriented correctly, it will be a more sound check. I will check the code to verify the normals aren't rearranged later and if it works I'll make the changes. Thanks for pointing it out

return false;
}
}
}
return true;
}

template<typename T>
void QuickHull<T>::createConvexHalfEdgeMesh() {
m_visibleFaces.clear();
Expand All @@ -113,6 +138,7 @@ namespace quickhull {

// Process faces until the face list is empty.
size_t iter = 0;
int incorrect_count=0;
while (!m_faceList.empty()) {
iter++;
if (iter == std::numeric_limits<size_t>::max()) {
Expand All @@ -135,8 +161,16 @@ namespace quickhull {
}

// Pick the most distant point to this triangle plane as the point to which we extrude
// int breakFlag=0;
const vec3& activePoint = m_vertexData[tf.m_mostDistantPoint];
// std::cout << "Active Point" << activePoint.x << " " << activePoint.y << " " << activePoint.z << std::endl;
const size_t activePointIndex = tf.m_mostDistantPoint;
// std::cout << activePointIndex << std::endl;
// if (activePointIndex==0)
// {
// breakFlag=1;
// }
tempPoints.push_back(activePointIndex);

// Find out the faces that have our active point on their positive side
// (these are the "visible faces"). The face on top of the stack of course is
Expand All @@ -161,6 +195,9 @@ namespace quickhull {
pvf.m_visibilityCheckedOnIteration = iter;
const T d = P.m_N.dotProduct(activePoint)+P.m_D;
if (d>0) {
// if (d>m_epsilon && d*d > m_epsilonSquared*P.m_sqrNLength) {
// if (d>0 && d*d > m_epsilonSquared*P.m_sqrNLength) {
// if (d>m_epsilon) {
pvf.m_isVisibleFaceOnCurrentIteration = 1;
pvf.m_horizonEdgesOnCurrentIteration = 0;
m_visibleFaces.push_back(faceData.m_faceIndex);
Expand Down Expand Up @@ -250,7 +287,9 @@ namespace quickhull {
A = horizonEdgeVertexIndices[0];
B = horizonEdgeVertexIndices[1];
C = activePointIndex;

// std::cout << "PtA" << m_vertexData[A].x << " " << m_vertexData[A].y << " " << m_vertexData[A].z << std::endl;
// std::cout << "PtB" << m_vertexData[B].x << " " << m_vertexData[B].y << " " << m_vertexData[B].z << std::endl;
// std::cout << "PtC" << m_vertexData[C].x << " " << m_vertexData[C].y << " " << m_vertexData[C].z << std::endl;
const size_t newFaceIndex = m_mesh.addFace();
m_newFaceIndices.push_back(newFaceIndex);

Expand All @@ -272,6 +311,7 @@ namespace quickhull {

const Vector3<T> planeNormal = mathutils::getTriangleNormal(m_vertexData[A],m_vertexData[B],activePoint);
newFace.m_P = Plane<T>(planeNormal,activePoint);
// std::cout << "N" << planeNormal.x << " " << planeNormal.y << " " << planeNormal.z << std::endl;
newFace.m_he = AB;

m_mesh.m_halfEdges[CA].m_opp = m_newHalfEdgeIndices[i>0 ? i*2-1 : 2*horizonEdgeCount-1];
Expand All @@ -292,6 +332,9 @@ namespace quickhull {
}
}
// The points are no longer needed: we can move them to the vector pool for reuse.
// for (auto disabledPt : *disabledPoints) {
// std::cout << "Disabled Point" << m_vertexData[disabledPt].x << " " << m_vertexData[disabledPt].y << " " << m_vertexData[disabledPt].z << std::endl;
// }
reclaimToIndexVectorPool(disabledPoints);
}

Expand All @@ -306,8 +349,27 @@ namespace quickhull {
}
}
}
// if (isMeshConvex(m_newFaceIndices))
// {
// // if (breakFlag==1)
// // {
// // m_indexVectorPool.clear();
// // break;
// // }
// // std::cout << "YAY\n";
// }
// else
// {
// incorrect_count+=1;
// // std::cout << "OOP\n";
// // m_indexVectorPool.clear();
// // break;
// }
}

if (!isMeshConvex(m_newFaceIndices))
{
std::cout << "Invalid Output" << "\n";
}
// Cleanup
m_indexVectorPool.clear();
}
Expand Down Expand Up @@ -400,6 +462,10 @@ namespace quickhull {
if (trianglePlane.isPointOnPositiveSide(m_vertexData[v[3]])) {
std::swap(v[0],v[1]);
}
tempPoints.push_back(v[0]);
tempPoints.push_back(v[1]);
tempPoints.push_back(v[2]);
tempPoints.push_back(v[3]);
return m_mesh.setup(v[0],v[1],v[2],v[3]);
}

Expand Down Expand Up @@ -431,6 +497,7 @@ namespace quickhull {
if (distToRay > maxD) {
maxD=distToRay;
maxI=i;
// Maybe here
}
}
if (maxD == m_epsilonSquared) {
Expand Down Expand Up @@ -483,7 +550,15 @@ namespace quickhull {
}

// Create a tetrahedron half edge mesh and compute planes defined by each triangle
m_mesh.setup(baseTriangle[0],baseTriangle[1],baseTriangle[2],maxI);
// std::cout << "Active Point Base" << m_vertexData[baseTriangle[0]].x << " " << m_vertexData[baseTriangle[0]].y << " " << m_vertexData[baseTriangle[0]].z << std::endl;
// std::cout << "Active Point Base" << m_vertexData[baseTriangle[1]].x << " " << m_vertexData[baseTriangle[1]].y << " " << m_vertexData[baseTriangle[1]].z << std::endl;
// std::cout << "Active Point Base" << m_vertexData[baseTriangle[2]].x << " " << m_vertexData[baseTriangle[2]].y << " " << m_vertexData[baseTriangle[2]].z << std::endl;
// std::cout << "Active Point Base" << m_vertexData[maxI].x << " " << m_vertexData[maxI].y << " " << m_vertexData[maxI].z << std::endl;
tempPoints.push_back(baseTriangle[0]);
tempPoints.push_back(baseTriangle[1]);
tempPoints.push_back(baseTriangle[2]);
tempPoints.push_back(maxI);
m_mesh.setup(baseTriangle[0],baseTriangle[1],baseTriangle[2],maxI);
for (auto& f : m_mesh.m_faces) {
auto v = m_mesh.getVertexIndicesOfFace(f);
const Vector3<T>& va = m_vertexData[v[0]];
Expand All @@ -492,6 +567,10 @@ namespace quickhull {
const Vector3<T> N1 = mathutils::getTriangleNormal(va, vb, vc);
const Plane<T> plane(N1,va);
f.m_P = plane;
// std::cout << "Pt1" << va.x << " " << va.y << " " << va.z << std::endl;
// std::cout << "Pt2" << vb.x << " " << vb.y << " " << vb.z << std::endl;
// std::cout << "Pt3" << vc.x << " " << vc.y << " " << vc.z << std::endl;
// std::cout << "N" << N1.x << " " << N1.y << " " << N1.z << std::endl;
}

// Finally we assign a face for each vertex outside the tetrahedron (vertices inside the tetrahedron have no role anymore)
Expand Down
17 changes: 17 additions & 0 deletions QuickHull.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@ namespace quickhull {
std::vector< std::unique_ptr<std::vector<size_t>> > m_disabledFacePointVectors;
std::vector<size_t> m_visibleFaces;
std::vector<size_t> m_horizonEdges;
std::vector<size_t> tempPoints;
struct FaceData {
size_t m_faceIndex;
size_t m_enteredFromHalfEdge; // If the face turns out not to be visible,
Expand Down Expand Up @@ -125,6 +126,7 @@ namespace quickhull {
// The public getConvexHull functions will setup a VertexDataSource object and call this
ConvexHull<FloatType> getConvexHull(const VertexDataSource<FloatType>& pointCloud,
bool CCW, bool useOriginalIndices, FloatType eps);
bool isMeshConvex(const std::vector<size_t>& newFaces);
public:
// Computes convex hull for a given point cloud.
// Params:
Expand Down Expand Up @@ -194,11 +196,26 @@ namespace quickhull {
template<typename T>
bool QuickHull<T>::addPointToFace(typename MeshBuilder<T>::Face& f, size_t pointIndex) {
const T D = mathutils::getSignedDistanceToPlane(m_vertexData[ pointIndex ],f.m_P);
// if (D==0)
// {
// std::cout << "0point" << m_vertexData[pointIndex].x << " " << m_vertexData[pointIndex].y << " " << m_vertexData[pointIndex].z << std::endl;
// }
// if (D>m_epsilon && D*D > m_epsilonSquared*f.m_P.m_sqrNLength) {
if (D>0 && D*D > m_epsilonSquared*f.m_P.m_sqrNLength) {
// if (D>m_epsilon) {
// std::cout << D << " " << m_epsilonSquared*f.m_P.m_sqrNLength << std::endl;
// std::cout << "Point" << m_vertexData[pointIndex].x << " " << m_vertexData[pointIndex].y << " " << m_vertexData[pointIndex].z << std::endl;
if (!f.m_pointsOnPositiveSide) {
f.m_pointsOnPositiveSide = std::move(getIndexVectorFromPool());
}
f.m_pointsOnPositiveSide->push_back( pointIndex );
// if (D == f.m_mostDistantPointDist) {
// std::cout << "Check";
// std::cout << "Prior" << m_vertexData[f.m_mostDistantPoint].x << " " << m_vertexData[f.m_mostDistantPoint].y << " " << m_vertexData[f.m_mostDistantPoint].z << std::endl;
// std::cout << "New" << m_vertexData[pointIndex].x << " " << m_vertexData[pointIndex].y << " " << m_vertexData[pointIndex].z << std::endl;
// std::cout << "Normal" << f.m_P.m_N.x << " " << f.m_P.m_N.y << " " << f.m_P.m_N.z << std::endl;

// }
if (D > f.m_mostDistantPointDist) {
f.m_mostDistantPointDist = D;
f.m_mostDistantPoint = pointIndex;
Expand Down
2 changes: 1 addition & 1 deletion Structs/Plane.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ namespace quickhull {

// Normal length squared
T m_sqrNLength;

Vector3<T> errorNormal;
bool isPointOnPositiveSide(const Vector3<T>& Q) const {
T d = m_N.dotProduct(Q)+m_D;
if (d>=0) return true;
Expand Down