Creating a surface from a 3D data set and calculating the area in python

Question

I have a numpy array containing 70 points in 3D. These points are measurements of the edges of a flame with an ellipse-kind shape. Every row represent a different point, and the columns are the x,y,z coordinates. The array size is (70,3).

I want to fit/create an ellipse to this 3D data and calculate/approximate the area of that ellipse in Python.

3D Plot

The data array is:

data = np.array([[
    0.14893939, 0.14893939, 0.14621212, 0.14621212, 0.14651515, 0.14651515, 0.14984848, 0.14984848,
    0.15015152, 0.15015152, 0.15045455, 0.15045455, 0.15075758, 0.15075758, 0.14439394, 0.14439394,
    0.1519697, 0.1519697, 0.14378788, 0.14378788, 0.14348485, 0.14348485, 0.15318182, 0.15318182,
    0.15318182, 0.15318182, 0.15318182, 0.15318182, 0.15318182, 0.15318182, 0.15318182, 0.15318182,
    0.15318182, 0.15318182, 0.15318182, 0.15318182, 0.15287879, 0.15287879, 0.15287879, 0.15287879,
    0.14378788, 0.14378788, 0.14378788, 0.14378788, 0.14409091, 0.14409091, 0.14409091, 0.14409091,
    0.1519697, 0.1519697, 0.15166667, 0.15166667, 0.15106061, 0.15106061, 0.15136364, 0.15136364,
    0.14469697, 0.14469697, 0.15015152, 0.15015152, 0.15045455, 0.15045455, 0.14530303, 0.14530303,
    0.14530303, 0.14560606, 0.14560606, 0.14560606, 0.1480303, 0.1480303], [
    0.14560606, 0.14560606, 0.14590909, 0.14590909, 0.14590909, 0.14590909, 0.14590909, 0.14590909,
    0.14590909, 0.14590909, 0.14590909, 0.14590909, 0.14590909, 0.14590909, 0.14651515, 0.14651515,
    0.14681818, 0.14681818, 0.1480303, 0.1480303, 0.14924242, 0.14924242, 0.14924242, 0.14924242,
    0.15045455, 0.15045455, 0.15075758, 0.15075758, 0.15106061, 0.15106061, 0.15136364, 0.15136364,
    0.15166667, 0.15166667, 0.1519697, 0.1519697, 0.15287879, 0.15287879, 0.15318182, 0.15318182,
    0.15409091, 0.15409091, 0.15439394, 0.15439394, 0.15469697, 0.15469697, 0.155, 0.155, 0.155,
    0.155, 0.15530303, 0.15530303, 0.15560606, 0.15560606, 0.15560606, 0.15560606, 0.15621212,
    0.15621212, 0.15621212, 0.15621212, 0.15621212, 0.15621212, 0.15651515, 0.15651515, 0.15651515,
    0.15651515, 0.15651515, 0.15651515, 0.15651515, 0.15651515], [
    0.01907071, 0.01921212, 0.01935354, 0.01949495, 0.01935354, 0.01949495, 0.01907071, 0.01921212,
    0.01907071, 0.01921212, 0.01907071, 0.01921212, 0.01907071, 0.01921212, 0.02062626, 0.02076768,
    0.01878788, 0.01892929, 0.02034343, 0.02048485, 0.02006061, 0.02020202, 0.01793939, 0.01808081,
    0.01765657, 0.01779798, 0.01765657, 0.01779798, 0.01765657, 0.01779798, 0.01765657, 0.01779798,
    0.01765657, 0.01779798, 0.01765657, 0.01779798, 0.01765657, 0.01779798, 0.01765657, 0.01779798,
    0.01935354, 0.01949495, 0.01935354, 0.01949495, 0.01935354, 0.01949495, 0.01935354, 0.01949495,
    0.01793939, 0.01808081, 0.01793939, 0.01808081, 0.01793939, 0.01808081, 0.01793939, 0.01808081,
    0.02006061, 0.02020202, 0.01822222, 0.01836364, 0.01822222, 0.01836364, 0.01963636, 0.01977778,
    0.01991919, 0.01963636, 0.01977778, 0.01991919, 0.01822222, 0.01836364]]).T

I tried to fit triangles to this data and calculated the area but this method overestimated the area and created a really messy surface. Thus, I am looking for a better way to create a surface and calculate the area.

Answer 1

I believe to have finally understood what the OP is looking for: an ellipse (the parametric curve with one degree of freedom) fitted to a set of 3D points.

One way to do this is to project by PCA (principal components analysis) the points to a plane (the principal plane of the set of points), and then fit an ellipse in that 2D projection. We can then get the area of the ellipse (which is the same in 2D as in 3D, since the PCA is just a rotation), and optionally project backwards the ellipse into 3 dimensions.

The first step is to find the principal 2D plane that best spans the points. As described above, we use PCA to achieve that. We could use sklearn.decomposition.PCA, but since we already have numpy imported and it is quite straightforward, we roll our own:

import numpy as np

class PCA:
    def __init__(self, x, k=2):
        center = x.mean(axis=0)
        x = x - center
        u, v = np.linalg.eig(x.T @ x)
        ix = np.argsort(-u)  # u are non-negative; order them descending
        # keep k dimensions
        u = u[ix[:k]]
        v = v[:, ix[:k]]  # axes are columns
        self.center = center
        self.u = u
        self.v = v
        
    def proj(self, x):
        return (x - self.center) @ self.v
    
    def invproj(self, y):
        # pseudo-inverse of v is np.linalg.inv(v.T @ v) @ v.T, but
        # since columns of v are ortho-normal, we can just use v.T
        return y @ self.v.T + self.center

