Search code examples
pythonscipygeometryplane

Area of coplanar points in 3D - Python

I got a set of coplanar points in 3D of coordinates (x,y,z) as well as the equation of the plan P they belong to. I would like to know how to get the area occupied by those points on the pane using Python. I have tried Convhull in 3d, but I get an error because all my points are indeed coplanar. Do you have any idea ?

Thanks !

Solution

  • First I tried to reproduce your problem by creating some random, coplanar points and trying to fit a convex hull:

    import numpy as np
from scipy.spatial import ConvexHull

# 4 random 3D points
points = np.random.rand(4, 3)
# A set of coplanar points where their third dimension is 1.
coplanar = np.array([np.array([x[0], x[1], 1.0]) for x in points])

hull = ConvexHull(coplanar)

    This indeed generates an error:

    Traceback (most recent call last):
  File "main.py", line 9, in <module>
    hull = ConvexHull(coplanar)
  File "qhull.pyx", line 2428, in scipy.spatial.qhull.ConvexHull.__init__
  File "qhull.pyx", line 357, in scipy.spatial.qhull._Qhull.__init__
scipy.spatial.qhull.QhullError: QH6154 Qhull precision error: Initial simplex is flat (facet 1 is cop
lanar with the interior point)

..... OMITTING SOME LINES ....

If the input is lower dimensional:
  - use 'QJ' to joggle the input and make it full dimensional
  - use 'Qbk:0Bk:0' to delete coordinate k from the input.  You should
    pick the coordinate with the least range.  The hull will have the
    correct topology.
  - determine the flat containing the points, rotate the points
    into a coordinate plane, and delete the other coordinates.
  - add one or more points to make the input full dimensional.

    As we can see, the underlying library (qhull) gives us some suggestions in case your data is lower-dimensional. As you said, you already know that your data is coplanar and that it could be projected into a plane and be represented by 2D points.

    An alternative to projecting your data, as suggested in the error message, is to use the option QJ to joggle the input (option doc). If I understood it correctly, this "joggle" makes your data "non coplanar" by introducing a random noise to your data, allowing the optimization to proceed.

    hull = ConvexHull(coplanar, qhull_options="QJ")
hull.area
> 0.3100618849870332

    In summary:

    • find the angle between the plane and the unit xy plane, and rotate your data if you want a precise answer
    • use the option "QJ" if you want a quick solution and high precision is not a problem (I actually don't know how far away the answer can go from the truth)

    Update

    From Qhull FAQ:

    The area is the area of the surface of the convex hull, while the volume is the total volume of the convex hull.

    In 2-d, the convex hull is a polygon. Its surface is the edges of a polygon. So in 2-d, the 'area' is the length of the polygon's edges, while the 'volume' is the area of the polygon.

    From these statements, I interpreted that you should use:

    • hull.area / 2 for 3D points (your original points with the option QJ)
    • hull.volume for 2D points (if you rotate your points and get rid of one dimension)

    To help clarify the usage of ConvexHull, I used a simpler example

    square2D = np.array([
    [0.0, 0.0],
    [2.0, 0.0],
    [2.0, 2.0],
    [0.0, 2.0]
])
hull = ConvexHull(square2D)
hull.area
> 8.0
hull.volume
> 4.0

    The results are consistent with the documentation.

    Now, to understand the effects of the option QJ, I just added one dimension to the previous points:

    coplanar = np.array([
    [0.0, 0.0, 0.0],
    [2.0, 0.0, 0.0],
    [2.0, 2.0, 0.0],
    [0.0, 2.0, 0.0]
])
hull = ConvexHull(coplanar, qhull_options="QJ")
hull.area
> 8.000000000306578
hull.volume
> 5.719619827513867e-11

    For this last example, I imagine the convex hull as an approximately "flat" 3D object (because of the introduced noise) with two big facets. This is why I thought that you could use hull.area / 2 in this case.