Do some quick functional density mapping for a sample voxel. Note this code assumes you are running off rocky.

easy_functional_density.py

import numpy as np
import os
#from CPAC.network_centrality import load
#from CPAC.cwas.subdist import norm_cols

import nibabel as nib
from scipy.sparse import coo_matrix, cs_graph_components

def load(datafile, template):
    
    """
    Method to read data from datafile and mask/parcellation unit
    and store the mask data, timeseries, affine matrix, mask type
    and scans. The output of this method is used by all other nodes.
    
    Parameters
    ----------
    datafile : string (nifti file)
        path to subject data file
    template : string (nifti file)
        path to mask/parcellation unit
        
    Returns
    -------
    timeseries_data: string (numpy npy file)
        path to file containing timeseries of the input data
    affine: string (numpy npy file)
        path to file containing affine matrix of teh input data
    mask_data: string (numpy npy file)
        path to file containing mask/parcellation unit matrix
    template_type: string 
        0 for mask, 1 for parcellation unit 
    scans: string (int)
        total no of scans in the input data
        
    Raises
    ------
    Exception
    """

    import os
    import nibabel as nib
    import numpy as np

    
    try:    
        if isinstance(datafile, list):
            img = nib.load(datafile[0])
        else:
            img = nib.load(datafile) 
        
        data = img.get_data().astype(np.float32)
        aff = img.get_affine()    
        mask = nib.load(template).get_data().astype(np.float32)  
        scans = data.shape[3]
        
    except:
        print "Error in loading images for graphs"
        raise
    
    
    
    if mask.shape != data.shape[:3]:
        raise Exception("Invalid Shape Error. mask and data file have"\
                        "different shape please check the voxel size of the two files")

    
    #check for parcellation
    nodes = np.unique(mask).tolist()
    nodes.sort()
    print "sorted nodes: ", nodes
    
    #extract timeseries
    if len(nodes)>2:
        flag=1
        for n in nodes:
            if n > 0:
                node_array = data[mask == n]
                avg = np.mean(node_array, axis =0)
                if flag:
                    timeseries = avg
                    flag=0
                else:
                    timeseries = np.vstack((timeseries, avg))
        #template_type is 1 for parcellation
        template_type = 1
    else:
        #template_type is 0 for mask
        template_type = 0
        mask = mask.astype('bool')
        timeseries = data[mask]
        
    #saving files
    img_name = os.path.splitext(os.path.splitext(os.path.basename(datafile))[0])[0] + '.npy'
    mask_name = os.path.splitext(os.path.splitext(os.path.basename(template))[0])[0]  + '.npy'
    
    mask_data = os.path.join(os.getcwd(), mask_name)
    timeseries_data = os.path.join(os.getcwd(),img_name)
    affine = os.path.join(os.getcwd(),'affine.npy')
        
    np.save(mask_data, mask)
    np.save(timeseries_data, timeseries)
    np.save(affine, aff)
    
    return timeseries_data, affine, mask_data, template_type, scans


def norm_cols(X):
    """
    Theano expression which centers and normalizes columns of X `||x_i|| = 1`
    """
    Xc = X - X.mean(0)
    return Xc/np.sqrt( (Xc**2.).sum(0) )


