def MRSvoxelStats(segfile, MRSfile=None, center=None, dim=None, subjID=None,
                  gareas=GAREAS,wareas=WAREAS,csfareas=CSFAREAS):
    """
    returns grey/white/CSF content within MRS voxel

    Parameters
    ----------
    segfile:  nifti file
        path to segmentation file with grey/white matter labels (freesurfer
        aseg file converted from mgz). 

    MRSfile: nifti file
        path to MRSfile of MRS voxel. provide either this or center + dim

    center : integer array
        [x,y,z] where x, y and z are the coordinates of the point of
        interest. Provide either MRSfile or center+dim 
    
    dim : float array
        dimensions of voxel in mm. Provide either MRSfile or center+dim

    subjID: string
        optional subject identifier. Defaults to nims scan number

    gareas, wareas, csfareas: arrays of integers
        arrays of freesurfer labels for gray, white, and csf areas
        respectively. Determines which areas are considered
        grey/white/csf. Some areas like brainstem were lumped with grey
        matter. Did not include hypointensities. 
   
    Returns
    -------

    raw number of voxels and proportions of grey, white, csf and
    non-grey-or-white-matter 

    Notes
    -----
    See http://miykael.github.io/nipype-beginner-s-guide/regionOfInterest.html

    """
    # subject ID
    if subjID==None:
        m=re.search('(\w+?)_*_',segfile)
        subjID=m.group(0)[:-1] 

    # get segmentation file
    aseg = nib.load(segfile)
    aseg_data = aseg.get_data().squeeze()
    aseg_aff = aseg.get_affine()
    segdir = os.path.dirname(segfile)
    segvoxdim=np.diagonal(aseg_aff)[:3]

    # get nifti file of MRS voxel if one is provided
    if MRSfile is not None:
        if center is not None or dim is not None:
            msg = 'provide EITHER MRSfile OR center and dim, not both!'
            raise ValueError(msg)

        mrs = nib.load(MRSfile)
        mrs_data = mrs.get_data().squeeze()
        mrs_aff = mrs.get_affine()        
        # This applies the concatenation of the transforms from mrs space to
        # the T1 space. [0,0,0] is the center of the MRS voxel:
        center = np.round(np.dot(np.dot(np.linalg.pinv(aseg_aff), mrs_aff),
                                 [0,0,0,1]))[:3].astype(int)
        
        dim=np.diagonal(mrs_aff)[:3]
    else: # no MRSfile
        if center==None or dim==None:
            msg = 'if no MRSfile is provided, provide center and '
            msg += 'dimensions of voxel'
            raise ValueError()    

    # calculate beginning corner of MRS voxel
    print('Creating mask with center: [%s, %s, %s]'%(center[0],
                                                    center[1],
                                                    str(center[2])))

    print('and dimensions: [%s, %s, %s] mm'%(dim[0], dim[1], dim[2]))
    
    # calculate roi mask with numpy
    # round up to nearest number of voxel units
    voxdim=np.zeros(3)
    voxdim[0]=np.ceil(np.abs(dim[0] / segvoxdim[0]))
    voxdim[1]=np.ceil(np.abs(dim[1] / segvoxdim[1]))
    voxdim[2]=np.ceil(np.abs(dim[2] / segvoxdim[2]))

    print ('MRS voxel dimensions in T1 voxel units: [%s, %s, %s]'%(voxdim[0],
                                                                  voxdim[1],
                                                                  voxdim[2]))

    # corners
    lcorner = [center[0] - voxdim[0] / 2,
               center[1] - voxdim[1] / 2,
               center[2] - voxdim[2] / 2]

    ucorner = [center[0] + voxdim[0] / 2,
               center[1] + voxdim[1] / 2,
               center[2] + voxdim[2] / 2]
    # create mask
    mdata = np.zeros(aseg_data.shape)
    for i in range(int(lcorner[0]),int(ucorner[0])):
        for j in range(int(lcorner[1]),int(ucorner[1])):
            for k in range(int(lcorner[2]),int(ucorner[2])):
                mdata[i,j,k]=1


    # calculate grey/white/csf from freesurfer labels
    gdata = np.zeros(aseg_data.shape)
    wdata = np.zeros(aseg_data.shape)
    csfdata = np.zeros(aseg_data.shape)

    for data, areas in zip([gdata, wdata, csfdata], [gareas, wareas, csfareas]):
        for area in areas:
            data[np.where(aseg_data==area)] = 1

    # multiply voxel ROI with seg data
    gmasked= mdata * gdata
    wmasked= mdata * wdata
    csfmasked= mdata * csfdata
    
    # extract stats from a given segmentation
    total = np.sum(mdata)
    white = np.sum(wmasked)
    grey = np.sum(gmasked)
    csf = np.sum(csfmasked)
    nongmwm = total - grey - white

    # proportions
    pWhite = float(white) / total
    pGrey = float(grey) / total
    pCSF = float(csf) / total 
    pNongmwm = float(nongmwm) / total

    return (total, grey, white, csf,nongmwm, pGrey, pWhite,pCSF, pNongmwm)