Unmask CIVET 2 data in python - CoBrALab/documentation GitHub Wiki

Vertex-wise CIVET outputs contain data at 40962 vertices for each hemisphere. A portion of these vertices are located on the midline, where CIVET estimates are unreliable. If you are conducting analyses such as PLS or NMF with CIVET data, you want to remove those invalid vertices prior to your analysis. Then, when mapping outputs back to CIVET files you want assign 0 to the invalid vertices. The code snippets below can help do that using the CIVET 2 mask file and some python code.

This requires python, and numpy. The below is written for python 3. It should work in python 2 as well, just edit or remove the print statements. The code below is not a straight plug and play. Edit it to point to your files. The code below is operating on only 1 subjects data (only loads one .txt file). If you are operating with multiple subjects, you may be concatenating each subjects data together and your dimensions will be slightly different (eg instead of (1,n_vertex), you'll have (n_subjects,n_vertex) perhaps). But the below should be adjustable for most cases.

import numpy as np

##Load data, mask

#load in your CT data (or SA or other vertex data)
left_ct = np.loadtxt('raw_ct.txt') 
print(np.shape(left_ct)) #should be (40962,)

#load in the CIVET2 left hemisphere mask
#this file has is 0 in the masked area of civet (midline), and 1 elsewhere
left_mask = np.loadtxt('/opt/quarantine/resources/CIVET/CIVET_2.0_mask_left_short.txt') #this is at the CIC
left_valid = np.where(left_mask==1) # list of valid indices in civet .txt file
left_invalid = np.where(left_mask==0) #list of invalid indices in civet .txt file
print(np.shape(left_valid)) #should be (1, 38561)

left_ct_valid = left_ct[left_valid] #create a new variable, left_ct_valid, that only contains valid (not midline) data
print(np.shape(left_ct_valid)) #should be (38561,)

##do stuff. pls, nmf etc, where you perform some analysis on your valid vertex data
#say you get out vertex_results from this, which should still have dimensions (38561,) or similar 
#number of columns can be different if you concatenated multiple subjects data together at somepoint, but num_vertices should be same
#vertex_results = vertex-wise results from your analysis


##now map back to original dimensions and write out

vertex_outarray = np.zeros((np.shape(left_mask)[0],1)) #initialize an array w 40962 rows. 1 row for every civet vertex, including midline
print(np.shape(vertex_outarray)) #should be (40962, 1)

valid_idx = 0 #counter to keep track of things
#cycle through every element of your left_mask file
#remember this has length 40962, 0 in masked area, 1 elsewhere
#if left_mask = 1, you're at a valid vertex, copy data to your out array and increment valid_idx
#if left_mask = 0, dont copy, out array maintains 0 at that position
for idx in range(0,np.size(left_mask)):
    if left_mask[idx] == 1:
        vertex_outarray[idx,0] = vertex_results[valid_idx]
        valid_idx +=1
        
np.savetxt('processed_ct.txt',vertex_outarray) #write out your results