Convex hull of 4 points

2020-02-08 08:44发布

问题:

I would like an algorithm to calculate the convex hull of 4 2D points. I have looked at the algorithms for the generalized problem, but I wonder if there is a simple solution for 4 points.

回答1:

Take three of the points, and determine whether their triangle is clockwise or counterclockwise::

triangle_ABC= (A.y-B.y)*C.x + (B.x-A.x)*C.y + (A.x*B.y-B.x*A.y)

For a right-handed coordinate system, this value will be positive if ABC is counterclockwise, negative for clockwise, and zero if they are collinear. But, the following will work just as well for a left-handed coordinate system, as the orientation is relative.

Compute comparable values for three triangles containing the fourth point:

triangle_ABD= (A.y-B.y)*D.x + (B.x-A.x)*D.y + (A.x*B.y-B.x*A.y)
triangle_BCD= (B.y-C.y)*D.x + (C.x-B.x)*D.y + (B.x*C.y-C.x*B.y)
triangle_CAD= (C.y-A.y)*D.x + (A.x-C.x)*D.y + (C.x*A.y-A.x*C.y)

If all three of {ABD,BCD,CAD} have the same sign as ABC, then D is inside ABC, and the hull is triangle ABC.

If two of {ABD,BCD,CAD} have the same sign as ABC, and one has the opposite sign, then all four points are extremal, and the hull is quadrilateral ABCD.

If one of {ABD,BCD,CAD} has the same sign as ABC, and two have the opposite sign, then the convex hull is the triangle with the same sign; the remaining point is inside it.

If any of the triangle values are zero, the three points are collinear and the middle point is not extremal. If all four points are collinear, all four values should be zero, and the hull will be either a line or a point. Beware of numerical robustness problems in these cases!

For those cases where ABC is positive:

ABC  ABD  BCD  CAD  hull
------------------------
 +    +    +    +   ABC
 +    +    +    -   ABCD
 +    +    -    +   ABDC
 +    +    -    -   ABD
 +    -    +    +   ADBC
 +    -    +    -   BCD
 +    -    -    +   CAD
 +    -    -    -   [should not happen]


回答2:

Or just use Jarvis march.



回答3:

Here's a more ad-hoc algorithm specific to 4 points:

  • Find the indices of the points with minimum-X, maximum-X, minimum-Y and maximum-Y and get the unique values from this. For example, the indices may be 0,2,1,2 and the unique values will be 0,2,1.
  • If there are 4 unique values, then the convex hull is made up of all the 4 points.
  • If there are 3 unique values, then these 3 points are definitely in the convex hull. Check if the 4th point lies within this triangle; if not, it is also part of the convex hull.
  • If there are 2 unique values, then these 2 points are on the hull. Of the other 2 points, the point that is further away from this line joining these 2 points is definitely on the hull. Do a triangle containment test to check if the other point is also in the hull.
  • If there is 1 unique value, then all the 4 points are co-incident.

Some computation is required if there are 4 points to order them correctly so as to avoid getting a bow-tie shape. Hmmm.... Looks like there are enough special cases to justify using a generalized algorithm. However, you could possibly tune this to run faster than a generalized algorithm.



回答4:

I did a proof of concept fiddle based on a crude version of the gift wrapping algorithm.

Not efficient in the general case, but enough for only 4 points.

function Point (x, y)
{
    this.x = x;
    this.y = y;
}

Point.prototype.equals = function (p)
{
    return this.x == p.x && this.y == p.y;
};

Point.prototype.distance = function (p)
{ 
    return Math.sqrt (Math.pow (this.x-p.x, 2) 
                    + Math.pow (this.y-p.y, 2));
};

function convex_hull (points)
{
    function left_oriented (p1, p2, candidate)
    {
        var det = (p2.x - p1.x) * (candidate.y - p1.y) 
                - (candidate.x - p1.x) * (p2.y - p1.y);
        if (det > 0) return true;  // left-oriented 
        if (det < 0) return false; // right oriented
        // select the farthest point in case of colinearity
        return p1.distance (candidate) > p1.distance (p2);
    }

    var N = points.length;
    var hull = [];

    // get leftmost point
    var min = 0;
    for (var i = 1; i != N; i++)
    {
        if (points[i].y < points[min].y) min = i;
    }
    hull_point = points[min];

    // walk the hull
    do
    {
        hull.push(hull_point);

        var end_point = points[0];
        for (var i = 1; i != N; i++) 
        {
            if (  hull_point.equals (end_point)
               || left_oriented (hull_point, 
                                 end_point, 
                                 points[i]))
            {
                end_point = points[i];
            }
        }
        hull_point = end_point;
    }
    /*
     * must compare coordinates values (and not simply objects)
     * for the case of 4 co-incident points
     */
    while (!end_point.equals (hull[0])); 
    return hull;
}

