def bench_load_trk():
    rng = np.random.RandomState(42)
    dtype = 'float32'
    NB_STREAMLINES = 5000
    NB_POINTS = 1000
    points = [rng.rand(NB_POINTS, 3).astype(dtype)
              for i in range(NB_STREAMLINES)]
    scalars = [rng.rand(NB_POINTS, 10).astype(dtype)
               for i in range(NB_STREAMLINES)]

    repeat = 10

    with InTemporaryDirectory():
        trk_file = "tmp.trk"
        tractogram = Tractogram(points, affine_to_rasmm=np.eye(4))
        TrkFile(tractogram).save(trk_file)

        streamlines_old = [d[0] - 0.5
                           for d in tv.read(trk_file, points_space="rasmm")[0]]
        mtime_old = measure('tv.read(trk_file, points_space="rasmm")', repeat)
        print("Old: Loaded {:,} streamlines in {:6.2f}".format(NB_STREAMLINES,
                                                               mtime_old))

        trk = nib.streamlines.load(trk_file, lazy_load=False)
        streamlines_new = trk.streamlines
        mtime_new = measure('nib.streamlines.load(trk_file, lazy_load=False)',
                            repeat)
        print("\nNew: Loaded {:,} streamlines in {:6.2}".format(NB_STREAMLINES,
                                                                mtime_new))
        print("Speedup of {:.2f}".format(mtime_old / mtime_new))
        for s1, s2 in zip(streamlines_new, streamlines_old):
            assert_array_equal(s1, s2)

    # Points and scalars
    with InTemporaryDirectory():

        trk_file = "tmp.trk"
        tractogram = Tractogram(points,
                                data_per_point={'scalars': scalars},
                                affine_to_rasmm=np.eye(4))
        TrkFile(tractogram).save(trk_file)

        streamlines_old = [d[0] - 0.5
                           for d in tv.read(trk_file, points_space="rasmm")[0]]

        scalars_old = [d[1]
                       for d in tv.read(trk_file, points_space="rasmm")[0]]
        mtime_old = measure('tv.read(trk_file, points_space="rasmm")', repeat)
        msg = "Old: Loaded {:,} streamlines with scalars in {:6.2f}"
        print(msg.format(NB_STREAMLINES, mtime_old))

        trk = nib.streamlines.load(trk_file, lazy_load=False)
        scalars_new = trk.tractogram.data_per_point['scalars']
        mtime_new = measure('nib.streamlines.load(trk_file, lazy_load=False)',
                            repeat)
        msg = "New: Loaded {:,} streamlines with scalars in {:6.2f}"
        print(msg.format(NB_STREAMLINES, mtime_new))
        print("Speedup of {:2f}".format(mtime_old / mtime_new))
        for s1, s2 in zip(scalars_new, scalars_old):
            assert_array_equal(s1, s2)

def fg_from_trk(trk_file, affine=None):
    """
    Read data from a trackvis .trk file and create a FiberGroup object
    according to the information in it.
    """

    # Generate right away, since we're going to do it anyway:
    read_trk = tv.read(trk_file, as_generator=False)
    fibers_trk = read_trk[0]

    # Per default read from the affine from the file header:
    if affine is not None:
        aff = affine
    else: 
        hdr = read_trk[1]
        aff= tv.aff_from_hdr(hdr)
        # If the header contains a bogus affine, we revert to np.eye(4), so we
        # don't get into trouble later:
        try:
            np.matrix(aff).getI()
        except np.linalg.LinAlgError:
            e_s = "trk file contains bogus header, reverting to np.eye(4)" 
            warnings.warn(e_s)
            aff = np.eye(4)

    fibers = []
    for f in fibers_trk:
        fibers.append(ozf.Fiber(np.array(f[0]).T,affine=aff))

    return ozf.FiberGroup(fibers, affine=aff)

