Skip to content

Dim3

This module implements a 3D tessellation and trimming algorithm in the Tessellation3D class. This module is part of the tessellation package and can be used for 3D orbit tessellation tasks.

It inherits from TessellationBase and includes methods for calculating tetrahedron side lengths and tetrahedron volumes. The Normalization nested class includes normalization methods. Additionally, the class offers a plotting function to visualize the tessellation.

Tessellation3D

Bases: TessellationBase

A class for the tessellation and trimming algorithm applied in 3 dimensions.

Source code in commensurability/tessellation/dim3.py 
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
class Tessellation3D(TessellationBase):
    """
    A class for the tessellation and trimming algorithm applied in 3 dimensions.
    """

    @staticmethod
    def simplex_sides(*vertices: np.ndarray) -> list:
        """
        Compute the side lengths of a 3D simplex defined by its vertices.

        Args:
            *vertices: The vertices of the simplex.

        Returns:
            list: List of side lengths.

        """
        v1, v2, v3, v4 = vertices
        return [
            linalg.norm(v2 - v1),
            linalg.norm(v3 - v1),
            linalg.norm(v4 - v1),
            linalg.norm(v3 - v2),
            linalg.norm(v4 - v2),
            linalg.norm(v4 - v3),
        ]

    @staticmethod
    def simplex_measure(*vertices: np.ndarray) -> float:
        """
        Compute the measure (volume) of a 3D simplex defined by its vertices.

        Args:
            *vertices: The vertices of the simplex.

        Returns:
            float: The volume of the simplex.

        """
        (x1, y1, z1), (x2, y2, z2), (x3, y3, z3), (x4, y4, z4) = vertices
        a1 = (x2 - x1) * ((y3 - y1) * (z4 - z1) - (y4 - y1) * (z3 - z1))
        a2 = (x3 - x1) * ((y4 - y1) * (z2 - z1) - (y2 - y1) * (z4 - z1))
        a3 = (x4 - x1) * ((y2 - y1) * (z3 - z1) - (y3 - y1) * (z2 - z1))
        return abs(a1 + a2 + a3) / 6

    class Normalization:
        """
        A class providing various methods for normalization in 3D.

        Methods:
            sphere: Compute the volume of a sphere containing the points.
            cylinder: Compute the volume of a cylinder containing the points.
            Rz_convexhull: Compute the volume of the convex hull in the R vs z plane rotated around the z-axis.
            convexhull: Compute the volume of the convex hull of the points.
            convexhull_rot4: Compute the volume of the convex hull of the points including copied rotations by 90 degrees four times.
            default: Default normalization method (convexhull_rot4).
        """

        points: np.ndarray

        def sphere(self) -> float:
            """
            Compute the volume of a sphere containing the points.

            Returns:
                float: Volume of the sphere.
            """
            r = linalg.norm(self.points, axis=1)
            return 4 / 3 * np.pi * np.max(r) ** 3

        def cylinder(self) -> float:
            """
            Compute the volume of a cylinder containing the points.

            Returns:
                float: Volume of the cylinder.
            """
            x, y, z = self.points.T
            return np.pi * np.max(x**2 + y**2) * (np.max(z) - np.min(z))

        def Rz_convexhull(self) -> float:
            """
            Compute the volume of the convex hull in the R vs z plane rotated around the z-axis.

            Returns:
                float: Volume of the convex hull after rotation.
            """
            x, y, z = self.points.T
            R = np.sqrt(x**2 + y**2)
            points = np.array([R, z]).T
            points = np.array([*points, [0, max(z)], [0, min(z)]])
            hull = spatial.ConvexHull(points)

            poly = points[hull.vertices]
            start = poly[0]
            poly = poly[1:] - start

            # NOTE: requires "pairwise" from itertools (recent python version)
            from itertools import pairwise

            centroids = [start + (t1 + t2) / 3 for t1, t2 in pairwise(poly)]
            areas = [np.linalg.det([t1, t2]) for t1, t2 in pairwise(poly)]

            centroid = np.sum([a * c for c, a in zip(centroids, areas)], axis=0)
            return 2 * np.pi * float(np.linalg.norm(centroid[:2]))

        def convexhull(self) -> float:
            """
            Compute the volume of the convex hull of the points.

            Returns:
                float: Volume of the convex hull.

            Note:
                This method may not work correctly for co-rotation cases.
            """
            hull = spatial.ConvexHull(self.points)
            return hull.volume

        def convexhull_rot4(self) -> float:
            """
            Compute the volume of the convex hull of the points including copied rotations by 90 degrees four times.

            Returns:
                float: Volume of the convex hull.
            """
            x, y, z = self.points.T
            r000 = np.array([+x, +y, z]).T
            r090 = np.array([-y, +x, z]).T
            r180 = np.array([-x, -y, z]).T
            r270 = np.array([+y, -x, z]).T

            points = np.array([*r000, *r090, *r180, *r270])
            hull = spatial.ConvexHull(points)
            return hull.volume

        default = convexhull_rot4

    @property
    def volume(self) -> float:
        """
        Alias for `measure`.

        Returns:
            float: Volume of the tessellation (same as measure).
        """
        return self.measure

    def plot(self, ax, plot_included=True, plot_removed=False, plot_points=True):
        """
        Plot the 3D tessellation.
        Included tetrahedra are green, excluded tetrahedra are red.
        Shared edges are blue.

        Args:
            ax (mpl_toolkits.mplot3d.axes3d.Axes3D): Matplotlib 3D axes.
            plot_included (bool, optional): Whether to plot included triangles (default True).
            plot_removed (bool, optional): Whether to plot removed triangles (default False).
            plot_points (bool, optional): Whether to plot points (default True).

        Raises:
            RuntimeError: If tessellation failed.
        """
        if self.tri is None:
            raise RuntimeError("Tessellation failed; cannot produce tessellation plot")

        X, Y, Z = self.points.T

        included_edges = set()
        excluded_edges = set()
        for simplex, included in zip(self.tri.simplices, self.mask):
            i1, i2, i3, i4 = simplex
            edges = included_edges if included else excluded_edges
            for tuplet in (
                (i1, i2),
                (i1, i3),
                (i1, i4),
                (i2, i3),
                (i2, i4),
                (i3, i4),
            ):
                edges.add(tuplet)

        plotted_objects = set()

        if plot_removed and plot_included:
            inex_edges = included_edges & excluded_edges
            included_edges -= inex_edges
            excluded_edges -= inex_edges

            inex_lines = [(self.points[i1], self.points[i2]) for i1, i2 in inex_edges]

            line_collection = a3.art3d.Poly3DCollection(inex_lines)
            line_collection.set_edgecolor("blue")
            line_collection.set_linewidths(0.2)
            ax.add_collection3d(line_collection)
            plotted_objects.add(line_collection)

        if plot_removed:
            excluded_lines = [(self.points[i1], self.points[i2]) for i1, i2 in excluded_edges]

            line_collection = a3.art3d.Poly3DCollection(excluded_lines)
            line_collection.set_edgecolor("red")
            line_collection.set_linewidths(0.2)
            ax.add_collection3d(line_collection)
            plotted_objects.add(line_collection)

        if plot_included:
            included_lines = [(self.points[i1], self.points[i2]) for i1, i2 in included_edges]

            line_collection = a3.art3d.Poly3DCollection(included_lines)
            line_collection.set_edgecolor("green")
            line_collection.set_linewidths(0.2)
            ax.add_collection3d(line_collection)
            plotted_objects.add(line_collection)

        if plot_points:
            line = ax.scatter(X, Y, Z, marker=".", color="black")
            plotted_objects.add(line)

        ax_lim = 1.1 * max(max(X), max(Y), max(Z))
        ax_lim = (-ax_lim, ax_lim)
        ax.set(xlim=ax_lim, ylim=ax_lim, zlim=ax_lim)

        return plotted_objects

