from __future__ import annotations
from dataclasses import dataclass
import numpy as np
@dataclass(frozen=True)
class VoxelRegion:
"""Summary of one connected voxel region."""
label: int
voxel_count: int
volume: float
surface_area: float
[docs]
@dataclass(frozen=True)
class ProbeAccessibleResult:
"""Probe-center accessible volume and surface summary."""
probe_radius: float
accessible_voxel_count: int
accessible_volume: float
accessible_surface_area: float
regions: list[VoxelRegion]
accessible_mask: np.ndarray
[docs]
class VoxelGridAnalysis:
"""Analyze connected voxel volumes and their surfaces."""
def __init__(self, voxel_grid):
self.voxel_grid = voxel_grid
@property
def grid(self):
return self.voxel_grid.to_numpy()
@property
def cell(self):
return self.voxel_grid.cell
@property
def gpts(self):
return self.voxel_grid.gpts
@property
def voxel_volume(self):
return abs(float(np.linalg.det(self.cell))) / float(np.prod(self.gpts))
[docs]
def mask(self, min_value=None, max_value=None, threshold=None, above=True):
"""Build a boolean mask from voxel values."""
if threshold is not None and (min_value is not None or max_value is not None):
raise ValueError("Specify either threshold or min_value/max_value, not both")
grid = self.grid
if threshold is not None:
return grid > threshold if above else grid < threshold
selected = np.ones(grid.shape, dtype=bool)
if min_value is not None:
selected &= grid >= min_value
if max_value is not None:
selected &= grid <= max_value
return selected
[docs]
def connected_components(self, selected, connectivity=1, periodic=True):
"""Label connected components in a boolean mask."""
try:
from skimage.measure import label
except ImportError as exc: # pragma: no cover - depends on optional dependency
raise ImportError(
"VoxelGridAnalysis requires scikit-image. Install it with "
"`pip install AtomVoxelizer[analysis]`."
) from exc
selected = np.asarray(selected, dtype=bool)
labels = label(selected, connectivity=connectivity)
if periodic:
labels = self._merge_periodic_labels(labels, selected)
return labels, int(labels.max())
[docs]
def region_volume(self, selected):
"""Return the volume represented by a boolean mask."""
return int(np.count_nonzero(selected)) * self.voxel_volume
[docs]
def surface_area(self, selected, periodic=True):
"""Estimate the surface area of a boolean region with marching cubes."""
try:
from skimage.measure import marching_cubes
except ImportError as exc: # pragma: no cover - depends on optional dependency
raise ImportError(
"VoxelGridAnalysis requires scikit-image. Install it with "
"`pip install AtomVoxelizer[analysis]`."
) from exc
selected = np.asarray(selected, dtype=bool)
if not np.any(selected):
return 0.0
if periodic and np.all(selected):
return 0.0
if periodic:
values = np.tile(selected.astype(np.float32), (3, 3, 3))
vertices, faces, _normals, _values = marching_cubes(values, level=0.5)
central_offset = self.gpts
centroids = vertices[faces].mean(axis=1)
in_central_cell = np.all((centroids >= central_offset) & (centroids < 2 * central_offset), axis=1)
faces = faces[in_central_cell]
vertices = vertices - central_offset
else:
values = np.pad(selected.astype(np.float32), 1, mode="constant", constant_values=0.0)
vertices, faces, _normals, _values = marching_cubes(values, level=0.5)
vertices = vertices - 1.0
if faces.size == 0:
return 0.0
real_vertices = self._index_vertices_to_real(vertices)
return self._mesh_surface_area(real_vertices, faces)
[docs]
def surface_area_voxel_faces(self, selected, periodic=True):
"""Estimate surface area by counting exposed voxel faces.
This is faster than marching cubes and avoids tiled periodic volumes. It
is a grid-face estimate rather than a smoothed triangular surface, so it
is most useful for fast convergence scans or large grids.
"""
selected = np.asarray(selected, dtype=bool)
if not np.any(selected):
return 0.0
if periodic and np.all(selected):
return 0.0
voxel_vectors = self.cell / self.gpts[:, None]
face_areas = np.array(
[
np.linalg.norm(np.cross(voxel_vectors[1], voxel_vectors[2])),
np.linalg.norm(np.cross(voxel_vectors[0], voxel_vectors[2])),
np.linalg.norm(np.cross(voxel_vectors[0], voxel_vectors[1])),
],
dtype=float,
)
area = 0.0
for axis, face_area in enumerate(face_areas):
if periodic:
exposed = selected != np.roll(selected, -1, axis=axis)
area += float(np.count_nonzero(exposed)) * face_area
else:
inner_a = [slice(None)] * selected.ndim
inner_b = [slice(None)] * selected.ndim
inner_a[axis] = slice(None, -1)
inner_b[axis] = slice(1, None)
exposed = selected[tuple(inner_a)] != selected[tuple(inner_b)]
area += float(np.count_nonzero(exposed)) * face_area
first = np.take(selected, 0, axis=axis)
last = np.take(selected, -1, axis=axis)
area += float(np.count_nonzero(first) + np.count_nonzero(last)) * face_area
return area
[docs]
def mesh_at_value(self, level, periodic=True, clip_periodic=True):
"""Return a marching-cubes mesh for a scalar voxel value.
``level`` is interpreted in the units stored in the voxel grid. For a
distance-mask grid this traces the surface at a fixed distance from the
nearest atom. Vertices are returned in real-space coordinates and faces
index into that vertex array. Periodic meshes are clipped to the primary
cell by default so boundary-crossing triangles are cut at the cell edge
instead of being wrapped across the cell.
"""
try:
from skimage.measure import marching_cubes
except ImportError as exc: # pragma: no cover - depends on optional dependency
raise ImportError(
"VoxelGridAnalysis requires scikit-image. Install it with "
"`pip install scikit-image`."
) from exc
level = float(level)
grid = self._finite_grid_for_level(level)
finite = grid[np.isfinite(grid)]
if finite.size == 0 or level < finite.min() or level > finite.max():
raise ValueError("level must lie within the finite voxel value range")
if periodic:
values = np.tile(grid.astype(np.float32), (3, 3, 3))
vertices, faces, _normals, _values = marching_cubes(values, level=level)
central_offset = self.gpts
centroids = vertices[faces].mean(axis=1)
in_central_cell = np.all((centroids >= central_offset) & (centroids < 2 * central_offset), axis=1)
faces = faces[in_central_cell]
vertices = vertices - central_offset
if clip_periodic:
vertices, faces = self._clip_mesh_to_index_cell(vertices, faces)
else:
vertices, faces, _normals, _values = marching_cubes(grid.astype(np.float32), level=level)
if faces.size == 0:
return np.empty((0, 3), dtype=float), np.empty((0, 3), dtype=int)
used = np.unique(faces)
remap = np.full(vertices.shape[0], -1, dtype=int)
remap[used] = np.arange(used.shape[0])
real_vertices = self._index_vertices_to_real(vertices[used])
return real_vertices, remap[faces]
[docs]
def surface_area_at_value(self, level, periodic=True, clip_periodic=True):
"""Estimate the area of a scalar isosurface at ``level``."""
vertices, faces = self.mesh_at_value(level=level, periodic=periodic, clip_periodic=clip_periodic)
if faces.size == 0:
return 0.0
return self._mesh_surface_area(vertices, faces)
[docs]
def analyze_regions(
self,
min_value=None,
max_value=None,
threshold=None,
above=True,
connectivity=1,
periodic=True,
surface_method="marching-cubes",
):
"""Return volume and marching-cubes area for each connected region."""
selected = self.mask(min_value=min_value, max_value=max_value, threshold=threshold, above=above)
return self.analyze_mask(
selected,
connectivity=connectivity,
periodic=periodic,
surface_method=surface_method,
)
[docs]
def analyze_mask(
self,
selected,
connectivity=1,
periodic=True,
surface_method="marching-cubes",
):
"""Return volume and surface area for connected regions in a boolean mask."""
if surface_method not in {"marching-cubes", "voxel-faces"}:
raise ValueError("surface_method must be 'marching-cubes' or 'voxel-faces'")
selected = np.asarray(selected, dtype=bool)
labels, label_count = self.connected_components(selected, connectivity=connectivity, periodic=periodic)
regions = []
for label_id in range(1, label_count + 1):
region_mask = labels == label_id
voxel_count = int(np.count_nonzero(region_mask))
if surface_method == "marching-cubes":
surface_area = self.surface_area(region_mask, periodic=periodic)
else:
surface_area = self.surface_area_voxel_faces(region_mask, periodic=periodic)
regions.append(
VoxelRegion(
label=label_id,
voxel_count=voxel_count,
volume=voxel_count * self.voxel_volume,
surface_area=surface_area,
)
)
return regions
[docs]
def probe_accessible_mask(self, positions, radii, probe_radius, write_grid=False):
"""Return voxels accessible to the center of a spherical probe.
``positions`` must have shape ``(N, 3)`` and ``radii`` must contain one
atom radius per position. A voxel is accessible when a probe center at
that voxel does not overlap any atom, so atoms are excluded with radius
``radii + probe_radius``.
The input voxel grid is used as the geometry template. By default, the
existing grid values are not modified. Set ``write_grid=True`` to store
the binary accessible mask in ``self.voxel_grid.grid`` as 1 for
accessible voxels and 0 for excluded voxels.
"""
positions, radii = self._validate_probe_inputs(positions, radii, probe_radius)
from .voxelgrid import VoxelGrid
work_grid = VoxelGrid(self.cell, gpts=self.gpts, dtype=np.uint8)
work_grid.grid.fill(1)
work_grid.set_spheres(positions, radii + float(probe_radius), value=0)
accessible = work_grid.grid.astype(bool, copy=True)
if write_grid:
self.voxel_grid.grid[...] = accessible.astype(self.voxel_grid.dtype, copy=False)
return accessible
[docs]
def analyze_probe_accessibility(
self,
positions,
radii,
probe_radius,
connectivity=1,
periodic=True,
surface_method="voxel-faces",
write_grid=False,
):
"""Analyze probe-center accessible volume and surface area.
This is a probe-center analysis: atom exclusion radii are inflated by
``probe_radius`` and accessible voxels are positions where the probe
center fits. Volumes therefore describe the region available to the
probe center. Surface areas describe the boundary of that accessible
center region.
"""
accessible = self.probe_accessible_mask(positions, radii, probe_radius, write_grid=write_grid)
regions = self.analyze_mask(
accessible,
connectivity=connectivity,
periodic=periodic,
surface_method=surface_method,
)
accessible_voxel_count = int(np.count_nonzero(accessible))
return ProbeAccessibleResult(
probe_radius=float(probe_radius),
accessible_voxel_count=accessible_voxel_count,
accessible_volume=accessible_voxel_count * self.voxel_volume,
accessible_surface_area=sum(region.surface_area for region in regions),
regions=regions,
accessible_mask=accessible,
)
[docs]
def probe_accessible_surface_area(
self,
positions,
radii,
probe_radius,
samples_per_atom=1000,
surface_radius_scale=1.0,
):
"""Estimate probe-accessible surface area by sampling inflated atom surfaces.
Points are placed deterministically on each sphere with a Fibonacci
sphere rule. A sampled point contributes to the area when it does not
overlap the inflated sphere around any other atom under periodic
boundary conditions.
``surface_radius_scale`` defaults to 1.0 for a hard-sphere contact
surface. PoreBlazer uses 1.122 for its nitrogen accessible surface area
calculation.
"""
positions, radii = self._validate_probe_inputs(positions, radii, probe_radius)
samples_per_atom = int(samples_per_atom)
surface_radius_scale = float(surface_radius_scale)
if samples_per_atom < 1:
raise ValueError("samples_per_atom must be at least 1")
if surface_radius_scale <= 0.0:
raise ValueError("surface_radius_scale must be positive")
directions = self._fibonacci_sphere(samples_per_atom)
surface_radii = surface_radius_scale * (radii + float(probe_radius))
total_area = 0.0
for atom_index, (position, surface_radius) in enumerate(zip(positions, surface_radii)):
points = position + directions * surface_radius
accessible = np.ones(samples_per_atom, dtype=bool)
for other_index, (other_position, other_radius) in enumerate(zip(positions, surface_radii)):
if other_index == atom_index:
continue
disp_frac = (points - other_position) @ self.voxel_grid.cell_inv
disp_frac -= np.round(disp_frac)
disp = disp_frac @ self.cell
dist2 = np.einsum("ij,ij->i", disp, disp)
accessible &= dist2 >= other_radius * other_radius
total_area += 4.0 * np.pi * surface_radius * surface_radius * float(np.mean(accessible))
return float(total_area)
[docs]
@staticmethod
def volume_angstrom3_to_cm3_per_g(volume_angstrom3, mass_amu):
"""Convert a cell/supercell volume from Angstrom^3 to cm^3/g."""
if mass_amu <= 0:
raise ValueError("mass_amu must be positive")
mass_g = float(mass_amu) * 1.66053906660e-24
return float(volume_angstrom3) * 1.0e-24 / mass_g
[docs]
@staticmethod
def area_angstrom2_to_m2_per_g(area_angstrom2, mass_amu):
"""Convert a cell/supercell area from Angstrom^2 to m^2/g."""
if mass_amu <= 0:
raise ValueError("mass_amu must be positive")
mass_g = float(mass_amu) * 1.66053906660e-24
return float(area_angstrom2) * 1.0e-20 / mass_g
@staticmethod
def _validate_probe_inputs(positions, radii, probe_radius):
positions = np.asarray(positions, dtype=np.float64)
radii = np.asarray(radii, dtype=np.float64)
probe_radius = float(probe_radius)
if positions.ndim != 2 or positions.shape[1] != 3:
raise ValueError("positions must have shape (N, 3)")
if radii.ndim != 1 or radii.shape[0] != positions.shape[0]:
raise ValueError("radii must have shape (N,)")
if np.any(radii < 0.0):
raise ValueError("radii must be non-negative")
if probe_radius < 0.0:
raise ValueError("probe_radius must be non-negative")
return positions, radii
@staticmethod
def _fibonacci_sphere(count):
indices = np.arange(count, dtype=float) + 0.5
z = 1.0 - 2.0 * indices / float(count)
theta = np.pi * (1.0 + np.sqrt(5.0)) * indices
radius = np.sqrt(np.maximum(0.0, 1.0 - z * z))
return np.column_stack((radius * np.cos(theta), radius * np.sin(theta), z))
def _index_vertices_to_real(self, vertices):
frac = vertices / self.gpts
return frac @ self.cell
def _finite_grid_for_level(self, level):
grid = np.asarray(self.grid, dtype=np.float32)
if np.all(np.isfinite(grid)):
return grid
finite = grid[np.isfinite(grid)]
if finite.size == 0:
return grid
high = max(float(finite.max()), float(level)) + 1.0
low = min(float(finite.min()), float(level)) - 1.0
clean = grid.copy()
clean[np.isposinf(clean)] = high
clean[np.isneginf(clean)] = low
return clean
def _clip_mesh_to_index_cell(self, vertices, faces):
clipped_vertices = []
clipped_faces = []
bounds = [(0.0, float(n)) for n in self.gpts]
for face in faces:
polygon = [vertices[int(vertex_index)].astype(float) for vertex_index in face]
for axis, (low, high) in enumerate(bounds):
polygon = self._clip_polygon_axis(polygon, axis, low, keep_greater=True)
polygon = self._clip_polygon_axis(polygon, axis, high, keep_greater=False)
if len(polygon) < 3:
break
if len(polygon) < 3:
continue
start = len(clipped_vertices)
clipped_vertices.extend(polygon)
for index in range(1, len(polygon) - 1):
clipped_faces.append((start, start + index, start + index + 1))
if not clipped_faces:
return np.empty((0, 3), dtype=float), np.empty((0, 3), dtype=int)
return np.asarray(clipped_vertices, dtype=float), np.asarray(clipped_faces, dtype=int)
@staticmethod
def _clip_polygon_axis(polygon, axis, boundary, keep_greater):
if not polygon:
return []
clipped = []
previous = polygon[-1]
previous_inside = previous[axis] >= boundary if keep_greater else previous[axis] <= boundary
for current in polygon:
current_inside = current[axis] >= boundary if keep_greater else current[axis] <= boundary
if current_inside != previous_inside:
delta = current[axis] - previous[axis]
if delta != 0.0:
t = (boundary - previous[axis]) / delta
clipped.append(previous + t * (current - previous))
if current_inside:
clipped.append(current)
previous = current
previous_inside = current_inside
return clipped
@staticmethod
def _merge_periodic_labels(labels, selected):
label_count = int(labels.max())
if label_count == 0:
return labels
parent = np.arange(label_count + 1)
def find(label_id):
while parent[label_id] != label_id:
parent[label_id] = parent[parent[label_id]]
label_id = parent[label_id]
return label_id
def union(a, b):
if a == 0 or b == 0:
return
root_a = find(int(a))
root_b = find(int(b))
if root_a != root_b:
parent[root_b] = root_a
for axis in range(labels.ndim):
first_labels = np.take(labels, 0, axis=axis)
last_labels = np.take(labels, -1, axis=axis)
first_selected = np.take(selected, 0, axis=axis)
last_selected = np.take(selected, -1, axis=axis)
for a, b in zip(first_labels[first_selected & last_selected], last_labels[first_selected & last_selected]):
union(a, b)
root_to_new_label = {}
next_label = 1
merged = np.zeros_like(labels)
for label_id in range(1, label_count + 1):
root = find(label_id)
if root not in root_to_new_label:
root_to_new_label[root] = next_label
next_label += 1
merged[labels == label_id] = root_to_new_label[root]
return merged
@staticmethod
def _mesh_surface_area(vertices, faces):
triangles = vertices[faces]
cross = np.cross(triangles[:, 1] - triangles[:, 0], triangles[:, 2] - triangles[:, 0])
return float(0.5 * np.linalg.norm(cross, axis=1).sum())
[docs]
@staticmethod
def mesh_surface_area(vertices, faces):
"""Return the area of a triangular mesh."""
return VoxelGridAnalysis._mesh_surface_area(np.asarray(vertices), np.asarray(faces, dtype=int))
__all__ = ["ProbeAccessibleResult", "VoxelGridAnalysis", "VoxelRegion"]
