README.md

vc_reader

Tools for working with Vesuvius Challenge data: discover scans and annotations, read and sample large CT volumes, and transform aligned image and label inputs on CPU or CUDA. The package is built for large multidimensional images and PyTorch training pipelines; its Zarr reader and transforms also work with your own data.

• Discover data: browse published metadata for samples, volumes, segments, ink labels and geometry, including scans at a requested resolution.
• Read and sample: fetch local, HTTP or S3 Zarr crops with bounded caches, batch reads, and sample rotated boxes or curved surfaces.
• Align targets: read CT with ink and supervision masks, sample registered labels into a volume grid, or project surface annotations into 3D.
• Process images and labels: apply CPU/CUDA geometry, intensity, noise and filter transforms to 2D or 3D inputs while keeping images and labels aligned. CudaLoader adds GPU transforms and prefetching to a PyTorch DataLoader.

Volume metadata

Opening a volume also gathers what is known about it, from the store's meta.json / metadata.json and the Open Data catalog, merged into one value per field (the store's own documents win). Remote stores do this on first open only; the result is kept with the volume in the disk cache. ZarrReader.open(url, refresh=True) reads both sources again and replaces it.

volume = ZarrReader.open(url, level=2, cache_dir=cache)
volume.voxel_size        # micrometers per voxel of the opened level
volume.voxel_size_zyx    # per axis, for non-cubic pyramids
volume.metadata          # level-0 voxel_size_um, energy_kev, detector_distance_mm, beamline,
                         # scan_name, acquisition_date, value_window, z_crop, sample_id,
                         # volume_id, scan_id, name, created_at, data_format, documents

Capability map

Check this map before writing code that reads, samples, interpolates, parses, caches or writes Vesuvius data, in any project in this repository. If a row fits, call it directly, with no local copy, wrapper or renamed equivalent. If a row nearly fits, extend it in C++ here, then bind it to Python. The map covers each area's main entry points, so also search include/vc_reader/ before concluding something is absent.

• C++ API: in namespace vcr, declared in include/vc_reader/<header>. Each row's header is shown in its C++ column.
• Implementation: relative to src/cpp/. Headers there are also included as vc_reader/native/.... .cu files are the CUDA kernels.
• Python: in src/vc_reader/, binding the same C++. _reader is the compiled extension module, vc_reader._reader.
• Tests: native tests are tests/cpp/test_native_<area>.cpp; Python tests are tests/test_*.py.
• Conventions: arrays and grids are C order and ZYX. Geometry and surface coordinates are XYZ. Surface vertices are valid or invalid; invalid vertices are never sampled.

Catalog and metadata

Need	C++ API	Implementation	Python
Load or cache the Open Data manifest; parse samples, scans, volumes, segments, artifacts, surfaces	catalog.hpp: catalog::loadCatalog, OpenDataManifest, OpenDataSample, OpenDataSegment, OpenDataArtifact	catalog.cpp, catalog_io.cpp	load, Catalog
Sample name aliases, public URL rewriting, cache-safe ids	catalog.hpp: catalog::sampleNameKey, communityAliases, resolveOpenDataUrl, openDataS3Url, safePathComponent	catalog.cpp	to_https, to_s3
Browse scrolls, volumes, segments and labels; query by resolution	browsing.hpp: catalog::Catalog, ScrollSource, Collection, Volume, Segment, Asset, atResolution	browsing.cpp, sources/resolution.cpp	ScrollSource, Collection, Volume, Segment, Asset
Volume-to-volume registration routes	registrations.hpp: catalog::RegistrationGraph, resolveVolumeTransform	registrations.cpp	ScrollSource, VolumeMapping
Affines: compose, invert, map normals, read transform.json, block-reduced level placement	affine.hpp: catalog::Affine	header only	Affine
Immutable source, storage, encoding and registration descriptors	sources.hpp: sources::Source, Storage, Encoding, Surface, SurfaceRegistration, VolumeRegistration	sources.cpp	sources.py
Volume facts from meta.json, metadata.json and the catalog: voxel size, energy, per-level spacing	zarr_metadata.hpp: zarr::storeMetadata, volumeMetadata, voxelSizeUm, levelVoxelSizeZYX	zarr/store_metadata.cpp	ZarrReader.metadata, .voxel_size, .voxel_size_zyx

Zarr, storage and caches

