Creating and combining masks from VTCs

I wanted to create this page for two reasons:

• demonstrate quickly the ease of access on xff objects (in this case BrainVoyager QX types)
• give users the ability to create an SPM-comparable mask when running GLMs

You should have all the VTCs you want to combine in a group analysis (RFX-GLM) at hand. In this example, I assume that they are listed in an MDM file, which is BrainVoyager's way of specifying data for a group analysis.

The algorithm is as follows:

• make settings (i.e. required intensity level, either as ratio or fixed value)
• load MDM file
• load the first (next) VTC file
• for first VTC, create a MSK (mask) object with the same spatial layout
• for subsequent VTC, check spatial layout
• compute the mean over time
• clear the VTC
• build mask for this VTC (threshold mean)
• if first VTC, copy mask to combined mask array
• combine mask with combined mask array
• continue with step 3 (next VTC)
• optionally perform required step for mask thresholding (if not simple binary combination)
• store combined mask and save under new filename
• clear mask object

create_msk_from_vtcs.m

% settings
% - intensity threshold: values < 10 are treated as relative threshold!
ithresh = 1.5;
 
% - combine as sum or mask
summask = true;
 
% - if combine as sum, threshold ?
subthresh = 0.75;
 
 
% load MDM file
mdm = xff('*.mdm', 'Please specify the MDM file to create a mask for...');
 
% loop over VTCS
for vc = 1:size(mdm.XTC_RTC, 1)
 
    % load first VTC
    vtc = xff(mdm.XTC_RTC{vc, 1});
    if ~isxff(vtc, 'vtc')
        clearxffobjects({vtc});
        error('Not a VTC file!');
    end
 
    % create mask
    if vc == 1
        msk = xff('new:msk');
        msk.Resolution = vtc.Resolution;
        msk.XStart = vtc.XStart;
        msk.XEnd = vtc.XEnd;
        msk.YStart = vtc.YStart;
        msk.YEnd = vtc.YEnd;
        msk.ZStart = vtc.ZStart;
        msk.ZEnd = vtc.ZEnd;
 
    % test layout
    else
        if msk.Resolution ~= vtc.Resolution || ...
           msk.XStart ~= vtc.XStart || ...
           msk.XEnd ~= vtc.XEnd || ...
           msk.YStart ~= vtc.YStart || ...
           msk.YEnd ~= vtc.YEnd || ...
           msk.ZStart ~= vtc.ZStart || ...
           msk.ZEnd ~= vtc.ZEnd
            vtc.ClearObject;
            msk.ClearObject;
            error(sprintf('VTC %d (%s) mismatches layout!', vc, mdm.XTC_RTC{vc, 1}));
        end
    end
 
    % compute the mean over time
    mvtc = squeeze(mean(vtc.VTCData));
 
    % clear VTC
    vtc.ClearObject;
 
    % threshold
    if ithresh < 10
        mask = (mvtc >= (ithresh .* mean(mvtc(:))));
    else
        mask = (mvtc >= ithresh);
    end
 
    % store or combine
    if vc == 1
        combined = mask;
    else
 
        % sum or and mask?
        if summask
            combined = combined + mask;
        else
            combined = combined & mask;
        end
    end
end
 
% combine across subjects required (summed mask)
if summask
 
    % subject threshold
    combined = (combined >= (subthresh * size(mdm.XTC_RTC, 1)));
end
 
% set in mask
msk.Mask = uint8(combined);
 
% save mask
msk.SaveAs;
 
% clear msk and mdm
msk.ClearObject;
mdm.ClearObject;

In case the MSK files/objects have already been created (e.g. as gray matter mask per subject), this is a way to average them. This will only work if the masks are in the same space, e.g. large TAL box!

• making settings (i.e. percentage threshold for combined mask)
• locating all mask files to be averaged
• loading all masks and adding the data
• thresholding the data
• saving new mask

combine_existing_msks.m

% settings (80% of masks must have a set voxel)
gthresh = 0.8;
 
% locate masks
%  - this could be enhanced by changing the pattern, e.g. using
%    'SK*_TAL*.msk' to locate only masks of subjects beginning with SK in TAL space
%  - additionally, the startfolder and depth could be altered
mskfiles = findfiles(pwd, '*.msk');
% alternative: mskfiles = findfiles([pwd '/SK*/VTC*/RUN*'], 'SK*TAL*.msk', 'depth=1');
 
% loop over masks
msk = [];
for mc = 1:numel(mskfiles)
 
    % clear old mask
    if ~isempty(msk)
        msk.ClearObject;
    end
 
    % load mask
    msk = xff(mskfiles{mc});
 
    % for first mask
    if mc == 1
 
        % copy data
        mask = uint16(msk.Mask);
 
    % otherwise
    else
 
        % add data
        mask = mask + uint16(msk.Mask);
    end
end
 
% threshold mask
mask = uint8(mask >= uint16(ceil(gthresh * numel(mskfiles))));
 
% store
msk.Mask = mask;
 
% save
msk.SaveAs;
 
% clear object
msk.ClearObject;