''' Testing trackvis module '''
from __future__ import division, print_function, absolute_import

from functools import partial

import numpy as np

from ..externals.six import BytesIO
from .. import trackvis as tv
from ..orientations import aff2axcodes
from ..volumeutils import native_code, swapped_code
from ..checkwarns import ErrorWarnings, IgnoreWarnings

from nose.tools import assert_true, assert_false, assert_equal, assert_raises

from numpy.testing import assert_array_equal, assert_array_almost_equal


def test_write():
    streams = []
    out_f = BytesIO()
    tv.write(out_f, [], {})
    assert_equal(out_f.getvalue(), tv.empty_header().tostring())
    out_f.truncate(0); out_f.seek(0)
    # Write something not-default
    tv.write(out_f, [], {'id_string':'TRACKb'})
    # read it back
    out_f.seek(0)
    streams, hdr = tv.read(out_f)
    assert_equal(hdr['id_string'], b'TRACKb')
    # check that we can pass none for the header
    out_f.truncate(0); out_f.seek(0)
    tv.write(out_f, [])
    out_f.truncate(0); out_f.seek(0)
    tv.write(out_f, [], None)
    # check that we check input values
    out_f.truncate(0); out_f.seek(0)
    assert_raises(tv.HeaderError,
           tv.write, out_f, [],{'id_string':'not OK'})
    assert_raises(tv.HeaderError,
           tv.write, out_f, [],{'version': 3})
    assert_raises(tv.HeaderError,
           tv.write, out_f, [],{'hdr_size': 0})


def test_write_scalars_props():
    # Test writing of scalar array with streamlines
    N = 6
    M = 2
    P = 4
    points = np.arange(N*3).reshape((N,3))
    scalars = np.arange(N*M).reshape((N,M)) + 100
    props = np.arange(P) + 1000
    # If scalars not same size for each point, error
    out_f = BytesIO()
    streams = [(points, None, None),
               (points, scalars, None)]
    assert_raises(tv.DataError, tv.write, out_f, streams)
    out_f.seek(0)
    streams = [(points, np.zeros((N,M+1)), None),
               (points, scalars, None)]
    assert_raises(tv.DataError, tv.write, out_f, streams)
    # Or if scalars different N compared to points
    bad_scalars = np.zeros((N+1,M))
    out_f.seek(0)
    streams = [(points, bad_scalars, None),
               (points, bad_scalars, None)]
    assert_raises(tv.DataError, tv.write, out_f, streams)
    # Similarly properties must have the same length for each streamline
    out_f.seek(0)
    streams = [(points, scalars, None),
               (points, scalars, props)]
    assert_raises(tv.DataError, tv.write, out_f, streams)
    out_f.seek(0)
    streams = [(points, scalars, np.zeros((P+1,))),
               (points, scalars, props)]
    assert_raises(tv.DataError, tv.write, out_f, streams)
    # If all is OK, then we get back what we put in
    out_f.seek(0)
    streams = [(points, scalars, props),
               (points, scalars, props)]
    tv.write(out_f, streams)
    out_f.seek(0)
    back_streams, hdr = tv.read(out_f)
    for actual, expected in zip(streams, back_streams):
        for a_el, e_el in zip(actual, expected):
            assert_array_equal(a_el, e_el)
    # Also so if the datatype of points, scalars is already float32 (github
    # issue #53)
    out_f.seek(0)
    streams = [(points.astype('f4'),
                scalars.astype('f4'),
                props.astype('f4'))]
    tv.write(out_f, streams)
    out_f.seek(0)
    back_streams, hdr = tv.read(out_f)
    for actual, expected in zip(streams, back_streams):
        for a_el, e_el in zip(actual, expected):
            assert_array_almost_equal(a_el, e_el)


def streams_equal(stream1, stream2):
    if not np.all(stream1[0] == stream2[0]):
        return False
    if stream1[1] is None:
        if not stream2[1] is None:
            return False
    if stream1[2] is None:
        if not stream2[2] is None:
            return False
    if not np.all(stream1[1] == stream2[1]):
        return False
    if not np.all(stream1[2] == stream2[2]):
        return False
    return True