# Borrowed from nipy.graph.graph
# https://github.com/nipy/nipy/blob/master/nipy/algorithms/graph/graph.py
def graph_3d_grid(xyz, k=18):
    """ Utility that computes the six neighbors on a 3d grid

    Parameters
    ----------
    xyz: array of shape (n_samples, 3); grid coordinates of the points
    k: neighboring system, equal to 6, 18, or 26

    Returns
    -------
    i, j, d 3 arrays of shape (E),
            where E is the number of edges in the resulting graph
            (i, j) represent the edges, d their weights
    """
    if np.size(xyz) == 0:
        return None
    lxyz = xyz - xyz.min(0)
    m = 3 * lxyz.max(0).sum() + 2

    # six neighbours
    n6 = [np.array([1, m, m ** 2]), np.array([m ** 2, 1, m]),
         np.array([m, m ** 2, 1])]

    # eighteen neighbours
    n18 = [np.array([1 + m, 1 - m, m ** 2]),
           np.array([1 + m, m - 1, m ** 2]),
           np.array([m ** 2, 1 + m, 1 - m]),
           np.array([m ** 2, 1 + m, m - 1]),
           np.array([1 - m, m ** 2, 1 + m]),
           np.array([m - 1, m ** 2, 1 + m])]

    # twenty-six neighbours
    n26 = [np.array([1 + m + m ** 2, 1 - m, 1 - m ** 2]),
           np.array([1 + m + m ** 2, m - 1, 1 - m ** 2]),
           np.array([1 + m + m ** 2, 1 - m, m ** 2 - 1]),
           np.array([1 + m + m ** 2, m - 1, m ** 2 - 1])]

    # compute the edges in each possible direction
    def create_edges(lxyz, nn, l1dist=1, left=np.array([]), right=np.array([]),
                     weights=np.array([])):
        q = 0
        for nn_row in nn:
            v1 = np.dot(lxyz, nn_row)
            o1 = np.argsort(v1)
            sv1 = v1[o1]
            nz = np.squeeze(np.nonzero(sv1[: - 1] - sv1[1:] == - l1dist))
            o1z, o1z1 = o1[nz], o1[nz + 1]
            left = np.hstack((left, o1z, o1z1))
            right = np.hstack((right, o1z1, o1z))
            q += 2 * np.size(nz)
        weights = np.hstack((weights, np.sqrt(l1dist) * np.ones(q)))
        return left, right, weights

    i, j, d = create_edges(lxyz, n6, 1.)
    if k >= 18:
        i, j, d = create_edges(lxyz, n18, 2, i, j, d)
    if k == 26:
        i, j, d = create_edges(lxyz, n26, 3, i, j, d)
    i, j = i.astype(np.int), j.astype(np.int)

    # reorder the edges to have a more standard order
    order = np.argsort(i + j * (len(i) + 1))
    i, j, d = i[order], j[order], d[order]
    return i, j, d

def cluster_data(img, thr, mask):
    """docstring for cluster_data"""
    xyz     = np.where(mask > 0)
    xyz_a   = np.array(xyz).T
    xyz_th  = xyz_a[img[xyz] > thr]
    labels  = clusterize(xyz_th)
    return labels

def clusterize(xyz_th, k=26):
    """docstring for clusterize"""
    nvoxs       = xyz_th.shape[0]
    i,j,d       = graph_3d_grid(xyz_th, k=k)    
    adj         = coo_matrix((d, (i,j)), shape=(nvoxs,nvoxs))
    nc, labels  = cs_graph_components(adj)
    return labels

thresh = 0.1
datafile = "/data/Projects/ABIDE_Initiative/CPAC/test_qp/Centrality_Working/resting_preproc_0051466_session_1/_mask_mask_abide_90percent_gm_4mm/_scan_rest_1_rest/resample_functional_to_template_0/residual_wtsimt_flirt.nii.gz"
maskfile = "/home2/data/Projects/ABIDE_Initiative/CPAC/abide/templates/masks/mask_abide_90percent_gm_4mm.nii.gz"

timeseries, aff, mask, template_type, ntpts = load(datafile, maskfile)

# For this example, calculate connectivity for 10 voxels each with the rest of the brain
j=0; i=9
timeseries = norm_cols(timeseries.T)
corr_matrix = np.dot(timeseries[:,j:i].T, timeseries) # corr_matrix is 10 voxs by nvoxs

# Let's get the functional density for one voxel (the 1st one)
cur_vox = 0
## we put the correlation values into 3D space
vox_map = np.zeros(mask.shape)
vox_map[mask] = corr_matrix[cur_vox,:]
## then we do a somewhat redundant step of getting clusters throughout brain
## in future can modify the function to only look within the neighborhood
labels = cluster_data(vox_map, thresh, mask) # if label -2 then isolate?
## density is the sum of all regions in the same cluster as our voxel
local_density = np.sum(labels==labels[cur_vox]) # I believe labels should only have nvoxs because we specified a mask earlier
                                                # should check this
print 'value -> ', local_density