CaviTracer - refactor for enhanced accuracy and speed of channel identification to be competitive with established SOTA tools by briza81 · Pull Request #2249 · prody/ProDy

briza81 · 2026-07-16T04:59:30Z

Accuracy/correctness improvements in the current branch over the original implementation and compared to CAVER3

1. Voronoi diagrams (VD) - the atoms are no longer treated as equal-radius points

The original ran a plain Delaunay/Voronoi tessellation on atom centers (and auto-excluded hetero/H atoms). Now there's a diagram argument with three modes:

"homogenized" (new default) - each atom is replaced by concentric shells of equal-radius balls (max_deviation controls fidelity), so an ordinary tessellation approximates the additively-weighted (Apollonius) diagram within the max_deviation tolerance (MolAxis/CAVER3 trick)
"simple" - the original point-based behavior.
"weighted" - the true Apollonius diagram via the new _vorpy_aw.py (third-party vorpy + optional numba kernel, spatial box decomposition, content-keyed on-disk .npz cache, weighted_cache=True). However, this is super slow, circa 100x over high-accuracy homogenized VD
-> Because atom radii are now properly treated, we no longer need to drop hetero/H atoms.

2. Width/bottleneck/Cost function- now they are sampled at the gates, not just the vertices

The original measured channel clearance only at Voronoi vertices (circumcenters), which are local maxima of width. Here we move to a per-edge gate clearance (min over the shared Delaunay face, closed-form point-to-segment) + integrating the clearance along the actual edge:

Dijkstra cost is now driven by the gate clearance, making the routing cost direction-symmetric. Edge clearance integral provides the ultimate accuracy for the edge cost.
Reported bottleneck now takes the min clearance on the edge.
Volume tube injects a gate knot midway between vertices so the radius profile dips where the channel is actually narrow.

3. limiting mesh sensitivity - avoiding layer-based metrices that change with VD density, which changes with max_deviation

migrate to physical metrics - distance/topological distance, instead of relying on layers, which are r1, r2 and max_deviation dependent
hence, depth is now measured in geodesic Ångström, not tetrahedron layers -> a multi-source Dijkstra along Voronoi edges (with a scale-invariant). Defaults recalibrated to match old behavior for both channels and surface cavities (min_depth 10→5 Å, etc.).

4. surface/cavity definition - local enclosure-based peel/erosion

The original erosion depended on the tessellation and on r1. Now, a boundary tetrahedron is stripped only while it's open (based on ray-marched directions that hit protein within range), so burial is measured against real atoms. r1 only decides where erosion starts.
Importantly, mouths absorb instead of conduct. Exit tetrahedra are no longer pre-sampled targets (thinned by sparsity); before the Dijkstra search, every mouth's outgoing edges are deleted so a surface tetrahedron, a probe of radius r2 can leave through becomes a sink (incoming edges are kept, outgoing edges are removed). This makes "a channel ends at its first surface contact" a hard constraint of the search rather than a cut applied afterward, avoiding the mesh-dependent failure where a genuine narrow interior corridor lost the argmin to a wide surface groove and vanished from the output. Now, every reachable mouth emits a candidate, and their deduplication decides which mouths are one opening: when two channels reach the same opening, we test if they have similar corridor/path to this opening (similarity within route_tolerance), and if so => two channels reaching the same exit via the same path -> we keep only the cheaper one. ->
Therefore, sparsity is not needed anymore to control against channel explosion. It is now an argument to fine-tune the minimal required distance between channel openings of distinct channels.

5. start-point handling

restrict_channels_to_start_point limits the search only to the anchored cavity
start_point_search - optimize the seed to the widest nearby tetrahedron of at least the same depth, which is accessible from the start point, instead of using the nearest one

6. performance optimization

Here we have several overlaps with #2241 by Matthew, to name a few:
vectorized tessellation filtering, surface erosion restricted to the boundary shell (O(n)→O(n^⅔) per pass), circumcenters recovered from the Delaunay paraboloid lifting (skips a second Qhull pass), one CSR weighted graph + a single scipy multi-target Dijkstra per cavity (replacing per-(seed, exit) heap Dijkstra), vectorized Simpson volume, and ~30× faster vectorized buildSparseGraph

7. output & reporting

channels now carry a per-channel cost (Dijkstra path weight) and curvature, are ordered by ascending cost (channel 0 = best tunnel), and each gets a REMARK line in the PQR
multi-channel PQR now has corrected CONNECT items, and separate resIDs - same as the PQR files for cavities had.
added per-stage LOGGER timing
local _warn() workaround for a global filter that was swallowing WARNINGs

8. Recalibrated defaults to losely agree with CAVER3 outputs