Need	C++ API	Implementation	Python
Open a Zarr v2/v3 array, OME-Zarr group or pyramid level (local, HTTP, S3, SFTP; online or offline)	zarr.hpp: zarr::openArray, openPyramidLevel, pyramidLevels, Group	zarr/open.cpp	reader/zarr_reader.py, reader/zarr_group.py: ZarrReader, open_zarr, ZarrGroup, open_group, open_node
Level placement in level-0 coordinates; a level ready for grid::sample	zarr.hpp: zarr::LevelPlacement, levelPlacement, declaredPlacement, gridLevel	zarr/open.cpp	ZarrGroup.placement
Parse .zarray / zarr.json	zarr_array_meta.hpp: zarr::parseArrayMeta, ArrayMeta	zarr/meta.cpp	reader/zarr_meta.py: ArrayMeta
Region reads, integer point reads, nearest and linear sampling, affine grids, cache-only sampling, single chunks	array.hpp: Array::read_many, read_points, sample_nearest, sample_linear, read_affine, sample_nearest_cached, sample_linear_cached, chunk, cached_chunk	array.hpp, chunk_io.hpp, zarr_storage.hpp, shard_index.hpp	ZarrReader
The same reads on the GPU; one device brick cache shared by several readers (e.g. a pyramid's levels)	cuda_reader.hpp: reader_cuda::Reader, reader_cuda::Pool	reader_cuda/reader.hpp, reader_cuda/kernels.cu	reader/cuda.py: CudaReader, CudaPool
Chunks computed on demand, read and sampled like a stored array, with a persistent cache	computed.hpp: computed::open	computed.cpp	_reader.open_computed
Credentials and access: anonymous, basic, signed S3, AWS profiles, bucket region	zarr.hpp: zarr::resolveAccess, s3Access, basicAccess, awsCredentials, bucketRegion, loadBasicAuth	zarr/access.cpp, http.hpp, sigv4.hpp	reader/auth.py
Service API requests: status and body of every final answer (4xx included), transport errors/429/5xx retried; http::Options::pathAsIs preserves dot segments in opaque identifier paths	http.hpp: http::Session::request	http_session.cpp, http.hpp (HttpClient::exchange)	—
List S3 prefixes; URL helpers	zarr.hpp: zarr::listChildren, listObjects, listPage, isRemote, httpUrl, s3Location, levelGroupLocation; s3_listing.hpp: s3::listUrl, parseListPage	zarr/access.cpp, s3_listing.cpp	ObjectStore, reader/auth.py
SFTP sessions, pooled fetches, directory listing; owned SSH forwards to hosted services	sftp.hpp: sftp::Session, Fetcher, Url, setAuthenticator, Forward, parseForwardUrl	sftp/	through any sftp:// location; SshForward
Prefill array chunks or download a remote prefix through the existing disk caches, with progress and cancellation	cache.hpp: cache::prefill, cache::downloadPrefix	cache.cpp	_reader.prefill, _reader.download_prefix, _reader.CacheCancellation
Decoded chunk memory cache	array.hpp: MemoryCache	memory_cache.hpp, lru.hpp	reader/memory_cache.py: MemoryCache
Disk caches: budgets, cache discovery, VC3D cache roots and leases, low-disk handling	zarr.hpp: zarr::openLayers, diskBudget, decodedCachePath, vc3dConfigDirectory, vc3dCacheRoot, selectedVc3dCacheRoot, vc3dVolumeId, createVc3dVolumeDir, Vc3dLease, LowDiskSpace	zarr/layers.cpp, disk_budget.cpp, disk_cache.hpp, write_behind.hpp	reader/disk_cache.py
Write arrays (v2/v3, shards, local or S3), attributes, groups, inventories	zarr_write.hpp: zarr::Writer, createGroup, removeNode, nodeChildren, readAttributes, writeAttributes, inventory, atomicWrite, compressBytes, decompressBytes, decompressZstdFrame	zarr/writer.cpp	reader/writer.py: ZarrWriter, create_group, zarr_inventory, atomic_write, compress_bytes, decompress_bytes, decompress_zstd_frame
Object store over a local directory or S3 prefix	zarr_write.hpp: zarr::Store	zarr/store.cpp	ObjectStore
Build pyramids and OME multiscales metadata	zarr_write.hpp: zarr::writePyramid, writePyramidLevel, createPyramidLevels, omeMultiscales	zarr/pyramid.cpp	write_pyramid, write_pyramid_level, create_pyramid_levels, ome_multiscales
Compress and decompress bytes	zarr_write.hpp: zarr::compressBytes, decompressBytes	zarr/writer.cpp, codecs.hpp, lzw.hpp	compress_bytes, decompress_bytes
Atomic and durable file writes; directory publication with rollback; snapshots	zarr.hpp: zarr::atomicWrite, Durability; publication.hpp: files::publishDirectory, replaceFile, recoverDirectory, snapshot	atomic_file.hpp, zarr/layers.cpp, publication.cpp	snapshot

Images and labelled sources

Need	C++ API	Implementation	Python
Read TIFF, PNG or JPEG, whole, by page or windowed; memory-map a TIFF	source_reader.hpp: sources::imageSource, readImage, readImagePage, tiffMapping, ImageFile	sources/image_file.cpp, sources/source_reader.cpp	images.py: read, open, memmap
Write TIFF, PNG or JPEG atomically; encode to memory	image_write.hpp: sources::writeTiff, TiffWriter, writeImage, encodePng, encodeJpeg	sources/image_write.cpp, sources/image_encode.cpp	images.py: write_tiff, TiffWriter; write_image, encode_png, encode_jpeg
Crops of any Zarr, image or stack source, with exact encodings and decode budgets	source_reader.hpp: sources::SourceReader, DecodedCache, EncodedFiles	sources/source_reader.cpp, sources/caches.cpp	reader/source_reader.py, reader/decoded_cache.py
Aligned multi-target reads in a reference grid; flattened image and label pairs	aligned_reader.hpp: sources::SampleReader, PairedReader	sources/aligned_reader.cpp	sampling.py: open_sample; paired.py: PairedReader

Surfaces and tifxyz

Need	C++ API	Implementation	Python
Read tifxyz: metadata, shape, mask, windows, whole surface; list a folder	tifxyz.hpp: tifxyz::read, Reader, readMetadata, parseMetadata, storedShape, readMask, list, vertexValid, fullExtent	tifxyz/tifxyz.cpp	tifxyz.py: read, read_metadata, stored_shape, read_mask, list_surfaces, full_extent
Write tifxyz with channels and rollback	tifxyz.hpp: tifxyz::write, WriteOptions, Channel	tifxyz/tifxyz.cpp	tifxyz.write
In-memory surface: bounds, retargeting to another level, positions, normals and tangents	tifxyz.hpp: tifxyz::Surface	tifxyz/tifxyz.cpp	tifxyz.Surface
Bilinear or bicubic Hermite interpolation of a vertex grid, with normals and tangents	grid.hpp: grid::Surface, resample, pointsAt, snapshot, Raster, Interpolation	grid/grid.cpp, grid/points.hpp	tifxyz.Surface.resample, .points_at
Sample a volume along a surface or plane: layers along the normal, pyramid fallback, cache-only reads, masks; scalar CZYX channel views (CPU and uint8 CUDA, shared packed chunk cache)	grid.hpp: grid::sample, Request, Level (channel for CZYX), Plane, Points, levelFromMap	grid/grid.cpp, grid/level_map.hpp, grid/kernels.cu	grid.sample_grid (every request field); GridLevel.channel; ZarrReader.cuda(channel=…); SampleReader.render_surface; prepared per-pixel XYZ/directions: ZarrReader.render_points
Reduce layers to one value per pixel: mean, max, min, median, min-abs, alpha, Beer-Lambert	composite.hpp: composite::reduce	composite/	_reader.composite_reduce
Render surface layers; project a label patch into CT	aligned_reader.hpp: sources::SampleReader, SurfaceProjection	sources/surface_reader.cpp, projection/	SampleReader.render_surface, .project_surface, .iter_project_surface

Resident arrays and surface transfer

Need	C++ API	Implementation	Python
Nearest/linear samples of resident batched 2D/3D float32 arrays, with value and coordinate gradients; CPU/CUDA	resident_sampling.hpp: sampling::Dense, sampleDense, backwardDense	sample_point.hpp, sampling_math.hpp, sampling/	resident.sample_array, torch.sample_array
Batched bicubic surface positions/tangents, linear or triangle fields, explicit validity; CPU/CUDA	resident_sampling.hpp: sampling::prepareSurfaces, sampleSurfaces	grid/derivatives.hpp, grid/points.hpp, sampling/	resident.GridField, torch.surface_derivatives, torch.surface_points
Pixel-centre field tiles with border replication; CPU/CUDA	resident_sampling.hpp: sampling::Tile, sampleFieldTile, sampleSurfaceTile	sampling_math.hpp, sampling/	resident.sample_field_tile
Label lookup through a surface correspondence, distance rejection and seam provenance; CPU/CUDA	sampling.hpp: projection::RasterRequest, transfer	projection/transfer.cpp, projection/kernels.cu	resident.SurfaceTransfer
Binary support swept through a surface normal band; CPU/CUDA	sampling.hpp: projection::MaskBand, projectMask	projection/common.hpp, grid/points.hpp, projection/mask.cpp, projection/kernels.cu	resident.project_mask

Resident arrays are already loaded host or device buffers. Their forward samples use the same voxel-lookup sampler as stored arrays; the adapter supplies a contiguous block, batch/channel layout, and boundary policy. Queries use array indices (ZYX or YX), without normalized Torch coordinates. zeros interpolates with the fill value outside the array, border replicates its boundary, and constant rejects out-of-bounds queries. Nearest ties round up. Coordinate gradients use the shared corner weights, including all axes at integer queries. Field tiles and correspondence UV maps share the clamped array interpolation used by transforms. Surface queries use row/column indices and require valid corners; bicubic Hermite derivatives are shared with tifxyz.Surface. GridField copies its input into fixed storage; CUDA copies remain resident. Torch sampling supports first-order gradients; CUDA value gradients use atomic addition and may differ slightly between runs. MPS tensors stage through CPU. Surface supervision queries do not support autograd.

Surface transfer uses bicubic target geometry, triangle source geometry matching the correspondence mesh, and linear UV maps to avoid overshooting seams. Categorical labels are selected by index. Surface-to-volume projection prefers nonzero labels over zero; between nonzero labels, the nearest surface sample wins, with sample ID breaking ties. A multichannel label is nonzero when any channel is nonzero, and all channels come from the same winning sample. Binary masks use the same surface interpolation, sampling-density policy, endpoint-inclusive normal layers (at most half a voxel apart), and nearest-voxel placement as label projection. Supplied normals allow isolated valid points and edges to contribute too. Binary occupancy accumulates by union using deterministic integer writes. All these operations use float32 on CPU and CUDA.

Geometry and analysis

Need	C++ API	Implementation	Python
Polyline queries (nearest, interpolate, plane and box hits) and voxel-set membership	sampling.hpp: geometry::Curves, Voxels	geometry/	geometry.py: PolylineSet, VoxelSet, transform_points
Low-level aligned and projection sampling plans	sampling.hpp: aligned::Plan, PointPlan, projection::Plan	aligned/, projection/	through SampleReader
Triangle slice labels, interpolated normals and surface frames (ordered DDA rasterization)	sampling.hpp: geometry::MeshSlice, raster_slice	geometry/triangle_mesh.hpp, geometry/kernels.cu	geometry.rasterize_mesh_slice
Sheet normals from CT	reader.hpp: sheetNormals, sheetNormalsCpu; cuda_reader.hpp: sheet_normal::sample_cuda	sheet_normal/, native_reader.cpp, symmetric_eigen.hpp	sheet_normal.py: sheet_normals, sample_sheet_normals
Shortest path over a 2D cost grid	pathfinding.hpp: pathfinding::shortestPath	pathfinding/	shortest_path
Connected components of a 2D/3D mask (4/8, 6/18/26-connected), numbered as scipy.ndimage.label; CPU (threaded, run-based) or CUDA	components.hpp: components::label, components::cuda::label	components/	geometry.label_components; transforms.ConnectedComponents
Fill holes of a 2D/3D mask (scipy.ndimage.binary_fill_holes), CPU or CUDA	components.hpp: components::fillHoles, components::cuda::fillHoles	components/	geometry.fill_holes
Binary YX skeletons (Zhang–Suen, Guo–Hall, combined), independent slice batches, host/device buffers	thinning.hpp: thinning::skeleton, cuda::skeleton, cuda::skeletonDevice	thinning/skeleton.cpp, thinning/kernels.cu	thinning.skeletonize, skeletonize_into; transforms.MedialSurface
Ordered distance-ridge traces, pruning, reusable CPU scratch; complete device-side tracing	thinning.hpp: thinning::trace, prune, cuda::traceDevice	thinning/trace.cpp, thinning/trace_cuda.cu	thinning.trace_ridges, trace_ridges_into, TraceWorkspace; transforms.RidgeTrace
Exact uint8 Manhattan distance to background in 2D/3D, saturated at 255	distance_transform.hpp: distance::manhattan, cuda::manhattanDevice	distance/	thinning.manhattan_distance

Display

Need	C++ API	Implementation	Python
Float to uint8/uint16 conversion, window/level, palette and overlay LUTs, ARGB blending	display.hpp: display::toUint8, toUint16, windowLevelLut, paletteLut, overlayLut, applyLut, blend, blendArgb	display/	display_to_uint8, display_to_uint16, window_level_lut, palette_lut, overlay_lut, apply_lut, blend, blend_argb
Volumetric slab compositing, including perspective	slab.hpp: display::compositeSlab, opacityLut, Slab	display/slab.cpp, display/slab_pixel.hpp	composite_slab, opacity_lut, SlabProjection
Views: the level for a zoom, fallback levels, a surface raster's pixels per voxel; progressive frames drawn from the caches with the reads they lacked	grid.hpp: grid::levelForView, fallbackLevels, voxelExtent, pixelsPerVoxel, Frame, Pending, sampleCached, load	grid/view.cpp	grid.py: level_for_view, fallback_levels, pixels_per_voxel, Pending, sample_grid, sample_grid_cached, load; levels from ZarrReader.grid_level
The render engine from C or any language with a C FFI: volumes, surfaces, frames, progressive reads, composites, display tables	vc_reader_c.h (vc_reader::c, libvc_reader_c)	c_api/c_api.cpp	—

Transforms and training

Need	C++ API	Implementation	Python
Augmentation and processing pipelines with seeded replay, on CPU or CUDA (catalog)	transforms.hpp: transforms::Pipeline, Op, Sample, Plan, save_plan, load_plan	transforms/pipeline.hpp, transforms/ops/, transforms/cpu/, transforms/cuda/	vc_reader.transforms
Device buffers, streams, resident pipelines and batches; how much freed memory the default CUDA pool keeps (all, by default)	cuda.hpp: transforms::cuda::begin_device, run_device, run_batch; cuda_available.hpp: cuda::available	cuda_buffer.hpp, cuda_stream.hpp, cuda_pool.hpp, transforms/cuda/device_run.hpp	cuda.py: DeviceArray, is_available, asarray, wrap, set_pool_release_threshold, pool_release_threshold
PyTorch GPU transforms and prefetching; masked normalization; differentiable structure tensors	transforms.hpp: StructureTensor (component cotangent via gradient)	transforms/	torch.py: CudaLoader, structure_tensor; direct calls A.Normalize(...).batch
Pavel/Gaussian structure components, directions; Euclidean distances; Hessian ridge responses	transforms.hpp: StructureTensor, DistanceTransform, Ridges	transforms/cpu/, transforms/cuda/	transforms.StructureTensor, DistanceTransform, Ridges
Coherence-enhancing diffusion of a single image or volume, on CPU or CUDA	transforms.hpp: CoherenceEnhancingDiffusion	transforms/cpu/coherence_enhancing_diffusion.hpp, transforms/cuda/coherence_enhancing_diffusion.cu	transforms.CoherenceEnhancingDiffusion

Shared primitives

Use these instead of writing another copy. They are header-only under src/cpp/ and included as vc_reader/native/<file>. They have no stability guarantee beyond this repository.

Need	File
Tie rules for rounding coordinates and values	rounding.hpp
Bilinear, trilinear and masked interpolation; corner weights	sampling_math.hpp, sample_point.hpp
Numpy typestrs and element types; encoding element values	element_type.hpp, element_value.hpp
Symmetric eigendecomposition	symmetric_eigen.hpp
LRU storage, thread pool, reused output buffers	lru.hpp, thread_pool.hpp, output_pool.hpp
Atomic file writes	atomic_file.hpp
Portable file I/O (open, pread, fstat, mmap, file locks, rename-over, temp dirs) on POSIX and Windows; headers must not include <windows.h>	platform_io.hpp (Windows bodies in platform_io.cpp)
OME storage transforms	storage_transform.hpp
Percentile, median, histogram, Fourier, elastic field and resize helpers shared by CPU and CUDA transforms	transforms/shared/

Guides: C++, discovery, source reading, aligned sampling, surface sampling, precision.

Quick start

See installation and notebook setup. This example discovers two image/ink-label pairs on PHerc0009B and PHercParis4, then trains one tiny ink classifier on two automatically selected tiles from each pair. It requires PyTorch, a CUDA-enabled build (selected automatically when a CUDA compiler is found), and an accessible CUDA device. Reads use the CUDA backend; DLPack shares the returned GPU buffers with PyTorch without downloading them. A small CPU scan of the label masks finds annotated tiles before any CT pixels are fetched. The 25 surface layers become input channels for a 2D CNN with batch size 1.

A tile is a rectangular crop of the flattened surface image. Its (y, x) is the top-left row and column in the selected image pyramid level, shared by the paired image, label and masks. These are not coordinates in the original 3D scan. The example walks a grid of tiles and keeps the first two containing supervised training pixels for each pair.

from itertools import product

import numpy as np
import torch
from torch import nn
from torch.nn import functional as F
from vc_reader import load, ScrollSource

catalog = load()
selections = [  # Sample and segment IDs from the catalog.
    ("PHerc0009B", "20250919125754"),
    ("PHercParis4", "20231210121321"),
]

# discover labeled regions using a lower resolution
LEVEL, TILE_SIZE = 3, 128 
pairs = []
for sample_id, segment_id in selections:
    scroll = ScrollSource(sample_id, catalog=catalog)
    label = scroll.segments[segment_id].labels.ink(label_set="20260918")
    surface_volume = label.surface_volumes.zarr()
    tiles = []
    # An empty CT layer slice reads only labels and masks while locating tiles.
    with label.open_paired(image=surface_volume, level=LEVEL, layers=np.s_[0:0]) as overview:
        height, width = overview.shape
        for y, x in product(range(0, height, TILE_SIZE), range(0, width, TILE_SIZE)):
            patch = overview[y:y + TILE_SIZE, x:x + TILE_SIZE]
            support = patch.label.valid & (patch.label.masks["supervision"] > 0)
            if "validation" in patch.label.masks:
                support &= patch.label.masks["validation"] == 0
            if support.any():
                tiles.append((y, x))  # Top-left image pixels, found from annotations.
                if len(tiles) == 2:
                    break
    if not tiles:
        raise ValueError(f"No supervised training tiles found for {sample_id}/{segment_id}")
    pairs.append((label, surface_volume, tiles))

# construct a tiny ex model
model = nn.Sequential(
    nn.Conv2d(25, 16, kernel_size=3, padding=1),
    nn.ReLU(),
    nn.Conv2d(16, 1, kernel_size=1),
).cuda().train()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)
steps = 0

