from typing import Literal, NamedTuple, Protocol
from pathlib import Path
import logging
import requests
import scipy.io
import h5py
import numpy as np
logger = logging.getLogger(__name__)
[docs]
class Points(NamedTuple):
ED: np.ndarray
ES: np.ndarray
ALL_NOT_SUPPORTED = (
"The PCA atlas of all 4,329 UK Biobank subjects (--all) is not supported, and there is "
"no plan to support it in the near future. It uses a 5352-node version of the template "
"with the aortic and pulmonary outflow tracts removed, which does not match the 5810-node "
"biv-me template used to build the surfaces."
)
# The atlas files are mirrored as assets on a GitHub release since the
# original download links at cardiacatlas.org are no longer reliable.
ATLAS_BASE_URL = "https://github.com/ComputationalPhysiology/ukb-atlas/releases/download/atlas-data"
ATLAS_FILES = {
False: (
"UKBRVLV.h5",
"f6145755b5bf0fd43c2a26a5d20a14a30c32bb657a3ee61120ac1eecd06bef08",
),
True: (
"UKBRVLV_ALL.h5",
"a116bbaeae7d31b42df5c7562bffd4bc81717d92f7868fd8c255828f4e3f98aa",
),
}
[docs]
def download_atlas(outdir: Path, all: bool = False) -> Path:
"""Download the UK Biobank atlas from the Cardiac Atlas Project.
Parameters
----------
outdir : Path
Directory to download the atlas to.
all : bool, optional
If true, download the PCA atlas derived from all 4,329 subjects
from the UK Biobank Study. If false, downlaod PCA atlas derived
from 630 healthy reference subjects from the UK Biobank Study
(see [1]_) by default False
Returns
-------
Path
Path to the downloaded file.
References
----------
.. [1] Petersen, Steffen E., et al. "Reference ranges for cardiac
structure and function using cardiovascular magnetic resonance
(CMR) in Caucasians from the UK Biobank population cohort.
" Journal of cardiovascular magnetic resonance 19.1 (2016): 18.
"""
import hashlib
outdir = Path(outdir)
outdir.mkdir(exist_ok=True)
name, sha256 = ATLAS_FILES[all]
url = f"{ATLAS_BASE_URL}/{name}"
path = outdir / name
if not path.exists():
logger.info(f"Downloading {url} to {path}. This may take a while.")
# Download to a temporary file so an interrupted download is not cached
tmp_path = path.with_suffix(".h5.part")
digest = hashlib.sha256()
with requests.get(url, stream=True, timeout=60) as response:
response.raise_for_status()
with open(tmp_path, "wb") as f:
for chunk in response.iter_content(chunk_size=1 << 20):
f.write(chunk)
digest.update(chunk)
if digest.hexdigest() != sha256:
tmp_path.unlink()
raise RuntimeError(
f"Checksum mismatch for {url}: expected {sha256}, got {digest.hexdigest()}"
)
tmp_path.replace(path)
logger.info("Done downloading.")
else:
logger.debug(f"File {path} already exists. Skipping download.")
return path
[docs]
def generate_points(
filename: Path,
mode: int = -1,
std: float = 1.5,
score: np.ndarray | None = None,
) -> Points:
"""Generate points from the UK Biobank atlas.
Parameters
----------
filename : Path
Path to the UK Biobank atlas file.
mode : int, optional
Mode to generate points from. If -1, generate points from the mean
shape. If between 0 and the number of modes, generate points from
the specified mode. By default -1
std : float, optional
Standard deviation to scale the mode by, by default 1.5
score : np.ndarray | None, optional
PCA scores to generate points from. If None, use the mode and std
parameters to generate points. By default None
Returns
-------
Points
Named tuple containing the end-diastolic (ED) and end-systolic (ES)
points.
"""
logger.info(f"Generating points from {filename}")
with h5py.File(filename, "r") as hdf:
S = compute_S(hdf, mode, std, score=score)
# First half is ED and second half is ES
N = S.shape[1] // 2
ed = np.reshape(S[0, :N], (-1, 3))
es = np.reshape(S[0, N:], (-1, 3))
return Points(
ED=ed,
ES=es,
)
[docs]
def generate_points_burns(
filename: Path,
mode: int = -1,
std: float = 1.5,
score: np.ndarray | None = None,
) -> Points:
"""Generate points from the Burns atlas.
Parameters
----------
filename : Path
Path to the Burns atlas file.
mode : int, optional
Mode to generate points from. If -1, generate points from the mean
shape. If between 0 and the number of modes, generate points from
the specified mode. By default -1
std : float, optional
Standard deviation to scale the mode by, by default 1.5
score : np.ndarray | None, optional
PCA scores to generate points from. If None, use the mode and std
parameters to generate points. By default None
Returns
-------
Points
Named tuple containing the end-diastolic (ED) and end-systolic (ES)
points.
"""
logger.info(f"Generating points from {filename} (Burns atlas)")
data = scipy.io.loadmat(filename)
hdf = data["pca200"][0, 0]
S = compute_S(hdf, mode, std, score)
N = S.shape[0] // 2
ed = np.reshape(S[:N, 0], (-1, 3))
es = np.reshape(S[N:, 0], (-1, 3))
return Points(
ED=ed,
ES=es,
)
[docs]
class AtlasFile(Protocol):
def __getitem__(self, item: Literal["MU", "COEFF", "LATENT"]) -> np.ndarray: ...
[docs]
def compute_S(
hdf: AtlasFile | h5py.File,
mode: int = -1,
std: float = 1.5,
score: np.ndarray | None = None,
) -> np.ndarray:
"""Compute the shape matrix S from the PCA atlas.
Parameters
----------
hdf : h5py.File
HDF5 file containing the PCA atlas.
mode : int, optional
Mode to generate points from. If -1, generate points from the mean
shape. If between 0 and the number of modes, generate points from
the specified mode. By default -1
std : float, optional
Standard deviation to scale the mode by, by default 1.5
score : np.ndarray | None, optional
PCA scores to generate points from. If None, use the mode and std
parameters to generate points. By default None
Returns
-------
np.ndarray
Shape matrix S.
"""
if score is None:
logger.info(f"Using mode {mode} and std {std}")
mu = np.transpose(hdf["MU"])
if mode == -1:
S = mu
else:
if mode < 0 or mode >= hdf["COEFF"].shape[0]:
raise ValueError(
f"Mode {mode} is out of bounds. Needs to be between "
f"0 and {hdf['COEFF'].shape[0] - 1}"
)
eigenvalue = hdf["LATENT"][0, mode]
eigenvector = hdf["COEFF"][mode, :]
S = np.transpose(hdf["MU"]) + (std * np.sqrt(eigenvalue) * eigenvector)
else:
num_scores = len(score)
latent = np.asarray(hdf["LATENT"]).flatten()[
:num_scores
] # works for (1, modes) and (modes, 1)
d = score * np.sqrt(latent) # (num_scores,)
coeff = hdf["COEFF"]
if coeff.shape[0] == np.asarray(hdf["LATENT"]).size: # ukb h5 format: COEFF is (modes, N)
eigvecs = coeff[:num_scores, :] # (num_scores, N)
else: # Burns format: COEFF is (N, modes)
eigvecs = coeff[:, :num_scores].T # (num_scores, N)
mu = np.transpose(hdf["MU"])
S = mu + np.matmul(d, eigvecs).reshape(mu.shape)
return S