function extr = calculateImageData(filenameList, voxelList, rad)

V = filenameList;

vox = voxelList;
nVoxel = size(vox,1);
nImage = numel(V);

% rad = 0;

for kImage=1:nImage
    roicenter = round(inv(V(kImage).mat)*[vox, ones(nVoxel,1)]');
    x = roicenter(1,:);
    y = roicenter(2,:);
    z = roicenter(3,:);


    for iVoxel = 1:nVoxel

        if rad==0
            x = roicenter(1,iVoxel);
            y = roicenter(2,iVoxel);
            z = roicenter(3,iVoxel);
        else
            tmp = spm_imatrix(V(kImage).mat);
            vdim = tmp(7:9);
            vxrad = ceil((rad*ones(1,3))./(ones(nVoxel,1)*vdim))';
            [x y z] = ndgrid(-vxrad(1,iVoxel):sign(vdim(1)):vxrad(1,iVoxel), ...
                -vxrad(2,iVoxel):sign(vdim(2)):vxrad(2,iVoxel), ...
                -vxrad(3,iVoxel):sign(vdim(3)):vxrad(3,iVoxel));
            sel = (x./vxrad(1,iVoxel)).^2 + (y./vxrad(2,iVoxel)).^2 + ...
                (z./vxrad(3,iVoxel)).^2 <= 1;
            x = roicenter(1,iVoxel)+x(sel(:));
            y = roicenter(2,iVoxel)+y(sel(:));
            z = roicenter(3,iVoxel)+z(sel(:));
        end;


        dat = spm_sample_vol(V(kImage), x, y, z,0);
        extr(kImage).dat(iVoxel)      = dat;
        extr(kImage).mean(iVoxel)     = nanmean(dat);
        extr(kImage).nvx(iVoxel)      = numel(dat);
    end;

end;
end