regarding morphing in matlab

Dear MNE users,

Just a follow up mail to my previous request. I would like to transform my
data from one decimated source space to another from within matlab. Below
are the steps I have adopted to do so. If anybody has previous experience
in doing similar things, kindly share your thoughts.

1/ I precomputed the morph maps from each surface to the common brain with
mne_make_morph_maps.

2/ I read in the morph maps using mne_read_morph_maps as follows:
[leftmap,rightmap] = mne_read_morph_map('my_brain', 'common_brain');

3/ I read in an inverse operator for the input surface 'my_brain' into the
structure 'inv_op'

4/ I read in a sample stc file for the target surface i.e. 'common_brain'
to get the vertex info. The stc files were read into 'stc_l' and 'stc_r'

5/ Next, I used the following lines of code to transform my data:
leftmapRelevant = leftmap(stc_l.vertices,inv_op.src(1).vertno);
rightmapRelevant = rightmap(stc_r.vertices,inv_op.src(2).vertno);
morphed_dataL = leftmapRelevant*dataL;
morphed_dataR = rightmapRelevant*dataR;

Note that stc_l.vertices contains all used vertices for the left
hemishpere of the target surface (i.e. 'common_brain') and
inv_op.src(1).vertno contains all used vertices for the left hemisphere of
the source surface (i.e. 'my_brain'). The same is true for stc_r.vertices
and inv_op.src(2).vertno respectively.

Intuition suggests that these steps would suffice, but I did some checks
and figured that most rows of leftmapRelevant and rightmapRelevant are
entirely zero.

I am using a morphgrade of 4 i.e. target surface consists of 2562 used
vertices per hemisphere. By comparison, the source surfaces are computed
with 5mm spacing and therefore consist of about 5000 used vertices per
hemisphere.

What am I missing?

Thank you for your time.
Best regards,
Pavan

Hi Pavan,
   I don't have much experience with morphing data in MATLAB. Given that,
the procedure looks ok to me. Did you write any of your morphed data to
an stc file to see what it looks like?

