function data = get_voxel_values3d(V,NumSlices,interp)
%V   = spm_vol(spm_get(Inf,'*.IMAGE',['Select an Image to spit out']));
[x,y] = ndgrid(1:V.dim(1),1:V.dim(2));
for SliceNumber=1:NumSlices
  % Print progress information in command window
%   str   = sprintf('Slice %3d/%-3d ',SliceNumber,NumSlices);
%   fprintf('%-40s: %15s',str,' ')       
  % Get the data for this slice
  data(:,:,SliceNumber) = subfunc1(V,x,y,SliceNumber,interp);

end



function data = subfunc1(V,x,y,SliceNumber,interp)
% Create a handle to the image volume.
% V = spm_vol('filename.img');

% Define which voxels you want to get at.  e.g., plane 20...
%[x,y,z] = ndgrid(1:V.dim(1),1:V.dim(2),1:V.dim(3));

%z=repmat(1:V.dim(3)
z     = ones(size(x))*SliceNumber;

% sample the data...
data = spm_sample_vol(V,x,y,z,interp);% default interp = 0;