reconstruction working

mu_map/util.py

 import math
+import random
 from typing import Optional
 
 import numpy as np
@@ -46,16 +47,16 @@ def grayscale_to_rgb(img: np.ndarray):
     return img.repeat(3).reshape((*img.shape, 3))
 
 
-import random
 def reconstruct(recon, mu_map=None, use_gpu=True, seed=42):
-    random.seed(42)
+    random.seed(seed)
 
+    print("Init scan")
     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]
-    print(f"NPixels {n_pixels}, {recon_t.shape}")
+    print(f"Set to NPixels {n_pixels}, {recon_t.shape}")
     # spect.set_n_pixels(n_pixels, n_pixels)
     spect.set_n_pixels(128, 128)
     spect.set_gantry_angular_positions(0.0, 360.0, 59)
@@ -67,6 +68,7 @@ def reconstruct(recon, mu_map=None, use_gpu=True, seed=42):
         padding_upper = math.floor(n_pixels - mu_map_t.shape[2])
         mu_map_t = np.pad(mu_map_t, [(0, 0), (0, 0), (padding_lower, padding_upper)])
 
+        print("Set attenuation map...")
         spect.set_attenuation(mu_map_t)
 
     print("Forward projection .....", end="\r")
@@ -75,6 +77,7 @@ def reconstruct(recon, mu_map=None, use_gpu=True, seed=42):
     spect.set_measurement(measurement.data)
     print("Forward projection DONE!")
 
+    print("init more parameters...")
     spect.set_pixel_size(4.8, 4.8)
     spect.set_radius(200.0)
     spect.set_psf(fwhm0_mm=5.0, depth_dependence=0.0001)
@@ -87,15 +90,52 @@ def reconstruct(recon, mu_map=None, use_gpu=True, seed=42):
 
     activity = activity.data
     activity = np.transpose(activity, (2, 1, 0))
+    _shape = activity.shape
     activity, _ = align_images(activity, recon)
+    print(f"Realign activity from {_shape} to {activity.shape}")
     return activity
 
 if __name__ == "__main__":
     import pydicom
+    import time
+
+    import cv2 as cv
     
     from mu_map.dataset.util import load_dcm_img
 
-    img = load_dcm_img("./data/second/images/0001-stress-recon_nac_nsc.dcm")
+    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))
+
+    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
 
-    reconstruct(img, mu_map)