def use_camino(mask,maskname='maskA',iteration='20'):

    img=nib.Nifti1Image(mask.astype(np.uint8),np.array([[1,0,0,-31.5],[0,1,0,-31.5],[0,0,1,-32.5],[0,0,0,1]]))
    nib.save(img,'/tmp/'+maskname+'.nii.gz')

    cmd='track -inputfile /home/eg309/Data/orbital_phantoms/dwi_dir/workflow/tractography/_subject_id_subject1/picopdfs_twoten/data_fit_pdfs.Bdouble -seedfile /tmp/'+maskname+'.nii.gz -iterations '+iteration+' -numpds 2  -inputmodel pico -outputfile data_fit_pdfs_tracked'
    pipe(cmd)
    cmd='camino_to_trackvis -i ./data_fit_pdfs_tracked -o data_fit_pdfs_tracked.trk -l 30 -d 64,64,64 -x 1.0,1.0,1.0 --voxel-order LAS'
    pipe(cmd)

    fname='data_fit_pdfs_tracked.trk'
    streams,hdr=tv.read(fname,points_space=None)
    tracks=[s[0] for s in streams]
    
    shape=(64,64,64)
    msk=np.load('/tmp/allmasks.npy')
    #"""
    tracks=[t+np.array([32.,32,32]) for t in tracks]

    tracksA=count_tracks_mask(tracks,shape,msk,1)        
    tracksB=count_tracks_mask(tracks,shape,msk,2)
    tracksC=count_tracks_mask(tracks,shape,msk,3)    
    tracksD=count_tracks_mask(tracks,shape,msk,4)    

    print 'Initial', 'A', 'B', 'C', 'D'
    print len(tracks), len(tracksA), len(tracksB), len(tracksC), len(tracksD)
    #"""
    #show_tracks(tracks)
    alpha=1.
    lw=2.
    bg=(1.,1.,1.,1)    

    """

def tractography_from_trackvis_file(filename):
    tracts_and_data, header = trackvis.read(filename, points_space='rasmm')

    tracts, scalars, properties = list(zip(*tracts_and_data))

    scalar_names = [n for n in header['scalar_name'] if len(n) > 0]

    #scalar_names_unique = []
    #scalar_names_subcomp = {}
    # for sn in scalar_names:
    #    if re.match('.*_[0-9]{2}', sn):
    #        prefix = sn[:sn.rfind('_')]
    #        if prefix not in scalar_names_unique:
    #            scalar_names_unique.append(prefix)
    #            scalar_names_subcomp[prefix] = int(sn[-2:])
    #        scalar_names_subcomp[prefix] = max(sn[-2:], scalar_names_subcomp[prefix])
    #    else:
    #        scalar_names_unique.append(sn)

    tracts_data = {}
    for i, sn in enumerate(scalar_names):
        if hasattr(sn, 'decode'):
            sn = sn.decode()
        tracts_data[sn] = [scalar[:, i][:, None] for scalar in scalars]

    affine = header['vox_to_ras']
    image_dims = header['dim']

    tr = Tractography(
        tracts, tracts_data,
        affine=affine, image_dims=image_dims
    )

    return tr

    def launch(self, data_file, region_volume=None):
        datatype = self._base_before_launch(data_file, region_volume)

        # note the streaming parsing, we do not load the dataset in memory at once
        tract_gen, hdr = trackvis.read(data_file, as_generator=True)

        vox2ras = _SpaceTransform(hdr)
        tract_start_indices = [0]
        tract_region = []

        # we process tracts in bigger chunks to optimize disk write costs
        for tract_bundle in chunk_iter(tract_gen, self.READ_CHUNK):
            tract_bundle = [tr[0] for tr in tract_bundle]

            for tr in tract_bundle:
                tract_start_indices.append(tract_start_indices[-1] + len(tr))
                if region_volume is not None:
                    tract_region.append(self._get_tract_region(tr[0]))

            vertices = numpy.concatenate(tract_bundle) # in voxel space
            vertices = vox2ras.transform(vertices)

            datatype.store_data_chunk("vertices", vertices, grow_dimension=0, close_file=False)

        datatype.tract_start_idx = tract_start_indices
        datatype.tract_region = numpy.array(tract_region, dtype=numpy.int16)
        return datatype

	def _run_interface(self, runtime):
		tracks, header = trk.read(self.inputs.in_file)
		if not isdefined(self.inputs.data_dims):
			data_dims = header['dim']
		else:
			data_dims = self.inputs.data_dims

		if not isdefined(self.inputs.voxel_dims):
			voxel_size = header['voxel_size']
		else:
			voxel_size = self.inputs.voxel_dims

		affine = header['vox_to_ras']

		streams = ((ii[0]) for ii in tracks)
		data = density_map(streams, data_dims, voxel_size)
		if data.max() < 2**15:
		   data = data.astype('int16')

		img = nb.Nifti1Image(data,affine)
		out_file = op.abspath(self.inputs.out_filename)
		nb.save(img, out_file)
		iflogger.info('Track density map saved as {i}'.format(i=out_file))
		iflogger.info('Data Dimensions {d}'.format(d=data_dims))
		iflogger.info('Voxel Dimensions {v}'.format(v=voxel_size))
		return runtime

def load_tractography(filename_tractography ):
    
  
    
    tractography, header = trackvis.read(filename_tractography)
    tractography = [streamline[0] for streamline in tractography]
    return  tractography

def get_streamlines():
    from nibabel import trackvis as tv
    from dipy.data import get_data

    fname = get_data('fornix')
    streams, hdr = tv.read(fname)
    streamlines = [i[0] for i in streams]
    return streamlines

def show(fname):

    streams, hdr = tv.read(fname)
    streamlines = [s[0] for s in streams]

    renderer = fvtk.ren() 
    fvtk_tubes = vtk_a.line(streamlines, opacity=0.2, linewidth=5)
    fvtk.add(renderer, fvtk_tubes)
    fvtk.show(renderer)

def bench_compress_streamlines():
    repeat = 10
    fname = get_data('fornix')
    streams, hdr = tv.read(fname)
    streamlines = [i[0] for i in streams]

    print("Timing compress_streamlines() in Cython ({0} streamlines)".format(len(streamlines)))
    cython_time = measure("compress_streamlines(streamlines)", repeat)
    print("Cython time: {0:.3}sec".format(cython_time))
    del streamlines

    fname = get_data('fornix')
    streams, hdr = tv.read(fname)
    streamlines = [i[0] for i in streams]
    python_time = measure("map(compress_streamlines_python, streamlines)", repeat)
    print("Python time: {0:.2}sec".format(python_time))
    print("Speed up of {0}x".format(python_time/cython_time))
    del streamlines

	def _run_interface(self, runtime):
		# Loading the ROI file
	    import nibabel as nib
	    import numpy as np
	    from dipy.tracking import utils

	    img = nib.load(self.inputs.ROI_file)
	    data = img.get_data()
	    affine = img.get_affine()

	    # Getting the FA file
	    img = nib.load(self.inputs.FA_file)
	    FA_data = img.get_data()
	    FA_affine = img.get_affine()

	    # Loading the streamlines
	    from nibabel import trackvis
	    streams, hdr = trackvis.read(self.inputs.trackfile,points_space='rasmm')
	    streamlines = [s[0] for s in streams]
	    streamlines_affine = trackvis.aff_from_hdr(hdr,atleast_v2=True)

	    # Checking for negative values
	    from dipy.tracking._utils import _mapping_to_voxel, _to_voxel_coordinates
	    endpoints = [sl[0::len(sl)-1] for sl in streamlines]
	    lin_T, offset = _mapping_to_voxel(affine, (1.,1.,1.))
	    inds = np.dot(endpoints, lin_T)
	    inds += offset
	    negative_values = np.where(inds <0)[0]
	    for negative_value in sorted(negative_values, reverse=True):
			del streamlines[negative_value]

	    # Constructing the streamlines matrix
	    matrix,mapping = utils.connectivity_matrix(streamlines=streamlines,label_volume=data,affine=streamlines_affine,symmetric=True,return_mapping=True,mapping_as_streamlines=True)
	    matrix[matrix < 10] = 0

	    # Constructing the FA matrix
	    dimensions = matrix.shape
	    FA_matrix = np.empty(shape=dimensions)

	    for i in range(0,dimensions[0]):
	        for j in range(0,dimensions[1]):
	            if matrix[i,j]:
	                dm = utils.density_map(mapping[i,j], FA_data.shape, affine=streamlines_affine)
            		FA_matrix[i,j] = np.mean(FA_data[dm>0])
	            else:
	                FA_matrix[i,j] = 0

	    FA_matrix[np.tril_indices(n=len(FA_matrix))] = 0
	    FA_matrix = FA_matrix.T + FA_matrix - np.diagonal(FA_matrix)

	    from nipype.utils.filemanip import split_filename
	    _, base, _ = split_filename(self.inputs.trackfile)
	    np.savetxt(base + '_FA_matrix.txt',FA_matrix,delimiter='\t')
	    return runtime

    def __init__(self, filename):
        """Load the TrackVis file """
        self._filename = filename

        self.tracks, self.fiberHeader = niv.read(filename)
        self.fiberCount = self.fiberHeader['n_count'].ravel()[0]
        
        print "number of fibers ->", self.fiberCount
        self.currentFiber = 0
        self.shape = tuple(self.fiberHeader['dim'])

        # Position at first fiber
        self._rewind()

def load_tracks(method="pmt"):
    from nibabel import trackvis as tv

    dname = "/home/eg309/Data/orbital_phantoms/dwi_dir/subject1/"

    if method == "pmt":
        fname = "/home/eg309/Data/orbital_phantoms/dwi_dir/workflow/tractography/_subject_id_subject1/cam2trk_pico_twoten/data_fit_pdfs_tracked.trk"
        streams, hdr = tv.read(fname, points_space="voxel")
        tracks = [s[0] for s in streams]
    if method == "dti":
        fname = dname + "dti_tracks.dpy"
    if method == "dsi":
        fname = dname + "dsi_tracks.dpy"
    if method == "gqs":
        fname = dname + "gqi_tracks.dpy"
    if method == "eit":
        fname = dname + "eit_tracks.dpy"
    if method in ["dti", "dsi", "gqs", "eit"]:
        dpr_linear = Dpy(fname, "r")
        tracks = dpr_linear.read_tracks()
        dpr_linear.close()

    if method != "pmt":
        tracks = [t - np.array([96 / 2.0, 96 / 2.0, 55 / 2.0]) for t in tracks if track_range(t, 100 / 2.5, 150 / 2.5)]
    tracks = [t for t in tracks if track_range(t, 100 / 2.5, 150 / 2.5)]

    print "final no of tracks ", len(tracks)
    qb = QuickBundles(tracks, 25.0 / 2.5, 18)
    # from dipy.viz import fvtk
    # r=fvtk.ren()
    # fvtk.add(r,fvtk.line(qb.virtuals(),fvtk.red))
    # fvtk.show(r)
    # show_tracks(tracks)#qb.exemplars()[0])
    # qb.remove_small_clusters(40)
    del tracks
    # load
    tl = TrackLabeler(qb, qb.downsampled_tracks(), vol_shape=None, tracks_line_width=3.0, tracks_alpha=1)

    # return tracks
    w = World()
    w.add(tl)
    # create window
    wi = Window(caption="Fos", bgcolor=(1.0, 1.0, 1.0, 1.0), width=1600, height=900)
    wi.attach(w)
    # create window manager
    wm = WindowManager()
    wm.add(wi)
    wm.run()

def load_or_create(subject, side, len_threshold=20, k=100, outdir='data_als/cache/', seed=0):
    filename = 'data_als/%d/tracks_dti_3M_linear.trk' % subject

    print "Loading", filename
    streamlines, header = read(filename)
    streamlines = np.array(streamlines, dtype=np.object)[:,0]

    # hd = md5(streamlines).hexdigest()
    # print "hexdigest:", hd

    filename_cst = 'data_als/%d/%d_corticospinal_%s_3M.pkl'
    filename_cst = filename_cst % (subject, subject_segmentation[subject], side)
    print "Loading CST", filename_cst
    cst_ids = np.load(filename_cst)
    # cst_streamlines = streamlines[cst_ids]

    print "Building the dissimilarity representation."
    try:
        filename_prototypes = outdir+'Pi_ids_%d_%s.npy' % (subject, side)
        print "Trying to load", filename_prototypes
        Pi_ids = np.load(filename_prototypes)
        print "Done."
    except IOError:
        print "Not found."
        print "Creating prototypes."
        lenghts = np.array([len(s) for s in streamlines])
        streamlines_long_ids = np.where(lenghts > len_threshold)[0] # using long streamlines heuristic
        distance = bundles_distances_mam
        np.random.seed(seed)
        Pi_ids = streamlines_long_ids[subset_furthest_first(streamlines[streamlines_long_ids], k=k, distance=distance)] # using long streamlines heuristic
        print "Saving", filename_prototypes
        np.save(filename_prototypes, Pi_ids)
        Pi = streamlines[Pi_ids]
        
    try:
        filename_dr = outdir+'dr_%d_%s.npy' % (subject, side)
        print "Trying to load", filename_dr
        dr = np.load(filename_dr)
        print "Done."
    except IOError:
        print "Not found."
        print "Computing the dissimilarity matrix."
        dr = bundles_distances_mam(streamlines, Pi).astype(np.float32)
        print "Saving", filename_dr
        np.save(filename_dr, dr.astype(np.float32))

    return streamlines, cst_ids, Pi_ids, dr

def load_PX_tracks():

    roi = "LH_premotor"

    dn = "/home/hadron/from_John_mon12thmarch"
    dname = "/extra_probtrackX_analyses/_subject_id_subj05_101_32/particle2trackvis_" + roi + "_native/"
    fname = dn + dname + "tract_samples.trk"
    from nibabel import trackvis as tv

    points_space = [None, "voxel", "rasmm"]
    streamlines, hdr = tv.read(fname, as_generator=True, points_space="voxel")
    tracks = [s[0] for s in streamlines]
    del streamlines
    # return tracks

    qb = QuickBundles(tracks, 25.0 / 2.5, 18)
    # tl=Line(qb.exemplars()[0],line_width=1)
    del tracks
    qb.remove_small_clusters(20)

    tl = TrackLabeler(qb, qb.downsampled_tracks(), vol_shape=None, tracks_line_width=3.0, tracks_alpha=1)

    # put the seeds together
    # seeds=np.vstack((seeds,seeds2))
    # shif the seeds
    # seeds=np.dot(mat[:3,:3],seeds.T).T + mat[:3,3]
    # seeds=seeds-shift
    # seeds2=np.dot(mat[:3,:3],seeds2.T).T + mat[:3,3]
    # seeds2=seeds2-shift
    # msk = Point(seeds,colors=(1,0,0,1.),pointsize=2.)
    # msk2 = Point(seeds2,colors=(1,0,.ppppp2,1.),pointsize=2.)
    w = World()
    w.add(tl)
    # w.add(msk)
    # w.add(msk2)
    # w.add(sl)
    # create window
    wi = Window(caption="Fos", bgcolor=(0.3, 0.3, 0.6, 1.0), width=1600, height=900)
    wi.attach(w)
    # create window manager
    wm = WindowManager()
    wm.add(wi)
    wm.run()

def getPointsFromTrack(filename):
	geo = hou.pwd().geometry()
	# Read in stream data 
	streams, hdr = trackvis.read(filename)

	streamlines = [s[0] for s in streams]

	# For each streamline add a curve to the geometry
	j = 0
	for stream in streamlines:		
		i = 0
		curve = geo.createNURBSCurve(len(stream))
		if hou.updateProgressAndCheckForInterrupt(int(float(j)/float(len(streamlines))*100)):
			break
		for vertex in curve.vertices():
		    vertex.point().setPosition((float(stream[i][0]),float(stream[i][1]),float(stream[i][2])))
		    i = i + 1
		    if hou.updateProgressAndCheckForInterrupt():
        		break
        j = j+1

def filterlength(dname, fdwi, ffa, ftrk, thr_length, show=False):
    
    fa_img = nib.load(ffa)
    fa = fa_img.get_data()
    affine = fa_img.get_affine()

    img = nib.load(fdwi)
    data = img.get_data()
    
    from nibabel import trackvis
    streams, hdr = trackvis.read(ftrk)
    streamlines = [s[0] for s in streams]
    
    # threshold on streamline length

    from dipy.tracking.utils import length
    lengths = list(length(streamlines))

    new_streamlines = [ s for s, l in zip(streamlines, lengths) if l > thr_length ] #3.5
    
    # info length streamlines

    print(len(streamlines))
    print(len(new_streamlines))

    print(max(length(streamlines)))
    print(min(length(streamlines)))

    print(max(length(new_streamlines)))
    print(min(length(new_streamlines)))
    
    # show new tracto

    new_streamlines = list(new_streamlines)
    new_lengths = list(length(new_streamlines))

    fnew_tractogram = dname + 'filteredtractogram.trk'
    save_trk_old_style(fnew_tractogram, new_streamlines, affine, fa.shape)

    if show:
        show_results(data, new_streamlines, fa, affine, opacity=0.6)

def compute_length_array(trkfile=None, streams=None, savefname = 'lengths.npy'):
    if streams is None and not trkfile is None:
        print("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
            print(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] = length(fib[0])
    
    # store length array
    np.save(savefname, leng)
    print("Store lengths array to: %s" % savefname)
    
    return leng

 def _run_interface(self, runtime):
     
     tracts, hdr = tv.read(self.inputs.trackvis_file, as_generator=True)
     
     self.stat_files_data = []
     self.max_maps_data = []
     self.mean_maps_data = []
     for stat_file in self.inputs.stat_files:
         fmri_nii = nb.load(stat_file)
         self.stat_files_data.append(fmri_nii.get_data())
     
         self.max_maps_data.append(np.zeros(fmri_nii.get_header().get_data_shape()))
         self.mean_maps_data.append(np.zeros(fmri_nii.get_header().get_data_shape()))
     
     hdr = hdr.copy()
     if isdefined(self.inputs.stat_labels) and len(self.inputs.stat_labels) == len(self.inputs.stat_files):
         for i, label in enumerate(self.inputs.stat_labels):
             hdr['property_name'][i] = ('max_%s'%label)[0:19]
             #hdr['property_name'][1+i*2] = 'stat_mean_%s'%label
     else:        
         for i in range(len(self.inputs.stat_files)):
             hdr['property_name'][i] = ('max%d'%i)[0:19]
             #hdr['property_name'][1+i*2] = 'stat_mean%d'%i
     
     tv.write(self.inputs.out_tracks, self._gen_annotate_tracts(tracts, hdr), hdr)
     
     if isdefined(self.inputs.stat_labels) and len(self.inputs.stat_labels) == len(self.inputs.stat_files):
         for i, label in enumerate(self.inputs.stat_labels):
             nb.save(nb.Nifti1Image(self.max_maps_data[i], fmri_nii.get_affine(), fmri_nii.get_header()), self.inputs.out_max_map_prefix + "_%s"%label + '.nii')
             nb.save(nb.Nifti1Image(self.mean_maps_data[i], fmri_nii.get_affine(), fmri_nii.get_header()), self.inputs.out_mean_map_prefix + "_%s"%label + '.nii')
     else:
         for i in range(len(self.inputs.stat_files)):
             nb.save(nb.Nifti1Image(self.max_maps_data[i], fmri_nii.get_affine(), fmri_nii.get_header()), self.inputs.out_max_map_prefix + str(i) + '.nii')
             nb.save(nb.Nifti1Image(self.mean_maps_data[i], fmri_nii.get_affine(), fmri_nii.get_header()), self.inputs.out_mean_map_prefix + str(i) + '.nii')
     
     del self.mean_maps_data
     del self.max_maps_data
     del self.stat_files_data
         
     return runtime

    def _run_interface(self, runtime):
        from numpy import min_scalar_type
        tracks, header = nbt.read(self.inputs.in_file)
        streams = ((ii[0]) for ii in tracks)

        if isdefined(self.inputs.reference):
            refnii = nb.load(self.inputs.reference)
            affine = refnii.get_affine()
            data_dims = refnii.get_shape()[:3]
            kwargs = dict(affine=affine)
        else:
            iflogger.warn(('voxel_dims and data_dims are deprecated'
                           'as of dipy 0.7.1. Please use reference '
                           'input instead'))

            if not isdefined(self.inputs.data_dims):
                data_dims = header['dim']
            else:
                data_dims = self.inputs.data_dims
            if not isdefined(self.inputs.voxel_dims):
                voxel_size = header['voxel_size']
            else:
                voxel_size = self.inputs.voxel_dims

            affine = header['vox_to_ras']
            kwargs = dict(voxel_size=voxel_size)

        data = density_map(streams, data_dims, **kwargs)
        data = data.astype(min_scalar_type(data.max()))
        img = nb.Nifti1Image(data, affine)
        out_file = op.abspath(self.inputs.out_filename)
        nb.save(img, out_file)

        iflogger.info(
            ('Track density map saved as {i}, size={d}, '
             'dimensions={v}').format(
                i=out_file,
                d=img.get_shape(),
                v=img.get_header().get_zooms()))
        return runtime