def load_ppm(): 
    # Read input files 
    for i in range(ntissues):
        fname = os.path.join(subdatadir, id_prior+'_EM_'+str(i)+'.nii')
        print(fname)
        if i == 0: 
            im = Image(load(fname))
            affine = im.affine
            data = np.zeros(list(im.shape)+[ntissues])
            data[:,:,:,0] = im.data
        else: 
            data[:,:,:,i] = Image(load(fname)).data
    
    #ppm = Image(data, affine)

    # Normalize and mask ppms 
    psum = data.sum(3)
    X,Y,Z = np.where(psum>0)
    for i in range(ntissues): 
        data[X,Y,Z,i] /= psum[X,Y,Z]
    mask = (X.astype('uint'), Y.astype('uint'), Z.astype('uint'))

    # For now, we may need to reduce the mask because the vem module
    # does not handle boundary conditions
    X,Y,Z = mask
    I = np.where((X>0)*(X<im.shape[0]-1)*(Y>0)*(Y<im.shape[1]-1)*(Z>0)*(Z<im.shape[2]-1))
    mask = X[I], Y[I], Z[I]
    
    return data, mask, affine 

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 series_from_mask(filenames, mask, dtype=np.float32, smooth=False):
    """ Read the time series from the given sessions filenames, using the mask.

        Parameters
        -----------
        filenames: list of 3D nifti file names, or 4D nifti filename.
            Files are grouped by session.
        mask: 3d ndarray
            3D mask array: true where a voxel should be used.
        smooth: False or float, optional
            If smooth is not False, it gives the size, in voxel of the
            spatial smoothing to apply to the signal.
        
        Returns
        --------
        session_series: ndarray
            3D array of time course: (session, voxel, time)
        header: header object
            The header of the first file.
    """
    assert len(filenames) != 0, ( 
        'filenames should be a file name or a list of file names, '
        '%s (type %s) was passed' % (filenames, type(filenames)))
    mask = mask.astype(np.bool)
    if isinstance(filenames, basestring):
        # We have a 4D nifti file
        data_file = load(filenames)
        header = data_file.get_header()
        series = data_file.get_data()
        affine = data_file.get_affine()[:3, :3]
        del data_file
        if isinstance(series, np.memmap):
            series = np.asarray(series).copy()
        if smooth:
            smooth_sigma = np.dot(linalg.inv(affine), np.ones(3))*smooth
            for this_volume in np.rollaxis(series, -1):
                this_volume[...] = ndimage.gaussian_filter(this_volume,
                                                        smooth_sigma)
        series = series[mask].astype(dtype)
    else:
        nb_time_points = len(list(filenames))
        series = np.zeros((mask.sum(), nb_time_points), dtype=dtype)
        for index, filename in enumerate(filenames):
            data_file = load(filename)
            data = data_file.get_data()
            if smooth:
                affine = data_file.get_affine()[:3, :3]
                smooth_sigma = np.dot(linalg.inv(affine), np.ones(3))*smooth
                data = ndimage.gaussian_filter(data, smooth_sigma)
                
            series[:, index] = data[mask].astype(dtype)
            # Free memory early
            del data
            if index == 0:
                header = data_file.get_header()

    return series, header

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 from_position_and_image(self, image_path, position):
        """
         Define the ROI as the set of  voxels of the image
         that is closest to the provided position

        Parameters
        -----------
        image_path: string, 
            the path of a label (discrete valued) image
        position: array of shape (3,)
            x, y, z position in the world space

        Notes
        -------
        everything could be performed in the image space
        """
        # check that the header is OK indeed
        self.check_header(image_path)

        # get the image data and find the best matching ROI
        nim = load(image_path)
        data = nim.get_data().astype(np.int)
        k = data.max()+1
        cent = np.array([np.mean(np.where(data==i),1) for i in range(k)])
        cent = np.hstack((cent,np.ones((k,1))))
        coord = np.dot(cent, self.affine.T)[:,:3]
        
        # find the best match
        dx = coord-position
        k = np.argmin(np.sum(dx**2,1))
        self.discrete = np.where(data==k)

