File: multi_subject_parcelation.py

package info (click to toggle)
nipy 0.1.2%2B20100526-2
  • links: PTS, VCS
  • area: main
  • in suites: squeeze
  • size: 11,992 kB
  • ctags: 13,434
  • sloc: python: 47,720; ansic: 41,334; makefile: 197
file content (66 lines) | stat: -rw-r--r-- 1,948 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
"""
This script contains a quick demo on  a multi'subject parcellation
on a toy 2D example.
Note how the middle parcels adapt to the individual configuration.
"""
print __doc__

import numpy as np
import nipy.neurospin.spatial_models.hierarchical_parcellation as hp
import nipy.neurospin.utils.simul_multisubject_fmri_dataset as simul
import nipy.neurospin.spatial_models.parcellation as fp

# step 1:  generate some synthetic data
nsubj = 10
dimx = 60
dimy = 60
pos = 3*np.array([[ 6,  7],
                  [10, 10],
                  [15, 10]])
ampli = np.array([5, 7, 6])
sjitter = 6.0
dataset = simul.surrogate_2d_dataset(nbsubj=nsubj, dimx=dimx, dimy=dimy, 
                                     pos=pos, ampli=ampli, width=10.0)
# dataset represents 2D activation images from nsubj subjects,
# with shape (dimx,dimy)

# step 2 : prepare all the information for the parcellation
nbparcel = 10
ref_dim = (dimx,dimy)
xy = np.array(np.where(dataset[0])).T
nvox = np.size(xy,0)
xyz = np.hstack((xy,np.zeros((nvox,1))))
	
ldata = np.reshape(dataset,(nsubj,dimx*dimy,1))
anat_coord = xy
mask = np.ones((nvox,nsubj)).astype('bool')
Pa = fp.Parcellation(nbparcel,xyz,mask-1)

# step 3 : run the algorithm
Pa =  hp.hparcel(Pa, ldata, anat_coord, mu = 3.0)
# note: play with mu to change the 'stiffness of the parcellation'
	
# step 4:  look at the results
Label =  np.array([np.reshape(Pa.label[:,s],(dimx,dimy))
                   for s in range(nsubj)])

import matplotlib.pylab as mp
mp.figure()
mp.title('Input data')
for s in range(nsubj):
    mp.subplot(2, 5, s+1)
    mp.imshow(dataset[s], interpolation='nearest')
    mp.axis('off')

mp.figure()
mp.title('Resulting parcels')
for s in range(nsubj):
    mp.subplot(2, 5, s+1)
    mp.imshow(Label[s], interpolation='nearest', vmin=-1, vmax=nbparcel)
    mp.axis('off')
mp.show()