def streamlist_equal(streamlist1, streamlist2):
    if len(streamlist1) != len(streamlist2):
        return False
    for s1, s2 in zip(streamlist1, streamlist2):
        if not streams_equal(s1, s2):
            return False
    return True


def test_round_trip():
    out_f = BytesIO()
    xyz0 = np.tile(np.arange(5).reshape(5,1), (1, 3))
    xyz1 = np.tile(np.arange(5).reshape(5,1) + 10, (1, 3))
    streams = [(xyz0, None, None), (xyz1, None, None)]
    tv.write(out_f, streams, {})
    out_f.seek(0)
    streams2, hdr = tv.read(out_f)
    assert_true(streamlist_equal(streams, streams2))
    # test that we can write in different endianness and get back same result,
    # for versions 1, 2 and not-specified
    for in_dict, back_version in (({},2),
                                  ({'version':2}, 2),
                                  ({'version':1}, 1)):
        for endian_code in (native_code, swapped_code):
            out_f.seek(0)
            tv.write(out_f, streams, in_dict, endian_code)
            out_f.seek(0)
            streams2, hdr = tv.read(out_f)
            assert_true(streamlist_equal(streams, streams2))
            assert_equal(hdr['version'], back_version)
    # test that we can get out and pass in generators
    out_f.seek(0)
    streams3, hdr = tv.read(out_f, as_generator=True)
    # check this is a generator rather than a list
    assert_true(hasattr(streams3, 'send'))
    # but that it results in the same output
    assert_true(streamlist_equal(streams, list(streams3)))
    # write back in
    out_f.seek(0)
    streams3, hdr = tv.read(out_f, as_generator=True)
    # Now we need a new file object, because we're still using the old one for
    # our generator
    out_f_write = BytesIO()
    tv.write(out_f_write, streams3, {})
    # and re-read just to check
    out_f_write.seek(0)
    streams2, hdr = tv.read(out_f_write)
    assert_true(streamlist_equal(streams, streams2))


