1
""" Read and write trackvis files (old interface)
2

3
See :mod:`nibabel.streamlines` for the new interface.
4

5
We will deprecate this, the old interface, in some future release.
6
"""
7 33
import warnings
8 33
import struct
9 33
import itertools
10

11 33
import numpy as np
12 33
import numpy.linalg as npl
13 33
from numpy.compat.py3k import asstr
14

15 33
from .volumeutils import (native_code, swapped_code, endian_codes, rec2dict)
16 33
from .openers import ImageOpener
17 33
from .orientations import aff2axcodes
18 33
from .affines import apply_affine
19 33
from .deprecated import deprecate_with_version
20

21 33
warnings.warn("The trackvis interface has been deprecated and will be removed "
22
              "in v4.0; please use the 'nibabel.streamlines' interface.",
23
              DeprecationWarning,
24
              stacklevel=2)
25

26
# Definition of trackvis header structure.
27
# See http://www.trackvis.org/docs/?subsect=fileformat
28
# See https://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html
29 33
header_1_dtd = [
30
    ('id_string', 'S6'),
31
    ('dim', 'h', 3),
32
    ('voxel_size', 'f4', 3),
33
    ('origin', 'f4', 3),
34
    ('n_scalars', 'h'),
35
    ('scalar_name', 'S20', 10),
36
    ('n_properties', 'h'),
37
    ('property_name', 'S20', 10),
38
    ('reserved', 'S508'),
39
    ('voxel_order', 'S4'),
40
    ('pad2', 'S4'),
41
    ('image_orientation_patient', 'f4', 6),
42
    ('pad1', 'S2'),
43
    ('invert_x', 'S1'),
44
    ('invert_y', 'S1'),
45
    ('invert_z', 'S1'),
46
    ('swap_xy', 'S1'),
47
    ('swap_yz', 'S1'),
48
    ('swap_zx', 'S1'),
49
    ('n_count', 'i4'),
50
    ('version', 'i4'),
51
    ('hdr_size', 'i4'),
52
]
53

54
# Version 2 adds a 4x4 matrix giving the affine transformation going
55
# from voxel coordinates in the referenced 3D voxel matrix, to xyz
56
# coordinates (axes L->R, P->A, I->S).  IF (0 based) value [3, 3] from
57
# this matrix is 0, this means the matrix is not recorded.
58 33
header_2_dtd = [
59
    ('id_string', 'S6'),
60
    ('dim', 'h', 3),
61
    ('voxel_size', 'f4', 3),
62
    ('origin', 'f4', 3),
63
    ('n_scalars', 'h'),
64
    ('scalar_name', 'S20', 10),
65
    ('n_properties', 'h'),
66
    ('property_name', 'S20', 10),
67
    ('vox_to_ras', 'f4', (4, 4)),  # new field for version 2
68
    ('reserved', 'S444'),
69
    ('voxel_order', 'S4'),
70
    ('pad2', 'S4'),
71
    ('image_orientation_patient', 'f4', 6),
72
    ('pad1', 'S2'),
73
    ('invert_x', 'S1'),
74
    ('invert_y', 'S1'),
75
    ('invert_z', 'S1'),
76
    ('swap_xy', 'S1'),
77
    ('swap_yz', 'S1'),
78
    ('swap_zx', 'S1'),
79
    ('n_count', 'i4'),
80
    ('version', 'i4'),
81
    ('hdr_size', 'i4'),
82
]
83

84
# Full header numpy dtypes
85 33
header_1_dtype = np.dtype(header_1_dtd)
86 33
header_2_dtype = np.dtype(header_2_dtd)
87

88
# affine to go from DICOM LPS to MNI RAS space
89 33
DPCS_TO_TAL = np.diag([-1, -1, 1, 1])
90

91

92 33
class HeaderError(Exception):
93
    """ Error in trackvis header
94
    """
95

96

97 33
class DataError(Exception):
98
    """ Error in trackvis data
99
    """
100

101

102 33
@deprecate_with_version('trackvis.read is deprecated; please use '
103
                        'nibabel.streamlines.load, instead.',
104
                        since='2.5.0', until='4.0.0')
