structural_mri_uncertainty/calculate_sigmap.py at main · InesGP/structural_mri_uncertainty · GitHub

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
import significantdigits
import os
import torchio as tio
import numpy as np
from significantdigits import Error, Method
import argparse
import brainload as bl

def calculate(folder, subject, ieee):

    t1s = []

    files = os.listdir(f"{folder}/{subject}")

    for f in files:
        if 'reg' in f:
            img = tio.ScalarImage(f"{folder}/{subject}/{f}")
            t1s.append(img.data.numpy().squeeze())

    mean = np.array(t1s)

    mean = np.sum(mean, axis=0)/len(mean)

    np.save(f"{folder}/{subject}/mean_{subject}.npy", mean)


    if ieee:
        t1s = []

        files = os.listdir(f"{folder}/{subject}/ieee")

        for f in files:
            img = tio.ScalarImage(f"{folder}/{subject}/ieee/{f}")
            t1s.append(img.data.numpy().squeeze())


    sig = significantdigits.significant_digits(
            array=t1s,
            reference=mean,
            base=2,
            axis=0,
            error=Error.Relative,
            method=Method.General,
        )

    if ieee: np.save(f"{folder}/{subject}/ieee/nonlinear_ieee.npy", sig)
    else: np.save(f"{folder}/{subject}/nonlinear_fuzzy.npy", sig)

def calculate_warp(folder, subject):
    warp_x = []
    warp_y = []
    warp_z = []

    files = os.listdir(f"{folder}/{subject}")

    for f in files:
        if 'm3z' in f:
            print(f"{folder}/{subject}/{f}")
            warp = bl.freesurferdata.read_m3z_file(f"{folder}/{subject}/{f}")
            warp_x.append(warp[0])
            warp_y.append(warp[1])
            warp_z.append(warp[2])

    mean_x = np.array(warp_x)
    mean_x = np.sum(mean_x, axis=0)/len(mean_x)

    mean_y = np.array(warp_y)
    mean_y = np.sum(mean_y, axis=0)/len(mean_y)

    mean_z = np.array(warp_z)
    mean_z = np.sum(mean_z, axis=0)/len(mean_z)


    sig_x = significantdigits.significant_digits(
            array=warp_x,
            reference=mean_x,
            base=2,
            axis=0,
            error=Error.Relative,
            method=Method.General,
        )

    sig_y = significantdigits.significant_digits(
            array=warp_y,
            reference=mean_y,
            base=2,
            axis=0,
            error=Error.Relative,
            method=Method.General,
        )

    sig_z = significantdigits.significant_digits(
            array=warp_z,
            reference=mean_z,
            base=2,
            axis=0,
            error=Error.Relative,
            method=Method.General,
        )

    np.save(f"{folder}/{subject}/nonlinear_fuzzy_x.npy", sig_x)
    np.save(f"{folder}/{subject}/nonlinear_fuzzy_y.npy", sig_y)
    np.save(f"{folder}/{subject}/nonlinear_fuzzy_z.npy", sig_z)

if __name__ == '__main__':
    parser = argparse.ArgumentParser(('Calculate significant digits for Freesurfer output'))

    parser.add_argument('-f', '--folder', required=True,
                        help=('Folder containing registered images'))
    parser.add_argument('-s', '--subject', required=True,
                        help='ID of the subject to use')
    parser.add_argument('-w', '--warp', action='store_true')
    parser.add_argument('--ieee', action='store_true',
                        help=('If set, IEEE significant digits are calculated'))

    args = parser.parse_args()

    if args.warp: calculate_warp(folder=args.folder, subject=args.subject)
    else: calculate(folder=args.folder, subject=args.subject, ieee=args.ieee)