def test_points_processing():
    # We may need to process points if they are in voxel or mm format
    out_f = BytesIO()
    def _rt(streams, hdr, points_space):
        # run round trip through IO object
        out_f.seek(0)
        tv.write(out_f, streams, hdr, points_space=points_space)
        out_f.seek(0)
        res0 = tv.read(out_f)
        out_f.seek(0)
        return res0, tv.read(out_f, points_space=points_space)
    n_pts = 5
    ijk0 = np.arange(n_pts * 3).reshape((n_pts,3)) / 2.0
    ijk1 = ijk0 + 20
    # Check with and without some scalars
    for scalars in ((None, None),
                    (np.arange(n_pts)[:, None],
                     np.arange(n_pts)[:, None] + 99)):
        vx_streams = [(ijk0, scalars[0], None), (ijk1, scalars[1], None)]
        vxmm_streams = [(ijk0 * [[2,3,4]], scalars[0], None),
                        (ijk1 * [[2,3,4]], scalars[1], None)]
        # voxmm is the default.  In this case we don't do anything to the
        # points, and we let the header pass through without further checks
        (raw_streams, hdr), (proc_streams, _) = _rt(vxmm_streams, {}, None)
        assert_true(streamlist_equal(raw_streams, proc_streams))
        assert_true(streamlist_equal(vxmm_streams, proc_streams))
        (raw_streams, hdr), (proc_streams, _) = _rt(vxmm_streams, {}, 'voxmm')
        assert_true(streamlist_equal(raw_streams, proc_streams))
        assert_true(streamlist_equal(vxmm_streams, proc_streams))
        # with 'voxels' as input, check for not all voxel_size == 0, warn if any
        # voxel_size == 0
        for hdr in ( # these cause read / write errors
                    # empty header has 0 voxel sizes
                    {},
                    {'voxel_size': [0,0,0]}, # the default
                    {'voxel_size': [-2,3,4]}, # negative not valid
                ):
            # Check error on write
            out_f.seek(0)
            assert_raises(tv.HeaderError,
                        tv.write, out_f, vx_streams, hdr, None, 'voxel')
            out_f.seek(0)
            # bypass write error and check read
            tv.write(out_f, vxmm_streams, hdr, None, points_space = None)
            out_f.seek(0)
            assert_raises(tv.HeaderError, tv.read, out_f, False, 'voxel')
        # There's a warning for any voxel sizes == 0
        hdr = {'voxel_size': [2, 3, 0]}
        with ErrorWarnings():
            assert_raises(UserWarning, _rt, vx_streams, hdr, 'voxel')
        # This should be OK
        hdr = {'voxel_size': [2, 3, 4]}
        (raw_streams, hdr), (proc_streams, _) = _rt(vx_streams, hdr, 'voxel')
        assert_true(streamlist_equal(vxmm_streams, raw_streams))
        assert_true(streamlist_equal(vx_streams, proc_streams))
        # Now we try with rasmm points.  In this case we need valid voxel_size,
        # and voxel_order, and vox_to_ras.  The voxel_order has to match the
        # vox_to_ras, and so do the voxel sizes
        aff = np.diag([2,3,4,1])
        # In this case the trk -> vx and vx -> mm invert each other
        rasmm_streams = vxmm_streams
        for hdr in ( # all these cause read and write errors for rasmm
            # Empty header has no valid affine
            {},
            # Error if ras_to_mm not defined (as in version 1)
            {'voxel_size': [2, 3, 4], 'voxel_order': 'RAS', 'version':1},
            # or it's all zero
            {'voxel_size': [2, 3, 4], 'voxel_order': 'RAS',
            'vox_to_ras': np.zeros((4,4))},
            # as it is by default
            {'voxel_size': [2, 3, 4], 'voxel_order': 'RAS'},
            # or the voxel_size doesn't match the affine
            {'voxel_size': [2, 2, 4], 'voxel_order': 'RAS',
            'vox_to_ras': aff},
            # or the voxel_order doesn't match the affine
            {'voxel_size': [2, 3, 4], 'voxel_order': 'LAS',
            'vox_to_ras': aff},
            ):
            # Check error on write
            out_f.seek(0)
            assert_raises(tv.HeaderError,
                        tv.write, out_f, rasmm_streams, hdr, None, 'rasmm')
            out_f.seek(0)
            # bypass write error and check read
            tv.write(out_f, vxmm_streams, hdr, None, points_space = None)
            out_f.seek(0)
            assert_raises(tv.HeaderError, tv.read, out_f, False, 'rasmm')
        # This should be OK
        hdr = {'voxel_size': [2, 3, 4], 'voxel_order': 'RAS',
            'vox_to_ras': aff}
        (raw_streams, hdr), (proc_streams, _) = _rt(rasmm_streams, hdr, 'rasmm')
        assert_true(streamlist_equal(vxmm_streams, raw_streams))
        assert_true(streamlist_equal(rasmm_streams, proc_streams))
        # More complex test to check matrix orientation
        fancy_affine = np.array([[0., -2, 0, 10],
                                [3, 0, 0, 20],
                                [0, 0, 4, 30],
                                [0, 0, 0, 1]])
        hdr = {'voxel_size': [3, 2, 4], 'voxel_order': 'ALS',
            'vox_to_ras': fancy_affine}
        def f(pts): # from vx to mm
            pts = pts[:,[1,0,2]] * [[-2,3,4]] # apply zooms / reorder
            return pts + [[10,20,30]] # apply translations
        xyz0, xyz1 = f(ijk0), f(ijk1)
        fancy_rasmm_streams = [(xyz0, scalars[0], None),
                               (xyz1, scalars[1], None)]
        fancy_vxmm_streams = [(ijk0 * [[3,2,4]], scalars[0], None),
                              (ijk1 * [[3,2,4]], scalars[1], None)]
        (raw_streams, hdr), (proc_streams, _) = _rt(
            fancy_rasmm_streams, hdr, 'rasmm')
        assert_true(streamlist_equal(fancy_vxmm_streams, raw_streams))
        assert_true(streamlist_equal(fancy_rasmm_streams, proc_streams))