def save_all_images(contrast, dim, mask_url, kargs):
    """
    idem savel_all, but the names are now all included in kargs
    """
    z_file = kargs["z_file"]
    t_file = kargs["t_file"]
    res_file = kargs["res_file"]
    con_file = kargs["con_file"]
    html_file = kargs["html_file"]
    mask = load(mask_url)
    mask_arr = mask.get_data()
    affine = mask.get_affine()
    shape = mask.get_shape()    
   
    # load the values
    t = contrast.stat()
    z = contrast.zscore()

    # saving the Z statistics map
    save_volume(shape, z_file, affine, mask_arr, z, "z_file")
    
    # Saving the t/F statistics map
    save_volume(shape, t_file, affine, mask_arr, t, "t_file")
    
    if int(dim) != 1:
        shape = (shape[0], shape[1], shape[2],int(dim)**2)
        contrast.variance = contrast.variance.reshape(int(dim)**2, -1)

    ## saving the associated variance map
    # fixme : breaks with F contrasts !
    if contrast.type == "t":
        save_volume(shape, res_file, affine, mask_arr,
                    contrast.variance)
    if int(dim) != 1:
        shape = (shape[0], shape[1], shape[2], int(dim))

    # writing the associated contrast structure
    # fixme : breaks with F contrasts !
    if contrast.type == "t":    
        save_volume(shape, con_file, affine, mask_arr,
                    contrast.effect)
 
    # writing the results as an html page
    if kargs.has_key("method"):
        method = kargs["method"]
    else:
        method = 'fpr'

    if kargs.has_key("threshold"):
        threshold = kargs["threshold"]
    else:
        threshold = 0.001

    if kargs.has_key("cluster"):
        cluster = kargs["cluster"]
    else:
        cluster = 0
    Results.ComputeResultsContents(z_file, mask_url, html_file,
                                   threshold=threshold, method=method,
                                   cluster=cluster)

    def set_discrete_feature_from_image(self, fid, image_path=None,
                                        image=None):
        """
        extract some discrete information from an image

        Parameters
        ----------
        fid: string, feature id
        image_path, string, optional
            input image path
        image, brfiti image path,
            input image

        Note that either image_path or image has to be provided
        """
        if image_path==None and image==None:
            raise ValueError, "one image needs to be provided"
        if image_path is not None:
            self.check_header(image_path)
            nim = load(image_path)
        if image is not None:
            nim = image
        data = nim.get_data()
        ldata = []
        for k in range(self.k):
            dk = self.xyz[k].T
            ldk = data[dk[0],dk[1],dk[2]]
            if np.size(ldk)==ldk.shape[0]:
                ldk = np.reshape(ldk,(np.size(ldk),1))
            ldata.append(ldk)
        self.set_discrete_feature(fid,ldata)

def load(filename):
    """Load an image from the given filename.

    Parameters
    ----------
    filename : string
        Should resolve to a complete filename path.

    Returns
    -------
    image : An `Image` object
        If successful, a new `Image` object is returned.

    See Also
    --------
    save_image : function for saving images
    fromarray : function for creating images from numpy arrays

    Examples
    --------

    >>> from nipy.io.api import load_image
    >>> from nipy.testing import anatfile
    >>> img = load_image(anatfile)
    >>> img.shape
    (33, 41, 25)
    """
    img = formats.load(filename)
    aff = img.get_affine()
    shape = img.get_shape()
    hdr = img.get_header()

    # Get info from NIFTI header, if present, to tell which axes are
    # which.  This is a NIFTI-specific kludge, that might be abstracted
    # out into the image backend in a general way.  Similarly for
    # getting zooms

    # axis_renames is a dictionary: dict([(int, str)])
    # that has keys in range(3)
    # the axes of the Image are renamed from 'ijk'
    # using these names

    try:
        axis_renames = hdr.get_axis_renames()
    except (TypeError, AttributeError):
        axis_renames = {}

    try:
        zooms = hdr.get_zooms()
    except AttributeError:
        zooms = np.ones(len(shape))

    # affine_transform is a 3-d transform

    affine_transform3d, affine_transform = \
        affine_transform_from_array(aff, 'ijk', pixdim=zooms[3:])
    img = Image(img.get_data(), affine_transform.renamed_domain(axis_renames))
    img.header = hdr
    return img