readers = [
    (label.open_paired(
        image=surface_volume, level=LEVEL, layers=np.s_[40:65], backend="cuda"), tiles)
    for label, surface_volume, tiles in pairs
]

# train 
for epoch in range(3):
    for reader, tiles in readers:
        for y, x in tiles:
            patch = reader[y:y + TILE_SIZE, x:x + TILE_SIZE]
            image = torch.from_dlpack(patch.image.data).float() / 255  # (1, 25, H, W)
            ink = (torch.from_dlpack(patch.label.data) > 0).float()    # (1, H, W)
            image_valid = torch.from_dlpack(patch.image.valid).bool().all(dim=1)
            label_valid = torch.from_dlpack(patch.label.valid).bool()
            supervision = torch.from_dlpack(patch.label.masks["supervision"]) > 0
            train_mask = image_valid & label_valid & supervision
            if "validation" in patch.label.masks:
                train_mask &= torch.from_dlpack(patch.label.masks["validation"]) == 0
            if not train_mask.any():
                continue

            optimizer.zero_grad(set_to_none=True)
            pixel_loss = F.binary_cross_entropy_with_logits(model(image)[0], ink, reduction="none")
            loss = (pixel_loss * train_mask).sum() / train_mask.sum()
            loss.backward()
            optimizer.step()
            steps += 1