105 33
def read(fileobj, as_generator=False, points_space=None, strict=True):
106
    """ Read trackvis file from `fileobj`, return `streamlines`, `header`
107

108
    Parameters
109
    ----------
110
    fileobj : string or file-like object
111
       If string, a filename; otherwise an open file-like object
112
       pointing to trackvis file (and ready to read from the beginning
113
       of the trackvis header data)
114
    as_generator : bool, optional
115
       Whether to return tracks as sequence (False, default) or as a generator
116
       (True).
117
    points_space : {None, 'voxel', 'rasmm'}, optional
118
        The coordinates in which you want the points in the *output*
119
        streamlines expressed.  If None, then return the points exactly as they
120
        are stored in the trackvis file. The points will probably be in
121
        trackvis voxmm space - see Notes for ``write`` function.  If 'voxel',
122
        we convert the points to voxel space simply by dividing by the recorded
123
        voxel size.  If 'rasmm' we'll convert the points to RAS mm space (real
124
        space). For 'rasmm' we check if the affine is set and matches the voxel
125
        sizes and voxel order.
126
    strict : {True, False}, optional
127
        If True, raise error on read for badly-formed file.  If False, let pass
128
        files with last track having too few points.
129

130
    Returns
131
    -------
132
    streamlines : sequence or generator
133
       Returns sequence if `as_generator` is False, generator if True.  Value
134
       is sequence or generator of 3 element sequences with elements:
135

136
       #. points : ndarray shape (N,3)
137
          where N is the number of points
138
       #. scalars : None or ndarray shape (N, M)
139
          where M is the number of scalars per point
140
       #. properties : None or ndarray shape (P,)
141
          where P is the number of properties
142

143
    hdr : structured array
144
       structured array with trackvis header fields
145

146
    Notes
147
    -----
148
    The endianness of the input data can be deduced from the endianness
149
    of the returned `hdr` or `streamlines`
150

151
    Points are in trackvis *voxel mm*.  Each track has N points, each with 3
152
    coordinates, ``x, y, z``, where ``x`` is the floating point voxel
153
    coordinate along the first image axis, multiplied by the voxel size for
154
    that axis.
155
    """
156 33
    fileobj = ImageOpener(fileobj)
157 33
    hdr_str = fileobj.read(header_2_dtype.itemsize)
158
    # try defaulting to version 2 format
159 33
    hdr = np.ndarray(shape=(),
160
                     dtype=header_2_dtype,
161
                     buffer=hdr_str)
162 33
    if hdr['id_string'].item()[:5] != b'TRACK':
163 0
        raise HeaderError('Expecting TRACK as first '
164
                          '5 characters of id_string')
165 33
    if hdr['hdr_size'] == 1000:
166 33
        endianness = native_code
167
    else:
168 33
        hdr = hdr.newbyteorder()
169 33
        if hdr['hdr_size'] != 1000:
170 0
            raise HeaderError(f"Invalid hdr_size of {hdr['hdr_size']}")
171 33
        endianness = swapped_code
172
    # Check version and adapt structure accordingly
173 33
    version = hdr['version']
174 33
    if version not in (1, 2):
175 0
        raise HeaderError('Reader only supports versions 1 and 2')
176 33
    if version == 1:  # make a new header with the same data
177 33
        hdr = np.ndarray(shape=(),
178
                         dtype=header_1_dtype,
179
                         buffer=hdr_str)
180 33
        if endianness == swapped_code:
181 33
            hdr = hdr.newbyteorder()
182
    # Do points_space checks
183 33
    _check_hdr_points_space(hdr, points_space)
184
    # prepare transforms for later use
185 33
    if points_space == 'voxel':
186 33
        zooms = hdr['voxel_size'][None, :].astype('f4')
187 33
    elif points_space == 'rasmm':
188 33
        zooms = hdr['voxel_size']
189 33
        affine = hdr['vox_to_ras']
190 33
        tv2vx = np.diag((1. / zooms).tolist() + [1])
191 33
        tv2mm = np.dot(affine, tv2vx).astype('f4')
192 33
    n_s = hdr['n_scalars']
193 33
    n_p = hdr['n_properties']
194 33
    f4dt = np.dtype(endianness + 'f4')
195 33
    pt_cols = 3 + n_s
196 33
    pt_size = int(f4dt.itemsize * pt_cols)
197 33
    ps_size = int(f4dt.itemsize * n_p)
198 33
    i_fmt = endianness + 'i'
199 33
    stream_count = hdr['n_count']
200 33
    if stream_count < 0:
201 0
        raise HeaderError('Unexpected negative n_count')
202

203 33
    def track_gen():
204
        # For case where there are no scalars or no properties
205 33
        scalars = None
206 33
        ps = None
207 33
        n_streams = 0
208
        # stream_count == 0 signals read to end of file
209 33
        n_streams_required = stream_count if stream_count != 0 else np.inf
210 33
        end_of_file = False
211 33
        while not end_of_file and n_streams < n_streams_required:
212 33
            n_str = fileobj.read(4)
213 33
            if len(n_str) < 4:
214 33
                break
215 33
            n_pts = struct.unpack(i_fmt, n_str)[0]
216
            # Check if we got as many bytes as we expect for these points
217 33
            exp_len = n_pts * pt_size
218 33
            pts_str = fileobj.read(exp_len)
219 33
            if len(pts_str) != exp_len:
220
                # Short of bytes, should we raise an error or continue?
221 33
                actual_n_pts = int(len(pts_str) / pt_size)
222 33
                if actual_n_pts != n_pts:
223 33
                    if strict:
224 33
                        raise DataError(f'Expecting {n_pts} points for stream '
225
                                        f'{n_streams}, found {actual_n_pts}')
226 33
                    n_pts = actual_n_pts
227 33
                    end_of_file = True
228
            # Cast bytes to points array