r2 1.25→0.9, min_depth 10→5, bottleneck 1→0.0, sparsity 15→1, plus the new diagram/max_deviation/min_enclosure args

More details, execution time measurements, and use cases for accuracy are available in the attached file: Refoactored_CaviTracer_comparisons.pdf

…Voronoi by optionally replacing each atom by equal-radius balls (rho = smallest vdW radius) on concentric shells, so an ordinary Delaunay/Voronoi tessellation approximates the additively weighted Voronoi diagram (approach used by MolAxis / CAVER 3) instead of discarding small atoms. Adds homogenize=True and max_deviation=0.2 to calcChannels plus homogenize_atoms / _fibonacci_sphere / _shell_point_count. Also, we can stop auto-excluding hetero/H atoms as it is no longer needed.

…zing ops in delete_simplices3d, delete_section,surface_layer and the neighbour remap (row membership via hashed void-view _rows_isin). Recover Voronoi vertices from the Delaunay paraboloid lifting (calc_circumcenters), skipping a second Qhull pass. Replace per-(seed,exit) heap Dijkstra with one CSR weighted graph + a single multi-target scipy.csgraph Dijkstra per cavity. Tube volume by vectorized composite Simpson instead of adaptive scipy.quad. Numerically equivalent for both channels and cavities; large speedup.

…ing of wall-clock timings for homogenization, tessellation, surface filtering, cavity detection and pathfinding, plus an overall total, via LOGGER.timeit/report. No behavioural change.

…hich is active when a start_point is set. False keeps the existing behavior (seed every cavity, search all channels); True restricts the search to the single cavity whose closest tetrahedron is globally nearest to start_point.

…n delete_simplices3d. The surface=True pass only ever deletes tetrahedra with a -1 neighbour, but the vectorized code computed np.linalg.norm over all n tetrahedra every erosion iteration and masked to the boundary afterwards (O(n) per pass). Compute the norm only on boundary rows (~n^(2/3)) and default interior tetrahedra to keep; the non-surface pass still needs all n and is unchanged. Numerically identical to the original (verified bit-identical simp/neigh/verti vs the unvectorized channels_original across 400 randomized Delaunay cases, both surface flags). Removes the last full-array cost that made erosion the crossover stage at large point counts (homogenize).

…e, order channels by ascending cost so channel 0 is the best tunnel. Write a REMARK line with length/bottleneck/curvature/cost before each channel in the output PQR file. Channel set and geometry unchanged.

…ion. 1. Surface truncation + de-duplication of channels. New calcChannels params truncate_at_surface=True and similarity=0.8, towards dijkstra. The Dijkstra cost only rewards width, so a cheapest path to a far exit can run through or past a nearer mouth. Each reconstructed seed->exit path is now cut at the first qualified surface mouth it enters -- a surface (exit) tetrahedron whose inscribed clearance is >= bottleneck. Channel properties are recalculated up to the truncated surface terminal. Truncated paths are de-duplicated by merging two when their exit spheres overlap (|Ti - Tj| < ri + rj) AND they share most of their route as a common prefix from the seed (>= similarity). Distinct mouths, and different corridors reaching one mouth (diverging early), are kept as separate tunnels. 2. Vectorized get_end_tetrahedra. The greedy sparse sampling of mouth tetrahedra (seed with the widest, then repeatedly add the widest tetra still >= sparsity from every selected one) was O(N_exit x M^2). This was replaced with a running min-distance array folded with one vectorized norm per pick, and inscribed radii (geometry-only) precomputed once. producing bit-identical results, while removing considerable load

…eighted")

…oat" by eroding the r1 surface inward round(r1-r2) layers with the r2 probe so the wide former-exterior shell a large r1 bridges over can't act as a low-cost path that collapsed channels at high r1. max_peel_depth caps it (None=uncapped, 0=off).

… results with Caver3 defaults

…ming function names to mixedCase and :param: to :arg:, line length limits

…onoi by implements the true additively-weighted (Apollonius) Voronoi diagram as an alternative to the "homogenized" approximation, via the new _vorpy_aw.py (third-party vorpy package + optional compiled numba kernel, spatial box decomposition, and a content-keyed on-disk cache).

…ions with preceived per-frame parallelism

…dence by perforimng the erosion with a local stop: a boundary tetrahedron is stripped only while it is open, i.e. fewer than min_enclosure of the directions leaving it meet protein within 15 A, ray-marched against the real atoms so burial does not depend on the tessellation. r1 then decides only where erosion starts, not where it stops, and is left capping the mouths. Also fix the dedup reporting one tunnel as a bundle of near-copies: a tunnel splays as it widens into its opening, and that fan was counted as corridor divergence. _routeCoverage now discards the shared opening, and cutting, matching and comparing use one definition of the opening instead of three.

