spm_vb_geom

 Create GMRF spatial precision matrix with geometric boundary conditions 
 FORMAT [D] = spm_vb_geom (vxyz)

 vxyz      list of neighbours of voxels to be analysed
 D         spatial precision matrix

 See equation in Woolrich et al. (ref 26 in paper VB2) 
 also described in dicussion in paper VB2

 %W% Will Penny %E%

0001 function [D] = spm_vb_geom (vxyz)
0002 % Create GMRF spatial precision matrix with geometric boundary conditions
0003 % FORMAT [D] = spm_vb_geom (vxyz)
0004 %
0005 % vxyz      list of neighbours of voxels to be analysed
0006 % D         spatial precision matrix
0007 %
0008 % See equation in Woolrich et al. (ref 26 in paper VB2)
0009 % also described in dicussion in paper VB2
0010 %
0011 % %W% Will Penny %E%
0012 
0013 voxels=size(vxyz,1);
0014 D=speye(voxels);
0015 
0016 for j=1:voxels,
0017     num_neighbors(j)=length(nonzeros(vxyz(j,:)));
0018 end
0019 
0020 for j=1:voxels,
0021     neighbors=nonzeros(vxyz(j,:));
0022     geom_means=sqrt(num_neighbors(j)*num_neighbors(neighbors));
0023     D(j,neighbors)=-1./geom_means;
0024 end