function coordinateStruct = computeCoordinates(imageStruct)
%imageStruct.(someName).images
%imageStruct.(someName).roiImages
coordinateStruct = struct;

subjectCellList = fieldnames(imageStruct);

for subject = 1:length(subjectCellList)
    subjectFieldName = cell2mat(subjectCellList(subject));
    % ROI Image Coordinates
    V  = imageStruct.(subjectFieldName).images;
    Vm = imageStruct.(subjectFieldName).roiImages;
    coordinateStruct.(subjectFieldName) = ...
        computeRoiImageCoordinates(V,Vm);
    % Parsed Voxel Definitions
    coordinateStruct.(subjectFieldName) = ...
        appendManualVoxelCoordinates(coordinateStruct.(subjectFieldName),V);
    
end
end

function coordinates = appendManualVoxelCoordinates(coordinateStruct,V)
    coordinates =coordinateStruct;
end

function coordinates = computeRoiImageCoordinates(V,Vm)
    nVolImage = length(V);
    nRoiImage = length(Vm);
    
    coordinates = cell(nVolImage,1);
 
    for iVolImage = 1:nVolImage
        coordinates{iVolImage} = [];
        for jRoiImage = 1:nRoiImage
            x = []; y = []; z = [];
            [x1 y1] = ndgrid(1:V(iVolImage).dim(1),1:V(iVolImage).dim(2));
            for p = 1:V(iVolImage).dim(3) % resample mask Vm(jRoiImage) in space of V(iVolImage)
                B = spm_matrix([0 0 -p 0 0 0 1 1 1]);
                M = inv(B*inv(V(iVolImage).mat)*Vm(jRoiImage).mat);
                msk = find(spm_slice_vol(Vm(jRoiImage),M,V(iVolImage).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;
            coordinates{iVolImage} = [coordinates{iVolImage}; x, y, z];
        end
    end
end