for reader, _ in readers:
    reader.close()

torch.save({
    "model_state_dict": model.state_dict(),
    "optimizer_state_dict": optimizer.state_dict(),
    "epochs_completed": epoch + 1,
    "steps": steps,
}, "ink_checkpoint.pt")
print(f"Completed {steps} optimizer updates; saved ink_checkpoint.pt")

Only supervised pixels with valid labels and valid CT in every input layer contribute to the loss. Published validation masks are excluded when present; PHercParis4 has no published validation mask in this selection. Readers stay open across epochs; call close() when training is finished to release their caches and resources. For full training and evaluation with a spatial holdout, see the 2.5D ink training notebook.

Each surface_volume is CT sampled along its segment's curved surface; its depth axis indexes surface layers. label.volume identifies the original 3D scan, which can be read in ZYX voxel coordinates through label.volume.open(level=0).

.surface_volumes.zarr() requires a unique matching render; multiple matches raise an error. Use filters such as artifact_id=... to select among them.

patch keeps each target's data, geometric validity and companion masks together.

Performance

Zarr reading and sampling:

Rows are sorted fastest first by aggregate Mvox/s across all workloads (total output voxels / total measured seconds). Each cell shows total seconds, MB/s, then Mvox/s.

Each workload covers at least 5,000 distinct 128³ storage chunks in nonzero regions of PHercParis4 level 2. Measured 2026-10-19 with eight reader threads. Total seconds sum all measured request durations across five repetitions after one warmup pass. Host decoded caches are disabled; CUDA chunks are cleared per request and host downloads are timed.

zarr-python and tensorstore used scipy for trilinear so the most direct raw read speed comps are the chunk based ones (batch and slab)

Reader	Batch, 8 × 128³	Plane, 512²	Surface, 96 × 512	Slab, 16 × 512²	Rotated box, 128³
vc_reader CPU	2.738 s 13,539.4 MB/s 13,539.4 Mvox/s	3.631 s 557.3 MB/s 139.3 Mvox/s	3.091 s 249.4 MB/s 62.3 Mvox/s	0.910 s 9,283.8 MB/s 9,283.8 Mvox/s	12.612 s 2,700.4 MB/s 675.1 Mvox/s
vc_reader CUDA	11.738 s 3,158.8 MB/s 3,158.8 Mvox/s	3.578 s 565.6 MB/s 141.4 Mvox/s	7.342 s 105.0 MB/s 26.2 Mvox/s	4.504 s 1,876.4 MB/s 1,876.4 Mvox/s	6.853 s 4,969.8 MB/s 1,242.4 Mvox/s
TensorStore	13.256 s 2,796.9 MB/s 2,796.9 Mvox/s	34.737 s 58.3 MB/s 14.6 Mvox/s	27.506 s 28.0 MB/s 7.0 Mvox/s	5.968 s 1,416.2 MB/s 1,416.2 Mvox/s	263.496 s 129.3 MB/s 32.3 Mvox/s
zarr-python	32.237 s 1,150.2 MB/s 1,150.2 Mvox/s	59.458 s 34.0 MB/s 8.5 Mvox/s	47.532 s 16.2 MB/s 4.1 Mvox/s	12.356 s 684.0 MB/s 684.0 Mvox/s	275.765 s 123.5 MB/s 30.9 Mvox/s

MB/s counts decimal output megabytes; Mvox/s counts output voxels. Box/slab outputs are uint8, and interpolated outputs are float32. Throughput uses the same total seconds shown above. See benchmark details for median latencies, coverage, request counts, and reproduction settings.

Read + normalize + 12 transforms:

Rows are sorted fastest first by aggregate Mvox/s across all workloads. Each cell shows total seconds, MB/s, then Mvox/s.

Reader	Batch, 8 × 128³	Plane, 512²	Surface, 96 × 512	Slab, 16 × 512²	Rotated box, 128³
vc_reader CUDA batched	4.818 s 6,156.1 MB/s 1,539.0 Mvox/s	0.747 s 542.1 MB/s 135.5 Mvox/s	1.585 s 97.3 MB/s 24.3 Mvox/s	2.074 s 3,260.3 MB/s 815.1 Mvox/s	1.687 s 4,038.8 MB/s 1,009.7 Mvox/s
vc_reader CUDA unbatched	5.376 s 5,517.4 MB/s 1,379.3 Mvox/s	0.759 s 533.0 MB/s 133.2 Mvox/s	1.565 s 98.5 MB/s 24.6 Mvox/s	2.147 s 3,149.5 MB/s 787.4 Mvox/s	1.666 s 4,089.2 MB/s 1,022.3 Mvox/s
vc_reader CPU	515.514 s 57.5 MB/s 14.4 Mvox/s	8.563 s 47.3 MB/s 11.8 Mvox/s	3.989 s 38.6 MB/s 9.7 Mvox/s	134.556 s 50.2 MB/s 12.6 Mvox/s	120.389 s 56.6 MB/s 14.1 Mvox/s