def test__check_hdr_points_space():
    # Test checking routine for points_space input given header
    # None or voxmm -> no checks, pass through
    assert_equal(tv._check_hdr_points_space({}, None), None)
    assert_equal(tv._check_hdr_points_space({}, 'voxmm'), None)
    # strange value for points_space -> ValueError
    assert_raises(ValueError,
                  tv._check_hdr_points_space, {}, 'crazy')
    # Input not in (None, 'voxmm', 'voxels', 'rasmm') - error
    # voxels means check voxel sizes present and not all 0.
    hdr = tv.empty_header()
    assert_array_equal(hdr['voxel_size'], [0,0,0])
    assert_raises(tv.HeaderError,
                  tv._check_hdr_points_space, hdr, 'voxel')
    # Negative voxel size gives error - because it is not what trackvis does,
    # and this not what we mean by 'voxmm'
    hdr['voxel_size'] = [-2, 3, 4]
    assert_raises(tv.HeaderError,
                  tv._check_hdr_points_space, hdr, 'voxel')
    # Warning here only
    hdr['voxel_size'] = [2, 3, 0]
    with ErrorWarnings():
        assert_raises(UserWarning,
                      tv._check_hdr_points_space, hdr, 'voxel')
    # This is OK
    hdr['voxel_size'] = [2, 3, 4]
    assert_equal(tv._check_hdr_points_space(hdr, 'voxel'), None)
    # rasmm - check there is an affine, that it matches voxel_size and
    # voxel_order
    # no affine
    hdr['voxel_size'] = [2, 3, 4]
    assert_raises(tv.HeaderError,
                  tv._check_hdr_points_space, hdr, 'rasmm')
    # still no affine
    hdr['voxel_order'] = 'RAS'
    assert_raises(tv.HeaderError,
                  tv._check_hdr_points_space, hdr, 'rasmm')
    # nearly an affine, but 0 at position 3,3 - means not recorded in trackvis
    # standard
    hdr['vox_to_ras'] = np.diag([2,3,4,0])
    assert_raises(tv.HeaderError,
                  tv._check_hdr_points_space, hdr, 'rasmm')
    # This affine doesn't match RAS voxel order
    hdr['vox_to_ras'] = np.diag([-2,3,4,1])
    assert_raises(tv.HeaderError,
                  tv._check_hdr_points_space, hdr, 'rasmm')
    # This affine doesn't match the voxel size
    hdr['vox_to_ras'] = np.diag([3,3,4,1])
    assert_raises(tv.HeaderError,
                  tv._check_hdr_points_space, hdr, 'rasmm')
    # This should be OK
    good_aff = np.diag([2,3,4,1])
    hdr['vox_to_ras'] = good_aff
    assert_equal(tv._check_hdr_points_space(hdr, 'rasmm'),
                 None)
    # Default voxel order of LPS assumed
    hdr['voxel_order'] = ''
    # now the RAS affine raises an error
    assert_raises(tv.HeaderError,
                  tv._check_hdr_points_space, hdr, 'rasmm')
    # this affine does have LPS voxel order
    good_lps = np.dot(np.diag([-1,-1,1,1]), good_aff)
    hdr['vox_to_ras'] = good_lps
    assert_equal(tv._check_hdr_points_space(hdr, 'rasmm'),
                 None)


def test_empty_header():
    for endian in '<>':
        for version in (1, 2):
            hdr = tv.empty_header(endian, version)
            assert_equal(hdr['id_string'], b'TRACK')
            assert_equal(hdr['version'], version)
            assert_equal(hdr['hdr_size'], 1000)
            assert_array_equal(
                hdr['image_orientation_patient'],
                [0,0,0,0,0,0])
    hdr = tv.empty_header(version=2)
    assert_array_equal(hdr['vox_to_ras'], np.zeros((4,4)))
    hdr_endian = tv.endian_codes[tv.empty_header().dtype.byteorder]
    assert_equal(hdr_endian, tv.native_code)