volume: float property

Alias for measure.

Returns:

Name Type Description
float float

Volume of the tessellation (same as measure).

Normalization

A class providing various methods for normalization in 3D.

Methods:

Name Description
sphere

Compute the volume of a sphere containing the points.
cylinder

Compute the volume of a cylinder containing the points.
Rz_convexhull

Compute the volume of the convex hull in the R vs z plane rotated around the z-axis.
convexhull

Compute the volume of the convex hull of the points.
convexhull_rot4

Compute the volume of the convex hull of the points including copied rotations by 90 degrees four times.
default

Default normalization method (convexhull_rot4).
Source code in commensurability/tessellation/dim3.py 
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
class Normalization:
    """
    A class providing various methods for normalization in 3D.

    Methods:
        sphere: Compute the volume of a sphere containing the points.
        cylinder: Compute the volume of a cylinder containing the points.
        Rz_convexhull: Compute the volume of the convex hull in the R vs z plane rotated around the z-axis.
        convexhull: Compute the volume of the convex hull of the points.
        convexhull_rot4: Compute the volume of the convex hull of the points including copied rotations by 90 degrees four times.
        default: Default normalization method (convexhull_rot4).
    """

    points: np.ndarray

    def sphere(self) -> float:
        """
        Compute the volume of a sphere containing the points.

        Returns:
            float: Volume of the sphere.
        """
        r = linalg.norm(self.points, axis=1)
        return 4 / 3 * np.pi * np.max(r) ** 3

    def cylinder(self) -> float:
        """
        Compute the volume of a cylinder containing the points.

        Returns:
            float: Volume of the cylinder.
        """
        x, y, z = self.points.T
        return np.pi * np.max(x**2 + y**2) * (np.max(z) - np.min(z))

    def Rz_convexhull(self) -> float:
        """
        Compute the volume of the convex hull in the R vs z plane rotated around the z-axis.

        Returns:
            float: Volume of the convex hull after rotation.
        """
        x, y, z = self.points.T
        R = np.sqrt(x**2 + y**2)
        points = np.array([R, z]).T
        points = np.array([*points, [0, max(z)], [0, min(z)]])
        hull = spatial.ConvexHull(points)

        poly = points[hull.vertices]
        start = poly[0]
        poly = poly[1:] - start

        # NOTE: requires "pairwise" from itertools (recent python version)
        from itertools import pairwise

        centroids = [start + (t1 + t2) / 3 for t1, t2 in pairwise(poly)]
        areas = [np.linalg.det([t1, t2]) for t1, t2 in pairwise(poly)]

        centroid = np.sum([a * c for c, a in zip(centroids, areas)], axis=0)
        return 2 * np.pi * float(np.linalg.norm(centroid[:2]))

    def convexhull(self) -> float:
        """
        Compute the volume of the convex hull of the points.

        Returns:
            float: Volume of the convex hull.

        Note:
            This method may not work correctly for co-rotation cases.
        """
        hull = spatial.ConvexHull(self.points)
        return hull.volume

    def convexhull_rot4(self) -> float:
        """
        Compute the volume of the convex hull of the points including copied rotations by 90 degrees four times.

        Returns:
            float: Volume of the convex hull.
        """
        x, y, z = self.points.T
        r000 = np.array([+x, +y, z]).T
        r090 = np.array([-y, +x, z]).T
        r180 = np.array([-x, -y, z]).T
        r270 = np.array([+y, -x, z]).T

        points = np.array([*r000, *r090, *r180, *r270])
        hull = spatial.ConvexHull(points)
        return hull.volume

    default = convexhull_rot4

