From dac43d7429262020e43537fb3593d4f66bd0c978 Mon Sep 17 00:00:00 2001 From: Stefan Date: Thu, 28 Apr 2022 10:11:10 +0200 Subject: [PATCH] added 10. in scripts which calculates froc with method Chris and optuna --- .gitignore | 4 +- scripts/10.optimize_thresholds_froc.py | 318 +++++++++++++++++++++++++ scripts/4.frocs.py | 3 +- scripts/5.Visualize_frocs.py | 16 +- scripts/calc_adc.py | 216 +++++++++-------- 5 files changed, 455 insertions(+), 102 deletions(-) create mode 100755 scripts/10.optimize_thresholds_froc.py diff --git a/.gitignore b/.gitignore index ccb8b4d..74a8f8c 100755 --- a/.gitignore +++ b/.gitignore @@ -148,4 +148,6 @@ cython_debug/ *.out sqliteDB/ *.db -/train_output/ \ No newline at end of file +/train_output/ +/umcglib/ +/k2s_umcg/ \ No newline at end of file diff --git a/scripts/10.optimize_thresholds_froc.py b/scripts/10.optimize_thresholds_froc.py new file mode 100755 index 0000000..2426144 --- /dev/null +++ b/scripts/10.optimize_thresholds_froc.py @@ -0,0 +1,318 @@ +from pickle import FALSE +from umcglib.froc import calculate_froc, plot_multiple_froc, partial_auc +from umcglib.binarize import dynamic_threshold +import optuna +import sqlite3 +from sfransen.utils_quintin import * +from os import path +from sfransen.DWI_exp.preprocessing_function import preprocess +import SimpleITK as sitk +import numpy as np +from sfransen.DWI_exp.callbacks import dice_coef +from sfransen.DWI_exp.losses import weighted_binary_cross_entropy +from tensorflow.keras.models import load_model +from tqdm import tqdm +from optuna.samplers import TPESampler +import argparse +import shutil +import os + +parser = argparse.ArgumentParser( + description='Calculate the froc metrics and store in froc_metrics.yml') +parser.add_argument('-experiment', + help='Title of experiment') +parser.add_argument('-series', '-s', + metavar='[series_name]', required=True, nargs='+', + help='List of series to include') +parser.add_argument('-fold', + default='', + help='List of series to include') +args = parser.parse_args() + +def does_table_exist(tablename: str, db_path: str): + + conn = sqlite3.connect(db_path) + c = conn.cursor() + + #get the count of tables with the name + c.execute(f''' SELECT count(name) FROM sqlite_master WHERE type='table' AND name='{tablename}' ''') + + does_exist = False + + #if the count is 1, then table exists + if c.fetchone()[0] == 1: + print(f"Table '{tablename}' exists.") + does_exist = True + else: + print(f"Table '{tablename}' does not exists.") + + #commit the changes to db + conn.commit() + #close the connection + conn.close() + return does_exist + +def load_or_create_study( + is_new_study: bool, + study_dir: str, +): + # Create an optuna if it does not exist. + storage = f"sqlite:///{study_dir}/{DB_FNAME}" + if is_new_study: + print(f"Creating a NEW study. With name: {storage}") + study = optuna.create_study(storage=storage, + study_name=study_dir, + direction='maximize', + sampler=TPESampler(n_startup_trials=N_STARTUP_TRIALS)) + else: + print(f"LOADING study {storage} from database file.") + study = optuna.load_study(storage=storage, + study_name=study_dir) + + return study + +def p_auc_froc_obj(trial, y_true_val, y_pred_val): + + dyn_thresh = trial.suggest_float('dyn_thresh', 0.0, 1.0) + min_conf = trial.suggest_float('min_conf', 0.0, 1.0) + + stats = calculate_froc(y_true=y_true_val, + y_pred=y_pred_val, + preprocess_func=dynamic_threshold, + dynamic_threshold_factor=dyn_thresh, + minimum_confidence=min_conf) + + sens, fpp = stats['sensitivity'], stats['fp_per_patient'] + p_auc_froc = partial_auc(sens, fpp, low=0.1, high=2.5) + + print(f"dyn_threshold: {dyn_thresh}, min_conf{min_conf}") + print(f"Trial {trial.number} pAUC FROC: {p_auc_froc}") + + return p_auc_froc + +def convert_np_to_list(flat_numpy_arr): + ans = [] + for elem in flat_numpy_arr: + ans.append(float(elem)) + return ans + +# >>>>>>>>> main <<<<<<<<<<<<< +# DB_FNAME = "calc_exp_t2_b1400_adc.db" +num_trials = 50 +N_STARTUP_TRIALS = 10 + +SERIES = args.series +series_ = '_'.join(SERIES) +EXPERIMENT = args.experiment + +DB_FNAME = f'{EXPERIMENT}_{series_}_{args.fold}.db' +MODEL_PATH = f'./../train_output/{EXPERIMENT}_{series_}/models/{EXPERIMENT}_{series_}.h5' +YAML_DIR = f'./../train_output/{EXPERIMENT}_{series_}' +DATA_DIR = "./../data/Nijmegen paths/" +TARGET_SPACING = (0.5, 0.5, 3) +INPUT_SHAPE = (192, 192, 24, len(SERIES)) +IMAGE_SHAPE = INPUT_SHAPE[:3] + +DATA_SPLIT_INDEX = read_yaml_to_dict(f'./../data/Nijmegen paths/train_val_test_idxs.yml') +TEST_INDEX = DATA_SPLIT_INDEX['val_set0'] +print("test test_index",TEST_INDEX[:5]) + +############ load data en preprocess / old method +# print(">>>>> read images <<<<<<<<<<") +# 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) + +# images = [] +# images_list = [] +# segmentations = [] + +# # 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) + +# images_list = np.transpose(images_list, (0, 2, 3, 4, 1)) +# print("shape of segmentations is",np.shape(segmentations)) +# 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, +# 'weighted_cross_entropy_fn':weighted_binary_cross_entropy +# } +# reconstructed_model = load_model(MODEL_PATH, custom_objects=dependencies) +# # reconstructed_model.summary(line_length=120) + +# # make predictions on all TEST_INDEX +# print(' >>>>>>> START prediction <<<<<<<<<') +# predictions_blur = reconstructed_model.predict(images_list, batch_size=1) + +# ############# preprocess ################# +# # preprocess predictions by removing the blur and making individual blobs +# print('>>>>>>>> START preprocess') + +# # def move_dims(arr): +# # # UMCG numpy dimensions convention: dims = (batch, width, heigth, depth) +# # # Joeran numpy dimensions convention: dims = (batch, depth, heigth, width) +# # arr = np.moveaxis(arr, 3, 1) +# # arr = np.moveaxis(arr, 3, 2) +# # return arr + +# # # Joeran has his numpy arrays ordered differently. +# # predictions_blur = move_dims(np.squeeze(predictions_blur)) +# # segmentations = move_dims(np.squeeze(segmentations)) + +# y_pred_val = np.squeeze(predictions_blur) +# y_true_val = segmentations +study_dir = f"./../sqliteDB/optuna_dbs" + +check_for_file = path.isfile(f"{study_dir}/{DB_FNAME}") +if check_for_file == False: + shutil.copyfile(f"{study_dir}/dyn_thres_min_conf_opt_OG.db", f"{study_dir}/{DB_FNAME}") + +table_exists = does_table_exist('trials', f"{study_dir}/{DB_FNAME}") +study = load_or_create_study(is_new_study=not table_exists, study_dir=study_dir) + +# # dyn_thresh = study.best_trial.params['dyn_thresh'] +# # min_conf = study.best_trial.params['min_conf'] +# dyn_thresh = 0.4 +# min_conf = 0.01 + +# # print("step 1:",np.shape(y_pred_val)) +# stats = calculate_froc(y_true=y_true_val, +# y_pred=y_pred_val, +# preprocess_func=dynamic_threshold, +# dynamic_threshold_factor=dyn_thresh, +# minimum_confidence=min_conf) + +# sens, fpp = stats['sensitivity'], stats['fp_per_patient'] +# p_auc = partial_auc(sens, fpp, low=0.1, high=2.5) + +# print(f"the p_auc with old setting is: {p_auc}" ) + +# # Try to find the best value for the dynamic threshold and min_confidence +# opt_func = lambda trail: p_auc_froc_obj(trail, y_true_val, y_pred_val) +# study.optimize(opt_func, n_trials=num_trials) + +dyn_thresh = study.best_trial.params['dyn_thresh'] +min_conf = study.best_trial.params['min_conf'] +print(f"done. best dyn_thresh: {dyn_thresh} . Best min_conf: {min_conf}") + +########## dump dict to yaml of best froc curve ############# + +######## gooi dit in functie ############### +DATA_SPLIT_INDEX = read_yaml_to_dict(f'./../data/Nijmegen paths/train_val_test_idxs.yml') +TEST_INDEX = DATA_SPLIT_INDEX['test_set0'] +print("test test_index",TEST_INDEX[:5]) + +############ load data en preprocess / old method +print(">>>>> read images <<<<<<<<<<") +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) + +images = [] +images_list = [] +segmentations = [] + +# 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) + +images_list = np.transpose(images_list, (0, 2, 3, 4, 1)) +print("shape of segmentations is",np.shape(segmentations)) +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, + 'weighted_cross_entropy_fn':weighted_binary_cross_entropy +} +reconstructed_model = load_model(MODEL_PATH, custom_objects=dependencies) +# reconstructed_model.summary(line_length=120) + +# make predictions on all TEST_INDEX +print(' >>>>>>> START prediction <<<<<<<<<') +predictions_blur = reconstructed_model.predict(images_list, batch_size=1) + +############# preprocess ################# +# preprocess predictions by removing the blur and making individual blobs +print('>>>>>>>> START preprocess') + +# def move_dims(arr): +# # UMCG numpy dimensions convention: dims = (batch, width, heigth, depth) +# # Joeran numpy dimensions convention: dims = (batch, depth, heigth, width) +# arr = np.moveaxis(arr, 3, 1) +# arr = np.moveaxis(arr, 3, 2) +# return arr + +# # Joeran has his numpy arrays ordered differently. +# predictions_blur = move_dims(np.squeeze(predictions_blur)) +# segmentations = move_dims(np.squeeze(segmentations)) + +y_pred_val = np.squeeze(predictions_blur) +y_true_val = segmentations +########### einde functie ############ + +stats = calculate_froc(y_true=y_true_val, + y_pred=y_pred_val, + preprocess_func=dynamic_threshold, + dynamic_threshold_factor=dyn_thresh, + minimum_confidence=min_conf) + + +subject_idxs = list(range(len(y_true_val))) +metrics = { + "num_patients": int(stats['num_patients']), + "auroc": int(stats['patient_auc']), + 'tpr': convert_np_to_list(stats['roc_tpr']), + 'fpr': convert_np_to_list(stats['roc_fpr']), + "roc_true": convert_np_to_list(stats['roc_patient_level_label'][s] for s in subject_idxs), + "roc_pred": convert_np_to_list(stats['roc_patient_level_conf'][s] for s in subject_idxs), + "num_lesions": int(stats['num_lesions']), + "thresholds": convert_np_to_list(stats['thresholds']), + "sensitivity": convert_np_to_list(stats['sensitivity']), + "FP_per_case": convert_np_to_list(stats['fp_per_patient']), + "precision": convert_np_to_list(stats['precision']), + "recall": convert_np_to_list(stats['recall']), + "AP": int(stats['average_precision']), +} +dump_dict_to_yaml(metrics, YAML_DIR, "froc_metrics_optuna_test", verbose=True) + + + + + + diff --git a/scripts/4.frocs.py b/scripts/4.frocs.py index a5d7afe..eae0e68 100755 --- a/scripts/4.frocs.py +++ b/scripts/4.frocs.py @@ -7,7 +7,8 @@ import multiprocessing from functools import partial import os from os import path - +from tqdm import tqdm +import argparse from sfransen.utils_quintin import * from sfransen.DWI_exp.helpers import * diff --git a/scripts/5.Visualize_frocs.py b/scripts/5.Visualize_frocs.py index 93a68fc..9fd153e 100755 --- a/scripts/5.Visualize_frocs.py +++ b/scripts/5.Visualize_frocs.py @@ -3,7 +3,7 @@ from sfransen.utils_quintin import * import matplotlib.pyplot as plt import argparse import matplotlib.ticker as tkr -from sfransen.FROC.p_auc import partial_auc +from umcglib.froc.p_auc import partial_auc parser = argparse.ArgumentParser( description='Visualise froc results') @@ -15,6 +15,9 @@ parser.add_argument('--experiment', '-s', metavar='[series_name]', required=True, nargs='+', help='List of series to include, must correspond with' + "path files in ./data/") +parser.add_argument('-yaml_metric', + help='List of series to include, must correspond with' + + "path files in ./data/") args = parser.parse_args() if args.comparison: @@ -24,6 +27,7 @@ else: colors = ['r','b','g','k','y','c'] plot_type = ['-','-','-','-','-','-'] +yaml_metric = args.yaml_metric experiments = args.experiment print(experiments) experiment_path = [] @@ -34,11 +38,11 @@ paufroc = [] fig = plt.figure(1) ax = fig.add_subplot(111) for idx in range(len(args.experiment)): - experiment_path = f'./../train_output/{experiments[idx]}/froc_metrics_focal_10.yml' + experiment_path = f'./../train_output/{experiments[idx]}/froc_metrics_{yaml_metric}.yml' experiment_metrics = read_yaml_to_dict(experiment_path) - pfroc = partial_auc(experiment_metrics["sensitivity"],experiment_metrics["FP_per_case"]) - paufroc.append(round(pfroc,2) + pfroc = partial_auc(experiment_metrics["sensitivity"],experiment_metrics["FP_per_case"],low=0.1, high=2.5) + paufroc.append(round(pfroc,2)) plt.plot(experiment_metrics["FP_per_case"], experiment_metrics["sensitivity"],color=colors[idx],linestyle=plot_type[idx]) ax.set(xscale="log") @@ -50,7 +54,7 @@ for idx in range(len(args.experiment)): ax.axes.xaxis.set_major_locator(tkr.FixedLocator([0,0.1,1,3])) for idx in range(len(args.experiment)): - experiment_path = f'./../train_output/{experiments[idx]}/froc_metrics_focal_10.yml' + experiment_path = f'./../train_output/{experiments[idx]}/froc_metrics_{yaml_metric}.yml' experiment_metrics = read_yaml_to_dict(experiment_path) auroc.append(round(experiment_metrics['auroc'],3)) @@ -86,4 +90,4 @@ plt.legend(experiments_auroc,loc='lower right') plt.xlabel('False positive rate') plt.ylabel('True positive rate') plt.grid() -plt.savefig(f"./../train_output/ROC_{args.saveas}.png", dpi=300) +plt.savefig(f"./../train_output/ROC_{args.saveas}.png", dpi=300) \ No newline at end of file diff --git a/scripts/calc_adc.py b/scripts/calc_adc.py index 97244e9..c76f631 100755 --- a/scripts/calc_adc.py +++ b/scripts/calc_adc.py @@ -1,12 +1,17 @@ import numpy as np import SimpleITK as sitk import matplotlib.pyplot as plt -######## load images ############# -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' +# ######## load images ############# +# 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_adc = '/data/pca-rad/datasets/anonymized_mri/only_nii_directory/110-M-110/2015NB/dADC_0-50-500-1000/702_.nii.gz' +path_b0 = '/data/pca-rad/datasets/anonymized_mri/only_nii_directory/110-M-110/2015NB/Diffusie/b-0/701_.nii.gz' +path_b800 = '/data/pca-rad/datasets/anonymized_mri/only_nii_directory/110-M-110/2015NB/Diffusie/b-800/701_.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' @@ -20,14 +25,14 @@ path_adc = '/data/pca-rad/datasets/radboud_new/pat0634/2016/dADC/nifti_image.nii # 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_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) +b0_img = sitk.ReadImage(path_b0, sitk.sitkFloat32) +b0 = sitk.GetArrayFromImage(b0_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) +# 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) @@ -39,21 +44,21 @@ def show_img(greyscale_img): fig.savefig(path, dpi=300, bbox_inches='tight') def calc_adc(b50, b400, b800): + "Calculates the adc based on b50, b400 and b800 DWI images/arrays." mean_dwi = (50 + 400 + 800) / 3 mean_si = np.divide(np.add(np.add(np.log(b50), np.log(b400)), np.log(b800)), 3) - 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_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 - mean_si = np.divide(np.add(np.log(b50), np.log(b800)), 2) +def calc_adc_1(b0,b800): + mean_dwi = (0 + 800) / 2 + mean_si = np.divide(np.add(np.log(b0), np.log(b800)), 2) - 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) + denominator = np.multiply((0 - mean_dwi), np.subtract(np.log(b0), mean_si)) + np.multiply((800 - mean_dwi), np.subtract(np.log(b800), mean_si)) + numerator = np.power((0 - mean_dwi), 2) + np.power((800 - 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 @@ -78,92 +83,115 @@ def calc_adc_3(b400,b800): 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))))) +def calc_high_b(b_value_high,b_value,b_image,adc): + """ + Calculates a high b-value image. + b_value_high = the requered b-value integer + b_value = the b_value integer used as reference image + b_value = the corresponding array + adc = the corresponding adc array + """" + high_b = np.multiply(b_image, np.exp(np.multiply(np.subtract(b_value,b_value_high), (np.divide(adc,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) +adc_0_800 = calc_adc_1(b0,b800) +high_b_1400_0 = calc_high_b(1400,0,b0,adc_0_800) -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) +adc_0_800 = sitk.GetImageFromArray(adc_0_800) +adc_0_800.CopyInformation(b0_img) +sitk.WriteImage(adc_0_800, f"../train_output/test.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) +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") -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') -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') +quit() +# 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) -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') +# 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) -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') +# 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) + +# 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') + +# 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)