It was fun :)



回答5:

I have written a fast implementation of comingstorm's answer using a lookup table. The case that all four points are colinear is not treated since my application does not need it. If the points are colinear the algorithm sets the first pointer point[0] to null. The hull contains of 3 points if point[3] is the null pointer, else the hull has 4 points. The hull is in counter-clockwise order for a coordinate system where the y-axis points to the top and the x-axis to the right.

const char hull4_table[] = {        
    1,2,3,0,1,2,3,0,1,2,4,3,1,2,3,0,1,2,3,0,1,2,4,0,1,2,3,4,1,2,4,0,1,2,4,0,
    1,2,3,0,1,2,3,0,1,4,3,0,1,2,3,0,0,0,0,0,0,0,0,0,2,3,4,0,0,0,0,0,0,0,0,0,
    1,4,2,3,1,4,3,0,1,4,3,0,2,3,4,0,0,0,0,0,0,0,0,0,2,3,4,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,2,4,3,0,0,0,0,0,0,0,0,0,1,2,4,0,1,3,4,0,1,2,4,0,1,2,4,0,
    0,0,0,0,0,0,0,0,1,4,3,0,0,0,0,0,0,0,0,0,0,0,0,0,1,3,4,0,0,0,0,0,0,0,0,0,
    1,4,2,0,1,4,2,0,1,4,3,0,1,4,2,0,0,0,0,0,0,0,0,0,2,3,4,0,0,0,0,0,0,0,0,0,
    0,0,0,0,0,0,0,0,2,4,3,0,0,0,0,0,0,0,0,0,2,4,3,0,1,3,4,0,1,3,4,0,1,3,2,4,
    0,0,0,0,0,0,0,0,2,4,3,0,0,0,0,0,0,0,0,0,1,3,2,0,1,3,4,0,1,3,2,0,1,3,2,0,
    1,4,2,0,1,4,2,0,1,4,3,2,1,4,2,0,1,3,2,0,1,3,2,0,1,3,4,2,1,3,2,0,1,3,2,0
};
struct Vec2i {
    int x, y;
};
typedef long long int64;    
inline int sign(int64 x) {
    return (x > 0) - (x < 0);
}
inline int64 orientation(const Vec2i& a, const Vec2i& b, const Vec2i& c) {
    return (int64)(b.x - a.x) * (c.y - b.y) - (b.y - a.y) * (c.x - b.x);
}

void convex_hull4(const Vec2i** points) {
    const Vec2i* p[5] = {(Vec2i*)0, points[0], points[1], points[2], points[3]};
    char abc = (char)1 - sign(orientation(*points[0], *points[1], *points[2]));
    char abd = (char)1 - sign(orientation(*points[0], *points[1], *points[3]));
    char cad = (char)1 - sign(orientation(*points[2], *points[0], *points[3]));
    char bcd = (char)1 - sign(orientation(*points[1], *points[2], *points[3]));

    const char* t = hull4_table + (int)4 * (bcd + 3*cad + 9*abd + 27*abc);
    points[0] = p[t[0]];
    points[1] = p[t[1]];
    points[2] = p[t[2]];
    points[3] = p[t[3]];
}


回答6:

Based on @comingstorm answer I created Swift solution:

func convexHull4(a: Pt, b: Pt, c: Pt, d: Pt) -> [LineSegment]? {

    let abc = (a.y-b.y)*c.x + (b.x-a.x)*c.y + (a.x*b.y-b.x*a.y)
    let abd = (a.y-b.y)*d.x + (b.x-a.x)*d.y + (a.x*b.y-b.x*a.y)
    let bcd = (b.y-c.y)*d.x + (c.x-b.x)*d.y + (b.x*c.y-c.x*b.y)
    let cad = (c.y-a.y)*d.x + (a.x-c.x)*d.y + (c.x*a.y-a.x*c.y)

    if  (abc > 0 && abd > 0 && bcd > 0 && cad > 0) ||
        (abc < 0 && abd < 0 && bcd < 0 && cad < 0) {
        //abc
        return [
            LineSegment(p1: a, p2: b),
            LineSegment(p1: b, p2: c),
            LineSegment(p1: c, p2: a)
        ]
    } else if   (abc > 0 && abd > 0 && bcd > 0 && cad < 0) ||
                (abc < 0 && abd < 0 && bcd < 0 && cad > 0) {
        //abcd
        return [
            LineSegment(p1: a, p2: b),
            LineSegment(p1: b, p2: c),
            LineSegment(p1: c, p2: d),
            LineSegment(p1: d, p2: a)
        ]
    } else if   (abc > 0 && abd > 0 && bcd < 0 && cad > 0) ||
                (abc < 0 && abd < 0 && bcd > 0 && cad < 0) {
        //abdc
        return [
            LineSegment(p1: a, p2: b),
            LineSegment(p1: b, p2: d),
            LineSegment(p1: d, p2: c),
            LineSegment(p1: c, p2: a)
        ]
    } else if   (abc > 0 && abd < 0 && bcd > 0 && cad > 0) ||
                (abc < 0 && abd > 0 && bcd < 0 && cad < 0) {
        //acbd
        return [
            LineSegment(p1: a, p2: c),
            LineSegment(p1: c, p2: b),
            LineSegment(p1: b, p2: d),
            LineSegment(p1: d, p2: a)
        ]
    } else if   (abc > 0 && abd > 0 && bcd < 0 && cad < 0) ||
                (abc < 0 && abd < 0 && bcd > 0 && cad > 0) {
        //abd
        return [
            LineSegment(p1: a, p2: b),
            LineSegment(p1: b, p2: d),
            LineSegment(p1: d, p2: a)
        ]
    } else if   (abc > 0 && abd < 0 && bcd > 0 && cad < 0) ||
                (abc < 0 && abd > 0 && bcd < 0 && cad > 0) {
        //bcd
        return [
            LineSegment(p1: b, p2: c),
            LineSegment(p1: c, p2: d),
            LineSegment(p1: d, p2: b)
        ]
    } else if   (abc > 0 && abd < 0 && bcd < 0 && cad > 0) ||
                (abc < 0 && abd > 0 && bcd > 0 && cad < 0) {
        //cad
        return [
            LineSegment(p1: c, p2: a),
            LineSegment(p1: a, p2: d),
            LineSegment(p1: d, p2: c)
        ]
    }

    return nil

}


回答7:

Based on comingstorm's solution I've created a C# solution which handles degenerate cases (e.g. 4 points form line or point).

https://gist.github.com/miyu/6e32e993d93d932c419f1f46020e23f0

  public static IntVector2[] ConvexHull3(IntVector2 a, IntVector2 b, IntVector2 c) {
     var abc = Clockness(a, b, c);
     if (abc == Clk.Neither) {
        var (s, t) = FindCollinearBounds(a, b, c);
        return s == t ? new[] { s } : new[] { s, t };
     }
     if (abc == Clk.Clockwise) {
        return new[] { c, b, a };
     }
     return new[] { a, b, c };
  }

  public static (IntVector2, IntVector2) FindCollinearBounds(IntVector2 a, IntVector2 b, IntVector2 c) {
     var ab = a.To(b).SquaredNorm2();
     var ac = a.To(c).SquaredNorm2();
     var bc = b.To(c).SquaredNorm2();
     if (ab > ac) {
        return ab > bc ? (a, b) : (b, c);
     } else {
        return ac > bc ? (a, c) : (b, c);
     }
  }

  // See https://stackoverflow.com/questions/2122305/convex-hull-of-4-points
  public static IntVector2[] ConvexHull4(IntVector2 a, IntVector2 b, IntVector2 c, IntVector2 d) {
     var abc = Clockness(a, b, c);

     if (abc == Clk.Neither) {
        var (s, t) = FindCollinearBounds(a, b, c);
        return ConvexHull3(s, t, d);
     }

     // make abc ccw
     if (abc == Clk.Clockwise) (a, c) = (c, a);

     var abd = Clockness(a, b, d);
     var bcd = Clockness(b, c, d);
     var cad = Clockness(c, a, d);

     if (abd == Clk.Neither) {
        var (s, t) = FindCollinearBounds(a, b, d);
        return ConvexHull3(s, t, c);
     }

     if (bcd == Clk.Neither) {
        var (s, t) = FindCollinearBounds(b, c, d);
        return ConvexHull3(s, t, a);
     }

     if (cad == Clk.Neither) {
        var (s, t) = FindCollinearBounds(c, a, d);
        return ConvexHull3(s, t, b);
     }

     if (abd == Clk.CounterClockwise) {
        if (bcd == Clk.CounterClockwise && cad == Clk.CounterClockwise) return new[] { a, b, c };
        if (bcd == Clk.CounterClockwise && cad == Clk.Clockwise) return new[] { a, b, c, d };
        if (bcd == Clk.Clockwise && cad == Clk.CounterClockwise) return new[] { a, b, d, c };
        if (bcd == Clk.Clockwise && cad == Clk.Clockwise) return new[] { a, b, d };
        throw new InvalidStateException();
     } else {
        if (bcd == Clk.CounterClockwise && cad == Clk.CounterClockwise) return new[] { a, d, b, c };
        if (bcd == Clk.CounterClockwise && cad == Clk.Clockwise) return new[] { d, b, c };
        if (bcd == Clk.Clockwise && cad == Clk.CounterClockwise) return new[] { a, d, c };
        // 4th state impossible
        throw new InvalidStateException();
     }
  }