229 33
            pts = np.ndarray(shape=(n_pts, pt_cols), dtype=f4dt,
230
                             buffer=pts_str)
231
            # Add properties
232 33
            if n_p:
233 33
                ps_str = fileobj.read(ps_size)
234 33
                ps = np.ndarray(shape=(n_p,), dtype=f4dt, buffer=ps_str)
235 33
            xyz = pts[:, :3]
236 33
            if points_space == 'voxel':
237 33
                xyz = xyz / zooms
238 33
            elif points_space == 'rasmm':
239 33
                xyz = apply_affine(tv2mm, xyz)
240 33
            if n_s:
241 33
                scalars = pts[:, 3:]
242 33
            yield (xyz, scalars, ps)
243 33
            n_streams += 1
244
        # Always close file if we opened it
245 33
        fileobj.close_if_mine()
246
        # Raise error if we didn't get as many streams as claimed
247 33
        if n_streams_required != np.inf and n_streams < n_streams_required:
248 33
            raise DataError(
249
                f'Expecting {stream_count} streamlines, found only {n_streams}')
250

251 33
    streamlines = track_gen()
252 33
    if not as_generator:
253 33
        streamlines = list(streamlines)
254 33
    return streamlines, hdr
255

256

257 33
@deprecate_with_version('trackvis.write is deprecated; please use '
258
                        'nibabel.streamlines.save, instead.',
259
                        since='2.5.0', until='4.0.0')
260 33
def write(fileobj, streamlines, hdr_mapping=None, endianness=None,
261
          points_space=None):
262
    """ Write header and `streamlines` to trackvis file `fileobj`
263

264
    The parameters from the streamlines override conflicting parameters
265
    in the `hdr_mapping` information.  In particular, the number of
266
    streamlines, the number of scalars, and the number of properties are
267
    written according to `streamlines` rather than `hdr_mapping`.
268

269
    Parameters
270
    ----------
271
    fileobj : filename or file-like
272
       If filename, open file as 'wb', otherwise `fileobj` should be an
273
       open file-like object, with a ``write`` method.
274
    streamlines : iterable
275
       iterable returning 3 element sequences with elements:
276

277
       #. points : ndarray shape (N,3)
278
          where N is the number of points
279
       #. scalars : None or ndarray shape (N, M)
280
          where M is the number of scalars per point
281
       #. properties : None or ndarray shape (P,)
282
          where P is the number of properties
283

284
       If `streamlines` has a ``len`` (for example, it is a list or a tuple),
285
       then we can write the number of streamlines into the header.  Otherwise
286
       we write 0 for the number of streamlines (a valid trackvis header) and
287
       write streamlines into the file until the iterable is exhausted.
288
       M - the number of scalars - has to be the same for each streamline in
289
       `streamlines`.  Similarly for P. See `points_space` and Notes for more
290
       detail on the coordinate system for ``points`` above.
291
    hdr_mapping : None, ndarray or mapping, optional
292
       Information for filling header fields.  Can be something
293
       dict-like (implementing ``items``) or a structured numpy array
294
    endianness : {None, '<', '>'}, optional
295
       Endianness of file to be written.  '<' is little-endian, '>' is
296
       big-endian.  None (the default) is to use the endianness of the
297
       `streamlines` data.
298
    points_space : {None, 'voxel', 'rasmm'}, optional
299
        The coordinates in which the points in the input streamlines are
300
        expressed.  If None, then assume the points are as you want them
301
        (probably trackvis voxmm space - see Notes).  If 'voxel', the points
302
        are in voxel space, and we will transform them to trackvis voxmm space.
303
        If 'rasmm' the points are in RAS mm space (real space).  We transform
304
        them to trackvis voxmm space.  If 'voxel' or 'rasmm' we insist that the
305
        voxel sizes and ordering are set to non-default values.  If 'rasmm' we
306
        also check if the affine is set and matches the voxel sizes
307

308
    Returns
309
    -------
310
    None
311

312
    Examples
313
    --------
314
    >>> from io import BytesIO
315
    >>> file_obj = BytesIO()
316
    >>> pts0 = np.random.uniform(size=(10,3))
317
    >>> pts1 = np.random.uniform(size=(10,3))
318
    >>> streamlines = ([(pts0, None, None), (pts1, None, None)])
319
    >>> write(file_obj, streamlines)
320
    >>> _ = file_obj.seek(0)  # returns 0
321
    >>> streams, hdr = read(file_obj)
322
    >>> len(streams)
323
    2
324

325
    If there are too many streamlines to fit in memory, you can pass an
326
    iterable thing instead of a list
327

328
    >>> file_obj = BytesIO()
329
    >>> def gen():
330
    ...     yield (pts0, None, None)
331
    ...     yield (pts0, None, None)
332
    >>> write(file_obj, gen())
333
    >>> _ = file_obj.seek(0)
334
    >>> streams, hdr = read(file_obj)
335
    >>> len(streams)
336
    2
337

338
    Notes
339
    -----
340
    Trackvis (the application) expects the ``points`` in the streamlines be in
341
    what we call *trackvis voxmm* coordinates.  If we have a point (x, y, z) in
342
    voxmm coordinates, and ``voxel_size`` has the voxel sizes for each of the 3
343
    dimensions, then x, y, z refer to mm in voxel space. Thus if i, j, k is a
344
    point in voxel coordinates, then ``x = i * voxel_size[0]; y = j *
345
    voxel_size[1]; z = k * voxel_size[2]``.   The spatial direction of x, y and
346
    z are defined with the "voxel_order" field.  For example, if the original
347
    image had RAS voxel ordering then "voxel_order" would be "RAS".  RAS here
348
    refers to the spatial direction of the voxel axes: "R" means that moving
349
    along first voxel axis moves from left to right in space, "A" -> second
350
    axis goes from posterior to anterior, "S" -> inferior to superior.  If
351
    "voxel_order" is empty we assume "LPS".
352

353
    This information comes from some helpful replies on the trackvis forum
354
    about `interpreting point coordiantes
355
    <http://trackvis.org/blog/forum/diffusion-toolkit-usage/interpretation-of-track-point-coordinates>`_
356
    """
