for f_idx in v2f[v_idx]:
for node_idx in v2fn[f_idx][v_idx]:
// Sum the signal from each node of the fiber in that voxel:
f_matrix_sig[keep_ct:keep_ct+n_bvecs] = \
f_matrix_sig[keep_ct:keep_ct+n_bvecs] + fiber_signal[f_idx][node_idx]
// For each fiber-voxel combination, we now store the row/column
// indices in the pre-allocated linear arrays
After Change
// In each voxel:
for v_idx in xrange(vox_coords.shape[0]):
// dbg:
if not np.mod(v_idx, 10000):
print("voxel %s"%(100*float(v_idx)/n_vox))
//For each fiber in that voxel:
