def feature_map(self, feature, imPath=None, pw=0.95):
     """
     Given a set of feature values, produce a feature map,
     assuming that one feature corresponds to one region
     
     Parameters
     ----------
     feature, array of shape (self.k) : the information to map
     imPath=None, string yielding the output image path
                  if not None
     pw=0.95: volume of the Gaussian ellipsoid associated with the ROIs
     
     Returns
     -------
     The image object
     """
     if np.size(feature)!=self.k:
         raise ValueError, 'Incompatible feature dimension'
     from nipy.io.imageformats import save, Nifti1Image          
             
     label = self.map_label(self.generate_coordinates(), pval=pw)
     label = np.reshape(label, self.shape)
     values = np.zeros(self.shape)
     values[label>-1] = feature[label[label>-1].astype(np.int)]
     wim = Nifti1Image(values, self.affine)
     wim.get_header()['descrip']='feature image'
     if imPath!=None:
        save(wim,imPath)
     
     return wim

    def make_image(self, path=None):
        """
        write a int image where the nonzero values are the ROIs

        Parameters
        ----------
        path: string, optional
            the desired image path

        Returns
        -------
        brifti image instance
        
        Note
        ----
        the background values are set to -1
        the ROIs values are set as [0..self.k-1]
        """
        if self.shape==None:
            raise ValueError, 'Need self.shape to be defined'

        data = -np.ones(self.shape,np.int)
        for k in range(self.k):
            dk = self.xyz[k].T
            data[dk[0], dk[1], dk[2]] = k

        wim =  Nifti1Image(data, self.affine)
        header = wim.get_header()
        header['descrip'] = "Multiple ROI image"
        if path!=None:
            save(wim, path)
        return wim

def ffx( maskImages, effectImages, varianceImages, resultImage=None):
    """
    Computation of the fixed effecst statistics

    Parameters
    ----------
    maskImages, string or list of strings
                the paths of one or several masks
                when several masks, the half thresholding heuristic is used
    effectImages, list of strings
                the paths ofthe effect images   
    varianceImages, list of strings
                    the paths of the associated variance images
    resultImage=None, string,
                 path of the result images

    Returns
    -------
    the computed values
    """
    # fixme : check that the images have same referntial
    # fixme : check that mask_Images is a list
    if len(effectImages)!=len(varianceImages):
        raise ValueError, 'Not the correct number of images'
    tiny = 1.e-15
    nsubj = len(effectImages)
    mask = intersect_masks(maskImages, None, threshold=0.5, cc=True)
    
    effects = []
    variance = []
    for s in range(nsubj):
        rbeta = load(effectImages[s])
        beta = rbeta.get_data()[mask>0]
        rbeta = load(varianceImages[s])
        varbeta = rbeta.get_data()[mask>0]
        
        effects.append(beta)
        variance.append(varbeta)
    
    effects = np.array(effects)
    variance = np.array(variance)
    effects[np.isnan(effects)] = 0
    effects[np.isnan(variance)] = 0
    variance[np.isnan(variance)] = tiny
    variance[variance==0] = tiny    
    
    t = effects/np.sqrt(variance)
    t = t.mean(0)*np.sqrt(nsubj)     
    #t = np.sum(effects/variance,0)/np.sum(1.0/np.sqrt(variance),0)

    nim = load(effectImages[0])
    affine = nim.get_affine()
    tmap = np.zeros(nim.get_shape())
    tmap[mask>0] = t
    tImage = Nifti1Image(tmap, affine)
    if resultImage!=None:
       save(tImage, resultImage)

    return tmap

def save_results(ppm, id_algo, noise): 
    savedir = os.path.join(baseres, subjects[s_idx])
    tag = id_algo
    if noise == 'laplace': 
        tag += '_laplace'
    for i in range(ntissues): 
        im = Image(ppm[:,:,:,i], affine)
        fname = id_posterior + '_' + regs[r_idx] + '_' + tag + '_' + str(i) + '.nii'
        save(Nifti1Image(im), os.path.join(savedir, fname))