Uses the exact same accepted requests and 5,000+ chunk coverage as the reader comparison. Measured 2026-10-19: one complete timed pass after one complete warmup pass per configuration. All outputs are float32; MB/s is four times Mvox/s. Timing includes reading, normalization, all twelve transforms, device batch assembly where applicable, and CUDA host downloads. Both CUDA modes use reused pinned host buffers prepared during warmup. Batched CUDA uses Pipeline.batch() per request: eight crops for ML and one otherwise. Unbatched CUDA and CPU apply transforms per crop. CPU measurements are retained from the matching full-coverage CPU run. See pipeline benchmark details for validation and reproduction settings, and the focused batching comparison for transform-only timings and crop-size scaling.

Examples and guides

The notebooks cover 2.5D and 3D ink training, general Zarr reading, CUDA transform, and mesh-segment rendering.

Topic	Guide
Metadata, resolution queries and custom data	Discovery and pairing
Source formats, storage and caches	Source reading
Registered labels and fields	Aligned sampling
Surface layers and projecting ink into CT	Surface sampling
Measurements and reproduction commands	Benchmarks

Transforms

Image transforms and normalization support channel-first 2D (C, Y, X) and 3D (C, Z, Y, X) arrays on CPU and CUDA. Use from vc_reader import transforms as A. Expand an entry for its measured configuration, timings and previews; the recipes are examples, not constructor defaults. Z/Y/X slice pairs show the original above the transformed output.

Use fixed parameters for processing, such as A.GaussianBlur(sigma=1.0), or sample parameters with A.Uniform, A.Normal and A.Choice for augmentation. Several constructors have random defaults; pipeline calls take an explicit seed for reproducibility, and Plan records sampled operations for replay.

Every transform can also be called on its own. A transform is random only in a Pipeline: called directly it always applies (p is ignored) and draws anything else from a fixed seed, so the same input gives the same output. Pass plain values for processing. An array (NumPy, Torch CPU/CUDA or DeviceArray) stands for the field the transform reads by default and the call returns the field it creates, else the transformed input; a sample dict works like a pipeline call. .batch does the same over (B, C, ...) batches, CUDA ones in one launch.

sdt = A.DistanceTransform(mode="outside")(mask)        # (1, Z, Y, X) label
st = A.StructureTensor(sigma=1, rho=1.5)(image)          # created field
z = A.Normalize("zscore", nonzero=False).batch(crops)    # (B, 1, Z, Y, X) CUDA
out = A.Normalize("minmax", masked=True).batch(
    {"image": x, "valid": v}, targets={"image": A.Image(valid="valid", refill=False), "valid": A.ValidMask()})

Operation family	Examples
Spatial	Affine, Flip, CropPad, Elastic, ForegroundCrop
Intensity	Normalize, Gamma, CLAHE, HistogramMatching
Filters	GaussianBlur, MedianFilter, OrientedBlur, UnsharpMask
Noise and acquisition effects	GaussianNoise, BiasField, GibbsRinging, ThickSlice
Masking and dropout	MaskImage, GridDropout, ObjectDropout
Labels and morphology	Morphology, ReplaceLabel, Regions, OneHot
Features	DistanceTransform, StructureTensor, MedialSurface
Layout and output construction	ToPseudo2D, FromPseudo2D, DeepSupervision

Gallery timings use 128³ float32 CT crops, one CPU thread, and medians of five calls after one warmup. CUDA uploaded includes upload and download; CUDA resident starts and ends on the GPU. Reading is excluded.

AdditiveNoise

Add independent noise to each voxel. Supports uniform, Gaussian, Laplace and centered beta distributions; this example uses Laplace noise.

A.AdditiveNoise(scale=0.05, distribution='laplace', mean=0)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
15.608 ms	1.313 ms	0.240 ms	11.9×	65.1×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Affine

Rotate, translate and scale the volume, applying the same geometry to images and labels.

A.Affine(degrees=5, translate=2, scale=0.95)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
68.183 ms	6.039 ms	0.608 ms	11.3×	112.1×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

BiasField

Multiply intensities by a smooth 3D gain field generated from a coarse Gaussian random lattice.

A.BiasField(std=0.3, scale=0.025)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
11.655 ms	1.442 ms	0.323 ms	8.1×	36.1×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

BlankRectangle

Fill randomly placed cuboids with a chosen value, simulating missing regions.

A.BlankRectangle(size=16, value=0, count=3)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
0.797 ms	1.205 ms	0.111 ms	0.7×	7.2×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

BrightnessAdd

Add a constant offset to every voxel in the selected image channels.

A.BrightnessAdd(shift=0.1)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
1.650 ms	1.329 ms	0.158 ms	1.2×	10.4×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

BrightnessGradient

Add a localized Gaussian brightness field. This example fixes its center at the middle of the volume; loc=(0, 1) samples a random position. The field is mean-centered by default.

A.BrightnessGradient(scale=12, loc=0.5, strength=0.3)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
2.587 ms	1.380 ms	0.184 ms	1.9×	14.1×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

BrightnessMultiply

Multiply every voxel in the selected image channels by a brightness factor.

A.BrightnessMultiply(factor=1.1)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
1.632 ms	1.359 ms	0.119 ms	1.2×	13.7×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

ChannelMisalignment

Independently shift, rotate and scale selected image channels to simulate registration errors. This example moves the second of two identical CT channels.

A.ChannelMisalignment(degrees=5, translate=2, scale=1, channels=(1,))

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
46.855 ms	2.591 ms	0.252 ms	18.1×	186.2×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Two input channels (second is a copy of the CT); only channel 1 moves. Both channels are included in timing.

CLAHE

Equalize local histograms in 3D tiles, limit their contrast, and interpolate between neighboring tiles to avoid seams.

A.CLAHE(clip_limit=2, tile_shape=(16, 16, 16), bins=256)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
52.663 ms	1.712 ms	0.549 ms	30.8×	96.0×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

ClipOutliers

Clip intensities at selected percentiles. This example then restores the original standard deviation.

A.ClipOutliers(lower=1, upper=99, p_retain_std=1)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
21.853 ms	1.757 ms	0.532 ms	12.4×	41.1×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Contrast

Scale intensity differences from the channel mean. The default range-preservation option clips results to the original range.

A.Contrast(factor=1.2)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
5.635 ms	1.297 ms	0.133 ms	4.3×	42.3×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

CropPad

Crop or pad all aligned targets to a requested ZYX shape. This example takes a centered 112³ crop from the 128³ input.

A.CropPad(shape=(112, 112, 112))

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
0.911 ms	1.932 ms	0.369 ms	0.5×	2.5×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Defocus

Average intensities within a spherical 3D neighborhood to simulate out-of-focus acquisition.

A.Defocus(radius=2)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
9.619 ms	1.944 ms	0.776 ms	4.9×	12.4×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Elastic

Warp the volume with a smooth random 3D B-spline displacement field, keeping images and labels aligned. The previews use interior slices so boundary padding does not dominate the view.

A.Elastic(sigma=32, magnitude=8, border=True)

CPU and CUDA construct and store elastic fields in float32. Ordinary crop geometry also uses float32, with wider arithmetic where coordinate range or numerical conditioning requires it. Output dtypes are preserved. Small coordinate differences can select a neighboring label at a nearest-neighbor boundary. Replay plans store float32 fields. See the precision policy and transform benchmarks.

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
85.995 ms	2.536 ms	0.696 ms	33.9×	123.6×

▶ Watch Elastic animation (WebM)