There are a few possibilities that occur to me:
(1) Are you sure 'most rows' of leftmapRelevant and rightmapRelevant are
infact zero? They ought to be sparse with probably 2 or 3 elements
non-zero in the entire row of 5000 (if that's your source space size) odd
elements..

(2) Depending on how you made the sample stc file that is in
'common_brain' space, it might have lost track of the vertex numbering.
Could you look at stc_l.vertices and stc_r.vertices to see if they make
sense? If they are simply [0:2561]' or [0:10241]' for instance, the
procedure of course will give you mostly zeros.

Hope it helps.

Regards,
Hari

Hi Hari (and others),

Thanks very much for the suggestions and apologies for the lengthy reply!

I did check that most rows are entirely zero (about 2300/2562). Secondly,
I did look at the vertex numbering and they seem to make sense. Also, I
visualized the result as stc and they do seem as sparse as the data matrix
suggests.

It seems to me that smoothing (or as Matti recommends to call it in
Section 8.3 of the manual v2.6.1 : smudging/blurring) is the root of the
difference between the mne_make_movie result and my MATLAB implemenation.

For example, I checked that the mne_make_movie result with only a single
smoothing iteration (i.e. --smooth 1) is nearly as sparse. Specifically,
in the resulting stc file 2372/2562 used vertices in stc_l.data are zero.

Here are some snapshots from a a single time frame of a single subject stc
file along with the same file morphed to a common brain with 3 different
smoothing factors (1, 2 and 5). Matti recommends a smoothing factor
between 4-7.

https://neuro.hut.fi/~pavan/temp/smudging_subject01-lh-lat.png
https://neuro.hut.fi/~pavan/temp/smudging01_common-lh-lat.png
https://neuro.hut.fi/~pavan/temp/smudging02_common-lh-lat.png
https://neuro.hut.fi/~pavan/temp/smudging05_common-lh-lat.png

As you might guess, the sparseness of the respective data matrices
decreases with increased smoothing.

I would appreciate your (or anyone else's) comments on whether the
smoothing assessment above is right and if so, what would your
recommendations be regarding the choice of the smoothing factor. Please
note that I am subjecting my morphed data to a substantial analysis
pipeline and it is not a last stage visualization of group data as is
typically the case.

Also, any comments regarding implementation of the smoothing/blurring
operation in MATLAB are welcome. Dealing with a large "vertex x vertex"
distance matrix, as well as the choice of the neighborhood parameter 'N_j'
as defined in Section 8.3 of the v2.6.1 manual are not straightforward to
me.

Thanks again,
Pavan

You're right about the problem being because of the sparseness of the map.
We had the same problem for utilizing the MNE morph maps. To solve this I
implemented a nearest neighbor interpolation onto the hires freesurfer mesh,
then passed the morphing through that operator onto the low resolution mesh
for another subject.

I've appended the code I use to do this mapping below. This code is a
custom function that assumes a bunch of things about where data lives so you
won't be able to run it. It also depends on a function from the Stanford
VISTASOFT tools. The function is nearpoints, a mex implementation of a
nearest neighbor mapping. Nearpoints is extremely fast so it works well for
this largish datasets, so it is a generically useful function. Hopefully
this code gives you the idea of what is going on.

The total transformation is a combination of mappings. I take the From
subject and map it onto the hi-res FS surface. I then apply the hi res
morphing of the FROM subject to the TO subject. And finally I take the hi
res to lo res mapping on the to subject:

The key line is this one:
totalTrans{iHemi} = toHi2Lo*(hi2hi{iHemi}*fromLo2Hi);

Justin Ales

function [mapMtx totalTrans] = makeDefaultCortexMorphMap(fromSub,toSub)
%function [mapMtx] = makeDefaultCortexMorphMap(fromSub,toSub)
%
% This function makes a morphing matrix that maps values from 1 subject to
another
% Freesurfer based MNE morphing matrices must have been computed already
% with command:
% mne_make_morph_maps --to skeri0055_fs4 --from skeri0001_fs4
%
%example usage:
%[mapMtx] = makeDefaultCortexMorphMap('skeri0001','skeri0055');

% $Log: makeDefaultCortexMorphMap.m,v $
% Revision 1.4 2009/06/26 21:12:30 ales
% mrSimScript does arbitrary waveforms
% skeriDefaultSimParameters has updated help
% makeDefaultCortexMorphMap returns an identity matrix for mapping an
individual on to itself
% makeForwardMatrixFromMne has some added help
%
% Revision 1.3 2009/05/29 16:58:14 ales
% Changed something. I don't remember what now. sorry.
%
% Revision 1.2 2009/05/21 17:04:23 ales
% Added auto log message comments
%

% Read stuff in
toSub = [toSub '_fs4'];
fromSub = [fromSub '_fs4'];

freesurfDir = getpref('freesurfer','SUBJECTS_DIR');

toSubSrcSpaceFile = fullfile(freesurfDir,toSub,'bem',[ toSub
'-ico-5p-src.fif']);
fromSubSrcSpaceFile = fullfile(freesurfDir,fromSub,'bem',[ fromSub
'-ico-5p-src.fif']);

if strcmp(fromSub,toSub)
    toSubSrc = mne_read_source_spaces(toSubSrcSpaceFile);
    mapMtx = speye(sum([toSubSrc.nuse])); %<- tricky use of [] to make a
matrix from structure
    totalTrans = [];
    return;
end

toSubSrc = mne_read_source_spaces(toSubSrcSpaceFile);
fromSubSrc = mne_read_source_spaces(fromSubSrcSpaceFile);

[hi2hi{1},hi2hi{2}] = mne_read_morph_map(fromSub,toSub,freesurfDir);

%% Make lowrez transfer matrix

for iHemi = 1:2,

fromDecIdx = double(fromSubSrc(iHemi).vertno); %Decimation index for fromSub
toDecIdx = double(toSubSrc(iHemi).vertno); %Decimation index for toSub

nHiFrom = double(fromSubSrc(iHemi).np); %Num vertices hires in from subject
nHiTo = double(toSubSrc(iHemi).np); %Num vertices hires in "to" subject

nLoTo = double(toSubSrc(iHemi).nuse);
nLoFrom = double(fromSubSrc(iHemi).nuse);

idxFrom =
nearpoints(fromSubSrc(iHemi).rr',fromSubSrc(iHemi).rr(fromDecIdx,:)');
idxTo = nearpoints(toSubSrc(iHemi).rr',toSubSrc(iHemi).rr(toDecIdx,:)');

%Sparse mapping matrix from low res to hi res surface using nearest
%neighbour interpolation.
fromLo2Hi = sparse(1:nHiFrom,idxFrom,ones(size(idxFrom)),nHiFrom,nLoFrom);
toLo2Hi = sparse(1:nHiTo,idxTo,ones(size(idxTo)),nHiTo,nLoTo);

%Indexing matrix from Hi res "to" subject to low res "to" subject
toHi2Lo = sparse(1:nLoTo,toDecIdx,ones(nLoTo,1),nLoTo,nHiTo);

%Fix if missing some end columns because of a zero mapping
[i j s] = find(hi2hi{iHemi});
hi2hi{iHemi} = sparse(i,j,s,nHiTo,nHiFrom);

totalTrans{iHemi} = toHi2Lo*(hi2hi{iHemi}*fromLo2Hi);

end

%Make the hemi parts into 1 big matrix suitable for our default cortex.
mapMtx = [totalTrans{1}
sparse(size(totalTrans{1},1),size(totalTrans{2},2)); ...
    sparse(size(totalTrans{2},1),size(totalTrans{1},2)) totalTrans{2}];

Hi Pavan,

A small suggestion, if you in Linux, you can also called MNE function from
matlab, use some thing like this.

command = ['mne_make_movie --stcin <stcINfile> --subject <subj>
--morph fsaverage --smooth 5 --stc <stcOUTfile>'];

[st,wt] = unix(command);

if st ~=0
error('ERROR : error in generating morph stc file')
wt
end

In the pipeline in matlab you can use mixture of matlab and mne code.

Sheraz Khan, Ph.D.

Research Fellow,
Department of Neurology /
Martinos Center for Biomedical Imaging
Massachusetts General Hospital /
Harvard Medical School
149 13th Street, CNY-10.023
Boston, MA 02129
fax 617-948-5966

Hi Pavan, Hari, Justin, and others,

Hi Hari (and others),

Thanks very much for the suggestions and apologies for the lengthy
reply!

I did check that most rows are entirely zero (about 2300/2562).
Secondly,
I did look at the vertex numbering and they seem to make sense.
Also, I
visualized the result as stc and they do seem as sparse as the data
matrix
suggests.

It seems to me that smoothing (or as Matti recommends to call it in
Section 8.3 of the manual v2.6.1 : smudging/blurring) is the root of
the
difference between the mne_make_movie result and my MATLAB
implemenation.

The smoothing or spreading operator is, indeed, essential. Here is how
the morphing works in mne_make_movie:

1. A morphing map is composed. Using the aligned spherical surfaces (?
h.sphere.reg) of the two subjects, you get a linear interpolation
(sparse) matrix which interpolates the values at the vertices of a
triangle of a "source" surface to get the value at a vertex of the
"destination" surface.

2. The vertices of interest on the "destination" surface are
determined. This is done by finding the vertices nearest to the
vertices of a recursively subdivided icosahedron on the sphere.reg
surfaces. Only the rows corresponding to these vertices are retained
in the morphing map.

3. Since the interpolation (morphing) matrix not only assumes but
requires that data are available on all vertices of the source
surface, the spreading operator is applied. This is the one step not
currently implemented in Matlab.

4. The smeared data are multiplied by the (restricted) morphing matrix.

Here is the recipe we came up with Alex Gramfort during lunch to make
the spreading operator matrix. The C code uses direct iteration
instead of forming this matrix explicitly. By "valid neigboring
vertex" I mean a vertex which is adjacent to the vertex of interest
AND belongs to the source space OR has been defined during previous
iterations of the spreading operator.

This is what you need to do:

0. Initialize the spreading operator as an nvert x nvert identity
matrix (S_0)
1. Create an nvert x nvert empty matrix (S)
2. On each row, insert ones to the columns corresponding to the
vertices adjacent to the vertex corresponding to the row.
3. Zero columns corresponding to vertices that are not valid.
4. Divide each row by the number of non-zero entries on that row.
4. Compute the spreading operator corresponding to the k'th smoothstep
as S_k = S*S_{k-1}

I hope this helps.

- Matti

Dear Matti, Justin, Hari, Dan, Sheraz, and others!

Thanks very much for clarifying and proposing a solution. I guess Justin's
interpolation method (nearpoints) and Matti's iterative method are
slightly different in how sparse they allow the 'destination' surface to
be. I have tested the nearpoints function and it seems to work really well
& fast! I will also test Matti's suggestion and report back in case I find
something unexpected.

Thanks again,
Pavan

Hi,

if you feel like making this accessible to the MNE community
it would be great to have it directly in the MNE matlab toolbox

https://github.com/mne-tools/mne-matlab

It would be another mne_* function. I can give you access to the
code repository, or you can just send me the m files if you
prefer.

Alex

Hi Alex!

Sure and happy to do so. Attached is a self explanatory script called
mne_create_decimated_morph_map.m (along with the nearpoints function
scraped off the Vistasoft tools as recommended by Justin).

The script is a just slightly cleaned up version of what Justin had sent:
so I suspect the permissions etc. would rest with him. I will try to
implement Matti's interpolation method when I have some more time.

Thanks again,
Pavan

Hi,

if you feel like making this accessible to the MNE community
it would be great to have it directly in the MNE matlab toolbox

GitHub - mne-tools/mne-matlab: MNE scripting with Matlab

It would be another mne_* function. I can give you access to the
code repository, or you can just send me the m files if you
prefer.

Alex
--
Alexandre Gramfort, PhD
gramfort at nmr.mgh.harvard.edu
Dept. of Radiology MGH Martinos Center / Harvard Medical School
http://www-sop.inria.fr/members/Alexandre.Gramfort/

Dear Matti, Justin, Hari, Dan, Sheraz, and others!

Thanks very much for clarifying and proposing a solution. I guess
Justin's
interpolation method (nearpoints) and Matti's iterative method are
slightly different in how sparse they allow the 'destination' surface to
be. I have tested the nearpoints function and it seems to work really
well
& fast! I will also test Matti's suggestion and report back in case I
find
something unexpected.

Thanks again,
Pavan

Hi Pavan, Hari, Justin, and others,

Hi Hari (and others),

Thanks very much for the suggestions and apologies for the lengthy
reply!

I did check that most rows are entirely zero (about 2300/2562).
Secondly,
I did look at the vertex numbering and they seem to make sense.
Also, I
visualized the result as stc and they do seem as sparse as the data
matrix
suggests.

It seems to me that smoothing (or as Matti recommends to call it in
Section 8.3 of the manual v2.6.1 : smudging/blurring) is the root of
the
difference between the mne_make_movie result and my MATLAB
implemenation.

The smoothing or spreading operator is, indeed, essential. Here is how
the morphing works in mne_make_movie:

1. A morphing map is composed. Using the aligned spherical surfaces (?
h.sphere.reg) of the two subjects, you get a linear interpolation
(sparse) matrix which interpolates the values at the vertices of a
triangle of a "source" surface to get the value at a vertex of the
"destination" surface.

2. The vertices of interest on the "destination" surface are
determined. This is done by finding the vertices nearest to the
vertices of a recursively subdivided icosahedron on the sphere.reg
surfaces. Only the rows corresponding to these vertices are retained
in the morphing map.

3. Since the interpolation (morphing) matrix not only assumes but
requires that data are available on all vertices of the source
surface, the spreading operator is applied. This is the one step not
currently implemented in Matlab.

4. The smeared data are multiplied by the (restricted) morphing matrix.

Here is the recipe we came up with Alex Gramfort during lunch to make
the spreading operator matrix. The C code uses direct iteration
instead of forming this matrix explicitly. By "valid neigboring
vertex" I mean a vertex which is adjacent to the vertex of interest
AND belongs to the source space OR has been defined during previous
iterations of the spreading operator.

This is what you need to do:

0. Initialize the spreading operator as an nvert x nvert identity
matrix (S_0)
1. Create an nvert x nvert empty matrix (S)
2. On each row, insert ones to the columns corresponding to the
vertices adjacent to the vertex corresponding to the row.
3. Zero columns corresponding to vertices that are not valid.
4. Divide each row by the number of non-zero entries on that row.
4. Compute the spreading operator corresponding to the k'th smoothstep
as S_k = S*S_{k-1}

I hope this helps.

- Matti

---------

Matti Hamalainen, Ph.D.
Athinoula A. Martinos Center for Biomedical Imaging
Massachusetts General Hospital

msh at nmr.mgh.harvard.edu

The information in this e-mail is intended only for the person to whom
it
is
addressed. If you believe this e-mail was sent to you in error and the
e-mail
contains patient information, please contact the Partners Compliance
HelpLine at
MyComplianceReport.com: Compliance and Ethics Reporting . If the e-mail was sent to you
in
error
but does not contain patient information, please contact the sender and
properly
dispose of the e-mail.

_______________________________________________
Mne_analysis mailing list
Mne_analysis at nmr.mgh.harvard.edu
Mne_analysis Info Page

-------------- next part --------------
A non-text attachment was scrubbed...
Name: mne_create_decimated_morph_map.tar.gz
Type: application/postscript
Size: 45844 bytes
Desc: not available
Url : http://mail.nmr.mgh.harvard.edu/pipermail/mne_analysis/attachments/20110317/85166304/attachment.eps

Hi Pavan,

Sure and happy to do so. Attached is a self explanatory script called
mne_create_decimated_morph_map.m (along with the nearpoints function
scraped off the Vistasoft tools as recommended by Justin).

great. Could this code be BSD licensed? including the vistasoft code?

The script is a just slightly cleaned up version of what Justin had sent:
so I suspect the permissions etc. would rest with him. I will try to
implement Matti's interpolation method when I have some more time.

what would be great would be give an example script that works
with the MNE sample dataset and the freesurfer average brain (fsaverage)

I'll try to do this unless someone volunteers.

Alex

Hi Alex,

great. Could this code be BSD licensed? including the vistasoft code?

This from the Vistasoft homepage:

http://white.stanford.edu/software/
"Unless otherwise noted, all our code is released under the GPL. The
authors of the code retain all copyrights. By downloading our code from
the links below or from our CVS repository, you agree to be bound by the
terms of the GPL."

I'll try to do this unless someone volunteers.

Thanks!

Hi,

The GPL license may constitute a problem...

- Matti

The GPL license may constitute a problem...

yes it is a problem as it means this code cannot be distributed with MNE.
A solution consist in replacing the 1 nearest neighbor with KD-Tree with a
BSD licensed one found in the matlab file exchange e.g.

Alex

Hi,

we've come up with Matti to a pure matlab solution to the morphing step.
It might require some more testing but preliminary results are fine.

See:

https://github.com/mne-tools/mne-matlab/commit/5c77f11650c7d248f51fe457de6a5764d5dfcc06

Typically you would do something like :

from = 'sample';
to = 'fsaverage';
src = mne_read_source_spaces([getenv('SUBJECTS_DIR'),'/', from,
'/bem/sample-oct-6-src.fif']);
data{1} = randn(4098, 10);
data{2} = randn(4098, 10);
dmap_sel = mne_morph_data(from,src,to,data,5);

the output dmap_sel is data from 'sample' in the space of fsaverage.

feed back welcome.

Alex