def test_get_affine():
    # Test get affine behavior, including pending deprecation
    hdr = tv.empty_header()
    # Using version 1 affine is not a good idea because is fragile and not
    # very useful. The default atleast_v2=None mode raises a FutureWarning
    with ErrorWarnings():
        assert_raises(FutureWarning, tv.aff_from_hdr, hdr)
    # testing the old behavior
    old_afh = partial(tv.aff_from_hdr, atleast_v2=False)
    # default header gives useless affine
    assert_array_equal(old_afh(hdr),
                       np.diag([0,0,0,1]))
    hdr['voxel_size'] = 1
    assert_array_equal(old_afh(hdr),
                       np.diag([0,0,0,1]))
    # DICOM direction cosines
    hdr['image_orientation_patient'] = [1,0,0,0,1,0]
    assert_array_equal(old_afh(hdr),
                       np.diag([-1,-1,1,1]))
    # RAS direction cosines
    hdr['image_orientation_patient'] = [-1,0,0,0,-1,0]
    assert_array_equal(old_afh(hdr),
                       np.eye(4))
    # translations
    hdr['origin'] = [1,2,3]
    exp_aff = np.eye(4)
    exp_aff[:3,3] = [-1,-2,3]
    assert_array_equal(old_afh(hdr),
                       exp_aff)
    # check against voxel order.  This one works
    hdr['voxel_order'] = ''.join(aff2axcodes(exp_aff))
    assert_equal(hdr['voxel_order'], b'RAS')
    assert_array_equal(old_afh(hdr), exp_aff)
    # This one doesn't
    hdr['voxel_order'] = 'LAS'
    assert_raises(tv.HeaderError, old_afh, hdr)
    # This one does work because the routine allows the final dimension to
    # be flipped to try and match the voxel order
    hdr['voxel_order'] = 'RAI'
    exp_aff = exp_aff * [[1,1,-1,1]]
    assert_array_equal(old_afh(hdr), exp_aff)
    # Check round trip case for flipped and unflipped, when positive voxels
    # only allowed.  This checks that the flipping heuristic works.
    flipped_aff = exp_aff
    unflipped_aff = exp_aff * [1,1,-1,1]
    for in_aff, o_codes in ((unflipped_aff, b'RAS'),
                            (flipped_aff, b'RAI')):
        hdr = tv.empty_header()
        tv.aff_to_hdr(in_aff, hdr, pos_vox=True, set_order=True)
        # Unset easier option
        hdr['vox_to_ras'] = 0
        assert_equal(hdr['voxel_order'], o_codes)
        # Check it came back the way we wanted
        assert_array_equal(old_afh(hdr), in_aff)
    # Check that the default case matches atleast_v2=False case
    with IgnoreWarnings():
        assert_array_equal(tv.aff_from_hdr(hdr), flipped_aff)
    # now use the easier vox_to_ras field
    hdr = tv.empty_header()
    aff = np.eye(4)
    aff[:3,:] = np.arange(12).reshape(3,4)
    hdr['vox_to_ras'] = aff
    # Pass v2 flag explicitly to avoid warnings
    assert_array_equal(tv.aff_from_hdr(hdr, atleast_v2=False), aff)
    # mappings work too
    d = {'version': 1,
        'voxel_size': np.array([1,2,3]),
        'image_orientation_patient': np.array([1,0,0,0,1,0]),
        'origin': np.array([10,11,12])}
    aff = tv.aff_from_hdr(d, atleast_v2=False)


