Matlab Tutorials

This page is for those interested in little Matlab tidbits I use frequently that might be of use to others? It assumes you already know a bit about Matlab (e.g., you know about strings, cell arrays, etc). In these little tutorials, I usually put comments up to describe each step of the code. Comments are in green text and usually start with the %. The actual code is in this font.

Subscribe to posts

Compute Effect Sizes

posted Feb 27, 2012 12:21 PM by Michael Lombardo [ updated Apr 22, 2012 12:53 PM ]

Recently I came across Marco del Giudice's paper on using Mahalanobis distance as a multivariate effect size, and comparing it to univariate effect sizes like Cohen's d. This prompted me to write up something similar for Matlab that will take multivariate data from 2 groups and compute both Cohen's d for univariate effects and Mahalanobis D for the multivariate effect size. In addition to this, the attached script has the option of implementing optimal shrinkage towards the diagonal sample covariance (e.g., Ledoit and Wolf) necessary for computing Mahalanobis D and will use bootstrapping for hypothesis testing on both univariate and multivariate effect sizes. For this to work in Matlab, you'll need the Statistics Toolbox as there is heavy reliance on various functions (e.g., nanmean, nancov).

| Attachments: ComputeEffectSizesBootstrap.m

Creating Masks and Applying them to Neuroimaging Data in Matlab

posted May 23, 2009 6:42 PM by Michael Lombardo [ updated Aug 12, 2009 3:01 PM ]

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);

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

Using Matlab to Create Unix Shell Scripts

posted May 21, 2009 5:24 PM by Michael Lombardo [ updated Aug 12, 2009 3:02 PM ]

One use of Matlab I find highly useful is in creating unix shell scripts. Since I'm more fluent in Matlab than in shell scripting, I essentially use Matlab to format and write my shell scripts, which can then be used in the unix shell. This is very helpful when you want to create a batching script on the go and don't have time to really code it out in a unix shell script, when it could be done in less time by just writing it in Matlab code.

In this example, I'm going to make a shell script that essentially batches my creation text files that are the mean timeseries of deep white matter seeds. To extract the timeseries from each seed I'm going to call on the fslmeants script which is an FSL utility that I find handy for ripping out the mean timeseries of seed regions of interest.

% Step 1. Make a file identifier that points to an open .sh script that you'll be writing your script to. You need this in order to start writing text into the shell script you are creating.

fid = fopen('CreateWhiteMatterTS.sh','wt');

% Step 2. Loop over a small set of subjects. While looping, create a string that will be called in unix to rip out the timeseries for each participant

sublist = {'12001','12002','12004'}; % My list of subject IDs to loop over

WhiteMatterSeedFilename = 'WhiteMatterSeed.nii'; % My white matter seed filename

% Now I loop over the three subjects, each time making a new input filename, output filename, and unix string with the correct subject ID

for isub = 1:length(sublist)

infile = sprintf('%s_RestingState.nii.gz',sublist{isub}); % This is the input file name for fslmeants to use

outfile = sprintf('%s_WhiteMatterTS.txt',sublist{isub}); % This is the output file name for fslmeants to use

unixstr = sprintf('fslmeants -i %s -o %s -m %s;',infile,outfile,WhiteMatterSeedFilename); % This is what I'm executing in the shell script.

fprintf(fid,'%s\n',unixstr); % This adds each line to the shell script

end % end my for loop

fclose(fid); % now close the shell script so I can't write anymore into it

% Step 3. Make sure the shell script has the proper read, write, and execute permissions and then run it from the Matlab prompt.

% Notice I put an ! before each unix command. This is how you execute unix commands in the unix shell while in Matlab.

!chmod u+rwx CreateWhiteMatterTS.sh % This will give the user read, write, and execute permissions on the shell script.

!./CreateWhiteMatterTS.sh % This will pass this string to the unix shell and execute the shell script

Reading and Writing Neuroimaging Data from Matlab

posted May 21, 2009 4:57 PM by Michael Lombardo [ updated Apr 19, 2010 7:21 AM ]

*** Note that text in green starting with a % is a comment in Matlab and is here to explain each step of the code in this font

Reading data into the Matlab workspace

% Step 1. Define filepath and filename of the data I want to read into the Matlab workspace

datadir = '/data/fMRI';

fname = '12001.nii';

Data2Read = fullfile(datadir,fname);

% What I've done in Step 1 is to construct a filename such as /data/fMRI/12001.nii from its filepath (the variable called datadir) and filename (the variable called fname). Both these variables are string variables and I've used the function fullfile to essentially concatenate these two strings together. The end result is a string variable called Data2Read, which is essentially a string that looks like this: /data/fMRI/12001.nii

% Step 2. Now pass Data2Read into the function spm_vol in order to read in the header information for the file. Note that the spm_vol function comes with the SPM toolbox in Matlab that is a standard software package for neuroimagers. Having SPM is essential for this part.

HeaderInfo = spm_vol(Data2Read);

% In Step 2, the header information for the data you want to read in is now in your workspace. This is required for Step 3, where you will use the header information to actually read in the entire data.

% Step 3. Now use spm_read_vols to read in the data

Data = spm_read_vols(HeaderInfo);

Your data is now loaded into the Matlab workspace.

Writing data out of Matlab into a NIFTI file

% Step 1. Take the header information from a previous file with similar dimensions and voxel sizes and change the filename in the header.

HeaderInfo.fname = 'NewData.nii'; % This is where you fill in the new filename

HeaderInfo.private.dat.fname = HeaderInfo.fname; % This just replaces the old filename in another location within the header.

% Step 2. Now use spm_write_vol to write out the new data. You need to give spm_write_vol the new header information and corresponding data matrix

spm_write_vol(HeaderInfo,Data); % where HeaderInfo is your header information for the new file, and Data is the image matrix corresponding to the image you'll be writing out.

1-4 of 4