Strong cubic B-spline deformation: original → deformed → original. Fixed Z/Y/X slices, seed 125, achieved RMS 5.36 voxels.

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Strong cubic B-spline example: 64-voxel knot spacing, requested RMS displacement 8 voxels, with border clamping.

Equalize

Equalize the histogram across the entire volume while retaining the channel intensity range.

A.Equalize(bins=256)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
19.669 ms	1.588 ms	0.369 ms	12.4×	53.4×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

FDA

Blend low-frequency amplitudes from a reference volume into the 3D Fourier spectrum while preserving the source phase.

A.FDA(reference=reference, blend=1, fraction=0.1)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
613.754 ms	3.463 ms	2.119 ms	177.2×	289.6×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Flip

Mirror the volume along selected axes. This example always flips X.

A.Flip(axes='x', p_axis=1)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
1.425 ms	1.869 ms	0.323 ms	0.8×	4.4×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

ForegroundCrop

Choose a random crop containing a sampled foreground voxel. All aligned targets follow the crop; empty masks fall back to random cropping.

A.ForegroundCrop(shape=(112, 112, 112))

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
2.387 ms	2.059 ms	0.531 ms	1.2×	4.5×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

ForegroundDropout

Place cuboids around sampled foreground voxels and erase only foreground inside those cuboids.

A.ForegroundDropout(size=(16, 16, 16), count=3)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
6.501 ms	1.885 ms	0.452 ms	3.4×	14.4×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Gamma

Apply a power-law intensity curve over the channel range, changing the relative brightness of dark and light structures.

A.Gamma(gamma=0.8)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
8.660 ms	1.420 ms	0.133 ms	6.1×	65.1×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

GaussianBlur

Smooth the volume with a separable 3D Gaussian filter.

full_width=False, truncate=4 uses radius floor(4*sigma + 0.5), matching SciPy Gaussian filter support. The default full_width=True interprets truncate as the full width. CPU and CUDA use the same convention.

A.GaussianBlur(sigma=1)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
6.254 ms	1.412 ms	0.229 ms	4.4×	27.4×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

GaussianNoise

Add independent zero-mean Gaussian noise with the specified standard deviation.

A.GaussianNoise(std=0.05)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
9.104 ms	1.255 ms	0.125 ms	7.3×	73.0×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

GibbsRinging

Truncate high frequencies in the 3D Fourier spectrum, creating ringing around sharp boundaries.

A.GibbsRinging(cutoff=0.7)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
425.044 ms	1.617 ms	0.485 ms	262.9×	875.9×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

GridDropout

Erase a periodic grid of cuboid holes with a randomly sampled offset.

A.GridDropout(unit_size=(16, 16, 16), ratio=0.4)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
39.059 ms	1.297 ms	0.146 ms	30.1×	267.9×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

HistogramMatching

Match the volume intensity histogram to a reference volume, with optional blending.

A.HistogramMatching(reference=reference, blend=1, bins=256)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
86.937 ms	3.000 ms	1.899 ms	29.0×	45.8×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Illumination

Apply a smooth multiplicative illumination field. Supports directional ramps, corner falloff and Gaussian illumination; this example uses a diagonal ramp.

A.Illumination(strength=0.3, mode='linear', direction=(1, 1, 1))

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
28.240 ms	1.355 ms	0.201 ms	20.8×	140.3×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Invert

Reflect intensities around the channel mean, reversing bright and dark structures.

A.Invert()

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
2.582 ms	1.389 ms	0.178 ms	1.9×	14.5×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

LocalContrast

Change contrast around a Gaussian-weighted local mean inside a smooth region.

A.LocalContrast(scale=12, loc=0.5, factor=1.2)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
4.167 ms	1.297 ms	0.171 ms	3.2×	24.4×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

LocalGamma

Blend a gamma-adjusted image into a localized Gaussian region.

A.LocalGamma(scale=12, loc=0.5, gamma=0.8)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
12.362 ms	1.308 ms	0.161 ms	9.5×	76.6×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

LocalSmoothing

Blend Gaussian smoothing into a localized region, leaving distant voxels largely unchanged.

A.LocalSmoothing(scale=12, loc=0.5, strength=0.7, sigma=1)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
7.336 ms	1.695 ms	0.414 ms	4.3×	17.7×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

LowResolution

Downsample and upsample the image to simulate reduced spatial resolution.

A.LowResolution(scale=0.5)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
6.298 ms	1.284 ms	0.127 ms	4.9×	49.6×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

MedianFilter

Replace each voxel with the median of its 3D neighborhood, suppressing isolated noise while retaining edges.

A.MedianFilter(size=3)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
15.783 ms	1.310 ms	0.149 ms	12.1×	105.9×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Morphology

Apply a morphological operation to labels. This example closes small gaps using dilation followed by erosion. footprint is ball (Euclidean radius), cross (L1 radius: SciPy's default structure with iterations=radius), cube (Chebyshev radius: SciPy's full-connectivity structure with iterations=radius) or repeated_cross. Outside the array counts as background, as SciPy's border_value=0; neutral_border=True keeps erosion from eating in from the edges.

A.Morphology(operation='close', radius=1)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
12.452 ms	0.732 ms	0.312 ms	17.0×	39.9×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

MotionBlur

Average samples along a 3D line to simulate motion during acquisition.

A.MotionBlur(length=5, direction=(1, 1, 1), samples=9)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
240.809 ms	2.594 ms	1.356 ms	92.8×	177.6×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

MultiplicativeNoise

Multiply each voxel by an independently sampled uniform random gain.

A.MultiplicativeNoise(low=0.8, high=1.2)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
8.368 ms	1.310 ms	0.149 ms	6.4×	56.3×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Normalize

Normalize intensities per channel. This example uses z-score normalization of nonzero voxels, keeping zeros unchanged; its display window includes negative values.

A.Normalize(mode='zscore')

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
6.726 ms	1.347 ms	0.245 ms	5.0×	27.5×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

ObjectDropout

Randomly remove whole connected foreground objects, clearing their image voxels and labels by default.

A.ObjectDropout(drop_probability=0.5)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
13.619 ms	2.345 ms	0.509 ms	5.8×	26.8×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

OrientedBlur

Smooth along a specified 3D direction and two perpendicular directions with separate blur widths. This example is aligned with Z.

A.OrientedBlur(sigma_along=2, sigma_across=0.5, axis=(1, 0, 0))

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
7.649 ms	1.874 ms	0.719 ms	4.1×	10.6×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

PixelDistributionAdaptation

Transfer reference intensity statistics. The standard mode shown here blends the reference mean and standard deviation; minmax mode transfers the range.

A.PixelDistributionAdaptation(reference=reference, blend=1, mode='standard')

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
18.454 ms	2.614 ms	1.448 ms	7.1×	12.7×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

PlasmaBrightnessContrast

Use a smooth multiscale 3D random field to vary brightness and contrast throughout the volume.

A.PlasmaBrightnessContrast(brightness=0.2, contrast=0.3, scale=0.025, persistence=0.5, octaves=5)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
31.285 ms	1.831 ms	0.677 ms	17.1×	46.2×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

PlasmaShadow

Darken the volume with a smooth multiscale 3D random shadow field.

A.PlasmaShadow(strength=0.5, scale=0.025, persistence=0.5, octaves=5)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
25.357 ms	1.690 ms	0.591 ms	15.0×	42.9×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

RandomBinaryOperator

Apply a selected binary morphological operation to image channels. This example dilates a thresholded foreground volume.

A.RandomBinaryOperator(channels=(0,), radius=1, operations=('dilate',))

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
10.243 ms	1.404 ms	0.241 ms	7.3×	42.6×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

RandomToneCurve

