diff --git a/mu_map/recon/project.py b/mu_map/recon/project.py
index dc92b54e6ca326b3b232f8ff36d3b3bc02d700cc..fbf76861a1c88337c85d8d45011bee74f074f545 100644
--- a/mu_map/recon/project.py
+++ b/mu_map/recon/project.py
@@ -1,25 +1,22 @@
-import argparse
import os
import tempfile
-from typing import Optional, Tuple
+from typing import Dict, Optional, Tuple
-import cv2 as cv
import numpy as np
from pydicom.dataset import FileDataset as DCMImage
import stir
-import stirextra
-# from vis import to_gray
from mu_map.file.interfile import (
load_interfile,
- KEY_DIM_1,
- KEY_DIM_2,
- KEY_DIM_3,
+ write_interfile,
KEY_SPACING_1,
- KEY_SPACING_2,
+ KEY_SPACING_3,
)
+"""
+A template for a header of a projection in INTERFILE format.
+"""
TEMPLATE_HEADER_PROJ = """
!INTERFILE :=
!imaging modality := nucmed
@@ -56,70 +53,121 @@ radius := {RADIUS}
def forward_project(
- recon: np.ndarray,
- spacing_slices: Optional[float],
- spacing_bins: Optional[float],
+ recon_header: Dict[str, str],
+ recon_img: np.ndarray,
n_slices: Optional[int] = None,
- n_bin: Optional[int] = None,
+ n_bins: Optional[int] = None,
n_projections: int = 120,
rotation: int = 360,
start_angle: int = 180,
radius: int = 150,
- verbose: bool = False,
-):
- n_slices = recon.shape[0] if n_slices is None else n_slices
- n_bins = recon.shape[1] if n_bins is None else n_bins
-
- # dir_tmp = tempfile.TemporaryDirectory()
- header_proj = HEADER_TEMPLATE_PROJ.replace("{SLICES}", str(n_slices))
+ **kwargs,
+) -> Tuple[Dict[str, str], np.ndarray]:
+ """
+ Forward project a reconstruction in INTERFILE format.
+
+ :param recon_header: the header of the reconstruction
+ :param recon_img: the reconstruction image
+ :param n_slices: the number of slices of the projection - defaults to the number of slices of the reconstruction
+ :param n_bins: the number of bins of the projection - defaults to the number of rows in the reconstruction
+ :param n_projections: the number of views used for the projection
+ :param rotation: the total rotation in degrees
+ :param start_angle: the start angle of the rotation
+ :param radius: the radius of the detector circle in mm
+ :return: an interfile image as a tuple of a header and a numpy array
+ """
+ n_slices = recon_img.shape[0] if n_slices is None else n_slices
+ n_bins = recon_img.shape[1] if n_bins is None else n_bins
+ # work around a bug in STIR with start_angle = 0
+ start_angle = 360 if start_angle == 0 else start_angle
+
+ # prepare writing temporary files to configure the projection with STIR
+ filename_recon = "recon.hv"
+ filename_proj = "projection"
+ filename_proj_data = f"{filename_proj}.s"
+ filename_proj_header = f"{filename_proj}.hs"
+
+ header_proj = TEMPLATE_HEADER_PROJ.replace("{DATA_FILE}", filename_proj_data)
+ header_proj = header_proj.replace("{SLICES}", str(n_slices))
header_proj = header_proj.replace("{BINS}", str(n_bins))
- header_proj = header_proj.replace("{SPACING_SLICES}", f"{spacing_slices:.4f}")
- header_proj = header_proj.replace("{SPACING_BINS}", f"{spacing_bins:.4f}")
+ header_proj = header_proj.replace("{SPACING_SLICES}", recon_header[KEY_SPACING_3])
+ header_proj = header_proj.replace("{SPACING_BINS}", recon_header[KEY_SPACING_1])
header_proj = header_proj.replace("{N_PROJECTIONS}", str(n_projections))
header_proj = header_proj.replace("{ROTATION}", str(rotation))
header_proj = header_proj.replace("{START_ANGLE}", str(start_angle))
header_proj = header_proj.replace("{RADIUS}", str(radius))
header_proj = header_proj.strip()
- print("Project header...")
- print(header_proj)
+ # write temporary files to configure the projection with STIR
+ dir_tmp = tempfile.TemporaryDirectory()
-def forward_project_dcm(dcm: DCMImage, recon: np.ndarray, **kwargs):
- spacing_slices = float(dcm.SliceThickness)
- spacing_bins = float(dcm.PixelSpacing[0])
- return forward_project(recon, spacing_slices=spacing_slices, spacing_bins=spacing_bins, **kwargs)
+ filename_proj_header = os.path.join(dir_tmp.name, filename_proj_header)
+ filename_recon = os.path.join(dir_tmp.name, filename_recon)
+ with open(filename_proj_header, mode="w") as f:
+ f.write(header_proj)
+ with open(os.path.join(dir_tmp.name, filename_proj_data), mode="wb") as f:
+ f.write(b"")
+ write_interfile(filename_recon, recon_header, recon_img)
+ projdata = stir.ProjData.read_from_file(filename_proj_header)
+ recon = stir.FloatVoxelsOnCartesianGrid.read_from_file(filename_recon)
-if __name__ == "__main__":
- from mu_map.file.dicom import load_dcm
- dcm, recon = load_dcm("tmp/0001-stress-recon_nac_nsc.dcm")
- forward_project_dcm(dcm, recon)
+ projmatrix = stir.ProjMatrixByBinUsingRayTracing()
+ projmatrix.set_up(projdata.get_proj_data_info(), recon)
+ forwardprojector = stir.ForwardProjectorByBinUsingProjMatrixByBin(projmatrix)
- exit(0)
+ projdataout = stir.ProjDataInMemory(
+ projdata.get_exam_info(), projdata.get_proj_data_info()
+ )
+
+ forwardprojector.set_up(projdataout.get_proj_data_info(), recon)
+ forwardprojector.forward_project(projdataout, recon)
+
+ filename_out = os.path.join(dir_tmp.name, "out.hs")
+ projdataout.write_to_file(filename_out)
+ return load_interfile(filename_out)
+
+
+if __name__ == "__main__":
+ import argparse
+
+ from mu_map.file.dicom import load_dcm
+ from mu_map.file.dicom_to_interfile import to_interfile
parser = argparse.ArgumentParser(
- description="Project a reconstruction with STIR",
+ description="Forward project a reconstruction with STIR",
formatter_class=argparse.ArgumentDefaultsHelpFormatter,
)
parser.add_argument(
"--recon",
type=str,
required=True,
- help="the reconstruction to project in INTERFILE format",
+ help="the reconstruction to project in INTERFILE or DICOM format",
)
parser.add_argument(
"--out", type=str, required=True, help="the projection output by this script"
)
parser.add_argument(
- "--n_projections", type=int, default=120, help="the number of performed projections"
+ "--n_slices",
+ type=int,
+ help="configure a different number of slices than in the reconstruction",
+ )
+ parser.add_argument(
+ "--n_projections",
+ type=int,
+ default=120,
+ help="the number of performed projections",
)
parser.add_argument(
"--rotation", type=int, default=360, help="the amount of rotation in degrees"
)
parser.add_argument(
- "--start_angle", type=int, default=360, help="the start angle for the projections"
+ "--start_angle",
+ type=int,
+ default=180,
+ help="the start angle for the projections",
)
parser.add_argument(
"--radius", type=int, default=150, help="the radius of the rotation"
@@ -127,53 +175,15 @@ if __name__ == "__main__":
parser.add_argument(
"-v", "--verbosity", type=int, default=0, help="configure the verbosity of STIR"
)
- parser.add_argument("--header", type=str, help="configure the verbosity of STIR")
args = parser.parse_args()
- args.out = os.path.splitext(args.out)[0]
- args.start_angle = 360 if args.start_angle == 0 else args.start_angle
stir.Verbosity_set(args.verbosity)
- tmp_dir = tempfile.TemporaryDirectory()
- if args.header is None:
+ _, ext = os.path.splitext(args.recon)
+ if ext == ".dcm":
+ header, image = to_interfile(*load_dcm(args.recon))
+ else:
header, image = load_interfile(args.recon)
- print(header[KEY_DIM_3])
- # header_proj = HEADER_TEMPLATE_PROJ.replace("{ROWS}", header[KEY_DIM_1])
- header_proj = HEADER_TEMPLATE_PROJ.replace("{ROWS}", str(33))
- header_proj = header_proj.replace("{COLUMNS}", header[KEY_DIM_2])
- header_proj = header_proj.replace("{SPACING_X}", header[KEY_SPACING_1])
- header_proj = header_proj.replace("{SPACING_Y}", header[KEY_SPACING_2])
- header_proj = header_proj.replace("{N_PROJECTIONS}", str(args.n_projections))
- header_proj = header_proj.replace("{ROTATION}", str(args.rotation))
- header_proj = header_proj.replace("{START_ANGLE}", str(args.start_angle))
- header_proj = header_proj.replace("{RADIUS}", str(args.radius))
- header_proj = header_proj.strip()
- print(f"Header:\n{header_proj}")
-
- filename = "proj"
- filename_header = f"{filename}.hs"
- filename_data = f"{filename}.s"
- header_proj = header_proj.replace("{DATA_FILE}", filename_data)
- print(header_proj)
-
- with open(os.path.join(tmp_dir.name, filename_header), mode="w") as f:
- f.write(header_proj)
- with open(os.path.join(tmp_dir.name, filename_data), mode="wb") as f:
- f.write(b"")
- args.header = os.path.join(tmp_dir.name, filename_header)
-
- projdata = stir.ProjData.read_from_file(args.header)
- recon = stir.FloatVoxelsOnCartesianGrid.read_from_file(args.recon)
-
- projmatrix = stir.ProjMatrixByBinUsingRayTracing()
- projmatrix.set_up(projdata.get_proj_data_info(), recon)
- forwardprojector = stir.ForwardProjectorByBinUsingProjMatrixByBin(projmatrix)
-
- projdataout = stir.ProjDataInMemory(
- projdata.get_exam_info(), projdata.get_proj_data_info()
- )
-
- forwardprojector.set_up(projdataout.get_proj_data_info(), recon)
- forwardprojector.forward_project(projdataout, recon)
- projdataout.write_to_file(args.out)
+ header, image = forward_project(header, image, **vars(args))
+ write_interfile(args.out, header, image)