% filenameList (e.g. as in SPM.xY.VY)
% voxelList in [x y z], 1 coordinate per row, untransformed

function extr = calculateImageData(filenameList, voxelList, imageOpts)
TRANSFORM_METHOD = 'image';

V = filenameList;

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

for kImage=1:nImage
    extr(kImage) = struct(...
        'dat',   repmat(NaN, [1 nVoxel]),...
        'mean',  repmat(NaN, [1 nVoxel]),...
        'sum',   repmat(NaN, [1 nVoxel]),...
        'nvx',   repmat(NaN, [1 nVoxel]),...
        'posmm', repmat(NaN, [3 nVoxel]),...
        'posvx', repmat(NaN, [3 nVoxel]));

    switch TRANSFORM_METHOD
        case 'single'
            roicenter = round(inv(V(kImage).mat)*[vox, ones(nVoxel,1)]');
            x = roicenter(1,:);
            y = roicenter(2,:);
            z = roicenter(3,:);
        case 'image'
            x = []; y = []; z = [];
            for iRoiFile = 1:nRoiFiles
                [x1 y1] = ndgrid(1:V(k).dim(1),1:V(k).dim(2));
                for p = 1:V(k).dim(3) % resample mask Vm(iRoiFile) in space of V(k)
                    B = spm_matrix([0 0 -p 0 0 0 1 1 1]);
                    M = inv(B*inv(V(k).mat)*Vm(iRoiFile).mat);
                    msk = find(spm_slice_vol(Vm(iRoiFile),M,V(k).dim(1:2),0));
                    if ~isempty(msk)
                        z1 = p*ones(size(msk(:)));
                        x = [x; x1(msk(:))];
                        y = [y; y1(msk(:))];
                        z = [z; z1];
                    end;
                end;
            end
    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