Apply a random monotone nonlinear intensity curve while preserving the channel range endpoints.

A.RandomToneCurve(scale=0.1)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
13.061 ms	1.232 ms	0.194 ms	10.6×	67.4×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

RemoveComponents

Randomly remove connected foreground components below a size threshold. This example modifies the label mask.

A.RemoveComponents(on='mask', max_fraction=0.25, p_fill=0)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
20.792 ms	1.335 ms	0.970 ms	15.6×	21.4×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

RicianNoise

Combine the image with two independent Gaussian noise fields to simulate magnitude-image Rician noise.

A.RicianNoise(std=0.05)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
21.196 ms	1.422 ms	0.185 ms	14.9×	114.8×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

RingingOvershoot

Apply a separable windowed-sinc filter in 3D, retaining edge overshoot and ringing.

A.RingingOvershoot(cutoff=0.7, radius=4)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
7.994 ms	5.418 ms	4.292 ms	1.5×	1.9×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Rotate90

Rotate aligned targets by multiples of 90 degrees. This example rotates the YX plane by one quarter turn.

A.Rotate90(planes=('yx',), turns=(1,))

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
2.052 ms	2.296 ms	0.587 ms	0.9×	3.5×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

SaltAndPepper

Replace selected voxels with the channel minimum or maximum, producing dark and bright impulse noise.

A.SaltAndPepper(amount=0.02, salt_fraction=0.5)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
9.965 ms	1.250 ms	0.206 ms	8.0×	48.3×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Sharpen

Enhance local edges with a discrete Laplacian sharpening filter. This example keeps overshoots instead of clipping them.

A.Sharpen(strength=0.25, p_clamp=0)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
4.557 ms	1.449 ms	0.174 ms	3.1×	26.3×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

SheetCompression

Use the CT intensity profile to compress sheet-like structures along a selected axis while transforming aligned targets together.

A.SheetCompression(strength=0.2, intensity_sigma=1, smoothing=2, axes='z')

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
41.401 ms	2.545 ms	0.887 ms	16.3×	46.7×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

ShotNoise

Sample Poisson photon counts from nonnegative intensities and rescale them, producing signal-dependent shot noise.

A.ShotNoise(photons=100)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
61.419 ms	2.374 ms	1.268 ms	25.9×	48.4×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

SliceIllumination

Reduce brightness across selected slices and wider bands to simulate slice-dependent acquisition defects.

A.SliceIllumination(defects=3, width=4, reduction=0.5, base_p=0.1, base_reduction=0.5)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
1.115 ms	1.277 ms	0.156 ms	0.9×	7.2×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Smear

Blend shifted copies of preceding slices into each slice to simulate directional smearing.

A.Smear(shift=(3, 0), alpha=0.25)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
5.003 ms	1.556 ms	0.461 ms	3.2×	10.8×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

StripeMasking

Erase periodic slabs along selected axes with random offsets; overlapping slabs form a stripe pattern through the volume.

A.StripeMasking(unit_size=(16, 16, 16), ratio=0.25, axes='zyx')

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
39.640 ms	1.319 ms	0.140 ms	30.1×	283.3×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

ThickSlice

Downsample and upsample along selected axes to simulate thicker acquisition slices. Labels follow the same slice mapping.

A.ThickSlice(scale=0.5, axes='z')

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
2.337 ms	5.592 ms	0.443 ms	0.4×	5.3×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Transpose

Permute spatial axes for all aligned targets. This example exchanges Y and X.

A.Transpose(order='zxy')

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
2.025 ms	2.199 ms	0.611 ms	0.9×	3.3×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

UnsharpMask

Add a scaled difference between the image and a Gaussian-smoothed copy to enhance fine structure.

A.UnsharpMask(strength=0.5, sigma=1)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
10.299 ms	1.632 ms	0.394 ms	6.3×	26.2×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

Vignetting

Attenuate intensities with a radial falloff in all three dimensions, darkening voxels farther from the center.

A.Vignetting(strength=0.5, radius=0.8)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
28.360 ms	1.411 ms	0.173 ms	20.1×	164.3×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

ZoomBlur

Average multiple 3D zoom states around a center point to create radial blur.

A.ZoomBlur(strength=0.15, samples=9)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
257.941 ms	3.091 ms	1.939 ms	83.5×	133.1×

Z slice · YX plane	Y slice · ZX plane	X slice · ZY plane

2D tracing companions

These operate on 128 × 128 inputs derived from the 128³ crop. Their timings are for 2D work and they have no Z/Y/X volume slice triplet.

StripGeometry

Warp a 2D strip with rotation and a smooth lateral offset. This companion example uses the central Z slice of the CT crop.

A.StripGeometry(degrees=5, smooth_offset=2, stride=16)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
1.092 ms	0.038 ms	0.042 ms	28.9×	25.8×

2D-only companion: central Z slice, not a cubic-volume timing.

SurfacePerturb

Add localized smooth displacements to a 2D coordinate surface. This companion uses a synthetic surface with CT-derived height; the last panel shows displacement magnitude.

A.SurfacePerturb(blobs=3, amplitude=1, points_affected=100, p=1)

CPU	CUDA uploaded	CUDA resident	Uploaded speedup	Resident speedup
0.527 ms	0.287 ms	0.269 ms	1.8×	2.0×

2D-only companion: synthetic coordinate surface with CT-derived height. Preview is displacement magnitude, not an anatomical surface.

Image analysis and training losses

transforms.StructureTensor(sigma=1, rho=1.5) returns Pavel structure tensor components in zz, zy, zx, yy, yx, xx order (2D: yy, yx, xx). Set rho=0 to omit component smoothing. vc_reader.torch.structure_tensor accepts (B,1,Y,X) or (B,1,Z,Y,X) tensors and supplies first-order autograd on CPU and CUDA using native forward and backward kernels; MPS tensors stage through the CPU kernels and retain autograd. Double backward is not supported. For the central-plane 2D Pavel convention, pass scale=3/8.

transforms.DistanceTransform(implicit_background=True) matches SciPy's all-one EDT convention: when no background feature exists, measure from (-1,0,...). The default instead assigns zero when the other class is absent.

For label ridge filtering, compose GaussianBlur(sigma=2, boundary="reflect", truncate=4, full_width=False), Normalize("minmax"), and Ridges(). Ridges uses second-order one-sided differences at boundaries and requires at least three samples on each spatial axis.

CoherenceEnhancingDiffusion processes a single channel-first (C,Y,X) image or (C,Z,Y,X) volume. Channels are independent, labels are unchanged, and the usual on, channels, p, and p_channel options apply. Both backends preserve float32/float64 input dtype. Volume buffers and eigendecomposition use that dtype; filter coefficients and CPU convolution accumulation follow StructureTensor.

from vc_reader import transforms as A

diffuse = A.Pipeline([
    A.CoherenceEnhancingDiffusion(sigma=1, rho=2, lambda_=0.01,
                                 exponent=2, alpha=0.01,
                                 step_size=0.1, num_steps=10),
], backend="cuda")  # or "cpu"
result = diffuse({"image": volume[None]}, seed=0)["image"][0]

volume is a contiguous floating-point (Z,Y,X) array. CUDA device inputs and Pipeline.batch also work; iterations keep their buffers on the device. The corresponding C++ operation is vcr::transforms::CoherenceEnhancingDiffusion in a Pipeline, with the threshold member named lambda.

Each iteration recomputes a Pavel structure tensor with nearest-edge filtering. sigma controls presmoothing and rho controls tensor integration, both in voxels. The diffusion tensor follows the CED eigenvalue mapping in Mirebeau et al., equation (11): mu_i = alpha + (1-alpha)*exp(-(lambda_/(l_max-l_i))**exponent), with mu_i=alpha for zero eigengap. lambda_ is in squared input-intensity units; inputs are not automatically normalized. Scaling intensities by k requires scaling lambda_ by k*k to retain the same diffusion strengths.

