2022-04-07 09:13:17 +02:00
|
|
|
import numpy as np
|
|
|
|
import SimpleITK as sitk
|
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
######## load images #############
|
2022-04-21 14:29:36 +02:00
|
|
|
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_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_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)
|
|
|
|
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)
|
2022-04-07 09:13:17 +02:00
|
|
|
|
|
|
|
def show_img(greyscale_img):
|
|
|
|
fig = plt.figure()
|
|
|
|
plt.imshow(greyscale_img)
|
|
|
|
plt.axis('on')
|
|
|
|
path = f"iets.png"
|
|
|
|
fig.savefig(path, dpi=300, bbox_inches='tight')
|
|
|
|
|
|
|
|
def calc_adc(b50, b400, b800):
|
|
|
|
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)
|
2022-04-21 14:29:36 +02:00
|
|
|
adc_with_zeros = np.divide(denominator, numerator) * -1000000
|
|
|
|
adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros))
|
|
|
|
return adc
|
2022-04-07 09:13:17 +02:00
|
|
|
|
|
|
|
def calc_adc_1(b50,b800):
|
|
|
|
mean_dwi = (50 + 800) / 2
|
|
|
|
mean_si = np.divide(np.add(np.log(b50), 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)
|
2022-04-21 14:29:36 +02:00
|
|
|
adc_with_zeros = np.divide(denominator, numerator) * -1000000
|
|
|
|
adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros))
|
|
|
|
return adc
|
2022-04-07 09:13:17 +02:00
|
|
|
|
|
|
|
def calc_adc_2(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)
|
2022-04-21 14:29:36 +02:00
|
|
|
adc_with_zeros = np.divide(denominator, numerator) * -1000000
|
|
|
|
adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros))
|
|
|
|
return adc
|
2022-04-07 09:13:17 +02:00
|
|
|
|
|
|
|
def calc_adc_3(b400,b800):
|
|
|
|
mean_dwi = (400 + 800) / 2
|
|
|
|
mean_si = np.divide(np.add(np.log(b400), np.log(b800)), 2)
|
|
|
|
|
|
|
|
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)
|
2022-04-21 14:29:36 +02:00
|
|
|
adc_with_zeros = np.divide(denominator, numerator) * -1000000
|
|
|
|
adc = np.clip(adc_with_zeros,0,np.amax(adc_with_zeros))
|
|
|
|
return adc
|
2022-04-07 09:13:17 +02:00
|
|
|
|
|
|
|
def calc_high_b(b_value_high,b_value,b_image,ADC_map):
|
2022-04-21 14:29:36 +02:00
|
|
|
high_b = np.multiply(b_image, np.exp(np.multiply(np.subtract(b_value,b_value_high), (np.divide(ADC_map,1000000)))))
|
2022-04-07 09:13:17 +02:00
|
|
|
return high_b
|
|
|
|
|
2022-04-21 14:29:36 +02:00
|
|
|
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_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_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)
|
|
|
|
# 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")
|
|
|
|
|
|
|
|
##########################################################################################
|
2022-04-07 09:13:17 +02:00
|
|
|
# b50 = sitk.GetImageFromArray(b50)
|
|
|
|
# sitk.WriteImage(b50, "./../train_output/adc_exp/b50_true.nii.gz")
|
|
|
|
|
|
|
|
# b400 = sitk.GetImageFromArray(b400)
|
|
|
|
# sitk.WriteImage(b400, "./../train_output/adc_exp/b400_true.nii.gz")
|
|
|
|
|
|
|
|
# b800 = sitk.GetImageFromArray(b800)
|
|
|
|
# sitk.WriteImage(b800, "./../train_output/adc_exp/b800_true.nii.gz")
|
|
|
|
|
|
|
|
# b1400 = sitk.GetImageFromArray(b1400)
|
|
|
|
# sitk.WriteImage(b1400,f"./../train_output/adc_exp/b1400_true.nii.gz")
|
|
|
|
|
|
|
|
# adc = sitk.GetImageFromArray(adc)
|
2022-04-21 14:29:36 +02:00
|
|
|
# sitk.WriteImage(adc, f"adc_true.nii.gz")
|