my_voxblock

 returns voxel block and modified mat file for pointlist

 $Id$

0001 function vblock = my_voxblock(pts, mat, vals)
0002 % returns voxel block and modified mat file for pointlist
0003 %
0004 % $Id$
0005 
0006 if isempty(vals)
0007   vals = ones(1, size(pts,2));
0008 end
0009 
0010 % dimensions of containing block
0011 st = min(pts, [], 2)';
0012 fn = max(pts, [], 2)';
0013 dsz = fn-st+1;
0014 
0015 % reset mat according to voxel block
0016 newM = mat * spm_matrix(st-1);
0017 
0018 % set voxels in block
0019 dat1 = zeros(dsz);
0020 dpts = pts - (st'-1) * ones(1, size(pts,2)); 
0021 dinds = dpts(1,:) + ...
0022     (dpts(2,:)-1) * dsz(1) + ...
0023     (dpts(3,:)-1) * dsz(1)*dsz(2);
0024 dat1(dinds) = vals;
0025 
0026 % add a shell of one voxel to isolate from interpolation 0s
0027 dat = zeros(dsz+2);
0028 dat(2:end-1,2:end-1,2:end-1) = dat1;
0029 mat = newM * spm_matrix([-1 -1 -1]);
0030 
0031 vblock = struct('dat', dat, 'mat', mat);