function [extr voxelcount] = calculateImageData(filenameList, coordlist,gFkt)

%center coordinates
vox = [];
nCoords = size(coordlist,1);
for iCoord = 1:nCoords
   vox = [vox ; coordlist(iCoord).coord];
end

%maximum extracted voxels over all images
voxelcount = 0;

V = filenameList;
nImage = numel(V);
nVoxel = nCoords;
for kImage=1:nImage
    roicenter = round(inv(V(kImage).mat)*[vox, ones(nVoxel,1)]');

    subvoxelcount = 0;
    for iVoxel = 1:nVoxel

        %radius
        rad = coordlist(iVoxel).rad;

        tmp = spm_imatrix(V(kImage).mat);
        vdim = tmp(7:9);
        vxrad = ceil((rad*ones(nVoxel,3))./(ones(nVoxel,1)*vdim))';
        
        if rad==0
            x = roicenter(1,iVoxel);
            y = roicenter(2,iVoxel);
            z = roicenter(3,iVoxel);
        else
            [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(:));
%             [x y z]
        end;

        dat = spm_sample_vol(V(kImage), x, y, z,0);
 
        %%Implement Spatial grouping here
        
        dat = gFkt(dat);
        
        % no grouping
        for iSubVoxel = 1:size(dat,1)
            subvoxelcount = subvoxelcount +1;
            extr(kImage).dat(subvoxelcount)    = dat(iSubVoxel);
        end
        
    end
    voxelcount = max(voxelcount,subvoxelcount);
end   

end