bone-contouring

SimpleITK-first bone contouring and mask generation for volumetric bone images.

Standard contouring: Matthias Walle. Optional U-Net method and trained model: Nathan J. Neeteson, Bryce A. Besler, Danielle E. Whittier and Steven K. Boyd.

Published HR-pQCT U-Net (optional)

Install pip install 'bone-contouring[unet]' to use the fixed-default published 61 µm radius/tibia model. The base package does not require PyTorch, scikit-image or aimio-py. No Bonelab or vtkbone is required.

bone-contouring unet /data/raw --output /data/results --device auto
ssh host 'bone-contouring unet /data/raw --output /data/results --device cuda'

from bone_contouring.unet import Segmenter
masks = Segmenter(device="auto").segment(density_zyx)  # calibrated mg HA/cm³

Devices: auto (CUDA → MPS → CPU), cpu, cuda, mps. Explicit unavailable devices fail clearly. The only inference control is device; published preprocessing, five-slice channels, model and morphology are unchanged. XCT-I and knee accuracy, and equivalence to scanner/IPL contours, are not established.

Post-processing runs on the CPU even when inference uses CUDA/MPS. Repeated erosions/dilations are batched into compiled SciPy iterations with the same 3D connectivity, border handling, padding and iteration counts; no smoothing or compartment defaults are reduced. Progress reports the four morphology stages and their total elapsed time. A read-only source-checkout benchmark compares masks voxel-for-voxel against the published implementation:

python benchmarks/benchmark_unet_postprocessing.py

First use downloads SHA256-verified radius_tibia_final.pth to the shared cache. Use HRPQCT_SEGMENTATION_MODEL_DIR to share a pre-provisioned cache, or HRPQCT_SEGMENTATION_WEIGHTS for a verified explicit weights file. Slicer stores weights in its application data HRpQCTSegmentation/models, beside MotionScore. Setup installs dependencies, not model weights.

The CLI writes native 0/127 full/trab/cort AIM compartment masks with original geometry and standard BoneContours names, individual provenance sidecars and a _UNET.json completion marker. It never overwrites existing artifacts and reuses one model per batch. Complete files are published atomically; interrupted cases without a valid completion marker are withheld by derivative discovery. These are not bone-tissue SEG masks. The normalized Slicer batch profile adds XCTII Laplace–Hamming tissue SEG and material labels in the same run, then loads full/trab/cort/SEG into one segmentation node. Completed U-Net compartments are reused without inference when only SEG/material is missing; imported/conflicting or incomplete compartments are preserved and blocked rather than overwritten. An existing tissue SEG is reused only with matching XCTII LH provenance/defaults; an incompatible or unverified SEG is preserved and reported as a conflict. The same one-case normalized recipe is available without Slicer:

bone-contouring unet /data/study/sub-001/ses-1/xct/sub-001_ses-1_voi-radiusleft_xct.AIM --output /data/study/derivatives/BoneContours/sub-001/ses-1/xct --dataset-root /data/study --device auto

LH uses the original AIM native intensities, not density-calibrated U-Net inputs. Site comes from the normalized radius/tibia VOI. Standalone raw-folder inference without --dataset-root remains compartment-only; existing standard batch profiles can also add missing SEG/material while preserving compartments.

Code is GPL-3.0-only from version 0.3.0; previously distributed MIT versions retain their original terms. See LICENSE and NOTICE for notices and upstream credit. The original inference draft also retains GPL-3.0-only; it need not be published as another PyPI package.

Method: Neeteson et al., Scientific Reports 13, 252 (2023). Original code: Bonelab/HR-pQCT-Segmentation. Published weights: Zenodo.

from bone_contouring import generate_masks_from_image, resolve_preset

masks = generate_masks_from_image(image, resolve_preset(modality="xct1", site="radius"))

Standard XCTI periosteal thresholds are 250 mg HA/cm³ for radius/tibia and 150 mg HA/cm³ for knee, in both resolved and shipped batch presets. These apply to calibrated density inputs, not the native-gray Laplace–Hamming tissue threshold. XCTII, endosteal thresholds, morphology, and smoothing are unchanged. Saved custom profiles retain their explicit thresholds; existing masks are not automatically regenerated when defaults change.

For manually configured images, use resolve_preset(modality="custom", site="none", segmentation="gauss") and set thresholds in the input image's units. Custom recipes do not apply XCTII radius/tibia threshold calibration. Users can save their own named profiles (for example, for micro-CT) for scene and batch use; there is no built-in micro-CT preset.

Dual threshold

The current compartment method is named Dual threshold (dual_threshold). The earlier buie, ipl_match, and stable_3d public methods are retired. Saved profiles selecting them raise a replacement error; they are never silently converted to another algorithm. The previous standard name remains an input alias for dual_threshold. Existing masks are not automatically regenerated.

Tissue segmentation

Gaussian tissue segmentation defaults to a single filter of the original density image with sigma 0.8 voxels and support 1 voxel (a finite 3×3×3 sampled Gaussian, reflected at image boundaries), followed by thresholds of 320 mg HA/cm³ in trabecular and 450 mg HA/cm³ in cortical compartments. Sigma is scaled by the smallest voxel spacing; it is not specified in millimetres. Contour prefilters and signed-distance smoothing are separate and do not pre-smooth the tissue input. Explicit custom-profile settings remain unchanged. Tissue SEG retains disconnected bone components; largest-component selection belongs to downstream FEA preparation and is never applied here. The legacy segmentation.keep_largest_component field is accepted but ignored, including in saved profiles. The separate 64-voxel minimum-component noise filter remains.

Experimental IPL script candidate

