scikit-image

1
import numpy as np
2
from scipy.special import elliprg
3

4

5
def ellipsoid(a, b, c, spacing=(1.0, 1.0, 1.0), levelset=False):
6
    """Generate ellipsoid for given semi-axis lengths.
7

8
    The respective semi-axis lengths are given along three dimensions in
9
    Cartesian coordinates. Each dimension may use a different grid spacing.
10

11
    Parameters
12
    ----------
13
    a : float
14
        Length of semi-axis along x-axis.
15
    b : float
16
        Length of semi-axis along y-axis.
17
    c : float
18
        Length of semi-axis along z-axis.
19
    spacing : 3-tuple of floats
20
        Grid spacing in three spatial dimensions.
21
    levelset : bool
22
        If True, returns the level set for this ellipsoid (signed level
23
        set about zero, with positive denoting interior) as np.float64.
24
        False returns a binarized version of said level set.
25

26
    Returns
27
    -------
28
    ellipsoid : (M, N, P) array
29
        Ellipsoid centered in a correctly sized array for given `spacing`.
30
        Boolean dtype unless `levelset=True`, in which case a float array is
31
        returned with the level set above 0.0 representing the ellipsoid.
32

33
    """
34
    if (a <= 0) or (b <= 0) or (c <= 0):
35
        raise ValueError('Parameters a, b, and c must all be > 0')
36

37
    offset = np.r_[1, 1, 1] * np.r_[spacing]
38

39
    # Calculate limits, and ensure output volume is odd & symmetric
40
    low = np.ceil(-np.r_[a, b, c] - offset)
41
    high = np.floor(np.r_[a, b, c] + offset + 1)
42

43
    for dim in range(3):
44
        if (high[dim] - low[dim]) % 2 == 0:
45
            low[dim] -= 1
46
        num = np.arange(low[dim], high[dim], spacing[dim])
47
        if 0 not in num:
48
            low[dim] -= np.max(num[num < 0])
49

50
    # Generate (anisotropic) spatial grid
51
    x, y, z = np.mgrid[
52
        low[0] : high[0] : spacing[0],
53
        low[1] : high[1] : spacing[1],
54
        low[2] : high[2] : spacing[2],
55
    ]
56

57
    if not levelset:
58
        arr = ((x / float(a)) ** 2 + (y / float(b)) ** 2 + (z / float(c)) ** 2) <= 1
59
    else:
60
        arr = ((x / float(a)) ** 2 + (y / float(b)) ** 2 + (z / float(c)) ** 2) - 1
61

62
    return arr
63

64

65
def ellipsoid_stats(a, b, c):
66
    """Calculate analytical volume and surface area of an ellipsoid.
67

68
    The surface area of an ellipsoid is given by
69

70
    .. math:: S=4\\pi b c R_G\\!\\left(1, \\frac{a^2}{b^2}, \\frac{a^2}{c^2}\\right)
71

72
    where :math:`R_G` is Carlson's completely symmetric elliptic integral of
73
    the second kind [1]_. The latter is implemented as
74
    :py:func:`scipy.special.elliprg`.
75

76
    Parameters
77
    ----------
78
    a : float
79
        Length of semi-axis along x-axis.
80
    b : float
81
        Length of semi-axis along y-axis.
82
    c : float
83
        Length of semi-axis along z-axis.
84

85
    Returns
86
    -------
87
    vol : float
88
        Calculated volume of ellipsoid.
89
    surf : float
90
        Calculated surface area of ellipsoid.
91

92
    References
93
    ----------
94
    .. [1] Paul Masson (2020). Surface Area of an Ellipsoid.
95
           https://analyticphysics.com/Mathematical%20Methods/Surface%20Area%20of%20an%20Ellipsoid.htm
96

97
    """
98
    if (a <= 0) or (b <= 0) or (c <= 0):
99
        raise ValueError('Parameters a, b, and c must all be > 0')
100

101
    # Volume
102
    vol = 4 / 3.0 * np.pi * a * b * c
103

104
    # Surface area
105
    surf = 3 * vol * elliprg(1 / a**2, 1 / b**2, 1 / c**2)
106

107
    return vol, surf
108