extract_voxel_values.m

% Extracts values from an image data file based on regions/clusters of
% voxels in mask marsbar or image file(s) or voxel.
%
% Usage: [Ym R info] = extract_voxel_values(mask,data)
%
% Ym    – matrix of output mean of values across voxels per image (row) per
%         region (column).
% R     – structure containing voxel info per region (R), per image (I):
%         R.I.Ya  – output vector of values of all voxels.
%         R.I.xyz – matrix indices of voxels in mask region/cluster.
%         R.I.mni – mni coordinates of voxels in mask region/cluster.
%         E.g. to get vector of value for region r, image i, after 
%         extraction, type R(r).I(i).Ya
% info  – structure containing region and image files used:
%         info.regions – string array of region files used.
%         info.images – string array of image files used.
% mask  – string array of mask image file name(s).
% data  – string array of image data file name(s). Can be any analyze or
%         nifti image that contains voxel values of interest, which can be
%         beta values, contrast values, t values etc.
%         e.g. [‘beta_0005.img’;’beta_0001.img’];
%
% If mask or data not defined, GUI asks to select files.
%
% Example usage:
% extract_voxel_values(‘L_Fusiform_roi.mat’,’beta_0005.img’)
%
% Dependencies: Marsbar (if *_roi.mat ROI mask files used).

% Modified from extract_voxel_values_old.m by Josh Goh 24 May 2013.
 
function [Ym R info] = extract_voxel_values(mask,data)
 
% Check input
if nargin<1
    mask = spm_select(Inf,’any’,’Select mask ROI files’,[],pwd);
    data = spm_select(Inf,’image’,’Select data file (*.img or *.nii)’,[],pwd);
end
 
info.regions = mask;
info.images  = data;
 
% Loop image
for i = 1:size(data,1)
    
    % Read image header
    V = spm_vol(deblank(data(i,:)));
    
    % Loop regions
    for r = 1:size(mask,1)
        
        % Get voxels from mask
        [~,~,e] = fileparts(deblank(mask(r,:)));
        switch e
            case {‘.mat’}
                roi = maroi(deblank(mask(r,:)));
                xyz = voxpts(roi,deblank(data(i,:)));
            case {‘.nii’,’.img’}
                maskdata = spm_read_vols(spm_vol(deblank(mask(r,:))));
                [x,y,z] = ind2sub(size(maskdata),find(maskdata));
                xyz = [x y z]’;
        end
        
        % Extract data
        Ya = spm_sample_vol(V,xyz(1,:),xyz(2,:),xyz(3,:),0);
        R(r).I(i).Ya = Ya;
        
        % Compute mean across voxels
        Ym(i,r) = mean(Ya(~isnan(Ya)));
        
        % Compute MNI coordinates
        R(r).I(i).mni = vox2mni(V.mat,xyz);
        R(r).I(i).xyz = xyz;
        
    end
end