p307800/fast-mri
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

working froc, changes in preprocessing

Browse Source
This commit is contained in:
Stefan 2022-04-21 14:29:36 +02:00
parent 10cdf43509
commit ee0a672f0b
11 changed files with 463 additions and 124 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
scripts/1.U-net_chris.py
View File

@ -13,6 +13,7 @@ from sfransen.utils_quintin import *
from sfransen.utils_quintin import * from sfransen.utils_quintin import *
from sfransen.DWI_exp import IntermediateImages, dice_coef from sfransen.DWI_exp import IntermediateImages, dice_coef
from sfransen.DWI_exp.preprocessing_function import preprocess from sfransen.DWI_exp.preprocessing_function import preprocess
from sfransen.DWI_exp.losses import weighted_binary_cross_entropy
import yaml import yaml
import numpy as np import numpy as np
from tqdm import tqdm 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 MODEL_SELECTION_DIRECTION = "min" # Change to 'max' if higher value is better
EARLY_STOPPING_METRIC = 'val_loss' EARLY_STOPPING_METRIC = 'val_loss'
EARLY_STOPPING_DIRECTION = 'min' 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 # Training configuration
# add metric ROC_AUC # add metric ROC_AUC
TRAINING_METRICS = ["binary_crossentropy", "binary_accuracy", dice_coef] 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) optimizer = Adam(learning_rate=INITIAL_LEARNING_RATE)
# Create folder structure in the output directory # 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) num_images = len(seg_paths)
# Read and preprocess each of the paths for each series, and the segmentations. # 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) img_s = {s: sitk.ReadImage(image_paths[s][img_idx], sitk.sitkFloat32)
for s in args.series} for s in args.series}
seg_s = sitk.ReadImage(seg_paths[img_idx], sitk.sitkFloat32) 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])}") print_(f"The shape of valid_data label = {np.shape(valid_data[1])}")
callbacks = [ callbacks = [
EarlyStopping(
monitor=EARLY_STOPPING_METRIC,
mode=EARLY_STOPPING_DIRECTION,
patience=EARLY_STOPPING,
verbose=2),
ModelCheckpoint( ModelCheckpoint(
filepath=path.join(PROJECT_DIR, "models", JOB_NAME + ".h5"), filepath=path.join(PROJECT_DIR, "models", JOB_NAME + ".h5"),
monitor=MODEL_SELECTION_METRIC, monitor=MODEL_SELECTION_METRIC,
@ -185,6 +188,11 @@ callbacks = [
save_best_only=True), save_best_only=True),
CSVLogger( CSVLogger(
filename=path.join(PROJECT_DIR, "logs", f"{JOB_NAME}.csv")), 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( IntermediateImages(
validation_set=valid_data, validation_set=valid_data,
sequences=args.series, sequences=args.series,

131
scripts/4.frocs.py
View File

@ -6,6 +6,8 @@ import numpy as np
import multiprocessing import multiprocessing
from functools import partial from functools import partial
import os import os
from os import path
from sfransen.utils_quintin import * from sfransen.utils_quintin import *
from sfransen.DWI_exp.helpers 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.blob_preprocess import *
from sfransen.FROC.cal_froc_from_np import * from sfransen.FROC.cal_froc_from_np import *
from sfransen.load_images import load_images_parrallel from sfransen.load_images import load_images_parrallel
from sfransen.DWI_exp.losses import weighted_binary_cross_entropy
parser = argparse.ArgumentParser( parser = argparse.ArgumentParser(
@ -36,7 +39,7 @@ SERIES = args.series
series_ = '_'.join(args.series) series_ = '_'.join(args.series)
EXPERIMENT = args.experiment 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_}' YAML_DIR = f'./../train_output/{EXPERIMENT}_{series_}'
IMAGE_DIR = f'./../train_output/{EXPERIMENT}_{series_}' IMAGE_DIR = f'./../train_output/{EXPERIMENT}_{series_}'
@ -54,8 +57,46 @@ N_CPUS = 12
########## load images in parrallel ############## ########## load images in parrallel ##############
print_(f"> Loading images into RAM...") 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 = {} image_paths = {}
for s in SERIES: for s in SERIES:
with open(path.join(DATA_DIR, f"{s}.txt"), 'r') as f: with open(path.join(DATA_DIR, f"{s}.txt"), 'r') as f:
image_paths[s] = [l.strip() for l in f.readlines()] 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()] seg_paths = [l.strip() for l in f.readlines()]
num_images = len(seg_paths) 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 = [] images = []
for s in SERIES: images_list = []
image_paths_seq = image_paths[s] segmentations = []
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 = 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 images_list = np.transpose(images_list, (0, 2, 3, 4, 1))
seg_paths_index = np.asarray(seg_paths)[TEST_INDEX]
data_list = pool.map(partial_seg,seg_paths_index) print('>>>>> size image_list nmr 2:', np.shape(images_list), '. equal to: (5, 192, 192, 24, 3)?')
segmentations = np.stack(data_list, axis=0)
########### load module ################## ########### load module ##################
print(' >>>>>>> LOAD MODEL <<<<<<<<<') print(' >>>>>>> LOAD MODEL <<<<<<<<<')
dependencies = { 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 = load_model(MODEL_PATH, custom_objects=dependencies)
# reconstructed_model.summary(line_length=120) # reconstructed_model.summary(line_length=120)
@ -127,16 +163,51 @@ predictions = zeros
# perform Froc # perform Froc
metrics = evaluate(y_true=segmentations, y_pred=predictions) 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 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 # save image nmr 3
img_s = sitk.GetImageFromArray(predictions_blur[3].squeeze()) 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()) 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()) img_s = sitk.GetImageFromArray(segmentations[3].squeeze())
sitk.WriteImage(img_s, f"{IMAGE_DIR}/segmentations_001.nii.gz") 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")

2
scripts/5.Visualize_frocs.py
View File

@ -28,7 +28,7 @@ experiment_path = []
experiment_metrics = {} experiment_metrics = {}
auroc = [] auroc = []
for idx in range(len(args.experiment)): 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) experiment_metrics = read_yaml_to_dict(experiment_path)
auroc.append(round(experiment_metrics['auroc'],3)) auroc.append(round(experiment_metrics['auroc'],3))

11
scripts/6.saliency_map.py
View File

@ -9,6 +9,8 @@ from sfransen.utils_quintin import *
from sfransen.DWI_exp import preprocess from sfransen.DWI_exp import preprocess
from sfransen.DWI_exp.helpers import * from sfransen.DWI_exp.helpers import *
from sfransen.DWI_exp.callbacks import dice_coef 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.blob_preprocess import *
from sfransen.FROC.cal_froc_from_np import * from sfransen.FROC.cal_froc_from_np import *
from sfransen.load_images import load_images_parrallel 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)) INPUT_SHAPE = (192, 192, 24, len(SERIES))
IMAGE_SHAPE = INPUT_SHAPE[:3] IMAGE_SHAPE = INPUT_SHAPE[:3]
froc_metrics = read_yaml_to_dict(f'{YAML_DIR}/froc_metrics.yml') # froc_metrics = read_yaml_to_dict(f'{YAML_DIR}/froc_metrics.yml')
top_10_idx = np.argsort(froc_metrics['roc_pred'])[-1 :] # 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') DATA_SPLIT_INDEX = read_yaml_to_dict('./../data/Nijmegen paths/train_val_test_idxs.yml')
TEST_INDEX = DATA_SPLIT_INDEX['val_set0'] 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] TEST_INDEX_image = [371]
N_CPUS = 12 N_CPUS = 12
@ -93,7 +95,8 @@ segmentations = np.stack(data_list, axis=0)
print(' >>>>>>> LOAD MODEL <<<<<<<<<') print(' >>>>>>> LOAD MODEL <<<<<<<<<')
dependencies = { 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 = load_model(MODEL_PATH, custom_objects=dependencies)
# reconstructed_model.layers[-1].activation = tf.keras.activations.linear # reconstructed_model.layers[-1].activation = tf.keras.activations.linear

4
scripts/8.Visualize_training.py
View File

@ -29,9 +29,9 @@ df = pd.read_csv(f'{folder_input}')
# read csv file # read csv file
for metric in df: for metric in df:
if not metric == 'epoch': if not metric == 'epoch':
# if metric == 'loss' or metric == 'val_loss': if metric == 'loss' or metric == 'val_loss':
plt.plot(df['epoch'], df[metric], label=metric) plt.plot(df['epoch'], df[metric], label=metric)
# plt.ylim(ymin=0,ymax=0.01) plt.ylim(ymin=0,ymax=0.01)

207
scripts/calc_adc.py
View File

@ -2,34 +2,34 @@ import numpy as np
import SimpleITK as sitk import SimpleITK as sitk
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
######## load images ############# ######## load images #############
# path_b50 = '/data/pca-rad/datasets/radboud_new/pat0351/2016/diffusie_cro/b-50/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/pat0351/2016/diffusie_cro/b-400/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/pat0351/2016/diffusie_cro/b-800/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/pat0351/2016/diffusie_cro/b-1400/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/pat0351/2016/dADC/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_b50 = 'X:/sfransen/train_output/adc_exp/b50_true.nii.gz'
# path_b400 = 'X:/sfransen/train_output/adc_exp/b400.nii.gz' # path_b400 = 'X:/sfransen/train_output/adc_exp/b400_true.nii.gz'
# path_b800 = 'X:/sfransen/train_output/adc_exp/b800.nii.gz' # path_b800 = 'X:/sfransen/train_output/adc_exp/b800_true.nii.gz'
# path_b1400 = 'X:/sfransen/train_output/adc_exp/b1400.nii.gz' # path_b1400 = 'X:/sfransen/train_output/adc_exp/b1400_true.nii.gz'
# path_adc = 'X:/sfransen/train_output/adc_exp/adc.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_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_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_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_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_adc = '/data/pca-rad/sfransen/train_output/adc_exp/adc_true.nii.gz'
b50 = sitk.ReadImage(path_b50, sitk.sitkFloat32) b50_img = sitk.ReadImage(path_b50, sitk.sitkFloat32)
b50 = sitk.GetArrayFromImage(b50) b50 = sitk.GetArrayFromImage(b50_img)
b400 = sitk.ReadImage(path_b400, sitk.sitkFloat32) b400_img = sitk.ReadImage(path_b400, sitk.sitkFloat32)
b400 = sitk.GetArrayFromImage(b400) b400 = sitk.GetArrayFromImage(b400_img)
b800 = sitk.ReadImage(path_b800, sitk.sitkFloat32) b800_img = sitk.ReadImage(path_b800, sitk.sitkFloat32)
b800 = sitk.GetArrayFromImage(b800) b800 = sitk.GetArrayFromImage(b800_img)
b1400 = sitk.ReadImage(path_b1400, sitk.sitkFloat32) b1400_img = sitk.ReadImage(path_b1400, sitk.sitkFloat32)
b1400 = sitk.GetArrayFromImage(b1400) b1400_original = sitk.GetArrayFromImage(b1400_img)
adc = sitk.ReadImage(path_adc, sitk.sitkFloat32) adc_img = sitk.ReadImage(path_adc, sitk.sitkFloat32)
adc = sitk.GetArrayFromImage(adc) adc_original = sitk.GetArrayFromImage(adc_img)
def show_img(greyscale_img): def show_img(greyscale_img):
fig = plt.figure() 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)) 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) numerator = np.power((50 - mean_dwi), 2) + np.power((400 - mean_dwi), 2) + np.power((800 - mean_dwi), 2)
adc = np.divide(denominator, numerator) adc_with_zeros = np.divide(denominator, numerator) * -1000000
return adc * -1000000 adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros))
return adc
def calc_adc_1(b50,b800): def calc_adc_1(b50,b800):
mean_dwi = (50 + 800) / 2 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)) 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) numerator = np.power((50 - mean_dwi), 2) + np.power((800 - mean_dwi), 2)
adc = np.divide(denominator, numerator) adc_with_zeros = np.divide(denominator, numerator) * -1000000
return adc * -1000000 adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros))
return adc
def calc_adc_2(b50,b400): def calc_adc_2(b50,b400):
mean_dwi = (50 + 400) / 2 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)) 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) numerator = np.power((50 - mean_dwi), 2) + np.power((400 - mean_dwi), 2)
adc = np.divide(denominator, numerator) adc_with_zeros = np.divide(denominator, numerator) * -1000000
return adc * -1000000 adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros))
return adc
def calc_adc_3(b400,b800): def calc_adc_3(b400,b800):
mean_dwi = (400 + 800) / 2 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)) 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) numerator = np.power((400 - mean_dwi), 2) + np.power((800 - mean_dwi), 2)
adc = np.divide(denominator, numerator) adc_with_zeros = np.divide(denominator, numerator) * -1000000
return adc * -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): 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 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) high_b_1400_50 = calc_high_b(1400,50,b50,adc_50_800)
sitk.WriteImage(high_b, f"./../train_output/adc_exp/b1400_true.nii.gz") 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_1400_400 = calc_high_b(1400,400,b400,adc_50_800)
high_b = sitk.GetImageFromArray(high_b) high_b_1400_800 = calc_high_b(1400,800,b800,adc_50_800)
sitk.WriteImage(high_b, f"./../train_output/adc_exp/b1400_ref_b50.nii.gz")
high_b = calc_high_b(1400,400,b400,adc) fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
high_b = sitk.GetImageFromArray(high_b) fig.suptitle('ADC calculated with b50 b400 b800')
sitk.WriteImage(high_b, f"./../train_output/adc_exp/b1400_ref_b400.nii.gz") 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) fig, (ax1, ax2, ax3) = plt.subplots(1, 3)
high_b = sitk.GetImageFromArray(high_b) fig.suptitle('ADC calculated with b50 b800')
sitk.WriteImage(high_b, f"./../train_output/adc_exp/b1400_ref_b800.nii.gz") 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) # b50 = sitk.GetImageFromArray(b50)
# sitk.WriteImage(b50, "./../train_output/adc_exp/b50_true.nii.gz") # 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) # b1400 = sitk.GetImageFromArray(b1400)
# sitk.WriteImage(b1400,f"./../train_output/adc_exp/b1400_true.nii.gz") # sitk.WriteImage(b1400,f"./../train_output/adc_exp/b1400_true.nii.gz")
# adc = sitk.GetImageFromArray(adc) # adc = sitk.GetImageFromArray(adc)
# sitk.WriteImage(adc, f"adc_true.nii.gz") # 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")

