demo script

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