357 33
    stream_iter = iter(streamlines)
358 33
    try:
359 33
        streams0 = next(stream_iter)
360 33
    except StopIteration:  # empty sequence or iterable
361
        # write header without streams
362 33
        hdr = _hdr_from_mapping(None, hdr_mapping, endianness)
363 33
        with ImageOpener(fileobj, 'wb') as fileobj:
364 33
            fileobj.write(hdr.tobytes())
365 33
        return
366 33
    if endianness is None:
367 33
        endianness = endian_codes[streams0[0].dtype.byteorder]
368
    # fill in a new header from mapping-like
369 33
    hdr = _hdr_from_mapping(None, hdr_mapping, endianness)
370
    # Try and get number of streams from streamlines.  If this is an iterable,
371
    # we don't have a len, so we write 0 for length.  The 0 is a valid trackvis
372
    # value with meaning - keep reading until you run out of data.
373 33
    try:
374 33
        n_streams = len(streamlines)
375 33
    except TypeError:  # iterable; we don't know the number of streams
376 33
        n_streams = 0
377 33
    hdr['n_count'] = n_streams
378
    # Get number of scalars and properties
379 33
    pts, scalars, props = streams0
380
    # calculate number of scalars
381 33
    if scalars is not None:
382 33
        n_s = scalars.shape[1]
383
    else:
384 33
        n_s = 0
385 33
    hdr['n_scalars'] = n_s
386
    # calculate number of properties
387 33
    if props is not None:
388 33
        n_p = props.size
389 33
        hdr['n_properties'] = n_p
390
    else:
391 33
        n_p = 0
392
    # do points_space checks
393 33
    _check_hdr_points_space(hdr, points_space)
394
    # prepare transforms for later use
395 33
    if points_space == 'voxel':
396 33
        zooms = hdr['voxel_size'][None, :].astype('f4')
397 33
    elif points_space == 'rasmm':
398 33
        zooms = hdr['voxel_size']
399 33
        affine = hdr['vox_to_ras']
400 33
        vx2tv = np.diag(zooms.tolist() + [1])
401 33
        mm2vx = npl.inv(affine)
402 33
        mm2tv = np.dot(vx2tv, mm2vx).astype('f4')
403
    # write header
404 33
    fileobj = ImageOpener(fileobj, mode='wb')
405 33
    fileobj.write(hdr.tobytes())
406
    # track preliminaries
407 33
    f4dt = np.dtype(endianness + 'f4')
408 33
    i_fmt = endianness + 'i'
409
    # Add back the read first streamline to the sequence
410 33
    for pts, scalars, props in itertools.chain([streams0], stream_iter):
411 33
        n_pts, n_coords = pts.shape
412 33
        if n_coords != 3:
413 0
            raise ValueError('pts should have 3 columns')
414 33
        fileobj.write(struct.pack(i_fmt, n_pts))
415 33
        if points_space == 'voxel':
416 33
            pts = pts * zooms
417 33
        elif points_space == 'rasmm':
418 33
            pts = apply_affine(mm2tv, pts)
419
        # This call ensures that the data are 32-bit floats, and that
420
        # the endianness is OK.
421 33
        if pts.dtype != f4dt:
422 33
            pts = pts.astype(f4dt)
423 33
        if n_s == 0:
424 33
            if not (scalars is None or len(scalars) == 0):
425 33
                raise DataError('Expecting 0 scalars per point')
426
        else:
427 33
            if scalars.shape != (n_pts, n_s):
428 33
                raise DataError(f'Scalars should be shape ({n_pts}, {n_s})')
429 33
            if scalars.dtype != f4dt:
430 33
                scalars = scalars.astype(f4dt)
431 33
            pts = np.c_[pts, scalars]
432 33
        fileobj.write(pts.tobytes())
433 33
        if n_p == 0:
