-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathInside.out_SC.dilation_MS.py
85 lines (65 loc) · 4 KB
/
Inside.out_SC.dilation_MS.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
# -*- coding: utf-8 -*-
"""
The purpose of this script is to compute the average metrix within a mask. The mask is growing outward at each iteration.
Algorithm:
- Generate spinal cord segmentation --> SCseg
- Manually create a 2-pixel mask at the center of the spinal canal.
- For each new iteration:
- The size of the mask increases by one pixel, using 2D dilation.
- We compute the intersection between the dilated mask and the SCseg --> MaskIntersect
- We compute the average voxel values within MaskIntersect
At each step, the mask is dilated by a ball of radius 1, then we substract
the previous step to get only a ring.
We then average the field over the ring
"""
from numpy import *
import os
import nibabel as nib
from skimage.morphology import binary_dilation, ball
###############################################################################
''' DATA LOAD '''
# To get speed and memory, data need to be unzipped format like *.nii, not .nii.gz
# a mgz file can be renamed mgh.gz and then unziped.
subject_list = ['WM.GM.test_MS073_R01_05_inf']
#['MS070_R01_02_inf','MS070_R01_02_sup','MS073_R01_05_inf','MS073_R01_05_sup','MS075_R01_07_inf','MS075_R01_07_sup','MS079_R01_11_inf','MS079_R01_11_sup','MS081_R01_13_inf','MS081_R01_13_sup','MS084_R01_16_inf','MS084_R01_16_sup','MS087_R01_19_inf','MS087_R01_19_sup','MS088_R01_20_inf','MS088_R01_20_mid','MS088_R01_20_sup','MS089_R01_29_inf','MS089_R01_29_sup','MS091_R01_22_inf','MS091_R01_22_sup','MS093_R01_24_inf','MS093_R01_24_sup','MS095_R01_26_inf','MS095_R01_26_sup','MS096_R01_27_inf','MS096_R01_27_sup','MS098_R01_30_inf','MS098_R01_30_sup','MS099_R01_31_inf','MS099_R01_31_sup','MS100_R01_32_inf','MS100_R01_32_sup','MS101_R01_33_inf','MS101_R01_33_sup','MS102_R01_34_inf','MS102_R01_34_sup','MS103_R01_35_inf','MS103_R01_35_sup','MS104_R01_36_inf','MS104_R01_36_sup','NIMMS01_inf','NIMMS01_mid','NIMMS01_sup','NIMMS02_inf','NIMMS02_sup','NIMMS03_inf','NIMMS03_sup','NIMMS04_inf','NIMMS04_sup','NIMMS05_inf','NIMMS05_sup','NIMMS06_inf','NIMMS06_sup','NIMMS07_inf','NIMMS07_sup','NIMMS08_inf','NIMMS08_sup','NIMMS09_inf','NIMMS09_sup','NIMMS10_inf','NIMMS10_sup','NIMMS11_inf','NIMMS11_sup']
print("Load SUVs")
for k in range(len(subject_list)):
subject = subject_list[k]
#recon = recon_list[k]
#subject_dir = "/autofs/cluster/mscat/users/sindhuja/data/3T/"+recon+"/mri"
path_files = "/Users/raouelletteiv/Research/7T_Spinal_Brain_Study/Data/Spinal_data/Spinal.Rings_outside_in/"+subject
SUV_filename = 't2s_seg_manual.nii.gz'
SUV = nib.load(os.path.join(path_files,SUV_filename)).get_data()
SUV_pos = SUV>0
print("Load region")
NAWM_filename = 't2s_seg_manual.nii.gz'
region = nib.load(os.path.join(path_files,NAWM_filename)).get_data() > 0
print("Load mask")
#LV_filename = 'LV.nii'
LV_filename = 't2s_centerline_expand.nii.gz'
raw_mask = nib.load(os.path.join(path_files,LV_filename)).get_data() >0
#
#
mask = raw_mask
del raw_mask
#
#
''' DILATIONS '''
print("dilations")
NDIL = 25 # Number of dilations
annot = nib.load(os.path.join(path_files,NAWM_filename))
annot.get_data()[:,:,:]=0
res = []
for ndil in range(NDIL):
mask_dil = binary_dilation(mask,ball(1)).astype('bool')
ring = (mask_dil ^ mask).astype('bool') & region
ring = ring & SUV_pos
av = sum(SUV[ring])/sum(ring)
mask = mask_dil.copy()
annot.get_data()[ring]=ndil+1
res = res + [[ndil,av]]
print(ndil, av)
annot_name = "/Users/raouelletteiv/Research/7T_Spinal_Brain_Study/Data/Spinal_data/Spinal.Rings_outside_in/"+subject+"/"+subject+"_INsideOut_2vox_dilation.nii"
nib.save(annot,annot_name) #to write the file with the dilations
outfilename = "/Users/raouelletteiv/Research/7T_Spinal_Brain_Study/Data/Spinal_data/Spinal.Rings_outside_in/"+subject+"/"+subject+"_INsideOut_2vox_dilation.csv"
savetxt(outfilename,array(res),delimiter=' ',header="ndil;average")