remove tomolab reconstruction

mu_map/util.py

@@ -3,7 +3,6 @@ import random
 from typing import Optional
 
 import numpy as np
-from tomolab.Reconstruction.SPECT import SPECT_Static_Scan
 
 from mu_map.dataset.util import align_images
 
@@ -46,87 +45,3 @@ def grayscale_to_rgb(img: np.ndarray):
     assert img.ndim == 2, f"grascale image has more than 2 dimensions {img.shape}"
     return img.repeat(3).reshape((*img.shape, 3))
 
-
-def reconstruct(recon, mu_map=None, iterations=10, use_gpu=True, seed=42):
-    random.seed(seed)
-
-    spect = SPECT_Static_Scan()
-    spect.set_use_gpu(use_gpu)
-
-    recon_t = np.transpose(recon, (2, 1, 0)).astype(np.float32)
-    n_pixels = recon_t.shape[0]
-
-    padding = (n_pixels - recon_t.shape[2]) / 2
-    padding_lower = math.ceil(padding)
-    padding_upper = math.floor(padding)
-    recon_t = np.pad(recon_t, [(0, 0), (0, 0), (padding_lower, padding_upper)])
-
-    spect.set_n_pixels(n_pixels, n_pixels)
-    spect.set_gantry_angular_positions(0.0, 360.0, 59)
-
-    measurement = spect.project(recon_t)
-    spect.set_measurement(measurement.data)
-
-    if mu_map is not None:
-        mu_map_t = np.transpose(mu_map, (2, 1, 0)).astype(np.float32)
-        mu_map_t = np.pad(mu_map_t, [(0, 0), (0, 0), (padding_lower, padding_upper)])
-        spect.set_attenuation(mu_map_t)
-
-    spect.set_pixel_size(4.8, 4.8)
-    spect.set_radius(200.0)
-    spect.set_psf(fwhm0_mm=5.0, depth_dependence=0.0001)
-
-    activity = spect.estimate_activity(
-        iterations=iterations, subset_size=16, subset_mode="random", method="EM"
-    )
-
-    activity = activity.data
-    activity = np.transpose(activity, (2, 1, 0))
-    activity, _ = align_images(activity, recon)
-    return activity
-
-if __name__ == "__main__":
-    import pydicom
-    import time
-
-    import cv2 as cv
-    
-    from mu_map.dataset.util import load_dcm_img
-
-    recon_nac = load_dcm_img("./data/second/images/0001-stress-recon_nac_nsc.dcm")
-    mu_map = load_dcm_img("./data/second/images/0001-stress-mu_map.dcm")
-    recon_nac, mu_map = align_images(recon_nac, mu_map)
-
-    since = time.time()
-    recon_ac = reconstruct(recon_nac)
-    took = time.time() - since
-    print(f"Reconstruction took {took:.3f}s")
-
-    print(f"            |   Max |  Mean |   Min")
-    print(f"  Recon nAC | {recon_nac.max():.3f} | {recon_nac.mean():.3f} | {recon_nac.min():.3f}")
-    print(f"  Recon  AC | {recon_ac.max():.3f} | {recon_ac.mean():.3f} | {recon_ac.min():.3f}")
-
-    def display(recon_nac, recon_ac, _slice):
-        img1 = to_grayscale(recon_nac[_slice], max_val=recon_nac.max())
-        img1 = cv.resize(img1, (1024, 1024))
-
-        img2 = to_grayscale(recon_ac[_slice], max_val=recon_ac.max())
-        img2 = cv.resize(img2, (1024, 1024))
-
-        s = np.full((1024, 10), 239, np.uint8)
-        return np.hstack((img1, s, img2))
-    exit(0)
-
-    i = 0
-    wname = "Test"
-    cv.namedWindow(wname, cv.WINDOW_NORMAL)
-    cv.resizeWindow(wname, 1600, 900)
-    while True:
-        cv.imshow(wname, display(recon_nac, recon_ac, i))
-
-        i = (i + 1) % recon_nac.shape[0]
-
-        key = cv.waitKey(100)
-        if key == ord("q"):
-            break
-