…ies. So far, sparsity fed two places, but they are never both live: the dedup opening floor (truncate_at_surface=True) and the getEndTetrahedra terminal spacing (False). Neither touches cavities - identical cavities at 1 and 15 - so calcSurfaceCavities now ignores it (still accepted, deprecated) and calcChannels keeps one sparsity=1. Also fixed calcSurfaceCavities which inherited min_enclosure=0.85 => peel strips shallow open regions, which is what a pocket is, hence we need to stop peeling there.

…to 0.70, and ENCLOSURE_RANGE 15 -> 25. Also warn when a structure has no hydrogens and r2 < 1.2: the void the missing H leave is then wide enough for the probe, and the interior percolates into a sponge => channel number explode. At r2 >= 1.2 protonated and X-ray input agree.

… tetrahedron by finding the seed as the widest tetrahedron of the same cavity within start_point_search (new arg, default 3 A) that is no shallower than the original closet tetrahedron and reachable from it through that neighbourhood. The cavity is still chosen by the anchor, never by the widened seed. Also report the seed actually used(vertex, distance, radius, depth) and the seed it replaced.

…ce of width — per-simplex vertex clearance (spline knots) and per-edge gate clearance (min over the shared Delaunay face, closed-form point-to-segment). buildSparseGraph fills both once, calculateRadiusSpline reads the vertex cache instead of recomputing per path; the edge map is inert here (nothing reads it yet), foundation for making the reported bottleneck and Dijkstra cost the paper's edge bottleneck radius rather than circumcenter-only sampling. Verified bitwise-identical channels.

…ce — calculateRadiusSpline now takes the path minimum over the per-edge gate clearances (shared-face min between consecutive circumcenters) via _pathBottleneck, reading the cached map and recomputing any missing edge, instead of min over the vertices which are local clearance maxima. Routing and volume untouched; the reported bottleneck can only tighten, so filterChannelsByBottleneck stops passing channels that are truly sub-threshold at a face it never sampled.

…red node's vertex clearance — buildSparseGraph now feeds d = shared-face gate (the cached edge bottleneck, symmetric) into l/(d^2+b), so the search stops preferring routes that are wide only at their circumcenters and narrow at a face they never measured. Cost is now direction-symmetric. Only the cost changes here; bottleneck and volume are computed as before, though rerouted channels naturally report different values.

…ects each edge gate clearance as an extra radius-spline knot midway between its two vertices, so the radius profile dips where the channel is actually narrow instead of interpolating only the wide circumcenters. _pathBottleneck becomes _pathGates returning the per-edge array, shared by the reported bottleneck (its min) and these knots. Centerline knots and spline domain unchanged, so length, bottleneck, routes and cap radii are identical; only volume moves

