/
xfmavg_scipy.py
executable file
·125 lines (102 loc) · 4.02 KB
/
xfmavg_scipy.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
#! /usr/bin/env python
import argparse
import subprocess
import traceback
import os
import numpy as np
# needed for matrix log and exp
import scipy.linalg
# needed to read and write XFM files
import minc
def do_cmd(cmds,verbose=False):
try:
if not verbose:
with open(os.devnull, "w") as fnull:
p=subprocess.Popen(cmds, stdout=fnull, stderr=subprocess.PIPE)
else:
p=subprocess.Popen(cmds, stderr=subprocess.PIPE)
(output,output_stderr)=p.communicate()
outvalue=p.wait()
except OSError:
raise Exception("ERROR: command {} Error:{}!\nMessage: {}\n{}".format(str(cmds),str(outvalue),output_stderr,traceback.format_exc()))
if not outvalue == 0:
raise Exception("ERROR: command {} failed {}!\nMessage: {}\n{}".format(str(cmds),str(outvalue),output_stderr,traceback.format_exc()))
return outvalue
def xfmavg(inputs, output, verbose=False):
# TODO: handl inversion flag correctly
all_linear=True
all_nonlinear=True
input_xfms=[]
for j in inputs:
x=minc.read_xfm(j)
if x[0].lin and len(x)==1 and (not x[0].inv):
# this is a linear matrix
input_xfms.append(x[0])
else:
all_linear &= False
# strip identity matrixes
nl = []
_identity = np.asmatrix(np.identity(4))
_eps = 1e-6
for i in x:
if i.lin:
if scipy.linalg.norm(_identity-i.trans)>_eps: # this is non-identity matrix
all_nonlinear&=False
else:
if i.inv:
raise Exception("Trying to average inverted nonlinear transform")
nl.append(i)
if len(nl)!=1:
all_nonlinear &= False
else:
input_xfms.append(nl[0])
if all_linear:
acc = np.asmatrix(np.zeros([4,4],dtype=np.complex))
for i in input_xfms:
acc+=scipy.linalg.logm(i.trans)
acc/=len(input_xfms)
out_xfm=[minc.xfm_entry(True,False,scipy.linalg.expm(acc).real)]
minc.write_xfm(output, out_xfm)
elif all_nonlinear:
input_grids=[]
for i in input_xfms:
input_grids.append(i.trans)
output_grid=output.rsplit('.xfm',1)[0]+'_grid_0.mnc'
cmds=['mincaverage','-clob']
cmds.extend(input_grids)
cmds.append(output_grid)
do_cmd(cmds,verbose=verbose)
out_xfm=[minc.xfm_entry(False,False,output_grid)]
minc.write_xfm(output, out_xfm)
else:
raise Exception("Mixed XFM files provided as input")
def parse_options():
parser = argparse.ArgumentParser(formatter_class=argparse.ArgumentDefaultsHelpFormatter,
description='Perform ')
parser.add_argument("--verbose",
action="store_true",
dest="verbose",
default=False,
help="Print verbose information" )
parser.add_argument("--clobber",
action="store_true",
dest="clobber",
default=False,
help="Overwrite output" )
parser.add_argument("inputs",
nargs='*',
help="Input xfm files")
parser.add_argument("output",
help="Output xfm file")
options = parser.parse_args()
return options
if __name__ == '__main__':
options = parse_options()
if options.inputs is not None and options.output is not None:
if not options.clobber and os.path.exists(options.output):
raise Exception("File {} exists! Run with --cloberr to overwrite".format(options.output))
xfmavg(options.inputs,options.output,verbose=options.verbose)
else:
print("Refusing to run without input data, run --help")
exit(1)
# kate: space-indent on; indent-width 4; indent-mode python;replace-tabs on;word-wrap-column 80