spm_vb_neighbors

 Create list of neighbors of voxels to be analysed
 FORMAT [vxyz] = spm_vb_neighbors (xyz)

 xyz      [Nvoxels x 3] list of voxel positions which are to be analysed

 vxyz     [Nvoxels x 4] list of neighbouring voxels
          vxyz(j,:)=[N1 N2 N3 0] means that there are only 3 neighbors
          of voxel j, and their numbers (ie. where they appear in the xyz list) 
          are N1, N2 and N3

 %W% Will Penny and Nelson Trujillo-Barreto %E%

0001 function [vxyz] = spm_vb_neighbors (xyz)
0002 % Create list of neighbors of voxels to be analysed
0003 % FORMAT [vxyz] = spm_vb_neighbors (xyz)
0004 %
0005 % xyz      [Nvoxels x 3] list of voxel positions which are to be analysed
0006 %
0007 % vxyz     [Nvoxels x 4] list of neighbouring voxels
0008 %          vxyz(j,:)=[N1 N2 N3 0] means that there are only 3 neighbors
0009 %          of voxel j, and their numbers (ie. where they appear in the xyz list)
0010 %          are N1, N2 and N3
0011 %
0012 % %W% Will Penny and Nelson Trujillo-Barreto %E%
0013 
0014 [voxels,Ndims]=size(xyz);
0015 % if voxels<Ndims
0016 %     disp('Error in neighbours: xyz wrong way round');
0017 %     return
0018 % end
0019 
0020 x=kron(xyz(:,1),ones(voxels,1));
0021 xx=repmat(xyz(:,1),voxels,1);
0022 yz=nonzeros(abs(x-xx));
0023 hx=min(yz);
0024 
0025 y=kron(xyz(:,2),ones(voxels,1));
0026 yy=repmat(xyz(:,2),voxels,1);
0027 xz=nonzeros(abs(y-yy));
0028 hy=min(yz);
0029 
0030 vxyz=zeros(voxels,4);
0031 
0032 for j=1:voxels,
0033 
0034     pxyz=repmat(xyz(j,:),voxels,1);
0035     dist=sqrt(sum((pxyz-xyz).^2,2));
0036     vx=find((dist<=hx)&(dist~=0));
0037     vy=find((dist<=hy)&(dist~=0));
0038     vxyzt=union(vx,vy);
0039     vxyz(j,1:length(vxyzt))=vxyzt';
0040     
0041 end;
0042 
0043