rename recon module to osem and add postfilter support

mu_map/recon/recon.py → mu_map/recon/osem.py

@@ -16,6 +16,7 @@ from mu_map.file.interfile import (
     KEY_SPACING_2,
     TEMPLATE_HEADER_IMAGE,
 )
+from mu_map.recon.filter import GaussianFilter
 
 
 TEMPLATE_RECON_PARAMS = """
@@ -65,21 +66,6 @@ OSMAPOSLParameters :=
   number of subiterations:= {N_SUBITERATIONS}
   Save estimates at subiteration intervals:= {SAVE_INTERVALS}
 
-  ; keywords that specify the filtering that occurs after every subiteration
-  ; warning: do not normally use together with a prior
-  inter-iteration filter subiteration interval := 4
-  inter-iteration filter type := Separable Gaussian
-  ;post-filter type := Separable Gaussian
-  separable gaussian filter parameters :=
-    x-dir filter fwhm (in mm) := 9.6
-    y-dir filter fwhm (in mm) := 9.6
-    z-dir filter fwhm (in mm) := 9.6
-    x-dir maximum kernel size := 128
-    y-dir maximum kernel size := 128
-    z-dir maximum kernel size := 33
-    Normalise filter to 1 := 1
-  end separable gaussian filter parameters :=
-
 END :=
 """
 
@@ -112,6 +98,7 @@ def reconstruct(
     init: Optional[Tuple[Dict[str, str], np.ndarray]] = None,
     n_subsets: Optional[int] = 4,
     n_iterations: Optional[int] = 10,
+    postfilter: Optional[GaussianFilter] = None,
     **kwargs,
 ):
     # sanitize parameters
@@ -155,6 +142,9 @@ def reconstruct(
     write_interfile(filename_init, *init)
     params = params.replace("{INIT_FILE}", filename_init)
 
+    if postfilter:
+        params = postfilter.insert(params)
+
     filename_params = os.path.join(dir_tmp.name, "OSEM_SPECT.par")
     with open(filename_params, mode="w") as f:
         f.write(params.strip())
@@ -162,7 +152,6 @@ def reconstruct(
     recon = stir.OSMAPOSLReconstruction3DFloat(filename_params)
     recon.reconstruct()
 
-    print(params)
     return load_interfile(filename_out)
 
 
@@ -209,6 +198,18 @@ if __name__ == "__main__":
         default=10,
         help="the number of iterations for OSEM reconstruction",
     )
+    parser.add_argument(
+        "--postfilter",
+        choices=["gaussian", "none"],
+        default="gaussian",
+        help="apply a postfilter to the reconstruction",
+    )
+    parser.add_argument(
+        "--postfilter_width",
+        type=float,
+        default=1.0,
+        help="the filter witdth for the postfilter is based on the spacing in the projection and can be modified with this factor",
+    )
     parser.add_argument(
         "-v", "--verbosity", type=int, default=0, help="configure the verbosity of STIR"
     )
@@ -230,11 +231,19 @@ if __name__ == "__main__":
     else:
         projection = load_as_interfile(args.projection)
 
+    postfilter = None
+    if args.postfilter == "gaussian":
+        postfilter = GaussianFilter(projection, args.postfilter_width)
+
     kwargs = vars(args)
     del kwargs["mu_map"]
     del kwargs["projection"]
-    print(kwargs)
-    header, image = reconstruct(projection, mu_map=mu_map, **kwargs)
-    image = image[:, :, ::-1]
+    del kwargs["postfilter"]
 
+    header, image = reconstruct(
+        projection, mu_map=mu_map, postfilter=postfilter, **kwargs
+    )
+    image = image[
+        :, :, ::-1
+    ]  # STIR creates reconstructions flipped on the x-axes, this is reverted here
     write_interfile(args.out, header, image)