[pymvpa] Using searchlight for quick spherical ROIs?

[Michael Hanke]

Hi Johan,

On Sun, Mar 15, 2009 at 03:43:51PM +0000, Johan Carlin wrote:
> Hi all,
> 
> I have just started using PyMVPA for some fMRI analysis. It has worked like
> a charm so far. :) Thanks very much for developing such a useful tool.

I'm glad you like it!

> I would like to run a clasifier on a number of ROIs that were defined in a
> univariate SPM analysis. The most obvious solution is to re-load the dataset
> with a new mask for each ROI. However, I wonder if there is a quicker
> option, using the Searchlight and the "center id" parameter?
> 
> The idea would be to start with a list of peak voxel coordinates for each
> ROI from the nifti image, map these to feature IDs somehow, and pass the
> list of ROI IDs to the Searchlight. The result would then be a (fast!)
> searchlight analysis, that runs a spherical classifier centered on each ROI
> in turn. The resulting classification map should only contain one voxel per
> ROI.
> 
> Has anyone attempted something similar? I can't quite figure out how to map
> the ROI coordinates in the original nifti to the feature ID that the
> classifier runs on.

I cannot look into the issue in detail right now, but with respect to
the "translation" problem I might provide some help. Let's look at the
following code. The first snippet loads the example dataset, which shall
for now be a replacement for your dataset -- the important bit is
however the type being 'NiftiDataset':

  >>> from mvpa.suite import *
  >>> attr = SampleAttributes(os.path.join(mvpa.pymvpa_dataroot,'attributes.txt'))
  >>> ds = NiftiDataset(samples=os.path.join(mvpa.pymvpa_dataroot,'bold.nii.gz'),
  ...                   labels=attr.labels, chunks=attr.chunks)


The dataset apparently has 1452 features

  >>> ds
  <Dataset / int16 1452 x 800 uniq: 12 chunks 9 labels>

If you now want to know which feature corresponds to a particular
coordinate, you can simply query the mapper instance of the dataset.

  # feature id for coordinate z=0 y=10 x=20
  >>> ds.mapper.getOutId((0,10,20))
  420

That should get you started. You can even go a bit further and compile
the list of neighbouring voxels yourself (again with the mappers help):

  >>> ds.mapper.getNeighbors(420, radius=5)
  [379, 380, 381, 419, 420, 421, 459, 460, 461]

which return you a list of feature ids all in a 5mm range around "420".


Does that help?


Michael

-- 
GPG key:  1024D/3144BE0F Michael Hanke
http://apsy.gse.uni-magdeburg.de/hanke
ICQ: 48230050