% filenameList as in SPM.xY.VY
% voxelList in [x y z], 1 coordinate per row

function extr = calculateImageData(filenameList, voxelList)
TRANSFORM_METHOD = 'single';
%function extr = calculateImageData(voxelList,des,smoothed)
% global USE_DRIVE_CHECK_HACK;

% dtype='PSTH';

% switch dtype 
%     case 'PSTH'
%         V=des.xY.VY;
%     case 'betas'
%         V=des.Vbeta;
% end;

% if USE_DRIVE_CHECK_HACK
% if V(1).fname(1)~='D'
%      for z=1:length(V) % Change Drive Letter - HACK!
%       	V(z).fname(1) = 'D';
%      end; 
% end
% end

% if (~smoothed)
%   for z=1:length(V) % Change smoothed Filename - HACK!
%       % D:....SUBJECTID\session\swfanders...
%       % D:....SUBJECTID\session\wfanders...
%       tmp = findstr(filesep,V(z).fname);
%       V(z).fname = strcat(V(z).fname(1:tmp(end)),V(z).fname(tmp(end)+2:end));
%   end;
% end

V = filenameList;

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

for kImage=1:nImage
% 	extr(kImage) = struct(...
%         'dat',  ...
%         'val',   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
    

	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