def test_aff_to_hdr():
    # The behavior is changing soon, change signaled by FutureWarnings
    # This is the call to get the old behavior
    old_a2h = partial(tv.aff_to_hdr, pos_vox=False, set_order=False)
    hdr = {'version': 1}
    affine = np.diag([1,2,3,1])
    affine[:3,3] = [10,11,12]
    old_a2h(affine, hdr)
    assert_array_almost_equal(tv.aff_from_hdr(hdr, atleast_v2=False), affine)
    # put flip into affine
    aff2 = affine.copy()
    aff2[:,2] *=-1
    old_a2h(aff2, hdr)
    # Historically we flip the first axis if there is a negative determinant
    assert_array_almost_equal(hdr['voxel_size'], [-1,2,3])
    assert_array_almost_equal(tv.aff_from_hdr(hdr, atleast_v2=False), aff2)
    # Test that default mode raises DeprecationWarning
    with ErrorWarnings():
        assert_raises(FutureWarning, tv.aff_to_hdr, affine, hdr)
        assert_raises(FutureWarning, tv.aff_to_hdr, affine, hdr, None, None)
        assert_raises(FutureWarning, tv.aff_to_hdr, affine, hdr, False, None)
        assert_raises(FutureWarning, tv.aff_to_hdr, affine, hdr, None, False)
    # And has same effect as above
    with IgnoreWarnings():
        tv.aff_to_hdr(affine, hdr)
    assert_array_almost_equal(tv.aff_from_hdr(hdr, atleast_v2=False), affine)
    # Check pos_vox and order flags
    for hdr in ({}, {'version':2}, {'version':1}):
        tv.aff_to_hdr(aff2, hdr, pos_vox=True, set_order=False)
        assert_array_equal(hdr['voxel_size'], [1, 2, 3])
        assert_false('voxel_order' in hdr)
        tv.aff_to_hdr(aff2, hdr, pos_vox=False, set_order=True)
        assert_array_equal(hdr['voxel_size'], [-1, 2, 3])
        assert_equal(hdr['voxel_order'], 'RAI')
        tv.aff_to_hdr(aff2, hdr, pos_vox=True, set_order=True)
        assert_array_equal(hdr['voxel_size'], [1, 2, 3])
        assert_equal(hdr['voxel_order'], 'RAI')
        if 'version' in hdr and hdr['version'] == 1:
            assert_false('vox_to_ras' in hdr)
        else:
            assert_array_equal(hdr['vox_to_ras'], aff2)


def test_tv_class():
    tvf = tv.TrackvisFile([])
    assert_equal(tvf.streamlines, [])
    assert_true(isinstance(tvf.header, np.ndarray))
    assert_equal(tvf.endianness, tv.native_code)
    assert_equal(tvf.filename, None)
    out_f = BytesIO()
    tvf.to_file(out_f)
    assert_equal(out_f.getvalue(), tv.empty_header().tostring())
    out_f.truncate(0); out_f.seek(0)
    # Write something not-default
    tvf = tv.TrackvisFile([], {'id_string':'TRACKb'})
    tvf.to_file(out_f)
    # read it back
    out_f.seek(0)
    tvf_back = tv.TrackvisFile.from_file(out_f)
    assert_equal(tvf_back.header['id_string'], b'TRACKb')
    # check that we check input values
    out_f.truncate(0); out_f.seek(0)
    assert_raises(tv.HeaderError,
                  tv.TrackvisFile,
                  [],{'id_string':'not OK'})
    assert_raises(tv.HeaderError,
                  tv.TrackvisFile,
                  [],{'version': 3})
    assert_raises(tv.HeaderError,
                  tv.TrackvisFile,
                  [],{'hdr_size':0})
    affine = np.diag([1,2,3,1])
    affine[:3,3] = [10,11,12]
    # affine methods will raise same warnings and errors as function
    with ErrorWarnings():
        assert_raises(FutureWarning, tvf.set_affine, affine)
        assert_raises(FutureWarning, tvf.set_affine, affine, None, None)
        assert_raises(FutureWarning, tvf.set_affine, affine, False, None)
        assert_raises(FutureWarning, tvf.set_affine, affine, None, False)
        assert_raises(FutureWarning, tvf.get_affine)
        assert_raises(FutureWarning, tvf.get_affine, None)
    tvf.set_affine(affine, pos_vox=True, set_order=True)
    aff = tvf.get_affine(atleast_v2=True)
    assert_array_almost_equal(aff, affine)
    # Test that we raise an error with an iterator
    assert_raises(tv.TrackvisFileError,
                  tv.TrackvisFile,
                  iter([]))


