Dear FreeSurfer experts,
I am following the LME tutorial for the spatiotemporal mass-univariate model on
Matlab 8.2. When I get to the point where the covariance estimates are
segmented into homogeneous regions I run into a problem. Nothing really happens
in Matlab anymore.
This is my command:
[lhRgs,lhRgMeans] = lme_mass_RgGrow (lhsphere, lhRe, lhTh0, lhcortex, 2, 95);
And this is what the Matlab output looks like:
connected to 12 workers.
Starting region growing algorithm
Computing seeds ...
I've tried this already last week and nothing happened for a few days.
When I terminate the command this is the Matlab output:
Operation terminated by user during kron (after line 36)
In lme_mass_RgGrow>homogenity (line 299)
if (sum(sum(DistRgParams <= kron(ones(1,nv),thr),1) == np) >= prc*nv/100)
&&...
In lme_mass_RgGrow>merge (line 320)
if homogenity(Re,Params,auxRgs,Rkn,nst,prc)
In lme_mass_RgGrow>split (line 269)
[Rgs,mrg] = merge(Rgs,Re,Params,mrgRgs(i),mrgRgs(j),nst,prc);
In lme_mass_RgGrow>seeds (line 208)
[Rgs,spltf(i),nnRgs] = split(Rgs,Re,Params,coord,splRgs(i),nst,prc);
In lme_mass_RgGrow>vtxwGrowing (line 77)
[Rgs1,Rgseed,Rgstd]=
seeds(Re(:,maskvtx),Params(:,maskvtx),coord(:,maskvtx),nst,prc);
In lme_mass_RgGrow (line 62)
Regions = vtxwGrowing(Re,Theta,sphcoord,AdjM,maskvtx,nst,prc);
I have then copied the function and added a few "fprintf" (see attachment) to
get an idea as to where matlab has problems and the output looks like this:
Starting region growing algorithm
Computing seeds ...
point5
point6
point7
point8
point7
point8
point7
point8
point7
point8
point7
point8
point7
point8
point7
point6
point7
point8
point7
point8
point7
point8
point7
point8
point7
point8
point7
[...]
If I understand that correctly, matlab gets hung up on functions 'function tf =
homogenity(Re,Params,Rgs,Rn,nst,prc)' and 'function [Rgs,tf] =
merge(Rgs,Re,Params,Rn1,Rn2,nst,prc)'. Does anybody have an idea what that
means and how it can be fixed?
How long is this step supposed to take? I have 81 subjects with 1 - 3 time
points.
Any help or suggestions would be much appreciated!
Cheers
Clara
--
Clara Kühn, Phd Candidate
Max-Planck-Institute for Human Cognitive and Brain Science
Department of Neuropsychology
Stephanstrasse 1A
04103 Leipzig, Germany
Phone: +49 341 - 9940 2271
Fax: +49 341 - 9940 2260
Web: www.cbs.mpg.de
E-Mail: cku...@cbs.mpg.defunction [Regions,RgMeans] = lme_mass_RgGrow(SphSurf,Re,Theta,maskvtx,nst,prc)
% [Regions,RgMeans] = lme_mass_RgGrow(SphSurf,Re,Theta,maskvtx,nst,prc)
%
% This function implements a region growing algorithm along the spherical
% surface to find homogeneous regions comprising locations with similar
% covariance components.
%
% Input
% SphSurf: Spherical surface, which represent a spherical coordinate system
% for the surface where the analysis is to be done. It is a structure in
% fs_read_surf format. SphSurf.tri = t x 3 matrix of triangle indices,
% 1-based, t=#triangles and SphSurf.coord = 3 x nv matrix of coordinates,
% nv=#vertices.
% Re: Matrix of residual errors at each location (nmxnv, nm total # of
% maps, nv #vertices.
% Theta: Matrix whose colums are estimators of the covariance components at.
% each location.
% maskvtx: Mask's vertices (1-based). Default [] (all vertices included).
% nst: Number of standard deviations of Theta inside a region for it to
% be considered as homogeneous. Default 2.
% prc: Percent of vertices inside an homogeneous region that are allowed to
% not meet the homogenity criteria. This parameter makes the algorithm
% robust to noise and outliers. Default 95%.
%
% Output
% Regions: 1 x nv segmentation vector containing a region number assigned
% to each vertice along the surface.
% RgMeans: Matrix assigning the mean of Theta within a region to each vertex
% belonging to that region.
%
% $Revision: 1.1.2.2 $ $Date: 2013/02/23 21:08:10 $
% Original Author: Jorge Luis Bernal Rusiel
% CVS Revision Info:
% $Author: nicks $
% $Date: 2013/02/23 21:08:10 $
% $Revision: 1.1.2.2 $
% References:
% Gonzalez, R.C., Woods, R.E. (2002). Digital Image Processing. 2nd Edition,
% New Jersey: Prentice Hall, Inc.
%
tic;
if nargin < 3
error('Too few inputs');
elseif nargin < 6
prc = 95;
if nargin < 5
nst = 2;
if nargin < 4
maskvtx = [];
end
end;
end;
nv = size(Theta,2);
if isempty(maskvtx)
maskvtx = 1:nv;
end;
AdjM = AdjMtx(SphSurf,maskvtx);
display(' ');
display('Starting region growing algorithm');
[phi,theta] = sphCoord(SphSurf.coord');
sphcoord = [phi,theta]';
Regions = vtxwGrowing(Re,Theta,sphcoord,AdjM,maskvtx,nst,prc);
[RgMeans,nRg] = lme_mass_RgMean(Regions,Theta);
display(' ');
display([num2str(nRg-1) ' homogeneous regions were obtained.']);
et = toc;
display(['Elapsed time is ' num2str(et/60) ' minutes.']);
end
function Rgs = vtxwGrowing(Re,Params,coord,Adj,maskvtx,nst,prc)
display('Computing seeds ...');
[Rgs1,Rgseed,Rgstd]= seeds(Re(:,maskvtx),Params(:,maskvtx),coord(:,maskvtx),nst,prc);
fprintf(2,'point1\n');
Rgseed = maskvtx(Rgseed);
ns = length(Rgseed);
display([num2str(ns) ' seeds were computed after splitting the surface']);
[np,nv] = size(Params);
Rgs = zeros(1,nv,'uint32');
Rgs(maskvtx) = Rgs1;
notIncl = true(1,nv);
notIncl(Rgs == 0) = false;
notIncl(Rgseed) = false;
maxRgsz = 90; %Do not form regions with size greater than maxRgsz
Rgvtxs = zeros(ns,maxRgsz,'uint32');
Rgvtxs(:,1) = Rgseed;
Rgvtxs_prvtxind = zeros(1,ns,'uint32');
Rgvtxs_lastind = ones(1,ns,'uint32');
Rgind = 1:ns;
thr = nst * Rgstd;
growing_ns = ns;
while (growing_ns > 0)
fprintf(2,'point2 while schleife\n');
i = 1;
display(' ');
display([num2str(growing_ns) ' seeds for growing']);
display(['Current maximum region size ' num2str(max(Rgvtxs_lastind)) ' vertices']);
while (i <= growing_ns)
fprintf(2,'point3 while schleife 2\n');
Rgvtxs_prvtxind(Rgind(i)) = Rgvtxs_prvtxind(Rgind(i)) + 1;
prvtx = Rgvtxs(Rgind(i),Rgvtxs_prvtxind(Rgind(i)));
RgParams = Params(:,Rgvtxs(Rgind(i),1:Rgvtxs_lastind(Rgind(i))));
mRgParams = mean(RgParams,2);
Rgsz = Rgvtxs_lastind(Rgind(i));
%processing the current vertex prvtx
adjvtxs = Adj(prvtx,Adj(prvtx,:) > 0);
%try to include only the adjacent vertices not included in other
%regions or in this region itself
adjvtxs = adjvtxs(notIncl(adjvtxs));
%try to include first the adjacent vertices with the param closest
%to the region mean param.
[~,loc] = min(sum(abs(Params(:,adjvtxs) - kron(ones(1,length(adjvtxs)),mRgParams)),1));
NewRgParams = [RgParams Params(:,adjvtxs(loc))];
mNewRgParams = mean(NewRgParams,2);
nRgnv = size(NewRgParams,2);
DistNewRgParams = abs(NewRgParams-kron(ones(1,nRgnv),mNewRgParams));
%an adjacent vertex is included in the region if the difference
%between its params and the region mean params is bellow the
%threshold thr or if the region still contains a percent
%prc of vertices bellow the threshold.
nadjv = length(adjvtxs); j = 1;
while (Rgsz < maxRgsz) && (j <= nadjv) && ...
sum(sum(DistNewRgParams <= kron(ones(1,nRgnv),thr(:,Rgind(i))),1) == np) >= prc*(nRgnv)/100
%In addition the vertex's residuals must be correlated above 0.5
%with any other vertex's residual in the region.
fprintf(2,'point3 while schleife 3\n');
if min(min(corrcoef([Re(:,Rgvtxs(Rgind(i),1:Rgvtxs_lastind(Rgind(i)))),Re(:,adjvtxs(loc))]))) >= 0.5
Rgvtxs_lastind(Rgind(i)) = Rgvtxs_lastind(Rgind(i)) + 1;
Rgsz = Rgsz + 1;
Rgvtxs(Rgind(i),Rgvtxs_lastind(Rgind(i))) = adjvtxs(loc);
notIncl(adjvtxs(loc)) = false;
RgParams = Params(:,Rgvtxs(Rgind(i),1:Rgvtxs_lastind(Rgind(i))));
mRgParams = mean(RgParams,2);
end;
adjvtxs = [adjvtxs(1:loc-1),adjvtxs(loc+1:end)];
[~,loc] = min(sum(abs(Params(:,adjvtxs) - kron(ones(1,length(adjvtxs)),mRgParams)),1));
NewRgParams = [RgParams Params(:,adjvtxs(loc))];
mNewRgParams = mean(NewRgParams,2);
nRgnv = size(NewRgParams,2);
DistNewRgParams = abs(NewRgParams-kron(ones(1,nRgnv),mNewRgParams));
j = j + 1;
end;
if (Rgvtxs_prvtxind(Rgind(i)) == Rgvtxs_lastind(Rgind(i))) || (Rgvtxs_lastind(Rgind(i)) == maxRgsz)
Rgind = [Rgind(1:i-1) Rgind(i+1:end)];
growing_ns = growing_ns - 1;
else
i = i + 1;
end;
end;
end;
for i=1:ns
Rgs(Rgvtxs(i,1:Rgvtxs_lastind(i))) = i;
end;
nni = sum(notIncl);
notInclv = find(notIncl == 1);
Rgs(notInclv) = ns+1:ns+nni;
%try to add holes to their most similar adjacent homogeneous region
display(' ');
display([num2str(nni) ' holes (unassigned vertices). Trying to form new'...
' regions among them or add them to their most similar adjacent'...
' region ...']);
fprintf(2,'point4\n');
for i=1:nni
hadjvtxs = Adj(notInclv(i),Adj(notInclv(i),:) > 0);
if ~isempty(hadjvtxs)
closestdistRParams = sum(abs(Params(:,notInclv(i)) - mean(Params(:,Rgs == Rgs(hadjvtxs(1))),2)));
nadjvtxs = length(hadjvtxs);
m = 1;
j = 2;
while j <= nadjvtxs
distRParams = sum(abs(Params(:,notInclv(i)) - mean(Params(:,Rgs == Rgs(hadjvtxs(j))),2)));
if distRParams < closestdistRParams
closestdistRParams = distRParams;
m = j;
end;
j = j+1;
end;
[Rgs,mrg] = merge(Rgs,Re,Params,Rgs(notInclv(i)),Rgs(hadjvtxs(m)),nst,prc);
if mrg
notIncl(notInclv(i)) = 0;
end;
end;
end;
nni = sum(notIncl);
if nni > 0
display([num2str(nni) ' holes (unassigned vertices) remained. Will be'...
' treated as independent regions of size 1.']);
else
display(' No holes (unassigned vertices) remained. ');
end;
end
function [Rgs,Rgseed,Rgstd] = seeds(Re,Params,coord,nst,prc)
[np,nv] = size(Params);
Rgs = ones(1,nv,'uint32');
spl = true;
fprintf(2,'point5\n');
while spl
splRgs = unique(Rgs);
ntospl = length(splRgs);
spltf = zeros(1,ntospl);
nRgs = ntospl;
for i=1:ntospl
[Rgs,spltf(i),nnRgs] = split(Rgs,Re,Params,coord,splRgs(i),nst,prc);
nRgs = nRgs + nnRgs;
end;
if sum(spltf) == 0
spl = false;
end;
end;
splRgs = unique(Rgs);
nRgs = length(splRgs);
Rgseed = zeros(1,nRgs);
Rgstd = zeros(np,nRgs);
for i=1:nRgs
Rgv = find(Rgs == splRgs(i));
RgParams = Params(:,Rgv);
Rgstd(:,i) = std(RgParams,0,2);
nv = length(Rgv);
spp = sum(coord(:,Rgv),2)/nv;
[~,loc] = min(sum(abs(coord(:,Rgv) - kron(ones(1,nv),spp)),1));
Rgseed(i) = Rgv(loc);
end;
end
function [Rgs,tf,nRnv] = split(Rgs,Re,Params,coord,Rnum,nst,prc)
fprintf(2,'point6\n');
if homogenity(Re,Params,Rgs,Rnum,nst,prc)
nRnv = 0;
tf = false;
else
%split the non-homogeneous region
Rvts = find(Rgs == Rnum);
nv = length(Rvts);
spp = sum(coord(:,Rvts),2)/nv;
nR = max(Rgs);
for i=1:nv
if (coord(1,Rvts(i)) <= spp(1)) && (coord(2,Rvts(i)) <= spp(2))
Rgs(Rvts(i)) = 1;
elseif (coord(1,Rvts(i)) <= spp(1)) && (coord(2,Rvts(i)) > spp(2))
Rgs(Rvts(i)) = 2;
elseif (coord(1,Rvts(i)) > spp(1)) && (coord(2,Rvts(i)) <= spp(2))
Rgs(Rvts(i)) = 3;
elseif (coord(1,Rvts(i)) > spp(1)) && (coord(2,Rvts(i)) > spp(2))
Rgs(Rvts(i)) = 4;
end;
end
Rnv = unique(Rgs(Rvts));
nRnv = length(Rnv) - 1;
for i=1:nRnv
Rgs(Rvts(Rgs(Rvts) == Rnv(i+1))) = nR+i;
end;
Rgs(Rvts(Rgs(Rvts) == Rnv(1))) = Rnum;
mrgRgs = [Rnum nR+1:nR+nRnv];
%merge homogeneous subregions
nmtomrg = length(mrgRgs);
i = 1;
while i <= nmtomrg
j = i + 1; mrg = false;
while (~mrg) && (j <= nmtomrg)
[Rgs,mrg] = merge(Rgs,Re,Params,mrgRgs(i),mrgRgs(j),nst,prc);
j = j + 1;
end;
if mrg
mrgRgs = [mrgRgs(1:j-2) mrgRgs(j:end)];
nmtomrg = nmtomrg - 1;
else
i = i + 1;
end;
end;
tf = true;
end;
end
function tf = homogenity(Re,Params,Rgs,Rn,nst,prc)
tf = false;
maxRgsz = 90;
fprintf(2,'point7\n');
%Regions with size greater than maxRgsz are not considered homogeneous
Rgv = find(Rgs == Rn);
if length(Rgv) <= maxRgsz
RgParams = Params(:,Rgv);
mRgParams = mean(RgParams,2);
stRgParams = std(RgParams,0,2);
thr = nst*stRgParams;
[np,nv] = size(RgParams);
DistRgParams = abs(RgParams-kron(ones(1,nv),mRgParams));
if (sum(sum(DistRgParams <= kron(ones(1,nv),thr),1) == np) >= prc*nv/100) &&...
min(min(corrcoef(Re(:,Rgv)))) >= 0.5
tf = true;
end;
end;
end
function [Rgs,tf] = merge(Rgs,Re,Params,Rn1,Rn2,nst,prc)
fprintf(2,'point8\n');
if Rn1 <= Rn2
Rkn = Rn1;
Ren = Rn2;
else
Rkn = Rn2;
Ren = Rn1;
end;
auxRgs = Rgs;
auxRgs(auxRgs == Ren) = Rkn;
if homogenity(Re,Params,auxRgs,Rkn,nst,prc)
Rgs = auxRgs;
tf = true;
else
tf = false;
end;
end
_______________________________________________
Freesurfer mailing list
Freesurfer@nmr.mgh.harvard.edu
https://mail.nmr.mgh.harvard.edu/mailman/listinfo/freesurfer
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.
