# mr-to-pct -- pseudo-CT from a T1-weighted MR image for transcranial
# ultrasound (TUS) acoustic simulation (Yaakub et al., Brain Stimulation
# 16(1), 2023): a 3D residual U-Net (MONAI) trained to synthesise a 100 keV
# CT. The output is on the input T1's grid, so BabelBrain takes it as an
# ordinary CT with "co-register CT to T1" off, and k-Wave/SimNIBS pipelines
# can use it the same way.
#
# Layout facts the spec depends on (checked against the pinned commit):
# - Upstream is a notebook, a script and a package named `utils` (netdef,
#   infer_funcs, infer_funcs_noants). A top-level `utils` cannot go into
#   site-packages, so the tree is installed privately in
#   %%{_datadir}/mr-to-pct and %%{_bindir}/mr-to-pct runs the script, which
#   imports utils relative to itself (Patch1).
# - The trained weights are not in the repository: they are on OSF
#   (osf.io/e7sz9, "mr-to-pct", CC BY 4.0) and ship here as Source1, next to
#   the script, where the script looks by default.
# - The T1 preparation step uses ANTsPy, which Fedora does not have; Patch1
#   adds the same steps in SimpleITK as the fallback.

%global commit      60473b73f57c97ffbf8baa175f7987e4a2ddfc15
%global shortcommit %(c=%{commit}; echo ${c:0:7})
%global snapdate    20250225
%global weights     pretrained_net_final_20220825.pth

Name:           mr-to-pct
Version:        1.0.0^%{snapdate}git%{shortcommit}
Release:        1%{?dist}
Summary:        Pseudo-CT from T1-weighted MRI for transcranial ultrasound planning

# Code: Apache-2.0 (LICENSE). Trained network weights (Source1): CC-BY-4.0,
# the licence of the OSF project that distributes them.
License:        Apache-2.0 AND CC-BY-4.0
URL:            https://github.com/sitiny/mr-to-pct
# main at 2025-02-25: v1.0.0 (2022-09) plus fixes to both inference modules
# and the script, then README changes.
Source0:        %{url}/archive/%{commit}/mr-to-pct-%{shortcommit}.tar.gz
# sha256 c2c1f38fbd7806796d738f4d520ca1f63ce172649613ff513cd9ae43788cfe18
Source1:        https://osf.io/download/c3w98/#/%{weights}
# MONAI removed AddChannel in 1.3; both inference modules imported it unused.
# Submitted upstream by the packager.
Patch0:         mr-to-pct-monai-1x-imports.patch
# Command-line options for the script (--weights, --prep/--no-prep,
# --device), imports relative to the script, and T1 preparation in SimpleITK
# when ANTsPy is missing. Submitted upstream by the packager.
Patch1:         mr-to-pct-script-cli.patch

BuildArch:      noarch
BuildRequires:  python3-devel
BuildRequires:  python3-torch
BuildRequires:  python3-monai
BuildRequires:  python3-numpy
BuildRequires:  python3-nibabel
BuildRequires:  python3-simpleitk
BuildRequires:  python3-matplotlib

Requires:       python3-torch
Requires:       python3-monai
Requires:       python3-numpy
Requires:       python3-nibabel
Requires:       python3-simpleitk
Requires:       python3-matplotlib
# The pseudo-CT is a CT input for these planners.
Suggests:       babelbrain

%description
mr-to-pct synthesises a pseudo-CT image from a T1-weighted MR image with a
trained 3D residual U-Net, for use as the skull model in acoustic
simulations of transcranial ultrasound stimulation (TUS) when no CT of the
participant is available. The pseudo-CT is written on the T1's grid, in
Hounsfield units, and can be given to BabelBrain as a CT input or used with
k-Wave and other solvers.

The trained network weights (CC BY 4.0) are included. Run
"mr-to-pct T1.nii pct.nii"; by default the T1 is bias-corrected and
head-masked first.

%prep
%autosetup -n mr-to-pct-%{commit} -p1
# Run by %%{_bindir}/mr-to-pct through python3; nothing is executed directly.
sed -i '1{/^#!/d}' mr-to-pct_script.py
echo "c2c1f38fbd7806796d738f4d520ca1f63ce172649613ff513cd9ae43788cfe18  %{SOURCE1}" | sha256sum -c -

