2D surface searchlight classification using freesurfer

[Andreas Schindler]

Hello

I just read the searchlight tutorial and wondered if there is maybe a
hidden option for using a 2D circular spotlight mask to process
cortical surface data instead of the default 3d volume sphere / cube
approach.
As far as I read, such a surface driven classification approach has
mainly two advantages.
On the one hand by reconstructing the cortical surface one only puts
functionally adjacent voxels within the same mask (which isn't
guaranteed within the volume as a 3D sphere might contain voxels of
two neighbouring gyri). On the other hand by discarding white matter,
classification gains speed and maybe even accuracy as many
uninformative voxels are gone.

If I got it right, given the flexibility of mvpa toolbox such an
approach shouldn't be difficult to implement. As, instead of volume
data one just could feed in a surface time series by e.g. using
freesurfer's matlab extension. Further, a 2D circle based dependency
mask could be created manually extending the create_adj_list function.

However, as I am rather new to the mvpa toolbox I wanted to ask if
someone maybe has already some experiences with surface based
spotlight approaches. Further, I wanted to ensure myself if it is
really that simple to feed surface instead of volume data into the
toolbox... or if might have missed a crucial difference between the
two. :-)

Many thanks for any kind of hint

best wishes

Andreas

[Greg Detre]

superb question. i had never thought of doing that.

here's my knee-jerk, off-the-cuff, from-the-hip,
shoot-first-ask-questions-later, it-looked-like-he-was-going-for-his-gun
response. i've already killed a handful of mvpa toolbox questioners by
accident though, so read with caution.

if you can get your surface into matlab as a flattened sheet, rather
than a spherical surface, then this should be fairly straightforward.
here's what i'd try:

% your functional data, as surface, i.e. 2D+T (spatial_dimension_1 x
spatial_dimension_2 x time)
func2 = rand(nRows,nCols,nTimepoints);

% turn that into a standard 1D+T pattern matrix, i.e. collapse all the
spatial dimensions together
func1 = reshape(func2, [nRows*nCols nTimepoints]);

% create a 3D mask. this is really just a 2D surface, but embedded in a
3D volume.
% N.B. make sure to put the singleton dimension first,
% otherwise matlab will squeeze this down to a 2D matrix implicitly
surfmask = ones([1 nRows nCols]);
subj = initset_object(subj,'mask','surfmask',surfmask);

% load in the new pattern, using 'surfmask' as the mask
subj = initset_object(subj,'pattern','func1',func1, 'masked_by','surfmask');

i haven't checked this carefully, and you should make sure that the
reshape command is ordering the rows correctly, but hopefully this gives
you the gist of a solution.

N.B. since i don't work with surface data too much myself, superficial
people should pipe up if you have a better plan!

g

> --
>
> You received this message because you are subscribed to the Google Groups "Princeton MVPA Toolbox for Matlab" group.
> To post to this group, send email to mvpa-t...@googlegroups.com.
> To unsubscribe from this group, send email to mvpa-toolbox...@googlegroups.com.
> For more options, visit this group at http://groups.google.com/group/mvpa-toolbox?hl=en.
>
>
>


--


---
Greg Detre
cell: 617 642 3902
email: gr...@gregdetre.co.uk
web: http://www.gregdetre.co.uk

[MS Al-Rawi]

Hello Andreas 

Nice idea .....but.... does freesurfer's matlab extension yields cortical surface represented as 2D plane shape(s), I dont know about tha!t. In fact, a surface is generally defined as a 3D volume (e.g., the spher's surface). Now, let's forget the complicated cortical surface fractal shape and assume that we want to transform the locations at the sphere's surface into a 2D plane. This requires the use of an interpolation algorithm which causes digitization error. Thus, if (hypothetically speaking) a 2D plane shape cortical surface do exist, it'll be contaminated with noise due to digitization errors and aliasing (see any ref. on sampling theory).

This brings us back to the other choice which is a 3D cortical surface represented by its original coordinate system. This cortical surface could be used directly with MVPA by masking out voxels not belonging to the 3D cortical surface (i.e., reading a mask that represent the cortical surface). And since the 3D cortical surface has many zeor valued voxels we can speedup the computation of searchlight easily (please see my other contribution 'A faster ANOVA' ) and I myself gained a very fast searchlight implementation by making use of this logic as well as the one proposed by Raj (Suggested small change to searchlight code, produces very big speed-up).

Going back to the original notion proposed by you which is using a 2D circular searchlight mask to define the neighbored voxels. I understand that what you want is to consider points that are neighbors at the cortical surface and you want to see what effect they might give. Sounds nice, however, if I am not mistaken, what I described above minimizes the chance of using a 2D plane shape cortical surface and what we really have is a 3D cortical surface, thus,  we cannot implement a 2D circular searchlight.

