visualize.py: harden void mesh (per-axis grid clamp, vertex cap, full failure fallback to atoms+beads)

visualize.py

@@ -106,16 +106,18 @@ _METAL_ELEMENTS = {
 def _compute_void_mesh(
     structure,
     probe_radius: float = 1.2,
-    grid_spacing: float = 0.7,
+    target_grid_spacing: float = 0.7,
+    n_per_axis_min: int = 16,
+    n_per_axis_max: int = 60,
+    max_vertices: int = 40000,
     vdw_eff: float = 1.5,
 ):
     """Compute a Mercury-style contact-surface void mesh for the unit cell.
 
     Rolls a `probe_radius` probe (default 1.2 A, helium-sized) through the
-    unit cell on a `grid_spacing` grid (default 0.7 A) and runs marching
-    cubes at the void boundary. Mirrors Mercury's Display > Voids > Contact
-    Surface calculation (CCDC, https://www.ccdc.cam.ac.uk/discover/blog/
-    how-to-search-visualize-and-analyse-mof-structures).
+    unit cell and runs marching cubes at the void boundary. Mirrors
+    Mercury's Display > Voids > Contact Surface (CCDC, see
+    https://www.ccdc.cam.ac.uk/discover/blog/how-to-search-visualize-and-analyse-mof-structures).
 
     Atoms are replicated to a 3x3x3 super-image so distance queries respect
     periodic boundaries. A constant effective vdW radius of 1.5 A is used
@@ -123,59 +125,93 @@ def _compute_void_mesh(
     Zn etc.) the per-element value is within 0.2 A of this mean and the
     visual is unchanged.
 
-    Returns (vertices_cart, faces) where
+    Robustness behaviour. Per-axis grid resolution is clamped to
+    [`n_per_axis_min`, `n_per_axis_max`] so very large unit cells do not
+    blow up compute or browser payload, and very small cells still get a
+    decent surface. Mesh vertex count is capped at `max_vertices`; if the
+    surface would exceed that, the function returns (None, None) and the
+    caller falls back to a no-mesh render (atom colouring plus attribution
+    beads). All anticipated failure modes (missing scikit-image, missing
+    scipy, degenerate lattice, all-void or all-solid grid, marching-cubes
+    failure, oversize mesh) return (None, None).
+
+    Returns
+    -------
+    (vertices_cart, faces) where
         vertices_cart: np.ndarray [V, 3] in Cartesian Angstroms
-        faces:        np.ndarray [F, 3] of vertex indices per triangle
-    or (None, None) if scikit-image is unavailable, the mesh is empty, or
-    marching cubes fails.
+        faces:         np.ndarray [F, 3] of vertex indices per triangle
+    or (None, None).
     """
+    # Optional dependencies; degrade silently if either is missing
     try:
         from skimage.measure import marching_cubes
-    except ImportError:
+    except Exception:
         return None, None
     try:
         from scipy.spatial import cKDTree
-    except ImportError:
-        return None, None
-
-    lattice = np.asarray(structure.lattice.matrix, dtype=float)
-    a, b, c = structure.lattice.abc
-    nx = max(int(np.ceil(a / grid_spacing)), 16)
-    ny = max(int(np.ceil(b / grid_spacing)), 16)
-    nz = max(int(np.ceil(c / grid_spacing)), 16)
-
-    fx = np.linspace(0.0, 1.0, nx, endpoint=False)
-    fy = np.linspace(0.0, 1.0, ny, endpoint=False)
-    fz = np.linspace(0.0, 1.0, nz, endpoint=False)
-    grid_frac = np.stack(np.meshgrid(fx, fy, fz, indexing='ij'), axis=-1)
-    grid_cart = grid_frac.reshape(-1, 3) @ lattice
-
-    atom_cart = np.array([site.coords for site in structure], dtype=float)
-    shifts = []
-    for da in (-1, 0, 1):
-        for db in (-1, 0, 1):
-            for dc in (-1, 0, 1):
-                shifts.append(da * lattice[0] + db * lattice[1] + dc * lattice[2])
-    images = np.vstack([atom_cart + s for s in shifts])
-
-    tree = cKDTree(images)
-    cutoff = vdw_eff + probe_radius
-    dist, _ = tree.query(grid_cart, k=1)
-    void = (dist > cutoff).reshape(nx, ny, nz).astype(np.float32)
-
-    if void.max() < 0.5 or void.min() > 0.5:
+    except Exception:
         return None, None
 
     try:
-        verts_idx, faces, _, _ = marching_cubes(
-            void, level=0.5, allow_degenerate=False
-        )
-    except (ValueError, RuntimeError):
-        return None, None
+        lattice = np.asarray(structure.lattice.matrix, dtype=float)
+        a, b, c = structure.lattice.abc
+        if min(a, b, c) <= 0.0:
+            return None, None
+
+        # Per-axis resolution clamped so any cell stays within the compute
+        # and payload budget. Tiny cells get more samples than nominal
+        # spacing would suggest; huge cells get coarser sampling.
+        nx = int(np.clip(np.ceil(a / target_grid_spacing),
+                         n_per_axis_min, n_per_axis_max))
+        ny = int(np.clip(np.ceil(b / target_grid_spacing),
+                         n_per_axis_min, n_per_axis_max))
+        nz = int(np.clip(np.ceil(c / target_grid_spacing),
+                         n_per_axis_min, n_per_axis_max))
+
+        fx = np.linspace(0.0, 1.0, nx, endpoint=False)
+        fy = np.linspace(0.0, 1.0, ny, endpoint=False)
+        fz = np.linspace(0.0, 1.0, nz, endpoint=False)
+        grid_frac = np.stack(np.meshgrid(fx, fy, fz, indexing='ij'), axis=-1)
+        grid_cart = grid_frac.reshape(-1, 3) @ lattice
+
+        atom_cart = np.array([site.coords for site in structure], dtype=float)
+        if atom_cart.shape[0] == 0:
+            return None, None
+
+        shifts = []
+        for da in (-1, 0, 1):
+            for db in (-1, 0, 1):
+                for dc in (-1, 0, 1):
+                    shifts.append(da * lattice[0] + db * lattice[1] + dc * lattice[2])
+        images = np.vstack([atom_cart + s for s in shifts])
+
+        tree = cKDTree(images)
+        cutoff = vdw_eff + probe_radius
+        dist, _ = tree.query(grid_cart, k=1)
+        void = (dist > cutoff).reshape(nx, ny, nz).astype(np.float32)
+
+        if void.max() < 0.5 or void.min() > 0.5:
+            return None, None
+
+        try:
+            verts_idx, faces, _, _ = marching_cubes(
+                void, level=0.5, allow_degenerate=False
+            )
+        except Exception:
+            return None, None
 
-    verts_frac = verts_idx / np.array([nx, ny, nz], dtype=float)
-    verts_cart = verts_frac @ lattice
-    return verts_cart.astype(float), faces.astype(int)
+        # Vertex cap so the JSON payload shipped to the iframe stays
+        # bounded even for unusual very-high-porosity structures
+        if verts_idx.shape[0] > max_vertices:
+            return None, None
+
+        verts_frac = verts_idx / np.array([nx, ny, nz], dtype=float)
+        verts_cart = verts_frac @ lattice
+        return verts_cart.astype(float), faces.astype(int)
+    except Exception:
+        # Catch-all so a malformed structure never breaks the rest of the
+        # rendering pipeline. Caller falls back to no-mesh render.
+        return None, None
 
 
 def _cluster_pore_vertices(
@@ -447,17 +483,23 @@ def create_3d_visualization(
         verts_cart, faces = None, None
 
     if verts_cart is not None and faces is not None and len(faces) > 0:
-        vertex_arr = [
-            {'x': float(v[0]), 'y': float(v[1]), 'z': float(v[2])}
-            for v in verts_cart
-        ]
-        face_arr = faces.flatten().tolist()
-        view.addCustom({
-            'vertexArr': vertex_arr,
-            'faceArr': face_arr,
-            'color': '#fbbf24',  # warm Mercury-style yellow
-            'opacity': 0.25,
-        })
+        try:
+            vertex_arr = [
+                {'x': float(v[0]), 'y': float(v[1]), 'z': float(v[2])}
+                for v in verts_cart
+            ]
+            face_arr = faces.flatten().tolist()
+            view.addCustom({
+                'vertexArr': vertex_arr,
+                'faceArr': face_arr,
+                'color': '#fbbf24',  # warm Mercury-style yellow
+                'opacity': 0.25,
+            })
+        except Exception:
+            # Mesh shape construction failed (e.g. py3Dmol API mismatch on
+            # some build). Skip the mesh and continue with the rest of the
+            # render so the user still sees atoms and attribution beads.
+            pass
 
     # Per-cavity attribution glyphs: small opaque beads at cluster centres.
     # Size scales with |attribution| so the strongest contributing cavity