From 2733e5d1f4477009fa1a045800f6dce456d33112 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Mon, 13 Nov 2023 17:27:28 -0800 Subject: [PATCH 01/27] first commit --- python/parallelproj_sim.py | 97 +++ python/utils.py | 1192 ++++++++++++++++++++++++++++++++++++ 2 files changed, 1289 insertions(+) create mode 100644 python/parallelproj_sim.py create mode 100644 python/utils.py diff --git a/python/parallelproj_sim.py b/python/parallelproj_sim.py new file mode 100644 index 0000000..318e695 --- /dev/null +++ b/python/parallelproj_sim.py @@ -0,0 +1,97 @@ +from __future__ import annotations + +import utils +import array_api_compat.numpy as np +import matplotlib.pyplot as plt +from array_api_compat import to_device + +# device variable (cpu or cuda) that determines whether calculations +# are performed on the cpu or cuda gpu + +dev = 'cpu' + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +#--- setup the scanner / LOR geometry --------------------------------------- +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- + +# setup a line of response descriptor that describes the LOR start / endpoints of +# a "narrow" clinical PET scanner with 9 rings +lor_descriptor = utils.DemoPETScannerLORDescriptor(np, + dev, + num_rings=4, + radial_trim=141) + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +#--- setup a simple 3D test image ------------------------------------------- +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- + +# image properties +voxel_size = (2.66, 2.66, 2.66) +num_trans = 160 +num_ax = 2 * lor_descriptor.scanner.num_modules + +# setup a box like test image +img_shape = (num_trans, num_trans, num_ax) +n0, n1, n2 = img_shape + +# setup an image containing a box +img = np.zeros(img_shape, dtype=np.float32, device=dev) +img[(n0 // 4):(3 * n0 // 4), (n1 // 4):(3 * n1 // 4), :] = 1 + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +#--- setup a non-TOF projector and project ---------------------------------- +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- + +projector = utils.RegularPolygonPETProjector(lor_descriptor, img_shape, + voxel_size) + +projector.tof = False # set this to True to get a time of flight projector + +img_fwd = projector(img) +back_img = projector.adjoint(img_fwd) + + +# get the start and end points of all sinogram bins +xstart, xend = lor_descriptor.get_lor_coordinates() + +# get the two dimensional indices of all sinogram bins +start_mods, end_mods, start_inds, end_inds = lor_descriptor.get_lor_indices() + +# convert the 2D indices to a 1D index +start_index = lor_descriptor.scanner.num_modules * start_mods + start_inds +end_index = lor_descriptor.scanner.num_modules * end_mods + end_inds + + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- + +# show the scanner geometry and one view in one sinogram plane +fig = plt.figure() +ax = fig.add_subplot(projection='3d') +lor_descriptor.scanner.show_lor_endpoints(ax) +lor_descriptor.show_views(ax, + views=np.asarray([lor_descriptor.num_views // 4], + device=dev), + planes=np.asarray( + [lor_descriptor.scanner.num_modules // 2], + device=dev)) +fig.tight_layout() +fig.show() + +fig2, ax2 = plt.subplots(1, 3, figsize=(15, 5)) +ax2[0].imshow(np.asarray(to_device(img[:, :, 3], 'cpu'))) +if projector.tof: + ax2[1].imshow(np.asarray(to_device(img_fwd[:, :, 4, 15], 'cpu'))) +else: + ax2[1].imshow(np.asarray(to_device(img_fwd[:, :, 4], 'cpu'))) +ax2[2].imshow(np.asarray(to_device(back_img[:, :, 3], 'cpu'))) +fig2.tight_layout() +fig2.show() \ No newline at end of file diff --git a/python/utils.py b/python/utils.py new file mode 100644 index 0000000..9470e75 --- /dev/null +++ b/python/utils.py @@ -0,0 +1,1192 @@ +from __future__ import annotations + +import enum +import abc +from dataclasses import dataclass +import array_api_compat.numpy as np +import numpy.typing as npt +import matplotlib.pyplot as plt +from mpl_toolkits.mplot3d.art3d import Line3DCollection + +from types import ModuleType +from array_api_compat import device, to_device, size + +import parallelproj + + +@dataclass +class TOFParameters: + """ + generic time of flight (TOF) parameters for a scanner with 385ps FWHM TOF + + num_tofbins: int + number of time of flight bins + tofbin_width: float + width of the TOF bin in spatial units (mm) + sigma_tof: float + standard deviation of Gaussian TOF kernel in spatial units (mm) + num_sigmas: float + number of sigmas after which TOF kernel is truncated + tofcenter_offset: float + offset of center of central TOF bin from LOR center in spatial units (mm) + """ + num_tofbins: int = 29 + tofbin_width: float = 13 * 0.01302 * 299.792 / 2 # 13 TOF "small" TOF bins of 0.01302[ns] * (speed of light / 2) [mm/ns] + sigma_tof: float = (299.792 / 2) * ( + 0.385 / 2.355) # (speed_of_light [mm/ns] / 2) * TOF FWHM [ns] / 2.355 + num_sigmas: float = 3. + tofcenter_offset: float = 0 + + +class SinogramSpatialAxisOrder(enum.Enum): + """order of spatial axis in a sinogram R (radial), V (view), P (plane)""" + + RVP = enum.auto() + """[radial,view,plane]""" + RPV = enum.auto() + """[radial,plane,view]""" + VRP = enum.auto() + """[view,radial,plane]""" + VPR = enum.auto() + """[view,plane,radial]""" + PRV = enum.auto() + """[plane,radial,view]""" + PVR = enum.auto() + """[plane,view,radial]""" + + +class PETScannerModule(abc.ABC): + + def __init__( + self, + xp: ModuleType, + dev: str, + num_lor_endpoints: int, + affine_transformation_matrix: npt.NDArray | None = None) -> None: + """abstract base class for PET scanner module + + Parameters + ---------- + xp: ModuleType + array module to use for storing the LOR endpoints + dev: str + device to use for storing the LOR endpoints + num_lor_endpoints : int + number of LOR endpoints in the module + affine_transformation_matrix : npt.NDArray | None, optional + 4x4 affine transformation matrix applied to the LOR endpoint coordinates, default None + if None, the 4x4 identity matrix is used + """ + + self._xp = xp + self._dev = dev + self._num_lor_endpoints = num_lor_endpoints + self._lor_endpoint_numbers = self.xp.arange(num_lor_endpoints, + device=self.dev) + + if affine_transformation_matrix is None: + self._affine_transformation_matrix = self.xp.eye(4, + device=self.dev) + else: + self._affine_transformation_matrix = affine_transformation_matrix + + @property + def xp(self) -> ModuleType: + """array module to use for storing the LOR endpoints""" + return self._xp + + @property + def dev(self) -> str: + """device to use for storing the LOR endpoints""" + return self._dev + + @property + def num_lor_endpoints(self) -> int: + """total number of LOR endpoints in the module + + Returns + ------- + int + """ + return self._num_lor_endpoints + + @property + def lor_endpoint_numbers(self) -> npt.NDArray: + """array enumerating all the LOR endpoints in the module + + Returns + ------- + npt.NDArray + """ + return self._lor_endpoint_numbers + + @property + def affine_transformation_matrix(self) -> npt.NDArray: + """4x4 affine transformation matrix + + Returns + ------- + npt.NDArray + """ + return self._affine_transformation_matrix + + @abc.abstractmethod + def get_raw_lor_endpoints(self, + inds: npt.NDArray | None = None) -> npt.NDArray: + """mapping from LOR endpoint indices within module to an array of "raw" world coordinates + + Parameters + ---------- + inds : npt.NDArray | None, optional + an non-negative integer array of indices, default None + if None means all possible indices [0, ... , num_lor_endpoints - 1] + + Returns + ------- + npt.NDArray + a 3 x len(inds) float array with the world coordinates of the LOR endpoints + """ + if inds is None: + inds = self.lor_endpoint_numbers + raise NotImplementedError + + def get_lor_endpoints(self, + inds: npt.NDArray | None = None) -> npt.NDArray: + """mapping from LOR endpoint indices within module to an array of "transformed" world coordinates + + Parameters + ---------- + inds : npt.NDArray | None, optional + an non-negative integer array of indices, default None + if None means all possible indices [0, ... , num_lor_endpoints - 1] + + Returns + ------- + npt.NDArray + a 3 x len(inds) float array with the world coordinates of the LOR endpoints including an affine transformation + """ + + raw_lor_endpoints = self.get_raw_lor_endpoints(inds) + + tmp = self.xp.ones((raw_lor_endpoints.shape[0], 4), device=self.dev) + tmp[:, :-1] = raw_lor_endpoints + + return (tmp @ self.affine_transformation_matrix.T)[:, :3] + + def show_lor_endpoints(self, + ax: plt.Axes, + annotation_fontsize: float = 0, + annotation_prefix: str = '', + annotation_offset: int = 0, + transformed: bool = True, + **kwargs) -> None: + """show the LOR coordinates in a 3D scatter plot + + Parameters + ---------- + ax : plt.Axes + 3D matplotlib axes + annotation_fontsize : float, optional + fontsize of LOR endpoint number annotation, by default 0 + annotation_prefix : str, optional + prefix for annotation, by default '' + annotation_offset : int, optional + number to add to crystal number, by default 0 + transformed : bool, optional + use transformed instead of raw coordinates, by default True + """ + + if transformed: + all_lor_endpoints = self.get_lor_endpoints() + else: + all_lor_endpoints = self.get_raw_lor_endpoints() + + # convert to numpy array + all_lor_endpoints = np.asarray(to_device(all_lor_endpoints, 'cpu')) + + ax.scatter(all_lor_endpoints[:, 0], all_lor_endpoints[:, 1], + all_lor_endpoints[:, 2], **kwargs) + + ax.set_box_aspect([ + ub - lb for lb, ub in (getattr(ax, f'get_{a}lim')() for a in 'xyz') + ]) + + ax.set_xlabel('x0') + ax.set_ylabel('x1') + ax.set_zlabel('x2') + + if annotation_fontsize > 0: + for i in self.lor_endpoint_numbers: + ax.text(all_lor_endpoints[i, 0], + all_lor_endpoints[i, 1], + all_lor_endpoints[i, 2], + f'{annotation_prefix}{i+annotation_offset}', + fontsize=annotation_fontsize) + + +class RegularPolygonPETScannerModule(PETScannerModule): + + def __init__( + self, + xp: ModuleType, + dev: str, + radius: float, + num_sides: int, + num_lor_endpoints_per_side: int, + lor_spacing: float, + ax0: int = 2, + ax1: int = 1, + affine_transformation_matrix: npt.NDArray | None = None) -> None: + """regular Polygon PET scanner module + + Parameters + ---------- + xp: ModuleType + array module to use for storing the LOR endpoints + device: str + device to use for storing the LOR endpoints + radius : float + inner radius of the regular polygon + num_sides: int + number of sides of the regular polygon + num_lor_endpoints_per_sides: int + number of LOR endpoints per side + lor_spacing : float + spacing between the LOR endpoints in the polygon direction + ax0 : int, optional + axis number for the first direction, by default 2 + ax1 : int, optional + axis number for the second direction, by default 1 + affine_transformation_matrix : npt.NDArray | None, optional + 4x4 affine transformation matrix applied to the LOR endpoint coordinates, default None + if None, the 4x4 identity matrix is used + """ + + self._radius = radius + self._num_sides = num_sides + self._num_lor_endpoints_per_side = num_lor_endpoints_per_side + self._ax0 = ax0 + self._ax1 = ax1 + self._lor_spacing = lor_spacing + super().__init__(xp, dev, num_sides * num_lor_endpoints_per_side, + affine_transformation_matrix) + + @property + def radius(self) -> float: + """inner radius of the regular polygon + + Returns + ------- + float + """ + return self._radius + + @property + def num_sides(self) -> int: + """number of sides of the regular polygon + + Returns + ------- + int + """ + return self._num_sides + + @property + def num_lor_endpoints_per_side(self) -> int: + """number of LOR endpoints per side + + Returns + ------- + int + """ + return self._num_lor_endpoints_per_side + + @property + def ax0(self) -> int: + """axis number for the first module direction + + Returns + ------- + int + """ + return self._ax0 + + @property + def ax1(self) -> int: + """axis number for the second module direction + + Returns + ------- + int + """ + return self._ax1 + + @property + def lor_spacing(self) -> float: + """spacing between the LOR endpoints in a module along the polygon + + Returns + ------- + float + """ + return self._lor_spacing + + # abstract method from base class to be implemented + def get_raw_lor_endpoints(self, + inds: npt.NDArray | None = None) -> npt.NDArray: + if inds is None: + inds = self.lor_endpoint_numbers + + side = inds // self.num_lor_endpoints_per_side + tmp = inds - side * self.num_lor_endpoints_per_side + tmp = self.xp.astype( + tmp, float) - (self.num_lor_endpoints_per_side / 2 - 0.5) + + phi = 2 * self.xp.pi * self.xp.astype(side, float) / self.num_sides + + lor_endpoints = self.xp.zeros((self.num_lor_endpoints, 3), + device=self.dev) + lor_endpoints[:, self.ax0] = self.xp.cos( + phi) * self.radius - self.xp.sin(phi) * self.lor_spacing * tmp + lor_endpoints[:, self.ax1] = self.xp.sin( + phi) * self.radius + self.xp.cos(phi) * self.lor_spacing * tmp + + return lor_endpoints + + +class ModularizedPETScannerGeometry: + """description of a PET scanner geometry consisting of LOR endpoint modules""" + + def __init__(self, modules: tuple[PETScannerModule]): + """ + Parameters + ---------- + modules : tuple[PETScannerModule] + a tuple of scanner modules + """ + + # member variable that determines whether we want to use + # a numpy or cupy array to store the array of all lor endpoints + self._modules = modules + self._num_modules = len(self._modules) + self._num_lor_endpoints_per_module = self.xp.asarray( + [x.num_lor_endpoints for x in self._modules], device=self.dev) + self._num_lor_endpoints = int( + self.xp.sum(self._num_lor_endpoints_per_module)) + + self.setup_all_lor_endpoints() + + def setup_all_lor_endpoints(self) -> None: + """calculate the position of all lor endpoints by iterating over + the modules and calculating the transformed coordinates of all + module endpoints + """ + + self._all_lor_endpoints_index_offset = self.xp.asarray([ + int(sum(self._num_lor_endpoints_per_module[:i])) + for i in range(size(self._num_lor_endpoints_per_module)) + ], + device=self.dev) + + self._all_lor_endpoints = self.xp.zeros((self._num_lor_endpoints, 3), + device=self.dev, + dtype=self.xp.float32) + + for i, module in enumerate(self._modules): + self._all_lor_endpoints[ + int(self._all_lor_endpoints_index_offset[i]):int( + self._all_lor_endpoints_index_offset[i] + + module.num_lor_endpoints), :] = module.get_lor_endpoints() + + self._all_lor_endpoints_module_number = [ + int(self._num_lor_endpoints_per_module[i]) * [i] + for i in range(self._num_modules) + ] + + self._all_lor_endpoints_module_number = self.xp.asarray( + [i for r in self._all_lor_endpoints_module_number for i in r], + device=self.dev) + + @property + def modules(self) -> tuple[PETScannerModule]: + """tuple of modules defining the scanner""" + return self._modules + + @property + def num_modules(self) -> int: + """the number of modules defining the scanner""" + return self._num_modules + + @property + def num_lor_endpoints_per_module(self) -> npt.NDArray: + """numpy array showing how many LOR endpoints are in every module""" + return self._num_lor_endpoints_per_module + + @property + def num_lor_endpoints(self) -> int: + """the total number of LOR endpoints in the scanner""" + return self._num_lor_endpoints + + @property + def all_lor_endpoints_index_offset(self) -> npt.NDArray: + """the offset in the linear (flattend) index for all LOR endpoints""" + return self._all_lor_endpoints_index_offset + + @property + def all_lor_endpoints_module_number(self) -> npt.NDArray: + """the module number of all LOR endpoints""" + return self._all_lor_endpoints_module_number + + @property + def all_lor_endpoints(self) -> npt.NDArray: + """the world coordinates of all LOR endpoints""" + return self._all_lor_endpoints + + @property + def xp(self) -> ModuleType: + """module indicating whether the LOR endpoints are stored as numpy or cupy array""" + return self._modules[0].xp + + @property + def dev(self) -> str: + return self._modules[0].dev + + def linear_lor_endpoint_index( + self, + module: npt.NDArray, + index_in_module: npt.NDArray, + ) -> npt.NDArray: + """transform the module + index_in_modules indices into a flattened / linear LOR endpoint index + + Parameters + ---------- + module : npt.NDArray + containing module numbers + index_in_module : npt.NDArray + containing index in modules + + Returns + ------- + npt.NDArray + the flattened LOR endpoint index + """ + # index_in_module = self._xp.asarray(index_in_module) + + return self.xp.take(self.all_lor_endpoints_index_offset, + module) + index_in_module + + def get_lor_endpoints(self, module: npt.NDArray, + index_in_module: npt.NDArray) -> npt.NDArray: + """get the coordinates for LOR endpoints defined by module and index in module + + Parameters + ---------- + module : npt.NDArray + the module number of the LOR endpoints + index_in_module : npt.NDArray + the index in module number of the LOR endpoints + + Returns + ------- + npt.NDArray | cpt.NDArray + the 3 world coordinates of the LOR endpoints + """ + return self.xp.take(self.all_lor_endpoints, + self.linear_lor_endpoint_index( + module, index_in_module), + axis=0) + + def show_lor_endpoints(self, + ax: plt.Axes, + show_linear_index: bool = True, + **kwargs) -> None: + """show all LOR endpoints in a 3D plot + + Parameters + ---------- + ax : plt.Axes + a 3D matplotlib axes + show_linear_index : bool, optional + annotate the LOR endpoints with the linear LOR endpoint index + **kwargs : keyword arguments + passed to show_lor_endpoints() of the scanner module + """ + for i, module in enumerate(self.modules): + if show_linear_index: + offset = np.asarray( + to_device(self.all_lor_endpoints_index_offset[i], 'cpu')) + prefix = f'' + else: + offset = 0 + prefix = f'{i},' + + module.show_lor_endpoints(ax, + annotation_offset=offset, + annotation_prefix=prefix, + **kwargs) + + +class RegularPolygonPETScannerGeometry(ModularizedPETScannerGeometry): + """description of a PET scanner geometry consisting stacked regular polygons""" + + def __init__(self, xp: ModuleType, dev: str, radius: float, num_sides: int, + num_lor_endpoints_per_side: int, lor_spacing: float, + num_rings: int, ring_positions: npt.NDArray, + symmetry_axis: int) -> None: + """ + Parameters + ---------- + xp: ModuleType + array module to use for storing the LOR endpoints + dev: str + device to use for storing the LOR endpoints + radius : float + radius of the scanner + num_sides : int + number of sides (faces) of each regular polygon + num_lor_endpoints_per_side : int + number of LOR endpoints in each side (face) of each polygon + lor_spacing : float + spacing between the LOR endpoints in each side + num_rings : int + the number of rings (regular polygons) + ring_positions : npt.NDArray + 1D array with the coordinate of the rings along the ring axis + symmetry_axis : int + the ring axis (0,1,2) + """ + + self._radius = radius + self._num_sides = num_sides + self._num_lor_endpoints_per_side = num_lor_endpoints_per_side + self._num_rings = num_rings + self._lor_spacing = lor_spacing + self._symmetry_axis = symmetry_axis + self._ring_positions = ring_positions + + if symmetry_axis == 0: + self._ax0 = 2 + self._ax1 = 1 + elif symmetry_axis == 1: + self._ax0 = 0 + self._ax1 = 2 + elif symmetry_axis == 2: + self._ax0 = 1 + self._ax1 = 0 + + modules = [] + + for ring in range(num_rings): + aff_mat = xp.eye(4, device=dev) + aff_mat[symmetry_axis, -1] = ring_positions[ring] + + modules.append( + RegularPolygonPETScannerModule( + xp, + dev, + radius, + num_sides, + num_lor_endpoints_per_side=num_lor_endpoints_per_side, + lor_spacing=lor_spacing, + affine_transformation_matrix=aff_mat, + ax0=self._ax0, + ax1=self._ax1)) + + modules = tuple(modules) + super().__init__(modules) + + self._all_lor_endpoints_index_in_ring = self.xp.arange( + self.num_lor_endpoints, device=dev + ) - self.all_lor_endpoints_ring_number * self.num_lor_endpoints_per_module[ + 0] + + @property + def radius(self) -> float: + """radius of the scanner""" + return self._radius + + @property + def num_sides(self) -> int: + """number of sides (faces) of each polygon""" + return self._num_sides + + @property + def num_lor_endpoints_per_side(self) -> int: + """number of LOR endpoints per side (face) in each polygon""" + return self._num_lor_endpoints_per_side + + @property + def num_rings(self) -> int: + """number of rings (regular polygons)""" + return self._num_rings + + @property + def lor_spacing(self) -> float: + """the spacing between the LOR endpoints in every side (face) of each polygon""" + return self._lor_spacing + + @property + def symmetry_axis(self) -> int: + """The symmetry axis. Also called axial (or ring) direction.""" + return self._symmetry_axis + + @property + def all_lor_endpoints_ring_number(self) -> npt.NDArray: + """the ring (regular polygon) number of all LOR endpoints""" + return self._all_lor_endpoints_module_number + + @property + def all_lor_endpoints_index_in_ring(self) -> npt.NDArray: + """the index withing the ring (regular polygon) number of all LOR endpoints""" + return self._all_lor_endpoints_index_in_ring + + @property + def num_lor_endpoints_per_ring(self) -> int: + """the number of LOR endpoints per ring (regular polygon)""" + return int(self._num_lor_endpoints_per_module[0]) + + @property + def ring_positions(self) -> npt.NDArray: + """the ring (regular polygon) positions""" + return self._ring_positions + + +class DemoPETScanner(RegularPolygonPETScannerGeometry): + + def __init__(self, + xp: ModuleType, + dev: str, + num_rings: int = 36, + symmetry_axis: int = 2) -> None: + + ring_positions = 5.32 * xp.arange( + num_rings, device=dev, dtype=xp.float32) + (xp.astype( + xp.arange(num_rings, device=dev) // 9, xp.float32)) * 2.8 + ring_positions -= 0.5 * xp.max(ring_positions) + super().__init__(xp, + dev, + radius=0.5 * (744.1 + 2 * 8.51), + num_sides=34, + num_lor_endpoints_per_side=16, + lor_spacing=4.03125, + num_rings=num_rings, + ring_positions=ring_positions, + symmetry_axis=symmetry_axis) + + +class PETLORDescriptor(abc.ABC): + """abstract base class to describe which modules / indices in modules of a + modularized PET scanner are in coincidence; defining geometrical LORs""" + + def __init__(self, scanner: ModularizedPETScannerGeometry) -> None: + """ + Parameters + ---------- + scanner : ModularizedPETScannerGeometry + a modularized PET scanner + """ + self._scanner = scanner + + @abc.abstractmethod + def get_lor_coordinates(self, + **kwargs) -> tuple[npt.ArrayLike, npt.ArrayLike]: + """return the start and end coordinates of all (or a subset of) LORs""" + raise NotImplementedError + + @property + def scanner(self) -> ModularizedPETScannerGeometry: + """the scanner for which coincidences are described""" + return self._scanner + + @property + def xp(self) -> ModuleType: + """array module to use for storing the LOR endpoints""" + return self.scanner.xp + + @property + def dev(self) -> str: + """device to use for storing the LOR endpoints""" + return self.scanner.dev + + +class RegularPolygonPETLORDescriptor(PETLORDescriptor): + + def __init__( + self, + scanner: RegularPolygonPETScannerGeometry, + radial_trim: int = 3, + max_ring_difference: int | None = None, + ) -> None: + """Coincidence descriptor for a regular polygon PET scanner where + we have coincidences within and between "rings (polygons of modules)" + The geometrical LORs can be sorted into a sinogram having a + "plane", "view" and "radial" axis. + + Parameters + ---------- + scanner : RegularPolygonPETScannerGeometry + a regular polygon PET scanner + radial_trim : int, optional + number of geometrial LORs to disregard in the radial direction, by default 3 + max_ring_difference : int | None, optional + maximim ring difference to consider for coincidences, by default None means + all ring differences are included + """ + + super().__init__(scanner) + + self._radial_trim = radial_trim + + if max_ring_difference is None: + self._max_ring_difference = self.scanner.num_rings - 1 + else: + self._max_ring_difference = max_ring_difference + + self._num_rad = (self.scanner.num_lor_endpoints_per_ring + + 1) - 2 * self._radial_trim + self._num_views = self.scanner.num_lor_endpoints_per_ring // 2 + + self._setup_plane_indices() + self._setup_view_indices() + + @property + def radial_trim(self) -> int: + """number of geometrial LORs to disregard in the radial direction""" + return self._radial_trim + + @property + def max_ring_difference(self) -> int: + """the maximum ring difference""" + return self._max_ring_difference + + @property + def num_planes(self) -> int: + """number of planes in the sinogram""" + return self._num_planes + + @property + def num_rad(self) -> int: + """number of radial elements in the sinogram""" + return self._num_rad + + @property + def num_views(self) -> int: + """number of views in the sinogram""" + return self._num_views + + @property + def start_plane_index(self) -> npt.NDArray: + """start plane for all planes""" + return self._start_plane_index + + @property + def end_plane_index(self) -> npt.NDArray: + """end plane for all planes""" + return self._end_plane_index + + @property + def start_in_ring_index(self) -> npt.NDArray: + """start index within ring for all views - shape (num_view, num_rad)""" + return self._start_in_ring_index + + @property + def end_in_ring_index(self) -> npt.NDArray: + """end index within ring for all views - shape (num_view, num_rad)""" + return self._end_in_ring_index + + def _setup_plane_indices(self) -> None: + """setup the start / end plane indices (similar to a Michelogram) + """ + self._start_plane_index = self.xp.arange(self.scanner.num_rings, + dtype=self.xp.int32, + device=self.dev) + self._end_plane_index = self.xp.arange(self.scanner.num_rings, + dtype=self.xp.int32, + device=self.dev) + + for i in range(1, self._max_ring_difference + 1): + tmp1 = self.xp.arange(self.scanner.num_rings - i, + dtype=self.xp.int16, + device=self.dev) + tmp2 = self.xp.arange(self.scanner.num_rings - i, + dtype=self.xp.int16, + device=self.dev) + i + + self._start_plane_index = self.xp.concat( + (self._start_plane_index, tmp1, tmp2)) + self._end_plane_index = self.xp.concat( + (self._end_plane_index, tmp2, tmp1)) + + self._num_planes = self._start_plane_index.shape[0] + + def _setup_view_indices(self) -> None: + """setup the start / end view indices + """ + n = self.scanner.num_lor_endpoints_per_ring + + m = 2 * (n // 2) + + self._start_in_ring_index = self.xp.zeros( + (self._num_views, self._num_rad), + dtype=self.xp.int32, + device=self.dev) + self._end_in_ring_index = self.xp.zeros( + (self._num_views, self._num_rad), + dtype=self.xp.int32, + device=self.dev) + + for view in np.arange(self._num_views): + self._start_in_ring_index[view, :] = ( + self.xp.concat( + (self.xp.arange(m) // 2, self.xp.asarray([n // 2]))) - + view)[self._radial_trim:-self._radial_trim] + self._end_in_ring_index[view, :] = ( + self.xp.concat( + (self.xp.asarray([-1]), -((self.xp.arange(m) + 4) // 2))) - + view)[self._radial_trim:-self._radial_trim] + + # shift the negative indices + self._start_in_ring_index = self.xp.where( + self._start_in_ring_index >= 0, self._start_in_ring_index, + self._start_in_ring_index + n) + self._end_in_ring_index = self.xp.where(self._end_in_ring_index >= 0, + self._end_in_ring_index, + self._end_in_ring_index + n) + + def get_lor_indices( + self, + views: None | npt.ArrayLike = None, + sinogram_order=SinogramSpatialAxisOrder.RVP + ) -> tuple[npt.ArrayLike, npt.ArrayLike, npt.ArrayLike, npt.ArrayLike]: + """return the start and end indices of all LORs / or a subset of views + + Parameters + ---------- + views : None | npt.ArrayLike, optional + the views to consider, by default None means all views + sinogram_order : SinogramSpatialAxisOrder, optional + the order of the sinogram axes, by default SinogramSpatialAxisOrder.RVP + + Returns + ------- + start_mods, end_mods, start_inds, end_inds + """ + + if views is None: + views = self.xp.arange(self.num_views, device=self.dev) + + # setup the module and in_module (in_ring) indices for all LORs in PVR order + start_inring_inds = self.xp.reshape( + self.xp.take(self.start_in_ring_index, views, axis=0), (-1, )) + end_inring_inds = self.xp.reshape( + self.xp.take(self.end_in_ring_index, views, axis=0), (-1, )) + + start_mods, start_inds = self.xp.meshgrid(self.start_plane_index, + start_inring_inds, + indexing='ij') + end_mods, end_inds = self.xp.meshgrid(self.end_plane_index, + end_inring_inds, + indexing='ij') + + # reshape to PVR dimensions (radial moving fastest, planes moving slowest) + sinogram_spatial_shape = (self.num_planes, views.shape[0], + self.num_rad) + start_mods = self.xp.reshape(start_mods, sinogram_spatial_shape) + end_mods = self.xp.reshape(end_mods, sinogram_spatial_shape) + start_inds = self.xp.reshape(start_inds, sinogram_spatial_shape) + end_inds = self.xp.reshape(end_inds, sinogram_spatial_shape) + + if sinogram_order is not SinogramSpatialAxisOrder.PVR: + if sinogram_order is SinogramSpatialAxisOrder.RVP: + new_order = (2, 1, 0) + elif sinogram_order is SinogramSpatialAxisOrder.RPV: + new_order = (2, 0, 1) + elif sinogram_order is SinogramSpatialAxisOrder.VRP: + new_order = (1, 2, 0) + elif sinogram_order is SinogramSpatialAxisOrder.VPR: + new_order = (1, 0, 2) + elif sinogram_order is SinogramSpatialAxisOrder.PRV: + new_order = (0, 2, 1) + + start_mods = self.xp.permute_dims(start_mods, new_order) + end_mods = self.xp.permute_dims(end_mods, new_order) + + start_inds = self.xp.permute_dims(start_inds, new_order) + end_inds = self.xp.permute_dims(end_inds, new_order) + + return start_mods, end_mods, start_inds, end_inds + + def get_lor_coordinates( + self, + views: None | npt.ArrayLike = None, + sinogram_order=SinogramSpatialAxisOrder.RVP + ) -> tuple[npt.ArrayLike, npt.ArrayLike]: + + """return the start and end coordinates of all LORs / or a subset of views + + Parameters + ---------- + views : None | npt.ArrayLike, optional + the views to consider, by default None means all views + sinogram_order : SinogramSpatialAxisOrder, optional + the order of the sinogram axes, by default SinogramSpatialAxisOrder.RVP + + Returns + ------- + xstart, xend : npt.ArrayLike + 2 dimensional floating point arrays containing the start and end coordinates of all LORs + """ + + start_mods, end_mods, start_inds, end_inds = self.get_lor_indices(views, sinogram_order) + sinogram_spatial_shape = start_mods.shape + + start_mods = self.xp.reshape(start_mods, (-1, )) + start_inds = self.xp.reshape(start_inds, (-1, )) + + end_mods = self.xp.reshape(end_mods, (-1, )) + end_inds = self.xp.reshape(end_inds, (-1, )) + + x_start = self.xp.reshape( + self.scanner.get_lor_endpoints(start_mods, start_inds), + sinogram_spatial_shape + (3, )) + x_end = self.xp.reshape( + self.scanner.get_lor_endpoints(end_mods, end_inds), + sinogram_spatial_shape + (3, )) + + return x_start, x_end + + def show_views(self, + ax: plt.Axes, + views: npt.ArrayLike, + planes: npt.ArrayLike, + lw: float = 0.2, + **kwargs) -> None: + """show all LORs of a single view in a given plane + + Parameters + ---------- + ax : plt.Axes + a 3D matplotlib axes + view : int + the view number + plane : int + the plane number + lw : float, optional + the line width, by default 0.2 + """ + + xs, xe = self.get_lor_coordinates( + views=views, sinogram_order=SinogramSpatialAxisOrder.RVP) + xs = self.xp.reshape(self.xp.take(xs, planes, axis=2), (-1, 3)) + xe = self.xp.reshape(self.xp.take(xe, planes, axis=2), (-1, 3)) + + p1s = np.asarray(to_device(xs, 'cpu')) + p2s = np.asarray(to_device(xe, 'cpu')) + + ls = np.hstack([p1s, p2s]).copy() + ls = ls.reshape((-1, 2, 3)) + lc = Line3DCollection(ls, linewidths=lw, **kwargs) + ax.add_collection(lc) + + +class DemoPETScannerLORDescriptor(RegularPolygonPETLORDescriptor): + + def __init__(self, + xp: ModuleType, + dev: str, + num_rings: int = 9, + radial_trim: int = 65, + max_ring_difference: int | None = None, + symmetry_axis: int = 2) -> None: + + scanner = DemoPETScanner(xp, + dev, + num_rings, + symmetry_axis=symmetry_axis) + + super().__init__(scanner, + radial_trim=radial_trim, + max_ring_difference=max_ring_difference) + + +class RegularPolygonPETProjector(parallelproj.LinearOperator): + + def __init__(self, + lor_descriptor: RegularPolygonPETLORDescriptor, + img_shape: tuple[int, int, int], + voxel_size: tuple[float, float, float], + img_origin: None | npt.ArrayLike = None, + views: None | npt.ArrayLike = None, + resolution_model: None | parallelproj.LinearOperator = None, + tof: bool = False): + """Regular polygon PET projector + + Parameters + ---------- + lor_descriptor : RegularPolygonPETLORDescriptor + descriptor of the LOR start / end points + img_shape : tuple[int, int, int] + shape of the image to be projected + voxel_size : tuple[float, float, float] + the voxel size of the image to be projected + img_origin : None | npt.ArrayLike, optional + the origin of the image to be projected, by default None + means that image is "centered" in the scanner + views : None | npt.ArrayLike, optional + sinogram views to be projected, by default None + means that all views are being projected + resolution_model : None | parallelproj.LinearOperator, optional + an image-based resolution model applied before forward projection, by default None + means an isotropic 4.5mm FWHM Gaussian smoothing is used + tof: bool, optional, default False + whether to use non-TOF or TOF projections + """ + + super().__init__() + self._dev = lor_descriptor.dev + + self._lor_descriptor = lor_descriptor + self._img_shape = img_shape + self._voxel_size = self.xp.asarray(voxel_size, + dtype=self.xp.float32, + device=self._dev) + + if img_origin is None: + self._img_origin = (-(self.xp.asarray( + self._img_shape, dtype=self.xp.float32, device=self._dev) / 2) + + 0.5) * self._voxel_size + else: + self._img_origin = self.xp.asarray(img_origin, + dtype=self.xp.float32, + device=self._dev) + + if views is None: + self._views = self.xp.arange(self._lor_descriptor.num_views, + device=self._dev) + else: + self._views = views + + if resolution_model is None: + self._resolution_model = parallelproj.GaussianFilterOperator( + self.in_shape, sigma=4.5 / (2.355 * self._voxel_size)) + else: + self._resolution_model = resolution_model + + self._xstart, self._xend = lor_descriptor.get_lor_coordinates( + views=self._views, sinogram_order=SinogramSpatialAxisOrder['RVP']) + + self._tof = tof + self._tof_parameters = TOFParameters() + + @property + def in_shape(self) -> tuple[int, int, int]: + return self._img_shape + + @property + def out_shape(self) -> tuple[int, int, int]: + if self.tof: + out_shape = (self._lor_descriptor.num_rad, self._views.shape[0], + self._lor_descriptor.num_planes, + self.tof_parameters.num_tofbins) + else: + out_shape = (self._lor_descriptor.num_rad, self._views.shape[0], + self._lor_descriptor.num_planes) + + return out_shape + + @property + def xp(self) -> ModuleType: + return self._lor_descriptor.xp + + @property + def tof(self) -> bool: + return self._tof + + @tof.setter + def tof(self, value: bool) -> None: + if not isinstance(value, bool): + raise ValueError('tof must be a boolean') + self._tof = value + + @property + def tof_parameters(self) -> TOFParameters: + return self._tof_parameters + + @tof_parameters.setter + def tof_parameters(self, value: TOFParameters) -> None: + if not isinstance(value, TOFParameters): + raise ValueError('tof_parameters must be a TOFParameters object') + self._tof_parameters = value + + def _apply(self, x): + """nonTOF forward projection of input image x including image based resolution model""" + + dev = device(x) + x_sm = self._resolution_model(x) + + if not self.tof: + x_fwd = parallelproj.joseph3d_fwd(self._xstart, self._xend, x_sm, + self._img_origin, + self._voxel_size) + else: + x_fwd = parallelproj.joseph3d_fwd_tof_sino( + self._xstart, self._xend, x_sm, self._img_origin, + self._voxel_size, self._tof_parameters.tofbin_width, + self.xp.asarray([self._tof_parameters.sigma_tof], + dtype=self.xp.float32, + device=dev), + self.xp.asarray([self._tof_parameters.tofcenter_offset], + dtype=self.xp.float32, + device=dev), self.tof_parameters.num_sigmas, + self.tof_parameters.num_tofbins) + + return x_fwd + + def _adjoint(self, y): + """nonTOF back projection of sinogram y""" + dev = device(y) + + if not self.tof: + y_back = parallelproj.joseph3d_back(self._xstart, self._xend, + self._img_shape, + self._img_origin, + self._voxel_size, y) + else: + y_back = parallelproj.joseph3d_back_tof_sino( + self._xstart, self._xend, self._img_shape, self._img_origin, + self._voxel_size, y, self._tof_parameters.tofbin_width, + self.xp.asarray([self._tof_parameters.sigma_tof], + dtype=self.xp.float32, + device=dev), + self.xp.asarray([self._tof_parameters.tofcenter_offset], + dtype=self.xp.float32, + device=dev), self.tof_parameters.num_sigmas, + self.tof_parameters.num_tofbins) + + return self._resolution_model.adjoint(y_back) + + + +def distributed_subset_order(n: int) -> list[int]: + """subset order that maximizes distance between subsets + + Parameters + ---------- + n : int + number of subsets + + Returns + ------- + list[int] + """ + l = [x for x in range(n)] + o = [] + + for i in range(n): + if (i % 2) == 0: + o.append(l.pop(0)) + else: + o.append(l.pop(len(l)//2)) + + return o + From 30d447156d0dac6da580cbf5820c4477dda5d725 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Mon, 13 Nov 2023 21:54:18 -0800 Subject: [PATCH 02/27] add parallelproj to env --- environment.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/environment.yml b/environment.yml index a61c4e7..cbd25bf 100644 --- a/environment.yml +++ b/environment.yml @@ -18,3 +18,4 @@ dependencies: - shellcheck>=0.8.0 - xtensor-fftw>=0.2.5 - xtensor>=0.24.2 + - parallelproj>=1.6.1 From ae4630ad97d6208377bbf1fd9d58b6f4f860def2 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Mon, 13 Nov 2023 22:21:06 -0800 Subject: [PATCH 03/27] add scipy dependency and correct plot --- environment.yml | 3 +- python/parallelproj_sim.py | 75 ++++++++++++++++++-------------------- 2 files changed, 38 insertions(+), 40 deletions(-) diff --git a/environment.yml b/environment.yml index cbd25bf..f1ed0e4 100644 --- a/environment.yml +++ b/environment.yml @@ -10,7 +10,7 @@ dependencies: - h5py>=3.7.0 - hdf5>=1.12.1 - howardhinnant_date>=3.0.1 - - ipykernel>=6.19.2 + - ipykernel>=6.19.2 - ninja>=1.11.0 - nlohmann_json>=3.11.2 - numpy>=1.24.3 @@ -19,3 +19,4 @@ dependencies: - xtensor-fftw>=0.2.5 - xtensor>=0.24.2 - parallelproj>=1.6.1 + - scipy diff --git a/python/parallelproj_sim.py b/python/parallelproj_sim.py index 318e695..82310fe 100644 --- a/python/parallelproj_sim.py +++ b/python/parallelproj_sim.py @@ -8,26 +8,25 @@ # device variable (cpu or cuda) that determines whether calculations # are performed on the cpu or cuda gpu -dev = 'cpu' +dev = "cpu" -#---------------------------------------------------------------------------- -#---------------------------------------------------------------------------- -#--- setup the scanner / LOR geometry --------------------------------------- -#---------------------------------------------------------------------------- -#---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- setup the scanner / LOR geometry --------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- # setup a line of response descriptor that describes the LOR start / endpoints of # a "narrow" clinical PET scanner with 9 rings -lor_descriptor = utils.DemoPETScannerLORDescriptor(np, - dev, - num_rings=4, - radial_trim=141) +lor_descriptor = utils.DemoPETScannerLORDescriptor( + np, dev, num_rings=2, radial_trim=141 +) -#---------------------------------------------------------------------------- -#---------------------------------------------------------------------------- -#--- setup a simple 3D test image ------------------------------------------- -#---------------------------------------------------------------------------- -#---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- setup a simple 3D test image ------------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- # image properties voxel_size = (2.66, 2.66, 2.66) @@ -40,16 +39,15 @@ # setup an image containing a box img = np.zeros(img_shape, dtype=np.float32, device=dev) -img[(n0 // 4):(3 * n0 // 4), (n1 // 4):(3 * n1 // 4), :] = 1 +img[(n0 // 4) : (3 * n0 // 4), (n1 // 4) : (3 * n1 // 4), :] = 1 -#---------------------------------------------------------------------------- -#---------------------------------------------------------------------------- -#--- setup a non-TOF projector and project ---------------------------------- -#---------------------------------------------------------------------------- -#---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- setup a non-TOF projector and project ---------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- -projector = utils.RegularPolygonPETProjector(lor_descriptor, img_shape, - voxel_size) +projector = utils.RegularPolygonPETProjector(lor_descriptor, img_shape, voxel_size) projector.tof = False # set this to True to get a time of flight projector @@ -68,30 +66,29 @@ end_index = lor_descriptor.scanner.num_modules * end_mods + end_inds -#---------------------------------------------------------------------------- -#---------------------------------------------------------------------------- -#---------------------------------------------------------------------------- -#---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- # show the scanner geometry and one view in one sinogram plane fig = plt.figure() -ax = fig.add_subplot(projection='3d') +ax = fig.add_subplot(projection="3d") lor_descriptor.scanner.show_lor_endpoints(ax) -lor_descriptor.show_views(ax, - views=np.asarray([lor_descriptor.num_views // 4], - device=dev), - planes=np.asarray( - [lor_descriptor.scanner.num_modules // 2], - device=dev)) +lor_descriptor.show_views( + ax, + views=np.asarray([lor_descriptor.num_views // 4], device=dev), + planes=np.asarray([lor_descriptor.scanner.num_modules // 2], device=dev), +) fig.tight_layout() fig.show() fig2, ax2 = plt.subplots(1, 3, figsize=(15, 5)) -ax2[0].imshow(np.asarray(to_device(img[:, :, 3], 'cpu'))) +ax2[0].imshow(np.asarray(to_device(img[:, :, 3], "cpu"))) if projector.tof: - ax2[1].imshow(np.asarray(to_device(img_fwd[:, :, 4, 15], 'cpu'))) + ax2[1].imshow(np.asarray(to_device(img_fwd[:, :, 0, 15], "cpu"))) else: - ax2[1].imshow(np.asarray(to_device(img_fwd[:, :, 4], 'cpu'))) -ax2[2].imshow(np.asarray(to_device(back_img[:, :, 3], 'cpu'))) + ax2[1].imshow(np.asarray(to_device(img_fwd[:, :, 0], "cpu"))) +ax2[2].imshow(np.asarray(to_device(back_img[:, :, 3], "cpu"))) fig2.tight_layout() -fig2.show() \ No newline at end of file +fig2.show() From b48ed6abedc328396450167b48c90e23c3b5fdb6 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Mon, 13 Nov 2023 22:45:53 -0800 Subject: [PATCH 04/27] add noise to the emission sinogram --- python/parallelproj_sim.py | 41 ++++++++++++++++++++++++++++---------- 1 file changed, 30 insertions(+), 11 deletions(-) diff --git a/python/parallelproj_sim.py b/python/parallelproj_sim.py index 82310fe..e775cd1 100644 --- a/python/parallelproj_sim.py +++ b/python/parallelproj_sim.py @@ -9,6 +9,8 @@ # are performed on the cpu or cuda gpu dev = "cpu" +expected_num_trues = 1e6 +np.random.seed(1) # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -30,7 +32,7 @@ # image properties voxel_size = (2.66, 2.66, 2.66) -num_trans = 160 +num_trans = 140 num_ax = 2 * lor_descriptor.scanner.num_modules # setup a box like test image @@ -48,22 +50,39 @@ # ---------------------------------------------------------------------------- projector = utils.RegularPolygonPETProjector(lor_descriptor, img_shape, voxel_size) - projector.tof = False # set this to True to get a time of flight projector -img_fwd = projector(img) -back_img = projector.adjoint(img_fwd) +# forward project the image +noise_free_sinogram = projector(img) + +# rescale the forward projection and image such that we get the expected number of trues +scale = expected_num_trues / np.sum(noise_free_sinogram) +noise_free_sinogram *= scale +img *= scale +# calculate the sensitivity image +sens_img = projector.adjoint(np.ones(noise_free_sinogram.shape, device=dev, dtype=np.float32)) -# get the start and end points of all sinogram bins -xstart, xend = lor_descriptor.get_lor_coordinates() +## get the start and end points of all sinogram bins +#xstart, xend = lor_descriptor.get_lor_coordinates() # get the two dimensional indices of all sinogram bins start_mods, end_mods, start_inds, end_inds = lor_descriptor.get_lor_indices() # convert the 2D indices to a 1D index -start_index = lor_descriptor.scanner.num_modules * start_mods + start_inds -end_index = lor_descriptor.scanner.num_modules * end_mods + end_inds +det_start_index = lor_descriptor.scanner.num_modules * start_mods + start_inds +det_end_index = lor_descriptor.scanner.num_modules * end_mods + end_inds + +# add poisson noise to the noise free sinogram +noisy_sinogram = np.random.poisson(noise_free_sinogram) + +# ravel the noisy sinogram and the detector start and end "index" sinograms +noisy_sinogram = noisy_sinogram.ravel() +det_start_index = det_start_index.ravel() +det_end_index = det_end_index.ravel() + +# generate listmode data from the noisy sinogram +inds = np.repeat(np.arange(noisy_sinogram.shape[0]), noisy_sinogram) # ---------------------------------------------------------------------------- @@ -86,9 +105,9 @@ fig2, ax2 = plt.subplots(1, 3, figsize=(15, 5)) ax2[0].imshow(np.asarray(to_device(img[:, :, 3], "cpu"))) if projector.tof: - ax2[1].imshow(np.asarray(to_device(img_fwd[:, :, 0, 15], "cpu"))) + ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0, 15], "cpu"))) else: - ax2[1].imshow(np.asarray(to_device(img_fwd[:, :, 0], "cpu"))) -ax2[2].imshow(np.asarray(to_device(back_img[:, :, 3], "cpu"))) + ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0], "cpu"))) +ax2[2].imshow(np.asarray(to_device(sens_img[:, :, 3], "cpu"))) fig2.tight_layout() fig2.show() From a4c72d189f2da8dca8c329d6a36dcdf4f454cd24 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Mon, 13 Nov 2023 23:22:55 -0800 Subject: [PATCH 05/27] add noise + event generation and plot of first 3 events --- python/parallelproj_sim.py | 69 ++++++++++++++++++++++++-------------- 1 file changed, 43 insertions(+), 26 deletions(-) diff --git a/python/parallelproj_sim.py b/python/parallelproj_sim.py index e775cd1..ea2e1b5 100644 --- a/python/parallelproj_sim.py +++ b/python/parallelproj_sim.py @@ -9,7 +9,7 @@ # are performed on the cpu or cuda gpu dev = "cpu" -expected_num_trues = 1e6 +expected_num_trues = 5e6 np.random.seed(1) # ---------------------------------------------------------------------------- @@ -63,27 +63,34 @@ # calculate the sensitivity image sens_img = projector.adjoint(np.ones(noise_free_sinogram.shape, device=dev, dtype=np.float32)) -## get the start and end points of all sinogram bins -#xstart, xend = lor_descriptor.get_lor_coordinates() +# get a lookup table that contains the world coordinates of all scanner detectors +scanner_lut = lor_descriptor.scanner.all_lor_endpoints # get the two dimensional indices of all sinogram bins start_mods, end_mods, start_inds, end_inds = lor_descriptor.get_lor_indices() -# convert the 2D indices to a 1D index -det_start_index = lor_descriptor.scanner.num_modules * start_mods + start_inds -det_end_index = lor_descriptor.scanner.num_modules * end_mods + end_inds +# generate two sinograms that contain the linearized detector start and end indices +sino_det_start_index = lor_descriptor.scanner.num_lor_endpoints_per_module[0] * start_mods + start_inds +sino_det_end_index = lor_descriptor.scanner.num_lor_endpoints_per_module[0] * end_mods + end_inds # add poisson noise to the noise free sinogram noisy_sinogram = np.random.poisson(noise_free_sinogram) # ravel the noisy sinogram and the detector start and end "index" sinograms noisy_sinogram = noisy_sinogram.ravel() -det_start_index = det_start_index.ravel() -det_end_index = det_end_index.ravel() +sino_det_start_index = sino_det_start_index.ravel() +sino_det_end_index = sino_det_end_index.ravel() + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- # generate listmode data from the noisy sinogram -inds = np.repeat(np.arange(noisy_sinogram.shape[0]), noisy_sinogram) +event_sino_inds = np.repeat(np.arange(noisy_sinogram.shape[0]), noisy_sinogram) +# shuffle the event sinogram indices +np.random.shuffle(event_sino_inds) +event_det_id_1 = sino_det_start_index[event_sino_inds] +event_det_id_2 = sino_det_end_index[event_sino_inds] # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -91,23 +98,33 @@ # ---------------------------------------------------------------------------- # show the scanner geometry and one view in one sinogram plane -fig = plt.figure() +fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(projection="3d") -lor_descriptor.scanner.show_lor_endpoints(ax) -lor_descriptor.show_views( - ax, - views=np.asarray([lor_descriptor.num_views // 4], device=dev), - planes=np.asarray([lor_descriptor.scanner.num_modules // 2], device=dev), -) +lor_descriptor.scanner.show_lor_endpoints(ax, show_linear_index=True, annotation_fontsize=4) + +# plot the LORs of the first n events +for i, col in enumerate(('red','green','blue')): + xstart = scanner_lut[event_det_id_1[i], :] + xend = scanner_lut[event_det_id_2[i], :] + + ax.plot( + [xstart[0], xend[0]], + [xstart[1], xend[1]], + [xstart[2], xend[2]], + color=col, + linewidth=1., + ) + fig.tight_layout() fig.show() - -fig2, ax2 = plt.subplots(1, 3, figsize=(15, 5)) -ax2[0].imshow(np.asarray(to_device(img[:, :, 3], "cpu"))) -if projector.tof: - ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0, 15], "cpu"))) -else: - ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0], "cpu"))) -ax2[2].imshow(np.asarray(to_device(sens_img[:, :, 3], "cpu"))) -fig2.tight_layout() -fig2.show() +# +#fig2, ax2 = plt.subplots(1, 3, figsize=(15, 5)) +#ax2[0].imshow(np.asarray(to_device(img[:, :, 3], "cpu"))) +#if projector.tof: +# ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0, 15], "cpu"))) +#else: +# ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0], "cpu"))) +#ax2[2].imshow(np.asarray(to_device(sens_img[:, :, 3], "cpu"))) +#fig2.tight_layout() +#fig2.show() +# \ No newline at end of file From ffc63f56c0d02d9892b7f9b1783cd5a3bb1c4a6a Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Tue, 14 Nov 2023 00:09:05 -0800 Subject: [PATCH 06/27] add LM recon (WIP) --- python/parallelproj_sim.py | 107 ++++++++++++++++++++++++++++--------- python/utils.py | 4 ++ 2 files changed, 87 insertions(+), 24 deletions(-) diff --git a/python/parallelproj_sim.py b/python/parallelproj_sim.py index ea2e1b5..149fe2d 100644 --- a/python/parallelproj_sim.py +++ b/python/parallelproj_sim.py @@ -1,15 +1,22 @@ +#TODO: - absolute scale of recon (maybe with sinogram recon) +# - additive MLEM + from __future__ import annotations +import parallelproj import utils import array_api_compat.numpy as np import matplotlib.pyplot as plt from array_api_compat import to_device +from scipy.ndimage import gaussian_filter # device variable (cpu or cuda) that determines whether calculations # are performed on the cpu or cuda gpu dev = "cpu" -expected_num_trues = 5e6 +expected_num_trues = 1e6 +num_iter = 2 +num_subsets = 20 np.random.seed(1) # ---------------------------------------------------------------------------- @@ -42,6 +49,7 @@ # setup an image containing a box img = np.zeros(img_shape, dtype=np.float32, device=dev) img[(n0 // 4) : (3 * n0 // 4), (n1 // 4) : (3 * n1 // 4), :] = 1 +img[(7*n0 // 16) : (9 * n0 // 16), (6*n1 // 16) : (8 * n1 // 16), :] = 2 # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -63,9 +71,6 @@ # calculate the sensitivity image sens_img = projector.adjoint(np.ones(noise_free_sinogram.shape, device=dev, dtype=np.float32)) -# get a lookup table that contains the world coordinates of all scanner detectors -scanner_lut = lor_descriptor.scanner.all_lor_endpoints - # get the two dimensional indices of all sinogram bins start_mods, end_mods, start_inds, end_inds = lor_descriptor.get_lor_indices() @@ -81,6 +86,9 @@ sino_det_start_index = sino_det_start_index.ravel() sino_det_end_index = sino_det_end_index.ravel() +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- convert the index sinograms in to listmode data ------------------------ # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -92,11 +100,61 @@ event_det_id_1 = sino_det_start_index[event_sino_inds] event_det_id_2 = sino_det_end_index[event_sino_inds] -# ---------------------------------------------------------------------------- -# ---------------------------------------------------------------------------- -# ---------------------------------------------------------------------------- -# ---------------------------------------------------------------------------- +print(f'number of events: {event_det_id_1.shape[0]}') + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +#---- convert LM detector ID arrays into PRD here --------------------------- +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- + +# get a lookup table that contains the world coordinates of all scanner detectors +# this is a 2D array of shape (num_detectors, 3) +scanner_lut = lor_descriptor.scanner.all_lor_endpoints +# +# +# +# +# +# + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +#---- read events back from PRD here ---------------------------------------- +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- + +# +# +# +# +# +# + +# hack until we have the reader / writer implemented +xstart = scanner_lut[event_det_id_1, :] +xend = scanner_lut[event_det_id_2, :] + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +#---- LM recon using the event detector IDs and the scanner LUT ------------- +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- + +recon = np.ones(img_shape, dtype=np.float32, device=dev) + +for it in range(num_iter): + for isub in range(num_subsets): + print(f'it {(it+1):03} / ss {(isub+1):03}', end='\r') + xs_sub = xstart[isub::num_subsets,:] + xe_sub = xend[isub::num_subsets,:] + exp = parallelproj.joseph3d_fwd(xs_sub, xe_sub, recon, projector.img_origin, voxel_size) + tmp = parallelproj.joseph3d_back(xs_sub, xe_sub, img_shape, projector.img_origin, voxel_size, 1/exp) + recon *= (tmp / (sens_img / num_subsets)) + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- # show the scanner geometry and one view in one sinogram plane fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(projection="3d") @@ -104,27 +162,28 @@ # plot the LORs of the first n events for i, col in enumerate(('red','green','blue')): - xstart = scanner_lut[event_det_id_1[i], :] - xend = scanner_lut[event_det_id_2[i], :] + xs = scanner_lut[event_det_id_1[i], :] + xe= scanner_lut[event_det_id_2[i], :] ax.plot( - [xstart[0], xend[0]], - [xstart[1], xend[1]], - [xstart[2], xend[2]], + [xs[0], xe[0]], + [xs[1], xe[1]], + [xs[2], xe[2]], color=col, linewidth=1., ) fig.tight_layout() fig.show() -# -#fig2, ax2 = plt.subplots(1, 3, figsize=(15, 5)) -#ax2[0].imshow(np.asarray(to_device(img[:, :, 3], "cpu"))) -#if projector.tof: -# ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0, 15], "cpu"))) -#else: -# ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0], "cpu"))) -#ax2[2].imshow(np.asarray(to_device(sens_img[:, :, 3], "cpu"))) -#fig2.tight_layout() -#fig2.show() -# \ No newline at end of file + +vmax = 1.2*img.max() +fig2, ax2 = plt.subplots(1, 4, figsize=(16, 4)) +ax2[0].imshow(np.asarray(to_device(img[:, :, 1], "cpu")), vmin = 0, vmax = vmax, cmap = 'Greys') +if projector.tof: + ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0, 15], "cpu")), cmap = 'Greys') +else: + ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0], "cpu")), cmap = 'Greys') +ax2[2].imshow(np.asarray(to_device(recon[:, :, 1], "cpu")), vmin = 0, vmax = vmax, cmap = 'Greys') +ax2[3].imshow(gaussian_filter(np.asarray(to_device(recon[:, :, 1], "cpu")), 1.5), vmin = 0, vmax = vmax, cmap = 'Greys') +fig2.tight_layout() +fig2.show() diff --git a/python/utils.py b/python/utils.py index 9470e75..9bc3a57 100644 --- a/python/utils.py +++ b/python/utils.py @@ -1118,6 +1118,10 @@ def tof_parameters(self, value: TOFParameters) -> None: raise ValueError('tof_parameters must be a TOFParameters object') self._tof_parameters = value + @property + def img_origin(self) -> npt.NDArray: + return self._img_origin + def _apply(self, x): """nonTOF forward projection of input image x including image based resolution model""" From f4ac538db2373cdb6e8dbed143dec879b780675c Mon Sep 17 00:00:00 2001 From: Fedor Date: Tue, 14 Nov 2023 17:02:20 +0000 Subject: [PATCH 07/27] added timestamps for the list of events --- python/parallelproj_sim.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/python/parallelproj_sim.py b/python/parallelproj_sim.py index 149fe2d..43fb7c8 100644 --- a/python/parallelproj_sim.py +++ b/python/parallelproj_sim.py @@ -81,7 +81,7 @@ # add poisson noise to the noise free sinogram noisy_sinogram = np.random.poisson(noise_free_sinogram) -# ravel the noisy sinogram and the detector start and end "index" sinograms +# ravel the noisy sinogram and the detector start and end "index" sinograms noisy_sinogram = noisy_sinogram.ravel() sino_det_start_index = sino_det_start_index.ravel() sino_det_end_index = sino_det_end_index.ravel() @@ -94,9 +94,27 @@ # generate listmode data from the noisy sinogram event_sino_inds = np.repeat(np.arange(noisy_sinogram.shape[0]), noisy_sinogram) + +# generate timestamps +acquisition_length = 20 # say, in mins +timestamps_in_lors = np.array([np.sort(np.random.uniform(0, acquisition_length, + size=noisy_sinogram[l])) for l in range(len(noisy_sinogram))]) + # shuffle the event sinogram indices np.random.shuffle(event_sino_inds) +# assign timestamps for events - need one forward run over event_sino_inds +timestamps_iter_table = np.zeros_like(noisy_sinogram, dtype=np.int64) # possibly many counts +timestamps_in_events = np.zeros_like(noisy_sinogram, dtype=np.float32) +for ind in event_sino_inds: # sorry, this is slow and ugly + timestamps_in_events[ind] = timestamps_in_lors[ind, timestamps_iter_table[ind]] + timestamps_iter_table[ind] += 1 + +# at this stage - lors are shuffled but timestamps are out of sequential order +# need to sort globally event_sino_inds according to timestamps +evend_sino_inds = event_sino_inds[np.argsort(timestamps_in_events)] + + event_det_id_1 = sino_det_start_index[event_sino_inds] event_det_id_2 = sino_det_end_index[event_sino_inds] From 343bda80c2c1b0f7b6c3cc8193da4b75565785fb Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Tue, 14 Nov 2023 09:17:21 -0800 Subject: [PATCH 08/27] correct use of resolution model --- python/parallelproj_sim.py | 45 +++++++++++++++++++++++--------------- python/utils.py | 17 ++++++++------ 2 files changed, 37 insertions(+), 25 deletions(-) diff --git a/python/parallelproj_sim.py b/python/parallelproj_sim.py index 149fe2d..d75308f 100644 --- a/python/parallelproj_sim.py +++ b/python/parallelproj_sim.py @@ -1,5 +1,4 @@ -#TODO: - absolute scale of recon (maybe with sinogram recon) -# - additive MLEM +#TODO: - additive MLEM from __future__ import annotations @@ -15,7 +14,7 @@ dev = "cpu" expected_num_trues = 1e6 -num_iter = 2 +num_iter = 3 num_subsets = 20 np.random.seed(1) @@ -28,7 +27,7 @@ # setup a line of response descriptor that describes the LOR start / endpoints of # a "narrow" clinical PET scanner with 9 rings lor_descriptor = utils.DemoPETScannerLORDescriptor( - np, dev, num_rings=2, radial_trim=141 + np, dev, num_rings=4, radial_trim=141 ) # ---------------------------------------------------------------------------- @@ -48,8 +47,8 @@ # setup an image containing a box img = np.zeros(img_shape, dtype=np.float32, device=dev) -img[(n0 // 4) : (3 * n0 // 4), (n1 // 4) : (3 * n1 // 4), :] = 1 -img[(7*n0 // 16) : (9 * n0 // 16), (6*n1 // 16) : (8 * n1 // 16), :] = 2 +img[(n0 // 4) : (3 * n0 // 4), (n1 // 4) : (3 * n1 // 4), 2:-2] = 1 +img[(7*n0 // 16) : (9 * n0 // 16), (6*n1 // 16) : (8 * n1 // 16), 2:-2] = 2. # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -57,7 +56,9 @@ # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- -projector = utils.RegularPolygonPETProjector(lor_descriptor, img_shape, voxel_size) +# setup a simple image-based resolution model using 4.5mm FWHM Gaussian smoothing +res_model = parallelproj.GaussianFilterOperator(img_shape, sigma=4.5 / (2.355 * np.asarray(voxel_size))) +projector = utils.RegularPolygonPETProjector(lor_descriptor, img_shape, voxel_size, resolution_model=res_model) projector.tof = False # set this to True to get a time of flight projector # forward project the image @@ -149,9 +150,15 @@ print(f'it {(it+1):03} / ss {(isub+1):03}', end='\r') xs_sub = xstart[isub::num_subsets,:] xe_sub = xend[isub::num_subsets,:] - exp = parallelproj.joseph3d_fwd(xs_sub, xe_sub, recon, projector.img_origin, voxel_size) - tmp = parallelproj.joseph3d_back(xs_sub, xe_sub, img_shape, projector.img_origin, voxel_size, 1/exp) - recon *= (tmp / (sens_img / num_subsets)) + + recon_sm = res_model(recon) + + exp = parallelproj.joseph3d_fwd(xs_sub, xe_sub, recon_sm, projector.img_origin, voxel_size) + ratio_back = parallelproj.joseph3d_back(xs_sub, xe_sub, img_shape, projector.img_origin, voxel_size, 1/exp) + + ratio_back_sm = res_model.adjoint(ratio_back) + + recon *= (ratio_back_sm / (sens_img / num_subsets)) #---------------------------------------------------------------------------- #---------------------------------------------------------------------------- @@ -177,13 +184,15 @@ fig.show() vmax = 1.2*img.max() -fig2, ax2 = plt.subplots(1, 4, figsize=(16, 4)) -ax2[0].imshow(np.asarray(to_device(img[:, :, 1], "cpu")), vmin = 0, vmax = vmax, cmap = 'Greys') -if projector.tof: - ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0, 15], "cpu")), cmap = 'Greys') -else: - ax2[1].imshow(np.asarray(to_device(noise_free_sinogram[:, :, 0], "cpu")), cmap = 'Greys') -ax2[2].imshow(np.asarray(to_device(recon[:, :, 1], "cpu")), vmin = 0, vmax = vmax, cmap = 'Greys') -ax2[3].imshow(gaussian_filter(np.asarray(to_device(recon[:, :, 1], "cpu")), 1.5), vmin = 0, vmax = vmax, cmap = 'Greys') +fig2, ax2 = plt.subplots(3, recon.shape[2], figsize=(recon.shape[2] * 2, 3 * 2)) +for i in range(recon.shape[2]): + ax2[0,i].imshow(np.asarray(to_device(img[:, :, i], "cpu")), vmin = 0, vmax = vmax, cmap = 'Greys') + ax2[1,i].imshow(np.asarray(to_device(recon[:, :, i], "cpu")), vmin = 0, vmax = vmax, cmap = 'Greys') + ax2[2,i].imshow(gaussian_filter(np.asarray(to_device(recon[:, :, i], "cpu")), 1.5), vmin = 0, vmax = vmax, cmap = 'Greys') + + ax2[0,i].set_title(f'ground truth sl {i+1}', fontsize = 'small') + ax2[1,i].set_title(f'LM recon {i+1}', fontsize = 'small') + ax2[2,i].set_title(f'LM recon smoothed {i+1}', fontsize = 'small') + fig2.tight_layout() fig2.show() diff --git a/python/utils.py b/python/utils.py index 9bc3a57..1953103 100644 --- a/python/utils.py +++ b/python/utils.py @@ -1066,11 +1066,7 @@ def __init__(self, else: self._views = views - if resolution_model is None: - self._resolution_model = parallelproj.GaussianFilterOperator( - self.in_shape, sigma=4.5 / (2.355 * self._voxel_size)) - else: - self._resolution_model = resolution_model + self._resolution_model = resolution_model self._xstart, self._xend = lor_descriptor.get_lor_coordinates( views=self._views, sinogram_order=SinogramSpatialAxisOrder['RVP']) @@ -1126,7 +1122,11 @@ def _apply(self, x): """nonTOF forward projection of input image x including image based resolution model""" dev = device(x) - x_sm = self._resolution_model(x) + + if self._resolution_model is not None: + x_sm = self._resolution_model(x) + else: + x_sm = x if not self.tof: x_fwd = parallelproj.joseph3d_fwd(self._xstart, self._xend, x_sm, @@ -1167,7 +1167,10 @@ def _adjoint(self, y): device=dev), self.tof_parameters.num_sigmas, self.tof_parameters.num_tofbins) - return self._resolution_model.adjoint(y_back) + if self._resolution_model is not None: + y_back = self._resolution_model.adjoint(y_back) + + return y_back From c2ae0b933cd6846ee3d9eb740f5998b65f99649c Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Tue, 14 Nov 2023 09:33:33 -0800 Subject: [PATCH 09/27] merge fedor's changes (commented our by now) --- python/parallelproj_sim.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/python/parallelproj_sim.py b/python/parallelproj_sim.py index f31fdfa..e04777a 100644 --- a/python/parallelproj_sim.py +++ b/python/parallelproj_sim.py @@ -96,24 +96,24 @@ # generate listmode data from the noisy sinogram event_sino_inds = np.repeat(np.arange(noisy_sinogram.shape[0]), noisy_sinogram) -# generate timestamps -acquisition_length = 20 # say, in mins -timestamps_in_lors = np.array([np.sort(np.random.uniform(0, acquisition_length, - size=noisy_sinogram[l])) for l in range(len(noisy_sinogram))]) +## generate timestamps +#acquisition_length = 20 # say, in mins +#timestamps_in_lors = np.array([np.sort(np.random.uniform(0, acquisition_length, +# size=noisy_sinogram[l])) for l in range(len(noisy_sinogram))]) # shuffle the event sinogram indices np.random.shuffle(event_sino_inds) -# assign timestamps for events - need one forward run over event_sino_inds -timestamps_iter_table = np.zeros_like(noisy_sinogram, dtype=np.int64) # possibly many counts -timestamps_in_events = np.zeros_like(noisy_sinogram, dtype=np.float32) -for ind in event_sino_inds: # sorry, this is slow and ugly - timestamps_in_events[ind] = timestamps_in_lors[ind, timestamps_iter_table[ind]] - timestamps_iter_table[ind] += 1 +## assign timestamps for events - need one forward run over event_sino_inds +#timestamps_iter_table = np.zeros_like(noisy_sinogram, dtype=np.int64) # possibly many counts +#timestamps_in_events = np.zeros_like(noisy_sinogram, dtype=np.float32) +#for ind in event_sino_inds: # sorry, this is slow and ugly +# timestamps_in_events[ind] = timestamps_in_lors[ind, timestamps_iter_table[ind]] +# timestamps_iter_table[ind] += 1 -# at this stage - lors are shuffled but timestamps are out of sequential order -# need to sort globally event_sino_inds according to timestamps -evend_sino_inds = event_sino_inds[np.argsort(timestamps_in_events)] +## at this stage - lors are shuffled but timestamps are out of sequential order +## need to sort globally event_sino_inds according to timestamps +#evend_sino_inds = event_sino_inds[np.argsort(timestamps_in_events)] event_det_id_1 = sino_det_start_index[event_sino_inds] From 5cb73d12191444f3a6367bd56692f2bb976f4907 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Tue, 14 Nov 2023 13:57:41 -0800 Subject: [PATCH 10/27] add TOF simulation example --- python/parallelproj_sim_tof.py | 232 +++++++++++++++++++++++++++++++++ 1 file changed, 232 insertions(+) create mode 100644 python/parallelproj_sim_tof.py diff --git a/python/parallelproj_sim_tof.py b/python/parallelproj_sim_tof.py new file mode 100644 index 0000000..9b8d891 --- /dev/null +++ b/python/parallelproj_sim_tof.py @@ -0,0 +1,232 @@ +#TODO: - additive MLEM + +from __future__ import annotations + +import parallelproj +import utils +import array_api_compat.numpy as np +import matplotlib.pyplot as plt +from array_api_compat import to_device +from scipy.ndimage import gaussian_filter + +# device variable (cpu or cuda) that determines whether calculations +# are performed on the cpu or cuda gpu + +dev = "cpu" +expected_num_trues = 1e6 +num_iter = 3 +num_subsets = 20 +np.random.seed(1) + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- setup the scanner / LOR geometry --------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- + +# setup a line of response descriptor that describes the LOR start / endpoints of +# a "narrow" clinical PET scanner with 9 rings +lor_descriptor = utils.DemoPETScannerLORDescriptor( + np, dev, num_rings=4, radial_trim=141 +) + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- setup a simple 3D test image ------------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- + +# image properties +voxel_size = (2.66, 2.66, 2.66) +num_trans = 140 +num_ax = 2 * lor_descriptor.scanner.num_modules + +# setup a box like test image +img_shape = (num_trans, num_trans, num_ax) +n0, n1, n2 = img_shape + +# setup an image containing a box +img = np.zeros(img_shape, dtype=np.float32, device=dev) +img[(n0 // 4) : (3 * n0 // 4), (n1 // 4) : (3 * n1 // 4), 2:-2] = 1 +img[(7*n0 // 16) : (9 * n0 // 16), (6*n1 // 16) : (8 * n1 // 16), 2:-2] = 2. + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- setup a non-TOF projector and project ---------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- + +# setup a simple image-based resolution model using 4.5mm FWHM Gaussian smoothing +res_model = parallelproj.GaussianFilterOperator(img_shape, sigma=4.5 / (2.355 * np.asarray(voxel_size))) +projector = utils.RegularPolygonPETProjector(lor_descriptor, img_shape, voxel_size, resolution_model=res_model) +projector.tof = True # set this to True to get a time of flight projector + +# forward project the image +noise_free_sinogram = projector(img) + +# rescale the forward projection and image such that we get the expected number of trues +scale = expected_num_trues / np.sum(noise_free_sinogram) +noise_free_sinogram *= scale +img *= scale +# +# calculate the sensitivity image +sens_img = projector.adjoint(np.ones(noise_free_sinogram.shape, device=dev, dtype=np.float32)) + +# get the two dimensional indices of all sinogram bins +start_mods, end_mods, start_inds, end_inds = lor_descriptor.get_lor_indices() + +# generate two sinograms that contain the linearized detector start and end indices +sino_det_start_index = lor_descriptor.scanner.num_lor_endpoints_per_module[0] * start_mods + start_inds +sino_det_end_index = lor_descriptor.scanner.num_lor_endpoints_per_module[0] * end_mods + end_inds + +# add poisson noise to the noise free sinogram +noisy_sinogram = np.random.poisson(noise_free_sinogram) + +## ravel the noisy sinogram and the detector start and end "index" sinograms +noisy_sinogram = noisy_sinogram.ravel() + +# repeat number of TOF bin times here +num_tof_bins = projector.tof_parameters.num_tofbins +sino_det_start_index = np.repeat(sino_det_start_index.ravel(), num_tof_bins) +sino_det_end_index = np.repeat(sino_det_end_index.ravel(), num_tof_bins) + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- convert the index sinograms in to listmode data ------------------------ +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- + +# generate listmode data from the noisy sinogram +event_sino_inds = np.repeat(np.arange(noisy_sinogram.shape[0]), noisy_sinogram) + +## generate timestamps +#acquisition_length = 20 # say, in mins +#timestamps_in_lors = np.array([np.sort(np.random.uniform(0, acquisition_length, +# size=noisy_sinogram[l])) for l in range(len(noisy_sinogram))]) + +# shuffle the event sinogram indices +np.random.shuffle(event_sino_inds) + +## assign timestamps for events - need one forward run over event_sino_inds +#timestamps_iter_table = np.zeros_like(noisy_sinogram, dtype=np.int64) # possibly many counts +#timestamps_in_events = np.zeros_like(noisy_sinogram, dtype=np.float32) +#for ind in event_sino_inds: # sorry, this is slow and ugly +# timestamps_in_events[ind] = timestamps_in_lors[ind, timestamps_iter_table[ind]] +# timestamps_iter_table[ind] += 1 + +## at this stage - lors are shuffled but timestamps are out of sequential order +## need to sort globally event_sino_inds according to timestamps +#evend_sino_inds = event_sino_inds[np.argsort(timestamps_in_events)] + +event_det_id_1 = sino_det_start_index[event_sino_inds] +event_det_id_2 = sino_det_end_index[event_sino_inds] +event_tof_bin = event_sino_inds % num_tof_bins + +print(f'number of events: {event_det_id_1.shape[0]}') + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +#---- convert LM detector ID arrays into PRD here --------------------------- +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- + +# get a lookup table that contains the world coordinates of all scanner detectors +# this is a 2D array of shape (num_detectors, 3) +scanner_lut = lor_descriptor.scanner.all_lor_endpoints + +# +# +# +# +# +# + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +#---- read events back from PRD here ---------------------------------------- +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- + +# +# +# +# +# +# + +# hack until we have the reader / writer implemented +xstart = scanner_lut[event_det_id_1, :] +xend = scanner_lut[event_det_id_2, :] +event_tof_bin = event_tof_bin +# tof bin width in mm +tofbin_width = projector.tof_parameters.tofbin_width +# sigma of the Gaussian TOF kernel in mm +sigma_tof = np.array([projector.tof_parameters.sigma_tof], dtype = np.float32) +tofcenter_offset = np.array([projector.tof_parameters.tofcenter_offset], dtype = np.float32) +nsigmas = projector.tof_parameters.num_sigmas + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +#---- LM recon using the event detector IDs and the scanner LUT ------------- +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- + +recon = np.ones(img_shape, dtype=np.float32, device=dev) + +for it in range(num_iter): + for isub in range(num_subsets): + print(f'it {(it+1):03} / ss {(isub+1):03}', end='\r') + xs_sub = xstart[isub::num_subsets,:] + xe_sub = xend[isub::num_subsets,:] + # in parallelproj the "0" TOF bin is the central TOF bin + # in PETSIRD the TOFbin number is non-negative + event_tof_bin_sub = event_tof_bin[isub::num_subsets] - num_tof_bins // 2 + + recon_sm = res_model(recon) + + exp = parallelproj.joseph3d_fwd_tof_lm(xs_sub, xe_sub, recon, projector.img_origin, voxel_size, tofbin_width, + sigma_tof, tofcenter_offset, nsigmas, event_tof_bin_sub) + + ratio_back = parallelproj.joseph3d_back_tof_lm(xs_sub, xe_sub, img_shape, projector.img_origin, + voxel_size, 1/exp, tofbin_width, sigma_tof, tofcenter_offset, nsigmas, event_tof_bin_sub) + + ratio_back_sm = res_model.adjoint(ratio_back) + + recon *= (ratio_back_sm / (sens_img / num_subsets)) + +#---------------------------------------------------------------------------- +#---------------------------------------------------------------------------- +# show the scanner geometry and one view in one sinogram plane +fig = plt.figure(figsize=(8, 8)) +ax = fig.add_subplot(projection="3d") +lor_descriptor.scanner.show_lor_endpoints(ax, show_linear_index=True, annotation_fontsize=4) + +# plot the LORs of the first n events +for i, col in enumerate(('red','green','blue')): + xs = scanner_lut[event_det_id_1[i], :] + xe= scanner_lut[event_det_id_2[i], :] + + ax.plot( + [xs[0], xe[0]], + [xs[1], xe[1]], + [xs[2], xe[2]], + color=col, + linewidth=1., + ) + +fig.tight_layout() +fig.show() + +vmax = 1.2*img.max() +fig2, ax2 = plt.subplots(3, recon.shape[2], figsize=(recon.shape[2] * 2, 3 * 2)) +for i in range(recon.shape[2]): + ax2[0,i].imshow(np.asarray(to_device(img[:, :, i], "cpu")), vmin = 0, vmax = vmax, cmap = 'Greys') + ax2[1,i].imshow(np.asarray(to_device(recon[:, :, i], "cpu")), vmin = 0, vmax = vmax, cmap = 'Greys') + ax2[2,i].imshow(gaussian_filter(np.asarray(to_device(recon[:, :, i], "cpu")), 1.5), vmin = 0, vmax = vmax, cmap = 'Greys') + + ax2[0,i].set_title(f'ground truth sl {i+1}', fontsize = 'small') + ax2[1,i].set_title(f'LM recon {i+1}', fontsize = 'small') + ax2[2,i].set_title(f'LM recon smoothed {i+1}', fontsize = 'small') + +fig2.tight_layout() +fig2.show() From 9c7453c644759f1feada7151bd71fe9267aaad2c Mon Sep 17 00:00:00 2001 From: Imraj-Singh <72553490+Imraj-Singh@users.noreply.github.com> Date: Tue, 14 Nov 2023 22:17:25 +0000 Subject: [PATCH 11/27] read write with yardl --- python/demo_read_reconstruct.py | 52 ++++++++++++++ python/demo_write.py | 121 ++++++++++++++++++++++++++++++++ python/io_yardl.py | 89 +++++++++++++++++++++++ 3 files changed, 262 insertions(+) create mode 100644 python/demo_read_reconstruct.py create mode 100644 python/demo_write.py create mode 100644 python/io_yardl.py diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py new file mode 100644 index 0000000..951e0a4 --- /dev/null +++ b/python/demo_read_reconstruct.py @@ -0,0 +1,52 @@ +# TODO: - absolute scale of recon (maybe with sinogram recon) +# - additive MLEM + +from __future__ import annotations + +import parallelproj +import array_api_compat.numpy as np +import matplotlib.pyplot as plt +from io_yardl import read_yardl + + +dev = "cpu" + +# image properties +voxel_size = (2.66, 2.66, 2.66) +img_shape = (128, 128, 4) +img_origin = [-168.91, -168.91, -3.9900002] + +event_det_id_1, event_det_id_2, scanner_lut = read_yardl("write_test.yardl") + +# hack until we have the reader / writer implemented + +xstart = scanner_lut[event_det_id_1, :] +xend = scanner_lut[event_det_id_2, :] + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# ---- LM recon using the event detector IDs and the scanner LUT ------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- + +recon = np.ones(img_shape, dtype=np.float32, device=dev) +img = np.tile(np.load("../SL.npy")[..., None], (1, 1, 4)) +num_iter = 2 +num_subsets = 20 +for it in range(num_iter): + for isub in range(num_subsets): + print(f"it {(it+1):03} / ss {(isub+1):03}", end="\r") + xs_sub = xstart[isub::num_subsets, :] + xe_sub = xend[isub::num_subsets, :] + exp = parallelproj.joseph3d_fwd(xs_sub, xe_sub, recon, img_origin, voxel_size) + tmp = parallelproj.joseph3d_back( + xs_sub, xe_sub, img_shape, img_origin, voxel_size, 1 / exp + ) + recon *= tmp + +fig, ax = plt.subplots(1, 3, figsize=(12, 4)) +fig.colorbar(ax[0].imshow(img[:, :, 2])) +fig.colorbar(ax[1].imshow(recon[:, :, 2])) +# difference +fig.colorbar(ax[2].imshow(recon[:, :, 2] - img[:, :, 2])) +fig.savefig("read_reconstruct_test.png") diff --git a/python/demo_write.py b/python/demo_write.py new file mode 100644 index 0000000..5943486 --- /dev/null +++ b/python/demo_write.py @@ -0,0 +1,121 @@ +# TODO: - absolute scale of recon (maybe with sinogram recon) +# - additive MLEM + +from __future__ import annotations + +import parallelproj +import utils +import array_api_compat.numpy as np +import matplotlib.pyplot as plt +from array_api_compat import to_device +from scipy.ndimage import gaussian_filter +from io_yardl import write_yardl, read_yardl + +# device variable (cpu or cuda) that determines whether calculations +# are performed on the cpu or cuda gpu + +dev = "cpu" +expected_num_trues = 1e6 +num_iter = 2 +num_subsets = 20 +np.random.seed(1) + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- setup the scanner / LOR geometry --------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- + +# setup a line of response descriptor that describes the LOR start / endpoints of +# a "narrow" clinical PET scanner with 9 rings +lor_descriptor = utils.DemoPETScannerLORDescriptor( + np, dev, num_rings=2, radial_trim=141 +) + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- setup a simple 3D test image ------------------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- + +# image properties +voxel_size = (2.66, 2.66, 2.66) +num_trans = 128 +num_ax = 2 * lor_descriptor.scanner.num_modules + +# setup a box like test image +img_shape = (num_trans, num_trans, num_ax) +n0, n1, n2 = img_shape + +# setup an image containing a box +img = np.tile(np.load("../SL.npy")[..., None], (1, 1, 4)) + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- setup a non-TOF projector and project ---------------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- + +projector = utils.RegularPolygonPETProjector(lor_descriptor, img_shape, voxel_size) +projector.tof = False # set this to True to get a time of flight projector + +# forward project the image +noise_free_sinogram = projector(img) + +# rescale the forward projection and image such that we get the expected number of trues +scale = expected_num_trues / np.sum(noise_free_sinogram) +noise_free_sinogram *= scale +img *= scale + +# calculate the sensitivity image +sens_img = projector.adjoint( + np.ones(noise_free_sinogram.shape, device=dev, dtype=np.float32) +) + +# get the two dimensional indices of all sinogram bins +start_mods, end_mods, start_inds, end_inds = lor_descriptor.get_lor_indices() + +# generate two sinograms that contain the linearized detector start and end indices +sino_det_start_index = ( + lor_descriptor.scanner.num_lor_endpoints_per_module[0] * start_mods + start_inds +) +sino_det_end_index = ( + lor_descriptor.scanner.num_lor_endpoints_per_module[0] * end_mods + end_inds +) + +# add poisson noise to the noise free sinogram +noisy_sinogram = np.random.poisson(noise_free_sinogram) + +# ravel the noisy sinogram and the detector start and end "index" sinograms +noisy_sinogram = noisy_sinogram.ravel() +sino_det_start_index = sino_det_start_index.ravel() +sino_det_end_index = sino_det_end_index.ravel() + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# --- convert the index sinograms in to listmode data ------------------------ +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- + +# generate listmode data from the noisy sinogram +event_sino_inds = np.repeat(np.arange(noisy_sinogram.shape[0]), noisy_sinogram) +# shuffle the event sinogram indices +np.random.shuffle(event_sino_inds) + +event_det_id_1 = sino_det_start_index[event_sino_inds] +event_det_id_2 = sino_det_end_index[event_sino_inds] + +print(f"number of events: {event_det_id_1.shape[0]}") + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- +# ---- convert LM detector ID arrays into PRD here --------------------------- +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- + +# get a lookup table that contains the world coordinates of all scanner detectors +# this is a 2D array of shape (num_detectors, 3) +scanner_lut = lor_descriptor.scanner.all_lor_endpoints + +# write and re-read +write_yardl(event_det_id_1, event_det_id_2, scanner_lut, output_file="write_test.yardl") diff --git a/python/io_yardl.py b/python/io_yardl.py new file mode 100644 index 0000000..6ba1532 --- /dev/null +++ b/python/io_yardl.py @@ -0,0 +1,89 @@ +from __future__ import annotations +import sys + +sys.path.append("/workspaces/LMSimRecon/PETSIRD/python") +import prd +import numpy as np + + +def write_yardl( + det_1_array: np.array[int], + det_2_array: np.array[int], + scanner_lut: np.array[float], + output_file: str | None = None, +) -> None: + """Writes a yardl list mode file from two detector arrays and a scanner lookup table. + + Args: + output_file (str): File to write to. + det_1_array (np.array[int]): Indices corresponding to detector 1. + det_2_array (np.array[int]): Indices corresponding to detector 2. + scanner_lut (np.array[float]): Lookup table for detector world coordinates. + """ + num_counts = len(det_1_array) + # create list of detectors + detectors = [] + for i in range(scanner_lut.shape[0]): + detectors.append( + prd.Detector( + id=i, x=scanner_lut[i, 0], y=scanner_lut[i, 1], z=scanner_lut[i, 2] + ) + ) + scanner_information = prd.ScannerInformation(detectors=detectors) + events = [] + for i in range(num_counts): + events.append( + prd.CoincidenceEvent( + detector_1_id=int(det_1_array[i]), detector_2_id=int(det_2_array[i]) + ) + ) + time_block = prd.TimeBlock(id=0, prompt_events=events) + if output_file is None: + with prd.BinaryPrdExperimentWriter(sys.stdout.buffer) as writer: + writer.write_header(prd.Header(scanner=scanner_information)) + writer.write_time_blocks((time_block,)) + else: + with prd.BinaryPrdExperimentWriter(output_file) as writer: + writer.write_header(prd.Header(scanner=scanner_information)) + writer.write_time_blocks((time_block,)) + + +def read_yardl(prd_file: str) -> tuple[np.array[int], np.array[int], np.array[float]]: + """Reads a yardl list mode file and returns two detector arrays and a scanner lookup table. + + Args: + prd_file (str): File to read from. + + Returns: + tuple[np.array[int], np.array[int], np.array[float]]: Two detector arrays and a scanner lookup table. + """ + with prd.BinaryPrdExperimentReader(prd_file) as reader: + # Read header and build lookup table + header = reader.read_header() + scanner_information = header.scanner + detectors = scanner_information.detectors + scanner_lut = np.zeros((len(detectors), 3)) + for i, detector in enumerate(detectors): + scanner_lut[i, 0] = detector.x + scanner_lut[i, 1] = detector.y + scanner_lut[i, 2] = detector.z + # Read events + detector_hits = [] + for time_blocks in reader.read_time_blocks(): + for event in time_blocks.prompt_events: + prompt = [event.detector_1_id, event.detector_2_id] + detector_hits.append(prompt) + det_1, det_2 = np.asarray(detector_hits).T + return det_1, det_2, scanner_lut + + +if __name__ == "__main__": + det_1 = np.asarray([0, 1, 2, 3, 4]) + det_2 = np.asarray([5, 6, 7, 8, 9]) + scanner_lut = np.random.rand(10, 3) + write_yardl(det_1, det_2, scanner_lut, output_file="test.yardl") + det_1_read, det_2_read, scanner_lut_read = read_yardl("test.yardl") + print(scanner_lut == scanner_lut_read) + print(np.isclose(scanner_lut, scanner_lut_read).all()) + print(scanner_lut.dtype) + print(scanner_lut_read.dtype) From a65f396a8973fa164e37475e63ec13fbad4b6613 Mon Sep 17 00:00:00 2001 From: Imraj-Singh <72553490+Imraj-Singh@users.noreply.github.com> Date: Tue, 14 Nov 2023 22:33:48 +0000 Subject: [PATCH 12/27] being a good dev ;) --- data/SL.npy | Bin 0 -> 65664 bytes python/demo_read_reconstruct.py | 2 +- python/demo_write.py | 4 ++-- python/io_yardl.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) create mode 100644 data/SL.npy diff --git a/data/SL.npy b/data/SL.npy new file mode 100644 index 0000000000000000000000000000000000000000..85b7aa89f9c34ae25b7dbd26ca81ea63b5a82a8a GIT binary patch literal 65664 zcmeI52UrwW*T)Sic7ru&Y$*0h-!gM9qNsGlf>^O)!3LrzWkqaYD1ri3?7eI3U82}~ zSJYrjEEuuHp!x1!-+f;`L|tGP7RBcocbIa{`Tx$nckZ-1(T<%QoZWsXG`f&q{eJzu z`+C)v)v7NHFsxh^zk7ZbgZLlr&!M@)Zj&EL$V`3eXEva24$!UjW-NqFa_JvuC6?(rF z-JXo2BRvbTyE^6B-J<2$YNyJ~z_R0q3G-Nnb=U^mVjt`a8?XhNu>H*g%%=kP{Fr8L zQFwvu%8n3Wne@6?HhL@Fo^8OkFYmz)>=?pM>G`pp%L3VqnF_X>PbcOmm~-mK3G-Nn zb=U^mVjt`a8?XhNuniyZl}|B!b77kqzPIq}5+oPRyhPg~3#rd@6K1vBlSPjXV0On8 zZ2Uh8X7;awz2{S*LgQ7Mur2n%zOVsXunF7n0blS5-`^b2d?4s*yUQjQ?7StLu7$}r z4_PnP4}MAkC)`wgD>tAz7YdrK}`7oUA_)+>g$rm@yq8~cwDR^KQ>~oSf=4q zDm$w;vo#N7HQy@OzyHaXCI8ibpkQZ*@;sOKXZC(W*_sUl*tmxRb6;%4hQ4gdjK|%j zZb#35m@toJSch$}E%w2_umM}JiLwA6@CBdnjTjIMVnS?)F<(R0(%^j$7BZ{crjqQ` z#RlR@^Q9E;=E`>K1hQFwf6V(;o|njb3f5s%AlrX%6mv*(VI?aYvR(U6ljFJw>bKj6 zDje%VUEVp<_#@6Va-;+4zqS7`VIIq{4%=W`?1Oz_1GZojw&4T5P!8Z5F(4MigxC-x zV%1WN`7*5F=0z;-FE){vUtNWoIW}TOu4CC)-d3ma9JJ(9!mkQ8Xqz9?N$boiG^oI= z_phUcY%oa%I#XOD2ikbao=%^zr;{p9*cSU>U)X>x*o1BPfG_xjZ^VFD5EEiUjEEI6 z=S$dH7PKrpND`7L70ippQ)|mkY*_;ZlX2$H`4n-J_t}?>X0Nspv%H^5u_XfO&R zeQ8hKPO6^ouayhff=$?l5BP#l_(lwf1u-Eu#E4iCGh){=ocSu8@pXbEX>DtYw|oA; z9-JJ(o}O2*7o79z6a~BM;?Gw0a%Ik!kJ3V;-gLH#BQ-y7Pp4F@<6p~pX6s=aKHv*J z;Ttg^7Q}?u5F=tm%!nNt@)dL~2)0_Zr^)!DB-=NBusH9oK69Py&*J#8Hh}Y;A}H9H zJs#{;)i=~J!I!RBIMO^d*7eNtp4l?^gm1)vSP&CpLyU+OF(YBm!t2hx zOa>KgC-W*=nD!NRV4jBp*l5mF?}dW>s_Vzf){(Jub~9<8w-dcMt-dxsvmC3k3^5=U z#Dv%oBVt9&h#eZB1)88u>l2WN0oIK-HubGDUN(G-E;Vs*W_9L%Y*WuCDOla|!&qO- z=k(7;cbb-_*>+x)yrErR67Y_{Cvmy%CmGko1>&NEk z7%?l$2wHFDkfrl?Rq~|CZE*bnF(YbC4fXreks|pk_sVjGGJx2j0a~C5+Mp3yp&8n98Mjusv%7HA ztig`L!nso0>9>KtY%{MTJMRawT~>|R{I~%We?@zb^~yYd?H3xL1)87@8le@Mp&bmg zDk-@b%%tnyGV@jws7U*fOp>Z#CUOPaV(rLgnncjm_xV5f%CZ41&;)JJ2(8cz?O>3L z$hF2*O5s?E;pJfQUP=en=(2)21S;51EuY-+!S6nxG9Dp%t2; z9Sp!iYZ8+Sft0({(A>v6QK)dCCfnFQkS+T|!PX~uvJQvW(A`!J8oh>lj^{$)ef0L{ z?CHj>j}zvx4C^#i)^Lpg+Mp3yp&8o204%^H7trUqi#Mx{jYsNVlnq=QObPeLu>HJG z6x3=IONzZntHwCf-x|unG=6Swz}w*ommTR@Syx&!zZV((>Gfg4JeFY{w!yZVC`BbfSvQyyUv-xQTpxf|Xohw$01Gey z8x7;nQ}c%!zcIbxY$TkwD8()>3SiUM^Xp=}?bw*|lgZ+S`dMG)tb@Bd(JzlYso}#& z3QvtC!>~BAii@LJ873^lI&6b&u@Clz4Rx0RXoY5I2LrGG6R-iJJWY}Yo!^T5ReJDT zIXcr~AnX2$w|zB+u`?alQ&mCz$NW^DH{%&jbf!!YCH)mc?$L1+n<3{J%K-kyHrN*X zU|-mPEp?Uw^aY?9+Q9%Uzyxf-NP{WJGt>K%cbM04e<(z*sn1-04P?E8c)#yfWp=e* zdrDFBG5;T)uga^PNNy2C0~g0;THiCvefni+53nz6z!q$(vkX8pw1WXyfC<=u5m@C} zqBP{Vbo2hQ=8wjZ;rUUl=0kqIH@7ppvSb38^XIJ8$yJ}{sj;sUO|lH4pI5|^s|s~K zecsd0V_(>SE!c!@b;b(qU;q|i0yba-RvJn}o|sM6pJA?7|E6&GxDjg>r(iu-j%BU= z9@F4rj%wwJ+>A7JAmgLMNN|oNw;beNSq@+WwqO&s)tP_T0t2uB6R-gzumZC@NtOnj z_R+JK^?&O{rN{WNVqD|HU>P&`r6<)=^Lf@YZqI9PXNoqCBBPm~_VJWC&FU{~!6t0O zhdRpu7=Q(sfDIUd6_{xt6}fKGq44jf(czPYl;Ta9#XJQoY2(9I`tK%JYxUM8^r@`8 z9H^YFKb5VaT3=7qa)3S|Y{E8tz?VA909b$t*nkmOff?B4dcxG~FtP8?((XfxQ@1XI zn72Z~reBaUxkoPwQuDF21g}Rkx^yOYohTAw)K~YFWdt^18$RGmo%sg~FaaAd0xK{B zJ2j^x*UdQ_CCEEHF2QOg1hV&MMzPwSdue7j_11mlBJrs`HGSeqK@qVko#!d@s_HM= z0N92P_)=&7VGB&a28_T8%)lgJy<5%idkztpzrH)gsk!c> zJ&qmVNJa}L)8-)c*L`Iffo=GJFLfUSzyxf-2&}*i?7&csDaln+@2^ymUE92fA_x03 z%g%mm`oc{#BgasBWF9{f#Kv9jm$g!!~@t7ksL_4FDT30xK{BJ21@Eq^Z%_ zKPEK# zaJv~@G2rj$@$KF{mIeM)2cyZ!^ATz@C&qR|UY@oshVI?fi7 z!S|)VblV^_eU&DbK5s+u&-k&EUpFhR>w+UPa3n+YZYwQSz-rtyQ%yPZ@J_4byO^`}gukQ8#?u!Rl ziN-gA4qDuINs?Eawopv+ZAo{2A43gtA^#1zj|ET0&{z%SAJ2_~6_|k?7=k63=1TI^ z;3Q;e5sB;Tjp7vFK`h`%09!mRmC~N5e$IEz!JbC1elFSvmXd!-_tg8H%v)BtBgW5`!S`A1s;R$vBpUsZ zoE=d_P#1o$LzQ+A$NkWjZOL)Fi>R>L^;ONgx6Lu5MSAjwhjplbyBI33!TukP0b5JN z(b$Jx)H_8rJ*sHVxD8l=8Q6g#Sb}M;Bu@=aX8TN+zxFT7inUd+q|{dIZhtqG&zF;r z+0$QsG7?PVg|SaYs&o$|c*<#ue;OPSPU&@}sMWC)t-&&|E-;4r1a+sK|9SWH3O`oa z+1b;fE7tVnZ7crXVmo@ehxg-nInJ0@RsY7d1uz3UFa%35RYPhtZUSTEHQ!T6ls<0F zJarYUN-aHhLs#{F|CWRHbn$gHD*7s1emDM}O7F{#z8r1-^F$IJ)xRyixgAXI12l6S zXvWLc7{3`bV|5pL!O#C!Ew`uS>n%yIU12is(n2h<^Dp6++e}eE&WtW?%6|S)PQVQ8 zzz{6KG(U_#SoW&5q`T!Q@$8;~Z28O4Y?JXak|wLp!qAc*|KxYYjl-|2+~+O6_o<1E zyQzG~+PUH&y|PrI=5$)3#`}NDzW<;Zany2n7}?EtqEbujNoSTZz1!JTG&?k0zI6Xp ziD4^ui{R-CWS-?GiW}y%Ci@(;0bmApU17$yDVwX zSO&lj48ami!B*o5%%Ov4*XIh`4>V;fA_JMu(l>M>i}(7nn}55CV@$dit|IT+I83b8 z?->Fa%ZY!fONNKWgv6{EJa z%TfNp4h+E(Ou;t2A8ng|aA%w8nzKe?H`A`HQN?kryiNiY;n(Iemg_J5o!$I{g~@NH zB(++v5<*sN7n=;dAdWL0Mln7)Zv(dQ`#$E8GpO#A-qhl8Y3gwBCqaMbn$Mn_Ty(Bz z{`0LGXsjoFAniZ{lYQBn^Uf^nj(e8Ne@^lb z9;H`CTPzR0C!du#Q|#+|T72=ulj;x4UK@aRKW%Cpjk+?8#$0fs;LxSw<8vouog>S9 zwGOHy|7Z`u5G=tIY{6LD@~(`hte2-OrEMfl9X5_VT5QdHW_QZ<{GSQ`qea;x2PFS4 zjuG|*t`k$6m7%pKqscKtKOlksEzJ2iVP1KTr~IA4mmTKs13i5=ou*FiMU@5^QfF^V zVR`X$pKbS*9Pq2I{DUD_f+^U7ae6=6G#@d%i(K$3!Uj79vdMpzV9y3;{aKGUDR#8x z$OdtD#q7Ta6XQDBBx+`NsU&2EU{GVJI4^D-O`P_>9H%abqg}E&aWrmu7+vqti|)0yAu1Uyb}OtSyVfPCg(}R|m47e@tk1~M<>A6V%@s^9Zb_I)b(+tYcsqSWJDvmEIR-1)C@aQcL6Erh|M zk4kkP_pC-#v=S=U%h~ZCOu!HzAS3&RvOCRrB>!u z)nETZ_EgHY8qIYHkOx`WezA|GN=&oY+E~miVptO(^bok)M z3G-Nnbu+4k(Auy*bi0fr_2Jj}tzOzv!0A}Ae#&W;kGr8J-WH>5&I?I)mRZ)HkrNn# zC76OO7^|W?FShwd&=Y$ea$(Z~eb~B6G4$e0*7INKh&{c(*ocZ;lL*IV<$P?M?Mp&h zoVmJGHh1tjvB4w<`cvAC(pJTj)yA!~^W63i6XvlD>#SWI$>H9|u^`g;y{LD#QY!Ue zqp*IlwaU-nw6ZdlR7-FZe;i<&IQ9zYlBrawa{@ zM*a`*HqWoF75(&JuW;2eLZ$ZhYq72PYokPeQ!Vkv6HD40=1g68hSTF?adhLx#|iUT z_FK5~|I62(!S>Wl?_Y6+-D%l$!&;wQ=lh)BS&^-zR|iZG-B0B_|6mA~U<$Ti{5ijQ zy&SQ<;brt(kG`zV&%yK}o6kskoM}U`ZQ6>jn$A|sxvqZSpUm{%Y?JGb(5GIH?Wxn8 zf%N!EH1+m~`!Hc1%X*jlUw;MVu>T-yvhk}YcHJMJr8u5>PLxF#doBJJmi6{8WBU(= zUeP{Kp6SHwu5_m}oMA>is{EaL$c~CdEfR~=@>lto`?a#~ z>Qz~Gs>?d@`Vm`NwcUXxm+wtd{yrBb%wt*lv1g_v@aGVIu1A}yJeBKqlYI4jU+Mp~ z$0eiPhxW&a@i?+$H1&aqnbNDkG$CG%wTH4{mg!%Is%4$(vesBNm zQL*%+Mp-)Ft~t20#f(F9gqzc(^oaXOpUdpdgC&@PEf|Bfw&YzIM^xfTiOmgDy872> z_LE=E@;@rkp0@ODLG$(qf@jGfm5;m1STg(V@bs8z=Z=r%hkvO{9rk}YFS~eIs5PrK znZBwhdR#iCa=#CLr32SWraL)_yVkd(#0+IBC;10cumxkV&P-!oE?aZ_q{KJLl=?r< zdH#=cZhewl(Xv5>#iF_~S>pdC1}`VhmTj`VA|8A9<@>T`zwth3Uqf2)Pm0`azNpeO zNJ_td?st>zT%IH@Z1CkWmX)%A`j7mBEf|Bf(x0~csx$v!^Wv@@4Y%7Rjyb85rS|tr zJS>u;CA;jTqM1!wO87$E>^Ou!e^o9`tQUDaOTMf2SY~>w({=fU%hkzao`88pzM>KUW7wt+b_;4%TU)|LW}jL30g0>3I>=lrIq0$G%rff7ko$aid!+ z*UEd{t3<)uGw=WX!u`%V)RJ~AI3qZX6|;2gv9Gp5DsKHneA&? zfA1~)x%yW>(yH#=GrCT8w-VQ=QTv36%i(tZB? zQdV?kw(X_KcGUAtg6RCaMV9*FndKhKrZ)+d+h0i(PhSipL!Nuw^WV;&?=x#KovK%F zPiK7WsCLzr;^httRL=in$qMs6v#f=!MvbTuziybBj_fW2L$Cx>u+`@KKe#v9l;5Ka zx-*LHnKP4aW#j%&Mo!Utc|Yy1I@C3}y0F=(W{&%Ty^3F#EexL_zT6Z_$N88}7}sa# z$T+f|H<@luaHM04Thq*cONndx)mG^od;i(_ErJZ5ON*?XF7|rbj<)kQE2D;-{0)X+ z38rAH&G&zm&wteZeH5$uV+=h~L;k_@(V4bn6YV70HoU7+PCi%nvD|q;NlAAHXQBOb zCpr-yON}?h(y|U=G_Ot<^5)NatW0f1HagFQ4w0(e@AcSYW4`5nb0IOzkPd89Z|+fU zzz{6KRNJ5bP(J^4?V1mpw{H>c{#`vh{N3%zZFn(iv@TM9#Zq78K6z&4pxBvMqaG!)NT{{hn;nkD>LN4Tg{mp1Y12|3aDxqZP9EWuRUpZ`=o|NZihF>IE85q5T) zdh!o8&+YB0_y%1X89GaTQf8Nx@MWn@V)gpv$;~UV}U7E=~$_*HTC75da`ya~pzhWD9Wiv8y zuQ$9q^}zH{jXrWy;hxxFhus(NCFP{w18=(MP5e_)L7@^J%-Cde4qlT-_#`t*aozxqbSwh4JU!-l!lg zyRxQu;Eq6P{2pcc*wj4_c3=pW`Q`m@<@gUy^~N%jeSgw}O6uzag6~>BdH!HSjf~cc zz1B{WW%c|s>EA!vycADwDa-F4v&+QXd z%r&6XUZ3u(?)-xt7=mSf8UH~!{!4rJfoz#Ub9S$6d$re*h5Q|Xr9Vo^NKsOZcn~GY z>hl3EjpoTZzIrc8mwlQ;oWpPA=kk_&r;CC^PubF#K^E2T_LeM6t}R!LYAIH_Y(Vvm zc^kyHQ+K{I+5F{4*HpWz7%VA)`1z=fgxDtm+{}gTsi)u>#T8XuwydWYT&y`U|N;; zy`TNwhMEr8CtkYVSbooS-dCT~e;++YYF)xl^b4`cwOnQ=|6m4oV3^;=e^ieDI@w|n zn|#@rC5s(0ZD+H?CF4HfyiWR0WK_?+{f9E}rFP)ynA2wEb??gmZCi=_!!mFCGuqFq zzrhUbz%aj!|Ee7SdDx@IEV(w1zfKVO{b#&Pb@N1fd}YPQGJtjf#{%>RGWH4Ymsw*p zv9-SJZ0t}`ceOP+x`(OoTp`2ZQ9RD5p zSAPFtgx_E`cY+a1eWST!fU+#0KQJuFj%GQV(dM&z#cj_v2rns2h_63W@Qs`zmhWDd zJYQ<|daW{+ocw|nn1Nk>{r&^W0KQ$)bW{`ecQFs<=n_nU=QLmM&=yqT?~e`YZAabX zThgv|mNe?FC9P;=L$3VVj;8L@|Bk|%4hD4u?KqxpVdUga^ne!&E6 zzzD3qmG3{J4Zyd|d-HFad-fg0PQI|2|Ut6D+Y&~>hO$*?^$ zj_%2NbPA^;YI;vFV_E3l&Vi=A97N@NPNmq!)2YPXNa}lJA`M*HgIv5jQJ0yMK1`FQ zccP=WeCcWDIO^vRLz5SEr;ciVZaHIY&+Squc{zX~$#lbYD z*cfJD)0jmq=uVw9P~U7vyHMihC~C*Q;jx>)<6o0c)A_eR!b3x7>j-~pVG{RYx)kD1 z7fy!K^9cSPcuBs^Gk&g`qPb%Lv_dnqg8^883D{)h^)3Eg?ovZ@AMZqA$@!Y>G#}HD zx&2OFgB-IxUTns?rTp6}(t(Cm^QJb%qNtK|CLL}*lXiCur-f_#QX^?Uin$&AVe0!^ zKiV-Th$egU?|j9^(2dj{beO*vsC+iwD_?-g}6G<8=VdxrY+w5 zJ8PYX(vzbx)TCw~+_k|5v_KQIK_j$6Gqiu} z`Td%{Qnpf~!WZ?VcOHsj;B7tEQyYnT@e%O(Ud;?_w8lVN5pbZ+K6`G;_Ykak<&^xaVG+0CZu-s@D5IZzL3p7C+G(syhe-F7; zRbvT*g{B_&R><6(%g|;cNA~`?H%rt_rNjsqYNO$6F`6g~h#9d%1GGRBv_T`Zs>+8KgC-W*=n7S2mU`?y}vE9R)GHci2bkx~_F7b0{P1IL)*hFlI5wRj>#10M6 z0!`4CRUY5t_SK?2O~x1HRcQQRvH4wnX8C+9QzW-#^7-*(@90d&H1sZ6PRamcLTrc; zu_9)~4h_%(P2XcKbE0{CoghhC+nVCio()b8ndlC;%I8M z_EhMay+*FpeXTs;{w;hX2E>Ax5F27dtcV%0fA9HIhyLK^MJ( zd|9!Zec65sJ$7i~G$7@I){yzxew1*_ksSH6KaF#6{RsQQ25iA5Y{Lh9!6$qp2E>Ax5F27tzw!c$=?PqI zciH5EowsDuwJ`bSA?wATgP&5#!|p8c;5cT@zmsfm)tR+=P?}kY3AB93c#7=MiSAT& zrfzx;bod?rwqGXKm9Y%#uno4wKG+vFU<)>38$RF*KH(cNAQr@wKE{Ii+#z6Q_};>= zOORYN^Ac@~ETs0&P1xp*p3E@HpIzKGj@_~x#@e2+W`{4;X1*cMsA}RhN(%gi47Y~- zzX{8*4%=W`?1Oz_1GZojw&4T5;1j-cM`eMRuMybuW16`|;RUiQJ3@qI((7W`=&f{n zwgH>ItOr}SW(ZSE7{{s)3t;{~`?H0IyxGnGZ^kZuoG_1NSch$}E%w2_umM}J3ES`i zUtf!-V13>Y!1#k(QoBl$(+3J8x7Kc+D<>L*`{oYwr(}tg9X)ji29BE;$h@zkeoF wezwOYxtVDl5fkRI4C}BBw#7c!H@kid_A4MzK%js?0f7Pn1q2ERXe)vL1M#3(-T(jq literal 0 HcmV?d00001 diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py index 951e0a4..8752a15 100644 --- a/python/demo_read_reconstruct.py +++ b/python/demo_read_reconstruct.py @@ -30,7 +30,7 @@ # ---------------------------------------------------------------------------- recon = np.ones(img_shape, dtype=np.float32, device=dev) -img = np.tile(np.load("../SL.npy")[..., None], (1, 1, 4)) +img = np.tile(np.load("../data/SL.npy")[..., None], (1, 1, 4)) num_iter = 2 num_subsets = 20 for it in range(num_iter): diff --git a/python/demo_write.py b/python/demo_write.py index 5943486..b2f146b 100644 --- a/python/demo_write.py +++ b/python/demo_write.py @@ -18,7 +18,7 @@ expected_num_trues = 1e6 num_iter = 2 num_subsets = 20 -np.random.seed(1) +np.random.seed(42) # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -48,7 +48,7 @@ n0, n1, n2 = img_shape # setup an image containing a box -img = np.tile(np.load("../SL.npy")[..., None], (1, 1, 4)) +img = np.tile(np.load("../data/SL.npy")[..., None], (1, 1, 4)) # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- diff --git a/python/io_yardl.py b/python/io_yardl.py index 6ba1532..f3ad203 100644 --- a/python/io_yardl.py +++ b/python/io_yardl.py @@ -1,7 +1,7 @@ from __future__ import annotations import sys -sys.path.append("/workspaces/LMSimRecon/PETSIRD/python") +sys.path.append("../PETSIRD/python") import prd import numpy as np From 8774b464b190bbc6026a1413698f50ce9f110c82 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Tue, 14 Nov 2023 23:34:30 +0000 Subject: [PATCH 13/27] improve demo scripts --- python/demo_read_reconstruct.py | 63 ++++++++++++++++++++++----------- python/demo_write.py | 42 ++++++++++++++++------ 2 files changed, 75 insertions(+), 30 deletions(-) diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py index 8752a15..ba05dda 100644 --- a/python/demo_read_reconstruct.py +++ b/python/demo_read_reconstruct.py @@ -1,10 +1,8 @@ -# TODO: - absolute scale of recon (maybe with sinogram recon) -# - additive MLEM - from __future__ import annotations import parallelproj import array_api_compat.numpy as np +from array_api_compat import to_device import matplotlib.pyplot as plt from io_yardl import read_yardl @@ -13,40 +11,65 @@ # image properties voxel_size = (2.66, 2.66, 2.66) -img_shape = (128, 128, 4) -img_origin = [-168.91, -168.91, -3.9900002] +img_shape = (128, 128, 8) +img_origin = [-168.91, -168.91, -9.31] + +num_iter = 2 +num_subsets = 20 -event_det_id_1, event_det_id_2, scanner_lut = read_yardl("write_test.yardl") +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- -# hack until we have the reader / writer implemented +event_det_id_1, event_det_id_2, scanner_lut = read_yardl("write_test.prd") xstart = scanner_lut[event_det_id_1, :] xend = scanner_lut[event_det_id_2, :] +# HACK: load the sensitivity image +sens_img = np.load("sensitivity_image.npy") + # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- # ---- LM recon using the event detector IDs and the scanner LUT ------------- # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- +res_model = parallelproj.GaussianFilterOperator( + img_shape, sigma=4.5 / (2.355 * np.asarray(voxel_size)) +) + recon = np.ones(img_shape, dtype=np.float32, device=dev) -img = np.tile(np.load("../data/SL.npy")[..., None], (1, 1, 4)) -num_iter = 2 -num_subsets = 20 + for it in range(num_iter): for isub in range(num_subsets): print(f"it {(it+1):03} / ss {(isub+1):03}", end="\r") xs_sub = xstart[isub::num_subsets, :] xe_sub = xend[isub::num_subsets, :] - exp = parallelproj.joseph3d_fwd(xs_sub, xe_sub, recon, img_origin, voxel_size) - tmp = parallelproj.joseph3d_back( + + recon_sm = res_model(recon) + + exp = parallelproj.joseph3d_fwd( + xs_sub, xe_sub, recon_sm, img_origin, voxel_size + ) + ratio_back = parallelproj.joseph3d_back( xs_sub, xe_sub, img_shape, img_origin, voxel_size, 1 / exp ) - recon *= tmp - -fig, ax = plt.subplots(1, 3, figsize=(12, 4)) -fig.colorbar(ax[0].imshow(img[:, :, 2])) -fig.colorbar(ax[1].imshow(recon[:, :, 2])) -# difference -fig.colorbar(ax[2].imshow(recon[:, :, 2] - img[:, :, 2])) -fig.savefig("read_reconstruct_test.png") + + ratio_back_sm = res_model.adjoint(ratio_back) + + recon *= ratio_back_sm / (sens_img / num_subsets) + + +# ---------------------------------------------------------------------------- +# ---------------------------------------------------------------------------- + +vmax = 0.055 +fig, ax = plt.subplots(1, recon.shape[2], figsize=(recon.shape[2] * 2, 2)) +for i in range(recon.shape[2]): + ax[i].imshow( + np.asarray(to_device(recon[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" + ) + ax[i].set_title(f"ground truth sl {i+1}", fontsize="small") + +fig.tight_layout() +fig.savefig("lm_reconstruction.png") diff --git a/python/demo_write.py b/python/demo_write.py index b2f146b..804cb34 100644 --- a/python/demo_write.py +++ b/python/demo_write.py @@ -1,6 +1,3 @@ -# TODO: - absolute scale of recon (maybe with sinogram recon) -# - additive MLEM - from __future__ import annotations import parallelproj @@ -8,8 +5,7 @@ import array_api_compat.numpy as np import matplotlib.pyplot as plt from array_api_compat import to_device -from scipy.ndimage import gaussian_filter -from io_yardl import write_yardl, read_yardl +from io_yardl import write_yardl # device variable (cpu or cuda) that determines whether calculations # are performed on the cpu or cuda gpu @@ -29,7 +25,7 @@ # setup a line of response descriptor that describes the LOR start / endpoints of # a "narrow" clinical PET scanner with 9 rings lor_descriptor = utils.DemoPETScannerLORDescriptor( - np, dev, num_rings=2, radial_trim=141 + np, dev, num_rings=4, radial_trim=141 ) # ---------------------------------------------------------------------------- @@ -48,7 +44,10 @@ n0, n1, n2 = img_shape # setup an image containing a box -img = np.tile(np.load("../data/SL.npy")[..., None], (1, 1, 4)) + +img = np.tile(np.load("../data/SL.npy")[..., None], (1, 1, num_ax)).astype(np.float32) +img[:, :, :2] = 0 +img[:, :, -2:] = 0 # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -56,7 +55,12 @@ # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- -projector = utils.RegularPolygonPETProjector(lor_descriptor, img_shape, voxel_size) +res_model = parallelproj.GaussianFilterOperator( + img_shape, sigma=4.5 / (2.355 * np.asarray(voxel_size)) +) +projector = utils.RegularPolygonPETProjector( + lor_descriptor, img_shape, voxel_size, resolution_model=res_model +) projector.tof = False # set this to True to get a time of flight projector # forward project the image @@ -117,5 +121,23 @@ # this is a 2D array of shape (num_detectors, 3) scanner_lut = lor_descriptor.scanner.all_lor_endpoints -# write and re-read -write_yardl(event_det_id_1, event_det_id_2, scanner_lut, output_file="write_test.yardl") +# write the data to PETSIRD +write_yardl(event_det_id_1, event_det_id_2, scanner_lut, output_file="write_test.prd") + +# HACK: write the sensitivity image to file +np.save("sensitivity_image.npy", sens_img) + +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- +# ----------------------------------------------------------------------------- + +vmax = 1.2 * img.max() +fig2, ax2 = plt.subplots(1, img.shape[2], figsize=(img.shape[2] * 2, 2)) +for i in range(img.shape[2]): + ax2[i].imshow( + np.asarray(to_device(img[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" + ) + ax2[i].set_title(f"ground truth sl {i+1}", fontsize="small") + +fig2.tight_layout() +fig2.savefig("simulated_phantom.png") From c09e5137b6f579f825fbb6db97e5ad8b00d02fb0 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Wed, 15 Nov 2023 00:03:41 +0000 Subject: [PATCH 14/27] clean up --- .gitignore | 3 +++ data/.gitignore | 2 ++ python/.gitignore | 2 ++ python/demo_read_reconstruct.py | 13 ++++++++----- python/demo_write.py | 20 +++++++++++++------- 5 files changed, 28 insertions(+), 12 deletions(-) create mode 100644 data/.gitignore create mode 100644 python/.gitignore diff --git a/.gitignore b/.gitignore index a0232e3..57df96b 100644 --- a/.gitignore +++ b/.gitignore @@ -5,3 +5,6 @@ ~$* *~ *.bak + +# user defined things +figs/* diff --git a/data/.gitignore b/data/.gitignore new file mode 100644 index 0000000..d5f7aab --- /dev/null +++ b/data/.gitignore @@ -0,0 +1,2 @@ +*.npy +*.prd diff --git a/python/.gitignore b/python/.gitignore new file mode 100644 index 0000000..8d221a4 --- /dev/null +++ b/python/.gitignore @@ -0,0 +1,2 @@ +*.png +*.npy diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py index ba05dda..5115f9b 100644 --- a/python/demo_read_reconstruct.py +++ b/python/demo_read_reconstruct.py @@ -5,7 +5,7 @@ from array_api_compat import to_device import matplotlib.pyplot as plt from io_yardl import read_yardl - +from pathlib import Path dev = "cpu" @@ -20,13 +20,13 @@ # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- -event_det_id_1, event_det_id_2, scanner_lut = read_yardl("write_test.prd") +event_det_id_1, event_det_id_2, scanner_lut = read_yardl("../data/write_test.prd") xstart = scanner_lut[event_det_id_1, :] xend = scanner_lut[event_det_id_2, :] # HACK: load the sensitivity image -sens_img = np.load("sensitivity_image.npy") +sens_img = np.load("../data/sensitivity_image.npy") # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -63,13 +63,16 @@ # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- +fig_dir = Path("../figs") +fig_dir.mkdir(exist_ok=True) + vmax = 0.055 fig, ax = plt.subplots(1, recon.shape[2], figsize=(recon.shape[2] * 2, 2)) for i in range(recon.shape[2]): ax[i].imshow( np.asarray(to_device(recon[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" ) - ax[i].set_title(f"ground truth sl {i+1}", fontsize="small") + ax[i].set_title(f"LM recon sl {i+1}", fontsize="small") fig.tight_layout() -fig.savefig("lm_reconstruction.png") +fig.savefig(fig_dir / "lm_reconstruction.png") diff --git a/python/demo_write.py b/python/demo_write.py index 804cb34..2d55c4c 100644 --- a/python/demo_write.py +++ b/python/demo_write.py @@ -6,6 +6,7 @@ import matplotlib.pyplot as plt from array_api_compat import to_device from io_yardl import write_yardl +from pathlib import Path # device variable (cpu or cuda) that determines whether calculations # are performed on the cpu or cuda gpu @@ -122,22 +123,27 @@ scanner_lut = lor_descriptor.scanner.all_lor_endpoints # write the data to PETSIRD -write_yardl(event_det_id_1, event_det_id_2, scanner_lut, output_file="write_test.prd") +write_yardl( + event_det_id_1, event_det_id_2, scanner_lut, output_file="../data/write_test.prd" +) # HACK: write the sensitivity image to file -np.save("sensitivity_image.npy", sens_img) +np.save("../data/sensitivity_image.npy", sens_img) # ----------------------------------------------------------------------------- # ----------------------------------------------------------------------------- # ----------------------------------------------------------------------------- +fig_dir = Path("../figs") +fig_dir.mkdir(exist_ok=True) + vmax = 1.2 * img.max() -fig2, ax2 = plt.subplots(1, img.shape[2], figsize=(img.shape[2] * 2, 2)) +fig, ax = plt.subplots(1, img.shape[2], figsize=(img.shape[2] * 2, 2)) for i in range(img.shape[2]): - ax2[i].imshow( + ax[i].imshow( np.asarray(to_device(img[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" ) - ax2[i].set_title(f"ground truth sl {i+1}", fontsize="small") + ax[i].set_title(f"ground truth sl {i+1}", fontsize="small") -fig2.tight_layout() -fig2.savefig("simulated_phantom.png") +fig.tight_layout() +fig.savefig(fig_dir / "simulated_phantom.png") From d2f9729cb7aab661e64b9d0af11a301779ee0c65 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Tue, 14 Nov 2023 20:17:19 -0800 Subject: [PATCH 15/27] clean up of PRD writer + support for TOF and energy IDS --- python/demo_read_reconstruct.py | 21 ++++-- python/demo_write.py | 27 +++++--- python/io_yardl.py | 118 +++++++++++++++++++++++--------- 3 files changed, 115 insertions(+), 51 deletions(-) diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py index 5115f9b..26b2724 100644 --- a/python/demo_read_reconstruct.py +++ b/python/demo_read_reconstruct.py @@ -8,25 +8,32 @@ from pathlib import Path dev = "cpu" +lm_data_dir: str = "../data/sim_LM_acq_1" +sens_img_file: str = "sensitivity_image.npy" +prd_file: str = "simulated_lm.prd" -# image properties +num_iter: int = 2 +num_subsets: int = 20 + +# hard coded input parameters voxel_size = (2.66, 2.66, 2.66) img_shape = (128, 128, 8) img_origin = [-168.91, -168.91, -9.31] -num_iter = 2 -num_subsets = 20 - # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- -event_det_id_1, event_det_id_2, scanner_lut = read_yardl("../data/write_test.prd") +event_det_id_1, event_det_id_2, scanner_lut = read_yardl( + str(Path(lm_data_dir) / prd_file) +) xstart = scanner_lut[event_det_id_1, :] xend = scanner_lut[event_det_id_2, :] -# HACK: load the sensitivity image -sens_img = np.load("../data/sensitivity_image.npy") +# HACK: write the sensitivity image to file +# this is currently needed since it is not agreed on how to store +# all valid detector pair combinations + attn / sens values in the PRD file +sens_img = np.load(Path(lm_data_dir) / sens_img_file).astype(np.float32) # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- diff --git a/python/demo_write.py b/python/demo_write.py index 2d55c4c..d931821 100644 --- a/python/demo_write.py +++ b/python/demo_write.py @@ -5,18 +5,20 @@ import array_api_compat.numpy as np import matplotlib.pyplot as plt from array_api_compat import to_device -from io_yardl import write_yardl +from io_yardl import write_prd_from_numpy_arrays from pathlib import Path -# device variable (cpu or cuda) that determines whether calculations -# are performed on the cpu or cuda gpu +dev: str = "cpu" +output_dir: str = "../data/sim_LM_acq_1" +output_sens_img_file: str = "sensitivity_image.npy" +output_prd_file: str = "simulated_lm.prd" +expected_num_trues: float = 1e6 -dev = "cpu" -expected_num_trues = 1e6 -num_iter = 2 -num_subsets = 20 np.random.seed(42) +# create the output directory +Path(output_dir).mkdir(exist_ok=True, parents=True) + # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- # --- setup the scanner / LOR geometry --------------------------------------- @@ -123,12 +125,17 @@ scanner_lut = lor_descriptor.scanner.all_lor_endpoints # write the data to PETSIRD -write_yardl( - event_det_id_1, event_det_id_2, scanner_lut, output_file="../data/write_test.prd" +write_prd_from_numpy_arrays( + event_det_id_1, + event_det_id_2, + scanner_lut, + output_file=str(Path(output_dir) / output_prd_file), ) # HACK: write the sensitivity image to file -np.save("../data/sensitivity_image.npy", sens_img) +# this is currently needed since it is not agreed on how to store +# all valid detector pair combinations + attn / sens values in the PRD file +np.save(Path(output_dir) / output_sens_img_file, sens_img) # ----------------------------------------------------------------------------- # ----------------------------------------------------------------------------- diff --git a/python/io_yardl.py b/python/io_yardl.py index f3ad203..3959a89 100644 --- a/python/io_yardl.py +++ b/python/io_yardl.py @@ -6,46 +6,96 @@ import numpy as np -def write_yardl( - det_1_array: np.array[int], - det_2_array: np.array[int], - scanner_lut: np.array[float], +def write_prd_from_numpy_arrays( + detector_1_id_array: np.array[int], + detector_2_id_array: np.array[int], + scanner_lut: np.array[float, float], + tof_idx_array: np.array[int] | None = None, + energy_1_idx_array: np.array[int] | None = None, + energy_2_idx_array: np.array[int] | None = None, output_file: str | None = None, + num_events: int | None = None, ) -> None: - """Writes a yardl list mode file from two detector arrays and a scanner lookup table. + """Write a PRD file from numpy arrays. Currently all into one time block - Args: - output_file (str): File to write to. - det_1_array (np.array[int]): Indices corresponding to detector 1. - det_2_array (np.array[int]): Indices corresponding to detector 2. - scanner_lut (np.array[float]): Lookup table for detector world coordinates. + Parameters + ---------- + detector_1_id_array : np.array[int] + array containing the detector 1 id for each event + detector_2_id_array : np.array[int] + array containing the detector 2 id for each event + scanner_lut : np.array[float, float] + a 2D float array of size (num_det, 3) containing the world + coordinates of each detector (in mm) + tof_idx_array : np.array[int] | None, optional + array containing the tof bin index of each event + energy_1_idx_array : np.array[int] | None, optional + array containing the energy 1 index of each event + energy_2_idx_array : np.array[int] | None, optional + array containing the energy 2 index of each event + output_file : str | None, optional + output file, if None write to stdout + num_events : int | None, optional + number of events to write, if None write all events """ - num_counts = len(det_1_array) - # create list of detectors - detectors = [] - for i in range(scanner_lut.shape[0]): - detectors.append( - prd.Detector( - id=i, x=scanner_lut[i, 0], y=scanner_lut[i, 1], z=scanner_lut[i, 2] - ) + + if num_events is None: + num_events = detector_1_id_array.shape[0] + + detector_list = [ + prd.Detector( + id=int(i), x=float(coords[0]), y=float(coords[1]), z=float(coords[2]) ) - scanner_information = prd.ScannerInformation(detectors=detectors) + for i, coords in enumerate(scanner_lut) + ] + + scanner_information = prd.ScannerInformation(detectors=detector_list) + events = [] - for i in range(num_counts): + for i in range(num_events): + det_id_1 = int(detector_1_id_array[i]) + det_id_2 = int(detector_2_id_array[i]) + + if tof_idx_array is not None: + tof_idx = int(tof_idx_array[i]) + else: + tof_idx = 0 + + if energy_1_idx_array is not None: + energy_1_idx = int(energy_1_idx_array[i]) + else: + energy_1_idx = 0 + + if energy_2_idx_array is not None: + energy_2_idx = int(energy_2_idx_array[i]) + else: + energy_2_idx = 0 + events.append( prd.CoincidenceEvent( - detector_1_id=int(det_1_array[i]), detector_2_id=int(det_2_array[i]) + detector_1_id=det_id_1, + detector_2_id=det_id_2, + tof_idx=tof_idx, + energy_1_idx=energy_1_idx, + energy_2_idx=energy_2_idx, ) ) + time_block = prd.TimeBlock(id=0, prompt_events=events) + if output_file is None: with prd.BinaryPrdExperimentWriter(sys.stdout.buffer) as writer: writer.write_header(prd.Header(scanner=scanner_information)) writer.write_time_blocks((time_block,)) else: - with prd.BinaryPrdExperimentWriter(output_file) as writer: - writer.write_header(prd.Header(scanner=scanner_information)) - writer.write_time_blocks((time_block,)) + if output_file.endswith(".ndjson"): + with prd.NDJsonPrdExperimentWriter(output_file) as writer: + writer.write_header(prd.Header(scanner=scanner_information)) + writer.write_time_blocks((time_block,)) + else: + with prd.BinaryPrdExperimentWriter(output_file) as writer: + writer.write_header(prd.Header(scanner=scanner_information)) + writer.write_time_blocks((time_block,)) def read_yardl(prd_file: str) -> tuple[np.array[int], np.array[int], np.array[float]]: @@ -77,13 +127,13 @@ def read_yardl(prd_file: str) -> tuple[np.array[int], np.array[int], np.array[fl return det_1, det_2, scanner_lut -if __name__ == "__main__": - det_1 = np.asarray([0, 1, 2, 3, 4]) - det_2 = np.asarray([5, 6, 7, 8, 9]) - scanner_lut = np.random.rand(10, 3) - write_yardl(det_1, det_2, scanner_lut, output_file="test.yardl") - det_1_read, det_2_read, scanner_lut_read = read_yardl("test.yardl") - print(scanner_lut == scanner_lut_read) - print(np.isclose(scanner_lut, scanner_lut_read).all()) - print(scanner_lut.dtype) - print(scanner_lut_read.dtype) +#if __name__ == "__main__": +# det_1 = np.asarray([0, 1, 2, 3, 4]) +# det_2 = np.asarray([5, 6, 7, 8, 9]) +# scanner_lut = np.random.rand(10, 3) +# write_yardl(det_1, det_2, scanner_lut, output_file="test.yardl") +# det_1_read, det_2_read, scanner_lut_read = read_yardl("test.yardl") +# print(scanner_lut == scanner_lut_read) +# print(np.isclose(scanner_lut, scanner_lut_read).all()) +# print(scanner_lut.dtype) +# print(scanner_lut_read.dtype) From 854532e3a734118f3747492a3124ea48f28de5b5 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Wed, 15 Nov 2023 06:22:06 -0800 Subject: [PATCH 16/27] continue to clean up reader --- python/demo_read_reconstruct.py | 12 +++-- python/io_yardl.py | 91 +++++++++++++++++++++++---------- 2 files changed, 70 insertions(+), 33 deletions(-) diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py index 26b2724..a6e99e4 100644 --- a/python/demo_read_reconstruct.py +++ b/python/demo_read_reconstruct.py @@ -4,7 +4,7 @@ import array_api_compat.numpy as np from array_api_compat import to_device import matplotlib.pyplot as plt -from io_yardl import read_yardl +from io_yardl import read_prd_to_numpy_arrays from pathlib import Path dev = "cpu" @@ -23,12 +23,14 @@ # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- -event_det_id_1, event_det_id_2, scanner_lut = read_yardl( - str(Path(lm_data_dir) / prd_file) +# TODO: used array structures with named columns + +event_attributes, scanner_lut = read_prd_to_numpy_arrays( + str(Path(lm_data_dir) / prd_file), read_tof=False, read_energy=False ) -xstart = scanner_lut[event_det_id_1, :] -xend = scanner_lut[event_det_id_2, :] +xstart = scanner_lut[event_attributes[:, 0], :] +xend = scanner_lut[event_attributes[:, 1], :] # HACK: write the sensitivity image to file # this is currently needed since it is not agreed on how to store diff --git a/python/io_yardl.py b/python/io_yardl.py index 3959a89..2dab17a 100644 --- a/python/io_yardl.py +++ b/python/io_yardl.py @@ -98,7 +98,9 @@ def write_prd_from_numpy_arrays( writer.write_time_blocks((time_block,)) -def read_yardl(prd_file: str) -> tuple[np.array[int], np.array[int], np.array[float]]: +def read_prd_to_numpy_arrays( + prd_file: str, read_tof: bool | None = None, read_energy: bool | None = None +) -> tuple[np.array[int], np.array[int], np.array[float]]: """Reads a yardl list mode file and returns two detector arrays and a scanner lookup table. Args: @@ -107,33 +109,66 @@ def read_yardl(prd_file: str) -> tuple[np.array[int], np.array[int], np.array[fl Returns: tuple[np.array[int], np.array[int], np.array[float]]: Two detector arrays and a scanner lookup table. """ + with prd.BinaryPrdExperimentReader(prd_file) as reader: # Read header and build lookup table header = reader.read_header() - scanner_information = header.scanner - detectors = scanner_information.detectors - scanner_lut = np.zeros((len(detectors), 3)) - for i, detector in enumerate(detectors): - scanner_lut[i, 0] = detector.x - scanner_lut[i, 1] = detector.y - scanner_lut[i, 2] = detector.z - # Read events - detector_hits = [] - for time_blocks in reader.read_time_blocks(): - for event in time_blocks.prompt_events: - prompt = [event.detector_1_id, event.detector_2_id] - detector_hits.append(prompt) - det_1, det_2 = np.asarray(detector_hits).T - return det_1, det_2, scanner_lut - - -#if __name__ == "__main__": -# det_1 = np.asarray([0, 1, 2, 3, 4]) -# det_2 = np.asarray([5, 6, 7, 8, 9]) -# scanner_lut = np.random.rand(10, 3) -# write_yardl(det_1, det_2, scanner_lut, output_file="test.yardl") -# det_1_read, det_2_read, scanner_lut_read = read_yardl("test.yardl") -# print(scanner_lut == scanner_lut_read) -# print(np.isclose(scanner_lut, scanner_lut_read).all()) -# print(scanner_lut.dtype) -# print(scanner_lut_read.dtype) + + # bool that decides whether the scanner has TOF and whether it is + # meaningful to read TOF + if read_tof is None: + read_tof: bool = len(header.scanner.tof_bin_edges) <= 1 + + # bool that decides whether the scanner has energy and whether it is + # meaningful to read energy + if read_energy is None: + read_energy: bool = len(header.scanner.energy_bin_edges) <= 1 + + # read the detector coordinate look up table + scanner_lut = np.array( + [[det.x, det.y, det.z] for det in header.scanner.detectors], + dtype=np.float32, + ) + + # loop over all time blocks and read all meaningful event attributes + for time_block in reader.read_time_blocks(): + if read_tof and read_energy: + event_attribute_list = [ + [ + e.detector_1_id, + e.detector_2_id, + e.tof_idx, + e.energy_1_idx, + e.energy_2_idx, + ] + for e in time_block.prompt_events + ] + elif read_tof and (not read_energy): + event_attribute_list = [ + [ + e.detector_1_id, + e.detector_2_id, + e.tof_idx, + ] + for e in time_block.prompt_events + ] + elif (not read_tof) and read_energy: + event_attribute_list = [ + [ + e.detector_1_id, + e.detector_2_id, + e.energy_1_idx, + e.energy_2_idx, + ] + for e in time_block.prompt_events + ] + else: + event_attribute_list = [ + [ + e.detector_1_id, + e.detector_2_id, + ] + for e in time_block.prompt_events + ] + + return np.array(event_attribute_list), scanner_lut From 8522e8ed2fa94af5ec594f5b28690d097def58a0 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Wed, 15 Nov 2023 17:13:21 +0000 Subject: [PATCH 17/27] make prd reader python array api agnostic --- python/demo_read_reconstruct.py | 18 ++++++++++-------- python/io_yardl.py | 15 +++++++++++---- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py index a6e99e4..d5d3d95 100644 --- a/python/demo_read_reconstruct.py +++ b/python/demo_read_reconstruct.py @@ -1,12 +1,14 @@ from __future__ import annotations import parallelproj -import array_api_compat.numpy as np from array_api_compat import to_device import matplotlib.pyplot as plt from io_yardl import read_prd_to_numpy_arrays from pathlib import Path +import array_api_compat.numpy as np +import numpy.array_api as xp + dev = "cpu" lm_data_dir: str = "../data/sim_LM_acq_1" sens_img_file: str = "sensitivity_image.npy" @@ -26,16 +28,16 @@ # TODO: used array structures with named columns event_attributes, scanner_lut = read_prd_to_numpy_arrays( - str(Path(lm_data_dir) / prd_file), read_tof=False, read_energy=False + str(Path(lm_data_dir) / prd_file), xp, dev, read_tof=False, read_energy=False ) -xstart = scanner_lut[event_attributes[:, 0], :] -xend = scanner_lut[event_attributes[:, 1], :] +xstart = xp.take(scanner_lut, event_attributes[:, 0], axis=0) +xend = xp.take(scanner_lut, event_attributes[:, 1], axis=0) # HACK: write the sensitivity image to file # this is currently needed since it is not agreed on how to store # all valid detector pair combinations + attn / sens values in the PRD file -sens_img = np.load(Path(lm_data_dir) / sens_img_file).astype(np.float32) +sens_img = xp.asarray(np.load(Path(lm_data_dir) / sens_img_file), device=dev) # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -44,10 +46,10 @@ # ---------------------------------------------------------------------------- res_model = parallelproj.GaussianFilterOperator( - img_shape, sigma=4.5 / (2.355 * np.asarray(voxel_size)) + img_shape, sigma=4.5 / (2.355 * xp.asarray(voxel_size)) ) -recon = np.ones(img_shape, dtype=np.float32, device=dev) +recon = xp.ones(img_shape, dtype=xp.float32, device=dev) for it in range(num_iter): for isub in range(num_subsets): @@ -79,7 +81,7 @@ fig, ax = plt.subplots(1, recon.shape[2], figsize=(recon.shape[2] * 2, 2)) for i in range(recon.shape[2]): ax[i].imshow( - np.asarray(to_device(recon[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" + xp.asarray(to_device(recon[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" ) ax[i].set_title(f"LM recon sl {i+1}", fontsize="small") diff --git a/python/io_yardl.py b/python/io_yardl.py index 2dab17a..5e9caf9 100644 --- a/python/io_yardl.py +++ b/python/io_yardl.py @@ -3,7 +3,9 @@ sys.path.append("../PETSIRD/python") import prd -import numpy as np + +import numpy.array_api as np +from types import ModuleType def write_prd_from_numpy_arrays( @@ -99,7 +101,11 @@ def write_prd_from_numpy_arrays( def read_prd_to_numpy_arrays( - prd_file: str, read_tof: bool | None = None, read_energy: bool | None = None + prd_file: str, + xp: ModuleType, + dev: str, + read_tof: bool | None = None, + read_energy: bool | None = None, ) -> tuple[np.array[int], np.array[int], np.array[float]]: """Reads a yardl list mode file and returns two detector arrays and a scanner lookup table. @@ -125,9 +131,10 @@ def read_prd_to_numpy_arrays( read_energy: bool = len(header.scanner.energy_bin_edges) <= 1 # read the detector coordinate look up table - scanner_lut = np.array( + scanner_lut = xp.asarray( [[det.x, det.y, det.z] for det in header.scanner.detectors], dtype=np.float32, + device=dev, ) # loop over all time blocks and read all meaningful event attributes @@ -171,4 +178,4 @@ def read_prd_to_numpy_arrays( for e in time_block.prompt_events ] - return np.array(event_attribute_list), scanner_lut + return xp.asarray(event_attribute_list, device=dev), scanner_lut From 7b15e92d8d41b8cea1d24d76e5773e586add28b7 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Wed, 15 Nov 2023 17:14:41 +0000 Subject: [PATCH 18/27] rename io module --- python/demo_read_reconstruct.py | 2 +- python/{io_yardl.py => prd_io.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename python/{io_yardl.py => prd_io.py} (100%) diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py index d5d3d95..2dce497 100644 --- a/python/demo_read_reconstruct.py +++ b/python/demo_read_reconstruct.py @@ -3,7 +3,7 @@ import parallelproj from array_api_compat import to_device import matplotlib.pyplot as plt -from io_yardl import read_prd_to_numpy_arrays +from prd_io import read_prd_to_numpy_arrays from pathlib import Path import array_api_compat.numpy as np diff --git a/python/io_yardl.py b/python/prd_io.py similarity index 100% rename from python/io_yardl.py rename to python/prd_io.py From 77552317d8dde7e712a12e25ce27d4455e18551b Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Wed, 15 Nov 2023 17:26:47 +0000 Subject: [PATCH 19/27] make read / recon example array api agnostic --- python/demo_read_reconstruct.py | 11 ++++++++++- python/prd_io.py | 2 +- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py index 2dce497..71fb8c3 100644 --- a/python/demo_read_reconstruct.py +++ b/python/demo_read_reconstruct.py @@ -7,9 +7,18 @@ from pathlib import Path import array_api_compat.numpy as np -import numpy.array_api as xp +#---------------------------------------------------------------- +#-- Choose you favorite python backend and device here ---------- +#---------------------------------------------------------------- + +# import numpy.array_api as xp +import array_api_compat.torch as xp dev = "cpu" + +#---------------------------------------------------------------- +#---------------------------------------------------------------- + lm_data_dir: str = "../data/sim_LM_acq_1" sens_img_file: str = "sensitivity_image.npy" prd_file: str = "simulated_lm.prd" diff --git a/python/prd_io.py b/python/prd_io.py index 5e9caf9..8c317a9 100644 --- a/python/prd_io.py +++ b/python/prd_io.py @@ -133,7 +133,7 @@ def read_prd_to_numpy_arrays( # read the detector coordinate look up table scanner_lut = xp.asarray( [[det.x, det.y, det.z] for det in header.scanner.detectors], - dtype=np.float32, + dtype=xp.float32, device=dev, ) From 3ecc19ddf5299587a3d12943d4d1dc5bc5ec1e63 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Thu, 16 Nov 2023 11:41:51 +0100 Subject: [PATCH 20/27] update environment.yml to have pytorch and array-api-compat --- environment.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/environment.yml b/environment.yml index f1ed0e4..2c95e6e 100644 --- a/environment.yml +++ b/environment.yml @@ -19,4 +19,5 @@ dependencies: - xtensor-fftw>=0.2.5 - xtensor>=0.24.2 - parallelproj>=1.6.1 - - scipy + - pytorch>=2.0 + - array-api-compat>=1.4 From d3282af8d5282d802aa57d766ee08a3cf22c522a Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Thu, 16 Nov 2023 11:49:10 +0100 Subject: [PATCH 21/27] correct import of prd module --- python/demo_write.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/demo_write.py b/python/demo_write.py index d931821..67ef5ab 100644 --- a/python/demo_write.py +++ b/python/demo_write.py @@ -5,7 +5,7 @@ import array_api_compat.numpy as np import matplotlib.pyplot as plt from array_api_compat import to_device -from io_yardl import write_prd_from_numpy_arrays +from prd_io import write_prd_from_numpy_arrays from pathlib import Path dev: str = "cpu" From 6dda15c3a581fcf16fc6030cfed1fe42794a5e2f Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Thu, 16 Nov 2023 12:11:42 +0100 Subject: [PATCH 22/27] clean up --- python/demo_write.py | 4 ++- python/prd_io.py | 77 ++++++++++++++++++++++++++++---------------- 2 files changed, 52 insertions(+), 29 deletions(-) diff --git a/python/demo_write.py b/python/demo_write.py index 67ef5ab..67f7831 100644 --- a/python/demo_write.py +++ b/python/demo_write.py @@ -112,7 +112,7 @@ event_det_id_1 = sino_det_start_index[event_sino_inds] event_det_id_2 = sino_det_end_index[event_sino_inds] -print(f"number of events: {event_det_id_1.shape[0]}") +print(f"number of simulated events: {event_det_id_1.shape[0]}") # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -131,11 +131,13 @@ scanner_lut, output_file=str(Path(output_dir) / output_prd_file), ) +print(f"saved PETSIRD LM file to {str(Path(output_dir) / output_prd_file)}") # HACK: write the sensitivity image to file # this is currently needed since it is not agreed on how to store # all valid detector pair combinations + attn / sens values in the PRD file np.save(Path(output_dir) / output_sens_img_file, sens_img) +print(f"saved sensitivity image to {str(Path(output_dir) / output_sens_img_file)}") # ----------------------------------------------------------------------------- # ----------------------------------------------------------------------------- diff --git a/python/prd_io.py b/python/prd_io.py index 8c317a9..5f16ebc 100644 --- a/python/prd_io.py +++ b/python/prd_io.py @@ -4,17 +4,17 @@ sys.path.append("../PETSIRD/python") import prd -import numpy.array_api as np +from collections.abc import Sequence from types import ModuleType def write_prd_from_numpy_arrays( - detector_1_id_array: np.array[int], - detector_2_id_array: np.array[int], - scanner_lut: np.array[float, float], - tof_idx_array: np.array[int] | None = None, - energy_1_idx_array: np.array[int] | None = None, - energy_2_idx_array: np.array[int] | None = None, + detector_1_id_array: Sequence[int], + detector_2_id_array: Sequence[int], + scanner_lut: Sequence[Sequence[float]], + tof_idx_array: Sequence[int] | None = None, + energy_1_idx_array: Sequence[int] | None = None, + energy_2_idx_array: Sequence[int] | None = None, output_file: str | None = None, num_events: int | None = None, ) -> None: @@ -22,18 +22,18 @@ def write_prd_from_numpy_arrays( Parameters ---------- - detector_1_id_array : np.array[int] + detector_1_id_array : Sequence[int] array containing the detector 1 id for each event - detector_2_id_array : np.array[int] + detector_2_id_array : Sequence[int] array containing the detector 2 id for each event - scanner_lut : np.array[float, float] + scanner_lut : Sequence[Sequence[float]] a 2D float array of size (num_det, 3) containing the world coordinates of each detector (in mm) - tof_idx_array : np.array[int] | None, optional + tof_idx_array : Sequence[int] | None, optional array containing the tof bin index of each event - energy_1_idx_array : np.array[int] | None, optional + energy_1_idx_array : Sequence[int] | None, optional array containing the energy 1 index of each event - energy_2_idx_array : np.array[int] | None, optional + energy_2_idx_array : Sequence[int] | None, optional array containing the energy 2 index of each event output_file : str | None, optional output file, if None write to stdout @@ -42,7 +42,9 @@ def write_prd_from_numpy_arrays( """ if num_events is None: - num_events = detector_1_id_array.shape[0] + n_events: int = len(detector_1_id_array) + else: + n_events: int = num_events detector_list = [ prd.Detector( @@ -54,7 +56,7 @@ def write_prd_from_numpy_arrays( scanner_information = prd.ScannerInformation(detectors=detector_list) events = [] - for i in range(num_events): + for i in range(n_events): det_id_1 = int(detector_1_id_array[i]) det_id_2 = int(detector_2_id_array[i]) @@ -106,16 +108,31 @@ def read_prd_to_numpy_arrays( dev: str, read_tof: bool | None = None, read_energy: bool | None = None, -) -> tuple[np.array[int], np.array[int], np.array[float]]: - """Reads a yardl list mode file and returns two detector arrays and a scanner lookup table. - - Args: - prd_file (str): File to read from. +) -> tuple[Sequence[Sequence[int]], Sequence[Sequence[float]]]: + """Read all time blocks of a PETSIRD listmode file - Returns: - tuple[np.array[int], np.array[int], np.array[float]]: Two detector arrays and a scanner lookup table. + Parameters + ---------- + prd_file : str + the PETSIRD listmode file + xp : ModuleType + the array backend module + dev : str + device used for the returned arrays + read_tof : bool | None, optional + read the TOF bin information of every event + default None means that is is auto determined + based on the scanner information (length of tof bin edges) + read_energy : bool | None, optional + read the energy information of every event + default None means that is is auto determined + based on the scanner information (length of energy bin edges) + + Returns + ------- + tuple[Sequence[Sequence[int]], Sequence[Sequence[float]]] + 2D array containing all event attributes, detector coordinate look up table """ - with prd.BinaryPrdExperimentReader(prd_file) as reader: # Read header and build lookup table header = reader.read_header() @@ -123,12 +140,16 @@ def read_prd_to_numpy_arrays( # bool that decides whether the scanner has TOF and whether it is # meaningful to read TOF if read_tof is None: - read_tof: bool = len(header.scanner.tof_bin_edges) <= 1 + r_tof: bool = len(header.scanner.tof_bin_edges) <= 1 + else: + r_tof = read_tof # bool that decides whether the scanner has energy and whether it is # meaningful to read energy if read_energy is None: - read_energy: bool = len(header.scanner.energy_bin_edges) <= 1 + r_energy: bool = len(header.scanner.energy_bin_edges) <= 1 + else: + r_energy = read_energy # read the detector coordinate look up table scanner_lut = xp.asarray( @@ -139,7 +160,7 @@ def read_prd_to_numpy_arrays( # loop over all time blocks and read all meaningful event attributes for time_block in reader.read_time_blocks(): - if read_tof and read_energy: + if r_tof and r_energy: event_attribute_list = [ [ e.detector_1_id, @@ -150,7 +171,7 @@ def read_prd_to_numpy_arrays( ] for e in time_block.prompt_events ] - elif read_tof and (not read_energy): + elif r_tof and (not r_energy): event_attribute_list = [ [ e.detector_1_id, @@ -159,7 +180,7 @@ def read_prd_to_numpy_arrays( ] for e in time_block.prompt_events ] - elif (not read_tof) and read_energy: + elif (not r_tof) and r_energy: event_attribute_list = [ [ e.detector_1_id, From 9103272d2576c30bfe2a64cfa6f926a62a43a531 Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Thu, 16 Nov 2023 12:47:03 +0100 Subject: [PATCH 23/27] bugfixes to get demo write array backend agnostic --- python/demo_read_reconstruct.py | 14 +++++----- python/demo_write.py | 40 ++++++++++++++++++--------- python/prd_io.py | 48 ++++++++++++++++++--------------- 3 files changed, 61 insertions(+), 41 deletions(-) diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py index 71fb8c3..f8ff7ce 100644 --- a/python/demo_read_reconstruct.py +++ b/python/demo_read_reconstruct.py @@ -8,16 +8,16 @@ import array_api_compat.numpy as np -#---------------------------------------------------------------- -#-- Choose you favorite python backend and device here ---------- -#---------------------------------------------------------------- +# ---------------------------------------------------------------- +# -- Choose you favorite array backend and device here -i--------- +# ---------------------------------------------------------------- + +import numpy.array_api as xp -# import numpy.array_api as xp -import array_api_compat.torch as xp dev = "cpu" -#---------------------------------------------------------------- -#---------------------------------------------------------------- +# ---------------------------------------------------------------- +# ---------------------------------------------------------------- lm_data_dir: str = "../data/sim_LM_acq_1" sens_img_file: str = "sensitivity_image.npy" diff --git a/python/demo_write.py b/python/demo_write.py index 67f7831..a68eed4 100644 --- a/python/demo_write.py +++ b/python/demo_write.py @@ -8,7 +8,17 @@ from prd_io import write_prd_from_numpy_arrays from pathlib import Path +# ---------------------------------------------------------------- +# -- Choose you favorite array backend and device here -i--------- +# ---------------------------------------------------------------- + +import numpy.array_api as xp + dev: str = "cpu" + +# ---------------------------------------------------------------- +# ---------------------------------------------------------------- + output_dir: str = "../data/sim_LM_acq_1" output_sens_img_file: str = "sensitivity_image.npy" output_prd_file: str = "simulated_lm.prd" @@ -28,7 +38,7 @@ # setup a line of response descriptor that describes the LOR start / endpoints of # a "narrow" clinical PET scanner with 9 rings lor_descriptor = utils.DemoPETScannerLORDescriptor( - np, dev, num_rings=4, radial_trim=141 + xp, dev, num_rings=4, radial_trim=141 ) # ---------------------------------------------------------------------------- @@ -48,7 +58,11 @@ # setup an image containing a box -img = np.tile(np.load("../data/SL.npy")[..., None], (1, 1, num_ax)).astype(np.float32) +img = xp.asarray( + np.tile(np.load("../data/SL.npy")[..., None], (1, 1, num_ax)), + device=dev, + dtype=xp.float32, +) img[:, :, :2] = 0 img[:, :, -2:] = 0 @@ -59,7 +73,7 @@ # ---------------------------------------------------------------------------- res_model = parallelproj.GaussianFilterOperator( - img_shape, sigma=4.5 / (2.355 * np.asarray(voxel_size)) + img_shape, sigma=4.5 / (2.355 * xp.asarray(voxel_size)) ) projector = utils.RegularPolygonPETProjector( lor_descriptor, img_shape, voxel_size, resolution_model=res_model @@ -70,13 +84,13 @@ noise_free_sinogram = projector(img) # rescale the forward projection and image such that we get the expected number of trues -scale = expected_num_trues / np.sum(noise_free_sinogram) +scale = expected_num_trues / float(xp.sum(noise_free_sinogram)) noise_free_sinogram *= scale img *= scale # calculate the sensitivity image sens_img = projector.adjoint( - np.ones(noise_free_sinogram.shape, device=dev, dtype=np.float32) + xp.ones(noise_free_sinogram.shape, device=dev, dtype=xp.float32) ) # get the two dimensional indices of all sinogram bins @@ -91,12 +105,14 @@ ) # add poisson noise to the noise free sinogram -noisy_sinogram = np.random.poisson(noise_free_sinogram) +noisy_sinogram = xp.asarray( + np.random.poisson(np.asarray(to_device(noise_free_sinogram, "cpu"))), device=dev +) # ravel the noisy sinogram and the detector start and end "index" sinograms -noisy_sinogram = noisy_sinogram.ravel() -sino_det_start_index = sino_det_start_index.ravel() -sino_det_end_index = sino_det_end_index.ravel() +noisy_sinogram = xp.reshape(noisy_sinogram, (noisy_sinogram.size,)) +sino_det_start_index = xp.reshape(sino_det_start_index, (sino_det_start_index.size,)) +sino_det_end_index = xp.reshape(sino_det_end_index, (sino_det_end_index.size,)) # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -136,7 +152,7 @@ # HACK: write the sensitivity image to file # this is currently needed since it is not agreed on how to store # all valid detector pair combinations + attn / sens values in the PRD file -np.save(Path(output_dir) / output_sens_img_file, sens_img) +np.save(Path(output_dir) / output_sens_img_file, np.asarray(to_device(sens_img, "cpu"))) print(f"saved sensitivity image to {str(Path(output_dir) / output_sens_img_file)}") # ----------------------------------------------------------------------------- @@ -146,11 +162,11 @@ fig_dir = Path("../figs") fig_dir.mkdir(exist_ok=True) -vmax = 1.2 * img.max() +vmax = 1.2 * xp.max(img) fig, ax = plt.subplots(1, img.shape[2], figsize=(img.shape[2] * 2, 2)) for i in range(img.shape[2]): ax[i].imshow( - np.asarray(to_device(img[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" + xp.asarray(to_device(img[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" ) ax[i].set_title(f"ground truth sl {i+1}", fontsize="small") diff --git a/python/prd_io.py b/python/prd_io.py index 5f16ebc..03cd4ca 100644 --- a/python/prd_io.py +++ b/python/prd_io.py @@ -4,17 +4,17 @@ sys.path.append("../PETSIRD/python") import prd -from collections.abc import Sequence +from numpy.array_api._array_object import Array from types import ModuleType def write_prd_from_numpy_arrays( - detector_1_id_array: Sequence[int], - detector_2_id_array: Sequence[int], - scanner_lut: Sequence[Sequence[float]], - tof_idx_array: Sequence[int] | None = None, - energy_1_idx_array: Sequence[int] | None = None, - energy_2_idx_array: Sequence[int] | None = None, + detector_1_id_array: Array, + detector_2_id_array: Array, + scanner_lut: Array, + tof_idx_array: Array | None = None, + energy_1_idx_array: Array | None = None, + energy_2_idx_array: Array | None = None, output_file: str | None = None, num_events: int | None = None, ) -> None: @@ -22,18 +22,18 @@ def write_prd_from_numpy_arrays( Parameters ---------- - detector_1_id_array : Sequence[int] + detector_1_id_array : Array array containing the detector 1 id for each event - detector_2_id_array : Sequence[int] + detector_2_id_array : Array array containing the detector 2 id for each event - scanner_lut : Sequence[Sequence[float]] + scanner_lut : Array a 2D float array of size (num_det, 3) containing the world coordinates of each detector (in mm) - tof_idx_array : Sequence[int] | None, optional + tof_idx_array : Array | None, optional array containing the tof bin index of each event - energy_1_idx_array : Sequence[int] | None, optional + energy_1_idx_array : Array | None, optional array containing the energy 1 index of each event - energy_2_idx_array : Sequence[int] | None, optional + energy_2_idx_array : Array | None, optional array containing the energy 2 index of each event output_file : str | None, optional output file, if None write to stdout @@ -42,16 +42,20 @@ def write_prd_from_numpy_arrays( """ if num_events is None: - n_events: int = len(detector_1_id_array) + n_events: int = detector_1_id_array.size else: n_events: int = num_events - detector_list = [ - prd.Detector( - id=int(i), x=float(coords[0]), y=float(coords[1]), z=float(coords[2]) + detector_list = [] + for i in range(scanner_lut.shape[0]): + detector_list.append( + prd.Detector( + id=int(i), + x=float(scanner_lut[i, 0]), + y=float(scanner_lut[i, 1]), + z=float(scanner_lut[i, 2]), + ) ) - for i, coords in enumerate(scanner_lut) - ] scanner_information = prd.ScannerInformation(detectors=detector_list) @@ -108,7 +112,7 @@ def read_prd_to_numpy_arrays( dev: str, read_tof: bool | None = None, read_energy: bool | None = None, -) -> tuple[Sequence[Sequence[int]], Sequence[Sequence[float]]]: +) -> tuple[Array, Array]: """Read all time blocks of a PETSIRD listmode file Parameters @@ -130,8 +134,8 @@ def read_prd_to_numpy_arrays( Returns ------- - tuple[Sequence[Sequence[int]], Sequence[Sequence[float]]] - 2D array containing all event attributes, detector coordinate look up table + tuple[Array, Array] + 2D array containing all event attributes, 2D array with detector coordinate look up table """ with prd.BinaryPrdExperimentReader(prd_file) as reader: # Read header and build lookup table From 39d8d52d5f5efc43a1b2f69b253e0549b6332a8e Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Thu, 16 Nov 2023 15:19:04 +0100 Subject: [PATCH 24/27] pass scanner_information as input to writer --- python/demo_write.py | 30 ++++++++++++++++++++++++------ python/prd_io.py | 31 ++++++------------------------- 2 files changed, 30 insertions(+), 31 deletions(-) diff --git a/python/demo_write.py b/python/demo_write.py index a68eed4..acc8e4c 100644 --- a/python/demo_write.py +++ b/python/demo_write.py @@ -1,5 +1,9 @@ from __future__ import annotations +import sys +sys.path.append("../PETSIRD/python") + +import prd import parallelproj import utils import array_api_compat.numpy as np @@ -9,7 +13,7 @@ from pathlib import Path # ---------------------------------------------------------------- -# -- Choose you favorite array backend and device here -i--------- +# -- Choose you favorite array backend and device here ----------- # ---------------------------------------------------------------- import numpy.array_api as xp @@ -140,11 +144,25 @@ # this is a 2D array of shape (num_detectors, 3) scanner_lut = lor_descriptor.scanner.all_lor_endpoints +# generate a list of detector coordinates of our scanner +detector_list = [] +for i in range(scanner_lut.shape[0]): + detector_list.append( + prd.Detector( + id=int(i), + x=float(scanner_lut[i, 0]), + y=float(scanner_lut[i, 1]), + z=float(scanner_lut[i, 2]), + ) + ) + +scanner_information = prd.ScannerInformation(detectors=detector_list) + # write the data to PETSIRD write_prd_from_numpy_arrays( event_det_id_1, event_det_id_2, - scanner_lut, + scanner_information, output_file=str(Path(output_dir) / output_prd_file), ) print(f"saved PETSIRD LM file to {str(Path(output_dir) / output_prd_file)}") @@ -165,10 +183,10 @@ vmax = 1.2 * xp.max(img) fig, ax = plt.subplots(1, img.shape[2], figsize=(img.shape[2] * 2, 2)) for i in range(img.shape[2]): - ax[i].imshow( - xp.asarray(to_device(img[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" - ) - ax[i].set_title(f"ground truth sl {i+1}", fontsize="small") + ax[i].imshow( + xp.asarray(to_device(img[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" + ) + ax[i].set_title(f"ground truth sl {i+1}", fontsize="small") fig.tight_layout() fig.savefig(fig_dir / "simulated_phantom.png") diff --git a/python/prd_io.py b/python/prd_io.py index 03cd4ca..8afa32c 100644 --- a/python/prd_io.py +++ b/python/prd_io.py @@ -11,12 +11,11 @@ def write_prd_from_numpy_arrays( detector_1_id_array: Array, detector_2_id_array: Array, - scanner_lut: Array, + scanner_information: prd.ScannerInformation, tof_idx_array: Array | None = None, energy_1_idx_array: Array | None = None, energy_2_idx_array: Array | None = None, output_file: str | None = None, - num_events: int | None = None, ) -> None: """Write a PRD file from numpy arrays. Currently all into one time block @@ -26,9 +25,9 @@ def write_prd_from_numpy_arrays( array containing the detector 1 id for each event detector_2_id_array : Array array containing the detector 2 id for each event - scanner_lut : Array - a 2D float array of size (num_det, 3) containing the world - coordinates of each detector (in mm) + scanner_information : prd.ScannerInformation + description of the scanner according to PETSIRD + (e.g. including all detector coordinates) tof_idx_array : Array | None, optional array containing the tof bin index of each event energy_1_idx_array : Array | None, optional @@ -37,30 +36,12 @@ def write_prd_from_numpy_arrays( array containing the energy 2 index of each event output_file : str | None, optional output file, if None write to stdout - num_events : int | None, optional - number of events to write, if None write all events """ - if num_events is None: - n_events: int = detector_1_id_array.size - else: - n_events: int = num_events - - detector_list = [] - for i in range(scanner_lut.shape[0]): - detector_list.append( - prd.Detector( - id=int(i), - x=float(scanner_lut[i, 0]), - y=float(scanner_lut[i, 1]), - z=float(scanner_lut[i, 2]), - ) - ) - - scanner_information = prd.ScannerInformation(detectors=detector_list) + num_events: int = detector_1_id_array.size events = [] - for i in range(n_events): + for i in range(num_events): det_id_1 = int(detector_1_id_array[i]) det_id_2 = int(detector_2_id_array[i]) From ff89fb622975fc3941471b41a39225aa42d8abfc Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Thu, 16 Nov 2023 16:38:46 +0100 Subject: [PATCH 25/27] add TOF data simulationa and reconstruction --- python/demo_read_reconstruct.py | 68 ++++++++++++++++++++++++++++----- python/demo_write.py | 60 ++++++++++++++++++++--------- python/prd_io.py | 19 +++------ 3 files changed, 108 insertions(+), 39 deletions(-) diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py index f8ff7ce..aa67b31 100644 --- a/python/demo_read_reconstruct.py +++ b/python/demo_read_reconstruct.py @@ -34,15 +34,36 @@ # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- -# TODO: used array structures with named columns +# read the LM file header and all event attributes +header, event_attributes = read_prd_to_numpy_arrays( + str(Path(lm_data_dir) / prd_file), xp, dev, read_tof=None, read_energy=False +) -event_attributes, scanner_lut = read_prd_to_numpy_arrays( - str(Path(lm_data_dir) / prd_file), xp, dev, read_tof=False, read_energy=False +# read the detector coordinates into a 2D lookup table +scanner_lut = xp.asarray( + [[det.x, det.y, det.z] for det in header.scanner.detectors], + dtype=xp.float32, + device=dev, ) xstart = xp.take(scanner_lut, event_attributes[:, 0], axis=0) xend = xp.take(scanner_lut, event_attributes[:, 1], axis=0) +# check if we have TOF data and generate the corresponding TOF parameters we need for the +# TOF joseph projector +if event_attributes.shape[1] == 3: + tof = True + event_tof_bin = event_attributes[:, 2] + num_tof_bins = header.scanner.tof_bin_edges.shape[0] - 1 + tofbin_width = header.scanner.tof_bin_edges[1] - header.scanner.tof_bin_edges[0] + sigma_tof = xp.asarray([header.scanner.tof_resolution], dtype=xp.float32) + tofcenter_offset = xp.asarray([0], dtype=xp.float32) + nsigmas = 3.0 + print(f"read {event_attributes.shape[0]} TOF events") +else: + tof = False + print(f"read {event_attributes.shape[0]} non-TOF events") + # HACK: write the sensitivity image to file # this is currently needed since it is not agreed on how to store # all valid detector pair combinations + attn / sens values in the PRD file @@ -68,12 +89,41 @@ recon_sm = res_model(recon) - exp = parallelproj.joseph3d_fwd( - xs_sub, xe_sub, recon_sm, img_origin, voxel_size - ) - ratio_back = parallelproj.joseph3d_back( - xs_sub, xe_sub, img_shape, img_origin, voxel_size, 1 / exp - ) + if tof: + event_tof_bin_sub = event_tof_bin[isub::num_subsets] - num_tof_bins // 2 + exp = parallelproj.joseph3d_fwd_tof_lm( + xs_sub, + xe_sub, + recon, + img_origin, + voxel_size, + tofbin_width, + sigma_tof, + tofcenter_offset, + nsigmas, + event_tof_bin_sub, + ) + + ratio_back = parallelproj.joseph3d_back_tof_lm( + xs_sub, + xe_sub, + img_shape, + img_origin, + voxel_size, + 1 / exp, + tofbin_width, + sigma_tof, + tofcenter_offset, + nsigmas, + event_tof_bin_sub, + ) + else: + exp = parallelproj.joseph3d_fwd( + xs_sub, xe_sub, recon_sm, img_origin, voxel_size + ) + ratio_back = parallelproj.joseph3d_back( + xs_sub, xe_sub, img_shape, img_origin, voxel_size, 1 / exp + ) ratio_back_sm = res_model.adjoint(ratio_back) diff --git a/python/demo_write.py b/python/demo_write.py index acc8e4c..4fc62c9 100644 --- a/python/demo_write.py +++ b/python/demo_write.py @@ -82,7 +82,7 @@ projector = utils.RegularPolygonPETProjector( lor_descriptor, img_shape, voxel_size, resolution_model=res_model ) -projector.tof = False # set this to True to get a time of flight projector +projector.tof = True # set this to True to get a time of flight projector # forward project the image noise_free_sinogram = projector(img) @@ -97,6 +97,13 @@ xp.ones(noise_free_sinogram.shape, device=dev, dtype=xp.float32) ) +# add poisson noise to the noise free sinogram +noisy_sinogram = xp.asarray( + np.random.poisson(np.asarray(to_device(noise_free_sinogram, "cpu"))), device=dev +) +# ravel the noisy sinogram and the detector start and end "index" sinograms +noisy_sinogram = xp.reshape(noisy_sinogram, (noisy_sinogram.size,)) + # get the two dimensional indices of all sinogram bins start_mods, end_mods, start_inds, end_inds = lor_descriptor.get_lor_indices() @@ -108,15 +115,18 @@ lor_descriptor.scanner.num_lor_endpoints_per_module[0] * end_mods + end_inds ) -# add poisson noise to the noise free sinogram -noisy_sinogram = xp.asarray( - np.random.poisson(np.asarray(to_device(noise_free_sinogram, "cpu"))), device=dev +# repeat number of TOF bin times here +num_tof_bins = projector.tof_parameters.num_tofbins + +sino_det_start_index = xp.reshape( + xp.stack(num_tof_bins * [sino_det_start_index], axis=-1), + sino_det_start_index.size * num_tof_bins, ) -# ravel the noisy sinogram and the detector start and end "index" sinograms -noisy_sinogram = xp.reshape(noisy_sinogram, (noisy_sinogram.size,)) -sino_det_start_index = xp.reshape(sino_det_start_index, (sino_det_start_index.size,)) -sino_det_end_index = xp.reshape(sino_det_end_index, (sino_det_end_index.size,)) +sino_det_end_index = xp.reshape( + xp.stack(num_tof_bins * [sino_det_end_index], axis=-1), + sino_det_end_index.size * num_tof_bins, +) # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -128,9 +138,12 @@ event_sino_inds = np.repeat(np.arange(noisy_sinogram.shape[0]), noisy_sinogram) # shuffle the event sinogram indices np.random.shuffle(event_sino_inds) +# convert event sino indices to xp array +event_sino_inds = xp.asarray(event_sino_inds, device=dev) -event_det_id_1 = sino_det_start_index[event_sino_inds] -event_det_id_2 = sino_det_end_index[event_sino_inds] +event_det_id_1 = xp.take(sino_det_start_index, event_sino_inds) +event_det_id_2 = xp.take(sino_det_end_index, event_sino_inds) +event_tof_bin = event_sino_inds % num_tof_bins print(f"number of simulated events: {event_det_id_1.shape[0]}") @@ -156,22 +169,35 @@ ) ) -scanner_information = prd.ScannerInformation(detectors=detector_list) +# setup the edges of all TOF bins +tof_bin_edges = ( + xp.arange(num_tof_bins + 1, dtype=xp.float32) - ((num_tof_bins + 1) / 2 - 0.5) +) * projector.tof_parameters.tofbin_width + +# setup the scanner information containing detector and TOF information +# WARNING: DEFINITION OF TOF RESOLUTION (sigma vs FWHM) not clear yet +scanner_information = prd.ScannerInformation( + model_name="DummyPET", + detectors=detector_list, + tof_bin_edges=np.asarray(to_device(tof_bin_edges, "cpu")), + tof_resolution=projector.tof_parameters.sigma_tof, +) # write the data to PETSIRD write_prd_from_numpy_arrays( event_det_id_1, event_det_id_2, scanner_information, + tof_idx_array=event_tof_bin, output_file=str(Path(output_dir) / output_prd_file), ) -print(f"saved PETSIRD LM file to {str(Path(output_dir) / output_prd_file)}") +print(f"wrote PETSIRD LM file to {str(Path(output_dir) / output_prd_file)}") # HACK: write the sensitivity image to file # this is currently needed since it is not agreed on how to store # all valid detector pair combinations + attn / sens values in the PRD file np.save(Path(output_dir) / output_sens_img_file, np.asarray(to_device(sens_img, "cpu"))) -print(f"saved sensitivity image to {str(Path(output_dir) / output_sens_img_file)}") +print(f"wrote sensitivity image to {str(Path(output_dir) / output_sens_img_file)}") # ----------------------------------------------------------------------------- # ----------------------------------------------------------------------------- @@ -183,10 +209,10 @@ vmax = 1.2 * xp.max(img) fig, ax = plt.subplots(1, img.shape[2], figsize=(img.shape[2] * 2, 2)) for i in range(img.shape[2]): - ax[i].imshow( - xp.asarray(to_device(img[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" - ) - ax[i].set_title(f"ground truth sl {i+1}", fontsize="small") + ax[i].imshow( + xp.asarray(to_device(img[:, :, i], "cpu")), vmin=0, vmax=vmax, cmap="Greys" + ) + ax[i].set_title(f"ground truth sl {i+1}", fontsize="small") fig.tight_layout() fig.savefig(fig_dir / "simulated_phantom.png") diff --git a/python/prd_io.py b/python/prd_io.py index 8afa32c..3836ddf 100644 --- a/python/prd_io.py +++ b/python/prd_io.py @@ -93,7 +93,7 @@ def read_prd_to_numpy_arrays( dev: str, read_tof: bool | None = None, read_energy: bool | None = None, -) -> tuple[Array, Array]: +) -> tuple[prd.types.Header, Array]: """Read all time blocks of a PETSIRD listmode file Parameters @@ -115,8 +115,8 @@ def read_prd_to_numpy_arrays( Returns ------- - tuple[Array, Array] - 2D array containing all event attributes, 2D array with detector coordinate look up table + tuple[prd.types.Header, Array] + PRD listmode file header, 2D array containing all event attributes """ with prd.BinaryPrdExperimentReader(prd_file) as reader: # Read header and build lookup table @@ -125,24 +125,17 @@ def read_prd_to_numpy_arrays( # bool that decides whether the scanner has TOF and whether it is # meaningful to read TOF if read_tof is None: - r_tof: bool = len(header.scanner.tof_bin_edges) <= 1 + r_tof: bool = len(header.scanner.tof_bin_edges) > 1 else: r_tof = read_tof # bool that decides whether the scanner has energy and whether it is # meaningful to read energy if read_energy is None: - r_energy: bool = len(header.scanner.energy_bin_edges) <= 1 + r_energy: bool = len(header.scanner.energy_bin_edges) > 1 else: r_energy = read_energy - # read the detector coordinate look up table - scanner_lut = xp.asarray( - [[det.x, det.y, det.z] for det in header.scanner.detectors], - dtype=xp.float32, - device=dev, - ) - # loop over all time blocks and read all meaningful event attributes for time_block in reader.read_time_blocks(): if r_tof and r_energy: @@ -184,4 +177,4 @@ def read_prd_to_numpy_arrays( for e in time_block.prompt_events ] - return xp.asarray(event_attribute_list, device=dev), scanner_lut + return header, xp.asarray(event_attribute_list, device=dev) From e3269c4ad2381f0f6eb6522848eac21e27d1b09b Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Thu, 16 Nov 2023 16:54:25 +0100 Subject: [PATCH 26/27] remove hard-coded voxel sizes / origin / shape in recon script --- data/.gitignore | 1 + python/demo_read_reconstruct.py | 13 ++++++------- python/demo_write.py | 9 +++++++-- 3 files changed, 14 insertions(+), 9 deletions(-) diff --git a/data/.gitignore b/data/.gitignore index d5f7aab..2cecbf9 100644 --- a/data/.gitignore +++ b/data/.gitignore @@ -1,2 +1,3 @@ *.npy +*.npz *.prd diff --git a/python/demo_read_reconstruct.py b/python/demo_read_reconstruct.py index aa67b31..d447255 100644 --- a/python/demo_read_reconstruct.py +++ b/python/demo_read_reconstruct.py @@ -20,17 +20,12 @@ # ---------------------------------------------------------------- lm_data_dir: str = "../data/sim_LM_acq_1" -sens_img_file: str = "sensitivity_image.npy" +sens_img_file: str = "sensitivity_image.npz" prd_file: str = "simulated_lm.prd" num_iter: int = 2 num_subsets: int = 20 -# hard coded input parameters -voxel_size = (2.66, 2.66, 2.66) -img_shape = (128, 128, 8) -img_origin = [-168.91, -168.91, -9.31] - # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- @@ -67,7 +62,11 @@ # HACK: write the sensitivity image to file # this is currently needed since it is not agreed on how to store # all valid detector pair combinations + attn / sens values in the PRD file -sens_img = xp.asarray(np.load(Path(lm_data_dir) / sens_img_file), device=dev) +sens_img_data = np.load(Path(lm_data_dir) / sens_img_file) +sens_img = xp.asarray(sens_img_data["sens_img"], device=dev) +img_shape = sens_img.shape +voxel_size = xp.asarray(sens_img_data["voxel_size"], device=dev) +img_origin = xp.asarray(sens_img_data["img_origin"], device=dev) # ---------------------------------------------------------------------------- # ---------------------------------------------------------------------------- diff --git a/python/demo_write.py b/python/demo_write.py index 4fc62c9..2378170 100644 --- a/python/demo_write.py +++ b/python/demo_write.py @@ -24,7 +24,7 @@ # ---------------------------------------------------------------- output_dir: str = "../data/sim_LM_acq_1" -output_sens_img_file: str = "sensitivity_image.npy" +output_sens_img_file: str = "sensitivity_image.npz" output_prd_file: str = "simulated_lm.prd" expected_num_trues: float = 1e6 @@ -196,7 +196,12 @@ # HACK: write the sensitivity image to file # this is currently needed since it is not agreed on how to store # all valid detector pair combinations + attn / sens values in the PRD file -np.save(Path(output_dir) / output_sens_img_file, np.asarray(to_device(sens_img, "cpu"))) +np.savez( + Path(output_dir) / output_sens_img_file, + sens_img=np.asarray(to_device(sens_img, "cpu")), + voxel_size=np.asarray(voxel_size), + img_origin=np.asarray(to_device(projector.img_origin, "cpu")), +) print(f"wrote sensitivity image to {str(Path(output_dir) / output_sens_img_file)}") # ----------------------------------------------------------------------------- From 9b1862ffee68d85d3c5a6637a8eb8af2dcd8ad0f Mon Sep 17 00:00:00 2001 From: Georg Schramm Date: Thu, 16 Nov 2023 16:56:17 +0100 Subject: [PATCH 27/27] rename demo script --- python/{demo_write.py => 00_simulate_lm_data.py} | 0 python/{demo_read_reconstruct.py => 01_reconstruct_lm_data.py} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename python/{demo_write.py => 00_simulate_lm_data.py} (100%) rename python/{demo_read_reconstruct.py => 01_reconstruct_lm_data.py} (100%) diff --git a/python/demo_write.py b/python/00_simulate_lm_data.py similarity index 100% rename from python/demo_write.py rename to python/00_simulate_lm_data.py diff --git a/python/demo_read_reconstruct.py b/python/01_reconstruct_lm_data.py similarity index 100% rename from python/demo_read_reconstruct.py rename to python/01_reconstruct_lm_data.py