Is there a way around this? I think there is. To consider neighbored points at the surface we need to generate (in fact we can) adjacent points at the surface using adjacency rules, for that you need to go to the book (digital image processing by Gonzalez and Woods- http://www.imageprocessingplace.com/). You can also forget the surface adjacency and use the one (3D adjacency) that MVPA provides that will give close results to surface based adjacency when the radius is small. 

Regards

Al-Rawi 

IEETA

Portugal

[Ben Singer]

Haven't done it, but have worked a bit with freesurfer meshes w/epi data at each node in a different context. An idea that comes to mind is to create a volume that is the bounding cube of the Freesurfer spherical mesh and then fill in the voxels within the spherical mask with the values they had in their original location pre-inflation. Then do the normal searchlight within this new 3D space (being sure not to waste time in the empty middle), so the adjacency is more like the 2D cortical adjacency, but with some distortion introduced by the inflation/warping etc. Since Freesurfer creates a spherical mesh with a 100mm radius there would be 200^3 1mm-cube voxels, quite a few! but using something between 1 and 3mm makes more sense given epi data resolution. That's lots of voxels, but you would create a spherical mask-- plus after you map epi data to each mesh node via SUMA's 3dVol2Surf at the 1mm native resolution, you could smooth the epi data on the sphere to a lower mesh resolution with SUMA's SurfSmooth. Then for each voxel in the spherical mask, assign a value from the mesh nodes that fall within it (which would be quite a few the lower the resolution you use). You'd have to do that manually since 3dSurfToVol would fill in the voxels in their original curvy locations.

There's lots of averaging going on here, though-- 3dVol2Surf averages (independently at each TR) voxels (usually 3 or so) to get the mesh node's value (though you can ask to use functions other than the mean)-- then you are averaging the nodes that fall within each of your new voxels. This could be avoided by instead getting the mapping between spherical voxels and original voxels.

Hope that makes sense-- if there's any sense to made from it, since I was mostly thinking out loud!

Ben

[Ben Singer]

Come to think of it you wouldn't need to do any operations on the epi data at all-- just generate a new coords file and mask, so each original voxel has a corresponding post-warp, on-the-sphere 3D coords, if it is a cortical voxel (not masked out). I might try this as an exercise if people think it would be worth doing, since I deal mostly with cortical voxels anyway.

[Andreas Schindler]

Hi Greg and Al-Rawi

First of all many thanks for your inspiring posts!
Meanwhile I also figured out that there does not seem to be an easy
way to just import freesurfer’s surfaces as 2D matrices into matlab.
The function allowing surface import unfortunately only provides a
column of vertex values where adjacency information of neighbouring
vertices is not preserved.

Apart from that, you are right A-Rawi and surfaces generally are
defined as 3D shapes in freesurfer. Nevertheless, as there are also
ways to flatten these spherical surfaces to so called flat patches
which really constitute 2D meshes, I thought it might be worth working
with these.
However, I have to admit that I didn’t think too much about possible
interpolation noise arising from such a flattening process.

For that reason I really like your approach of masking the functional
volumes by a cortical surface mask as might for example be provided by
freesurfer’s mri_surf2vol function. However, by traveling through
such a masked volume using a sphere spotlight there still remains the
“danger” of sampling across neighbouring gyri (which contain voxels
distant from each other on the cortical surface).

A sort of compromise might be the solution I recently found suggested
by Soon et al.

http://www.sciencedirect.com/science?_ob=ArticleURL&_udi=B6WNP-4X3PHYG-NN&_user=29041&_rdoc=1&_fmt=&_orig=search&_sort=d&_docanchor=&view=c&_searchStrId=1139727845&_rerunOrigin=google&_acct=C000003178&_version=1&_urlVersion=0&_userid=29041&md5=4a2c847736fa8d5f3df1d5168518a819

As far as I understand they are using the inflated surface to define
their searchlight ROIs and project them back into the EPI volume
space. As the inflated surface isn’t flat they use a circular shaped
spotlight with a geodesic approximated radius (even though I haven’t
figured out how to realise that with freesurfer). That way I think one
cannot exclude interpolation errors arising from inflation (though
there are surely less than as for entire flattening) but at least can
avoid the aforementioned sampling of neighbouring gyri.

Looking forward to your comments / suggestions

Best wishes

Andreas

On Dec 17, 4:42 pm, MS Al-Rawi <rawi...@yahoo.com> wrote:
> Hello Andreas
>
> Nice idea .....but.... does freesurfer's matlab extension yields cortical surface represented as 2D plane shape(s), I dont know about tha!t. In fact, a surface is generally defined as a 3D volume (e.g., the spher's surface). Now, let's forget the complicated cortical surface fractal shape and assume that we want to transform the locations at the sphere's surface into a 2D plane. This requires the use of an interpolation algorithm which causes digitization error. Thus, if (hypothetically speaking) a 2D plane shape cortical surface do exist, it'll be contaminated with noise due to digitization errors and aliasing (see any ref. on sampling theory).
>
> This brings us back to the other choice which is a 3D cortical surface represented by its original coordinate system. This cortical surface could be used directly with MVPA by masking out voxels not belonging to the 3D cortical surface (i.e., reading a mask that represent the cortical surface). And since the 3D cortical surface has many zeor valued voxels we can speedup the computation of searchlight easily (please see my other contribution 'A faster ANOVA' ) and I myself gained a very fast searchlight implementation by making use of this logic as well as the one proposed by Raj (Suggested small change to searchlight code, produces very big speed-up).
>
> Going back to the original notion proposed by you which is using a 2D circular searchlight mask to define the neighbored voxels. I understand that what you want is to consider points that are neighbors at the cortical surface and you want to see what effect they might give. Sounds nice, however, if I am not mistaken, what I described above minimizes the chance of using a 2D plane shape cortical surface and what we really have is a 3D cortical surface, thus,  we cannot implement a 2D circular searchlight.
>

