

In [1]:

    
## Boiler plate code common to many notebooks.  See the TestFilesCommonCode.ipynb for details
from __future__ import print_function
%run TestFilesCommonCode.ipynb









    



SimpleITK Version: 0.9.0.dev2151-g2a716
Compiled: Dec 23 2014 21:10:53



In [3]:

    
FA_filename='/scratch/TESTS/BS/DWIWorkflow_Nipype/6_RunOnACroppedFile/19568/Outputs/FA.nrrd'
t2_filename = '/Shared/sinapse/CACHE/20141001_PREDICTHD_long_Results/PHD_024/1155/19568/TissueClassify/t2_average_BRAINSABC.nii.gz'
labels_filename = '/Shared/sinapse/CACHE/20141001_PREDICTHD_long_Results/PHD_024/1155/19568/TissueClassify/complete_brainlabels_seg.nii.gz'

FA = sitk.ReadImage(FA_filename)
T2 = sitk.ReadImage(t2_filename)
labels = sitk.ReadImage(labels_filename)

myshow(T2)
myshow(FA)
myshow(sitk.LabelToRGB(labels))



In [6]:

    
size=labels.GetSize()
myshow(sitk.LabelToRGB(labels[size[0]//2,:,::-1]))
size=FA.GetSize()
myshow(sitk.Expand(FA[size[0]//2,:,:],[5,5,5]))

FA and label maps are in the same physical space, but they have different "sizes". Also, FA has some missing brain regions, e.g. it does not cover the CSF from the label map.



In [53]:

    
# Resample label map to FA image space
resFilt = sitk.ResampleImageFilter()
resFilt.SetReferenceImage(FA)
resFilt.SetOutputPixelType(labels.GetPixelIDValue())
resFilt.SetInterpolator(1) # Nearest neighbour interpolation
labels_res = resFilt.Execute(labels)

size=labels_res.GetSize()
myshow(sitk.Expand(sitk.LabelToRGB(labels_res[size[0]//2,:,:]),[5,5,5]))

Now compute the mean of FA on region of interest:

WM -> 1

SURFGM -> 2

CSF -> 4

VB -> 5

CRBLGM -> 11

CRBLWM -> 12

BASAL -> 19

GLOBUS -> 23

THALAMUS -> 24

HIPPOCAMPUS -> 25

AIR -> 0



In [110]:

    
labelID = 24



In [111]:

    
statFilter = sitk.LabelStatisticsImageFilter()
statFilter.Execute(FA, labels_res)

print('mean:',statFilter.GetMean(labelID))
print('std:',statFilter.GetSigma(labelID))
print('max:',statFilter.GetMaximum(labelID))
print('min:',statFilter.GetMinimum(labelID))
print('median:',statFilter.GetMedian(labelID))
print('Num of voxels:',statFilter.GetCount(labelID))

voxelSize=FA.GetSpacing()[0]*FA.GetSpacing()[1]*FA.GetSpacing()[2]
effective_volume=statFilter.GetCount(labelID) * voxelSize
print('effective_volume:',effective_volume)

#print(statFilter.GetNumberOfLabels())









    



mean: 0.3125944811
std: 0.196391885731
max: 0.887071133766
min: 0.0
median: 0.282317265868
Num of voxels: 3943
effective_volume: 31544.0

FA spacing is 2mm, so the volume of each voxel in FA lattice is 8 mm^3

CSF_mean = 0.056

effective_CSF_Volume = 20430 x 8 = 163440 mm^3

WM_mean = 0.278

effective_WM_Volume = 71550 x 8 = 572400 mm^3

SURFGM_mean = 0.094

effective_SURFGM_Volume = 42086 x 8 = 336688 mm^3

effective_GLOBUS_Volume = 1709 x 8 = 13672 mm^3

effective_THALAMUS_Volume = 3943 x 8 = 31544 mm^3

effective_HIPPOCAMPUS_Volume = 2399 x 8 = 19192 mm^3



In [112]:

    
# Now resample FA to label map image space
resFilt = sitk.ResampleImageFilter()
resFilt.SetReferenceImage(labels)
FAres = resFilt.Execute(FA)

size=FAres.GetSize()
myshow(sitk.Expand(FAres[size[0]//2,:,::-1],[1,1,1]))



In [113]:

    
statFilter = sitk.LabelStatisticsImageFilter()
statFilter.Execute(FAres, labels)

print('mean:',statFilter.GetMean(labelID))
print('std:',statFilter.GetSigma(labelID))
print('Num of voxels:',statFilter.GetCount(labelID))

voxelSize=labels.GetSpacing()[0]*labels.GetSpacing()[1]*labels.GetSpacing()[2]
total_volume=statFilter.GetCount(labelID) * voxelSize
print('total_volume:',total_volume)
print('confidence_coeficient:',effective_volume/total_volume)

#print(statFilter.GetNumberOfLabels())









    



mean: 0.308409052079
std: 0.177767487138
Num of voxels: 31514
total_volume: 31514.0
confidence_coeficient: 1.00095195786

Label map spacing is 1mm, so the volume of each voxel in Label map lattice is 1 mm^3

total_CSF_Volume = 186805 x 1 = 186805 mm^3

CSF_confidence_coeficient = effective_CSF_Volume/total_CSF_Volume = 163440/186805 = 0.874

total_WM_Volume = 576212 x 1 = 576212 mm^3

WM_confidence_coeficient = effective_WM_Volume/total_WM_Volume = 572400/576212 = 0.993

total_SURFGM_Volume = 340779 x 1 = 340779 mm^3

SURFGM_confidence_coeficient =

effective_SURFGM_Volume/total_SURFGM_Volume = 336688/340779 = 0.987

For GLOBUS, THALAMUS and HIPPOCAMPUS, the confidence coeficient should be 1 because there is no missing part for these regions, and for each of the, the whole volume has contributed in computing the average FA. As expected computed confidence coeficient is close to 1 for these regions, and the mismatch has happened because of resampling error when label map is resampled to FA image space

total_GLOBUS_Volume = 13701 x 1 = 13701 mm^3

GLOBUS_confidence_coeficient =

effective_GLOBUS_Volume/total_GLOBUS_Volume = 13672/13701 = 0.997

total_THALAMUS_Volume = 31514 x 1 = 31514 mm^3

THALAMUS_confidence_coeficient =

effective_THALAMUS_Volume/total_THALAMUS_Volume = 31544/31514 = 1.0009

total_HIPPOCAMPUS_Volume = 19427 x 1 = 19427 mm^3

HIPPOCAMPUS_confidence_coeficient =

effective_HIPPOCAMPUS_Volume/total_HIPPOCAMPUS_Volume = 19192/19427 = 0.987



In [ ]: