From 03bd157664892c383e165ddb83f57c08dc969420 Mon Sep 17 00:00:00 2001 From: Stefan Date: Fri, 13 May 2022 13:30:23 +0200 Subject: [PATCH] start working on simulating averages & the miccai challenge: h5.gz > npy, build hists in parrallel --- scripts/12.simulate_averages.py | 42 +++++ scripts/miccai_calc_histograms.py | 71 ++++++++ scripts/miccai_rss_recon.py | 48 ++++++ scripts/miccai_rss_recon_undersampled.py | 66 ++++++++ scripts/test.py | 94 +++++++++++ scripts/tmp.py | 199 ++++++++++++++++++++--- 6 files changed, 496 insertions(+), 24 deletions(-) create mode 100755 scripts/12.simulate_averages.py create mode 100755 scripts/miccai_calc_histograms.py create mode 100755 scripts/miccai_rss_recon.py create mode 100755 scripts/miccai_rss_recon_undersampled.py create mode 100755 scripts/test.py diff --git a/scripts/12.simulate_averages.py b/scripts/12.simulate_averages.py new file mode 100755 index 0000000..797d65d --- /dev/null +++ b/scripts/12.simulate_averages.py @@ -0,0 +1,42 @@ +import numpy as np +import SimpleITK as sitk + +#paths +# path_b800 = '/data/pca-rad/sfransen/UMCG_009.nii.gz' +# path_seg_noise = '/data/pca-rad/sfransen/seg_umcg_0003.nii.gz' + +path_b800 = '/data/pca-rad/sfransen/b800_003.nii.gz' +path_seg_noise = '/data/pca-rad/sfransen/segment_003.nii.gz' +#load image/seg_signal/seg_noise +b800_img = sitk.ReadImage(path_b800, sitk.sitkFloat32) +b800_arr = sitk.GetArrayFromImage(b800_img) + +noise_img = sitk.ReadImage(path_seg_noise, sitk.sitkFloat32) +noise_arr = sitk.GetArrayFromImage(noise_img) + +noise_list = np.matrix.flatten(np.multiply(noise_arr,b800_arr)) +#calculate standard deviation +values = [i for i in np.matrix.flatten(np.multiply(noise_arr,b800_arr)) if i > 0] +std_old = np.std(values) +print('std old',std_old) +A1 = 12 +A2 = 4 +#calculate new std +std_new = np.multiply(std_old,np.divide(np.sqrt(A1),np.sqrt(A2))) +print('std new',std_new) +# std_2 = np.divide(std_1,np.sqrt(2)) +std_noise = np.sqrt(np.abs(np.power(std_old,2)-np.power(std_new,2))) +# std_noise = np.multiply(std_1,std_2) +print('std_noise',std_noise) +noise = np.random.normal(0,std_noise,np.shape(b800_arr)) +b800_arr_noise = np.add(b800_arr,noise) + +values_check = [i for i in np.matrix.flatten(np.multiply(noise_arr,b800_arr_noise)) if i > 0] +std_new_check = np.std(values_check) +print('std new check',std_new_check) + +b800_img_2 = sitk.GetImageFromArray(b800_arr_noise) +b800_img_2.CopyInformation(b800_img) +sitk.WriteImage(b800_img_2, f"../test_4.nii.gz") +b800_img = sitk.ReadImage(path_b800, sitk.sitkFloat32) +sitk.WriteImage(b800_img, f"../test_true.nii.gz") diff --git a/scripts/miccai_calc_histograms.py b/scripts/miccai_calc_histograms.py new file mode 100755 index 0000000..be0f585 --- /dev/null +++ b/scripts/miccai_calc_histograms.py @@ -0,0 +1,71 @@ +from traceback import FrameSummary +import numpy as np +import SimpleITK as sitk +from os import listdir +from pip import main +from scipy.fftpack import fftshift, ifftshift, ifftn +from umcglib.utils import apply_parallel +import time +import h5py +import matplotlib.pyplot as plt + + +mypath = f'/data/pca-rad/datasets/miccai_2022/K2S_MICCAI2022_GRP/train/data/TBrecon1/train/untarred' +files = [f for f in listdir(mypath)] + +dict = { + 'background': [], + 'femoral_cartilage': [], + 'tibial_cartilage': [], + 'patellar_cartilage': [], + 'femur': [], + 'tibia': [], + 'patella': [] +} + +for file_idx in range(300): + seg = [] + patient_id = files[file_idx] + seg = np.load(f'{mypath}/{patient_id}/{patient_id}_seg.npy') + # # 0: background; 1: femoral cartilage; 2: tibial cartilage; 3: patellar cartilage; 4: femur; 5: tibia; 6: patella. + + dict['background'].append(np.sum(seg == 0)) + dict['femoral_cartilage'].append(np.sum(seg == 1)) + dict['tibial_cartilage'].append(np.sum(seg == 2)) + dict['patellar_cartilage'].append(np.sum(seg == 3)) + dict['femur'].append(np.sum(seg == 4)) + dict['tibia'].append(np.sum(seg == 5)) + dict['patella'].append(np.sum(seg == 6)) + + print('done:',file_idx,' ',f'{mypath}/{patient_id}/{patient_id}_seg.npy') + +np.save('../dict.npy',dict) + +for keys in dict: + data = dict[f'keys'] + plt.hist(data) + plt.title(f"historgram of {keys}") + plt.savefig(f"../{keys}.png", dpi=300) + print('done ',keys) + + + +# print(f.keys()) +# print(np.shape(f['us_mask'][()])) +# print(type(f['us_mask'][()])) + + + + + + + +# quit() +# seg = sitk.GetImageFromArray(np.squeeze(seg)) +# sitk.WriteImage(seg, f"../test_seg_2.nii.gz") +# img_s = sitk.GetImageFromArray(np.squeeze(image)) +# # img_s.CopyInformation(seg) +# sitk.WriteImage(img_s, f"../test_image_3.nii.gz") +# img_s = sitk.GetImageFromArray(np.squeeze(coil_image)) +# # img_s.CopyInformation(seg) +# sitk.WriteImage(img_s, f"../test_last_coil_image_2.nii.gz") \ No newline at end of file diff --git a/scripts/miccai_rss_recon.py b/scripts/miccai_rss_recon.py new file mode 100755 index 0000000..77dd636 --- /dev/null +++ b/scripts/miccai_rss_recon.py @@ -0,0 +1,48 @@ +import numpy as np +import SimpleITK as sitk +from os import listdir +from pip import main +from scipy.fftpack import fftshift, ifftshift, ifftn +from umcglib.utils import apply_parallel +import time + +mypath = f'/data/pca-rad/datasets/miccai_2022/K2S_MICCAI2022_GRP/train/data/TBrecon1/train/untarred' +files = [f for f in listdir(mypath)] +# save_file = files[1] +# seg = np.load(f'{mypath}/{save_file}/{save_file}_seg.npy') +# kspace = np.load(f'{mypath}/{save_file}/{save_file}_ksp.npy') + +def k2p_to_img(kspace_single_coil): + kspace_padded = np.pad(np.squeeze(kspace_single_coil[:,:,:]),[(128,128),(128,128),(0,0)],'constant') + coil_image = np.abs(fftshift(fftshift(ifftn(ifftshift(kspace_padded)),axes=0),axes=1)) + return coil_image + +def img_from_ksp( + kspace: np.array): + kspace = np.transpose(kspace, (2, 0, 1, 3)) + coil_images = apply_parallel(kspace, k2p_to_img, 12) + image = np.round(np.sqrt(np.sum(np.power(coil_images,2),0))) + image = image[:,:,2:-2] + image = image[:,:,::-1] + return image + +for file_idx in range(167,300): + image_recon = [] + kspace = [] + save_file = [] + save_file = files[file_idx] + kspace = np.load(f'{mypath}/{save_file}/{save_file}_ksp.npy') + image_recon = img_from_ksp(kspace) + np.save(f'{mypath}/{save_file}/{save_file}_rss_recon.npy',image_recon) + print('done:',file_idx,' ',f'{mypath}/{save_file}/{save_file}_rss_recon.npy') + + +# quit() +# seg = sitk.GetImageFromArray(np.squeeze(seg)) +# sitk.WriteImage(seg, f"../test_seg_2.nii.gz") +# img_s = sitk.GetImageFromArray(np.squeeze(image)) +# # img_s.CopyInformation(seg) +# sitk.WriteImage(img_s, f"../test_image_3.nii.gz") +# img_s = sitk.GetImageFromArray(np.squeeze(coil_image)) +# # img_s.CopyInformation(seg) +# sitk.WriteImage(img_s, f"../test_last_coil_image_2.nii.gz") \ No newline at end of file diff --git a/scripts/miccai_rss_recon_undersampled.py b/scripts/miccai_rss_recon_undersampled.py new file mode 100755 index 0000000..7b16244 --- /dev/null +++ b/scripts/miccai_rss_recon_undersampled.py @@ -0,0 +1,66 @@ +import numpy as np +import SimpleITK as sitk +from os import listdir +from pip import main +from scipy.fftpack import fftshift, ifftshift, ifftn +from umcglib.utils import apply_parallel +import time +import h5py + +mypath = f'/data/pca-rad/datasets/miccai_2022/K2S_MICCAI2022_GRP/train/data/TBrecon1/train/untarred' +files = [f for f in listdir(mypath)] +KSPACE_U_PATH = f'/data/pca-rad/datasets/miccai_2022/K2S_MICCAI2022_GRP/train/data/TBrecon1/train/R_8_mask.h5' +f = h5py.File(KSPACE_U_PATH,'r') +kspace_mask = f['us_mask'][()] + +def k2p_to_img(kspace_single_coil): + kspace_padded = np.pad(np.squeeze(kspace_single_coil[:,:,:]),[(128,128),(128,128),(0,0)],'constant') + coil_image = np.abs(fftshift(fftshift(ifftn(ifftshift(kspace_padded)),axes=0),axes=1)) + return coil_image + +def img_from_ksp( + kspace: np.array): + kspace_masked = np.multiply(kspace, kspace_mask) + kspace_masked = np.transpose(kspace_masked, (2, 0, 1, 3)) + coil_images = apply_parallel(kspace_masked, k2p_to_img, 12) + image = np.round(np.sqrt(np.sum(np.power(coil_images,2),0))) + image = image[:,:,2:-2] + image = image[:,:,::-1] + return image.astype(np.float32) + +for file_idx in range(187,300): + image_recon = [] + kspace = [] + save_file = [] + save_file = files[file_idx] + kspace = np.load(f'{mypath}/{save_file}/{save_file}_ksp.npy') + image_recon = img_from_ksp(kspace) + np.save(f'{mypath}/{save_file}/{save_file}_rss_recon_usampled.npy',image_recon) + print('done:',file_idx,' ',f'{mypath}/{save_file}/{save_file}_rss_recon_usampled.npy') + + # img_s = sitk.GetImageFromArray(np.squeeze(image_recon)) + # sitk.WriteImage(img_s, f"../job_test_undersampled_{save_file}.nii.gz") + + + + + +# print(f.keys()) +# print(np.shape(f['us_mask'][()])) +# print(type(f['us_mask'][()])) + + + + + + + +# quit() +# seg = sitk.GetImageFromArray(np.squeeze(seg)) +# sitk.WriteImage(seg, f"../test_seg_2.nii.gz") +# img_s = sitk.GetImageFromArray(np.squeeze(image)) +# # img_s.CopyInformation(seg) +# sitk.WriteImage(img_s, f"../test_image_3.nii.gz") +# img_s = sitk.GetImageFromArray(np.squeeze(coil_image)) +# # img_s.CopyInformation(seg) +# sitk.WriteImage(img_s, f"../test_last_coil_image_2.nii.gz") \ No newline at end of file diff --git a/scripts/test.py b/scripts/test.py new file mode 100755 index 0000000..344dd7f --- /dev/null +++ b/scripts/test.py @@ -0,0 +1,94 @@ +from multiprocessing import set_start_method +import numpy as np +from umcglib.utils import apply_parallel,print_stats_np +from os import listdir +import matplotlib.pyplot as plt + +# mypath = f'/data/pca-rad/datasets/miccai_2022/K2S_MICCAI2022_GRP/train/data/TBrecon1/train/untarred' +# patient_id = 'pat00123' +# seg = np.load(f'{mypath}/{patient_id}/{patient_id}_seg.npy').astype(int) +# print(np.shape(seg)) +# print_stats_np(seg,'seg_sdgsd') + +def load_seg(file,mypath): + patient_id = file + print(patient_id) + seg = np.load(f'{mypath}/{patient_id}/{patient_id}_seg.npy').astype(int) + background = np.sum(seg == 0) + femoral_cartilage = np.sum(seg == 1) + tibial_cartilage = np.sum(seg == 2) + patellar_cartilage = np.sum(seg == 3) + femur = np.sum(seg == 4) + tibia = np.sum(seg == 5) + patella = np.sum(seg == 6) + output = [background,femoral_cartilage,tibial_cartilage,patellar_cartilage,femur,tibia,patella] + # print('background',background) + return output + +def calculate_hist(segs): + background = np.sum(segs == 0) + femoral_cartilage = np.sum(segs == 1) + tibial_cartilage = np.sum(segs == 2) + patellar_cartilage = np.sum(segs == 3) + femur = np.sum(segs == 4) + tibia = np.sum(segs == 5) + patella = np.sum(segs == 6) + output = [background,femoral_cartilage,tibial_cartilage,patellar_cartilage,femur,tibia,patella] + # print('background',background) + return output + +if __name__ == '__main__': + path = f'/data/pca-rad/datasets/miccai_2022/K2S_MICCAI2022_GRP/train/data/TBrecon1/train/untarred' + files = [f for f in listdir(path)] + print('files length',len(files)) + # For multiprocessing + set_start_method("spawn") + print('spawn done') + # segs_list = apply_parallel(files, load_seg, 4, mypath = path) + segs_list = apply_parallel(files, load_seg, 12, mypath = path) + print('done loading 50') + # segs_list_100 = apply_parallel(files[50:100], load_seg, 4, mypath = path) + # print('done loading 100') + # segs_list_150 = apply_parallel(files[100:150], load_seg, 4, mypath = path) + # print('done loading 150') + # segs_list_200 = apply_parallel(files[150:200], load_seg, 4, mypath = path) + # print('done loading 200') + # segs_list_250 = apply_parallel(files[200:250], load_seg, 4, mypath = path) + # print('done loading 250') + # segs_list_300 = apply_parallel(files[250:300], load_seg, 4, mypath = path) + print('done loading 300') + # segs_list.extend(segs_list_100) + # segs_list.extend(segs_list_150) + # segs_list.extend(segs_list_200) + # segs_list.extend(segs_list_250) + # segs_list.extend(segs_list_300) + print('segs_list is done') + + output = np.stack(segs_list, axis=0) + # print('load done',np.shape(segs_n)) + + # output = apply_parallel(segs_list,calculate_hist,4) + print('hist calc done',np.shape(output)) + # print(np.stack(output)[:,0]) + + dict = { + 'background': np.stack(output)[:,0], + 'femoral_cartilage': np.stack(output)[:,1], + 'tibial_cartilage': np.stack(output)[:,2], + 'patellar_cartilage': np.stack(output)[:,3], + 'femur': np.stack(output)[:,4], + 'tibia': np.stack(output)[:,5], + 'patella': np.stack(output)[:,6] + } + + np.save('../dict.npy',dict) + + for keys in dict: + plt.figure() + data = dict[f'{keys}'] + plt.hist(data) + plt.title(f"historgram of {keys}") + plt.savefig(f"../{keys}.png", dpi=300) + print('done ',keys) + + diff --git a/scripts/tmp.py b/scripts/tmp.py index caccff1..26b7c44 100755 --- a/scripts/tmp.py +++ b/scripts/tmp.py @@ -1,27 +1,178 @@ -from scipy import stats +from traceback import FrameSummary import numpy as np -# siemens_froc = [1.68,1.81,1.44,1.55] -# b400_froc = [3.4,3.93,2.82,] -# b800_froc = [1.58,1.99,1.36,1.6] - -# siemens_roc = [0.782, 0.732, 0.775, 0.854] -b400_roc = [0.746, 0.814, 0.789, 0.763] -b800_roc = [0.786, 0.731, 0.67, 0.782] - -# stat_test = stats.wilcoxon(siemens_froc,b800_froc,alternative='less') -# print('froc stats siemens > b400',stat_test) -# print(' Mean and std siemens froc:', np.mean(siemens_froc),'+-',np.std(siemens_froc)) -# print(' Mean and std b400 froc:', np.mean(b400_froc),'+-',np.std(b400_froc)) -# print(' Mean and std b800 froc:', np.mean(b800_froc),'+-',np.std(b800_froc)) - -# print(' Mean and std siemens roc:', np.mean(siemens_roc),'+-',np.std(siemens_roc)) -print(' Mean and std b400 roc:', np.mean(b400_roc),'+-',np.std(b400_roc)) -print(' Mean and std b800 roc:', np.mean(b800_roc),'+-',np.std(b800_roc)) +import SimpleITK as sitk +from os import listdir +from pip import main +from scipy.fftpack import fftshift, ifftshift, ifftn +from umcglib.utils import apply_parallel, print_stats_np +import time +import h5py +import matplotlib.pyplot as plt +from multiprocessing import set_start_method +import gzip +import os -# The test has been introduced in [4]. Given n independent samples (xi, yi) from a bivariate distribution -# (i.e. paired samples), it computes the differences di = xi - yi. One assumption of the test is that the -# differences are symmetric, see [2]. The two-sided test has the null hypothesis that the median of the -# differences is zero against the alternative that it is different from zero. The one-sided test has the -# null hypothesis that the median is positive against the alternative that it is negative -# (alternative == 'less'), or vice versa (alternative == 'greater.'). \ No newline at end of file +mypath = f'/data/pca-rad/datasets/miccai_2022/K2S_MICCAI2022_GRP/train/data/TBrecon1/train/untarred' +files = [f for f in listdir(mypath)] + +for file_idx in range(188): + image_recon = [] + save_file = [] + save_file = files[file_idx] + image_recon = np.load(f'{mypath}/{save_file}/{save_file}_rss_recon_usampled.npy').astype(np.float32) + np.save(f'{mypath}/{save_file}/{save_file}_rss_recon_usampled.npy',image_recon) + print('done:',file_idx,' ',f'{mypath}/{save_file}/{save_file}_rss_recon_usampled.npy') + + + + + + + + + +quit() +mypath = f'/data/pca-rad/datasets/miccai_2022/K2S_MICCAI2022_GRP/train/data/TBrecon1/train/untarred' +patient_id = 'pat00335' +seg = np.load(f'{mypath}/{patient_id}/{patient_id}_seg.npy') +print_stats_np(seg,'seg_') +image = np.load(f'{mypath}/{patient_id}/{patient_id}_rss_recon.npy') +print_stats_np(image,'image_') +image_u = np.load(f'{mypath}/{patient_id}/{patient_id}_rss_recon_usampled.npy') +print_stats_np(image_u,'image_u_') + +seg = sitk.GetImageFromArray(np.squeeze(seg)) +sitk.WriteImage(seg, f"../test_seg_335.nii.gz") +image = sitk.GetImageFromArray(np.squeeze(image)) +sitk.WriteImage(image, f"../test_image_335.nii.gz") +image_u = sitk.GetImageFromArray(np.squeeze(image_u)) +sitk.WriteImage(image_u, f"../test_image_u_335.nii.gz") +print(np.shape(seg)) + + + + +# ROOT_H5_GZ = r'../../../../../scratch/p290820/MICCAI/TBrecon-01-02-00335.h5.gz' +# SAVE_TO_DIR = r'/data/pca-rad/datasets/miccai_2022/K2S_MICCAI2022_GRP/train/data/TBrecon1/train/untarred/pat00335' +# unzipped = gzip.open(ROOT_H5_GZ,'rb') +# h5_file = h5py.File(unzipped,'r') # Keys: kspace, kspace_info, seg, seg_info + +# # ksp_n = h5_file['kspace'][()] #np.complex64 +# seg_n = h5_file['seg'][()].astype(np.float32) #np.float32 + +# # fname_ksp = os.path.join(patient_path, f"pat{pat_num}_ksp.npy") +# fname_seg = os.path.join(SAVE_TO_DIR, f"pat00335_seg_test.npy") +# # np.save(fname_ksp, ksp_n) +# np.save(fname_seg, seg_n) + + + +quit() + +def inladen(patient_id,mypath): + seg = np.load(f'{mypath}/{patient_id}/{patient_id}_seg.npy').astype(int) + return seg + + +def load_seg(file,mypath): + patient_id = file + print(patient_id) + seg = inladen(patient_id,mypath) + return seg + +def calculate_hist(segs): + background = np.sum(segs == 0) + femoral_cartilage = np.sum(segs == 1) + tibial_cartilage = np.sum(segs == 2) + patellar_cartilage = np.sum(segs == 3) + femur = np.sum(segs == 4) + tibia = np.sum(segs == 5) + patella = np.sum(segs == 6) + output = [background,femoral_cartilage,tibial_cartilage,patellar_cartilage,femur,tibia,patella] + # print('background',background) + return output + +if __name__ == '__main__': + path = f'/data/pca-rad/datasets/miccai_2022/K2S_MICCAI2022_GRP/train/data/TBrecon1/train/untarred' + files = [f for f in listdir(path)] + print('files length',len(files)) + # For multiprocessing + set_start_method("spawn") + print('spawn done') + segs_list = apply_parallel(files[:100], load_seg, 4, mypath = path) + print('done loading 100') + + segs_list_200 = apply_parallel(files[100:200], load_seg, 4, mypath = path) + print('done loading 200') + + segs_list_300 = apply_parallel(files[200:300], load_seg, 4, mypath = path) + print('done loading 300') + + segs_list.extend(segs_list_200) + segs_list.extend(segs_list_300) + + segs_n = np.stack(segs_list, axis=0) + print('load done',np.shape(segs_n)) + + output = apply_parallel(segs_list,calculate_hist,4) + print('hist calc done',np.shape(np.stack(output))) + # print(np.stack(output)[:,0]) + + dict = { + 'background': np.stack(output)[:,0], + 'femoral_cartilage': np.stack(output)[:,0], + 'tibial_cartilage': np.stack(output)[:,0], + 'patellar_cartilage': np.stack(output)[:,0], + 'femur': np.stack(output)[:,0], + 'tibia': np.stack(output)[:,0], + 'patella': np.stack(output)[:,0] + } + + # for file_idx in range(300): + # seg = [] + # patient_id = files[file_idx] + # seg = np.load(f'{mypath}/{patient_id}/{patient_id}_seg.npy') + # # # 0: background; 1: femoral cartilage; 2: tibial cartilage; 3: patellar cartilage; 4: femur; 5: tibia; 6: patella. + + + # dict['background'].append(np.sum(seg == 0)) + # dict['femoral_cartilage'].append(np.sum(seg == 1)) + # dict['tibial_cartilage'].append(np.sum(seg == 2)) + # dict['patellar_cartilage'].append(np.sum(seg == 3)) + # dict['femur'].append(np.sum(seg == 4)) + # dict['tibia'].append(np.sum(seg == 5)) + # dict['patella'].append(np.sum(seg == 6)) + + # print('done:',file_idx,' ',f'{mypath}/{patient_id}/{patient_id}_seg.npy') + + np.save('../dict.npy',dict) + + for keys in dict: + data = dict[f'{keys}'] + plt.hist(data) + plt.title(f"historgram of {keys}") + plt.savefig(f"../{keys}.png", dpi=300) + print('done ',keys) + + + + # print(f.keys()) + # print(np.shape(f['us_mask'][()])) + # print(type(f['us_mask'][()])) + + + + + + + + # quit() + # seg = sitk.GetImageFromArray(np.squeeze(seg)) + # sitk.WriteImage(seg, f"../test_seg_2.nii.gz") + # img_s = sitk.GetImageFromArray(np.squeeze(image)) + # # img_s.CopyInformation(seg) + # sitk.WriteImage(img_s, f"../test_image_3.nii.gz") + # img_s = sitk.GetImageFromArray(np.squeeze(coil_image)) + # # img_s.CopyInformation(seg) + # sitk.WriteImage(img_s, f"../test_last_coil_image_2.nii.gz") \ No newline at end of file