def intersect_masks(input_masks, output_filename=None, threshold=0.5, cc=True):
    """
    Given a list of input mask images, generate the output image which
    is the the threshold-level intersection of the inputs 

    
    Parameters
    ----------
    input_masks: list of strings or ndarrays
        paths of the input images nsubj set as len(input_mask_files), or
        individual masks.
    output_filename, string:
        Path of the output image, if None no file is saved.
    threshold: float within [0, 1], optional
        gives the level of the intersection.
        threshold=1 corresponds to keeping the intersection of all
        masks, whereas threshold=0 is the union of all masks.
    cc: bool, optional
        If true, extract the main connected component
        
    Returns
    -------
    grp_mask, boolean array of shape the image shape
    """  
    grp_mask = None 

    for this_mask in input_masks:
        if isinstance(this_mask, basestring):
            # We have a filename
            this_mask = load(this_mask).get_data()
        if grp_mask is None:
            grp_mask = this_mask.copy().astype(np.int)
        else:
            grp_mask += this_mask
    
    grp_mask = grp_mask>(threshold*len(input_masks))
    if np.any(grp_mask>0) and cc:
        grp_mask = largest_cc(grp_mask)
    
    if output_filename is not None:
        if isinstance(input_masks[0], basestring):
            nim = load(input_masks[0]) 
            header = nim.get_header()
            affine = nim.get_affine()
        else:
            header = dict()
            affine = np.eye(4)
        header['descrip'] = 'mask image'
        output_image = nifti1.Nifti1Image(grp_mask.astype(np.uint8),
                                            affine=affine,
                                            header=header,
                                         )
        save(output_image, output_filename)

    return grp_mask>0

 def to_image(self, path=None):
     """
     Write itself as an image, and returns it
     """
     data = np.zeros(self.shape).astype(np.int8)
     data[self.ijk[:,0], self.ijk[:,1], self.ijk[:,2]] = 1
     nim = Nifti1Image(data, self.affine)
     nim.get_header()['descrip'] = 'mask image'
     if path is not None:
         save(nim, path)
     return nim

def save(filename, obj):
    """ Save an nipy image object to a file.
    """
    obj = as_volume_img(obj, copy=False)
    hdr = imageformats.Nifti1Header()
    for key, value in obj.metadata.iteritems():
        if key in hdr:
            hdr[key] = value
    img = imageformats.Nifti1Image(obj.get_data(), 
                                   obj.affine,
                                   header=hdr)
    imageformats.save(img, filename)

def parcellation_output_with_paths(Pa, mask_images, group_path, indiv_path):
    """
    Function that produces images that describe the spatial structure
    of the parcellation.  It mainly produces label images at the group
    and subject level
    
    Parameters
    ----------
    Pa : Parcellation instance that describes the parcellation
    mask_images: list of images paths that define the mask
    coord: array of shape (nvox,3) that contains(approximated)
           MNI-coordinates of the brain mask voxels considered in the
           parcellation process
    group_path, string, path of the group-level parcellation image
    indiv_path, list of strings, paths of the individual parcellation images    
    
    fixme
    -----
    the referential-defining information should be part of the Pa instance
    """
    nsubj = Pa.nb_subj
    mxyz = Pa.ijk
    
    # write the template image
    tlabs = Pa.group_labels
    rmask = load(mask_images[0])
    ref_dim = rmask.get_shape()
    grid_size = np.prod(ref_dim)
    affine = rmask.get_affine()
    
    Label = np.zeros(ref_dim)
    Label[Pa.ijk[:,0],Pa.ijk[:,1],Pa.ijk[:,2]]=tlabs+1
    
    wim = Nifti1Image (Label, affine)
    hdr = wim.get_header()
    hdr['descrip'] = 'group_level Label image obtained from a \
                     parcellation procedure'
    save(wim, group_path)
    
    # write subject-related stuff
    for s in range(nsubj):
        # write the images
        labs = Pa.label[:,s]
        Label = np.zeros(ref_dim).astype(np.int)
        Label[Pa.ijk[:,0],Pa.ijk[:,1],Pa.ijk[:,2]]=labs+1
        wim = Nifti1Image (Label, affine)
        hdr = wim.get_header()
        hdr['descrip'] = 'individual Label image obtained \
                         from a parcellation procedure'
        save(wim, indiv_path[s])

