diff --git a/python/recon_block_scanner_listmode_data.py b/python/recon_block_scanner_listmode_data.py index 986e04a..c78eaed 100644 --- a/python/recon_block_scanner_listmode_data.py +++ b/python/recon_block_scanner_listmode_data.py @@ -166,7 +166,7 @@ def transform_BoxShape( assert num_tofbins % 2 == 1, "Number of TOF bins must be odd" # %% # calculate the sensitivity image - +print("Calculating sensitivity image") # we loop through the symmetric group ID look up table to see which module pairs # are in coincidence @@ -185,7 +185,7 @@ def transform_BoxShape( sgid = header.scanner.detection_efficiencies.module_pair_sgidlut[i, j] if sgid >= 0: - print(i, j, sgid) + print(f"mod1 {i:03}, mod2 {j:03}, SGID {sgid:03}", end="\r") start_det_el = det_element_center_list[i] end_det_el = det_element_center_list[j] @@ -200,8 +200,6 @@ def transform_BoxShape( ) proj.tof_parameters = tof_params - # TODO: add TOF parameters - # get the module pair efficiencies - asumming that we only use 1 energy bin module_pair_eff = ( header.scanner.detection_efficiencies.module_pair_efficiencies_vector[ @@ -212,7 +210,6 @@ def transform_BoxShape( start_el_eff = xp.repeat(det_el_efficiencies[i], len(end_det_el), axis=0) end_el_eff = xp.tile(det_el_efficiencies[j], (len(start_det_el))) - # TODO loop over TOF! for tofbin in xp.arange(-(num_tofbins // 2), num_tofbins // 2 + 1): # print(tofbin) proj.event_tofbins = xp.full( @@ -220,6 +217,8 @@ def transform_BoxShape( ) sens_img += proj.adjoint(start_el_eff * end_el_eff * module_pair_eff) +print("") + # for some reason we have to divide the sens image by the number of TOF bins # right now unclear why that is sens_img /= num_tofbins @@ -227,46 +226,64 @@ def transform_BoxShape( # %% # read all coincidence events +print("Reading LM events") -if False: - num_prompts = 0 - event_counter = 0 - num_tof_bins = header.scanner.number_of_tof_bins() +num_prompts = 0 +event_counter = 0 +num_tof_bins = header.scanner.number_of_tof_bins() - xstart = [] - xend = [] - tof_bin = [] - effs = [] +xstart = [] +xend = [] +tof_bin = [] +effs = [] - for i_time_block, time_block in enumerate(reader.read_time_blocks()): - if isinstance(time_block, petsird.TimeBlock.EventTimeBlock): - num_prompts += len(time_block.value.prompt_events) +for i_time_block, time_block in enumerate(reader.read_time_blocks()): + if isinstance(time_block, petsird.TimeBlock.EventTimeBlock): + num_prompts += len(time_block.value.prompt_events) - for i_event, event in enumerate(time_block.value.prompt_events): - event_mods_and_els = get_module_and_element( - header.scanner.scanner_geometry, event.detector_ids - ) + for i_event, event in enumerate(time_block.value.prompt_events): + event_mods_and_els = get_module_and_element( + header.scanner.scanner_geometry, event.detector_ids + ) + + event_start_coord = det_element_center_list[event_mods_and_els[0].module][ + event_mods_and_els[0].el + ] + xstart.append(event_start_coord) - event_start_coord = element_centers[ - event_mods_and_els[0].module, event_mods_and_els[0].el - ] - xstart.append(event_start_coord) + event_end_coord = det_element_center_list[event_mods_and_els[1].module][ + event_mods_and_els[1].el + ] + xend.append(event_end_coord) - event_end_coord = element_centers[ - event_mods_and_els[1].module, event_mods_and_els[1].el - ] - xend.append(event_end_coord) + # get the event efficiencies + effs.append(get_detection_efficiency(header.scanner, event)) + # get the signed event TOF bin (0 is the central bin) + tof_bin.append(event.tof_idx - num_tof_bins // 2) - # get the event efficiencies - effs.append(get_detection_efficiency(header.scanner, event)) - # get the signed event TOF bin (0 is the central bin) - tof_bin.append(event.tof_idx - num_tof_bins // 2) + event_counter += 1 + +reader.close() + +xstart = xp.asarray(xstart, device=dev) +xend = xp.asarray(xend, device=dev) +effs = xp.asarray(effs, device=dev) +tof_bin = xp.asarray(tof_bin, device=dev) + + +# %% +# run a LM OSEM recon - event_counter += 1 +num_iter = 50 +recon = xp.ones(img_shape, dtype="float32") - reader.close() +proj = parallelproj.ListmodePETProjector(xstart, xend, img_shape, voxel_size) +proj.tof_parameters = tof_params +proj.event_tofbins = tof_bin - xstart = xp.asarray(xstart, device=dev) - xend = xp.asarray(xend, device=dev) - effs = xp.asarray(effs, device=dev) - tof_bin = xp.asarray(tof_bin, device=dev) +for i in range(num_iter): + print(f"it {(i +1):03} / {num_iter:03}") + lm_exp = effs * proj(res_model(recon)) + print(lm_exp.min()) + tmp = res_model(proj.adjoint(effs / lm_exp)) + recon *= tmp / sens_img diff --git a/python/simulate_block_scanner_listmode_data.py b/python/simulate_block_scanner_listmode_data.py index e8c78af..8a41232 100755 --- a/python/simulate_block_scanner_listmode_data.py +++ b/python/simulate_block_scanner_listmode_data.py @@ -220,6 +220,12 @@ def module_pair_eff_from_sgd(i_sgd: int) -> float: det_el_efficiencies = 0.2 + 2 * xp.astype( xp.random.rand(scanner.num_modules, lor_desc.num_lorendpoints_per_block), "float32" ) +# multiply the det el eff. of the first module by 3 to introduce more variation +det_el_efficiencies[0, :] *= 3 + +# divide the det el eff. of the last module by 3 to introduce more variation +det_el_efficiencies[-1, :] /= 3 + # simulate a few dead crystals det_el_efficiencies[det_el_efficiencies < 0.21] = 0 @@ -257,7 +263,6 @@ def module_pair_eff_from_sgd(i_sgd: int) -> float: # and normalization are ignored) ones_back_tof = fwd_op.adjoint(xp.ones(fwd_op.out_shape, dtype=xp.float32, device=dev)) -xp.save("ones_back_tof.npy", ones_back_tof) print(ones_back_tof.shape) # %% @@ -421,7 +426,6 @@ def module_pair_eff_from_sgd(i_sgd: int) -> float: tofbin=event_tof_bin, ) - xp.save("lm_back.npy", lm_back) lm_back_non_tof = parallelproj.joseph3d_back( xstart=scanner.get_lor_endpoints(event_start_block, event_start_el), xend=scanner.get_lor_endpoints(event_end_block, event_end_el), @@ -430,7 +434,6 @@ def module_pair_eff_from_sgd(i_sgd: int) -> float: voxsize=proj.voxel_size, img_fwd=xp.ones(num_events, dtype=xp.float32, device=dev), ) - xp.save("lm_back_non_tof.npy", lm_back_non_tof) vi = pv.ThreeAxisViewer([histo_back, lm_back, histo_back - lm_back])