91
scripts/calc_adc_b1400.py Executable file
View File

@ -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)

1
src/sfransen/DWI_exp/__init__.py
View File

@ -3,3 +3,4 @@ from .callbacks import *
from .helpers import * from .helpers import *
from .preprocessing_function import * from .preprocessing_function import *
from .unet import * from .unet import *
from .losses import *

78
src/sfransen/DWI_exp/losses.py Executable file
View File

@ -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

5
src/sfransen/FROC/cal_froc_from_np.py
View File

@ -187,7 +187,7 @@ def preprocess_softmax(softmax: np.ndarray,
def evaluate( def evaluate(
y_true: np.ndarray, y_true: np.ndarray,
y_pred: np.ndarray, y_pred: np.ndarray,
min_overlap=0.02, min_overlap=0.10,
overlap_func: str = 'DSC', overlap_func: str = 'DSC',
case_confidence: str = 'max', case_confidence: str = 'max',
multiple_lesion_candidates_selection_criteria='overlap', multiple_lesion_candidates_selection_criteria='overlap',
@ -382,7 +382,8 @@ def counts_from_lesion_evaluations(
TP[i] = tp TP[i] = tp
else: else:
# extend FROC curve to infinity # 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 FP[i] = np.inf
return TP, FP, thresholds, num_lesions return TP, FP, thresholds, num_lesions

31
src/sfransen/load_images.py
View File

@ -1,6 +1,7 @@
from typing import List from typing import List
import SimpleITK as sitk import SimpleITK as sitk
from sfransen.DWI_exp.helpers import * from sfransen.DWI_exp.helpers import *
import numpy as np
def load_images_parrallel( def load_images_parrallel(
image_paths: str, image_paths: str,
@ -12,19 +13,31 @@ def load_images_parrallel(
#resample #resample
mri_tra_s = resample(img_s, mri_tra_s = resample(img_s,
min_shape=target_shape, min_shape=(s+1 for s in target_shape),
method=sitk.sitkNearestNeighbor, method=sitk.sitkNearestNeighbor,
new_spacing=target_space) new_spacing=target_space)
#center crop #center crop
mri_tra_s = center_crop(mri_tra_s, shape=target_shape) mri_tra_s = center_crop(mri_tra_s, shape=target_shape)
#normalize #normalize
if seq != 'seg': # if seq != 'seg':
filter = sitk.NormalizeImageFilter() # filter = sitk.NormalizeImageFilter()
mri_tra_s = filter.Execute(mri_tra_s) # mri_tra_s = filter.Execute(mri_tra_s)
else: # else:
filter = sitk.BinaryThresholdImageFilter() # filter = sitk.BinaryThresholdImageFilter()
filter.SetLowerThreshold(1.0) # filter.SetLowerThreshold(1.0)
mri_tra_s = filter.Execute(mri_tra_s) # mri_tra_s = filter.Execute(mri_tra_s)
return sitk.GetArrayFromImage(mri_tra_s).T # 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