# =============================================================================
# This file is part of LION library
# License : BSD-3
#
# Author : Ander Biguri
# Modifications: -
# =============================================================================
from __future__ import annotations
# math/science imports
"""CT operator, noise, and attenuation-unit utilities."""
import numpy as np
import tomosipo as ts
import torch
import torchvision.transforms.functional as TF
# AItomotools imports
from LION.CTtools.ct_geometry import Geometry
from LION.operators import CTProjectionOp
[docs]
def from_HU_to_normal(img):
"""
Converts image in Hounsfield Units (air-> -1000, bone->500) into a [0-1] image.
Comercial scanners use a piecewise linear function. Check STIR for real values. (https://raw.githubusercontent.com/UCL/STIR/85cc1940c297b1749cf44a9fba937d7cefdccd47/src/utilities/share/ct_slopes.json)
"""
if isinstance(img, np.ndarray):
return np.minimum(np.maximum((img.astype(np.float32) + 1000) / 3000, 0), 1)
elif isinstance(img, torch.Tensor):
return torch.clip((img.float() + 1000) / 3000, min=0, max=1)
else:
raise NotImplementedError
[docs]
def from_normal_to_HU(img):
"""Convert a LION-normalised CT image in ``[0, 1]`` to HU.
This is the inverse of :func:`from_HU_to_normal` on its represented
Hounsfield-unit interval, mapping 0 to -1000 HU and 1 to 2000 HU.
Values are not clipped so reconstruction overshoots remain visible.
"""
if isinstance(img, np.ndarray):
return 3000 * img.astype(np.float32) - 1000
if isinstance(img, torch.Tensor):
return 3000 * img.float() - 1000
raise NotImplementedError
[docs]
def from_HU_to_mu(img):
"""
Converts image in Hounsfield Units (air-> -1000, bone->500) into linear attenuation coefficient (air-> 0.0012,
bone->1.52 g/cm^3). Approximate.
Comercial scanners use a piecewise linear function. Check STIR for real values. (https://raw.githubusercontent.com/UCL/STIR/85cc1940c297b1749cf44a9fba937d7cefdccd47/src/utilities/share/ct_slopes.json)
"""
if isinstance(img, np.ndarray):
return np.maximum(
((1.52 - 0.0012) / (500 + 1000)) * (img.astype(np.float32) + 1000) + 0.0012,
0,
)
elif isinstance(img, torch.Tensor):
return torch.clip(
((1.52 - 0.0012) / (500 + 1000)) * (img.float() + 1000) + 0.0012, min=0.0
)
else:
raise NotImplementedError
[docs]
def sinogram_add_noise(
proj,
I0=1000,
sigma=5,
sigma_blur=0.3015,
ks_value=3,
flat_field=None,
dark_field=None,
enable_gradients=False,
):
"""
Wraper for _sinogram_add_noise to support gradients
"""
if enable_gradients:
sino = _sinogram_add_noise(
proj, I0, sigma, sigma_blur, ks_value, flat_field, dark_field
)
else:
with torch.no_grad():
sino = _sinogram_add_noise(
proj.detach(), I0, sigma, sigma_blur, ks_value, flat_field, dark_field
)
return sino
def _sinogram_add_noise(
proj,
I0=1000,
sigma=5,
sigma_blur=0.3015,
ks_value=3,
flat_field=None,
dark_field=None,
):
"""
Adds realistic noise to sinograms.
- Poisson noise, with I0 counts in a scanner with no sample (bigger value==less noise)
- Gaussian noise of zero mean and sigma std
- Detector crosstalk in % of the signal of adjacent pixels.
- Add a flat_field to add even more realistic noise (computed at non-corrected flat fields)
"""
dev = torch.cuda.current_device()
if torch.is_tensor(proj):
istorch = True
dev = proj.get_device()
if dev == -1:
dev = torch.device("cpu")
elif isinstance(proj, np.ndarray):
# all good
istorch = False
proj = torch.from_numpy(proj).cuda(dev)
else:
raise ValueError("numpy or torch tensor expected")
if dark_field is None:
dark_field = torch.zeros(proj.shape, device=dev)
if flat_field is None:
flat_field = torch.ones(proj.shape, device=dev)
max_val = torch.amax(
proj
) # alternatively the highest power of 2 close to this value, but lets leave it as is.
if max_val <= 0:
max_val = 1.0
Im = I0 * torch.exp(-proj / max_val)
# Uncorrect the flat fields
Im = Im * (flat_field - dark_field) + dark_field
# Add Poisson noise
Im = torch.poisson(Im)
# Detector cross talk
ks = int(sigma_blur * ks_value) * 2 + 1
if ks < 3:
ks = 3
Im = TF.gaussian_blur(Im, kernel_size=[ks, ks], sigma=[sigma_blur, sigma_blur])
# Electronic noise:
Im = Im + sigma * torch.randn(Im.shape, device=dev)
Im[Im <= 0] = 1e-6
# Correct flat fields
Im = (Im - dark_field) / (flat_field - dark_field)
proj = -torch.log(Im / I0) * max_val
proj[proj < 0] = 0
if istorch:
return proj
else:
return proj.cpu().detach().numpy()
[docs]
def from_HU_to_material_id(img):
"""
Converts an image in Hounsfield units into a material index
May require some image filtering preprocessing
"""
materials = img
materials[img < -950] = 0 # air
materials[(img > -950) & (img < -750)] = 1 # lung
materials[(img > -750) & (img < -150)] = 2 # bronqui
materials[(img > -150) & (img < -0)] = 3 # fat
materials[(img > 0) & (img < 150)] = 4 # muscle
materials[(img > 150) & (img < 300)] = 5 # bone marrow
materials[img > 300] = 6 # bone
materials = materials.astype(np.uint8)
return materials
[docs]
def make_operator(geometry: Geometry) -> CTProjectionOp:
"""Construct a differentiable LION CT projection operator.
Parameters
----------
geometry : Geometry
Fan- or parallel-beam acquisition geometry.
Returns
-------
CTProjectionOp
LION wrapper around the matching tomosipo operator.
"""
if not isinstance(geometry, Geometry):
raise ValueError(
"Input geometry is not of class LION.CTtools.ct_geometry.Geometry"
)
if geometry.mode == "fan":
vg = ts.volume(
shape=geometry.image_shape, size=geometry.image_size, pos=geometry.image_pos
)
pg = ts.cone(
angles=geometry.angles,
shape=geometry.detector_shape,
size=geometry.detector_size,
src_orig_dist=geometry.dso,
src_det_dist=geometry.dsd,
)
elif geometry.mode == "parallel":
vg = ts.volume(
shape=geometry.image_shape, size=geometry.image_size, pos=geometry.image_pos
)
pg = ts.parallel(
angles=geometry.angles,
shape=geometry.detector_shape,
size=geometry.detector_size,
)
else:
raise ValueError("Geometry mode not understood, has to be 'fan' or 'parallel'")
A = ts.operator(vg, pg)
A = CTProjectionOp(A)
return A
[docs]
def forward_projection(
image: np.ndarray | torch.Tensor,
geometry: Geometry,
backend: str = "tomosipo",
) -> torch.Tensor:
"""
Produces a noise free forward projection, given np.array image, a size (in real world units), a sinogram shape and size,
distances from source to detector DSD and distance from source to object DSO.
May support other backends than tomosipo
"""
if backend != "tomosipo":
raise ValueError("Only tomosipo backend for CT supported")
# You can add other backends here
import tomosipo as ts
device = torch.cuda.current_device()
if isinstance(image, np.ndarray):
image = torch.from_numpy(image).float().to(device)
if len(image.shape) == 3:
if image.shape[0] > 1: # there is no reason to have this constraint
raise ValueError("Image must be 2D")
elif len(image.shape) == 2:
image = torch.unsqueeze(image, 0)
else:
raise ValueError("Image must be 2D")
vg = ts.volume(shape=geometry.image_shape, size=geometry.image_size)
pg = ts.cone(
angles=geometry.angles,
shape=geometry.detector_shape,
size=geometry.detector_size,
src_orig_dist=geometry.dso,
src_det_dist=geometry.dsd,
)
A = ts.operator(vg, pg)
sino = A(image)[0]
return sino