…edron layers, so min_depth/max_depth/weighted_mouth_depth stop drifting with mesh density (a layer is one tetrahedron thick, so the old BFS layer count inflated as max_deviation shrank. Depth is now a multi-source Dijkstra from the opening along Voronoi edges; edges touching a near-flat tetrahedron (runaway circumcenter) are dropped via a scale-invariant R/L>5 flatness flag. Recalibrated defaults against old behavior: min_depth 10->5 A, weighted_mouth_depth 4->2.5 A, calcSurfaceCavities 2/3->1.5/2.5 A. Deepest cavity now stable across the mesh; also fixes the param-file Depth [A] column, which wrote layer counts.

…cy from array ops over the (N, deg) neighbour table instead of the Python double loop, eliminating the per-tetra frozenset intersection and per-edge _edgeBottleneck dispatch. Directed edges come straight off neighbors.ravel(); shared Delaunay faces via a convention-agnostic broadcast intersection (still not scipy's opposite-vertex rule, so it holds for the weighted diagram); all gate clearances in one batched _edgeBottleneckBatch, with the twin-tetrahedron guard applied row-wise. Symmetry preserved by driving each edge from its lower-index endpoint. Verified across the max_deviation 0.2→0.01 sweep: identical CSR structure and channels (paths/length/bottleneck/volume/cost), weights agree to 2e-13 (einsum vs BLAS reduction order), ~30x faster on the graph build.

… — saveChannelsToPdb numbered ATOM serials globally but emitted CONECT with local indices (range(1, samples)), so every channel after the first re-bonded channel 0's atoms and left its own unbonded, and all channels shared FIL T 1. Combined output now bonds each channel over its own global serials (range(atom_index, atom_index+samples-1)) with no CONECT spanning channels, and gives each channel its own residue (FIL T%4d), matching saveCavitiesToPdb.

…e-integral Dijkstra cost, as the default for homogenized/simple; keep legacy l/(d^2+b) as 'bottleneck' and the default for weighted ('integral'+'weighted' errors). The old cost charges the whole edge at its narrowest gate (MOLE-style), not additive under subdivision, so routing drifts as max_deviation coarsens. Integrating r(t)^-z along the edge (profile-integral) is additive and mesh-invariant. Trapezoid over a fixed 0.3A arclength grid plus a forced node at the exact gate; reported bottleneck unchanged. Cost-only change. At worst, only negligible timing differences.

briza81 · 2026-07-16T05:06:39Z

@karolamik13 and @MatthewLicht pls have a look, more details over email.

@jamesmkrieger Karolina suggested asking if you could review this, too thank you

jamesmkrieger · 2026-07-16T08:15:42Z

Thank @briza81, from the summary these sound like great changes and I will take a look at the details soon

jamesmkrieger

Generally looks good. There are a few minor things of a more stylistic nature that I don't need to review again. I trust you know what you're doing on the details and it generally looks clear and well commented and documented

Addresses the stylistic review comments on prody#2249, plus the further typos an aspell pass over the file's comments and docstrings turned up. Comments and docstrings only; the one code change is a rewrap of list comperhension with identical code: - getVmdModel: "data a / nd uses VMD" was split mid-word across lines. - "developement" -> "development" in the three Open3D install notes - showChannels: wrap the create_sphere list comprehension more readably - calcChannels: drop the blank line between the diagram description and its :type: field. - calcChannels: generalize the no-hydrogen rationale from "an X-ray structure carries no hydrogens" to experimental structures generally, X-ray and cryo-EM alike, since neither resolves H except at the very highest resolutions. - Additional typos and formatting issues corrected, mostly in the older docstrings, predating this PR.

This should address the remaining review comments on prody#2249. Pure refactor, no behaviour change. The commented-out scalar sphereFit becomes a real vectorized method returning a "probe fits" boolean mask, with an optional `rows` mask that keeps the surface pass restricted to the boundary shell (~n^(2/3) rows). deleteSimplices3d now calls it: `should_delete = fits` when eroding the surface, `~fits` when carving the interior. Both methods get docstrings.

briza81 · 2026-07-16T12:22:46Z

@jamesmkrieger thank you for feedback, all should be resolved now
@karolamik13 should be ready for you

Docs changes (CaviTracer name update, Jan Brezovsky is added as an author of the code)

karolamik13 · 2026-07-17T16:21:22Z

After improvements from Jan and Matthew, calcChannelSurfaceOverlaps() becomes very slow due to the significantly larger number of channels generated by the new approach. Therefore, I added options to speed up calcChannelSurfaceOverlaps() and calcSurfaceCavitiesOverlaps() by enabling max_proc and parallel calculations.

jamesmkrieger · 2026-07-17T17:17:31Z

sounds good

karolamik13 · 2026-07-18T16:24:51Z

@briza81 @MatthewLicht

input atoms filtering from TODO list is now added

…msInputComposition is added to provide information about atoms composition that is used to compute channels. Waters are not excluded from selection in calcChannels() but other components will be used. The user is informed that other components are taken into account.

briza81 added 27 commits July 7, 2026 17:32

channels: per-stage LOGGER timing in calcChannels with unified report…

0504b5b

…ing of wall-clock timings for homogenization, tessellation, surface filtering, cavity detection and pathfinding, plus an overall total, via LOGGER.timeit/report. No behavioural change.

channels: cleaning comments

d53d53f

channels: record per-channel cost (Dijkstra path weight) and curvatur…

f30fcc7

…e, order channels by ascending cost so channel 0 is the best tunnel. Write a REMARK line with length/bottleneck/curvature/cost before each channel in the output PQR file. Channel set and geometry unchanged.

channels: rename homogenize bool → diagram ("homogenized"/"simple"/"w…

64271bd

…eighted")

channels: update of calcChannels parameters pre-optimized to resemble…

1c9f29d

… results with Caver3 defaults

channels: update of style to better align with styleguide.rst by rena…

63dfdd7

…ming function names to mixedCase and :param: to :arg:, line length limits

channels: cleaning unused imports

2cdbfa2

channels: additively weighted VD are run in serial to avoid complicat…

b9ca70d

…ions with preceived per-frame parallelism

briza81 marked this pull request as draft July 16, 2026 05:01

briza81 marked this pull request as ready for review July 16, 2026 05:02

jamesmkrieger reviewed Jul 16, 2026

View reviewed changes