From e167679de2617a0603478423a395a585f53ba70b Mon Sep 17 00:00:00 2001
From: Tamino Huxohl <thuxohl@techfak.uni-bielefeld.de>
Date: Thu, 17 Nov 2022 16:07:18 +0100
Subject: [PATCH] add script for reconstruction
---
mu_map/recon/recon.py | 228 ++++++++++++++++++++++++++++++++++++++++++
1 file changed, 228 insertions(+)
create mode 100644 mu_map/recon/recon.py
diff --git a/mu_map/recon/recon.py b/mu_map/recon/recon.py
new file mode 100644
index 0000000..1d48d6b
--- /dev/null
+++ b/mu_map/recon/recon.py
@@ -0,0 +1,228 @@
+import os
+from typing import Dict, Optional, Tuple
+import tempfile
+
+import numpy as np
+import stir
+
+from mu_map.file import load_as_interfile
+from mu_map.file.interfile import (
+ parse_interfile_header_str,
+ load_interfile,
+ write_interfile,
+ KEY_DIM_1,
+ KEY_DIM_2,
+ KEY_SPACING_1,
+ KEY_SPACING_2,
+ TEMPLATE_HEADER_IMAGE,
+)
+
+
+TEMPLATE_RECON_PARAMS = """
+OSMAPOSLParameters :=
+
+ objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData
+ PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+ input file := {PROJECTION}
+
+ projector pair type := Matrix
+ Projector Pair Using Matrix Parameters :=
+ Matrix type := SPECT UB
+ Projection Matrix By Bin SPECT UB Parameters:=
+ ; width of PSF
+ maximum number of sigmas:= 2.0
+
+ ;PSF type of correction { 2D // 3D // Geometrical }
+ psf type:= Geometrical
+ ; next 2 parameters define the PSF. They are ignored if psf_type is "Geometrical"
+ ; These values are mostly dependent on your collimator.
+ ; the PSF is modelled as a Gaussian with sigma dependent on the distance from the collimator
+ ; sigma_at_depth = collimator_slope * depth_in_cm + collimator sigma 0(cm)
+ collimator slope := 0.0163
+ collimator sigma 0(cm) := 0.1466
+
+ ;Attenuation correction { Simple // Full // No }
+ attenuation type := {ATTENUATION_TYPE}
+ ;Values in attenuation map in cm-1
+ attenuation map := {ATTENUATION_MAP}
+
+ ;Mask properties { Cylinder // Attenuation Map // Explicit Mask // No}
+ mask type := {MASK_TYPE}
+ mask file := {MASK_FILE}
+
+ ; if next variable is set to 0, only a single view is kept in memory
+ keep all views in cache:=1
+
+ End Projection Matrix By Bin SPECT UB Parameters:=
+ End Projector Pair Using Matrix Parameters :=
+ end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:=
+
+ initial estimate := {INIT_FILE}
+ output filename prefix := {OUTPUT_PREFIX}
+
+ number of subsets:= {N_SUBSETS}
+ number of subiterations:= {N_SUBITERATIONS}
+ Save estimates at subiteration intervals:= {SAVE_INTERVALS}
+
+ ; keywords that specify the filtering that occurs after every subiteration
+ ; warning: do not normally use together with a prior
+ ;inter-iteration filter subiteration interval := 4
+ ;inter-iteration filter type := Separable Gaussian
+ post-filter type := Separable Gaussian
+ separable gaussian filter parameters :=
+ x-dir filter fwhm (in mm) := 6
+ y-dir filter fwhm (in mm) := 6
+ z-dir filter fwhm (in mm) := 6
+ x-dir maximum kernel size := 129
+ y-dir maximum kernel size := 129
+ z-dir maximum kernel size := 31
+ Normalise filter to 1 := 1
+ end separable gaussian filter parameters :=
+
+END :=
+"""
+
+
+def uniform_estimate(projection: Tuple[Dict[str, str], np.ndarray]):
+ header_proj, image_proj = projection
+
+ image = np.ones(
+ (image_proj.shape[1], image_proj.shape[2], image_proj.shape[2]), np.float32
+ )
+
+ offset = -0.5 * image_proj.shape[2] * float(header_proj[KEY_SPACING_1])
+ header = TEMPLATE_HEADER_IMAGE.replace("{ROWS}", str(image.shape[2]))
+ header = header.replace("{COLUMNS}", str(image.shape[1]))
+ header = header.replace("{SLICES}", str(image.shape[0]))
+ header = header.replace("{SPACING_X}", header_proj[KEY_SPACING_1])
+ header = header.replace("{SPACING_Y}", header_proj[KEY_SPACING_1])
+ header = header.replace("{SPACING_Z}", header_proj[KEY_SPACING_2])
+ header = header.replace("{OFFSET_X}", f"{offset:.4f}")
+ header = header.replace("{OFFSET_Y}", f"{offset:.4f}")
+ header = parse_interfile_header_str(header)
+
+ return header, image
+
+
+def reconstruct(
+ projection: Tuple[Dict[str, str], np.ndarray],
+ mu_map: Optional[Tuple[Dict[str, str], np.ndarray]] = None,
+ mask: Optional[Tuple[Dict[str, str], np.ndarray]] = None,
+ init: Optional[Tuple[Dict[str, str], np.ndarray]] = None,
+ n_subsets: Optional[int] = 4,
+ n_iterations: Optional[int] = 10,
+):
+ # sanitize parameters
+ n_subiterations = n_subsets * n_iterations
+ save_intervals = n_subiterations
+
+ dir_tmp = tempfile.TemporaryDirectory()
+ filename_projection = os.path.join(dir_tmp.name, "projection.hv")
+ write_interfile(filename_projection, *projection)
+ params = TEMPLATE_RECON_PARAMS.replace("{PROJECTION}", filename_projection)
+
+ output_prefix = os.path.join(dir_tmp.name, "out")
+ filename_out = f"{output_prefix}_{save_intervals}.hv"
+ params = params.replace("{OUTPUT_PREFIX}", output_prefix)
+ params = params.replace("{N_SUBSETS}", str(n_subsets))
+ params = params.replace("{N_SUBITERATIONS}", str(n_subiterations))
+ params = params.replace("{SAVE_INTERVALS}", str(save_intervals))
+
+ if mu_map is not None:
+ filename_mu_map = os.path.join(dir_tmp.name, "mu_map.hv")
+ write_interfile(filename_mu_map, *mu_map)
+ params = params.replace("{ATTENUATION_TYPE}", "Full")
+ params = params.replace("{ATTENUATION_MAP}", filename_mu_map)
+ else:
+ params = params.replace("{ATTENUATION_TYPE}", "No")
+ params = params.replace("attenuation map", ";attenuation map")
+
+ if mask is not None:
+ filename_mask = os.path.join(dir_tmp.name, "mask.hv")
+ write_interfile(filename_mask, *mask)
+ params = params.replace("{MASK_TYPE}", "Explicit Mask")
+ params = params.replace("{MASK_FILE}", filename_mask)
+ else:
+ params = params.replace("mask file", ";mask file")
+ params = params.replace(
+ "{MASK_TYPE}", "Attenuation Map" if mu_map is not None else "No"
+ )
+
+ init = uniform_estimate(projection) if init is None else init
+ filename_init = os.path.join(dir_tmp.name, "init.hv")
+ write_interfile(filename_init, *init)
+ params = params.replace("{INIT_FILE}", filename_init)
+
+ filename_params = os.path.join(dir_tmp.name, "OSEM_SPECT.par")
+ with open(filename_params, mode="w") as f:
+ f.write(params.strip())
+
+ recon = stir.OSMAPOSLReconstruction3DFloat(filename_params)
+ recon.reconstruct()
+
+ print(params)
+ return load_interfile(filename_out)
+
+
+if __name__ == "__main__":
+ import argparse
+
+ from mu_map.file import load_as_interfile
+ from mu_map.recon.project import forward_project
+
+ parser = argparse.ArgumentParser(
+ description="Reconstruct a projection or another reconstruction",
+ formatter_class=argparse.ArgumentDefaultsHelpFormatter,
+ )
+ parser.add_argument(
+ "--projection", type=str, help="a projection in INTERFILE format"
+ )
+ parser.add_argument(
+ "--recon",
+ type=str,
+ help="a reconstruction in DICOM or INTERFILE format - if specified the projection will be overwritten",
+ )
+ parser.add_argument(
+ "--mu_map",
+ type=str,
+ help="a mu map for attenuation correction in DICOM or INTERFILE format",
+ )
+ parser.add_argument(
+ "--out", type=str, help="the filename to store the reconstruction"
+ )
+ parser.add_argument(
+ "--n_subsets",
+ type=int,
+ default=4,
+ help="the number of subsets for OSEM reconstruction",
+ )
+ parser.add_argument(
+ "--n_iterations",
+ type=int,
+ default=10,
+ help="the number of iterations for OSEM reconstruction",
+ )
+ parser.add_argument(
+ "-v", "--verbosity", type=int, default=0, help="configure the verbosity of STIR"
+ )
+
+ args = parser.parse_args()
+ assert (
+ args.projection is not None or args.recon is not None
+ ), "You have to specify either a projection or a reconstruction"
+ stir.Verbosity_set(args.verbosity)
+
+ mu_map = load_as_interfile(args.mu_map) if args.mu_map else None
+ mu_map_slices = None if mu_map is None else mu_map[1].shape[0]
+
+ if args.recon:
+ recon = load_as_interfile(args.recon)
+ projection = forward_project(*recon, n_slices=mu_map_slices)
+ else:
+ projection = load_as_interfile(args.projection)
+
+ header, image = reconstruct(projection, mu_map=mu_map)
+ image = image[:, :, ::-1]
+
+ write_interfile(args.out, header, image)
--
GitLab