434 33
            if not (props is None or len(props) == 0):
435 33
                raise DataError('Expecting 0 properties per point')
436
        else:
437 33
            if props.size != n_p:
438 33
                raise DataError(f'Properties should be size {n_p}')
439 33
            if props.dtype != f4dt:
440 33
                props = props.astype(f4dt)
441 33
            fileobj.write(props.tobytes())
442 33
    fileobj.close_if_mine()
443

444

445 33
def _check_hdr_points_space(hdr, points_space):
446
    """ Check header `hdr` for consistency with transform `points_space`
447

448
    Parameters
449
    ----------
450
    hdr : ndarray
451
        trackvis header as structured ndarray
452
    points_space : {None, 'voxmm', 'voxel', 'rasmm'
453
        nature of transform that we will (elsewhere) apply to streamlines
454
        paired with `hdr`.  None or 'voxmm' means pass through with no futher
455
        checks.  'voxel' checks for all ``hdr['voxel_sizes'] being <= zero
456
        (error) or any being zero (warning).  'rasmm' checks for presence of
457
        non-zeros affine in ``hdr['vox_to_ras']``, and that the affine therein
458
        corresponds to ``hdr['voxel_order']`` and ''hdr['voxel_sizes']`` - and
459
        raises an error otherwise.
460

461
    Returns
462
    -------
463
    None
464

465
    Notes
466
    -----
467
    """
468 33
    if points_space is None or points_space == 'voxmm':
469 33
        return
470 33
    if points_space == 'voxel':
471 33
        voxel_size = hdr['voxel_size']
472 33
        if np.any(voxel_size < 0):
473 33
            raise HeaderError(f'Negative voxel sizes {voxel_size} not '
474
                              'valid for voxel - voxmm conversion')
475 33
        if np.all(voxel_size == 0):
476 33
            raise HeaderError('Cannot convert between voxels and voxmm when '
477
                              '"voxel_sizes" all 0')
478 33
        if np.any(voxel_size == 0):
479 33
            warnings.warn(f'zero values in "voxel_size" - {voxel_size}')
480 33
        return
481 33
    elif points_space == 'rasmm':
482 33
        try:
483 33
            affine = hdr['vox_to_ras']
484 33
        except ValueError:
485 33
            raise HeaderError('Need "vox_to_ras" field to get '
486
                              'affine with which to convert points; '
487
                              'this is present for headers >= version 2')
488 33
        if np.all(affine == 0) or affine[3, 3] == 0:
489 33
            raise HeaderError('Need non-zero affine to convert between '
490
                              'rasmm points and voxmm')
491 33
        zooms = hdr['voxel_size']
492 33
        aff_zooms = np.sqrt(np.sum(affine[:3, :3]**2, axis=0))
493 33
        if not np.allclose(aff_zooms, zooms):
494 33
            raise HeaderError(f'Affine zooms {aff_zooms} differ '
495
                              f'from voxel_size field value {zooms}')
496 33
        aff_order = ''.join(aff2axcodes(affine))
497 33
        voxel_order = asstr(hdr['voxel_order'].item())
498 33
        if voxel_order == '':
499 33
            voxel_order = 'LPS'  # trackvis default
500 33
        if not voxel_order == aff_order:
501 33
            raise HeaderError(f'Affine implies voxel_order {aff_order} '
502
                              f'but header voxel_order is {voxel_order}')
503
    else:
504 33
        raise ValueError(f'Painfully confusing "points_space" value of "{points_space}"')
505

506

507 33
def _hdr_from_mapping(hdr=None, mapping=None, endianness=native_code):
508
    """ Fill `hdr` from mapping `mapping`, with given endianness """
509 33
    if hdr is None:
510
        # passed a valid mapping as header?  Copy and return
511 33
        if isinstance(mapping, np.ndarray):
512 33
            test_dtype = mapping.dtype.newbyteorder('=')
513 33
            if test_dtype in (header_1_dtype, header_2_dtype):
514 33
                return mapping.copy()
515
        # otherwise make a new empty header.   If no version specified,
516
        # go for default (2)
517 33
        if mapping is None:
518 33
            version = 2
519
        else:
520 33
            version = mapping.get('version', 2)
521 33
        hdr = empty_header(endianness, version)
522 33
    if mapping is None:
523 33
        return hdr
524 33
    if isinstance(mapping, np.ndarray):
525 0
        mapping = rec2dict(mapping)
526 33
    for key, value in mapping.items():
527 33
        hdr[key] = value
528
    # check header values
529 33
    if hdr['id_string'].item()[:5] != b'TRACK':
530 33
        raise HeaderError('Expecting TRACK as first '
531
                          '5 characaters of id_string')
532 33
    if hdr['version'] not in (1, 2):
533 0
        raise HeaderError('Reader only supports version 1')
534 33
    if hdr['hdr_size'] != 1000:
535 33
        raise HeaderError('hdr_size should be 1000')
536 33
    return hdr
537

538