Rz_convexhull()

Compute the volume of the convex hull in the R vs z plane rotated around the z-axis.

Returns:

Name Type Description
float float

Volume of the convex hull after rotation.
Source code in commensurability/tessellation/dim3.py 
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
def Rz_convexhull(self) -> float:
    """
    Compute the volume of the convex hull in the R vs z plane rotated around the z-axis.

    Returns:
        float: Volume of the convex hull after rotation.
    """
    x, y, z = self.points.T
    R = np.sqrt(x**2 + y**2)
    points = np.array([R, z]).T
    points = np.array([*points, [0, max(z)], [0, min(z)]])
    hull = spatial.ConvexHull(points)

    poly = points[hull.vertices]
    start = poly[0]
    poly = poly[1:] - start

    # NOTE: requires "pairwise" from itertools (recent python version)
    from itertools import pairwise

    centroids = [start + (t1 + t2) / 3 for t1, t2 in pairwise(poly)]
    areas = [np.linalg.det([t1, t2]) for t1, t2 in pairwise(poly)]

    centroid = np.sum([a * c for c, a in zip(centroids, areas)], axis=0)
    return 2 * np.pi * float(np.linalg.norm(centroid[:2]))

convexhull()

Compute the volume of the convex hull of the points.

Returns:

Name Type Description
float float