def test_conversion():

    brifti_obj = imageformats.load(data_file)
    vol_img = as_volume_img(data_file)
    yield nose.tools.assert_equals, as_volume_img(vol_img), \
                    vol_img
    yield nose.tools.assert_equals, as_volume_img(brifti_obj), \
                    vol_img

def load_image(image_path, mask_path=None ):
    """ Return an array of image data masked by mask data 

    Parameters
    ----------
    image_path string or list of strings 
               that yields the data of interest
    mask_path=None: string that yields the mask path

    Returns
    -------
    image_data a data array that can be 1, 2, 3  or 4D 
               depending on chether mask==None or not
               and on the length of the times series
    """
    # fixme : do some check
    if mask_path !=None:
       rmask = load(mask_path)
       shape = rmask.get_shape()[:3]
       mask = np.reshape(rmask.get_data(),shape)
    else:
        mask = None

    image_data = []
    
    if hasattr(image_path, '__iter__'):
       if len(image_path)==1:
          image_path = image_path[0]

    if hasattr(image_path, '__iter__'):
       for im in image_path:
           if mask is not None:
               temp = np.reshape(load(im).get_data(),shape)[mask>0,:]    
           else:
                temp = np.reshape(load(im).get_data(),shape) 
           image_data.append(temp)
       image_data = np.array(image_data).T
    else:
        image_data = load(image_path).get_data()
        if mask != None:
            image_data = image_data[mask>0,:]
    
    return image_data

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 read_ppms():
    """
    Open PPMs (White Matter, Gray Matter, CSF, Rest)
    """
    Pdict = {}
    for tissue in tissues: 
        fname = 'out'+tissue+'_100.img'
        im = brifti.load(os.path.join(datadir, fname))
        Pdict[tissue] = im.get_data()/1000.
    return Pdict

    def from_binary_image(self, image_path):
        """
        Take all the <>0 sites of the image as the ROI

        Parameters
        -----------
        image_path: string
            the path of an image
        """
        self.check_header(image_path)
        nim = load(image_path)
        self.discrete = np.where(nim.get_data())

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]
        beta.append(temp)
        rvar = load(var_images[s])
        temp = (rvar.get_data())[mask]
        varbeta.append(temp)

    VarFunctional = np.array(varbeta).T
    Functional = np.array(beta).T
    Functional[np.isnan(Functional)] = 0
    VarFunctional[np.isnan(VarFunctional)] = 0
    VarFunctional = np.maximum(VarFunctional, tiny)
    return Functional,  VarFunctional

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 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 set_feature_from_image(self, fid, image_path):
        """
        extract some roi-related information from an image

        Parameters
        -----------
        fid: string
            feature id
        image: string
            image path
        """
        self.check_header(image_path)
        nim = load(image_path)  
        data = nim.get_data()
        self.set_feature(fid,data)

 def get_anat(cls):
     filename = find_mni_template()
     if cls.anat is None:
         if filename is None:
             raise OSError('Cannot find template file T1_brain.nii.gz '
                     'required to plot anatomy, see the nipy documentation '
                     'installaton section for how to install template files.')
         anat_im = load(filename)
         anat = anat_im.get_data()
         anat = anat.astype(np.float)
         anat_mask = ndimage.morphology.binary_fill_holes(anat > 0)
         anat = np.ma.masked_array(anat, np.logical_not(anat_mask))
         cls.anat_sform = anat_im.get_affine()
         cls.anat = anat
         cls.anat_max = anat.max()
     return cls.anat, cls.anat_sform, cls.anat_max

    def from_labelled_image(self, image_path, label):
        """
        Define the ROI as the set of  voxels of the image
        that have the pre-defined label

        Parameters
        -----------
        image_path: ndarray
            a label (discrete valued) image
        label: int
            the desired label
        """
        self.check_header(image_path)
        nim = load(image_path)
        data = nim.get_data()
        self.discrete = np.where(data==label)

def save_masked_volume(data, mask_url, path, descrip=None):
    """
    volume saving utility for masked volumes
    
    Parameters
    ----------
    data, array of shape(nvox) data to be put in the volume
    mask_url, string, the mask path
    path string, output image path
    descrip = None, a string descibing what the image is
    """
    rmask = load(mask_url)
    mask = rmask.get_data()
    shape = rmask.get_shape()
    affine = rmask.get_affine()
    save_volume(shape, path, affine, mask, data, descrip)