For the ellipse fit, we adapt a bit the excellent solution here. We further define conversions from implicit to parametric forms.

def ellipse_solve2D(p):
    # inspiration: https://stackoverflow.com/a/47881806/758174
    # 
    # Fit a 2D ellipse to an array of 2D points (p: (n, 2)).
    # This uses ordinary least squares to find best-fit coefficients
    # of the implicit equation:
    #
    #   a * x**2 + b * x * y + c * y**2 + d * x + e * y + f = 0, with f = -1
    x, y = p.T
    a = np.c_[x**2, x*y, y**2, x, y]
    b = np.ones_like(x)
    a, b, c, d, e = np.linalg.lstsq(a, b, rcond=None)[0]
    f = -1
    return a, b, c, d, e, f

def ellipse_convert_coeffs(A, B, C, D, E, F):
    # convert implicit ellipse equation:
    #
    #   a * x**2 + b * x * y + c * y**2 + d * x + e * y + f = 0
    #
    # into the canonical form:
    #
    #   x**2 / a**2 + y**2 / b**2 = 1
    #
    # with rotation by theta and shift center to (x0, y0)
    #
    # source: https://en.wikipedia.org/wiki/Ellipse#General_ellipse
    den = B**2 - 4*A*C
    z = np.sqrt((A - C)**2 + B**2)
    term = 2 * (A * E**2 + C * D**2 - B * D * E + den * F)
    a = -np.sqrt(term * (A + C + z)) / den
    b = -np.sqrt(term * (A + C - z)) / den
    x0 = (2 * C * D - B * E) / den
    y0 = (2 * A * E - B * D) / den
    theta = np.arctan((C - A - np.sqrt((A - C)**2 + B**2)) / B) \
        if B != 0 else 0 if A < C else np.pi / 2
    return a, b, x0, y0, theta

def ellipse_convert_vects(a, b, x0, y0, theta):
    # convert ellipse from canonical form to parametric:
    #
    # p = f0 + f1 cos(t) + f2 sin(t) | t: 0 .. 2 pi
    #
    # s.t. f1 and f2 are orthogonal (not generally required).
    f0 = np.array([x0, y0])
    ct, st = np.cos(theta), np.sin(theta)
    R = np.array([[ct, -st], [st, ct]])
    f1, f2 = np.array([[a, 0], [0, b]]) @ R
    return f0, f1, f2

Now, let's use these functions on the data provided in the question:

>>> data.shape
(70, 3)

pca = PCA(data, k=2)
p2d = pca.proj(data)
a, b, x0, y0, theta = ellipse_convert_coeffs(*ellipse_solve2D(p2d))
area = np.pi * a * b

>>> a, b
(0.005820995952691145, 0.005327126913750171)

>>> area
9.741822529039281e-05

Plot the 2D projection of the points, along with the ellipse found:

f0, f1, f2 = ellipse_convert_vects(a, b, x0, y0, theta)

# generate 2D points on the ellipse
t = np.linspace(0, 2*np.pi, 200)[:, None]
ellipse = f0 + f1 * np.cos(t) + f2 * np.sin(t)

fig, ax = plt.subplots()
ax.axline(f0, f0+f1, color='k', linewidth=1)
ax.axline(f0, f0+f2, color='k', linewidth=1)
ax.plot(*ellipse.T, 'r', linewidth=1)
ax.plot(*p2d.T, '.')
ax.set_title(f'area = {area:.3g}')
ax.set_aspect('equal')
plt.tight_layout()
plt.show()

Notes

1. It is interesting to observe that the fitted ellipse is not centered at (0,0) on the PCA projection: this is due to the fact that the original data is not uniformly distributed along the ellipse. Thus, the center of gravity of these points (corresponding to (0,0) in the PCA projection), is not necessarily at the center of the best fit ellipse. This shows that fitting explicitly the ellipse will yield a better fit than simply using the eigenvalues from the PCA.

2. The vectors f1 and f2 are two conjugate diameters of the ellipse. As noted on the Wikipedia page, they generally don't have to be orthogonal. But, in our case, we chose to make them orthogonal by construction.

Inverse projection (back to 3D)

Let us project back to 3D the ellipse found above, along with line segments representing the ellipse's principal axes f1 and f2:

import plotly.graph_objects as go

# 2D "axes" (using f1, f2, which from the above are orthogonal,
# with norm a and b, respectively)
axis_0 = f0 + np.c_[-f1, f1].T
axis_1 = f0 + np.c_[-f2, f2].T

# projection back to 3D
axis_0_ = pca.invproj(axis_0)
axis_1_ = pca.invproj(axis_1)
ellipse_ = pca.invproj(ellipse)

# 3D plot
fig = go.Figure([
    go.Scatter3d(
        **dict(zip('xyz', data.T)), mode='markers',
        marker=dict(size=5, symbol='circle', color=-data[:, 2], opacity=1),
    ),
    go.Scatter3d(
        **dict(zip('xyz', ellipse_.T)), mode='lines',
        line=dict(color='red', width=2),
    ),
    go.Scatter3d(
        **dict(zip('xyz', axis_0_.T)), mode='lines',
        line=dict(color='black', width=2),
    ),
    go.Scatter3d(
        **dict(zip('xyz', axis_1_.T)), mode='lines',
        line=dict(color='black', width=2),
    ),
])
fig.update_layout(showlegend=False)
fig.show()

If you run the code, the above is an interactive plotly 3D view that you can play with to inspect the fit in 3D.

Answer 2

You can use Delaunay tessellation method from Scipy to triangulate the surface. Not sure if this is what you tried that "overestimated the area and created a really messy surface" but let me know if you need a proof of concept code