Vasculature.XCAT turns XCAT's surfaces, meshes and label volumes into the domains, inlet trees and blockers grow takes.
src/xcat/classify.jl
tissue_numbers(mask) -> Array{Int32, 3}
Each voxel's number among the true voxels of mask (in findall order), 0 elsewhere: the tissue voxel numbers flow and phantom use.
nearest_tissue_voxel(numbers, k) -> Int32
The number of the tissue voxel at grid index k (a tuple), or of the nearest one within 2 voxels (0: none).
imaged_segments(tree) -> (indices, starts, ends, radii)
The imaged (base) segments of a grown tree (vtree): their indices in the tree, end points (3 × n, mm) and radii (mm).
nearest_segment(segments, p) -> (index, distance)
Of segments (imaged_segments), the one whose surface is nearest to p (mm): its index in the tree and the distance from p to its surface (mm; negative inside the lumen).
nearest_segment(index::SegmentIndex, p) -> (index, distance)
As nearest_segment over the indexed segments, for a point within the index's reach of a lumen; elsewhere (0, Inf32).
SegmentIndex
A bucket grid over a tree's imaged_segments for fast nearest-lumen queries near the vessels: each cell lists the segments whose lumen, widened by the index's reach, touches it. Build it with segment_index.
segment_index(segments; reach, cell = 4.0) -> SegmentIndex
An index over segments (imaged_segments) for nearest_segment queries at points within reach mm of a lumen; cell is the bucket size (mm).
closing(mask, r) -> Array{Bool, 3}
Morphological closing of a 3-D mask by a ball of r voxels: dilate, then erode, each from an exact squared Euclidean distance transform (Felzenszwalb & Huttenlocher, one line at a time).
great_vessel_transit(labels, origin, h; cell_mm = 1.0, svc_mm = 70) -> NamedTuple
Which great vessel each vessel voxel belongs to, and where it lies along it, on a cell_mm grid: its tree (0 off the vessels) and the blood (mL) between it and the heart. labels is XCAT's label volume with voxel (1, 1, 1) centred at origin (mm), spacing h. The trees grow breadth first through XCAT's vessel labels from where they meet the heart, by connection rather than by mesh:
the aorta (label 28) and the systemic arteries joined to it, from the LV blood;
the pulmonary arteries (label 25 and those in the lungs, 125), from the RV blood, never stepping beside the aorta, so they cannot run into it where the two touch; grown first, so the aortic tree takes the rest;
the superior vena cava, from the upper half of the RA's veins: its first svc_mm (the SVC itself, about 70 mm long [mahlon2007_svc_ct_length]) is tree 3, and the veins beyond it are the right (6) or left (7) brachiocephalic tree, by side (the liver's is the right);
the inferior vena cava, from the lower half of the RA's veins;
the pulmonary veins (label 24 and those in the lungs, 124), from the LA blood.
The venous trees grow at once, each voxel joining the first to reach it. A voxel's volume is its tree's volume nearer the heart: a cross-section of the pulmonary or venous trees cuts every parallel branch, which together carry the whole flow, so all of them count; in the aortic tree only the aorta counts, since its branches run beside it, not in its path. Past the arch's top the descending aorta carries its own, smaller flow, so it is timed from the arch's volume (apex_ml) on; a brachiocephalic voxel is timed from the SVC's top (confluence_ml). Returns (; origin_mm, h, tree, volume, descending, apex\_ml, confluence\_ml, tree\_ml, ostia\_ml); ostia\_ml is the aorta's volume where the coronary arteries (label 26) leave it.
scale_tree(transit, tree, k) -> NamedTuple
transit (great_vessel_transit) with the blood along tree tree scaled by k: each of its voxels keeps its place along the vessel, and the blood upstream of it is k times XCAT's. For a phantom whose drawn vessel is narrower than the blood the body model that times it holds (XCAT vmale_50's aorta). The marks on that tree (apex_ml, ostia_ml on the aorta, confluence_ml on the superior vena cava) scale with it; tree_ml stays XCAT's.
vessel_volume_at(transit, p, t) -> Float32
Where the great-vessel voxel at p (mm) lies along tree t of transit (great_vessel_transit): the blood (mL) between it and the heart, interpolated trilinearly over the tree's cells around it (smooth along the vessel), else its nearest cell on the tree's value.
vessel_tree_at(transit, p, candidates) -> UInt8
The great-vessel tree the voxel at p (mm) is on: of candidates, the one most of the cells around it are on (trilinear weights), else its nearest cell's, else the first candidate.
great_vessel_pool(transit, p, t) -> (VOX_POOL, pool, ml)
The body pool (XCAT_BODY_POOLS) of a voxel at p (mm) on great-vessel tree t, and the blood (mL) along its vessel that times it: the aorta (5) or, past the arch's top, the descending aorta (27, timed from the arch's end at its own flow); the pulmonary arteries (6); the superior (8) and inferior (11) venae cavae; the pulmonary veins (7); the right (9) and left (10) brachiocephalic veins, timed from the SVC's top.
src/xcat/heart.jl
Centerline
Polyline X (3×K, mm) with radius r (K, mm), sampled from a tube surface, running from the vessel's proximal end to its distal end.
centerline(s::NurbsSurface; step = 1.0) -> Centerline
Centerline of an XCAT tube surface: ring centroids and mean ring radii along v (u runs around the tube), resampled every step mm. The collapsed cap rings at both ends are dropped.
BaseTree
Binary base tree for one coronary root: node 1 is the ostium; parent[k] (0 for the ostium), positions X (3×n, mm), radii r and the vessel each node came from.
base_forest(lines; roots, attach_tol = 3.0) -> Vector{BaseTree}
Joins centerlines into one binary tree per root. roots names the vessels that start at an ostium (e.g. ["dias_lad1", "dias_rca1"]). Every other vessel attaches at the node nearest to its proximal end on the closest already-placed vessel (within attach_tol mm of that vessel's surface); its proximal end is the end nearer to that tree. A node may carry at most two children, so a third branch moves to the next node downstream.
write_base_gxl(path, t::BaseTree)
Writes the base tree in the GXL layout import_gxl reads: node 0 is the ostium ("root node"), leaves are "terminal node", and each edge carries the child node's radius.
heart_masks(S, g; phase = "dias", res = (160, 160)) -> NamedTuple of Array{Bool,3}
LV and RV myocardium and cavities of an XCAT heart on grid g: LV myocardium = inside lv_4 minus inside lv_0, plus the papillary muscles; RV myocardium = inside rvm minus rv minus the LV; LV cavity = inside lv_0 minus papillary muscles; RV cavity = inside rv minus the LV.
label_masks(L, origin, h, g; roi = trues(g.dims), lv = 15, rv = 16, lvpool = 19, rvpool = 20,
pools = (19, 20, 21, 22, 24, 25)) -> NamedTuple of Array{Bool,3}
Heart masks on the coarse grid g from a fine XCAT label volume L (voxel size h mm, voxel (1,1,1) centred at origin). Each coarse voxel looks at the fine voxels whose centres fall inside it:
lvmyo / rvmyo: at least half of them are LV / RV myocardium (the larger count wins);
lvcav / rvcav: at least half are LV / RV blood pool;
pool: any of them is blood (the four chamber pools, great arteries and veins by default);
blood: the fraction of them that is blood;
domain: myocardium with no blood at all, so a vessel grown inside it cannot enter a chamber.
XCAT reuses some label values for distant organs (e.g. 16 for the kidneys), so roi limits the muscle and cavity masks to the heart region; pool covers the whole grid. The default label values are those of the XCAT heart phantom.
vessel_tube(t::BaseTree, g; margin = 1.5) -> Array{Bool,3}
Voxels of g within r + margin mm of the base tree's centerline. Epicardial arteries run in the fat outside the myocardium; adding this tube to a tree's growth domain lets new segments connect to them without leaving the domain.
vessel_bridges(t::BaseTree, myo, g; every = 2.0, reach = 15.0, margin = 1.5) -> Array{Bool,3}
Epicardial arteries run in the fat a few mm outside the myocardium. Every every mm along the centerline, a capsule (radius r + margin) joins the vessel to the nearest myocardial voxel within reach mm, so new branches can pass from the vessel into the wall, as real penetrating branches do, without the growth domain having to include the whole fat layer.
src/xcat/labels.jl
label_mask(L, origin, h, g::VoxelGrid, labels; fraction = 0.5) -> Array{Bool, 3}
Cells of the (coarser) grid g in which at least fraction of the label volume's voxels carry one of labels. L is a label volume with voxel (1, 1, 1) centred at origin (mm) and spacing h (mm); a cell takes the voxels whose centres lie inside it.
label_blockers(L, origin, h, g::VoxelGrid, labels; factor = 2) -> (mask, origin_mm, cell_mm)
The blocker mask over the box of grid g, in the form grow and penetrating_stubs take: cells of factor × factor × factor label voxels (0.4 mm for XCAT's 0.2 mm), a cell blocked (1) when any of its voxels carries one of labels (chambers, great vessels, airways). L, origin and h are as for label_mask.
src/xcat/nurbs.jl
NurbsSurface
Cubic B-spline surface with control points P[:, j, i]: j along u (fastest in the file, nu points) and i along v (nv points); knot vectors u (length nu + 4) and v (length nv + 4). Coordinates in mm.
read_nrb(path) -> Dict{String, NurbsSurface}
Reads every surface in an XCAT .nrb file (for example vmale_50_heart.nrb), keyed by name ("dias_lad1", "sys_lv_0", ...).
write_nrb(path, surfaces)
Writes surfaces in XCAT .nrb format, in the given order (a vector of NurbsSurface, or a Dict, written sorted by name).
evaluate(s, t_u, t_v) -> NTuple{3, Float64}
Point on the surface at parameters (t_u, t_v), each within its knot range.
tessellate(s; nu = 64, nv = 64) -> (points::Matrix{Float64} (3×K), triangles::Matrix{Int} (3×T))
Samples the surface on a regular parameter grid and triangulates it. Degenerate triangles (collapsed caps) are dropped.
src/xcat/stubs.jl
penetrating_stubs(t::BaseTree, g::VoxelGrid, tissues; scale = 1, spacing = 2.8scale, reach = 6scale, depth = 1.5scale,
domain = nothing, forbid = nothing, blockers = nothing, radius = 0.1) -> (tree, report)
Adds stubs, short branches from the vessels of t into the tissue they run beside, about every spacing mm along each vessel. A stub runs from a node of t to the nearest voxel of tissues and depth mm on into it; its tip is a new end of the tree (radius radius mm until growth sizes it by flow).
The step decides where stubs are needed, so it can be run on any inlet tree:
Sites lie about spacing mm apart along each vessel: the node nearest one spacing past the previous site, a side branch of t (an imaged branch counts as a stub that is already there) or the root. The first node of every branch is a site. Sample the centrelines at spacing or finer (centerline(s; step = spacing)).
A site gets no stub when it already has two children (the tree is binary); when no voxel of tissues lies within reach mm; when the tip lies outside the grid, outside domain or in forbid; or when the stub crosses blockers.
So a densely imaged tree gets few stubs, a bare trunk one per spacing, and stretches far from the tissue (e.g. a left main running free in fat) none. To grow from the imaged tree alone, do not call it.
Where the vessel runs through the tissue (its nearest tissue voxel lies within its own radius, as for a femoral artery inside a soft-tissue domain), the stub leaves sideways, perpendicular to the vessel and radius_k + depth mm long, turning by the golden angle from one such stub to the next; of eight directions the first whose tip is inside and clear is taken.
tissues is one Bool mask on g or several (the heart's LV and RV walls); where several are within reach, successive stubs alternate between them. domain (default: the union of tissues) and forbid are Bool masks on g; blockers is (mask, origin_mm, cell_mm) as for grow, a fine mask no stub may cross (chambers, great vessels). report counts the sites and what became of them: (; sites, stubs, sideways, no_tissue, outside, blocked).
Lengths: the heart's are the defaults at scale = 1: spacing 2.8 mm (the pig LAD's mean segment length; RCA 3.2, Kassab et al. 1993), reach 6 mm, depth 1.5 mm. Like every other length set for the heart, they scale with the organ (docs/deviations.md §16): pass scale = size_scale(spec) and sample the centrelines at 2.8 * scale mm. An organ whose feeding arteries' branch spacing is measured passes that as spacing.
src/xcat/supply_gaps.jl
fill_supply_gaps(parent, X, r, domain, lower; r_resolved, h = 1.0, max_gap_mm = nothing, min_gap_ml = nothing, qexp = 7 / 3,
step_mm = 2.0, forbidden = nothing, region = nothing, free_mm = 3.0, max_iter = 200)
-> (; parent, X, r, added, gaps)
parent (0 = root), X (3 × n node positions, mm) and r (radii, mm) describe the tree as nodes; segment k runs from parent[k] to k. domain is the organ mask on a grid whose voxel (1, 1, 1) is centred at lower (mm) with spacing h. A gap is a connected set (6-neighbour) of domain voxels farther than max_gap_mm from every segment, of at least min_gap_ml. Both follow from the imaging's resolution by default: the most tissue a vessel thinner than r_resolved (the smallest radius the imaging resolves, mm; required, since a grown or flow-resized tree's smallest radius says nothing about the imaging) can supply is V_max = V_organ · (r_resolved / r_root)^qexp, so a region larger than V_max that is farther than G = (3 V_max / 4π)^(1/3) from the imaged tree must be fed by a vessel the imaging would have shown, and one is missing (min_gap_ml = V_max, max_gap_mm = G). Its branch ends at the gap voxel nearest the gap's centroid and starts at the nearest node whose radius is at least the new branch's radius and whose straight path to it obeys the organ's blockers: no point of it lies in forbidden (a Bool mask on the domain grid: chambers, airways, ...), and every point farther than free_mm from the start lies in the domain and, with a region label map (left / right lung, lobes, ...), in the target's own region. Nodes are placed every step_mm. A start node that already has two children is moved to the midpoint of its parent segment (no node gets three children). Radius: r_root × (V_gap / V_domain)^(1/qexp). Returns the new arrays, one entry per added branch in added, the gap volumes found in each round in gaps, and in unfilled the gaps no node can reach without crossing a blocker or a region boundary (an organ region with no imaged vessel of its own).
src/xcat/tubes.jl
read_mesh_blocks(path; names = nothing, offset = (0, 0, 0)) -> Dict{String, Vector{Tuple{Matrix{Float64}, Matrix{Int}}}}
The triangle blocks of an XCAT .raw mesh file, by name: each block as (points, triangles), points 3 × n (mm, plus offset) and triangles 3 × m (indices into the points), the form voxelize takes. A name heads every block of its surface, so a vessel system drawn as many tubes gives one entry per tube. names keeps only those surfaces. offset moves the mesh into another frame (XCAT writes its meshes offset mm from its NURBS frame; its log gives it).
tube_centerline(points; name = "tube", step = 1.0, max_nodes = 400) -> Centerline
The centreline of one tube mesh, given as its vertices (3 × n, mm): slices across its principal axis about step mm apart, each slice's centroid a node. XCAT's mesh tubes have one calibre, so every node takes the tube's median slice radius (the end slices are partial rings). A tube too short to slice is its two ends at its axis' extremes.
tube_forest(lines; toward, join_mm = 6.0, min_nodes = 12) -> Vector{BaseTree}
Joins tube centrelines into trees, largest tube first: a tube attaches by its nearer end to the nearest node of an already placed tube at least half its radius, when that end lies within the two radii plus join_mm of it (a node that already has a child is split at the midpoint of its parent segment, so the tree stays binary); otherwise it starts a new tree, whose root is its end nearer toward (mm; the hilum or the heart). Trees of fewer than min_nodes nodes (stray single tubes) are dropped. Each tree's nodes run from its root, parents first.
graft(trunk::Centerline, trees; name = trunk.name) -> BaseTree
One tree: the trunk (e.g. the main pulmonary artery from its valve) with trees joined to it, in order of their roots' distance from the trunk's distal end, each at the nearest node already in the tree (the trunk's or an earlier tree's). The trunk runs from the end farther from every tree root; its rings narrower than half its widest (a valve sinus or dome collapsing at an end) are dropped. A node that already has two children is split at the midpoint of its parent segment. The root edge takes the trunk's median ring radius, robust to the sinus and the bifurcation dome at its two ends.
drain(trees, surface; root_radius, name = "veins") -> BaseTree
One draining tree: the largest of trees takes the others, each joined by its root at its nearest node, and a straight segment carries the joined tree from its root to the nearest point of surface (3 × n points, mm: the chamber it drains into), which becomes the root, of radius root_radius (mm).