import SimpleITK as sitk
import numpy as np
import matplotlib.pyplot as plt

img = sitk.ReadImage("Thorax1.0B40f_Downsampled_by_2.nrrd")

# Abstand zwischen Voxelmittelpunkten
xSp, ySp, zSp = img.GetSpacing()
volumePerVoxel = xSp * ySp * zSp

print(f"1) Die Voxel sind Quader mit den verschiedenen Seitenlängen: X: {xSp}, Y: {ySp}, Z: {zSp} \n")

# Menge an Voxeln in eine Richtung
sizeVector = img.GetSize()
xSi, ySi, zSi = sizeVector

print(f"2) Das Bild hat in den drei Raumrichtungen {xSi * ySi * zSi} Voxel. \n")
print(f"3) Die Methode GetSize() gibt den Typ tuple zurück. \n")
print(f"4) Der Bildausschnitt in die drei Raumrichtungen ist X: {xSp*xSi} mm, Y: {ySp*ySi} mm, Z: {zSp*zSi} mm\n")

# CT-Bild der Lunge erstellen
seedList = [(120,140,100)]
lungs = sitk.ConnectedThreshold(img, seedList , -1100, -200)
sitk.WriteImage(lungs, "lungs.nrrd")

# Volumen der Lunge berechnen
lungsArray = np.array(sitk.GetArrayFromImage(lungs))
voxelInLungs = lungsArray[lungsArray > 0].size * volumePerVoxel
lungsInLitre = voxelInLungs / 1000000
#print(f"{lungsInLitre} L \n")

# Histogramm des CT-Bilds
thoraxArray = np.array(sitk.GetArrayFromImage(img))
thoraxArray1D = thoraxArray.flatten()
plt.hist(thoraxArray1D, bins=100)
plt.title("Helligkeitswerte eines Thorax")
plt.xlabel("Helligkeitswert in Hounsfield-Einheiten (HU)")
plt.ylabel("Häufigkeit in Millionen (mio)")
plt.show()
plt.clf()

# Alle Helligkeitswerte der Lunge separieren
lungsArray1D = lungsArray.flatten()
lungBrightness = thoraxArray1D[lungsArray1D > 0]

# Histogram aller Helligkeitswerte in der Lunge 
plt.hist(lungBrightness, bins=100)
plt.title("Helligkeitswerte einer Lunge")
plt.xlabel("Helligkeitswert in Hounsfield-Einheiten (HU)")
plt.ylabel("Häufigkeit")
plt.show()
plt.clf()

print("5) In der Lunge treten die Helligkeitswerte ~ -1024 bis -200 auf.")
print("   Bei ~ -900 ist ein peak an Helligkeit zu sehen. \n")