Metadata-Version: 2.4
Name: bone-contouring
Version: 0.5.0
Summary: SimpleITK-first bone contour and mask generation for volumetric bone images
Author: Matthias Walle
License-Expression: GPL-3.0-only
Requires-Python: >=3.11
Description-Content-Type: text/markdown
License-File: LICENSE
License-File: NOTICE
License-File: LICENSES/MIT-previous.txt
Requires-Dist: bone-imaging-derivatives>=0.1.7
Requires-Dist: numpy<3.0,>=1.26
Requires-Dist: SimpleITK>=2.3
Requires-Dist: scipy<2.0,>=1.11
Provides-Extra: unet
Requires-Dist: torch>=2.2; extra == "unet"
Requires-Dist: scikit-image<0.26,>=0.24; extra == "unet"
Requires-Dist: aimio-py>=0.1.8; extra == "unet"
Provides-Extra: geodesic
Requires-Dist: hrpqct-geodesic-contour>=0.1.2; extra == "geodesic"
Provides-Extra: test
Requires-Dist: pytest>=8; extra == "test"
Requires-Dist: pytest-cov>=7; extra == "test"
Dynamic: license-file

<p align="center">
  <img src="resources/bone-contouring.png" alt="bone-contouring icon" width="240">
</p>

# 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.

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

```python
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:

```bash
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:

```bash
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)](https://doi.org/10.1038/s41598-022-27350-0).
Original code: [Bonelab/HR-pQCT-Segmentation](https://github.com/Bonelab/HR-pQCT-Segmentation).
Published weights: [Zenodo](https://zenodo.org/records/14755838).

```python
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.

```python
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.

```python
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.

| Default | Radius | Tibia | Knee adaptation |
| --- | --- | --- | --- |
| Cortical seed threshold (mg HA/cm³) | 500–3000 | 500–3000 | 150–3000 |
| Seed Gaussian sigma/support (voxels) | 2 / 3 | 2 / 3 | 2 / 3 |
| Minimum XY cortical rim (voxels) | 6 | 6 | 6 |
| Final 3D closing radius (voxels) | 30 | 50 | 36 |
| Initial corner-component minimum (voxels) | 800 | 200000 | 200000 |

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.