You'll need to implement this boilerplate for your vector type:

  // relative to screen coordinates, so top left origin, x+ right, y+ down.
  // clockwise goes from origin to x+ to x+/y+ to y+ to origin, like clockwise if
  // you were to stare at a clock on your screen
  //
  // That is, if you draw an angle between 3 points on your screen, the clockness of that
  // direction is the clockness this would return.
  public enum Clockness {
     Clockwise = -1,
     Neither = 0,
     CounterClockwise = 1
  }
  public static Clockness Clockness(IntVector2 a, IntVector2 b, IntVector2 c) => Clockness(b - a, b - c);
  public static Clockness Clockness(IntVector2 ba, IntVector2 bc) => Clockness(ba.X, ba.Y, bc.X, bc.Y);
  public static Clockness Clockness(cInt ax, cInt ay, cInt bx, cInt by, cInt cx, cInt cy) => Clockness(bx - ax, by - ay, bx - cx, by - cy);
  public static Clockness Clockness(cInt bax, cInt bay, cInt bcx, cInt bcy) => (Clockness)Math.Sign(Cross(bax, bay, bcx, bcy));


回答8:

here is a complete analysis for the problem and efficient ruby code (minimizes the number of compares)

#   positions for d:
#
#           abc > 0                               abc < 0
#         (+-+- doesn't exist)               (-+-+ doesn't exist)
#
#
#                     |        /       ---+ \   --++    | -+++
#                     |       /        bdc   \  acbd    | acd
#                     | +-++ /                \         |
#                     | abd /         ---------A--------B---------
#                     |    /                    \  --+- |
#                     |   /                      \ acb  |
#                     |  /                        \     |
#                     | /                          \    |
#                     |/                  ----      \   | -++-
#                     C                   adcb       \  | acdb
#                    /|                               \ |
#                   / |                                \|
#         ++++     /  |                                 C
#         abcd    /   |                                 |\
#                /    | +--+                            | \
#               /     | abdc                            |  \
#              / ++-+ |                                 |   \
#             /  abc  |                                 |    \
#   ---------A--------B---------                        |     \
#     +++-  /         |                                 |      \
#     bcd  /   ++--   | +---                            | -+--  \
#         /    adbc   | adc                             | adb    \
#
#   or as table
#
#   ++++ abcd       -+++ acd
#   +++- bcd        -++- acdb
#   ++-+ abc        -+-+ XXXX
#   ++-- adbc       -+-- adb
#   +-++ abd        --++ acbd
#   +-+- XXXX       --+- acb
#   +--+ abdc       ---+ bdc
#   +--- adc        ---- adcb
#
# if there are some collinear points, the hull will be nil (for the moment)
#
def point_point_point_orientation(p, q, r)
  (q.x - p.x) * (r.y - q.y) - (q.y - p.y) * (r.x - q.x)
end



def convex_hull_4_points(a, b, c, d)
  abc = point_point_point_orientation(a, b, c)
  if abc.zero?
    # todo
    return nil
  end

  bcd = point_point_point_orientation(b, c, d)
  if bcd.zero?
    # todo
    return nil
  end

  cda = point_point_point_orientation(c, d, a)
  if cda.zero?
    # todo
    return nil
  end

  dab = point_point_point_orientation(d, a, b)
  if dab.zero?
    # todo
    return nil
  end

  if abc.positive?
    if bcd.positive?
      if cda.positive?
        if dab.positive?
          [a, b, c, d] # ++++
        else
          [b, c, d]    # +++-
        end
      else
        if dab.positive?
          [a, b, c]    # ++-+
        else
          [a, d, b, c] # ++--
        end
      end
    else
      if cda.positive?
        if dab.positive?
          [a, b, d]    # +-++
        else
          raise        # +-+-
        end
      else
        if dab.positive?
          [a, b, d, c] # +--+
        else
          [a, d, c]    # +---
        end
      end
    end
  else
    if bcd.positive?
      if cda.positive?
        if dab.positive?
          [a, c, d]    # -+++
        else
          [a, c, d, b] # -++-
        end
      else
        if dab.positive?
          raise        # -+-+
        else
          [a, d, b]    # -+--
        end
      end
    else
      if cda.positive?
        if dab.positive?
          [a, c, b, d] # --++
        else
          [a, c, b]    # --+-
        end
      else
        if dab.positive?
          [b, d, c]    # ---+
        else
          [a, d, c, b] # ----
        end
      end
    end
  end
end