def ffx_from_stat( maskImages, statImages, resultImage=None):
    """
    Computation of the fixed effects statistics from statistic
    

    Parameters
    ----------
    maskImages, string or list of strings
                the paths of one or several masks
                when several masks, the half thresholding heuristic is used
    statImages, list of strings
                the paths ofthe statitsic images   
    resultImage=None, string,
                 path of the result images

    Returns
    -------
    the computed values
    """
    # fixme : check that the images have same referntial
    # fixme : check that mask_Images is a list
    nsubj = len(statImages)
    mask = intersect_masks(maskImages, None, threshold=0.5, cc=True)
    
    t = []
    for s in range(nsubj):
        rbeta = load(statImages[s])
        beta = rbeta.get_data()[mask>0]            
        t.append(beta)
    
    t = np.array(t)
    t[np.isnan(t)] = 0
    t = t.mean(0)*np.sqrt(nsubj)     

    nim = load(statImages[0])
    affine = nim.get_affine()
    tmap = np.zeros(nim.get_shape())
    tmap[mask>0] = t
    tImage = Nifti1Image(tmap, affine)
    if resultImage!=None:
       save(tImage,resultImage)

    return tmap

def mask_parcellation(mask_images, nb_parcel, output_image=None):
    """
    Performs the parcellation of a certain mask

    Parameters
    ----------
    mask_images: list of strings,
                 paths of the mask images that define the common space.
    nb_parcel: int,
               number of desired parcels
    output_image: string, optional
                   path of the output image
                   
    Returns
    -------
    wim: Nifti1Imagine instance,  the resulting parcellation
    """
    from ..mask import intersect_masks

    # compute the group mask
    affine = load(mask_images[0]).get_affine()
    shape = load(mask_images[0]).get_shape()
    mask = intersect_masks(mask_images, threshold=0)>0
    ijk = np.where(mask)
    ijk = np.array(ijk).T
    nvox = ijk.shape[0]

    # Get and cluster  coordinates 
    ijk = np.hstack((ijk,np.ones((nvox,1))))
    coord = np.dot(ijk, affine.T)[:,:3]
    cent, tlabs, J = kmeans(coord, nb_parcel)
        
    # Write the results
    label = -np.ones(shape)
    label[mask]= tlabs
    wim = Nifti1Image(label, affine)
    wim.get_header()['descrip'] = 'Label image in %d parcels'%nb_parcel    
    if output_image is not None:
        save(wim, output_image)
    return wim

def save_volume(shape, path, affine, mask=None, data=None, descrip=None):
    """
    volume saving utility for masked volumes
    
    Parameters
    ----------
    shape, tupe of dimensions of the data
    path, string, output image path
    affine, transformation of the grid to a coordinate system
    mask=None, binary mask used to reduce the volume size
    data=None data to be put in the volume
    descrip=None, a string descibing what the image is

    Fixme
    -----
    Handle the case where data is multi-dimensional
    """
    volume = np.zeros(shape)
    if mask== None: 
       print "Could not write the image: no mask"
       return

    if data == None:
       print "Could not write the image: no data"
       return

    if np.size(data.shape) == 1:
        volume[mask > 0] = data
    else:
        for i in range(data.shape[0]):
            volume[i][mask[0] > 0] = data[i]

    wim = Nifti1Image(volume, affine)
    if descrip !=None:
        wim.get_header()['descrip']=descrip
    save(wim, path)

