diff --git a/src/test1.py b/src/test1.py new file mode 100644 index 0000000..f670f7d --- /dev/null +++ b/src/test1.py @@ -0,0 +1,103 @@ +import sys + +from PIL import Image +from scipy.stats import chi2_contingency + +import numpy as np +import SimpleITK as sitk + + +CT = '/mnt/1218/Public/dataset/Spine/CTSpine1K/data/liver/liver_169.nii.gz' +CT = '/mnt/1218/Public/dataset/Spine/CTSpine1K/data/liver/liver_167.nii.gz' +OUT = '/mnt/1248/xfr/Downloads/out.nii.gz' + +WINDOW=1800 +LEVEL=400 + +def rescale_and_write(image_array, output_image, new_size=(512,512)): + max_val = np.max(image_array) + min_val = np.min(image_array) + scaled_array = ((image_array - min_val) / (max_val - min_val) * 255).astype(np.uint8) + im = Image.fromarray(scaled_array).resize(new_size) + im.save(output_image) + +def rescale_cast(image): + a = sitk.RescaleIntensity(image) + return sitk.Cast(a, sitk.sitkUInt8) + + +img = sitk.ReadImage(CT) +img = sitk.Clamp(img, sitk.sitkFloat32, LEVEL-WINDOW/2, LEVEL+WINDOW/2) +img = sitk.RescaleIntensity(img) +# img = sitk.RescaleIntensity(img, 0, 65535) + + +size = img.GetSize() + +p0 = sitk.MeanProjection(img, 0)[0,:,:] +rescale_and_write(sitk.GetArrayFromImage(p0), "/mnt/1248/xfr/Downloads/p0.jpg") + +p1 = sitk.MeanProjection(img, 1)[:,0,:] +rescale_and_write(sitk.GetArrayFromImage(p1), "/mnt/1248/xfr/Downloads/p1.jpg") + +# too slow +# b0 = [p0] * size[0] +# b0 = sitk.JoinSeries(b0) +# sitk.WriteImage(b0, "/mnt/1248/xfr/Downloads/b0.nii.gz") + +arr = sitk.GetArrayFromImage(img) + +# Compute means +means = arr.mean(axis=1, keepdims=True) +# Broadcast to full shape +result1 = np.tile(means, (1, arr.shape[1], 1)) +b1 = sitk.GetImageFromArray(result1) +b1.CopyInformation(img) +# sitk.WriteImage(b1, "/mnt/1248/xfr/Downloads/b1.nii.gz") + +# Compute means +means = arr.mean(axis=2, keepdims=True) +# Broadcast to full shape +result2 = np.tile(means, (1, 1, arr.shape[2])) +b2 = sitk.GetImageFromArray(result2) +b2.CopyInformation(img) +# sitk.WriteImage(b2, "/mnt/1248/xfr/Downloads/b2.nii.gz") + +# b3=sitk.Multiply(b1, b2) +result3 = result1*result2 +b3 = sitk.GetImageFromArray(result3) +b3.CopyInformation(img) + +sitk.WriteImage(rescale_cast(img), OUT) +sitk.WriteImage(rescale_cast(b1), "/mnt/1248/xfr/Downloads/b1.nii.gz") +sitk.WriteImage(rescale_cast(b2), "/mnt/1248/xfr/Downloads/b2.nii.gz") +sitk.WriteImage(rescale_cast(b3), "/mnt/1248/xfr/Downloads/b3.nii.gz") + + +exit() + + + + + +array_zyx = sitk.GetArrayFromImage(img) +print(array_zyx.shape) + +slice = array_zyx[array_zyx.shape[0]//2]+1 +print(slice.shape) + +rescale_and_write(slice, "/mnt/1248/xfr/Downloads/im1.jpg") + + +# Example observed contingency table +observed_data = np.array([[30, 20], + [10, 40]]) + +observed_data = slice + +chi2_stat, p_value, dof, expected_counts = chi2_contingency(observed_data) + +# print("Expected Counts:") +# print(expected_counts) + +rescale_and_write(expected_counts, "/mnt/1248/xfr/Downloads/im2.jpg")