Creating Masks and Applying them to Neuroimaging Data in Matlab

Post date: May 24, 2009 1:42:46 AM

This tutorial demonstrates how to implement simple procedures to mask images. This is like creating an ROI in a specific region and is useful especially in instances where you want to rip something out of an image (e.g., rip out part of an Atlas as your ROI) or when you want to rip out all voxels whose values are above a certain threshold.

%=====================================================

% Example 1: Ripping out BA 10 from a Brodmann Area Atlas

% Step 1: Load in the Atlas

BrodmannAtlas = spm_read_vols(spm_vol('/Users/mvlombardo/Documents/BrainAtlases/brodmann.img')); % Reads in the brodmann atlas into a variable called BrodmannAtlas

% Step 2: Make a mask that extracts BA 10. I know that the Brodmann atlas has its voxels labeled based on the Brodmann Area number. I simply have to extract all voxels whose value is 10

BA10 = ismember(BrodmannAtlas,10);

% ismember is a crucial function. Give it a matrix in the first argument and tell it the value you are looking for in that matrix within the second argument.

% Note that ismember takes BrodmannAtlas, assigns a logical value of "true" (or 1) to voxels whose value is 10, and assigns all other voxels a logical value of false (or 0).

% Step 3: Write out BA 10 as its own image.

% Load the header information of another file of similar dimensions and voxel sizes

V = spm_vol('/Users/mvlombardo/Documents/BrainAtlases/brodmann.img');

% Change the name in the header

V.fname = 'BA10.nii';

V.private.dat.fname = V.fname;

% Write out the new header and data

spm_write_vol(V,BA10);

%=====================================================

% Example 2: Extracting all voxels whose value is above a specified t-threshold

% Step 1: Load spmT*.img which is a map containing t-values from a statistical contrast

TMap = spm_read_vols(spm_vol('/Users/mvlombardo/Documents/fMRI/Contrast1/spmT_0001.img'));

tthresh = 3.1768423; % This is the corresponding t-threshold I'll need to extract all voxels whose value is greater than this.

% Step 2: Extract all voxels in TMap that are greater than tthresh

TMap_suprathresh = TMap>=tthresh;

% Note that TMap_suprathresh is a matrix of the same size as TMap filled with 0 or 1, where voxels greater than tthresh are 1, otherwise they are 0

% Step 3: Write out TMap_suprathresh as its own image.

V = spm_vol('/Users/mvlombardo/Documents/fMRI/Contrast1/spmT_0001.img');

V.fname = 'tmap_suprathresh.nii';

V.private.dat.fname = V.fname;

spm_write_vol(V,TMap_suprathresh);

%=====================================================

% Example 3: Making a conjunction map between two suprathreshold t-maps

% Read in both t-maps

TMap1 = spm_read_vols(spm_vol('/Users/mvlombardo/Documents/fMRI/Contrast1/spmT_0001.img'));

TMap2 = spm_read_vols(spm_vol('/Users/mvlombardo/Documents/fMRI/Contrast1/spmT_0002.img'));

% Define t-threshold

tthresh = 3.1768423;

% Mask out suprathreshold voxels from both t-maps

TMap1_suprathresh = TMap1>=tthresh;

TMap2_suprathresh = TMap2>=tthresh;

% Make conjunction map as the logical AND of both suprathreshold t-maps

ConjunctionMap = TMap1_suprathresh & TMap2_suprathresh;

% Write out the new conjunction map

V = spm_vol('/Users/mvlombardo/Documents/fMRI/Contrast1/spmT_0001.img');

V.fname = 'ConjunctionMap_TMap1_AND_TMap2.nii';

V.private.dat.fname = V.fname;

spm_write_vol(V,ConjunctionMap);

%=====================================================