Volume of the convex hull.
Note

This method may not work correctly for co-rotation cases.

Source code in commensurability/tessellation/dim3.py 
124
125
126
127
128
129
130
131
132
133
134
135
def convexhull(self) -> float:
    """
    Compute the volume of the convex hull of the points.

    Returns:
        float: Volume of the convex hull.

    Note:
        This method may not work correctly for co-rotation cases.
    """
    hull = spatial.ConvexHull(self.points)
    return hull.volume

convexhull_rot4()

Compute the volume of the convex hull of the points including copied rotations by 90 degrees four times.

Returns:

Name Type Description
float float

Volume of the convex hull.
Source code in commensurability/tessellation/dim3.py 
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
def convexhull_rot4(self) -> float:
    """
    Compute the volume of the convex hull of the points including copied rotations by 90 degrees four times.

    Returns:
        float: Volume of the convex hull.
    """
    x, y, z = self.points.T
    r000 = np.array([+x, +y, z]).T
    r090 = np.array([-y, +x, z]).T
    r180 = np.array([-x, -y, z]).T
    r270 = np.array([+y, -x, z]).T

    points = np.array([*r000, *r090, *r180, *r270])
    hull = spatial.ConvexHull(points)
    return hull.volume

cylinder()

Compute the volume of a cylinder containing the points.

Returns:

Name Type Description
float float

Volume of the cylinder.
Source code in commensurability/tessellation/dim3.py 
88
89
90
91
92
93
94
95
96
def cylinder(self) -> float:
    """
    Compute the volume of a cylinder containing the points.

    Returns:
        float: Volume of the cylinder.
    """
    x, y, z = self.points.T
    return np.pi * np.max(x**2 + y**2) * (np.max(z) - np.min(z))

sphere()

Compute the volume of a sphere containing the points.

Returns:

Name Type Description
float float

Volume of the sphere.
Source code in commensurability/tessellation/dim3.py 
78
79
80
81
82
83
84
85
86
def sphere(self) -> float:
    """
    Compute the volume of a sphere containing the points.

    Returns:
        float: Volume of the sphere.
    """
    r = linalg.norm(self.points, axis=1)
    return 4 / 3 * np.pi * np.max(r) ** 3

plot(ax, plot_included=True, plot_removed=False, plot_points=True)

Plot the 3D tessellation. Included tetrahedra are green, excluded tetrahedra are red. Shared edges are blue.

Parameters:

Name Type Description Default
ax Axes3D

Matplotlib 3D axes.

 required
plot_included bool

Whether to plot included triangles (default True).

 True
plot_removed bool

Whether to plot removed triangles (default False).

 False
plot_points bool

Whether to plot points (default True).

 True

Raises:

Type Description
RuntimeError

