diff --git a/.github/workflows/python-package-min.yml b/.github/workflows/python-package-min.yml index 0920b9e..1404d59 100644 --- a/.github/workflows/python-package-min.yml +++ b/.github/workflows/python-package-min.yml @@ -17,7 +17,7 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.10", "3.11", "3.12", "3.13"] + python-version: ["3.11", "3.12", "3.13"] steps: - uses: actions/checkout@v2 diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index d84dede..110d48a 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -16,7 +16,7 @@ jobs: strategy: fail-fast: false matrix: - python-version: ["3.10", "3.11", "3.12", "3.13"] + python-version: ["3.11", "3.12", "3.13"] steps: - uses: actions/checkout@v2 diff --git a/doc/requirements.txt b/doc/requirements.txt index 50c4f28..42933f1 100644 --- a/doc/requirements.txt +++ b/doc/requirements.txt @@ -18,3 +18,4 @@ kaleido ipywidgets<8,>=7.6.0 sphinxcontrib-video anywidget # for plotly v6+ +scikit-image # for 3D implicit plotting diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 6af21fc..ed812af 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -2,6 +2,13 @@ Changelog ========== +v4.2.0 +====== + +Implemented 3D volumetric implicit plotting with K3D-Jupyter, Plotly +and Matplotlib. + + v4.1.2 ====== diff --git a/requirements.txt b/requirements.txt index 32d64e8..33b0f2c 100644 --- a/requirements.txt +++ b/requirements.txt @@ -25,5 +25,6 @@ psutil av tqdm kaleido +scikit-image # mayavi>=4.8.0 # PyQt5>=5.15.7 diff --git a/setup.py b/setup.py index d4cb2a0..d53c1f9 100755 --- a/setup.py +++ b/setup.py @@ -23,6 +23,7 @@ def readme(): "vtk", # needed for streamlines in k3d "control>=0.10.0", "imageio", "imageio-ffmpeg", "psutil", "av", "tqdm", # to save animations + "scikit-image", # for 3D implicit plotting # mayavi-related # "mayavi>=4.8.0", # "PyQt5>=5.15.7", diff --git a/spb/_version.py b/spb/_version.py index 13ffcf4..0fd7811 100644 --- a/spb/_version.py +++ b/spb/_version.py @@ -1 +1 @@ -__version__ = "4.1.2" +__version__ = "4.2.0" diff --git a/spb/backends/k3d/k3d.py b/spb/backends/k3d/k3d.py index 82509aa..dd7fd6d 100644 --- a/spb/backends/k3d/k3d.py +++ b/spb/backends/k3d/k3d.py @@ -6,12 +6,13 @@ from spb.backends.k3d.renderers import ( Line3DRenderer, Vector3DRenderer, ComplexRenderer, SurfaceRenderer, Implicit3DRenderer, + Implicit3DVoxelRenderer, GeometryRenderer, Arrow3DRenderer ) from spb.series import ( Parametric3DLineSeries, ComplexParametric3DLineSeries, List3DSeries, Vector3DSeries, SliceVector3DSeries, - RiemannSphereSeries, Implicit3DSeries, + RiemannSphereSeries, Implicit3DSeries, Implicit3DVoxelSeries, ComplexDomainColoringSeries, ComplexSurfaceSeries, SurfaceOver2DRangeSeries, ParametricSurfaceSeries, PlaneSeries, Geometry3DSeries, Arrow3DSeries @@ -65,6 +66,7 @@ class K3DBackend(Plot): Vector3DSeries: Vector3DRenderer, SliceVector3DSeries: Vector3DRenderer, Implicit3DSeries: Implicit3DRenderer, + Implicit3DVoxelSeries: Implicit3DVoxelRenderer, ComplexDomainColoringSeries: ComplexRenderer, ComplexSurfaceSeries: ComplexRenderer, RiemannSphereSeries: ComplexRenderer, diff --git a/spb/backends/k3d/renderers/__init__.py b/spb/backends/k3d/renderers/__init__.py index 99e9cf4..cd24cf9 100644 --- a/spb/backends/k3d/renderers/__init__.py +++ b/spb/backends/k3d/renderers/__init__.py @@ -1,6 +1,6 @@ from spb.backends.k3d.renderers.line3d import Line3DRenderer from spb.backends.k3d.renderers.surface import SurfaceRenderer -from spb.backends.k3d.renderers.implicit3d import Implicit3DRenderer +from spb.backends.k3d.renderers.implicit3d import Implicit3DRenderer, Implicit3DVoxelRenderer from spb.backends.k3d.renderers.vector3d import Vector3DRenderer from spb.backends.k3d.renderers.complex import ComplexRenderer from spb.backends.k3d.renderers.geometry import GeometryRenderer @@ -8,6 +8,7 @@ __all__ = [ "Line3DRenderer", "SurfaceRenderer", "Implicit3DRenderer", + "Implicit3DVoxelRenderer", "Vector3DRenderer", "ComplexRenderer", "GeometryRenderer", "Arrow3DRenderer" ] diff --git a/spb/backends/k3d/renderers/implicit3d.py b/spb/backends/k3d/renderers/implicit3d.py index d4beefc..66c2557 100644 --- a/spb/backends/k3d/renderers/implicit3d.py +++ b/spb/backends/k3d/renderers/implicit3d.py @@ -32,3 +32,43 @@ class Implicit3DRenderer(Renderer): draw_update_map = { _draw_implicit3d_helper: _update_implicit3d_helper } + + +def _draw_implicit3dvoxel_helper(renderer, data): + p, s = renderer.plot, renderer.series + np = p.np + + # NOTE: here I'm going to use voxels. K3D also exposes the volume function: + # https://k3d-jupyter.org/reference/factory.volume.html + # however, at the time of writing this, it doesn't scale the data to the + # appropriate boundaries, hence I'll only use voxels. + + mesh_x, mesh_y, mesh_z, f = data + xmin, xmax = s.start_x, s.end_x + ymin, ymax = s.start_y, s.end_y + zmin, zmax = s.start_z, s.end_z + color = p._convert_to_int(next(p._cl)) + + a = dict( + bounds=[xmin, xmax, ymin, ymax, zmin, zmax], + color_map=[color], + outlines=False + ) + kw = p.merge({}, a, s.rendering_kw) + plt_vox = p.k3d.voxels( + f.transpose(2,1,0).astype(np.uint8), + **kw + ) + p._fig += plt_vox + + return plt_vox + + +def _update_implicit3dvoxel_helper(renderer, data, handle): + raise NotImplementedError + + +class Implicit3DVoxelRenderer(Renderer): + draw_update_map = { + _draw_implicit3dvoxel_helper: _update_implicit3dvoxel_helper + } diff --git a/spb/backends/matplotlib/matplotlib.py b/spb/backends/matplotlib/matplotlib.py index ccd69e8..1819e62 100644 --- a/spb/backends/matplotlib/matplotlib.py +++ b/spb/backends/matplotlib/matplotlib.py @@ -10,7 +10,7 @@ NyquistRenderer, Arrow2DRendererFancyArrowPatch, Arrow3DRendererFancyArrowPatch, RootLocusRenderer, SGridLineRenderer, ZGridLineRenderer, NGridLineRenderer, MCirclesRenderer, PoleZeroRenderer, - NicholsLineRenderer + NicholsLineRenderer, Implicit3DRenderer, Implicit3DVoxelRenderer ) from spb.series import ( LineOver1DRangeSeries, List2DSeries, Parametric2DLineSeries, @@ -25,7 +25,8 @@ Arrow2DSeries, Arrow3DSeries, RootLocusSeries, SGridLineSeries, ZGridLineSeries, SystemResponseSeries, ColoredSystemResponseSeries, PoleZeroSeries, NGridLineSeries, PoleZeroWithSympySeries, - MCirclesSeries, HLineSeries, VLineSeries + MCirclesSeries, HLineSeries, VLineSeries, Implicit3DSeries, + Implicit3DVoxelSeries ) from sympy.external import import_module from packaging import version @@ -97,6 +98,8 @@ class MatplotlibBackend(Plot): MCirclesSeries: MCirclesRenderer, PoleZeroSeries: PoleZeroRenderer, PoleZeroWithSympySeries: PoleZeroRenderer, + Implicit3DSeries: Implicit3DRenderer, + Implicit3DVoxelSeries: Implicit3DVoxelRenderer } def __init__(self, *args, **kwargs): @@ -125,6 +128,7 @@ def __init__(self, *args, **kwargs): self.Line3DCollection = self.mpl_toolkits.mplot3d.art3d.Line3DCollection self.Path3DCollection = self.mpl_toolkits.mplot3d.art3d.Path3DCollection self.Axes3D = self.mpl_toolkits.mplot3d.Axes3D + self.Poly3DCollection = self.mpl_toolkits.mplot3d.art3d.Poly3DCollection kwargs["_library"] = "matplotlib" kwargs.setdefault("colormaps", [ diff --git a/spb/backends/matplotlib/renderers/__init__.py b/spb/backends/matplotlib/renderers/__init__.py index 46d808c..080d794 100644 --- a/spb/backends/matplotlib/renderers/__init__.py +++ b/spb/backends/matplotlib/renderers/__init__.py @@ -21,6 +21,7 @@ from spb.backends.matplotlib.renderers.mcircles import MCirclesRenderer from spb.backends.matplotlib.renderers.polezero import PoleZeroRenderer from spb.backends.matplotlib.renderers.nichols import NicholsLineRenderer +from spb.backends.matplotlib.renderers.implicit3d import Implicit3DRenderer, Implicit3DVoxelRenderer __all__ = [ "Line2DRenderer", "Line3DRenderer", "ContourRenderer", "SurfaceRenderer", @@ -30,5 +31,6 @@ "Arrow2DRendererQuivers", "Arrow2DRendererFancyArrowPatch", "Arrow3DRendererFancyArrowPatch", "RootLocusRenderer", "SGridLineRenderer", "ZGridLineRenderer", "NGridLineRenderer", "MCirclesRenderer", - "PoleZeroRenderer", "NicholsLineRenderer" + "PoleZeroRenderer", "NicholsLineRenderer", + "Implicit3DRenderer", "Implicit3DVoxelRenderer" ] diff --git a/spb/backends/matplotlib/renderers/implicit3d.py b/spb/backends/matplotlib/renderers/implicit3d.py new file mode 100644 index 0000000..f2e3f93 --- /dev/null +++ b/spb/backends/matplotlib/renderers/implicit3d.py @@ -0,0 +1,87 @@ +from spb.backends.matplotlib.renderers.renderer import MatplotlibRenderer +from spb.backends.utils import compute_streamtubes +from sympy.external import import_module + + +def _draw_implicit3d_isosurface_helper(renderer, data): + skimage = import_module("skimage", import_kwargs={'fromlist': ['measure']}) + if skimage is None: + raise ModuleNotFoundError( + "`skimage` is a mandatory dependency for 3D implicit plots," + " which is currently missing.") + + p, s = renderer.plot, renderer.series + mesh_x, mesh_y, mesh_z, f = data + verts, faces, normals, values = skimage.measure.marching_cubes(f, 0) + + x = mesh_x[:, 0, 0] + y = mesh_y[0, :, 0] + z = mesh_z[0, 0, :] + dx = abs(x[0] - x[1]) + dy = abs(y[0] - y[1]) + dz = abs(z[0] - z[1]) + verts[:, 0] = x[0] + verts[:, 0] * dx + verts[:, 1] = y[0] + verts[:, 1] * dy + verts[:, 2] = z[0] + verts[:, 2] * dz + + mkw = dict(linewidths=0.25, edgecolor="k", facecolors=next(p._cl)) + kw = p.merge({}, mkw, s.rendering_kw) + + # Fancy indexing: `verts[faces]` to generate a collection of triangles + mesh = p.Poly3DCollection(verts[faces], **kw) + p.ax.add_collection3d(mesh) + + return [mesh] + + +def _update_implicit3d_isosurface_helper(renderer, data, handle): + raise NotImplementedError + + +class Implicit3DRenderer(MatplotlibRenderer): + _sal = True + draw_update_map = { + _draw_implicit3d_isosurface_helper: _update_implicit3d_isosurface_helper + } + + +def _draw_implicit3d_volume_helper(renderer, data): + skimage = import_module("skimage", import_kwargs={'fromlist': ['measure']}) + if skimage is None: + raise ModuleNotFoundError( + "`skimage` is a mandatory dependency for 3D implicit plots," + " which is currently missing.") + + p, s = renderer.plot, renderer.series + mesh_x, mesh_y, mesh_z, f = data + x_edges = mesh_x[:, 0, 0] + y_edges = mesh_y[0, :, 0] + z_edges = mesh_z[0, 0, :] + dx = abs(x_edges[1] - x_edges[0]) + dy = abs(y_edges[1] - y_edges[0]) + dz = abs(z_edges[1] - z_edges[0]) + + x_min, x_max = float(s.range_x[1]), float(s.range_x[2]) + y_min, y_max = float(s.range_y[1]), float(s.range_y[2]) + z_min, z_max = float(s.range_z[1]), float(s.range_z[2]) + x_edges = p.np.linspace(x_min - dx / 2, x_max + dx / 2, s.n1) + y_edges = p.np.linspace(y_min - dy / 2, y_max + dy / 2, s.n2) + z_edges = p.np.linspace(z_min - dz / 2, z_max + dz / 2, s.n3) + x, y, z = p.np.meshgrid(x_edges, y_edges, z_edges, indexing="ij") + + mkw = dict(linewidths=0.25, edgecolor="k", facecolors=next(p._cl)) + kw = p.merge({}, mkw, s.rendering_kw) + mesh = p.ax.voxels(x, y, z, f, **kw) + + return [mesh] + + +def _update_implicit3d_volume_helper(renderer, data, handle): + raise NotImplementedError + + +class Implicit3DVoxelRenderer(MatplotlibRenderer): + _sal = True + draw_update_map = { + _draw_implicit3d_volume_helper: _update_implicit3d_volume_helper + } diff --git a/spb/backends/plotly/plotly.py b/spb/backends/plotly/plotly.py index a411d6c..04c6be9 100644 --- a/spb/backends/plotly/plotly.py +++ b/spb/backends/plotly/plotly.py @@ -8,14 +8,15 @@ from spb.backends.plotly.renderers import ( Line2DRenderer, Line3DRenderer, Vector2DRenderer, Vector3DRenderer, ComplexRenderer, ContourRenderer, SurfaceRenderer, Implicit3DRenderer, - GeometryRenderer, GenericRenderer, HVLineRenderer, Arrow2DRenderer + GeometryRenderer, GenericRenderer, HVLineRenderer, Arrow2DRenderer, + Implicit3DVoxelRenderer ) from spb.series import ( LineOver1DRangeSeries, List2DSeries, Parametric2DLineSeries, ColoredLineOver1DRangeSeries, AbsArgLineSeries, ComplexPointSeries, Parametric3DLineSeries, ComplexParametric3DLineSeries, List3DSeries, Vector2DSeries, Vector3DSeries, SliceVector3DSeries, - RiemannSphereSeries, Implicit3DSeries, + RiemannSphereSeries, Implicit3DSeries, Implicit3DVoxelSeries, ComplexDomainColoringSeries, ComplexSurfaceSeries, ContourSeries, SurfaceOver2DRangeSeries, ParametricSurfaceSeries, PlaneSeries, Geometry2DSeries, Geometry3DSeries, GenericDataSeries, @@ -70,6 +71,7 @@ class PlotlyBackend(Plot): Vector3DSeries: Vector3DRenderer, SliceVector3DSeries: Vector3DRenderer, Implicit3DSeries: Implicit3DRenderer, + Implicit3DVoxelSeries: Implicit3DVoxelRenderer, ComplexDomainColoringSeries: ComplexRenderer, ComplexSurfaceSeries: ComplexRenderer, RiemannSphereSeries: ComplexRenderer, diff --git a/spb/backends/plotly/renderers/__init__.py b/spb/backends/plotly/renderers/__init__.py index ddb8f08..4c9f57f 100644 --- a/spb/backends/plotly/renderers/__init__.py +++ b/spb/backends/plotly/renderers/__init__.py @@ -3,7 +3,7 @@ from spb.backends.plotly.renderers.contour import ContourRenderer from spb.backends.plotly.renderers.surface import SurfaceRenderer from spb.backends.plotly.renderers.geometry import GeometryRenderer -from spb.backends.plotly.renderers.implicit3d import Implicit3DRenderer +from spb.backends.plotly.renderers.implicit3d import Implicit3DRenderer, Implicit3DVoxelRenderer from spb.backends.plotly.renderers.vector2d import Vector2DRenderer from spb.backends.plotly.renderers.vector3d import Vector3DRenderer from spb.backends.plotly.renderers.complex import ComplexRenderer @@ -13,7 +13,7 @@ __all__ = [ "Line2DRenderer", "Line3DRenderer", "ContourRenderer", "SurfaceRenderer", - "GeometryRenderer", "Implicit3DRenderer", "Vector2DRenderer", - "Vector3DRenderer", "ComplexRenderer", "GenericRenderer", "HVLineRenderer", - "Arrow2DRenderer" + "GeometryRenderer", "Implicit3DRenderer", "Implicit3DVoxelRenderer", + "Vector2DRenderer", "Vector3DRenderer", "ComplexRenderer", + "GenericRenderer", "HVLineRenderer", "Arrow2DRenderer" ] diff --git a/spb/backends/plotly/renderers/implicit3d.py b/spb/backends/plotly/renderers/implicit3d.py index 461c908..db2f7f3 100644 --- a/spb/backends/plotly/renderers/implicit3d.py +++ b/spb/backends/plotly/renderers/implicit3d.py @@ -1,5 +1,6 @@ +from sympy.external import import_module from spb.backends.base_renderer import Renderer - +import numpy as np def _draw_implicit3d_helper(renderer, data): p, s = renderer.plot, renderer.series @@ -33,3 +34,232 @@ class Implicit3DRenderer(Renderer): draw_update_map = { _draw_implicit3d_helper: _update_implicit3d_helper } + + +def _draw_implicit3dvoxel_helper(renderer, data): + p, s = renderer.plot, renderer.series + np = p.np + + # NOTE: plotly exposes the volume plots: + # https://plotly.com/python/3d-volume-plots/ + # However, it is very difficult to use them in this context, where people + # can write (f(x) > 0) | (g(x) < 2). Hence, I'll just use voxels. + + x_edges = np.linspace(float(s.range_x[1]), float(s.range_x[2]), s.n1) + y_edges = np.linspace(float(s.range_y[1]), float(s.range_y[2]), s.n2) + z_edges = np.linspace(float(s.range_z[1]), float(s.range_z[2]), s.n3) + mesh_x, mesh_y, mesh_z, f = data + + col = next(p._cl) + colorscale = [[0, col], [1, col]] + + mesh = voxel_mesh( + f, x_edges, y_edges, z_edges, color=col, **s.rendering_kw) + p._fig.add_trace(mesh) + p._colorbar_counter += 1 + return len(p._fig.data) - 1 + + + +def _update_implicit3dvoxel_helper(renderer, data, handle): + raise NotImplementedError + + +class Implicit3DVoxelRenderer(Renderer): + draw_update_map = { + _draw_implicit3dvoxel_helper: _update_implicit3dvoxel_helper + } + + +def voxel_mesh( + mask, x_edges, y_edges, z_edges, color="steelblue", opacity=1.0, + flatshading=True, cull_hidden_faces=True, dedupe_vertices=True, + **mesh_kwargs +): + """ + Build a single go.Mesh3d trace rendering a cube for every + True entry in `mask`. + + Disclaimer: this function was built with Claude Sonnet 4.6 and ChatGPT 5.5. + + Parameters + ---------- + mask : ndarray of bool, shape (nx, ny, nz) + Mask[i, j, k] = True means draw a cube occupying the cell between + x_edges[i]:x_edges[i+1], y_edges[j]:y_edges[j+1], z_edges[k]:z_edges[k+1]. + x_edges, y_edges, z_edges : ndarray + 1D arrays of cell boundaries (length nx+1, ny+1, nz+1). + color : str + Single color for all voxels (string, e.g. "steelblue" or "#1f77b4"). + opacity : float + flatshading : bool + Keeps each cube face looking flat/blocky rather than smooth-shaded. + On its own this is not always sufficient -- see the `lighting` note below. + cull_hidden_faces : bool, default True. + Skips faces shared by two filled voxels (interior faces nobody will + ever see). This is the standard voxel-engine trick and massively + reduces triangle count for solid blobs. + dedupe_vertices : bool, default True. + Merges duplicate (x, y, z) points so vertices shared between adjacent + faces are stored once instead of once per face. Doesn't change triangle + count or appearance, but typically cuts the serialized figure size + (and browser memory) by an order of magnitude or more at high + resolution, since Mesh3d's vertex arrays are the dominant cost. + Implemented as a single vectorized np.unique pass, + and -- importantly -- vertices are only ever built for faces that + survive culling in the first place, so the array being deduplicated + is already small; this avoids a return-to-Python-loop-speed dedup + pass on a much larger, mostly-discarded array. + + One side effect: with dedupe on, many faces now share exact vertices + at voxel edges (by construction -- that's the point of deduping). + Plotly's WebGL renderer computes per-vertex normals by averaging the + normals of every triangle touching a vertex, then has a separate + epsilon-based merge step for near-identical normals before flat + shading is applied. At a voxel edge, two perpendicular faces meet at + a shared vertex whose averaged normal points diagonally between + them, and that merge step can still blend lighting across the edge + even with flatshading=True, showing up as a faint gradient/halo right + on the seams instead of a crisp line. Passing + lighting=dict(vertexnormalsepsilon=0, facenormalsepsilon=0) disables + that merge step entirely, which is why it's the default below -- it + forces every triangle to use its own exact face normal with no + cross-triangle blending. Override by passing your own `lighting=...` + in mesh_kwargs if you want Plotly's default behavior back. + **mesh_kwargs : passed through to go.Mesh3d (e.g. lighting, etc.) + + Returns + ------- + go.Mesh3d trace. + """ + np = import_module("numpy") + plotly = import_module( + 'plotly', + import_kwargs={'fromlist': ['graph_objects', 'figure_factory']}, + warn_not_installed=True, + min_module_version='5.0.0') + go = plotly.graph_objects + + # 8 corners of a unit cube, in a fixed order, indexed by the face definitions below + _UNIT_CORNERS = np.array([ + [0, 0, 0], + [1, 0, 0], + [1, 1, 0], + [0, 1, 0], + [0, 0, 1], + [1, 0, 1], + [1, 1, 1], + [0, 1, 1], + ], dtype=float) + + # Outward-facing quads for each of the 6 faces (4 corner indices into + # _UNIT_CORNERS, in winding order), plus which neighbor cell (delta i,j,k) to + # check when deciding whether that face is hidden. + # Corner numbering: 0=(0,0,0) 1=(1,0,0) 2=(1,1,0) 3=(0,1,0) 4=(0,0,1) 5=(1,0,1) 6=(1,1,1) 7=(0,1,1) + _FACES = [ + ((-1, 0, 0), [0, 4, 7, 3]), # x- face + ((1, 0, 0), [1, 2, 6, 5]), # x+ face + ((0, -1, 0), [0, 1, 5, 4]), # y- face + ((0, 1, 0), [3, 7, 6, 2]), # y+ face + ((0, 0, -1), [0, 3, 2, 1]), # z- face + ((0, 0, 1), [4, 5, 6, 7]), # z+ face + ] + # Two triangles covering a quad's 4 local corners (0,1,2,3), fan-split from corner 0 + _QUAD_TRIS_LOCAL = np.array([[0, 1, 2], [0, 2, 3]], dtype=int) + + + idx = np.argwhere(mask) + n = len(idx) + if n == 0: + return go.Mesh3d(x=[], y=[], z=[], i=[], j=[], k=[], color=color, opacity=opacity) + + nx, ny, nz = mask.shape + + if not cull_hidden_faces: + face_keep = np.ones((n, 6), dtype=bool) + else: + face_keep = np.empty((n, 6), dtype=bool) + for f, (delta, _) in enumerate(_FACES): + ni = idx[:, 0] + delta[0] + nj = idx[:, 1] + delta[1] + nk = idx[:, 2] + delta[2] + in_bounds = ( + (ni >= 0) & (ni < nx) & + (nj >= 0) & (nj < ny) & + (nk >= 0) & (nk < nz) + ) + neighbor_filled = np.zeros(n, dtype=bool) + if in_bounds.any(): + neighbor_filled[in_bounds] = mask[ni[in_bounds], nj[in_bounds], nk[in_bounds]] + # keep this face only if there's no filled neighbor on that side + face_keep[:, f] = ~neighbor_filled + + # Lower/upper bounds of each occupied cell + x0 = x_edges[idx[:, 0]] + x1 = x_edges[idx[:, 0] + 1] + y0 = y_edges[idx[:, 1]] + y1 = y_edges[idx[:, 1] + 1] + z0 = z_edges[idx[:, 2]] + z1 = z_edges[idx[:, 2] + 1] + + lo = np.stack([x0, y0, z0], axis=1) # (n, 3) + span = np.stack([x1 - x0, y1 - y0, z1 - z0], axis=1) # (n, 3) + + # Build vertices/triangles ONLY for faces that survived culling. This is + # the key difference from a naive "build all 8 corners of every cube, + # then dedupe": a fully-interior cube touching filled neighbors on all 6 + # sides contributes 0 rows here instead of 8 rows that np.unique would + # later have to sort through and discard. Each kept face gets exactly 4 + # corner vertices (still duplicated across adjacent faces at this stage -- + # the dedupe step below merges those). + vert_chunks = [] + tri_chunks = [] + running_offset = 0 + for f, (_, quad_corners) in enumerate(_FACES): + keep = face_keep[:, f] + if not keep.any(): + continue + corners = _UNIT_CORNERS[quad_corners] # (4, 3) + # (m, 4, 3): one quad of 4 vertices per kept face + quad_verts = lo[keep][:, None, :] + span[keep][:, None, :] * corners[None, :, :] + m = quad_verts.shape[0] + vert_chunks.append(quad_verts.reshape(-1, 3)) + + offs = (np.arange(m) * 4 + running_offset)[:, None, None] # (m, 1, 1) + tris = _QUAD_TRIS_LOCAL[None, :, :] + offs # (m, 2, 3) + tri_chunks.append(tris.reshape(-1, 3)) + running_offset += m * 4 + + if vert_chunks: + verts = np.concatenate(vert_chunks, axis=0) + all_tris = np.concatenate(tri_chunks, axis=0) + else: + verts = np.empty((0, 3)) + all_tris = np.empty((0, 3), dtype=int) + + if dedupe_vertices and len(verts) > 0: + # Merge identical (x, y, z) rows in one vectorized pass. inverse[i] is + # the index into the deduped array for original vertex i, so remapping + # the triangle indices through `inverse` repoints every triangle at + # the shared vertex without ever touching vertices in a Python loop. + # Because we only built vertices for kept faces above, every row here + # is referenced by at least one triangle -- no separate "drop unused + # vertices" pass is needed. + verts, inverse = np.unique(verts, axis=0, return_inverse=True) + inverse = inverse.reshape(-1) # numpy >=2.0 returns a column vector here + all_tris = inverse[all_tris] + + # Disable Plotly's epsilon-based normal-merging so flatshading actually + # produces crisp edges at shared voxel-face vertices (see docstring). + # Respect an explicit lighting= the caller passed in mesh_kwargs instead + # of silently overriding it. + mesh_kwargs.setdefault("lighting", dict(vertexnormalsepsilon=0, facenormalsepsilon=0)) + + return go.Mesh3d( + x=verts[:, 0], y=verts[:, 1], z=verts[:, 2], + i=all_tris[:, 0], j=all_tris[:, 1], k=all_tris[:, 2], + color=color, + opacity=opacity, + flatshading=flatshading, + **mesh_kwargs, + ) diff --git a/spb/graphics/functions_3d.py b/spb/graphics/functions_3d.py index 39cd3f7..60b7611 100644 --- a/spb/graphics/functions_3d.py +++ b/spb/graphics/functions_3d.py @@ -1,11 +1,15 @@ from sympy import ( - pi, Symbol, sin, cos, sqrt, atan2, Tuple, Plane + pi, Symbol, sin, cos, sqrt, atan2, Tuple, Plane, Rel, + StrictGreaterThan, GreaterThan, StrictLessThan, LessThan, + Equality ) +from sympy.logic.boolalg import BooleanFunction from spb.doc_utils.docstrings import _PARAMS from spb.doc_utils.ipython import modify_graphics_series_doc from spb.series import ( Parametric3DLineSeries, SurfaceOver2DRangeSeries, ParametricSurfaceSeries, - Implicit3DSeries, List3DSeries, ComplexSurfaceBaseSeries, PlaneSeries + Implicit3DSeries, List3DSeries, ComplexSurfaceBaseSeries, PlaneSeries, + Implicit3DVoxelSeries ) from spb.utils import ( _create_missing_ranges, _preprocess_multiple_ranges, @@ -780,7 +784,7 @@ def implicit_3d( rendering_kw=None, **kwargs ): """ - Plots an isosurface of a function. + Plots an isosurface or a volumetric plot of a function. Notes ===== @@ -792,6 +796,7 @@ def implicit_3d( 3. To plot ``f(x, y, z) = c`` either write ``expr = f(x, y, z) - c`` or pass the appropriate keyword to ``rendering_kw``. Read the backends documentation to find out the available options. + 4. Volumetric plots will use voxels. Returns ======= @@ -819,6 +824,19 @@ def implicit_3d( implicit_3d(x**4 + y**4 + z**4 - (x**2 + y**2 + z**2 - 0.3), (x, -2, 2), (y, -2, 2), (z, -2, 2)), backend=PB) + + .. k3d-screenshot:: + :camera: 4.3, -5.82, 4.95, 0.4, -0.25, -0.67, -0.32, 0.5, 0.8 + + from sympy import * + from spb import * + var("x, y, z") + graphics( + implicit_3d( + Eq(sin(x) + cos(y) + sin(y) * cos(z) + sin(z) * cos(x), 0), + (x, -2*pi, 2*pi), (y, -2*pi, 2*pi), (z, -2*pi, 2*pi), n=100), + backend=KB + ) Visualize the isocontours from `isomin=0` to `isomax=2` by providing a ``rendering_kw`` dictionary: @@ -836,6 +854,27 @@ def implicit_3d( backend=PB ) + Volumetric plots: + + .. k3d-screenshot:: + :camera: 10.29, -15.33, -9.3, 0, 0, 0, 0.28, -0.41, 0.87 + + graphics( + implicit_3d(x*y*z < 1, (x, -5, 5), (y, -5, 5), (z, -5, 5), n=150), + backend=KB + ) + + .. k3d-screenshot:: + :camera: -5.2, -1.86, 2.06, 0, 0.87, 0, 0.23, 0.16, 0.96 + + graphics( + implicit_3d( + (x**2 + y**2 + z**2 >= 1) & (x**2 + y**2 + z**2 <= 3) & (y >= 0), + (x, -2, 2), (y, -2, 2), (z, -2, 2), n=150), + backend=KB + ) + + See Also ======== @@ -861,10 +900,15 @@ def implicit_3d( "visualization might be flipped." ) + if isinstance(expr, Equality): + expr = expr.lhs - expr.rhs + ranges = _preprocess_multiple_ranges( [expr], [range_x, range_y, range_z], 3, params) - s = Implicit3DSeries( - expr, *ranges, label, rendering_kw=rendering_kw, **kwargs) + types = (BooleanFunction, StrictGreaterThan, GreaterThan, + StrictLessThan, LessThan) + cl = Implicit3DVoxelSeries if isinstance(expr, types) else Implicit3DSeries + s = cl(expr, *ranges, label, rendering_kw=rendering_kw, **kwargs) return [s] diff --git a/spb/plot_functions/functions_3d.py b/spb/plot_functions/functions_3d.py index d79750a..802170a 100644 --- a/spb/plot_functions/functions_3d.py +++ b/spb/plot_functions/functions_3d.py @@ -782,6 +782,25 @@ def plot3d_implicit(*args, **kwargs): (expr2, range_x2, range_y2, range_z2, rendering_kw2 [opt]), **kwargs)` + Volumetric plots: + + .. k3d-screenshot:: + :camera: 10.29, -15.33, -9.3, 0, 0, 0, 0.28, -0.41, 0.87 + + plot3d_implicit( + x*y*z < 1, (x, -5, 5), (y, -5, 5), (z, -5, 5), + n=150, backend=KB + ) + + .. k3d-screenshot:: + :camera: -5.2, -1.86, 2.06, 0, 0.87, 0, 0.23, 0.16, 0.96 + + plot3d_implicit( + (x**2 + y**2 + z**2 >= 1) & (x**2 + y**2 + z**2 <= 3) & (y >= 0), + (x, -2, 2), (y, -2, 2), (z, -2, 2), + n=150, backend=KB + ) + Notes ===== 1. the number of discretization points is crucial as the algorithm will @@ -792,6 +811,7 @@ def plot3d_implicit(*args, **kwargs): 3. To plot ``f(x, y, z) = c`` either write ``expr = f(x, y, z) - c`` or pass the appropriate keyword to ``rendering_kw``. Read the backends documentation to find out the available options. + 4. Volumetric plots will use voxels. Examples ======== diff --git a/spb/series/__init__.py b/spb/series/__init__.py index 20659c1..b721f4c 100644 --- a/spb/series/__init__.py +++ b/spb/series/__init__.py @@ -43,6 +43,7 @@ ParametricSurfaceSeries, ImplicitSeries, Implicit3DSeries, + Implicit3DVoxelSeries, PlaneSeries, Geometry2DSeries, Geometry3DSeries, @@ -90,6 +91,7 @@ ParametricSurfaceSeries, ImplicitSeries, Implicit3DSeries, + Implicit3DVoxelSeries, PlaneSeries, Geometry2DSeries, Geometry3DSeries, diff --git a/spb/series/evaluator.py b/spb/series/evaluator.py index 1def886..61924a8 100644 --- a/spb/series/evaluator.py +++ b/spb/series/evaluator.py @@ -572,6 +572,14 @@ class GridEvaluator(_LambdifierMixin): numpy functions don't like complex-type arguments, and appropriate actions should be taken.""" ) + at_cell_center = param.Boolean(False, doc=""" + Consider a range that should be discretized with N number of + discretization points, for example from -2 to 2 with N=5. + + * if ``at_cell_center=False``: will produce [-2, -1, 0, 1, 2, 3] + * if ``at_cell_center=True``: will produce [-1.5, -0.5, 0.5, 1.5]. + Note that the number of points is N-1. This is particularly useful + when dealing with volumetric discretization, such in voxels.""") def _create_discretized_domain(self): """ @@ -598,6 +606,9 @@ def _create_discretized_domain(self): only_integers=needs_integer_discr ) + if self.at_cell_center: + d = (d[:-1] + d[1:]) / 2 + if ( (not self.series.force_real_eval) and (not needs_integer_discr) and diff --git a/spb/series/series_2d_3d.py b/spb/series/series_2d_3d.py index 529fdb3..0647dee 100644 --- a/spb/series/series_2d_3d.py +++ b/spb/series/series_2d_3d.py @@ -2223,7 +2223,8 @@ def __init__(self, expr, range_x, range_y, range_z, label="", **kwargs): kwargs["range_y"] = range_y kwargs["range_z"] = range_z kwargs["_range_names"] = ["range_x", "range_y", "range_z"] - kwargs.setdefault("evaluator", GridEvaluator(series=self)) + kwargs.setdefault("evaluator", GridEvaluator( + series=self, at_cell_center=False)) super().__init__(**kwargs) self.evaluator.set_expressions() self.var_x, self.start_x, self.end_x = self.range_x @@ -2275,6 +2276,14 @@ def _apply_transform(self, *args): return t(x, self.tx), t(y, self.ty), t(z, self.tz), f +@modify_parameterized_doc() +class Implicit3DVoxelSeries(Implicit3DSeries): + def __init__(self, expr, range_x, range_y, range_z, label="", **kwargs): + kwargs.setdefault("evaluator", GridEvaluator( + series=self, at_cell_center=True)) + super().__init__(expr, range_x, range_y, range_z, label, **kwargs) + + @modify_parameterized_doc() class PlaneSeries(SurfaceBaseSeries): """Represents a plane in a 3D domain.""" diff --git a/tests/backends/make_tests.py b/tests/backends/make_tests.py index be537e0..0b92494 100644 --- a/tests/backends/make_tests.py +++ b/tests/backends/make_tests.py @@ -7,7 +7,7 @@ graphics, arrow_2d, arrow_3d, plot_root_locus, plot_pole_zero, ngrid, sgrid, zgrid, mcircles, surface, surface_parametric, line, root_locus, contour, list_2d, line_parametric_2d, hline, vline, - line_parametric_3d + line_parametric_3d, implicit_3d ) from spb.series import ( SurfaceOver2DRangeSeries, ParametricSurfaceSeries, LineOver1DRangeSeries, @@ -1110,3 +1110,29 @@ def make_test_colorbar_ticks_line_parametric_3d( ), backend=B, show=False ) + + +def make_test_implicit_3d_voxels_trisurf(B): + a = ( + (x ** 2 + y ** 2 + z ** 2 < 25) + & (x > 0) + & (x + 5 * y ** 2 > 4) + & (x ** 2 + y - z ** 3 > 4) + ) + b = ( + (y + z > -1) + & (y ** 2 < 2) + & (x ** 1 - y ** 2 < -5) + ) + + return graphics( + implicit_3d( + a | b, + (x, -5, 5), + (y, -5, 5), + (z, -5, 5), + n1=5, n2=7, n3=9 + ), + backend=B, + show=False + ) diff --git a/tests/backends/test_base_backend.py b/tests/backends/test_base_backend.py index bf64fda..fdb3bc6 100644 --- a/tests/backends/test_base_backend.py +++ b/tests/backends/test_base_backend.py @@ -160,6 +160,8 @@ def test_plot_sum(): second_col = np.array(p3.colorloop[1]) if line_col.dtype == second_col.dtype == " 1 diff --git a/tests/backends/test_plotly.py b/tests/backends/test_plotly.py index e7f7c66..0ee1e2d 100644 --- a/tests/backends/test_plotly.py +++ b/tests/backends/test_plotly.py @@ -96,6 +96,7 @@ make_test_colorbar_ticks_surface, make_test_colorbar_ticks_line_parametric_2d, make_test_colorbar_ticks_line_parametric_3d, + make_test_implicit_3d_voxels_trisurf, ) @@ -1981,3 +1982,11 @@ def test_colorbar_ticks_line_parametric_3d(): assert np.allclose( p2.fig.data[0].line.colorbar.tickvals, expected_tick_vals) + + +def test_implicit_3d_voxels_trisurf(): + p1 = make_test_implicit_3d_voxels_trisurf(PB) + assert len(p1.fig.data) == 1 + assert isinstance(p1.fig.data[0], go.Mesh3d) + mesh_x, mesh_y, mesh_z, f = p1[0].get_data() + assert all(t.shape == (4,6,8) for t in [mesh_x, mesh_y, mesh_z, f])