539 33
@deprecate_with_version('empty_header is deprecated; please use '
540
                        'nibabel.streamlines.TrkFile.create_empty_header, instead.',
541
                        since='2.5.0', until='4.0.0')
542 33
def empty_header(endianness=None, version=2):
543
    """ Empty trackvis header
544

545
    Parameters
546
    ----------
547
    endianness : {'<','>'}, optional
548
       Endianness of empty header to return. Default is native endian.
549
    version : int, optional
550
       Header version.  1 or 2.  Default is 2
551

552
    Returns
553
    -------
554
    hdr : structured array
555
       structured array containing empty trackvis header
556

557
    Examples
558
    --------
559
    >>> hdr = empty_header()
560
    >>> print(hdr['version'])
561
    2
562
    >>> hdr['id_string'].item() == b'TRACK'
563
    True
564
    >>> endian_codes[hdr['version'].dtype.byteorder] == native_code
565
    True
566
    >>> hdr = empty_header(swapped_code)
567
    >>> endian_codes[hdr['version'].dtype.byteorder] == swapped_code
568
    True
569
    >>> hdr = empty_header(version=1)
570
    >>> print(hdr['version'])
571
    1
572

573
    Notes
574
    -----
575
    The trackvis header can store enough information to give an affine
576
    mapping between voxel and world space.  Often this information is
577
    missing.  We make no attempt to fill it with sensible defaults on
578
    the basis that, if the information is missing, it is better to be
579
    explicit.
580
    """
581 33
    if version == 1:
582 33
        dt = header_1_dtype
583 33
    elif version == 2:
584 33
        dt = header_2_dtype
585
    else:
586 33
        raise HeaderError('Header version should be 1 or 2')
587 33
    if endianness:
588 33
        dt = dt.newbyteorder(endianness)
589 33
    hdr = np.zeros((), dtype=dt)
590 33
    hdr['id_string'] = 'TRACK'
591 33
    hdr['version'] = version
592 33
    hdr['hdr_size'] = 1000
593 33
    return hdr
594

595

596 33
@deprecate_with_version('aff_from_hdr is deprecated; please use '
597
                        'nibabel.streamlines.trk.get_affine_trackvis_to_rasmm, instead.',
598
                        since='2.5.0', until='4.0.0')
599 33
def aff_from_hdr(trk_hdr, atleast_v2=True):
600
    """ Return voxel to mm affine from trackvis header
601

602
    Affine is mapping from voxel space to Nifti (RAS) output coordinate
603
    system convention; x: Left -> Right, y: Posterior -> Anterior, z:
604
    Inferior -> Superior.
605

606
    Parameters
607
    ----------
608
    trk_hdr : mapping
609
       Mapping with trackvis header keys ``version``. If ``version == 2``, we
610
       also expect ``vox_to_ras``.
611
    atleast_v2 : None or bool
612
        If None, currently defaults to False.  This will change to True in
613
        future versions.  If True, require that there is a valid 'vox_to_ras'
614
        affine, raise HeaderError otherwise.  If False, look for valid
615
        'vox_to_ras' affine, but fall back to best guess from version 1 fields
616
        otherwise.
617

618
    Returns
619
    -------
620
    aff : (4,4) array
621
       affine giving mapping from voxel coordinates (affine applied on
622
       the left to points on the right) to millimeter coordinates in the
623
       RAS coordinate system
624

625
    Notes
626
    -----
627
    Our initial idea was to try and work round the deficiencies of the version
628
    1 format by using the DICOM orientation fields to store the affine.  This
629
    proved difficult in practice because trackvis (the application) doesn't
630
    allow negative voxel sizes (needed for recording axis flips) and sets the
631
    origin field to 0. In future, we'll raise an error rather than try and
632
    estimate the affine from version 1 fields
633
    """
634 33
    if trk_hdr['version'] == 2:
635 33
        aff = trk_hdr['vox_to_ras']
636 33
        if aff[3, 3] != 0:
637 33
            return aff
638 33
        if atleast_v2:
639 33
            raise HeaderError('Requiring version 2 affine and this affine is '
640
                              'not valid')
641
    # Now we are in the dark world of the DICOM fields.  We might have made
642
    # this one ourselves, in which case the origin might be set, and it might
643
    # have negative voxel sizes
644 33
    aff = np.eye(4)
645
    # The IOP field has only two of the three columns we need
646 33
    iop = trk_hdr['image_orientation_patient'].reshape(2, 3).T
647
    # R might be a rotation matrix (and so completed by the cross product of
648
    # the first two columns), or it might be an orthogonal matrix with negative
649
    # determinant. We try pure rotation first
650 33
    R = np.c_[iop, np.cross(*iop.T)]
651 33
    vox = trk_hdr['voxel_size']
652 33
    aff[:3, :3] = R * vox
653 33
    aff[:3, 3] = trk_hdr['origin']
654 33
    aff = np.dot(DPCS_TO_TAL, aff)
655
    # Next we check against the 'voxel_order' field if present and not empty.
656 33
    try:
657 33
        voxel_order = asstr(trk_hdr['voxel_order'].item())
658 33
    except (KeyError, ValueError):
659 33
        voxel_order = ''
660 33
    if voxel_order == '':
661 33
        return aff
662
    # If the voxel_order conflicts with the affine by one flip, this may have
663
    # been a negative determinant affine saved with positive voxel sizes
664 33
    exp_order = ''.join(aff2axcodes(aff))
665 33
    if voxel_order != exp_order:
666
        # If first pass doesn't match, try flipping the (estimated) third
667
        # column
668 33
        aff[:, 2] *= -1
669 33
        exp_order = ''.join(aff2axcodes(aff))
670 33
        if voxel_order != exp_order:
671 33
            raise HeaderError('Estimate of header affine does not match '
672
                              f'voxel_order of {exp_order}')
673 33
    return aff
674

675

676 33
@deprecate_with_version('aff_to_hdr is deprecated; please use the '
677
                        'nibabel.streamlines.TrkFile.affine_to_rasmm property, instead.',
678
                        since='2.5.0', until='4.0.0')
679 33
def aff_to_hdr(affine, trk_hdr, pos_vox=True, set_order=True):
680
    """ Set affine `affine` into trackvis header `trk_hdr`
681

682
    Affine is mapping from voxel space to Nifti RAS) output coordinate
683
    system convention; x: Left -> Right, y: Posterior -> Anterior, z:
684
    Inferior -> Superior.  Sets affine if possible, and voxel sizes, and voxel
685
    axis ordering.
686

687
    Parameters
688
    ----------
689
    affine : (4,4) array-like
690
       Affine voxel to mm transformation
691
    trk_hdr : mapping
692
       Mapping implementing __setitem__
693
    pos_vos : None or bool
694
        If None, currently defaults to False - this will change in future
695
        versions of nibabel.  If False, allow negative voxel sizes in header to
696
        record axis flips.  Negative voxels cause problems for trackvis (the
697
        application).  If True, enforce positive voxel sizes.
698
    set_order : None or bool
699
        If None, currently defaults to False - this will change in future
700
        versions of nibabel.  If False, do not set ``voxel_order`` field in
701
        `trk_hdr`.  If True, calculcate ``voxel_order`` from `affine` and set
702
        into `trk_hdr`.
703

704
    Returns
705
    -------
706
    None
707

708
    Notes
709
    -----
710
    version 2 of the trackvis header has a dedicated field for the nifti RAS
711
    affine. In theory trackvis 1 has enough information to store an affine,
712
    with the fields 'origin', 'voxel_size' and 'image_orientation_patient'.
713
    Unfortunately, to be able to store any affine, we'd need to be able to set
714
    negative voxel sizes, to encode axis flips. This is because
715
    'image_orientation_patient' is only two columns of the 3x3 rotation matrix,
716
    and we need to know the number of flips to reconstruct the third column
717
    reliably.  It turns out that negative flips upset trackvis (the
718
    application).  The application also ignores the origin field, and may not
719
    use the 'image_orientation_patient' field.
720
    """
721 33
    try:
722 33
        version = trk_hdr['version']
723 33
    except (KeyError, ValueError):  # dict or structured array
724 33
        version = 2
725 33
    if version == 2:
726 33
        trk_hdr['vox_to_ras'] = affine
727 33
    if set_order:
728 33
        trk_hdr['voxel_order'] = ''.join(aff2axcodes(affine))
729
    # Now on dodgy ground with DICOM fields in header
730
    # RAS to DPCS output
731 33
    affine = np.dot(DPCS_TO_TAL, affine)
732 33
    trans = affine[:3, 3]
733
    # Get zooms
734 33
    RZS = affine[:3, :3]
735 33
    zooms = np.sqrt(np.sum(RZS * RZS, axis=0))
736 33
    RS = RZS / zooms
737
    # If you said we could, adjust zooms to make RS correspond (below) to a
738
    # true rotation matrix.  We need to set the sign of one of the zooms to
739
    # deal with this.  Trackvis (the application) doesn't like negative zooms
740
    # at all, so you might want to disallow this with the pos_vox option.
741 33
    if not pos_vox and npl.det(RS) < 0:
742 33
        zooms[0] *= -1
743 33
        RS[:, 0] *= -1
744
    # retrieve rotation matrix from RS with polar decomposition.
745
    # Discard shears because we cannot store them.
746 33
    P, S, Qs = npl.svd(RS)
747 33
    R = np.dot(P, Qs)
748
    # it's an orthogonal matrix
749 33
    assert np.allclose(np.dot(R, R.T), np.eye(3))
750
    # set into header
751 33
    trk_hdr['origin'] = trans
752 33
    trk_hdr['voxel_size'] = zooms
753 33
    trk_hdr['image_orientation_patient'] = R[:, 0:2].T.ravel()
754

755

756 33
class TrackvisFileError(Exception):
757
    """ Error from TrackvisFile class
758
    """
759

760

