diff --git a/documentation/tutorial/visualization.py b/documentation/tutorial/visualization.py index 7010636..0168171 100644 --- a/documentation/tutorial/visualization.py +++ b/documentation/tutorial/visualization.py @@ -38,6 +38,7 @@ from medpy.metric import binary import ntpath +from scipy.ndimage import label def superimpose_segmentation_images(pet_gt_prd_display, file_name, logzscore=None): @@ -176,8 +177,6 @@ def read_predicted_images(path: str = None): pred[pred>0.5] =1 pred[pred<0.5] = 0 - print(pet.shape) - for coronal_sagittal in range(2): if len(gt) and len(pred): pet_gt_prd_display = [pet[coronal_sagittal], gt[coronal_sagittal], pred[coronal_sagittal]] diff --git a/src/LFBNet/losses/__pycache__/__init__.cpython-36.pyc b/src/LFBNet/losses/__pycache__/__init__.cpython-36.pyc deleted file mode 100644 index 570bdd7..0000000 Binary files a/src/LFBNet/losses/__pycache__/__init__.cpython-36.pyc and /dev/null differ diff --git a/src/LFBNet/losses/__pycache__/losses.cpython-36.pyc b/src/LFBNet/losses/__pycache__/losses.cpython-36.pyc deleted file mode 100644 index faba8c2..0000000 Binary files a/src/LFBNet/losses/__pycache__/losses.cpython-36.pyc and /dev/null differ diff --git a/src/LFBNet/network_architecture/__pycache__/__init__.cpython-36.pyc b/src/LFBNet/network_architecture/__pycache__/__init__.cpython-36.pyc deleted file mode 100644 index d2e85ac..0000000 Binary files a/src/LFBNet/network_architecture/__pycache__/__init__.cpython-36.pyc and /dev/null differ diff --git a/src/LFBNet/network_architecture/__pycache__/get_conv_blocks.cpython-36.pyc b/src/LFBNet/network_architecture/__pycache__/get_conv_blocks.cpython-36.pyc deleted file mode 100644 index 548c98a..0000000 Binary files a/src/LFBNet/network_architecture/__pycache__/get_conv_blocks.cpython-36.pyc and /dev/null differ diff --git a/src/LFBNet/network_architecture/__pycache__/lfbnet.cpython-36.pyc b/src/LFBNet/network_architecture/__pycache__/lfbnet.cpython-36.pyc deleted file mode 100644 index d27037c..0000000 Binary files a/src/LFBNet/network_architecture/__pycache__/lfbnet.cpython-36.pyc and /dev/null differ diff --git a/src/LFBNet/postprocessing/__init__.py b/src/LFBNet/postprocessing/__init__.py index 4fa414a..f4e9742 100644 --- a/src/LFBNet/postprocessing/__init__.py +++ b/src/LFBNet/postprocessing/__init__.py @@ -1,4 +1,4 @@ from __future__ import absolute_import from src import * -from . import * +from .postprocessing import * diff --git a/src/LFBNet/postprocessing/postprocessing.py b/src/LFBNet/postprocessing/postprocessing.py new file mode 100644 index 0000000..3008a13 --- /dev/null +++ b/src/LFBNet/postprocessing/postprocessing.py @@ -0,0 +1,68 @@ +""" Script to preprocess a given FDG PET image in .nii. +""" +# Import libraries +import glob +import os + +import numpy as np +from numpy import ndarray +from numpy.random import seed +from tqdm import tqdm +from typing import List, Tuple +import matplotlib.pyplot as plt +import warnings + +from skimage.transform import resize +import scipy.ndimage +import nibabel as nib +from scipy.ndimage import label + +# seed random number generator +seed(1) + + +def remove_outliers_in_sagittal(predicted): + """ when the sagittal image based segmentation does not have corresponding image in coronal remove it + """ + sagittal = np.squeeze(predicted[0, ...]) + coronal = np.squeeze(predicted[1, ...]) + + try: + sagittal = np.squeeze(sagittal) + except: + pass + + try: + coronal = np.squeeze(coronal) + except: + pass + + binary_mask = sagittal.copy() + binary_mask[binary_mask >= 0.5] = 1 + binary_mask[binary_mask < 0.5] = 0 + coronal[coronal >= 0.5] = 1 + coronal[coronal < 0.5] = 0 + + labelled_mask, num_labels = label(binary_mask) + # Let us now remove all the small regions, i.e., less than the specified mm. + print(binary_mask.shape) + print(labelled_mask.shape) + print(coronal.shape) + for get_label in range(num_labels): + refined_mask = np.zeros(binary_mask.shape) + refined_mask[labelled_mask == get_label] = 1 + x, y = np.nonzero(refined_mask) + x1, y1 = np.max(x), np.min(y) + x2, y2 = np.min(x), np.max(y) + refined_mask[:, y1-5:y2+5] = 1 + if (refined_mask * coronal).sum() == 0: + binary_mask[labelled_mask == get_label] = 0 + + predicted[0, ...] = np.expand_dims(binary_mask, axis=-1) + + return predicted + + +# Read .nii files using itk +if __name__ == '__main__': + print("remove wrong segmentation on sagittal views") diff --git a/src/LFBNet/utilities/compute_surrogate_features.py b/src/LFBNet/utilities/compute_surrogate_features.py index 2fd5670..9fe048c 100644 --- a/src/LFBNet/utilities/compute_surrogate_features.py +++ b/src/LFBNet/utilities/compute_surrogate_features.py @@ -119,7 +119,7 @@ def get_features(mask_to_compute_feature_on: ndarray = None): @staticmethod def compute_features_in_physical_space( - sagittal: dict = None, coronal: dic = None, voxel_size=None + sagittal: dict = None, coronal: dict = None, voxel_size=None ): """ Compute features in physical space. Basically multiply the given TMV or dissemination by the voxel space. diff --git a/read_3D_nifti_mask_image_compute_TMTV_Dmax_save_as_csv_file.py b/src/LFBNet/utilities/read_3D_nifti_mask_image_compute_TMTV_Dmax_save_as_csv_file.py similarity index 96% rename from read_3D_nifti_mask_image_compute_TMTV_Dmax_save_as_csv_file.py rename to src/LFBNet/utilities/read_3D_nifti_mask_image_compute_TMTV_Dmax_save_as_csv_file.py index 049ada3..6ce8ce1 100644 --- a/read_3D_nifti_mask_image_compute_TMTV_Dmax_save_as_csv_file.py +++ b/src/LFBNet/utilities/read_3D_nifti_mask_image_compute_TMTV_Dmax_save_as_csv_file.py @@ -1,117 +1,120 @@ -""" Script to compute metabolic tumor volume (TMTV) and dissemination from a given 3D mask images. -1. read .nii files -2. read the pixel spacing -3. compute the total metabolic tumor volume (TMTV) -4. Compute the lesion dissemination (Dmax) -5. calculate TMTV and Dmax in physical spacing, using the pixel spacing -6. save in .CSV as columns of: patient name (anonymous), pixel_spaing, TMTV, Dmax -""" - -# import important libraries -import os -import glob -from tqdm import tqdm -import argparse -from numpy.random import seed -seed(1) - -import nibabel as nib -import numpy as np -import csv - -# library for dmax computation -from skimage.measure import label, regionprops -from skimage import data, util -from scipy.spatial import distance - -import pathlib - - - -# function to write csv file -def write_to_csv_file(array, output_path, file_name="csv"): - """ - :param array: array that consists rows and columns to be saved to csv file - :param output_path: The directory to save csv file - :param file_name: Name of the csv file - :return: saved file_name.csv in the older output_path - """ - array = np.array(array) - file_name = str(output_path) + '/' + str(file_name) + '.csv' - - with open(file_name, 'w', newline='') as output: - output_data = csv.writer(output, delimiter=",") - for row in range(array.shape[0]): - output_data.writerow(array[row][:]) - - print("saved at: ", file_name) - - -# function to read .nii and compute biomarker values -def read_nii_mask_save_csv_tmtv_dmax(input_path, output_path): - """ - :param input_path: Path to the directory that consists the directory for .nii files - :param output_path: The directory to save csv file, after computing the TMTV and Dmax, pixel spacing - :return: read .nii, compute TMTV, Dmax - """ - case_ids = os.listdir(input_path) - print("Total number of cases to read: 0.1%d", len(case_ids)) - name_x_y_z_TMTV_dmax = [["ID", "X", "Y", "Z", 'TMTV', "Dmax"]] - - - for n, case_name in tqdm(enumerate(case_ids), total=len(case_ids)): - path_img_nii = str(input_path) + "/" + str(case_name) + "/" + str(case_name) + ".nii.gz" - try: - # Read .nii files - gt = nib.load(path_img_nii) - res_pet = gt.header.get_zooms() - gt = np.asanyarray(gt.dataobj) - - # Compute TMTV - def compute_TMTV(gt_): - gt_[gt_ > 0] = 1 - gt_[gt_ <= 0] = 0 - return np.sum(gt_ > 0) - - #compute dmax - def compute_dmax(gt_): - # label images - gt_ = util.img_as_ubyte(gt_) > 0 - gt_= label(gt_, connectivity=gt_.ndim) - props = regionprops(gt_) - - dist_all = [] - for k in range(len(props)): - # physical space - a = np.multiply(np.array(props[k].centroid), np.array(res_pet)) - - # compute the distance between all other centroids: - if len(props) >1: - for kk in range(len(props)): - b = np.multiply(np.array(props[kk].centroid), np.array(res_pet)) - dist = distance.euclidean(a, b) - dist_all.append(dist) - else: - dist_all.append(0) - return np.max(dist_all) - - tmtv = compute_TMTV(gt.copy()) - dmax = compute_dmax(gt.copy()) - name_x_y_z_TMTV_dmax.append([str(case_name), res_pet[0], res_pet[1], res_pet[2], - tmtv * res_pet[0] * res_pet[1] * res_pet[2], dmax]) - except: - print(f"Error reading {path_img_nii}") - continue - - write_to_csv_file(name_x_y_z_TMTV_dmax, output_path, file_name="vienna_data_xyz_tmtv_dmax") - print('Total number of patients correctly read and their volume calculated: ', len(name_x_y_z_TMTV_dmax) - 1) - print("Done !!") - - -if __name__ == "__main__": - # We assume the .nii file name and the folder name are the same - parser = argparse.ArgumentParser(description="script to read nii files and compute TMTV and Dmax") - parser.add_argument("--input_dir", dest='input_dir', type=pathlib.Path, help="Input directory path to .nii files") - parser.add_argument("--output_dir", dest='output_dir', type=pathlib.Path, help='output directory path') - args = parser.parse_args() - read_nii_mask_save_csv_tmtv_dmax(args.input_dir, args.output_dir) +""" Script to compute metabolic tumor volume (TMTV) and dissemination from a given 3D mask images. +1. read .nii files +2. read the pixel spacing +3. compute the total metabolic tumor volume (TMTV) +4. Compute the lesion dissemination (Dmax) +5. calculate TMTV and Dmax in physical spacing, using the pixel spacing +6. save in .CSV as columns of: patient name (anonymous), pixel_spaing, TMTV, Dmax +""" + +# import important libraries +import os +import glob +from tqdm import tqdm +import argparse +from numpy.random import seed +seed(1) + +import nibabel as nib +import numpy as np +import csv + +# library for dmax computation +from skimage.measure import label, regionprops +from skimage import data, util +from scipy.spatial import distance + +import pathlib + + + +# function to write csv file +def write_to_csv_file(array, output_path, file_name="csv"): + """ + :param array: array that consists rows and columns to be saved to csv file + :param output_path: The directory to save csv file + :param file_name: Name of the csv file + :return: saved file_name.csv in the older output_path + """ + array = np.array(array) + file_name = str(output_path) + '/' + str(file_name) + '.csv' + + with open(file_name, 'w', newline='') as output: + output_data = csv.writer(output, delimiter=",") + for row in range(array.shape[0]): + output_data.writerow(array[row][:]) + + print("saved at: ", file_name) + + +# function to read .nii and compute biomarker values +def read_nii_mask_save_csv_tmtv_dmax(input_path, output_path): + """ + :param input_path: Path to the directory that consists the directory for .nii files + :param output_path: The directory to save csv file, after computing the TMTV and Dmax, pixel spacing + :return: read .nii, compute TMTV, Dmax + """ + case_ids = os.listdir(input_path) + print("Total number of cases to read: 0.1%d", len(case_ids)) + name_x_y_z_TMTV_dmax = [["ID", "X", "Y", "Z", 'TMTV', "Dmax"]] + + + for n, case_name in tqdm(enumerate(case_ids), total=len(case_ids)): + path_img_nii = str(input_path) + "/" + str(case_name) + "/gt/" + path_img_nii = glob.glob(path_img_nii + "/*.nii.gz")[0] + try: + # Read .nii files + gt = nib.load(path_img_nii) + res_pet = gt.header.get_zooms() + gt = np.asanyarray(gt.dataobj) + gt[gt > 0] = 1 + gt[gt <= 0] = 0 + + # Compute TMTV + def compute_TMTV(gt_): + gt_[gt_ > 0] = 1 + gt_[gt_ <= 0] = 0 + return np.sum(gt_ > 0) + + #compute dmax + def compute_dmax(gt_): + # label images + gt_ = util.img_as_ubyte(gt_) > 0 + gt_= label(gt_, connectivity=gt_.ndim) + props = regionprops(gt_) + + dist_all = [] + for k in range(len(props)): + # physical space + a = np.multiply(np.array(props[k].centroid), np.array(res_pet)) + + # compute the distance between all other centroids: + if len(props) >1: + for kk in range(len(props)): + b = np.multiply(np.array(props[kk].centroid), np.array(res_pet)) + dist = distance.euclidean(a, b) + dist_all.append(dist) + else: + dist_all.append(0) + return np.max(dist_all) + + tmtv = compute_TMTV(gt.copy()) + dmax = compute_dmax(gt.copy()) + name_x_y_z_TMTV_dmax.append([str(case_name), res_pet[0], res_pet[1], res_pet[2], + tmtv * res_pet[0] * res_pet[1] * res_pet[2], dmax]) + except: + print(f"Error reading {path_img_nii}") + continue + + write_to_csv_file(name_x_y_z_TMTV_dmax, output_path, file_name="data_xyz_tmtv_dmax") + print('Total number of patients correctly read and their volume calculated: ', len(name_x_y_z_TMTV_dmax) - 1) + print("Done !!") + + +if __name__ == "__main__": + # We assume the .nii file name and the folder name are the same + parser = argparse.ArgumentParser(description="script to read nii files and compute TMTV and Dmax") + parser.add_argument("--input_dir", dest='input_dir', type=pathlib.Path, help="Input directory path to .nii files") + parser.add_argument("--output_dir", dest='output_dir', type=pathlib.Path, help='output directory path') + args = parser.parse_args() + read_nii_mask_save_csv_tmtv_dmax(args.input_dir, args.output_dir) diff --git a/src/__pycache__/__init__.cpython-36.pyc b/src/__pycache__/__init__.cpython-36.pyc deleted file mode 100644 index af3f7d9..0000000 Binary files a/src/__pycache__/__init__.cpython-36.pyc and /dev/null differ diff --git a/src/run/__pycache__/__init__.cpython-36.pyc b/src/run/__pycache__/__init__.cpython-36.pyc deleted file mode 100644 index bdd7def..0000000 Binary files a/src/run/__pycache__/__init__.cpython-36.pyc and /dev/null differ diff --git a/src/run/__pycache__/parse_argument.cpython-36.pyc b/src/run/__pycache__/parse_argument.cpython-36.pyc deleted file mode 100644 index 9e7a2c5..0000000 Binary files a/src/run/__pycache__/parse_argument.cpython-36.pyc and /dev/null differ diff --git a/src/run/__pycache__/trainer.cpython-36.pyc b/src/run/__pycache__/trainer.cpython-36.pyc deleted file mode 100644 index edc9ee5..0000000 Binary files a/src/run/__pycache__/trainer.cpython-36.pyc and /dev/null differ diff --git a/src/run/trainer.py b/src/run/trainer.py index a926015..04a7d2e 100644 --- a/src/run/trainer.py +++ b/src/run/trainer.py @@ -51,7 +51,7 @@ from src.LFBNet.losses import losses from src.LFBNet.preprocessing import save_nii_images from src.LFBNet.utilities import train_valid_paths - +from src.LFBNet.postprocessing import remove_outliers_in_sagittal # choose cuda gpu CUDA_VISIBLE_DEVICES = 1 os.environ["CUDA_VISIBLE_DEVICES"] = "1" @@ -804,6 +804,7 @@ def evaluation_test( 'sensitivity': binary.sensitivity(predicted, ground_truth)} predicted = self.model.combine_and_train.predict([input_image, feedback_latent]) + predicted = remove_outliers_in_sagittal(predicted) save_nii_images( [predicted, ground_truth, input_image], identifier=str(case_name), name=[case_name + "_predicted", case_name + "_ground_truth", case_name + "_pet"], @@ -831,6 +832,7 @@ def prediction(self, input_image: ndarray = None, case_name: str = None): feedback_latent = self.model.feedback_latent.predict(predicted) # feedback: hf predicted = self.model.combine_and_train.predict([input_image, feedback_latent]) + predicted = remove_outliers_in_sagittal(predicted) save_nii_images( image=[predicted, input_image], identifier=str(case_name), name=[case_name + "_predicted", case_name + "_pet"],