%build
# Nothing to build.

%install
mkdir -p %{buildroot}%{_datadir}/%{name}
cp -pr mr-to-pct_script.py utils %{buildroot}%{_datadir}/%{name}/
install -p -m 0644 %{SOURCE1} %{buildroot}%{_datadir}/%{name}/%{weights}
%py_byte_compile %{python3} %{buildroot}%{_datadir}/%{name}

mkdir -p %{buildroot}%{_bindir}
cat > %{buildroot}%{_bindir}/%{name} <<EOF
#!/usr/bin/bash
exec %{python3} %{_datadir}/%{name}/mr-to-pct_script.py "\$@"
EOF
chmod 0755 %{buildroot}%{_bindir}/%{name}

%check
root=%{buildroot}%{_datadir}/%{name}
export HOME=$PWD/home MPLBACKEND=Agg; mkdir -p "$HOME"
%{python3} "$root/mr-to-pct_script.py" --help | grep -q -- '--no-prep'

# The network builds, the shipped weights load into it key for key, and one
# patch of the inference runs to a finite, CT-like output. A full volume is
# about ten minutes on a CPU; tests/mr-to-pct.bats runs that.
cd "$root"
%{python3} - <<'EOF'
import numpy as np, torch
from utils.netdef import ShuffleUNet
from utils.infer_funcs_noants import do_inference_3D
net = ShuffleUNet(dimensions=3, in_channels=1, out_channels=1,
                  channels=(32, 64, 128, 256, 512), strides=(2, 2, 2, 2))
res = net.load_state_dict(torch.load("pretrained_net_final_20220825.pth", map_location="cpu"))
assert not res.missing_keys and not res.unexpected_keys, res
net.eval()
# A 256x256x40 slab: a bright ellipsoid "head" in a dark background, scaled
# to the network's [-1, 1] input range. step_size 8 keeps it to a few patches.
z, y, x = np.mgrid[0:256, 0:256, 0:40].astype(np.float32)
head = ((z - 128) / 90) ** 2 + ((y - 128) / 110) ** 2 < 1
vol = np.where(head, 0.4, -1.0).astype(np.float32)
with torch.no_grad():
    out = do_inference_3D(net, vol, torch.device("cpu"), step_size=8)
assert out.shape == vol.shape and np.isfinite(out).all(), (out.shape, np.isfinite(out).mean())
print("inference ok: output range %.1f .. %.1f" % (out.min(), out.max()))
EOF

# The SimpleITK preparation (Patch1) masks a synthetic head.
%{python3} - <<'EOF'
import os, tempfile, numpy as np, SimpleITK as sitk
from utils.prep_sitk import prep_t1_sitk
d = tempfile.mkdtemp()
z, y, x = np.mgrid[0:64, 0:64, 0:64]
head = ((z - 32) ** 2 + (y - 32) ** 2 + (x - 32) ** 2) < 22 ** 2
rng = np.random.default_rng(0)
arr = (np.where(head, 600.0, 20.0) * (1 + 0.3 * x / 64) + rng.normal(0, 5, head.shape)).astype(np.float32)
f = os.path.join(d, "t1.nii"); sitk.WriteImage(sitk.GetImageFromArray(arr), f)
out = sitk.GetArrayFromImage(sitk.ReadImage(prep_t1_sitk(f)))
m = out > 0
covered, grown = (m & head).sum() / head.sum(), m.sum() / head.sum()
# The mask is dilated and smoothed on purpose, so it is a little larger.
assert covered > 0.99 and grown < 1.6, (covered, grown)
print("prep ok: head covered %.3f, mask %.2fx the head" % (covered, grown))
EOF

%files
%license LICENSE
%doc README.md CITATION.cff mr-to-pct_infer.ipynb
%{_bindir}/%{name}
%{_datadir}/%{name}/

%changelog
* Wed Oct 14 2026 Morgan Hough <morgan.hough@gmail.com> - 1.0.0^20250225git60473b7-1
- Initial package, with the trained weights from OSF. Two patches submitted
  upstream: drop MONAI's removed AddChannel import, and command-line options
  with a SimpleITK T1 preparation when ANTsPy is missing.