The explicit stencil uses forward differences and their negative adjoint, with zero outward flux. It conserves total intensity (up to roundoff) and is L2 stable for the enforced 0 < step_size <= 1/6; alpha must be in (0,1]. num_steps=0 is identity. This stencil does not guarantee a maximum principle, so small intensity overshoots are possible; results are not clipped. Voxel spacing is assumed isotropic and equal to one. The complete volume and scratch buffers must fit memory; this API does not tile or stream diffusion.

Validate CED with the existing build and Python environment:

AGENTS_AGENT_MODE=1 CCACHE_DISABLE=1 cmake --build vc_reader/build-native-cuda --target _transforms test_native_transforms -j 2
AGENTS_AGENT_MODE=1 PYTEST_DISABLE_PLUGIN_AUTOLOAD=1 VC_READER_BUILD_DIR="$PWD/vc_reader/build-native-cuda" vesuvius/.venv/bin/python -m pytest vc_reader/tests/test_coherence_enhancing_diffusion.py vc_reader/tests/test_image_processing.py
ctest --test-dir vc_reader/build-native-cuda -R '^test_native_transforms$' --output-on-failure

geometry.rasterize_mesh_slice intersects XYZ triangles with a Z plane, rounds DDA samples to even, and returns labels plus optional normals and frames. Later triangles overwrite labels; normal/frame sums retain triangle order on both CPU and CUDA. Affine.load_transform_json reads registration matrices; geometry.transform_points(..., vectors=True, normalize=True) maps unit vectors using a supplied linear transform.

To validate these APIs against an existing local build (no installation):

AGENTS_AGENT_MODE=1 CCACHE_DISABLE=1 cmake --build vc_reader/build-native-cuda --target _transforms _reader _catalog -j 2
AGENTS_AGENT_MODE=1 PYTEST_DISABLE_PLUGIN_AUTOLOAD=1 VC_READER_BUILD_DIR="$PWD/vc_reader/build-native-cuda" vesuvius/.venv/bin/python -m pytest vc_reader/tests/test_image_processing.py

Run these commands from the monorepo root, substituting your existing CMake build and Python environment. Tests cover forward values, gradients, EDT edge cases, mesh labels/frames, and the bundled PHerc0139 CT crop when available.

Native C++ API

Readers, sampling, geometry, and the full transform catalog are available to C++20 clients through vc_reader::reader and vc_reader::transforms, including CUDA backends. Python and C++ use the same implementations. Configure with VC_READER_BUILD_PYTHON=OFF to build without Python or nanobind, and select VC_READER_ENABLE_CUDA=ON or OFF for the desired backend.

See the C++ guide for installation, target and header coverage, array/cache configuration, transform pipelines and replay, device-buffer ownership, streams, and a complete CPU/CUDA CT-normal example.

C API

libvc_reader_c (CMake target vc_reader::c, header <vc_reader/vc_reader_c.h>) exposes the render engine to C and to any language with a C FFI (Rust, C#, Swift, Go, Node, ctypes). It opens a volume's pyramid with its caches, chooses levels for a view, samples planes and vertex-grid or tifxyz surfaces with layers and composites, draws frames progressively (vcr_sample_cached now, vcr_pending_load on a worker, then draw again), and maps samples to display pixels. Handles are opaque, outputs go to caller buffers, every failure is a vcr_status with vcr_last_error(), and only vcr_* symbols are exported. A frontend's loop:

vcr_request_init(&req); req.volume = volume;
req.start_level = vcr_level_for_view(volume, zoom, 1.0f);
req.fallback_levels = vcr_fallback_levels(volume, req.start_level, w, h, zoom, 5);
vcr_sample_cached(&frame, &req, pending, &result);   /* draw result now */
if (!vcr_pending_empty(pending))                      /* on a worker thread: */
    vcr_pending_load(pending, redraw, ui);            /* then sample again */

tests/c/test_c_api.c builds as strict C99 and exercises each path, on the GPU too when one is available.

Native catalog

vc_reader::catalog provides the same Open Data parser and volume-registration resolver used by Python. Include <vc_reader/catalog.hpp> to load a catalog, browse its samples, scans, volumes, segments and artifacts, or build a RegistrationGraph. loadCatalog accepts a cache file and refresh policy; its result reports cached fallback and a warning when refresh fails.

The catalog target requires nlohmann_json, libcurl and zlib; it does not require Python, Qt, OpenCV, CUDA or the volume reader. Configure with VC_READER_BUILD_READER=OFF and VC_READER_ENABLE_CUDA=OFF for catalog-only use. Python requires the compiled _catalog extension but no NumPy or array backend for metadata browsing. Native regression tests are enabled by VC_READER_BUILD_NATIVE_TESTS=ON.

Skeletons, traces, and connected components

from vc_reader import thinning, transforms as A
from vc_reader.geometry import label_components

# uint8 contiguous (Y,X), or any leading dimensions batching independent YX masks.
skeleton = thinning.skeletonize(mask, method="guo_hall", backend="cuda")
trace_image, xy_paths = thinning.trace_ridges(mask, backend="cuda", return_traces=True)
labels, count = label_components(mask, connectivity=8)  # CUDA masks give CUDA labels

# Keep intermediate arrays on CUDA. Device inputs also produce device outputs.
pipeline = A.Pipeline([
    A.Morphology(on="mask", operation="close", radius=1),
    A.ConnectedComponents(on="mask"),
    A.RidgeTrace(on="mask"),
    A.MedialSurface(on="mask", method="guo_hall", tube=False, close=False),
], backend="cuda")
result = pipeline({"mask": mask[None]}, seed=0)
# mask_components, mask_traces, and mask_skel are added to the result.

Skeleton methods are zhang_suen, guo_hall, and zhang_suen_guo_hall. All nonzero input pixels are foreground; standalone skeleton outputs are 0/255. The standalone default preserves outer pixels; preserve_border=False treats outside pixels as background. MedialSurface defaults to the latter and retains input label values at skeleton pixels, with its existing optional morphology. Skeletons and ridge traces operate in YX planes; a stack is not a 3D skeleton. Components and Manhattan distances operate across the complete 2D or 3D grid. Component labels are int32 1..count in scan order of each component's first voxel, as scipy.ndimage.label numbers them, with zero for background. ConnectedComponents labels each channel independently (connectivity, 0 for full).

RidgeTrace/trace_ridges return visit-index pixels (seeds are 1, indices saturate at 255). return_traces=True also returns ordered XY polylines. min_length prunes polylines by pixel length and rasterizes survivors as 255. CPU TraceWorkspace buffers can be reused, with one workspace per thread. CUDA performs distance, seeds, components, and tracing on the device. Independent components run in parallel; within each component seed order is preserved. Short speculative path prefixes reduce serial work and are accepted only where the ordered visited mask allows them. Large single components can limit acceleration. Returning polylines copies the final paths to the host; optional pruning then runs on the host.

For resident arrays, skeletonize_into(input, output, backend="cuda") and trace_ridges_into(input, output) use preallocated CUDA uint8 buffers. Calls complete on vc_reader's stream; finish external writes before calling. The transforms pipeline keeps intermediate images on device, including in Pipeline.batch; convergence/count checks can synchronize small scalars. CUDA initialization is a one-time cost, so repeated processing should keep a process/pipeline alive. See the thinning benchmark for reproducible latency and reference comparisons.