How to use the nibabel.trackvis function in nibabel

gfa = np.where(np.isnan(gfa), 0., gfa)
    ttc = ThresholdTissueClassifier(gfa, .2)

    # Create around N seeds
    seeds = utils.seeds_from_mask(gfa > .25, 2, affine=affine)
    seeds = seeds[::len(seeds) // N + 1]

    # Create streamline generator
    streamlines = LocalTracking(csapeaks, ttc, seeds, affine, .5, max_cross=1)
    trk_streamlines = utils.move_streamlines(streamlines,
                                             input_space=affine,
                                             output_space=trackvis_affine)
    trk = ((streamline, None, None) for streamline in trk_streamlines)

    # Save streamlines
    nib.trackvis.write(filename, trk, hdr)

import sys
import os.path as op
from fos import *
import fos.util

import numpy as np

from PySide.QtGui import QApplication

import nibabel as nib
a=nib.trackvis.read( op.join(op.dirname(__file__), "data", "tracks300.trk") )
g=np.array(a[0], dtype=np.object)
trk = [tr[0] for tr in a[0]]
#trk = trk[:10]
#g=np.array(trk, dtype=np.object)
#g=g[:200]

def prepare_tracks(trk):
    pos = []
    con = []
    cons = [] # an id for each fiber!
    offset = 0
    for i, f in enumerate(trk):
        fiblen = len(f)
        conarr = np.vstack( (np.array(range(fiblen - 1)), np.array(range(1,fiblen)) )).T.ravel()
        conarr += offset
        con.append( conarr )

def bench_quickbundles():
    dtype = "float32"
    repeat = 10
    nb_points = 12

    streams, hdr = nib.trackvis.read(get_fnames('fornix'))
    fornix = [s[0].astype(dtype) for s in streams]
    fornix = streamline_utils.set_number_of_points(fornix, nb_points)

    # Create eight copies of the fornix to be clustered (one in each octant).
    streamlines = []
    streamlines += [s + np.array([100, 100, 100], dtype) for s in fornix]
    streamlines += [s + np.array([100, -100, 100], dtype) for s in fornix]
    streamlines += [s + np.array([100, 100, -100], dtype) for s in fornix]
    streamlines += [s + np.array([100, -100, -100], dtype) for s in fornix]
    streamlines += [s + np.array([-100, 100, 100], dtype) for s in fornix]
    streamlines += [s + np.array([-100, -100, 100], dtype) for s in fornix]
    streamlines += [s + np.array([-100, 100, -100], dtype) for s in fornix]
    streamlines += [s + np.array([-100, -100, -100], dtype) for s in fornix]

    # The expected number of clusters of the fornix using threshold=10 is 4.
    threshold = 10.

def save_fibers(oldhdr, oldfib, fname, indices):
    """ Stores a new trackvis file fname using only given indices """

    hdrnew = oldhdr.copy()

    outstreams = []
    for i in indices:
        outstreams.append( oldfib[i] )

    n_fib_out = len(outstreams)
    hdrnew['n_count'] = n_fib_out

    log.info("Writing final no orphan fibers: %s" % fname)
    nibabel.trackvis.write(fname, outstreams, hdrnew)

def _run_interface(self, runtime):    
        filename = os.path.basename(self.inputs.gibbs_output)
    
        dx, dy, dz = get_data_dims(self.inputs.gibbs_input)
        vx, vy, vz = get_vox_dims(self.inputs.gibbs_input)
        image_file = nib.load(self.inputs.gibbs_input)
        affine = image_file.get_affine()
        import numpy as np
        #Reads MITK tracks
        fib, hdr = nib.trackvis.read(self.inputs.gibbs_output)
        trk_header = nib.trackvis.empty_header()
        trk_header['dim'] = [dx,dy,dz]
        trk_header['voxel_size'] = [vx,vy,vz]
        trk_header['origin'] = [0 ,0 ,0]
        axcode = nib.orientations.aff2axcodes(affine)
        if axcode[0] != str(hdr['voxel_order'])[0]:
            flip_x = -1
        else:
            flip_x = 1
        if axcode[1] != str(hdr['voxel_order'])[1]:
            flip_y = -1
        else:
            flip_y = 1
        if axcode[2] != str(hdr['voxel_order'])[2]:
            flip_z = -1
        else:

def save_fibers(oldhdr, oldfib, fname, indices):
    """ Stores a new trackvis file fname using only given indices """

    hdrnew = oldhdr.copy()

    outstreams = []
    for i in indices:
        outstreams.append( oldfib[i] )

    n_fib_out = len(outstreams)
    hdrnew['n_count'] = n_fib_out

    print("Writing final no orphan fibers: %s" % fname)
    nibabel.trackvis.write(fname, outstreams, hdrnew)

import nibabel as nib

hdr = nib.trackvis.empty_header()
hdr['voxel_size'] = (2., 2., 2.)
hdr['voxel_order'] = 'LAS'
hdr['dim'] = csapeaks.gfa.shape[:3]

"""
Save the streamlines.
"""

csa_streamlines_trk = ((sl, None, None) for sl in csa_streamlines)

csa_sl_fname = 'csa_prob_streamline.trk'

nib.trackvis.write(csa_sl_fname, csa_streamlines_trk, hdr)

"""
Visualize the streamlines with fvtk (python vtk is required).
"""

from dipy.viz import fvtk
from dipy.viz.colormap import line_colors

r = fvtk.ren()

fvtk.add(r, fvtk.line(csa_streamlines, line_colors(csa_streamlines)))

print('Saving illustration as csa_prob_tracks.png')

fvtk.record(r, n_frames=1, out_path='csa_prob_tracks.png', size=(600, 600))

trk_header['voxel_size'] = [r_vx, r_vy, r_vz]

            affine = np.dot(affine, np.diag(1. / np.array([vx, vy, vz, 1])))
            transformed_streamlines = transform_to_affine(
                streamlines, trk_header, affine)

            aff = affine_from_fsl_mat_file(xfm, [vx, vy, vz],
                                           [r_vx, r_vy, r_vz])
            iflogger.info(aff)

            axcode = aff2axcodes(reg_affine)
            trk_header['voxel_order'] = axcode[0] + axcode[1] + axcode[2]

            final_streamlines = move_streamlines(transformed_streamlines, aff)
            trk_tracks = ((ii, None, None) for ii in final_streamlines)
            trk.write(out_filename, trk_tracks, trk_header)
            iflogger.info('Saving transformed Trackvis file as %s',
                          out_filename)
            iflogger.info('New TrackVis Header:')
            iflogger.info(trk_header)
        else:
            iflogger.info(
                'Applying transformation from scanner coordinates to %s',
                self.inputs.image_file)
            axcode = aff2axcodes(affine)
            trk_header['voxel_order'] = axcode[0] + axcode[1] + axcode[2]
            trk_header['vox_to_ras'] = affine
            transformed_streamlines = transform_to_affine(
                streamlines, trk_header, affine)
            trk_tracks = ((ii, None, None) for ii in transformed_streamlines)
            trk.write(out_filename, trk_tracks, trk_header)
            iflogger.info('Saving Trackvis file as %s', out_filename)

def compute_length_array(trkfile=None, streams=None, savefname = 'lengths.npy'):
    if streams is None and not trkfile is None:
        log.info("Compute length array for fibers in %s" % trkfile)
        streams, hdr = tv.read(trkfile, as_generator = True)
        n_fibers = hdr['n_count']
        if n_fibers == 0:
            msg = "Header field n_count of trackfile %s is set to 0. No track seem to exist in this file." % trkfile
            log.error(msg)
            raise Exception(msg)
    else:
        n_fibers = len(streams)
        
    leng = np.zeros(n_fibers, dtype = np.float)
    for i,fib in enumerate(streams):
        leng[i] = util.length(fib[0])
    
    # store length array
    lefname = op.join(gconf.get_cmp_fibers(), savefname)
    np.save(lefname, leng)
    log.info("Store lengths array to: %s" % lefname)

del streams#,hdr

if not os.path.isfile(C_fname):

    print 'Starting LARCH ...'
    tim=time.clock()
    C,atracks=tl.larch(tracks,[50.**2,20.**2,5.**2],True,True)
    #tracks=[tm.downsample(t,3) for t in tracks]
    #C=pf.local_skeleton_clustering(tracks,20.)
    print 'Done in total of ',time.clock()-tim,'seconds.'

    print 'Saving result...'
    pkl.save_pickle(C_fname,C)
    
    streams=[(i,None,None)for i in atracks]
    tv.write(appr_fname,streams,hdr)

else:

    print 'Loading result...'
    C=pkl.load_pickle(C_fname)

skel=[]
for c in C:
    skel.append(C[c]['repz'])
    
print 'Showing dataset after clustering...'
r=fos.ren()
fos.clear(r)
colors=np.zeros((len(skel),3))
for (i,s) in enumerate(skel):