def Parcellation_output(Pa, mask_images, learning_images, coord, nbru, 
                        verbose=1,swd = "/tmp"):
    """
    Function that produces images that describe the spatial structure
    of the parcellation.  It mainly produces label images at the group
    and subject level
    
    Parameters
    ----------
    Pa : Parcellation instance that describes the parcellation
    mask_images: list of images paths that define the mask
    learning_images: list of float images containing the input data
    coord: array of shape (nvox,3) that contains(approximated)
           MNI-coordinates of the brain mask voxels considered in the
           parcellation process
    nbru: list of subject ids
    verbose=1 : verbosity level
    swd = '/tmp': write directory
    
    Results
    -------
    Pa: the updated Parcellation instance
    """
    nsubj = Pa.nb_subj
    Pa.set_subjects(nbru)
    
    # write the template image
    tlabs = Pa.group_labels
    LabelImage = os.path.join(swd,"template_parcel.nii") 
    rmask = load(mask_images[0])
    ref_dim = rmask.get_shape()
    affine = rmask.get_affine()
    
    Label = np.zeros(ref_dim)
    Label[Pa.ijk[:,0],Pa.ijk[:,1],Pa.ijk[:,2]]=tlabs+1
    
    wim = Nifti1Image (Label, affine)
    hdr = wim.get_header()
    hdr['descrip'] = 'group_level Label image obtained from a \
                     parcellation procedure'
    save(wim, LabelImage)
    
    # write subject-related stuff
    Jac = []
    if Pa.isfield('jacobian'):
        Jac = Pa.get_feature('jacobian')
        Jac = np.reshape(Jac,(Pa.k,nsubj))
        
    for s in range(nsubj):
        # write the images
        labs = Pa.label[:,s]
        LabelImage = os.path.join(swd,"parcel%s.nii" % nbru[s])
        JacobImage = os.path.join(swd,"jacob%s.nii" % nbru[s])      

        Label = np.zeros(ref_dim).astype(np.int)
        Label[Pa.ijk[:,0],Pa.ijk[:,1],Pa.ijk[:,2]]=labs+1
        wim = Nifti1Image (Label, affine)
        hdr = wim.get_header()
        hdr['descrip'] = 'individual Label image obtained \
                         from a parcellation procedure'
        save(wim, LabelImage)

        if ((verbose)&(np.size(Jac)>0)):
            Label = np.zeros(ref_dim)
            Label[Pa.ijk[:,0],Pa.ijk[:,1],Pa.ijk[:,2]]=Jac[labs,s]
            wim = Nifti1Image (Label, affine)
            hdr = wim.get_header()
            hdr['descrip'] = 'image of the jacobian of the deformation \
                              associated with the parcellation'
            save(wim, JacobImage)       

    return Pa

grp_mask = Nifti1Image(mask, load(mask_images[0]).get_affine())
ijk = np.array(np.where(mask)).T
nvox = np.sum(mask)

# output dir
b_smooth = True
if b_smooth:
    print "smoothed data"
    threshold_path = 'volume_threshold_smooth.con'
    swd = '/data/home/virgile/virgile_internship/group_analysis/smoothed_FWHM5'
else:
    print "unsmoothed data"
    threshold_path = 'volume_threshold.con'
    swd = '/data/home/virgile/virgile_internship/group_analysis/smoothed_FWHM0'

save(grp_mask, op.join(swd,'grp_mask.nii'))

################################################################
# Load the effects and variance
################################################################

def load_images(con_images, var_images):
    """
    """
    nsubj = len(con_images)
    beta = []
    varbeta = []
    tiny = 1.e-15
    for s in range(nsubj): 
        rbeta = load(con_images[s])
        temp = (rbeta.get_data())[mask]

########################################

paradigm = dm.EventRelatedParadigm(conditions, onsets)
X, names = dm.dmtx_light(frametimes, drift_model='Cosine', hfcut=128,
               hrf_model=hrf_model, paradigm=paradigm, add_regs=motion, add_reg_names=add_reg_names)


#######################################
# Get the FMRI data
#######################################

fmri_data = surrogate_4d_dataset(shape=shape, n_scans=n_scans)[0]

