NeuroElf wiki

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
====== Creating and combining masks from VTCs ======

===== Motivation =====
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

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

===== Algorithm =====
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 [[xff - MSK format|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

==== Code ====
<code matlab 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;</code>

===== Alternative  algorithm =====
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

==== Code ====
<code matlab 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;</code>