> Is there a way around this? I think there is. To consider neighbored points at the surface we need to generate (in fact we can) adjacent points at the surface using adjacency rules, for that you need to go to the book (digital image processing by Gonzalez and Woods-http://www.imageprocessingplace.com/). You can also forget the surface adjacency and use the one (3D adjacency) that MVPA provides that will give close results to surface based adjacency when the radius is small.
>
> Regards
>
> Al-Rawi
> IEETA
> Portugal

[Greg Detre]

if anyone is interested in directly manipulating the MVPA toolbox
internal data structures to define which voxels are adjacent to which,
you can create your own ADJ_LIST (see CREATE_ADJ_LIST.M and what it
calls/creates). then you could just feed your own adjacency matrix in,
and you'd have complete control over which voxels would be considered
neighboring.

i'd be happy to talk this over in more detail.

g

> errors and aliasing (see any ref. on sampling theory*).*

> <mailto:mvpa-t...@googlegroups.com>.

> To unsubscribe from this group, send email to
> mvpa-toolbox...@googlegroups.com

> <mailto:mvpa-toolbox%2Bunsu...@googlegroups.com>.

> For more options, visit this group at
> http://groups.google.com/group/mvpa-toolbox?hl=en.
>
>
>
> --
>
> You received this message because you are subscribed to the Google
> Groups "Princeton MVPA Toolbox for Matlab" group.
> To post to this group, send email to mvpa-t...@googlegroups.com.
> To unsubscribe from this group, send email to
> mvpa-toolbox...@googlegroups.com.
> For more options, visit this group at
> http://groups.google.com/group/mvpa-toolbox?hl=en.

[MS Al-Rawi]

Hi 

Ben .....

I got your idea, you are saying that it is possible to wrap (warp or fold) the 3D surface into a 2D plane like surface. But, you said, "There's lots of averaging going on here"  to avoid having multi-layer 2D like surfaces after wrapping, I guess. This leaves us with the same interpolation and aliasing error(s), please correct if I am wrong.

Andreas ..you noted

> For that reason I really ......(which contain voxels  distant from each other on the cortical surface).

You are definite about that, thus, I urge you to consider the other option, and it is very easy to implement. Just go read Chapter2 (2.2.2. Adjacency, Connectivity, Regions, and Boundaries) of digital image processing by Gonzalez and Woods 3d ed- http://www.imageprocessingplace.com/  .. and you will see that it is easy to develop an algorithm to ensure that the neighborhood voxels on the cortical surface contain non-distant voxels from each other. If you succeed with this algorithm, you just have to replace its function (program you write for Adjacency cortical surface) instead of the 3D sphere of the MVPA searchlight.

Good luck

Al-Rawi

[Ben Singer]

Hi Al-Rawi,

Sort of... what I was suggesting really boiled down to creating a new 3D coords file and using existing adjacency/searchlight code. A quick and dirty way of doing something like the Soon et al paper. In both cases you want to project back the nodes on the spherical freesurfer surface to their position in the original volume. That's easy-- freesurfer surfaces are all in node-to-node correspondence, so row 1 is node1's 3D coordinates, in the original folded surface file (the "white" file) and in the inflated sphere file ("sphere" or "sphere.reg" file). Now find which voxel the node came from. If you assign the "sphere.reg" 3D coordinates to that voxel (actually there will be multiple voxels due to cortical thickness, but never mind that), you are giving the voxel coordinates in a standard anatomical cortical template, and you can combine data across subjects before doing the adjacency. The 3D searchlight intersecting with in essence a 2D surface isn't as nice as the 2D searchlight of Soon et al, but it lets you get up and running with only a new coords file.

Ben

[Ben Singer]

(you will also need to align the epi volume with the high-res anatomical "surfvol" volume for that subject, that was used by freesurfer to create the surface-- the description below will give coords in the original surfvol but for the epi voxel coords you'll need to run suma's @AlignToExperiment script)