[Mne_analysis] regarding morphing in matlab
Justin Ales
justin.ales at gmail.com
Mon Mar 14 14:30:29 EDT 2011
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}];
On Mon, Mar 14, 2011 at 9:09 AM, Pavan Ramkumar <pavan at neuro.hut.fi> wrote:
> 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
>
> > 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
> >
> > On Mon, March 14, 2011 7:00 am, Pavan Ramkumar wrote:
> >> 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
> >>
> >>> Dear MNE users,
> >>>
> >>> I would like to morph my data from a single subject into a common
> >>> brain.
> >>> However, instead of writing the result out as an stc file (which takes
> >>> up
> >>> too much space) I would like to do the morphing from within matlab, and
> >>> then subsequently do some postprocessing on the morphed data before I
> >>> write out stc files.
> >>>
> >>> Does anybody know where to find the vertex numbers corresponding to the
> >>> --morphgrade parameter used in mne_make_movie? Specifically, for
> >>> morphgrade = 5, where do I find the indices of the 10242 vertices that
> >>> are
> >>> selectively stored into the stc file by mne_make_movie? As far as I
> >>> looked, the manual does not explicitly give this information.
> >>>
> >>> Please advice.
> >>>
> >>> Thanks in advance,
> >>> Pavan
> >>>
> >>
> >> _______________________________________________
> >> Mne_analysis mailing list
> >> Mne_analysis at nmr.mgh.harvard.edu
> >> https://mail.nmr.mgh.harvard.edu/mailman/listinfo/mne_analysis
> >>
> >>
> >>
> >
> >
> > --
> > Hari Bharadwaj
> >
> >
> > 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
> > http://www.partners.org/complianceline . 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
> https://mail.nmr.mgh.harvard.edu/mailman/listinfo/mne_analysis
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.nmr.mgh.harvard.edu/pipermail/mne_analysis/attachments/20110314/0a837edf/attachment.html
More information about the Mne_analysis
mailing list