diff --git a/scripts/1.U-net_chris.py b/scripts/1.U-net_chris.py index 0243dbe..de8ab83 100755 --- a/scripts/1.U-net_chris.py +++ b/scripts/1.U-net_chris.py @@ -13,6 +13,7 @@ from sfransen.utils_quintin import * from sfransen.utils_quintin import * from sfransen.DWI_exp import IntermediateImages, dice_coef from sfransen.DWI_exp.preprocessing_function import preprocess +from sfransen.DWI_exp.losses import weighted_binary_cross_entropy import yaml import numpy as np from tqdm import tqdm @@ -69,11 +70,18 @@ MODEL_SELECTION_METRIC = 'val_loss' MODEL_SELECTION_DIRECTION = "min" # Change to 'max' if higher value is better EARLY_STOPPING_METRIC = 'val_loss' EARLY_STOPPING_DIRECTION = 'min' +# MODEL_SELECTION_METRIC = 'weighted_binary_cross_entropy' +# MODEL_SELECTION_DIRECTION = "min" # Change to 'max' if higher value is better +# EARLY_STOPPING_METRIC = 'weighted_binary_cross_entropy' +# EARLY_STOPPING_DIRECTION = 'min' # Training configuration # add metric ROC_AUC TRAINING_METRICS = ["binary_crossentropy", "binary_accuracy", dice_coef] -loss = BinaryFocalLoss(gamma=FOCAL_LOSS_GAMMA) +# loss = BinaryFocalLoss(gamma=FOCAL_LOSS_GAMMA) +weight_for_0 = 0.05 +weight_for_1 = 0.95 +loss = weighted_binary_cross_entropy({0: weight_for_0, 1: weight_for_1}) optimizer = Adam(learning_rate=INITIAL_LEARNING_RATE) # Create folder structure in the output directory @@ -115,7 +123,7 @@ with open(path.join(DATA_DIR, f"seg.txt"), 'r') as f: num_images = len(seg_paths) # Read and preprocess each of the paths for each series, and the segmentations. -for img_idx in tqdm(range(num_images)): # [:40]):for less images +for img_idx in tqdm(range(num_images)): #[:40]): #for less images img_s = {s: sitk.ReadImage(image_paths[s][img_idx], sitk.sitkFloat32) for s in args.series} seg_s = sitk.ReadImage(seg_paths[img_idx], sitk.sitkFloat32) @@ -166,11 +174,6 @@ print_(f"The shape of valid_data input = {np.shape(valid_data[0])}") print_(f"The shape of valid_data label = {np.shape(valid_data[1])}") callbacks = [ - EarlyStopping( - monitor=EARLY_STOPPING_METRIC, - mode=EARLY_STOPPING_DIRECTION, - patience=EARLY_STOPPING, - verbose=2), ModelCheckpoint( filepath=path.join(PROJECT_DIR, "models", JOB_NAME + ".h5"), monitor=MODEL_SELECTION_METRIC, @@ -185,6 +188,11 @@ callbacks = [ save_best_only=True), CSVLogger( filename=path.join(PROJECT_DIR, "logs", f"{JOB_NAME}.csv")), + EarlyStopping( + monitor=EARLY_STOPPING_METRIC, + mode=EARLY_STOPPING_DIRECTION, + patience=EARLY_STOPPING, + verbose=2), IntermediateImages( validation_set=valid_data, sequences=args.series, diff --git a/scripts/4.frocs.py b/scripts/4.frocs.py index c1053db..a5d7afe 100755 --- a/scripts/4.frocs.py +++ b/scripts/4.frocs.py @@ -6,6 +6,8 @@ import numpy as np import multiprocessing from functools import partial import os +from os import path + from sfransen.utils_quintin import * from sfransen.DWI_exp.helpers import * @@ -14,6 +16,7 @@ from sfransen.DWI_exp.callbacks import dice_coef #from sfransen.FROC.blob_preprocess import * from sfransen.FROC.cal_froc_from_np import * from sfransen.load_images import load_images_parrallel +from sfransen.DWI_exp.losses import weighted_binary_cross_entropy parser = argparse.ArgumentParser( @@ -36,7 +39,7 @@ SERIES = args.series series_ = '_'.join(args.series) EXPERIMENT = args.experiment -MODEL_PATH = f'./../train_output/{EXPERIMENT}_{series_}/models/{EXPERIMENT}_{series_}_dice.h5' +MODEL_PATH = f'./../train_output/{EXPERIMENT}_{series_}/models/{EXPERIMENT}_{series_}.h5' YAML_DIR = f'./../train_output/{EXPERIMENT}_{series_}' IMAGE_DIR = f'./../train_output/{EXPERIMENT}_{series_}' @@ -54,8 +57,46 @@ N_CPUS = 12 ########## load images in parrallel ############## print_(f"> Loading images into RAM...") -# read paths from txt +# # read paths from txt +# image_paths = {} +# for s in SERIES: +# with open(path.join(DATA_DIR, f"{s}.txt"), 'r') as f: +# image_paths[s] = [l.strip() for l in f.readlines()] +# with open(path.join(DATA_DIR, f"seg.txt"), 'r') as f: +# seg_paths = [l.strip() for l in f.readlines()] +# num_images = len(seg_paths) + +# # create pool of workers +# pool = multiprocessing.Pool(processes=N_CPUS) +# partial_images = partial(load_images_parrallel, +# seq = 'images', +# target_shape=IMAGE_SHAPE, +# target_space = TARGET_SPACING) +# partial_seg = partial(load_images_parrallel, +# seq = 'seg', +# target_shape=IMAGE_SHAPE, +# target_space = TARGET_SPACING) + +# #load images +# images = [] +# for s in SERIES: +# image_paths_seq = image_paths[s] +# image_paths_index = np.asarray(image_paths_seq)[TEST_INDEX][:5] +# data_list = pool.map(partial_images,image_paths_index) +# data = np.stack(data_list, axis=0) +# images.append(data) + +# images_list = np.transpose(images, (1, 2, 3, 4, 0)) +# print('>>>>> size image_list nmr 1:', np.shape(images_list)) + +# #load segmentations +# seg_paths_index = np.asarray(seg_paths)[TEST_INDEX][:5] +# data_list = pool.map(partial_seg,seg_paths_index) +# segmentations = np.stack(data_list, axis=0) + +########## test with old method ############# image_paths = {} + for s in SERIES: with open(path.join(DATA_DIR, f"{s}.txt"), 'r') as f: image_paths[s] = [l.strip() for l in f.readlines()] @@ -63,38 +104,33 @@ with open(path.join(DATA_DIR, f"seg.txt"), 'r') as f: seg_paths = [l.strip() for l in f.readlines()] num_images = len(seg_paths) -# create pool of workers -pool = multiprocessing.Pool(processes=N_CPUS) -partial_images = partial(load_images_parrallel, - seq = 'images', - target_shape=IMAGE_SHAPE, - target_space = TARGET_SPACING) -partial_seg = partial(load_images_parrallel, - seq = 'seg', - target_shape=IMAGE_SHAPE, - target_space = TARGET_SPACING) - -#load images images = [] -for s in SERIES: - image_paths_seq = image_paths[s] - image_paths_index = np.asarray(image_paths_seq)[TEST_INDEX] - data_list = pool.map(partial_images,image_paths_index) - data = np.stack(data_list, axis=0) - images.append(data) +images_list = [] +segmentations = [] -images_list = np.transpose(images, (1, 2, 3, 4, 0)) +# Read and preprocess each of the paths for each series, and the segmentations. +for img_idx in tqdm(range(len(TEST_INDEX))): #[:40]): #for less images + # print('images number',[TEST_INDEX[img_idx]]) + img_s = {f'{s}': sitk.ReadImage(image_paths[s][TEST_INDEX[img_idx]], sitk.sitkFloat32) for s in SERIES} + seg_s = sitk.ReadImage(seg_paths[TEST_INDEX[img_idx]], sitk.sitkFloat32) + img_n, seg_n = preprocess(img_s, seg_s, + shape=IMAGE_SHAPE, spacing=TARGET_SPACING) + for seq in img_n: + images.append(img_n[f'{seq}']) + images_list.append(images) + images = [] + segmentations.append(seg_n) -#load segmentations -seg_paths_index = np.asarray(seg_paths)[TEST_INDEX] -data_list = pool.map(partial_seg,seg_paths_index) -segmentations = np.stack(data_list, axis=0) +images_list = np.transpose(images_list, (0, 2, 3, 4, 1)) + +print('>>>>> size image_list nmr 2:', np.shape(images_list), '. equal to: (5, 192, 192, 24, 3)?') ########### load module ################## print(' >>>>>>> LOAD MODEL <<<<<<<<<') dependencies = { - 'dice_coef': dice_coef + 'dice_coef': dice_coef, + 'weighted_cross_entropy_fn':weighted_binary_cross_entropy } reconstructed_model = load_model(MODEL_PATH, custom_objects=dependencies) # reconstructed_model.summary(line_length=120) @@ -127,16 +163,51 @@ predictions = zeros # perform Froc metrics = evaluate(y_true=segmentations, y_pred=predictions) -dump_dict_to_yaml(metrics, YAML_DIR, "froc_metrics", verbose=True) +dump_dict_to_yaml(metrics, YAML_DIR, "froc_metrics_focal_10", verbose=True) ############## save image as example ################# +# save image nmr 2 +img_s = sitk.GetImageFromArray(predictions_blur[2].squeeze()) +sitk.WriteImage(img_s, f"{IMAGE_DIR}/predictions_blur_002_old.nii.gz") + +img_s = sitk.GetImageFromArray(predictions[2].squeeze()) +sitk.WriteImage(img_s, f"{IMAGE_DIR}/predictions_002_old.nii.gz") + +img_s = sitk.GetImageFromArray(segmentations[2].squeeze()) +sitk.WriteImage(img_s, f"{IMAGE_DIR}/segmentations_002_old.nii.gz") + +img_s = sitk.GetImageFromArray(np.transpose(images_list[2,:,:,:,0].squeeze())) +sitk.WriteImage(img_s, f"{IMAGE_DIR}/t2_002_old.nii.gz") + # save image nmr 3 img_s = sitk.GetImageFromArray(predictions_blur[3].squeeze()) -sitk.WriteImage(img_s, f"{IMAGE_DIR}/predictions_blur_001.nii.gz") +sitk.WriteImage(img_s, f"{IMAGE_DIR}/predictions_blur_003_old.nii.gz") img_s = sitk.GetImageFromArray(predictions[3].squeeze()) -sitk.WriteImage(img_s, f"{IMAGE_DIR}/predictions_001.nii.gz") +sitk.WriteImage(img_s, f"{IMAGE_DIR}/predictions_003_old.nii.gz") img_s = sitk.GetImageFromArray(segmentations[3].squeeze()) -sitk.WriteImage(img_s, f"{IMAGE_DIR}/segmentations_001.nii.gz") \ No newline at end of file +sitk.WriteImage(img_s, f"{IMAGE_DIR}/segmentations_003_old.nii.gz") + +img_s = sitk.GetImageFromArray(np.transpose(images_list[3,:,:,:,0].squeeze())) +sitk.WriteImage(img_s, f"{IMAGE_DIR}/t2_003_old.nii.gz") + +img_s = sitk.GetImageFromArray(np.transpose(images_list[3,:,:,:,1].squeeze())) +sitk.WriteImage(img_s, f"{IMAGE_DIR}/highb_003_old.nii.gz") + +# save image nmr 3 +img_s = sitk.GetImageFromArray(predictions_blur[4].squeeze()) +sitk.WriteImage(img_s, f"{IMAGE_DIR}/predictions_blur_004_old.nii.gz") + +img_s = sitk.GetImageFromArray(predictions[4].squeeze()) +sitk.WriteImage(img_s, f"{IMAGE_DIR}/predictions_004_old.nii.gz") + +img_s = sitk.GetImageFromArray(segmentations[4].squeeze()) +sitk.WriteImage(img_s, f"{IMAGE_DIR}/segmentations_004_old.nii.gz") + +img_s = sitk.GetImageFromArray(np.transpose(images_list[4,:,:,:,0].squeeze())) +sitk.WriteImage(img_s, f"{IMAGE_DIR}/t2_004_old.nii.gz") + +img_s = sitk.GetImageFromArray(np.transpose(images_list[4,:,:,:,1].squeeze())) +sitk.WriteImage(img_s, f"{IMAGE_DIR}/highb_004_old.nii.gz") \ No newline at end of file diff --git a/scripts/5.Visualize_frocs.py b/scripts/5.Visualize_frocs.py index a48b2c2..1c70233 100755 --- a/scripts/5.Visualize_frocs.py +++ b/scripts/5.Visualize_frocs.py @@ -28,7 +28,7 @@ experiment_path = [] experiment_metrics = {} auroc = [] for idx in range(len(args.experiment)): - experiment_path = f'./../train_output/{experiments[idx]}/froc_metrics.yml' + experiment_path = f'./../train_output/{experiments[idx]}/froc_metrics_focal_10.yml' experiment_metrics = read_yaml_to_dict(experiment_path) auroc.append(round(experiment_metrics['auroc'],3)) diff --git a/scripts/6.saliency_map.py b/scripts/6.saliency_map.py index 40265be..a7730bc 100755 --- a/scripts/6.saliency_map.py +++ b/scripts/6.saliency_map.py @@ -9,6 +9,8 @@ from sfransen.utils_quintin import * from sfransen.DWI_exp import preprocess from sfransen.DWI_exp.helpers import * from sfransen.DWI_exp.callbacks import dice_coef +from sfransen.DWI_exp.losses import weighted_binary_cross_entropy + from sfransen.FROC.blob_preprocess import * from sfransen.FROC.cal_froc_from_np import * from sfransen.load_images import load_images_parrallel @@ -40,13 +42,13 @@ TARGET_SPACING = (0.5, 0.5, 3) INPUT_SHAPE = (192, 192, 24, len(SERIES)) IMAGE_SHAPE = INPUT_SHAPE[:3] -froc_metrics = read_yaml_to_dict(f'{YAML_DIR}/froc_metrics.yml') -top_10_idx = np.argsort(froc_metrics['roc_pred'])[-1 :] +# froc_metrics = read_yaml_to_dict(f'{YAML_DIR}/froc_metrics.yml') +# top_10_idx = np.argsort(froc_metrics['roc_pred'])[-1 :] DATA_SPLIT_INDEX = read_yaml_to_dict('./../data/Nijmegen paths/train_val_test_idxs.yml') TEST_INDEX = DATA_SPLIT_INDEX['val_set0'] -TEST_INDEX_top10 = [TEST_INDEX[i] for i in top_10_idx] +# TEST_INDEX_top10 = [TEST_INDEX[i] for i in top_10_idx] TEST_INDEX_image = [371] N_CPUS = 12 @@ -93,7 +95,8 @@ segmentations = np.stack(data_list, axis=0) print(' >>>>>>> LOAD MODEL <<<<<<<<<') dependencies = { - 'dice_coef': dice_coef + 'dice_coef': dice_coef, + 'weighted_cross_entropy_fn':weighted_binary_cross_entropy } reconstructed_model = load_model(MODEL_PATH, custom_objects=dependencies) # reconstructed_model.layers[-1].activation = tf.keras.activations.linear diff --git a/scripts/8.Visualize_training.py b/scripts/8.Visualize_training.py index 69f9aa5..5ce8d87 100755 --- a/scripts/8.Visualize_training.py +++ b/scripts/8.Visualize_training.py @@ -29,9 +29,9 @@ df = pd.read_csv(f'{folder_input}') # read csv file for metric in df: if not metric == 'epoch': - # if metric == 'loss' or metric == 'val_loss': - plt.plot(df['epoch'], df[metric], label=metric) - # plt.ylim(ymin=0,ymax=0.01) + if metric == 'loss' or metric == 'val_loss': + plt.plot(df['epoch'], df[metric], label=metric) + plt.ylim(ymin=0,ymax=0.01) diff --git a/scripts/calc_adc.py b/scripts/calc_adc.py index 7017274..97244e9 100755 --- a/scripts/calc_adc.py +++ b/scripts/calc_adc.py @@ -2,34 +2,34 @@ import numpy as np import SimpleITK as sitk import matplotlib.pyplot as plt ######## load images ############# -# path_b50 = '/data/pca-rad/datasets/radboud_new/pat0351/2016/diffusie_cro/b-50/nifti_image.nii.gz' -# path_b400 = '/data/pca-rad/datasets/radboud_new/pat0351/2016/diffusie_cro/b-400/nifti_image.nii.gz' -# path_b800 = '/data/pca-rad/datasets/radboud_new/pat0351/2016/diffusie_cro/b-800/nifti_image.nii.gz' -# path_b1400 = '/data/pca-rad/datasets/radboud_new/pat0351/2016/diffusie_cro/b-1400/nifti_image.nii.gz' -# path_adc = '/data/pca-rad/datasets/radboud_new/pat0351/2016/dADC/nifti_image.nii.gz' +path_b50 = '/data/pca-rad/datasets/radboud_new/pat0634/2016/diffusie_cro/b-50/nifti_image.nii.gz' +path_b400 = '/data/pca-rad/datasets/radboud_new/pat0634/2016/diffusie_cro/b-400/nifti_image.nii.gz' +path_b800 = '/data/pca-rad/datasets/radboud_new/pat0634/2016/diffusie_cro/b-800/nifti_image.nii.gz' +path_b1400 = '/data/pca-rad/datasets/radboud_new/pat0634/2016/diffusie_cro/b-1400/nifti_image.nii.gz' +path_adc = '/data/pca-rad/datasets/radboud_new/pat0634/2016/dADC/nifti_image.nii.gz' -# path_b50 = 'X:/sfransen/train_output/adc_exp/b50.nii.gz' -# path_b400 = 'X:/sfransen/train_output/adc_exp/b400.nii.gz' -# path_b800 = 'X:/sfransen/train_output/adc_exp/b800.nii.gz' -# path_b1400 = 'X:/sfransen/train_output/adc_exp/b1400.nii.gz' -# path_adc = 'X:/sfransen/train_output/adc_exp/adc.nii.gz' +# path_b50 = 'X:/sfransen/train_output/adc_exp/b50_true.nii.gz' +# path_b400 = 'X:/sfransen/train_output/adc_exp/b400_true.nii.gz' +# path_b800 = 'X:/sfransen/train_output/adc_exp/b800_true.nii.gz' +# path_b1400 = 'X:/sfransen/train_output/adc_exp/b1400_true.nii.gz' +# path_adc = 'X:/sfransen/train_output/adc_exp/adc_true.nii.gz' -path_b50 = '/data/pca-rad/sfransen/train_output/adc_exp/b50_true.nii.gz' -path_b400 = '/data/pca-rad/sfransen/train_output/adc_exp/b400_true.nii.gz' -path_b800 = '/data/pca-rad/sfransen/train_output/adc_exp/b800_true.nii.gz' -path_b1400 = '/data/pca-rad/sfransen/train_output/adc_exp/b1400_true.nii.gz' -path_adc = '/data/pca-rad/sfransen/train_output/adc_exp/adc_calc_b50_b400_b800.nii.gz' +# path_b50 = '/data/pca-rad/sfransen/train_output/adc_exp/b50_true.nii.gz' +# path_b400 = '/data/pca-rad/sfransen/train_output/adc_exp/b400_true.nii.gz' +# path_b800 = '/data/pca-rad/sfransen/train_output/adc_exp/b800_true.nii.gz' +# path_b1400 = '/data/pca-rad/sfransen/train_output/adc_exp/b1400_true.nii.gz' +# path_adc = '/data/pca-rad/sfransen/train_output/adc_exp/adc_true.nii.gz' -b50 = sitk.ReadImage(path_b50, sitk.sitkFloat32) -b50 = sitk.GetArrayFromImage(b50) -b400 = sitk.ReadImage(path_b400, sitk.sitkFloat32) -b400 = sitk.GetArrayFromImage(b400) -b800 = sitk.ReadImage(path_b800, sitk.sitkFloat32) -b800 = sitk.GetArrayFromImage(b800) -b1400 = sitk.ReadImage(path_b1400, sitk.sitkFloat32) -b1400 = sitk.GetArrayFromImage(b1400) -adc = sitk.ReadImage(path_adc, sitk.sitkFloat32) -adc = sitk.GetArrayFromImage(adc) +b50_img = sitk.ReadImage(path_b50, sitk.sitkFloat32) +b50 = sitk.GetArrayFromImage(b50_img) +b400_img = sitk.ReadImage(path_b400, sitk.sitkFloat32) +b400 = sitk.GetArrayFromImage(b400_img) +b800_img = sitk.ReadImage(path_b800, sitk.sitkFloat32) +b800 = sitk.GetArrayFromImage(b800_img) +b1400_img = sitk.ReadImage(path_b1400, sitk.sitkFloat32) +b1400_original = sitk.GetArrayFromImage(b1400_img) +adc_img = sitk.ReadImage(path_adc, sitk.sitkFloat32) +adc_original = sitk.GetArrayFromImage(adc_img) def show_img(greyscale_img): fig = plt.figure() @@ -44,8 +44,9 @@ def calc_adc(b50, b400, b800): denominator = np.multiply((50 - mean_dwi), np.subtract(np.log(b50), mean_si)) + np.multiply((400 - mean_dwi), np.subtract(np.log(b400), mean_si)) + np.multiply((800 - mean_dwi), np.subtract(np.log(b800), mean_si)) numerator = np.power((50 - mean_dwi), 2) + np.power((400 - mean_dwi), 2) + np.power((800 - mean_dwi), 2) - adc = np.divide(denominator, numerator) - return adc * -1000000 + adc_with_zeros = np.divide(denominator, numerator) * -1000000 + adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros)) + return adc def calc_adc_1(b50,b800): mean_dwi = (50 + 800) / 2 @@ -53,8 +54,9 @@ def calc_adc_1(b50,b800): denominator = np.multiply((50 - mean_dwi), np.subtract(np.log(b50), mean_si)) + np.multiply((800 - mean_dwi), np.subtract(np.log(b800), mean_si)) numerator = np.power((50 - mean_dwi), 2) + np.power((800 - mean_dwi), 2) - adc = np.divide(denominator, numerator) - return adc * -1000000 + adc_with_zeros = np.divide(denominator, numerator) * -1000000 + adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros)) + return adc def calc_adc_2(b50,b400): mean_dwi = (50 + 400) / 2 @@ -62,8 +64,9 @@ def calc_adc_2(b50,b400): denominator = np.multiply((50 - mean_dwi), np.subtract(np.log(b50), mean_si)) + np.multiply((400 - mean_dwi), np.subtract(np.log(b400), mean_si)) numerator = np.power((50 - mean_dwi), 2) + np.power((400 - mean_dwi), 2) - adc = np.divide(denominator, numerator) - return adc * -1000000 + adc_with_zeros = np.divide(denominator, numerator) * -1000000 + adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros)) + return adc def calc_adc_3(b400,b800): mean_dwi = (400 + 800) / 2 @@ -71,29 +74,121 @@ def calc_adc_3(b400,b800): denominator = np.multiply((400 - mean_dwi), np.subtract(np.log(b400), mean_si)) + np.multiply((800 - mean_dwi), np.subtract(np.log(b800), mean_si)) numerator = np.power((400 - mean_dwi), 2) + np.power((800 - mean_dwi), 2) - adc = np.divide(denominator, numerator) - return adc * -1000000 + adc_with_zeros = np.divide(denominator, numerator) * -1000000 + adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros)) + return adc def calc_high_b(b_value_high,b_value,b_image,ADC_map): - high_b = np.multiply(b_image, np.log(np.multiply(np.subtract(b_value,b_value_high), ADC_map))) + high_b = np.multiply(b_image, np.exp(np.multiply(np.subtract(b_value,b_value_high), (np.divide(ADC_map,1000000))))) return high_b +adc_50_400_800 = calc_adc(b50,b400,b800) +adc_50_800 = calc_adc_1(b50,b800) +adc_50_400 = calc_adc_2(b50,b400) +adc_400_800 = calc_adc_3(b400,b800) -high_b = sitk.GetImageFromArray(b1400) -sitk.WriteImage(high_b, f"./../train_output/adc_exp/b1400_true.nii.gz") +high_b_1400_50 = calc_high_b(1400,50,b50,adc_50_800) +high_b_1400_all = calc_high_b(1400,50,b50,adc_50_400_800) -high_b = calc_high_b(1400,50,b50,adc) -high_b = sitk.GetImageFromArray(high_b) -sitk.WriteImage(high_b, f"./../train_output/adc_exp/b1400_ref_b50.nii.gz") +high_b_1400_400 = calc_high_b(1400,400,b400,adc_50_800) +high_b_1400_800 = calc_high_b(1400,800,b800,adc_50_800) -high_b = calc_high_b(1400,400,b400,adc) -high_b = sitk.GetImageFromArray(high_b) -sitk.WriteImage(high_b, f"./../train_output/adc_exp/b1400_ref_b400.nii.gz") +fig, (ax1, ax2, ax3) = plt.subplots(1, 3) +fig.suptitle('ADC calculated with b50 b400 b800') +ax1.imshow(adc_50_400_800[:][:][13],cmap='gray') +ax1.set_title('calculated b50 b400 b800') +ax2.imshow(adc_original[:][:][13],cmap='gray') +ax2.set_title('original') +error_map = np.subtract(adc_original[:][:][13],adc_50_800[:][:][13]) +ax3.imshow(error_map,cmap='gray') +ax3.set_title('error map') +path = f"ADC_634.png" +fig.savefig(path, dpi=300, bbox_inches='tight') -high_b = calc_high_b(1400,800,b800,adc) -high_b = sitk.GetImageFromArray(high_b) -sitk.WriteImage(high_b, f"./../train_output/adc_exp/b1400_ref_b800.nii.gz") +fig, (ax1, ax2, ax3) = plt.subplots(1, 3) +fig.suptitle('ADC calculated with b50 b800') +ax1.imshow(adc_50_800[:][:][13],cmap='gray') +ax1.set_title('calculated b50 b800') +ax2.imshow(adc_original[:][:][13],cmap='gray') +ax2.set_title('original') +error_map = np.subtract(adc_original[:][:][13],adc_50_800[:][:][13]) +ax3.imshow(error_map,cmap='gray') +ax3.set_title('error map') +path = f"ADC_634_1.png" +fig.savefig(path, dpi=300, bbox_inches='tight') +fig, (ax1, ax2, ax3) = plt.subplots(1, 3) +fig.suptitle('Difference between ADC calculation') +ax1.imshow(adc_50_400_800[:][:][13],cmap='gray') +ax1.set_title('calculated b50 b400 b800') +ax2.imshow(adc_50_800[:][:][13],cmap='gray') +ax2.set_title('calculated b50 b800') +error_map = np.subtract(adc_50_800[:][:][13],adc_50_400_800[:][:][13]) +ax3.imshow(error_map,cmap='gray') +ax3.set_title('error map') +path = f"ADC_634_different_calculation_1.png" +fig.savefig(path, dpi=300, bbox_inches='tight') + +fig, (ax1, ax2, ax3) = plt.subplots(1, 3) +fig.suptitle('Difference between ADC calculation') +ax1.imshow(adc_50_800[:][:][13],cmap='gray') +ax1.set_title('calculated b50 b800') +ax2.imshow(adc_50_400[:][:][13],cmap='gray') +ax2.set_title('calculated b50 b400') +error_map = np.subtract(adc_50_800[:][:][13],adc_50_400[:][:][13]) +ax3.imshow(error_map,cmap='gray') +ax3.set_title('error map') +path = f"ADC_634_different_calculation_2.png" +fig.savefig(path, dpi=300, bbox_inches='tight') + +fig, (ax1, ax2, ax3) = plt.subplots(1, 3) +fig.suptitle('High B calculated with b50 reference and ADC from b50&b800') +ax1.imshow(high_b_1400_50[:][:][13],cmap='gray') +ax1.set_title('calculated') +ax2.imshow(b1400_original[:][:][13],cmap='gray') +ax2.set_title('original') +error_map = np.subtract(b1400_original[:][:][13],high_b_1400_50[:][:][13]) +ax3.imshow(error_map,cmap='gray') +ax3.set_title('error map') +path = f"HighB_b50_b800_634_1.png" +fig.savefig(path, dpi=300, bbox_inches='tight') + +fig, (ax1, ax2, ax3) = plt.subplots(1, 3) +fig.suptitle('High B calculated with b50 reference and ADC from b50&b400&b800') +ax1.imshow(high_b_1400_50[:][:][13],cmap='gray') +ax1.set_title('calculated') +ax2.imshow(b1400_original[:][:][13],cmap='gray') +ax2.set_title('original') +error_map = np.subtract(b1400_original[:][:][13],high_b_1400_50[:][:][13]) +ax3.imshow(error_map,cmap='gray') +ax3.set_title('error map') +path = f"HighB_b50_b800_634_2.png" +fig.savefig(path, dpi=300, bbox_inches='tight') + +# adc_50_400_800 = sitk.GetImageFromArray(adc_50_400_800) +# adc_50_400_800.CopyInformation(b50_img) +# sitk.WriteImage(adc_50_400_800, f"../train_output/adc_exp/adc_copied_with_b50.nii.gz") +# adc_50_800 = sitk.GetImageFromArray(adc_50_800) +# adc_50_800.CopyInformation(adc_img) +# sitk.WriteImage(adc_50_800, f"../train_output/adc_exp/adc_copied_with_adc.nii.gz") +# adc_50_400 = sitk.GetImageFromArray(adc_50_400) +# sitk.WriteImage(adc_50_400, f"../train_output/adc_exp/output/adc_calc_b50_b400.nii.gz") +# adc_400_800 = sitk.GetImageFromArray(adc_400_800) +# sitk.WriteImage(adc_400_800, f"../train_output/adc_exp/output/adc_calc_b400_b800.nii.gz") + +# high_b_1400_50 = sitk.GetImageFromArray(high_b_1400_50) +# high_b_1400_50.CopyInformation(b50_img) +# sitk.WriteImage(high_b_1400_50, f"../train_output/adc_exp/hb_copied_with_b50.nii.gz") +# high_b_1400_50_1 = sitk.GetImageFromArray(high_b_1400_50) +# high_b_1400_50_1.CopyInformation(adc_img) +# sitk.WriteImage(high_b_1400_50_1, f"../train_output/adc_exp/hb_copied_with_adc.nii.gz") + +# high_b_1400_400 = sitk.GetImageFromArray(high_b_1400_400) +# sitk.WriteImage(high_b_1400_400, f"../train_output/adc_exp/output/b1400_ref_b400.nii.gz") +# high_b_1400_800 = sitk.GetImageFromArray(high_b_1400_800) +# sitk.WriteImage(high_b_1400_800, f"../train_output/adc_exp/output/b1400_ref_b800.nii.gz") + +########################################################################################## # b50 = sitk.GetImageFromArray(b50) # sitk.WriteImage(b50, "./../train_output/adc_exp/b50_true.nii.gz") @@ -106,27 +201,5 @@ sitk.WriteImage(high_b, f"./../train_output/adc_exp/b1400_ref_b800.nii.gz") # b1400 = sitk.GetImageFromArray(b1400) # sitk.WriteImage(b1400,f"./../train_output/adc_exp/b1400_true.nii.gz") - # adc = sitk.GetImageFromArray(adc) -# sitk.WriteImage(adc, f"adc_true.nii.gz") - -# adc = calc_adc(b50,b400,b800) -# print("calculated with 3 adc shape:",adc.shape) -# adc = sitk.GetImageFromArray(adc) -# sitk.WriteImage(adc, f"adc_calc_b50_b400_b800.nii.gz") - -# adc = calc_adc_1(b50,b800) -# print("calculated with 2 adc shape:",adc.shape) -# adc = sitk.GetImageFromArray(adc) -# sitk.WriteImage(adc, f"adc_calc_b50_b800.nii.gz") - -# adc = calc_adc_2(b50,b400) -# print("calculated with 2 adc shape:",adc.shape) -# adc = sitk.GetImageFromArray(adc) -# sitk.WriteImage(adc, f"adc_calc_b50_b400.nii.gz") - -# adc = calc_adc_3(b400,b800) -# print("calculated with 2 adc shape:",adc.shape) -# adc = sitk.GetImageFromArray(adc) -# sitk.WriteImage(adc, f"adc_calc_b400_b800.nii.gz") - +# sitk.WriteImage(adc, f"adc_true.nii.gz") \ No newline at end of file diff --git a/scripts/calc_adc_b1400.py b/scripts/calc_adc_b1400.py new file mode 100755 index 0000000..c7edc63 --- /dev/null +++ b/scripts/calc_adc_b1400.py @@ -0,0 +1,91 @@ +import os +from os import path +from tqdm import tqdm +import numpy as np +import SimpleITK as sitk + +def calc_adc(b50,b400): + mean_dwi = (50 + 400) / 2 + mean_si = np.divide(np.add(np.log(b50), np.log(b400)), 2) + + denominator = np.multiply((50 - mean_dwi), np.subtract(np.log(b50), mean_si)) + np.multiply((400 - mean_dwi), np.subtract(np.log(b400), mean_si)) + numerator = np.power((50 - mean_dwi), 2) + np.power((400 - mean_dwi), 2) + adc_with_zeros = np.divide(denominator, numerator) * -1000000 + adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros)) + return adc + +def calc_high_b(b_value_high,b_value,b_image,ADC_map): + high_b = np.multiply(b_image, np.exp(np.multiply(np.subtract(b_value,b_value_high), (np.divide(ADC_map,1000000))))) + return high_b + +def append_new_line(file_name, text_to_append): + """Append given text as a new line at the end of file""" + # Open the file in append & read mode ('a+') + with open(file_name, "a+") as file_object: + # Move read cursor to the start of file. + file_object.seek(0) + # If file is not empty then append '\n' + data = file_object.read(100) + if len(data) > 0: + file_object.write("\n") + # Append text at the end of file + file_object.write(text_to_append) + + + +series = ['b50','b400'] +# series = ['t2'] + +DATA_DIR = "./../data/Nijmegen paths/" + +# Read the image paths from the data directory. +# Texts files are expected to have the name "[series_name].txt" +images, image_paths = {s: [] for s in series}, {} +segmentations = [] + +for s in series: + with open(path.join(DATA_DIR, f"{s}.txt"), 'r') as f: + image_paths[s] = [l.strip() for l in f.readlines()] +with open(path.join(DATA_DIR, f"seg.txt"), 'r') as f: + seg_paths = [l.strip() for l in f.readlines()] +num_images = len(seg_paths) + +for img_idx in tqdm(range(num_images)): + img_b50 = sitk.ReadImage(image_paths['b50'][img_idx], sitk.sitkFloat32) + arr_b50 = sitk.GetArrayFromImage(img_b50) + img_b400 = sitk.ReadImage(image_paths['b400'][img_idx], sitk.sitkFloat32) + arr_b400 = sitk.GetArrayFromImage(img_b400) + # img_b800 = sitk.ReadImage(image_paths['b800'][img_idx], sitk.sitkFloat32) + # arr_b800 = sitk.GetArrayFromImage(img_b800) + # img_t2 = sitk.ReadImage(image_paths['t2'][img_idx], sitk.sitkFloat32) + # img_t2 = sitk.GetArrayFromImage(img_t2) + # img_seg = sitk.ReadImage(seg_paths[img_idx], sitk.sitkFloat32) + # img_seg = sitk.GetArrayFromImage(img_seg) + + adc = calc_adc(arr_b50,arr_b400) + high_b = calc_high_b(1400,50,arr_b50,adc) + + adc = sitk.GetImageFromArray(adc) + high_b = sitk.GetImageFromArray(high_b) + # t2 = sitk.GetImageFromArray(img_t2) + # seg = sitk.GetImageFromArray(img_seg) + + patient_ID = os.path.split(os.path.dirname(os.path.dirname(os.path.dirname(os.path.dirname(image_paths['b50'][img_idx])))))[1] + path_adc = f'./../data/adc_calc_2/{patient_ID}.nii.gz' + path_high_b = f'./../data/b1400_calc_2/{patient_ID}.nii.gz' + # patient_ID = os.path.split(os.path.dirname(os.path.dirname(os.path.dirname(image_paths['t2'][img_idx]))))[1] + # path_t2 = f'./../data/t2_calc/{patient_ID}.nii.gz' + # path_seg = f'./../data/seg_calc/{patient_ID}.nii.gz' + + adc.CopyInformation(img_b50) + high_b.CopyInformation(img_b50) + sitk.WriteImage(adc, path_adc) + sitk.WriteImage(high_b, path_high_b) + # sitk.WriteImage(t2, path_t2) + # sitk.WriteImage(seg, path_seg) + + append_new_line('adccalc2.txt',path_adc) + append_new_line('b1400calc2.txt',path_high_b) + # append_new_line('t2calc.txt',path_t2) + # append_new_line('segcalc.txt',path_seg) + diff --git a/src/sfransen/DWI_exp/__init__.py b/src/sfransen/DWI_exp/__init__.py index 6a0390a..10452f3 100755 --- a/src/sfransen/DWI_exp/__init__.py +++ b/src/sfransen/DWI_exp/__init__.py @@ -3,3 +3,4 @@ from .callbacks import * from .helpers import * from .preprocessing_function import * from .unet import * +from .losses import * diff --git a/src/sfransen/DWI_exp/losses.py b/src/sfransen/DWI_exp/losses.py new file mode 100755 index 0000000..e6e746b --- /dev/null +++ b/src/sfransen/DWI_exp/losses.py @@ -0,0 +1,78 @@ +import tensorflow as tf +from tensorflow.image import ssim +from tensorflow.keras import backend as K + +def weighted_binary_cross_entropy(weights: dict, from_logits: bool = False): + ''' + Return a function for calculating weighted binary cross entropy + It should be used for multi-hot encoded labels + + # Example + y_true = tf.convert_to_tensor([1, 0, 0, 0, 0, 0], dtype=tf.int64) + y_pred = tf.convert_to_tensor([0.6, 0.1, 0.1, 0.9, 0.1, 0.], dtype=tf.float32) + weights = { + 0: 1., + 1: 2. + } + # with weights + loss_fn = get_loss_for_multilabels(weights=weights, from_logits=False) + loss = loss_fn(y_true, y_pred) + print(loss) + # tf.Tensor(0.6067193, shape=(), dtype=float32) + + # without weights + loss_fn = get_loss_for_multilabels() + loss = loss_fn(y_true, y_pred) + print(loss) + # tf.Tensor(0.52158177, shape=(), dtype=float32) + + # Another example + y_true = tf.convert_to_tensor([[0., 1.], [0., 0.]], dtype=tf.float32) + y_pred = tf.convert_to_tensor([[0.6, 0.4], [0.4, 0.6]], dtype=tf.float32) + weights = { + 0: 1., + 1: 2. + } + # with weights + loss_fn = get_loss_for_multilabels(weights=weights, from_logits=False) + loss = loss_fn(y_true, y_pred) + print(loss) + # tf.Tensor(1.0439969, shape=(), dtype=float32) + + # without weights + loss_fn = get_loss_for_multilabels() + loss = loss_fn(y_true, y_pred) + print(loss) + # tf.Tensor(0.81492424, shape=(), dtype=float32) + + @param weights A dict setting weights for 0 and 1 label. e.g. + { + 0: 1. + 1: 8. + } + For this case, we want to emphasise those true (1) label, + because we have many false (0) label. e.g. + [ + [0 1 0 0 0 0 0 0 0 1] + [0 0 0 0 1 0 0 0 0 0] + [0 0 0 0 1 0 0 0 0 0] + ] + + + + @param from_logits If False, we apply sigmoid to each logit + @return A function to calcualte (weighted) binary cross entropy + ''' + assert 0 in weights + assert 1 in weights + + def weighted_cross_entropy_fn(y_true, y_pred): + tf_y_true = tf.cast(y_true, dtype=y_pred.dtype) + tf_y_pred = tf.cast(y_pred, dtype=y_pred.dtype) + + weights_v = tf.where(tf.equal(tf_y_true, 1), weights[1], weights[0]) + ce = K.binary_crossentropy(tf_y_true, tf_y_pred, from_logits=from_logits) + loss = K.mean(tf.multiply(ce, weights_v)) + return loss + + return weighted_cross_entropy_fn \ No newline at end of file diff --git a/src/sfransen/FROC/cal_froc_from_np.py b/src/sfransen/FROC/cal_froc_from_np.py index 11c4e6f..e357b52 100755 --- a/src/sfransen/FROC/cal_froc_from_np.py +++ b/src/sfransen/FROC/cal_froc_from_np.py @@ -187,7 +187,7 @@ def preprocess_softmax(softmax: np.ndarray, def evaluate( y_true: np.ndarray, y_pred: np.ndarray, - min_overlap=0.02, + min_overlap=0.10, overlap_func: str = 'DSC', case_confidence: str = 'max', multiple_lesion_candidates_selection_criteria='overlap', @@ -382,7 +382,8 @@ def counts_from_lesion_evaluations( TP[i] = tp else: # extend FROC curve to infinity - TP[i] = TP[-2] + TP[i] = TP[-2] #note: aangepast stefan 11-04-2022 + # TP[i] = TP[-1] FP[i] = np.inf return TP, FP, thresholds, num_lesions diff --git a/src/sfransen/load_images.py b/src/sfransen/load_images.py index 96b6935..a348bd7 100755 --- a/src/sfransen/load_images.py +++ b/src/sfransen/load_images.py @@ -1,6 +1,7 @@ from typing import List import SimpleITK as sitk from sfransen.DWI_exp.helpers import * +import numpy as np def load_images_parrallel( image_paths: str, @@ -12,19 +13,31 @@ def load_images_parrallel( #resample mri_tra_s = resample(img_s, - min_shape=target_shape, + min_shape=(s+1 for s in target_shape), method=sitk.sitkNearestNeighbor, new_spacing=target_space) #center crop mri_tra_s = center_crop(mri_tra_s, shape=target_shape) #normalize - if seq != 'seg': - filter = sitk.NormalizeImageFilter() - mri_tra_s = filter.Execute(mri_tra_s) - else: - filter = sitk.BinaryThresholdImageFilter() - filter.SetLowerThreshold(1.0) - mri_tra_s = filter.Execute(mri_tra_s) + # if seq != 'seg': + # filter = sitk.NormalizeImageFilter() + # mri_tra_s = filter.Execute(mri_tra_s) + # else: + # filter = sitk.BinaryThresholdImageFilter() + # filter.SetLowerThreshold(1.0) + # mri_tra_s = filter.Execute(mri_tra_s) - return sitk.GetArrayFromImage(mri_tra_s).T \ No newline at end of file + # return sitk.GetArrayFromImage(mri_tra_s).T + + # Return sitk.Image instead of numpy np.ndarray. + + ### method trained in Unet + img_n = sitk.GetArrayFromImage(mri_tra_s).T + + if seq != 'seg': + image_return = (img_n - np.mean(img_n)) / ( 2* np.std(img_n)) + else: + image_return = np.clip(img_n, 0., 1.) + + return image_return