761 33
class TrackvisFile(object):
762
    """ Convenience class to encapsulate trackvis file information
763

764
    Parameters
765
    ----------
766
    streamlines : sequence
767
       sequence of streamlines.  This object does not accept generic iterables
768
       as input because these can be consumed and make the object unusable.
769
       Please use the function interface to work with generators / iterables
770
    mapping : None or mapping
771
       Mapping defining header attributes
772
    endianness : {None, '<', '>'}
773
       Set here explicit endianness if required.  Endianness otherwise inferred
774
       from `streamlines`
775
    filename : None or str, optional
776
       filename
777
    points_space : {None, 'voxel', 'rasmm'}, optional
778
        Space in which streamline points are expressed in memory.  Default
779
        (None) means streamlines contain points in trackvis *voxmm* space
780
        (voxel positions * voxel sizes).  'voxel' means points are in voxel
781
        space (and need to be multiplied by voxel size for saving in file).
782
        'rasmm' mean the points are expressed in mm space according to the
783
        affine.  See ``read`` and ``write`` function docstrings for more
784
        detail.
785
    affine : None or (4,4) ndarray, optional
786
        Affine expressing relationship of voxels in an image to mm in RAS mm
787
        space. If 'points_space' is not None, you can use this to give the
788
        relationship between voxels, rasmm and voxmm space (above).
789
    """
790

791 33
    @deprecate_with_version('TrackvisFile is deprecated; please use '
792
                            'nibabel.streamlines.TrkFile, instead.',
793
                            since='2.5.0', until='4.0.0')
794 33
    def __init__(self,
795
                 streamlines,
796
                 mapping=None,
797
                 endianness=None,
798
                 filename=None,
799
                 points_space=None,
800
                 affine=None,
801
                 ):
802 33
        try:
803 33
            n_streams = len(streamlines)
804 33
        except TypeError:
805 33
            raise TrackvisFileError('Need sequence for streamlines input')
806 33
        self.streamlines = streamlines
807 33
        if endianness is None:
808 33
            if n_streams > 0:
809 33
                pts0 = streamlines[0][0]
810 33
                endianness = endian_codes[pts0.dtype.byteorder]
811
            else:
812 33
                endianness = native_code
813 33
        self.header = _hdr_from_mapping(None, mapping, endianness)
814 33
        self.endianness = endianness
815 33
        self.filename = filename
816 33
        self.points_space = points_space
817 33
        if affine is not None:
818 33
            self.set_affine(affine, pos_vox=True, set_order=True)
819

820 33
    @classmethod
821 33
    def from_file(klass, file_like, points_space=None):
822 33
        streamlines, header = read(file_like, points_space=points_space)
823 33
        filename = file_like if isinstance(file_like, str) else None
824 33
        return klass(streamlines, header, None, filename, points_space)
825

826 33
    def to_file(self, file_like):
827 33
        write(file_like, self.streamlines, self.header, self.endianness,
828
              points_space=self.points_space)
829 33
        self.filename = file_like if isinstance(file_like, str) else None
830

831 33
    def get_affine(self, atleast_v2=True):
832
        """ Get affine from header in object
833

834
        Returns
835
        -------
836
        aff : (4,4) ndarray
837
            affine from header
838
        atleast_v2 : None or bool, optional
839
            See ``aff_from_hdr`` docstring for detail.  If True, require valid
840
            affine in ``vox_to_ras`` field of header.
841

842
        Notes
843
        -----
844
        This method currently works for trackvis version 1 headers, but we
845
        consider it unsafe for version 1 headers, and in future versions of
846
        nibabel we will raise an error for trackvis headers < version 2.
847
        """
848 33
        return aff_from_hdr(self.header, atleast_v2)
849

850 33
    def set_affine(self, affine, pos_vox=True, set_order=True):
851
        """ Set affine `affine` into trackvis header
852

853
        Affine is mapping from voxel space to Nifti RAS) output coordinate
854
        system convention; x: Left -> Right, y: Posterior -> Anterior, z:
855
        Inferior -> Superior.  Sets affine if possible, and voxel sizes, and
856
        voxel axis ordering.
857

858
        Parameters
859
        ----------
860
        affine : (4,4) array-like
861
            Affine voxel to mm transformation
862
        pos_vos : None or bool, optional
863
            If None, currently defaults to False - this will change in future
864
            versions of nibabel.  If False, allow negative voxel sizes in
865
            header to record axis flips.  Negative voxels cause problems for
866
            trackvis (the application).  If True, enforce positive voxel sizes.
867
        set_order : None or bool, optional
868
            If None, currently defaults to False - this will change in future
869
            versions of nibabel.  If False, do not set ``voxel_order`` field in
870
            `trk_hdr`.  If True, calculcate ``voxel_order`` from `affine` and
871
            set into `trk_hdr`.
872

873
        Returns
874
        -------
875
        None
876
        """
877 33
        return aff_to_hdr(affine, self.header, pos_vox, set_order)

Read our documentation on viewing source code .

Loading