Dear Torben (and Kevin who also reported this bug recently),
This is a problem introduced shortly before the SPM8 release due to
our attempt to remove the dependency on optimization toolbox. Putting
the two attached files in external/fieldtrip/private should solve the
problem. Please let me know if it doesn't. The files will be included
in the next SPM8 update.
Best,
Vladimir
On Thu, Apr 30, 2009 at 9:41 AM, Torben Ellegaard Lund
<[log in to unmask]> wrote:
> Dear List
>
> I have been analysing the MMN example dataset as described in the spm
> manual. When I get to the forward modelling I get the following error after
> having selected “EEG-BEM” in response to the “Which EEG head model?”
> question:
>
>
>
>
> using the pre-specified triangulated boundaries
> determining source compartment (1)
> determining skin compartment (3)
> not using the isolated source approach
> ??? Attempted to access pvec(:,2); index out of bounds because
> size(pvec)=[0,5].
>
> Error in ==> fitsphere at 64
> C = -0.5*pvec(2:4) / pvec(1);
>
> Error in ==> triangle4pt at 62
> [center,radius] = fitsphere(sph_pnt);
>
> Error in ==> prepare_bemmodel at 167
> vol = triangle4pt(vol);
>
> Error in ==> ft_prepare_bemmodel at 11
> [varargout{1:nargout}] = funhandle(varargin{:});
>
> Error in ==> spm_eeg_inv_forward at 77
> vol = ft_prepare_bemmodel(cfg, vol);
>
> Error in ==> spm_eeg_inv_forward_ui at 46
> D = spm_eeg_inv_forward(D);
>
> Error in ==> spm_eeg_inv_imag_api>Forward_Callback at 87
> handles.D = spm_eeg_inv_forward_ui(handles.D);
>
> Error in ==> spm_eeg_inv_imag_api at 53
> feval(varargin{:}); % FEVAL switchyard
>
> ??? Error while evaluating uicontrol Callback
>
>>>
>>>
>
>
>
> Does anybody have any suggestions?
>
>
> Best
> Torben
>
>
>
> Torben Ellegaard Lund
> Assistant Professor, PhD
> The Danish National Research Foundation's Center of Functionally Integrative
> Neuroscience (CFIN)
> Aarhus University
> Aarhus University Hospital
> Building 10G, 5th floor
> Noerrebrogade 44
> 8000 Aarhus C
> Denmark
> Phone: +4589494380
> Fax: +4589494400
> http://www.cfin.au.dk
> [log in to unmask]
>
function vol = triangle4pt(vol)
% FORMAT vol = triangle4pt(vol)
%
% Takes the volume model and estimates the 4th point of each triangle of
% each mesh.
%
% In each vol.bnd sub-structure, a field '.pnt4' is added. The '.pnt4'
% field is a Ntri*3 matrix, with the coordinates of a point for each
% triangle in the meshed surface.
%
% Explanations:
% The point is that for some BEM, specifically 'solid angle', calculation
% it is necessary to estimate the local curvature of the true surface which
% is approximated by the flat triangle. One way to proceed is to use
% "close by" vertices to estimate the overall area's curvature.
% A more elegant(?) way uses a 4th point for each triangle: the "centroid"
% of the triangle is simply pusehd away from the triangle surface to fix
% the local surface curvature (assuming the surface is smooth enough).
% This 4th point is thus hovering above/under the triangle and can be used
% to fit a sphere on the triangle in a realistic way.
%
% Method:
% - The 4th point can/could be defined at the tessalation stage, based on
% the anatomical images directly.
% - With any model, the curvature can be estimated/approximated by looking
% at the vertices around the triangle considered and fit a sphere on
% those few vertices, assuming the surface is smooth enough
% The latter option is the one followed here.
% The extra-vertices considered here are those 3 which are linked to the
% triangle by 2 edges.
%__________________________________________________________________________
%
% written by Christophe Phillips, 2009/01/19
% Cyclotron Research Centre, University of li?ge, belgium
% $Log: triangle4pt.m,v $
% Revision 1.4 2009/04/30 16:59:49 vlalit
% Bug fix for the problem of too flat mesh surfaces as suggested by Christophe
%
% Revision 1.3 2009/04/01 13:26:05 roboos
% use Taubin's method of sphere fitting (fitsphere) instead of the iterative implementation by Guido Nolte
%
Ns = length(vol.bnd);
for ii=1:Ns % treat each mesh one at a time
tri = vol.bnd(ii).tri;
pnt = vol.bnd(ii).pnt;
Nt = size(tri,1);
pnt4 = zeros(Nt,3);
for jj=1:Nt % treat each triangle on a t a time
lt1 = find(tri(:,1)==tri(jj,1) | tri(:,2)==tri(jj,1) | ...
tri(:,3)==tri(jj,1));
lt2 = find(tri(:,1)==tri(jj,2) | tri(:,2)==tri(jj,2) | ...
tri(:,3)==tri(jj,2));
lt3 = find(tri(:,1)==tri(jj,3) | tri(:,2)==tri(jj,3) | ...
tri(:,3)==tri(jj,3));
lt = [intersect(lt1,lt2); intersect(lt2,lt3); intersect(lt3,lt1)];
lt(lt==jj) = [];
% list of 3 directly surrounding triangles
lv = tri(lt,:);
lv = setxor(lv(:),tri(jj,:));
% list of 3 voxels connected by 2 edges to the jj_th triangle.
sph_pnt = pnt([tri(jj,:) lv],:);
[center,radius] = fitsphere(sph_pnt);
% best fitting sphere radius & centre, for the 6 points chosen
pnt_c = sum(pnt(tri(jj,:),:))/3;
% centroid of the triangle treated
if isfinite(radius)
tmp = pnt_c-center;
vd = tmp/norm(tmp);
% unit vector from sphere center in direction of middle triangle
pnt4(jj,:) = center+vd*radius ;
% projection of the centroid, at 'radius' from the center
else
% radius at infinite, i.e. flat surface
pnt4(jj,:) = pnt_c;
end
end
vol.bnd(ii).pnt4 = pnt4;
end
return
% figure, plot3(pnt(:,1),pnt(:,2),pnt(:,3),'.')
% axis equal, hold on
% plot3(pnt(tri(jj,:),1),pnt(tri(jj,:),2),pnt(tri(jj,:),3),'*')
% plot3(pnt(lv,1),pnt(lv,2),pnt(lv,3),'o')
% plot3(pnt4(jj,1),pnt4(jj,2),pnt4(jj,3),'rs')
function [C,R] = fitsphere(pnt)
% FITSPHERE fits the centre and radius of a sphere to a set of points
% using Taubin's method.
%
% Use as
% [center,radius] = fitsphere(pnt)
% where
% pnt = Nx3 matrix with the Carthesian coordinates of the surface points
% and
% center = the center of the fitted sphere
% radius = the radius of the fitted sphere
% Copyright (C) 2009, Jean Daunizeau (for SPM)
%
% $Log: fitsphere.m,v $
% Revision 1.3 2009/04/30 16:59:36 vlalit
% Bug fix for the problem of too flat mesh surfaces as suggested by Christophe
%
% Revision 1.2 2009/04/01 12:37:12 roboos
% updated documentation
%
% Revision 1.1 2009/04/01 12:16:14 roboos
% new function, see email from Guillaume Flandin from 23 March
%
x = pnt(:,1);
y = pnt(:,2);
z = pnt(:,3);
% Make sugary one and zero vectors
l = ones(length(x),1);
O = zeros(length(x),1);
% Make design mx
D = [(x.*x + y.*y + z.*z) x y z l];
Dx = [2*x l O O O];
Dy = [2*y O l O O];
Dz = [2*z O O l O];
% Create scatter matrices
M = D'*D;
N = Dx'*Dx + Dy'*Dy + Dz'*Dz;
% Extract eigensystem
[v, evalues] = eig(M);
evalues = diag(evalues);
Mrank = sum(evalues > eps*5*norm(M));
if (Mrank == 5)
% Full rank -- min ev corresponds to solution
% Minverse = v'*diag(1./evalues)*v;
[v,evalues] = eig(inv(M)*N);
[dmin,dminindex] = max(diag(evalues));
pvec = v(:,dminindex(1))';
else
% Rank deficient -- just extract nullspace of M
pvec = null(M)';
[m,n] = size(pvec);
if m > 1
pvec = pvec(1,:);
end
end
if isempty(pvec)
C = [NaN NaN NaN];
R = Inf;
else
% Convert to (R,C)
C = -0.5*pvec(2:4) / pvec(1);
R = sqrt(sum(C*C') - pvec(5)/pvec(1));
end
% if nargout == 1,
% if pvec(1) < 0
% pvec = -pvec;
% end
% C = pvec;
% else
% C = -0.5*pvec(2:4) / pvec(1);
% R = sqrt(sum(C*C') - pvec(5)/pvec(1));
% end
|