def test_tvfile_io():
    # Test reading and writing tracks with file class
    out_f = BytesIO()
    ijk0 = np.arange(15).reshape((5,3)) / 2.0
    ijk1 = ijk0 + 20
    vx_streams = [(ijk0, None, None), (ijk1, None, None)]
    vxmm_streams = [(ijk0 * [[2,3,4]], None, None),
                    (ijk1 * [[2,3,4]], None, None)]
    # Roundtrip basic
    tvf = tv.TrackvisFile(vxmm_streams)
    tvf.to_file(out_f)
    out_f.seek(0)
    tvf2 = tv.TrackvisFile.from_file(out_f)
    assert_equal(tvf2.filename, None)
    assert_true(streamlist_equal(vxmm_streams, tvf2.streamlines))
    assert_equal(tvf2.points_space, None)
    # Voxel points_space
    tvf = tv.TrackvisFile(vx_streams, points_space='voxel')
    out_f.seek(0)
    # No voxel size - error
    assert_raises(tv.HeaderError, tvf.to_file, out_f)
    out_f.seek(0)
    # With voxel size, no error, roundtrip works
    tvf.header['voxel_size'] = [2,3,4]
    tvf.to_file(out_f)
    out_f.seek(0)
    tvf2 = tv.TrackvisFile.from_file(out_f, points_space='voxel')
    assert_true(streamlist_equal(vx_streams, tvf2.streamlines))
    assert_equal(tvf2.points_space, 'voxel')
    out_f.seek(0)
    # Also with affine specified
    tvf = tv.TrackvisFile(vx_streams, points_space='voxel',
                          affine=np.diag([2,3,4,1]))
    tvf.to_file(out_f)
    out_f.seek(0)
    tvf2 = tv.TrackvisFile.from_file(out_f, points_space='voxel')
    assert_true(streamlist_equal(vx_streams, tvf2.streamlines))
    # Fancy affine test
    fancy_affine = np.array([[0., -2, 0, 10],
                             [3, 0, 0, 20],
                             [0, 0, 4, 30],
                             [0, 0, 0, 1]])
    def f(pts): # from vx to mm
        pts = pts[:,[1,0,2]] * [[-2,3,4]] # apply zooms / reorder
        return pts + [[10,20,30]] # apply translations
    xyz0, xyz1 = f(ijk0), f(ijk1)
    fancy_rasmm_streams = [(xyz0, None, None), (xyz1, None, None)]
    # Roundtrip
    tvf = tv.TrackvisFile(fancy_rasmm_streams, points_space='rasmm')
    out_f.seek(0)
    # No affine
    assert_raises(tv.HeaderError, tvf.to_file, out_f)
    out_f.seek(0)
    # With affine set, no error, roundtrip works
    tvf.set_affine(fancy_affine, pos_vox=True, set_order=True)
    tvf.to_file(out_f)
    out_f.seek(0)
    tvf2 = tv.TrackvisFile.from_file(out_f, points_space='rasmm')
    assert_true(streamlist_equal(fancy_rasmm_streams, tvf2.streamlines))
    assert_equal(tvf2.points_space, 'rasmm')
    out_f.seek(0)
    # Also when affine given in init
    tvf = tv.TrackvisFile(fancy_rasmm_streams, points_space='rasmm',
                          affine=fancy_affine)
    tvf.to_file(out_f)
    out_f.seek(0)
    tvf2 = tv.TrackvisFile.from_file(out_f, points_space='rasmm')
    assert_true(streamlist_equal(fancy_rasmm_streams, tvf2.streamlines))