If tessellation failed.
Source code in commensurability/tessellation/dim3.py 
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def plot(self, ax, plot_included=True, plot_removed=False, plot_points=True):
    """
    Plot the 3D tessellation.
    Included tetrahedra are green, excluded tetrahedra are red.
    Shared edges are blue.

    Args:
        ax (mpl_toolkits.mplot3d.axes3d.Axes3D): Matplotlib 3D axes.
        plot_included (bool, optional): Whether to plot included triangles (default True).
        plot_removed (bool, optional): Whether to plot removed triangles (default False).
        plot_points (bool, optional): Whether to plot points (default True).

    Raises:
        RuntimeError: If tessellation failed.
    """
    if self.tri is None:
        raise RuntimeError("Tessellation failed; cannot produce tessellation plot")

    X, Y, Z = self.points.T

    included_edges = set()
    excluded_edges = set()
    for simplex, included in zip(self.tri.simplices, self.mask):
        i1, i2, i3, i4 = simplex
        edges = included_edges if included else excluded_edges
        for tuplet in (
            (i1, i2),
            (i1, i3),
            (i1, i4),
            (i2, i3),
            (i2, i4),
            (i3, i4),
        ):
            edges.add(tuplet)

    plotted_objects = set()

    if plot_removed and plot_included:
        inex_edges = included_edges & excluded_edges
        included_edges -= inex_edges
        excluded_edges -= inex_edges

        inex_lines = [(self.points[i1], self.points[i2]) for i1, i2 in inex_edges]

        line_collection = a3.art3d.Poly3DCollection(inex_lines)
        line_collection.set_edgecolor("blue")
        line_collection.set_linewidths(0.2)
        ax.add_collection3d(line_collection)
        plotted_objects.add(line_collection)

    if plot_removed:
        excluded_lines = [(self.points[i1], self.points[i2]) for i1, i2 in excluded_edges]

        line_collection = a3.art3d.Poly3DCollection(excluded_lines)
        line_collection.set_edgecolor("red")
        line_collection.set_linewidths(0.2)
        ax.add_collection3d(line_collection)
        plotted_objects.add(line_collection)

    if plot_included:
        included_lines = [(self.points[i1], self.points[i2]) for i1, i2 in included_edges]

        line_collection = a3.art3d.Poly3DCollection(included_lines)
        line_collection.set_edgecolor("green")
        line_collection.set_linewidths(0.2)
        ax.add_collection3d(line_collection)
        plotted_objects.add(line_collection)

    if plot_points:
        line = ax.scatter(X, Y, Z, marker=".", color="black")
        plotted_objects.add(line)

    ax_lim = 1.1 * max(max(X), max(Y), max(Z))
    ax_lim = (-ax_lim, ax_lim)
    ax.set(xlim=ax_lim, ylim=ax_lim, zlim=ax_lim)

    return plotted_objects

simplex_measure(*vertices) staticmethod

Compute the measure (volume) of a 3D simplex defined by its vertices.

Parameters:

Name Type Description Default
*vertices ndarray

The vertices of the simplex.

 ()

Returns:

Name Type Description
float float

The volume of the simplex.
Source code in commensurability/tessellation/dim3.py 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
@staticmethod
def simplex_measure(*vertices: np.ndarray) -> float:
    """
    Compute the measure (volume) of a 3D simplex defined by its vertices.

    Args:
        *vertices: The vertices of the simplex.

    Returns:
        float: The volume of the simplex.

    """
    (x1, y1, z1), (x2, y2, z2), (x3, y3, z3), (x4, y4, z4) = vertices
    a1 = (x2 - x1) * ((y3 - y1) * (z4 - z1) - (y4 - y1) * (z3 - z1))
    a2 = (x3 - x1) * ((y4 - y1) * (z2 - z1) - (y2 - y1) * (z4 - z1))
    a3 = (x4 - x1) * ((y2 - y1) * (z3 - z1) - (y3 - y1) * (z2 - z1))
    return abs(a1 + a2 + a3) / 6

simplex_sides(*vertices) staticmethod

Compute the side lengths of a 3D simplex defined by its vertices.

Parameters:

Name Type Description Default
*vertices ndarray

The vertices of the simplex.

 ()

Returns:

Name Type Description
list list

List of side lengths.
Source code in commensurability/tessellation/dim3.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
@staticmethod
def simplex_sides(*vertices: np.ndarray) -> list:
    """
    Compute the side lengths of a 3D simplex defined by its vertices.

    Args:
        *vertices: The vertices of the simplex.

    Returns:
        list: List of side lengths.

    """
    v1, v2, v3, v4 = vertices
    return [
        linalg.norm(v2 - v1),
        linalg.norm(v3 - v1),
        linalg.norm(v4 - v1),
        linalg.norm(v3 - v2),
        linalg.norm(v4 - v2),
        linalg.norm(v4 - v3),
    ]