2
$\begingroup$

I implemented the $O(n^2)$ convex hull algorithm. That is:

  • Find a triangle known to be in the hull (by finding the lowest point, a point connected to it in the 2D convex hull, and the point that makes the widest angled triangle relative to all other triangles with thsoe 2 points).

  • Add the 3 edges and the triangle normal to q queue.

  • While the queue is not empty

    • pop and edge and the normal to its last seen triangle from the queue.
    • for each point, compute the dot between the last seen normal and the normal produced by the 2 edge points, pick the point with the largest dot (smallest angle).
    • for the 2 new edges, if they have not been seen before, push them and the new triangle's normal onto the queue.

This would be a rust implementation:

pub fn convex_hull(points: &Vec<Vec3>) -> Vec<usize>
{
    // The lowest point is part of the 3D convex hull.
    let lowest_id = points
        .iter()
        .enumerate()
        .min_by(|(_, p1), (_, p2)| p1.z.partial_cmp(&p2.z).unwrap())
        .map(|(id, _)| id)
        .unwrap();

    // The lowest point is guarateed to be connected to one of the points in the 2D convex
    // hull.
    let mut best_angle = 1.0;
    let mut best_angle_id = 0;
    for (id, p2) in points.iter().enumerate()
    {
        if id == lowest_id
        {
            continue;
        }

        let mut e1 = p2 - points[lowest_id];
        e1.y = 0.0;
        e1 = e1.normalize();

        let dir = Vec3::new(0.0, 0.0, 1.0);
        let dot = e1.dot(&dir);
        if dot < best_angle
        {
            best_angle = dot;
            best_angle_id = id;
        }
    }

    // Pick a candidate to initialize the next part of the algorithm. Should be distinct
    // from the points we already found.
    let mut candidate = 0;
    while candidate == lowest_id || candidate == best_angle_id
    {
        candidate += 1;
    }

    let get_normal = |i, j, k| -> Vec3 {
        return cross(&(points[j] - points[i]), &(points[k] - points[i])).normalize();
    };

    // Find the thrid point by finding a points such that all other points are in one side
    // of it.
    for i in 0..points.len()
    {
        if i == lowest_id || i == best_angle_id
        {
            continue;
        }

        let normal = get_normal(lowest_id, best_angle_id, candidate);
        if normal.dot(&(points[i] - points[lowest_id])) < 0.0
        {
            candidate = i;
        }
    }

    // Add the first triangle to the set.
    let mut face_ids = vec![lowest_id, best_angle_id, candidate];

    // Test and fix dingint order if necessary.
    let normal = get_normal(face_ids[0], face_ids[1], face_ids[2]).normalize();
    let test = normal.dot(&Vec3::new(0.0, 0.0, -1.0));
    if test < 0.0
    {
        face_ids.swap(0, 2);
    }

    let normal = get_normal(face_ids[0], face_ids[1], face_ids[2]).normalize();
    let mut edge_queue = vec![
        (face_ids[1], face_ids[0], normal),
        (face_ids[2], face_ids[1], normal),
        (face_ids[0], face_ids[2], normal),
    ];

    let mut seen_edges = HashSet::new();
    seen_edges.insert((face_ids[1], face_ids[0]));
    seen_edges.insert((face_ids[2], face_ids[1]));
    seen_edges.insert((face_ids[0], face_ids[2]));
    while !edge_queue.is_empty()
    {
        let (p1_id, p2_id, current_normal) = edge_queue.pop().unwrap();

        let mut best_point = 0;
        let mut best_dot = -1.0;
        let mut best_normal = Vec3::zeros();
        // Find the point that produces the most coplanar triangle with the other
        // triangle in the current edge.
        for point_id in 0..points.len()
        {
            let normal = get_normal(p1_id, p2_id, point_id);
            let dot = normal.dot(&current_normal);

            if dot > best_dot
            {
                best_point = point_id;
                best_dot = dot;
                best_normal = normal;
            }
        }

        face_ids.push(p1_id);
        face_ids.push(p2_id);
        face_ids.push(best_point);

        // Insert points in the opposite winding order to maintian orientation.
        let was_seen = !seen_edges.insert((best_point, p2_id));
        if !was_seen
        {
            edge_queue.push((best_point, p2_id, best_normal));
        }

        let was_seen = !seen_edges.insert((p1_id, best_point));
        if !was_seen
        {
            edge_queue.push((p1_id, best_point, best_normal));
        }
    }

    face_ids
}

However, although this works in theory, it greatly struggles with co-planarity. Let me show you. On a randomly distributed point cloud I get a fine CH:

enter image description here

But if I try it with a set where lots of points are coplanar I get something horrible:

enter image description here

enter image description here

The general shape of the CH is correct, but the topology I am getting out is just awful. How should I handle coplanarity?

$\endgroup$
3
  • $\begingroup$ I don't disagree, the question is however how to generate a properly topologized mesh in these cases. A bruteforcy way would be to remove and retriangulate but that is quite inelegant. $\endgroup$
    – Makogan
    Feb 19, 2023 at 21:54
  • $\begingroup$ Where there are four or more coplanar points on the hull, they form a face that is not a triangle but a convex polygon. I can't really follow your explanation, and I assume that your approach is gift wrapping. A possible solution is to accept all extreme vertices that are coplanar with the current edge and compute their convex hull (a planar problem). If your data structure only supports triangles, triangulate the face. $\endgroup$
    – user16034
    Feb 19, 2023 at 21:55
  • $\begingroup$ @YvesDaoust The current outpout is given by gift wrapping. On regions that are highly coplanar it generates the degenerate geomrtry I showed. I ned to modify the algorithm to not do that. "A possible solution is to accept all extreme vertices that are coplanar with the current edge and compute their convex hull (a planar problem)" This is potentially the solution. $\endgroup$
    – Makogan
    Feb 19, 2023 at 22:44

2 Answers 2

2
$\begingroup$

If You find a co-planar triangle candidate, it has to be one that does not contain any other points in that plane. I am not familiar with rust so here is a pseudo-code attempt to enhance your code:

if dot > best_dot || (
       dot == 0 &&
  best_dot == 0 && // <- should always be true at this point
  (triangle (p1_id, p2_id, best_point) contains point_id)
)
{
  best_point = point_id;
  best_dot = dot;
  best_normal = normal;
}

To test if a triangle $(a,b,c)$ contains a point $p$ on a plane with a normal vector $n$. You could think of $(a,b,c)$ as part of a triangular prism that is infinitely long in the direction of $n$. Or in other words check if $p$ is on the right side of the following three planes given by three points each:

  • (b,a,b+n)
  • (c,b,c+n)
  • (a,c,a+n)

Floating Point Rounding be a Harsh Mistress

While this would theoretically solve Your problem, there is another big headache in practice: Numeric accuracy. For this algorithm to work robustly, all Your dot products and co-planar and in-triangle testing have to be exact. You either could:

$\endgroup$
0
$\begingroup$

The geometric solution in case of coplanar points is to form the convex hull of the coplanar points in their plane, giving a convex polygon. The internal sites are discarded.

Unfortunately, I cannot tell you how to update your algorithm.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.