/
hurst_rs_script_final.py
144 lines (114 loc) · 5.33 KB
/
hurst_rs_script_final.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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#!/usr/bin/env python3
import sys
if len(sys.argv) > 1:
name = sys.argv[1]
else:
name = input("Enter name:")
print(name)
import multiprocessing
from joblib import Parallel, delayed
import os
import nibabel as nib
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
#import FracTool_Final
import nolds
#load sample slice and convert to np array
example_slice = name
slice_img = nib.load(example_slice, mmap=False)
#get header of image
n1_header = slice_img.header
#get TR from header
TR = n1_header.get_zooms()
TR = np.asarray(TR)
TR = TR[3]
slice_array = slice_img.get_fdata()
[N1, N2, slice_num, N3] = slice_array.shape
row = np.arange(0,N1)
column = np.arange(0,N2)
def FracTool_voxel(i, j):
'''This function returns the Hurst coefficient and class of each voxel in the slice'''
global TR
rawbold = (slice_sq[i,j])
result_fractool = FracTool_Final.FracTool(rawbold, TR) #run Fractool on each signal in signal array
H_val = result_fractool[1] #result[1] of FracTool is Hurst value
#uncomment Class_val if you want a 3D image of signal class
#Class_val = result[0] #result[0] of FracTool is Class
return result_fractool #result is array of 2 elements
def dfa_voxel(i,j):
rawbold = (slice_sq[i,j])
result_dfa = nolds.dfa(rawbold)
return result_dfa
def hurst_voxel(i,j):
rawbold = (slice_sq[i,j])
result_hurst = nolds.hurst_rs(rawbold)
return result_hurst
def corr_dim_voxel(i,j):
rawbold = (slice_sq[i,j])
result_corr_dim = nolds.corr_dim(rawbold, 1)
return result_corr_dim
def sampen_voxel(i,j):
rawbold = (slice_sq[i,j])
result_sampen = nolds.sampen(rawbold)
return result_sampen
#multiprocessing for loop
num_cores = multiprocessing.cpu_count()
#base will be the bottom of stacking 2D slice arrays on top to create 3D image
base_fractool = np.zeros((N1,N2))
base_dfa = np.zeros((N1,N2))
base_hurst = np.zeros((N1,N2))
base_corr_dim = np.zeros((N1,N2))
base_sampen = np.zeros((N1,N2))
#for loop iterates through every slice of brain and stacks slices on top of each other
for x in range(slice_num):
#change 4th index in slice_array to change the number of timepoints
#ex 0:700 is first 700, : is the whole signal
slice_sq = slice_array[:,:,x,:]
#output_fractool = Parallel(n_jobs=num_cores)(delayed(FracTool_voxel)(i,j) for i in row for j in column)
#output_fractool = np.array(output_fractool) #convert output into numpy array
#output_fractool = output_fractool.astype(np.float64) #convert type object list elements to type float
#Hurst_matrix_fractool = output_fractool[:,1] #second column of result output array of result is Hurst
#Hurst_matrix_fractool = Hurst_matrix_fractool.reshape(N1,N2)
#output_dfa = Parallel(n_jobs=num_cores)(delayed(dfa_voxel)(i,j) for i in row for j in column)
#output_dfa = np.array(output_dfa) #convert output into numpy array
#output_dfa = output_dfa.astype(np.float64) #convert type object list elements to type float
#Hurst_matrix_dfa = output_dfa.reshape(N1,N2)
output_hurst = Parallel(n_jobs=num_cores)(delayed(hurst_voxel)(i,j) for i in row for j in column)
output_hurst = np.array(output_hurst) #convert output into numpy array
output_hurst = output_hurst.astype(np.float64) #convert type object list elements to type float
Hurst_matrix_hurst = output_hurst.reshape(N1,N2)
#output_corr_dim = Parallel(n_jobs=num_cores)(delayed(corr_dim_voxel)(i,j) for i in row for j in column)
#output_corr_dim = np.array(output_corr_dim) #convert output into numpy array
#output_corr_dim = output_corr_dim.astype(np.float64) #convert type object list elements to type float
#Hurst_matrix_corr_dim = output_corr_dim.reshape(N1,N2)
#output_sampen = Parallel(n_jobs=num_cores)(delayed(sampen_voxel)(i,j) for i in row for j in column)
#output_sampen = np.array(output_sampen) #convert output into numpy array
#output_sampen = output_sampen.astype(np.float64) #convert type object list elements to type float
#Hurst_matrix_sampen = output_sampen.reshape(N1,N2)
#base_fractool = np.dstack((base_fractool,Hurst_matrix_fractool))
#base_dfa = np.dstack((base_dfa,Hurst_matrix_dfa))
base_hurst = np.dstack((base_hurst,Hurst_matrix_hurst))
#base_corr_dim = np.dstack((base_corr_dim,Hurst_matrix_corr_dim))
#base_sampen = np.dstack((base_sampen,Hurst_matrix_sampen))
#save nifti file of base to the Outputs folder in FractalDimension
#name = name.rstrip("nii.gz")
#base_fractool = base_fractool[:,:,1:slice_num+1]
#ni_img = nib.nifti1.Nifti1Pair(base_fractool, None, n1_header)
#nib.save(ni_img, f'{name}_Hurst_Fractool.nii.gz')
#name = name.rstrip("nii.gz")
#base_dfa = base_dfa[:,:,1:slice_num+1]
#ni_img = nib.nifti1.Nifti1Pair(base_dfa, None, n1_header)
#nib.save(ni_img, f'{name}_Hurst_DFA.nii.gz')
name = name.rstrip("nii.gz")
base_hurst = base_hurst[:,:,1:slice_num+1]
ni_img = nib.nifti1.Nifti1Pair(base_hurst, None, n1_header)
nib.save(ni_img, f'{name}_Hurst_HurstRS.nii.gz')
#name = name.rstrip("nii.gz")
#base_corr_dim = base_corr_dim[:,:,1:slice_num+1]
#ni_img = nib.nifti1.Nifti1Pair(base_corr_dim, None, n1_header)
#nib.save(ni_img, f'{name}_Hurst_corrdim.nii.gz')
#name = name.rstrip("nii.gz")
#base_sampen = base_sampen[:,:,1:slice_num+1]
#ni_img = nib.nifti1.Nifti1Pair(base_sampen, None, n1_header)
#nib.save(ni_img, f'{name}_Hurst_sampen.nii.gz')