# if you want to save it as an image
data_file = op.join(swd,'fmri_data.nii')
save(fmri_data, data_file)

########################################
# Perform a GLM analysis
########################################

# GLM fit
Y = fmri_data.get_data()
model = "ar1"
method = "kalman"
glm = GLM.glm()
mp.pcolor(X)
mp.show()
glm.fit(Y.T, X, method=method, model=model)
#explained = np.dot(X,glm.beta.reshape(X.shape[1],-1)).reshape(Y.T.shape).T
#residuals = Y - explained 

label = -np.ones(F.V)
nroi = hroi.NROI_from_field(F, affine, shape, xyz, 0, threshold, smin)
bmap = -np.zeros(F.V)
if nroi!=None:
    idx = nroi.discrete_features['index']
    for k in range(nroi.k):
        label[idx[k]] = k

# saving the blob image,i. e. a label image 
wlabel = -2*np.ones(shape)
wlabel[data!=0] = label
blobPath = os.path.join(swd,"blob.nii")

wim = Nifti1Image(wlabel, affine)
wim.get_header()['descrip'] = 'blob image extracted from %s'%InputImage
save(wim, blobPath)

# --- 2.b take blob labelled "1" as an ROI
roi = DiscreteROI( affine=affine, shape=shape)
roi.from_labelled_image(blobPath, 1)
roiPath2 = os.path.join(swd, "roi_blob_1.nii")
roi.make_image(roiPath2)

# --- 2.c take the blob closest to 'position as an ROI'
roiPath3 = os.path.join(swd, "blob_closest_to_%d_%d_%d.nii")%\
           (position[0], position[1], position[2])
roi.from_position_and_image(blobPath, np.array(position))
roi.make_image(roiPath3)

# --- 2.d make a set of ROIs from all the blobs
mroi = MultipleROI( affine=affine, shape=shape)

#########################################
# Estimate the contrasts
#########################################

print 'Computing contrasts...'
for index, contrast_id in enumerate(contrasts):
    print '  Contrast % 2i out of %i: %s' % (index+1, 
                                             len(contrasts), contrast_id)
    lcontrast = my_glm.contrast(contrasts[contrast_id])
    # 
    contrast_path = op.join(swd, '%s_z_map.nii'% contrast_id)
    write_array = mask_array.astype(np.float)
    write_array[mask_array] = lcontrast.zscore()
    contrast_image = Nifti1Image(write_array, fmri_image.get_affine() )
    save(contrast_image, contrast_path)
    affine = fmri_image.get_affine()

    
    vmax = max(-write_array.min(), write_array.max())
    plot_map(write_array, affine, 
             cmap=cm.cold_hot, 
             vmin=-vmax,
             vmax=vmax,
             anat=None,
             figure=10,
             threshold=2.5)
    pylab.savefig(op.join(swd, '%s_z_map.png' % contrast_id))
    pylab.clf()
    

R.set_source_fov(fixed_npoints=64**3)

# Make Gaussian spline transform instance
spacing = 16
slices = [slice(0,s.stop,s.step*spacing) for s in R._slices]
cp = np.mgrid[slices]
cp = np.rollaxis(cp, 0, 4)

# Start with an affine registration
A0 = Affine()
A = R.optimize(A0)
###A = Affine()

# Save affinely transformed target  
Jt = J.transform(A, reference=I)
save(asNifti1Image(Jt), 'affine_anubis_to_ammon.nii')

# Then add control points...
T0 = SplineTransform(I, cp, sigma=20., grid_coords=True, affine=A)