Select inner_contour="ipl" for the IPL candidate translating STEP_1 of the supplied IPL_UPAT_CALGARY_EVAL_XT2_NOREG.COM (Steven K. Boyd, Danielle E. Whittier). This fixed candidate is opt-in, for XCTII radius/tibia only. The shared standard uses a tunable adaptation described below. The script consumes a supplied outer GOBJ, so this candidate uses the full ROI from the separately selected outer-contour method.

params = resolve_preset(modality="xct2", site="tibia", segmentation="gauss",
                        outer_contour="standard", inner_contour="ipl")
masks = generate_masks_from_image(density_image, params)

The fixed recipe uses a sigma-2/support-3 cortical seed at 500–3000 mg HA/cm³, inversion/rank extraction, a six-voxel XY peel, 3D erosion–component–dilation at distance 3, close/open at 15, the literal corner-number filters, then close at 30 (radius) or 50 (tibia), axial 50–100% component cleanup and complementary CORT/TRAB masks. The final corner filter keeps 1–800 voxels even after the site-specific minimum of 800/200000; that apparent script inconsistency is preserved. No standard signed-distance smoothing or hole-fill fallback is added. The final TRAB = all - cleaned CORT can reassign discarded peripheral cortical fragments to TRAB outside the previous peel; metadata reports those voxels.

This is not validated native IPL equivalence. The implementation currently uses voxel-centre Euclidean balls, background outside XY, and edge-replicated terminal slices in Z during morphology, cropped back afterward. This avoids artificial endosteal scan-end caps; Z is not cropped at intermediate ROI bounds. The large 3D opening/closing distances are unchanged; no extra smoothing is added. Native stack-end handling, Scanco metric 11, GOBJ peel rasterization, Gaussian support/margins, bounding-box margins and component ties remain unverified. Metadata records these assumptions and intermediate counts/timings. Tissue SEG retains the separately configured package method; IPL STEP_2 compartment-boundary filtering is not emulated here. The STEP_1 recipe is commented out in the supplied script. Validate against native outputs generated with that step enabled and the same outer ROI before using this experimental candidate as a scanner replacement.

Dual threshold implementation

dual_threshold uses the existing repaired compartment sequence across XCTI and XCTII, for radius, tibia, and knee, including shipped profiles and Slicer. This naming cleanup does not change the numerical method. Input must represent one target bone. Native-mask comparisons currently cover XCTII radius/tibia, not XCTI or knee; native IPL equivalence is not established.

params = resolve_preset(
    modality="xct2", site="radius", segmentation="gauss",
    outer_contour="dual_threshold", inner_contour="dual_threshold",
)
masks = generate_masks_from_image(density_image, params)

The outer stage independently filters and thresholds calibrated density using the effective preset settings, selects the largest 6-connected bone component, dilates in XY, fills axial holes before erosion, then erodes with the same radius. Completing ordinary closing before filling can reopen a narrow bridge and leave marrow connected to the exterior. XCTII radius/tibia defaults remain sigma 0.8 voxel-relative and threshold 320 mg HA/cm³; other presets retain their existing density thresholds, Gaussian sigmas, and outer closing radii. It deliberately omits the legacy pre-fill opening: that operation can break a thin shell and prevent marrow filling. The ignored outer controls are periosteal_open_radius, fill_holes, and use_adaptive_threshold.

The inner stage uses the supplied STEP_1 order: a finite Gaussian cortical seed, inversion/rank selection inside the six-voxel XY-peeled outer ROI, 3D erosion and dilation by 3, closing/opening by 15, the documented corner filters, a final site-specific close, and axial 50–100% component cleanup. Morphology uses voxel-centre Euclidean balls with XY background and edge-replicated scan ends in Z; intermediate cropping never introduces a false Z boundary.

These voxel settings are shared by XCTI and XCTII; they are not physically resolution-independent. Knee retains its 36-voxel final close and uses the tibia corner rule, an unvalidated adaptation for one target bone. Standard knee outer threshold is 150 in both scanner presets. inner.endosteal_threshold, gaussian_sigma, peel and trabecular_close_radius remain effective custom controls. Legacy inner kernel/adaptive/Buie and signed-distance smoothing settings do not drive the new inner sequence. There is no additional signed-distance smoothing of TRAB. The outer contour retains its existing physical boundary regularization.

After the script-style cleanup, standard reapplies the configured XY peel; CORT is exactly FULL minus TRAB. This intentionally differs from the fixed candidate's literal final complement, which may reassign small peripheral fragments outside the rim. Set inner.peel=0 to disable the constraint. The rim is a compartment minimum, not measured cortical thickness or an original Buie requirement. Native metric 11, coefficients, margins, and border handling remain unverified. XCTI and knee transfer require further validation.

Metadata includes effective settings, timings, seed sizes, and advisory QA for empty/fragmented slices, axial holes, crop contact, area jumps, and adjacent boundary jumps. Warnings neither trigger a synthetic fallback nor establish anatomical correctness; true anatomy, motion, or a cropped input can trigger them. The absence of warnings is also not proof of correctness. Inspect the input and contours, especially when trabecular ROI reaches the outer boundary.

Use calibrated density and the image API to force regeneration. Existing batch mask reuse is unchanged; existing masks are not automatically regenerated. Settings hashes and provenance include standard_algorithm=dual_threshold_v1, the inner sequence revision and effective settings for the replacement standard. Tissue-segmentation methods/settings are unchanged; their contour-support override no longer drives this standard's envelopes. Select --method dual_threshold in benchmarks/benchmark_buie.py to produce separate results. Adjacent-boundary metrics measure envelope movement in mm (including anatomy/motion), with omitted empty pairs counted; they complement, rather than replace, native overlap and cortical-volume bias.