|
1 | 1 | #!python |
2 | 2 | """PAR/REC to NIfTI converter |
3 | 3 | """ |
4 | | -from __future__ import division, print_function, absolute_import |
5 | 4 |
|
6 | | -from optparse import OptionParser, Option |
7 | | -import numpy as np |
8 | | -import numpy.linalg as npl |
9 | | -import sys |
10 | | -import os |
11 | | -import csv |
12 | | -import nibabel |
13 | | -import nibabel.parrec as pr |
14 | | -from nibabel.parrec import one_line |
15 | | -from nibabel.mriutils import calculate_dwell_time, MRIError |
16 | | -import nibabel.nifti1 as nifti1 |
17 | | -from nibabel.filename_parser import splitext_addext |
18 | | -from nibabel.volumeutils import fname_ext_ul_case |
19 | | -from nibabel.orientations import (io_orientation, inv_ornt_aff, |
20 | | - apply_orientation) |
21 | | -from nibabel.affines import apply_affine, from_matvec, to_matvec |
22 | | - |
23 | | - |
24 | | -def get_opt_parser(): |
25 | | - # use module docstring for help output |
26 | | - p = OptionParser( |
27 | | - usage="%s [OPTIONS] <PAR files>\n\n" % sys.argv[0] + __doc__, |
28 | | - version="%prog " + nibabel.__version__) |
29 | | - p.add_option( |
30 | | - Option("-v", "--verbose", action="store_true", dest="verbose", |
31 | | - default=False, |
32 | | - help="""Make some noise.""")) |
33 | | - p.add_option( |
34 | | - Option("-o", "--output-dir", action="store", type="string", |
35 | | - dest="outdir", default=None, |
36 | | - help=one_line("""Destination directory for NIfTI files. |
37 | | - Default: current directory."""))) |
38 | | - p.add_option( |
39 | | - Option("-c", "--compressed", action="store_true", |
40 | | - dest="compressed", default=False, |
41 | | - help="Whether to write compressed NIfTI files or not.")) |
42 | | - p.add_option( |
43 | | - Option("-p", "--permit-truncated", action="store_true", |
44 | | - dest="permit_truncated", default=False, |
45 | | - help=one_line( |
46 | | - """Permit conversion of truncated recordings. Support for |
47 | | - this is experimental, and results *must* be checked |
48 | | - afterward for validity."""))) |
49 | | - p.add_option( |
50 | | - Option("-b", "--bvs", action="store_true", dest="bvs", default=False, |
51 | | - help=one_line( |
52 | | - """Output bvals/bvecs files in addition to NIFTI |
53 | | - image."""))) |
54 | | - p.add_option( |
55 | | - Option("-d", "--dwell-time", action="store_true", default=False, |
56 | | - dest="dwell_time", |
57 | | - help=one_line( |
58 | | - """Calculate the scan dwell time. If supplied, the magnetic |
59 | | - field strength should also be supplied using |
60 | | - --field-strength (default 3). The field strength must be |
61 | | - supplied because it is not encoded in the PAR/REC |
62 | | - format."""))) |
63 | | - p.add_option( |
64 | | - Option("--field-strength", action="store", type="float", |
65 | | - dest="field_strength", |
66 | | - help=one_line( |
67 | | - """The magnetic field strength of the recording, only needed |
68 | | - for --dwell-time. The field strength must be supplied |
69 | | - because it is not encoded in the PAR/REC format."""))) |
70 | | - p.add_option( |
71 | | - Option("-i", "--volume-info", action="store_true", dest="vol_info", |
72 | | - default=False, |
73 | | - help=one_line( |
74 | | - """Export .PAR volume labels corresponding to the fourth |
75 | | - dimension of the data. The dimension info will be stored in |
76 | | - CSV format with the first row containing dimension labels |
77 | | - and the subsequent rows (one per volume), the corresponding |
78 | | - indices. Only labels that vary along the 4th dimension are |
79 | | - exported (e.g. for a single volume structural scan there |
80 | | - are no dynamic labels and no output file will be created). |
81 | | - """))) |
82 | | - p.add_option( |
83 | | - Option("--origin", action="store", dest="origin", default="scanner", |
84 | | - help=one_line( |
85 | | - """Reference point of the q-form transformation of the NIfTI |
86 | | - image. If 'scanner' the (0,0,0) coordinates will refer to |
87 | | - the scanner's iso center. If 'fov', this coordinate will be |
88 | | - the center of the recorded volume (field of view). Default: |
89 | | - 'scanner'."""))) |
90 | | - p.add_option( |
91 | | - Option("--minmax", action="store", nargs=2, dest="minmax", |
92 | | - help=one_line( |
93 | | - """Mininum and maximum settings to be stored in the NIfTI |
94 | | - header. If any of them is set to 'parse', the scaled data is |
95 | | - scanned for the actual minimum and maximum. To bypass this |
96 | | - potentially slow and memory intensive step (the data has to |
97 | | - be scaled and fully loaded into memory), fixed values can be |
98 | | - provided as space-separated pair, e.g. '5.4 120.4'. It is |
99 | | - possible to set a fixed minimum as scan for the actual |
100 | | - maximum (and vice versa). Default: 'parse parse'."""))) |
101 | | - p.set_defaults(minmax=('parse', 'parse')) |
102 | | - p.add_option( |
103 | | - Option("--store-header", action="store_true", dest="store_header", |
104 | | - default=False, |
105 | | - help=one_line( |
106 | | - """If set, all information from the PAR header is stored in |
107 | | - an extension ofthe NIfTI file header. Default: off"""))) |
108 | | - p.add_option( |
109 | | - Option("--scaling", action="store", dest="scaling", default='dv', |
110 | | - help=one_line( |
111 | | - """Choose data scaling setting. The PAR header defines two |
112 | | - different data scaling settings: 'dv' (values displayed on |
113 | | - console) and 'fp' (floating point values). Either one can be |
114 | | - chosen, or scaling can be disabled completely ('off'). Note |
115 | | - that neither method will actually scale the data, but just |
116 | | - store the corresponding settings in the NIfTI header, unless |
117 | | - non-uniform scaling is used, in which case the data is |
118 | | - stored in the file in scaled form. Default: 'dv'"""))) |
119 | | - p.add_option( |
120 | | - Option('--keep-trace', action="store_true", dest='keep_trace', |
121 | | - default=False, |
122 | | - help=one_line("""Do not discard the diagnostic Philips DTI |
123 | | - trace volume, if it exists in the data."""))) |
124 | | - p.add_option( |
125 | | - Option("--overwrite", action="store_true", dest="overwrite", |
126 | | - default=False, |
127 | | - help=one_line("""Overwrite file if it exists. Default: |
128 | | - False"""))) |
129 | | - p.add_option( |
130 | | - Option("--strict-sort", action="store_true", dest="strict_sort", |
131 | | - default=False, |
132 | | - help=one_line("""Use additional keys in determining the order |
133 | | - to sort the slices within the .REC file. This may be necessary |
134 | | - for more complicated scans with multiple echos, |
135 | | - cardiac phases, ASL label states, etc."""))) |
136 | | - return p |
137 | | - |
138 | | - |
139 | | -def verbose(msg, indent=0): |
140 | | - if verbose.switch: |
141 | | - print("%s%s" % (' ' * indent, msg)) |
142 | | - |
143 | | - |
144 | | -def error(msg, exit_code): |
145 | | - sys.stderr.write(msg + '\n') |
146 | | - sys.exit(exit_code) |
147 | | - |
148 | | - |
149 | | -def proc_file(infile, opts): |
150 | | - # figure out the output filename, and see if it exists |
151 | | - basefilename = splitext_addext(os.path.basename(infile))[0] |
152 | | - if opts.outdir is not None: |
153 | | - # set output path |
154 | | - basefilename = os.path.join(opts.outdir, basefilename) |
155 | | - |
156 | | - # prep a file |
157 | | - if opts.compressed: |
158 | | - verbose('Using gzip compression') |
159 | | - outfilename = basefilename + '.nii.gz' |
160 | | - else: |
161 | | - outfilename = basefilename + '.nii' |
162 | | - if os.path.isfile(outfilename) and not opts.overwrite: |
163 | | - raise IOError('Output file "%s" exists, use --overwrite to ' |
164 | | - 'overwrite it' % outfilename) |
165 | | - |
166 | | - # load the PAR header and data |
167 | | - scaling = 'dv' if opts.scaling == 'off' else opts.scaling |
168 | | - infile = fname_ext_ul_case(infile) |
169 | | - pr_img = pr.load(infile, |
170 | | - permit_truncated=opts.permit_truncated, |
171 | | - scaling=scaling, |
172 | | - strict_sort=opts.strict_sort) |
173 | | - pr_hdr = pr_img.header |
174 | | - affine = pr_hdr.get_affine(origin=opts.origin) |
175 | | - slope, intercept = pr_hdr.get_data_scaling(scaling) |
176 | | - if opts.scaling != 'off': |
177 | | - verbose('Using data scaling "%s"' % opts.scaling) |
178 | | - # get original scaling, and decide if we scale in-place or not |
179 | | - if opts.scaling == 'off': |
180 | | - slope = np.array([1.]) |
181 | | - intercept = np.array([0.]) |
182 | | - in_data = pr_img.dataobj.get_unscaled() |
183 | | - out_dtype = pr_hdr.get_data_dtype() |
184 | | - elif not np.any(np.diff(slope)) and not np.any(np.diff(intercept)): |
185 | | - # Single scalefactor case |
186 | | - slope = slope.ravel()[0] |
187 | | - intercept = intercept.ravel()[0] |
188 | | - in_data = pr_img.dataobj.get_unscaled() |
189 | | - out_dtype = pr_hdr.get_data_dtype() |
190 | | - else: |
191 | | - # Multi scalefactor case |
192 | | - slope = np.array([1.]) |
193 | | - intercept = np.array([0.]) |
194 | | - in_data = np.array(pr_img.dataobj) |
195 | | - out_dtype = np.float64 |
196 | | - # Reorient data block to LAS+ if necessary |
197 | | - ornt = io_orientation(np.diag([-1, 1, 1, 1]).dot(affine)) |
198 | | - if np.all(ornt == [[0, 1], |
199 | | - [1, 1], |
200 | | - [2, 1]]): # already in LAS+ |
201 | | - t_aff = np.eye(4) |
202 | | - else: # Not in LAS+ |
203 | | - t_aff = inv_ornt_aff(ornt, pr_img.shape) |
204 | | - affine = np.dot(affine, t_aff) |
205 | | - in_data = apply_orientation(in_data, ornt) |
206 | | - |
207 | | - bvals, bvecs = pr_hdr.get_bvals_bvecs() |
208 | | - if not opts.keep_trace: # discard Philips DTI trace if present |
209 | | - if bvecs is not None: |
210 | | - bad_mask = np.logical_and(bvals != 0, (bvecs == 0).all(axis=1)) |
211 | | - if bad_mask.sum() > 0: |
212 | | - pl = 's' if bad_mask.sum() != 1 else '' |
213 | | - verbose('Removing %s DTI trace volume%s' |
214 | | - % (bad_mask.sum(), pl)) |
215 | | - good_mask = ~bad_mask |
216 | | - in_data = in_data[..., good_mask] |
217 | | - bvals = bvals[good_mask] |
218 | | - bvecs = bvecs[good_mask] |
219 | | - |
220 | | - # Make corresponding NIfTI image |
221 | | - nimg = nifti1.Nifti1Image(in_data, affine, pr_hdr) |
222 | | - nhdr = nimg.header |
223 | | - nhdr.set_data_dtype(out_dtype) |
224 | | - nhdr.set_slope_inter(slope, intercept) |
225 | | - |
226 | | - if 'parse' in opts.minmax: |
227 | | - # need to get the scaled data |
228 | | - verbose('Loading (and scaling) the data to determine value range') |
229 | | - if opts.minmax[0] == 'parse': |
230 | | - nhdr['cal_min'] = in_data.min() * slope + intercept |
231 | | - else: |
232 | | - nhdr['cal_min'] = float(opts.minmax[0]) |
233 | | - if opts.minmax[1] == 'parse': |
234 | | - nhdr['cal_max'] = in_data.max() * slope + intercept |
235 | | - else: |
236 | | - nhdr['cal_max'] = float(opts.minmax[1]) |
237 | | - |
238 | | - # container for potential NIfTI1 header extensions |
239 | | - if opts.store_header: |
240 | | - # dump the full PAR header content into an extension |
241 | | - with open(infile, 'rb') as fobj: # contents must be bytes |
242 | | - hdr_dump = fobj.read() |
243 | | - dump_ext = nifti1.Nifti1Extension('comment', hdr_dump) |
244 | | - nhdr.extensions.append(dump_ext) |
245 | | - |
246 | | - verbose('Writing %s' % outfilename) |
247 | | - nibabel.save(nimg, outfilename) |
248 | | - |
249 | | - # write out bvals/bvecs if requested |
250 | | - if opts.bvs: |
251 | | - if bvals is None and bvecs is None: |
252 | | - verbose('No DTI volumes detected, bvals and bvecs not written') |
253 | | - elif bvecs is None: |
254 | | - verbose('DTI volumes detected, but no diffusion direction info was' |
255 | | - 'found. Writing .bvals file only.') |
256 | | - with open(basefilename + '.bvals', 'w') as fid: |
257 | | - # np.savetxt could do this, but it's just a loop anyway |
258 | | - for val in bvals: |
259 | | - fid.write('%s ' % val) |
260 | | - fid.write('\n') |
261 | | - else: |
262 | | - verbose('Writing .bvals and .bvecs files') |
263 | | - # Transform bvecs with reorientation affine |
264 | | - orig2new = npl.inv(t_aff) |
265 | | - bv_reorient = from_matvec(to_matvec(orig2new)[0], [0, 0, 0]) |
266 | | - bvecs = apply_affine(bv_reorient, bvecs) |
267 | | - with open(basefilename + '.bvals', 'w') as fid: |
268 | | - # np.savetxt could do this, but it's just a loop anyway |
269 | | - for val in bvals: |
270 | | - fid.write('%s ' % val) |
271 | | - fid.write('\n') |
272 | | - with open(basefilename + '.bvecs', 'w') as fid: |
273 | | - for row in bvecs.T: |
274 | | - for val in row: |
275 | | - fid.write('%s ' % val) |
276 | | - fid.write('\n') |
277 | | - |
278 | | - # export data labels varying along the 4th dimensions if requested |
279 | | - if opts.vol_info: |
280 | | - labels = pr_img.header.get_volume_labels() |
281 | | - if len(labels) > 0: |
282 | | - vol_keys = list(labels.keys()) |
283 | | - with open(basefilename + '.ordering.csv', 'w') as csvfile: |
284 | | - csvwriter = csv.writer(csvfile, delimiter=',') |
285 | | - csvwriter.writerow(vol_keys) |
286 | | - for vals in zip(*[labels[k] for k in vol_keys]): |
287 | | - csvwriter.writerow(vals) |
288 | | - |
289 | | - # write out dwell time if requested |
290 | | - if opts.dwell_time: |
291 | | - try: |
292 | | - dwell_time = calculate_dwell_time( |
293 | | - pr_hdr.get_water_fat_shift(), |
294 | | - pr_hdr.get_echo_train_length(), |
295 | | - opts.field_strength) |
296 | | - except MRIError: |
297 | | - verbose('No EPI factors, dwell time not written') |
298 | | - else: |
299 | | - verbose('Writing dwell time (%r sec) calculated assuming %sT ' |
300 | | - 'magnet' % (dwell_time, opts.field_strength)) |
301 | | - with open(basefilename + '.dwell_time', 'w') as fid: |
302 | | - fid.write('%r\n' % dwell_time) |
303 | | - # done |
304 | | - |
305 | | - |
306 | | -def main(): |
307 | | - parser = get_opt_parser() |
308 | | - (opts, infiles) = parser.parse_args() |
309 | | - |
310 | | - verbose.switch = opts.verbose |
311 | | - |
312 | | - if opts.origin not in ['scanner', 'fov']: |
313 | | - error("Unrecognized value for --origin: '%s'." % opts.origin, 1) |
314 | | - if opts.dwell_time and opts.field_strength is None: |
315 | | - error('Need --field-strength for dwell time calculation', 1) |
316 | | - |
317 | | - # store any exceptions |
318 | | - errs = [] |
319 | | - for infile in infiles: |
320 | | - verbose('Processing %s' % infile) |
321 | | - try: |
322 | | - proc_file(infile, opts) |
323 | | - except Exception as e: |
324 | | - errs.append('%s: %s' % (infile, e)) |
325 | | - |
326 | | - if len(errs): |
327 | | - error('Caught %i exceptions. Dump follows:\n\n %s' |
328 | | - % (len(errs), '\n'.join(errs)), 1) |
329 | | - else: |
330 | | - verbose('Done') |
| 5 | +from nibabel.parrec2nii_cmd import main |
331 | 6 |
|
332 | 7 |
|
333 | 8 | if __name__ == '__main__': |
|
0 commit comments