"""
# Test 1
s = R.eval(T0)
sa = R.eval(T0.affine)
assert_almost_equal(s, sa)

# Test 2
T = SplineTransform(I, cp, sigma=5., grid_coords=True, affine=A)
T.param += 1.
s0 = R.eval(T0)
s = R.eval(T)

def Parcellation_based_analysis(Pa, test_images, numbeta, swd="/tmp", 
                                    DMtx=None, verbose=1, method_id=0):
    """
    This function computes parcel averages and RFX at the parcel-level

    Parameters
    ----------
    Pa Parcellation instance that is updated in this function
    test_images: double list of paths of functional images used 
                 as input to for inference. 
                 Normally these are contrast images.
                 double list is 
                 [number of subjects [number of contrasts]]
    numbeta: list of int of the associated ids
    swd='/tmp': write directory
    DMtx=None: array od shape (nsubj,ncon) 
               a design matrix for second-level analyses 
              (not implemented yet)
    verbose=1: verbosity level
    method_id = 0: an id of the method used.
              This is useful to compare the outcome of different 
              Parcellation+RFX  procedures

    Results
    -------
    Pa: the updated Parcellation instance
    """
    nsubj = Pa.nb_subj
    mxyz = Pa.ijk.T
    mask = Pa.label>-1
    nbeta = len(numbeta)
    
    # 1. read the test data
    # fixme: Check that everybody is in the same referential
    Test = []
    for s in range(nsubj):
        beta = []
        lxyz = mxyz[:,mask[:,s]]
        lxyz = np.array(lxyz)

        for b in range(nbeta):
            # the raw contrast images   
            rbeta = load(test_images[s][b])
            temp = rbeta.get_data()
            temp = temp[lxyz[0,:],lxyz[1,:],lxyz[2,:]]
            temp = np.reshape(temp, np.size(temp))
            beta.append(temp)
            temp[np.isnan(temp)]=0 ##

        beta = np.array(beta)
        Test.append(beta.T) 

    # 2. compute the parcel-based stuff
    # and make inference inference (RFX,...)

    prfx = np.zeros((Pa.k,nbeta))
    vinter = np.zeros(nbeta)
    for b in range(nbeta):
        unitest = [np.reshape(Test[s][:,b],(np.size(Test[s][:,b]),1)) \
                  for s in range(nsubj)]
        cname = 'contrast_%04d'%(numbeta[b])
        Pa.make_feature(unitest, cname)
        prfx[:,b] =  np.reshape(Pa.PRFX(cname,1),Pa.k)
        vinter[b] = Pa.variance_inter(cname)

    vintra = Pa.variance_intra(Test)

    if verbose:
        print 'average intra-parcel variance', vintra
        print 'average intersubject variance', vinter.mean()
            
    # 3. Write the stuff
    # write RFX images
    ref_dim = rbeta.get_shape()
    affine = rbeta.get_affine()
    tlabs = Pa.group_labels

    # write the prfx images
    for b in range(len(numbeta)):
        RfxImage = os.path.join(swd,"prfx_%s_%d.nii" % (numbeta[b],method_id))
        if ((verbose)&(np.size(prfx)>0)):
            rfx_map = np.zeros(ref_dim)
            rfx_map[Pa.ijk[:,0],Pa.ijk[:,1],Pa.ijk[:,2]] = prfx[tlabs,b]
            wim = Nifti1Image (rfx_map, affine)
            hdr = wim.get_header()
            hdr['descrip'] = 'parcel-based eandom effects image (in z-variate)'
            save(wim, RfxImage)     
        
    return Pa

#.........这里部分代码省略.........
    if method not in ['ward','gkm','ward_and_gkm','kmeans']:
        raise ValueError, 'unknown method'
    if nn not in [6,18,26]:
        raise ValueError, 'nn should be 6,18 or 26'
    nbeta = len(betas)
    
    # step 1: load the data ----------------------------
    #1.1 the mask image
    nim = load(MaskImage)
    ref_dim =  nim.get_shape()
    affine = nim.get_affine()
    mask = nim.get_data()
    xyz = np.array(np.where(mask>0)).T
    nvox = xyz.shape[0]

    if method is not 'kmeans':
        # 1.2 get the main cc of the graph 
        # to remove the small connected components
        g = fg.WeightedGraph(nvox)
        g.from_3d_grid(xyz.astype(np.int),nn)
        
        aux = np.zeros(g.V).astype('bool')
        imc = g.main_cc()
        aux[imc]= True
        if np.sum(aux)==0:
            raise ValueError, "empty mask. Cannot proceed"
        g = g.subgraph(aux)
        lmask = np.zeros(ref_dim)
        lmask[xyz[:,0],xyz[:,1],xyz[:,2]]=aux
        xyz = xyz[aux,:]
        nvox = xyz.shape[0]

    # 1.3 from vox to mm
    xyz2 = np.hstack((xyz,np.ones((nvox,1))))
    coord = np.dot(xyz2, affine.T)[:,:3]

    # 1.4 read the functional data
    beta = []
    for b in range(nbeta):
        rbeta = load(betas[b])
        lbeta = rbeta.get_data()
        lbeta = lbeta[lmask>0]
        beta.append(lbeta)
	
    beta = np.array(beta).T

    #step 2: parcel the data ---------------------------
    feature = np.hstack((beta, mu*coord/np.std(coord)))
    if method is not 'kmeans':
        g = ff.Field(nvox, g.edges, g.weights, feature)

    if method=='kmeans':
        cent, u, J = kmeans(feature, nbparcel)

    if method=='ward':
        u, J0 = g.ward(nbparcel)

    if method=='gkm':
        seeds = np.argsort(np.random.rand(g.V))[:nbparcel]
        seeds, u, J1 = g.geodesic_kmeans(seeds)

    if method=='ward_and_gkm':
        w,J0 = g.ward(nbparcel)
        seeds, u, J1 = g.geodesic_kmeans(label=w)

    lpa = Parcellation(nbparcel, xyz, np.reshape(u,(nvox,1)))
    if verbose:
        pi = np.reshape(lpa.population(), nbparcel)
        vi = np.sum(lpa.var_feature_intra([beta])[0], 1)
        vf = np.dot(pi,vi)/nvox
        va =  np.dot(pi,np.sum(lpa.var_feature_intra([coord])[0],1))/nvox
        print nbparcel, "functional variance", vf, "anatomical variance",va


    # step3:  write the resulting label image
    Label = -np.ones(ref_dim,'int16')
    Label[lmask>0] = u

    if fullpath is not None:
        LabelImage = fullpath
    elif write_dir is not None:
        if method=='kmeans':
            LabelImage = os.path.join(write_dir,"parcel_kmeans.nii")
        if method=='ward':
            LabelImage = os.path.join(write_dir,"parcel_wards.nii")
        elif method=='gkm':
            LabelImage = os.path.join(write_dir,"parcel_gkmeans.nii")
        elif method=='ward_and_gkm':
            LabelImage = os.path.join(write_dir,"parcel_wgkmeans.nii")
    else:
        LabelImage = None
    
    if LabelImage is not None:
        wim = Nifti1Image(Label, affine)
        hdr = wim.get_header()
        hdr['descrip'] = 'Intra-subject parcellation image'
        save(wim, LabelImage)
        print "Wrote the parcellation images as %s" %LabelImage

    return lpa, Label

"""
print __doc__

import numpy as np
import os
from nipy.io.imageformats import load, save, Nifti1Image 
from nipy.neurospin.graph.field import Field
import get_data_light
import tempfile
data_dir = get_data_light.get_it()

# paths
swd = tempfile.mkdtemp()
input_image = os.path.join(data_dir, 'spmT_0029.nii.gz')
mask_image = os.path.join(data_dir, 'mask.nii.gz')

mask = load(mask_image).get_data()>0
ijk = np.array(np.where(mask)).T
nvox = ijk.shape[0]
data = load(input_image).get_data()[mask]
image_field = Field(nvox)
image_field.from_3d_grid(ijk, k=6)
image_field.set_field(data)
u = image_field.ward(100)

label_image = os.path.join(swd, 'label.nii')
wdata = mask - 1
wdata[mask] = u
save(Nifti1Image(wdata, load(mask_image).get_affine()), label_image)
print "Label image written in %s"  % label_image