From 702f42088c3d106d57aa53d506d98d75dafb3289 Mon Sep 17 00:00:00 2001 From: econdami Date: Tue, 9 Jul 2024 19:05:49 +0200 Subject: [PATCH] Add Make_AIF brick to Perfdsc pipeline --- mia_processes/bricks/tools/tools.py | 106 ++++++++++++++----- mia_processes/pipelines/Perfusion/perfdsc.py | 21 +++- 2 files changed, 100 insertions(+), 27 deletions(-) diff --git a/mia_processes/bricks/tools/tools.py b/mia_processes/bricks/tools/tools.py index 693aabfb..9cd09bf5 100644 --- a/mia_processes/bricks/tools/tools.py +++ b/mia_processes/bricks/tools/tools.py @@ -39,6 +39,7 @@ # Other imports import csv +import json import os import re import shutil @@ -2168,7 +2169,7 @@ def __init__(self): def bol_ar_time(self, data): """Compute bolus arrival time - :param data: + :param data: Value of voxels selected during dynamics :returns: the bolus arrival time """ window_size = 8 @@ -2193,6 +2194,28 @@ def bol_ar_time(self, data): return t0 + def convert_to_native_types(self, data): + """Convert data to list of native Python types + + Useful, for example, for serializing with json. + + :param data: The data to be converted to a native Python type + """ + + if isinstance(data, np.ndarray): + return data.tolist() + + if isinstance(data, list): + return [self.convert_to_native_types(item) for item in data] + + if isinstance(data, np.int64): + return int(data) + + if isinstance(data, np.float64): + return float(data) + + return data + def list_outputs(self, is_plugged=None): """Dedicated to the initialisation step of the brick. The main objective of this method is to produce the outputs of the @@ -2220,7 +2243,7 @@ def list_outputs(self, is_plugged=None): _, file_name = os.path.split(self.func_file) file_name_no_ext, _ = os.path.splitext(file_name) self.outputs["aif_file"] = os.path.join( - self.output_directory, file_name_no_ext + "_aif.mat" + self.output_directory, file_name_no_ext + "_aif.json" ) else: @@ -2404,40 +2427,56 @@ def run_process_mia(self): for i in range(nb_vox_cand): # calculation of arrival time t0 = self.bol_ar_time(tmp_data[sorting[i], :]) - initslop = delta_base_min.flatten()[sorting[i]] / ( - idx_min.flatten()[sorting[i]] - t0 - ) - score[i] = (delta_base_min.flatten()[sorting[i]] * initslop) / ( - half_width_peak.flatten()[sorting[i]] * t0 + initslop = delta_base_min.flatten(order="F")[sorting[i]] / ( + idx_min.flatten(order="F")[sorting[i]] - t0 ) + # TODO: Do we need * (t0 + 1) because python is 0-based indexing ? + score[i] = ( + delta_base_min.flatten(order="F")[sorting[i]] * initslop + ) / (half_width_peak.flatten(order="F")[sorting[i]] * t0) sorted_score = np.argsort(score)[::-1] + # first-pass curve averaged over the 5 best voxels determined using + # the score. aif = np.mean(tmp_data[sorting[sorted_score[:nb_vox]], :], axis=0) - roi = np.zeros((nrow, ncol, nslices), dtype=bool) - - for i in range(nb_vox): - roi.flat[sorting[sorted_score[i]]] = True + # these lines are unnecessary: + # roi = np.zeros((nrow, ncol, nslices), dtype=bool) + # + # for i in range(nb_vox): + # roi.flat[sorting[sorted_score[i]]] = True # i, j, k correspond to the indexes in the brain (3D) for # the best results - i, j, k = np.where(roi) + i, j, k = np.unravel_index( + sorting[sorted_score[:nb_vox]], (nrow, ncol, nslices), order="F" + ) for v in range(nb_vox): warn = 0 + # In scores, the first 4 columns are: + # score value / row index / column index / slice index. scores[v][0:4] = [score[sorted_score[v]], i[v], j[v], k[v]] + # Then, from the sixth to the ninth column, the results of 4 + # tests are displayed. - # Pre-bolus baseline test - std_basepre = np.std(tmp_data[sorting[sorted_score[v]], 1:7]) - mean_basepre = np.mean(tmp_data[sorting[sorted_score[v]], 1:7]) + # 1/ Pre-bolus baseline test: indicates whether the baseline is + # too noisy + std_basepre = np.std( + tmp_data[sorting[sorted_score[v]], 1:6], ddof=1 + ) + mean_basepre = np.mean(tmp_data[sorting[sorted_score[v]], 1:6]) if std_basepre >= mean_basepre / 10: warn += 1 - scores[v][4 + warn] = ( - "The pre-bolus baseline of the voxel " "is too noisy" - ) - - # Post-bolus baseline test - std_basepost = np.std(tmp_data[sorting[sorted_score[v]], -6:]) + scores[v][ + 4 + warn + ] = "The pre-bolus baseline of the voxel is too noisy" + + # 2/ Post-bolus baseline test: indicates whether the baseline is + # too noisy + std_basepost = np.std( + tmp_data[sorting[sorted_score[v]], -6:], ddof=1 + ) mean_basepost = np.mean(tmp_data[sorting[sorted_score[v]], -6:]) if std_basepost >= mean_basepost / 10: @@ -2446,22 +2485,39 @@ def run_process_mia(self): "The post-bolus baseline of the voxel " "is too noisy" ) - # t0 point test + # 3/ t0 point test: indicates whether the voxel value at the time + # of bolus arrival (t0) is greater than 11/10th of the baseline + # mean t0 = self.bol_ar_time(tmp_data[sorting[sorted_score[v]], :]) if tmp_data[sorting[sorted_score[v]], t0] >= 1.1 * mean_basepre: warn += 1 scores[v][4 + warn] = "The voxel value at t0 is too high" - # Pre-bolus baseline length test + # 4/ Pre-bolus baseline length test if t0 < 8: warn += 1 scores[v][4 + warn] = "The voxel baseline is too short" + # Finally, the fifth column shows the number of warnings scores[v][4] = warn + + # Saving results ################################################### - # TODO: May be we don't need to save in matlab file format? - io.savemat(self.aif_file, {"aif": aif, "scores": scores}) + with open(self.aif_file, "w") as f: + json.dump( + { + "aif": self.convert_to_native_types(aif), + "scores": self.convert_to_native_types(scores), + }, + f, + ) + + # To load the data back: + # with open('pathTofile.json', 'r') as f: + # data = json.load(f) + # aif = np.array(data['aif']) # or just aif = data['aif'] + # scores = data['scores'] print(f"Make_AIF brick: {self.aif_file} created!") diff --git a/mia_processes/pipelines/Perfusion/perfdsc.py b/mia_processes/pipelines/Perfusion/perfdsc.py index 41956a5a..8287cb28 100644 --- a/mia_processes/pipelines/Perfusion/perfdsc.py +++ b/mia_processes/pipelines/Perfusion/perfdsc.py @@ -52,6 +52,9 @@ def pipeline_definition(self): self.nodes["1_spatial_preprocessing"].process.nodes[ "realign" ].process.register_to_mean = False + # We can do the same with: + # self.nodes["1_spatial_preprocessing"].process.nodes[ + # "realign"].set_plug_value("register_to_mean", False) self.nodes["1_spatial_preprocessing"].process.nodes[ "normalize12_1" ].process.write_voxel_sizes = [1.0, 1.0, 1.0] @@ -87,6 +90,9 @@ def pipeline_definition(self): "threshold_2": True, "resample1": True, } + self.add_process( + "make_aif", "mia_processes.bricks.tools.tools.Make_AIF" + ) # links self.export_parameter( @@ -112,16 +118,25 @@ def pipeline_definition(self): "coregistered_source", is_optional=False, ) + self.add_link( + "1_spatial_preprocessing.smoothed_func->" "make_aif.func_file" + ) # Only done for normalize12_2 activation! self.export_parameter( "1_spatial_preprocessing", "normalized_anat", is_optional=False ) self.export_parameter("2_spatial_mask", "mask_003", is_optional=False) + self.export_parameter("make_aif", "aif_file", is_optional=False) # parameters order - self.reorder_traits( - ("mask_003", "anat_file", "func_files", "coregistered_source") + ( + "mask_003", + "anat_file", + "func_files", + "coregistered_source", + "aif_file", + ) ) # nodes positions @@ -130,6 +145,7 @@ def pipeline_definition(self): "2_spatial_mask": (42.0, -76.0), "outputs": (337.00279291855725, -76.0), "inputs": (-548.5783216389599, -167.0), + "make_aif": (43.1492517677367, -270.66348836125746), } # nodes dimensions @@ -138,6 +154,7 @@ def pipeline_definition(self): "2_spatial_mask": (200.15625, 145.0), "outputs": (135.8125, 110.0), "inputs": (92.79751223929541, 110.0), + "make_aif": (117.0, 75.0), } self.do_autoexport_nodes_parameters = False