I recently sent out a code to do VBM without the user interface. Though this code seemed to work
for brains without lesions, it didn't work when a lesion mask was used. I have since edited it and
now it should work universally. I am also including the function that do_VBM calls to coregister
the pre and post scans. Feedback is appreciated. I also have a quick question. Does anyone
have a code to crop some slices off the bottom of an image? Some of my pre/post scans go to
different levels in the cerebellum. I want to make them the same so my difference data is not as
confounded. Also, if anyone has codes to do statistics on VBM data, I would very much appreciate
them. Thanks.

function do_VBM(subj_folder)
%
%Performs VBM analysis in SPM5. It outputs files c1*, wc1* for grey matter,
%c2* and wc2* for white matter, and c3* for CSF. You can then run stats on
%wc1* using SPM or SnPM.
%
%INPUT: subj_folder - name of subject folder containing Images_Pre and
% Images_Post
%

%-----------***Set Defaults and Paths***-----------------

templatename= '/Users/neurolabII/matlab/codes_for_imaging/spm5/canonical/avg152T1.nii';
root='/Volumes/CITherapy$/Neuroimaging_Backup/Neurolab_II/matlab/codes_for_imaging/
spm5/tpm/';
filename_pre=fullfile(subj_folder, 'Images_Pre/T1/Analyze/image_1.img');
filename_post=fullfile(subj_folder, 'Images_Post/T1/Analyze/rimage_1.img');

spm_defaults
global defaults
defaults.analyze.flip=0;
%make voxels 1mm^3
%defaults.normalise.write.vox = [1 1 1];
%rmfield(defaults.normalize.write, 'bb');


%----------------***Coregister the pre scan to the post scan***--------------------%

pre_post_coregister(subj_folder, 1);


%------------*******Segment the pre scan******---------------


fprintf('Segmenting the Pre scan \n')
%define the inputs to spm_preproc
opts.tpm = char(fullfile(root,'grey.nii'),...
               fullfile(root,'white.nii'),...
               fullfile(root,'csf.nii')); %use the default segmented templates
opts.ngaus = [3 2 2 5]; %number of Gaussians. Should I change this to [2 2 2 4]?
opts.warpreg = 1; %if normalized images appear distorted, increase this to 10
opts.warpco = 25; %decrease for more detailed deformations, but this will increase run time
opts.biasreg = 0.01; %medium bias regularization
opts.biasfwhm = 75; %higher for images with less bias, keep for now...
opts.regtype = 'mni'; %when using the mni templates .... otherwise select 'subj'?
opts.fudge = 5; %what does this do?
opts.samp = 3; %smaller values sample more of the data, but make the script run slower
mask=fullfile(subj_folder, 'Images_Pre/T1/Analyze/mlimage_1.img');
opts.msk = mask; %mask the lesions

%define the writing options
opts2.biascor =1; %write bias corrected images
% opts2.GM =[1 1 1]; %write modulated, unmodulated normalized, and native GM
% opts2.WM =[1 1 1]; %write modulated, unmodulated normalized, and native WM
opts2.GM =[0 0 1];
opts2.WM =[0 0 1];
opts2.CSF =[0 0 1]; %write CSF in native space only

%get the image information
V=spm_vol(filename_pre);

% get segmentation information
results = spm_preproc(V, opts);

% convert segmentation information to sn-files
[po,pin] = spm_prep2sn(results);

% write segmented data
spm_preproc_write(po,opts2);

%save the segmentation information
[pth fn ext] = fileparts(filename_pre);
save(fullfile(pth, [fn '_seg_sn.mat']), '-struct', 'po');
save(fullfile(pth, [fn '_seg_inv_sn.mat']), '-struct', 'pin');

%define the images to normalize
grey=fullfile(subj_folder, 'Images_Pre/T1/Analyze/c1image_1.img');
white=fullfile(subj_folder, 'Images_Pre/T1/Analyze/c2image_1.img');
Vgrey=spm_vol(grey);
Vwhite=spm_vol(white);

%write unmodulated normalized
spm_write_sn(Vgrey,po);
spm_write_sn(Vwhite,po);

%------------*******Segment the post scan******---------------

fprintf('Segmenting the Post scan \n')

%change the options, so only native segmented images are written
% opts2.GM =[0 0 1];
% opts2.WM =[0 0 1];

%get the image information
V2=spm_vol(filename_post);

% get segmentation information
results = spm_preproc(V2, opts);

% convert segmentation information to sn-files
[po2,pin2] = spm_prep2sn(results);

% write segmented data
spm_preproc_write(po2,opts2);

%save the segmentation information
[pth2 fn2 ext2] = fileparts(filename_post);
save(fullfile(pth2, [fn2 '_seg_sn.mat']), '-struct', 'po2');
save(fullfile(pth2, [fn2 '_seg_inv_sn.mat']), '-struct', 'pin2');

%------***Write the normalized and modulated normalized files***------
% for the post scan according to *_seg_sn.mat
% files from the pre scan

%replace pre scan with post scan in po to normalize pre scan with post scan
%parameters
% po.VF=po2.VF;
%
% %change options to only write normalized modulated and unmodulated
% opts2.GM =[1 1 0];
% opts2.WM =[1 1 0];
% opts2.CSF=[0 0 0]; %don't write CSF
%
% %write the images
% spm_preproc_write(po,opts2);


%define the images to normalize
grey=fullfile(subj_folder, 'Images_Post/T1/Analyze/c1rimage_1.img');
white=fullfile(subj_folder, 'Images_Post/T1/Analyze/c2rimage_1.img');
Vgrey=spm_vol(grey);
Vwhite=spm_vol(white);

%write unmodulated normalized
spm_write_sn(Vgrey,po);
spm_write_sn(Vwhite,po);

%write modulated normalized
%defaults.normalise.write.preserve = 1;
%spm_write_sn(Vgrey,po,defaults.normalise.write);



%------------------------***Smooth the Images***---------------------

%smooth the normalized grey matter images using an 10mm Gaussian kernel
%create name strings for the images that need to be smoothed
smoothed_grey_post=fullfile(pth2, ['sm_wc1' fn2 ext2]);
smoothed_grey_pre=fullfile(pth, ['sm_wc1' fn ext]);
pre_grey = fullfile(pth, ['wc1' fn ext]);
post_grey = fullfile(pth2, ['wc1' fn2 ext2]);
%smooth pre
Vpre_grey=spm_vol(pre_grey);
spm_smooth(Vpre_grey, smoothed_grey_pre, 10);
%smooth post
Vpost_grey=spm_vol(post_grey);
spm_smooth(Vpost_grey, smoothed_grey_post, 10);


This function calls:

function pre_post_coregister(subj_folder,which_scan)
%
%Coregisters the post scan to the pre scan (6 parameter affine transformation)
%
%INPUTS:
% subj_folder - full name of subject folder, including path
% which_scan - 1 for pre (default), 2 for post
%
%OUTPUTS:
% Resliced image saved as rimage_1.img

%set the default
if nargin<2
    which_scan=1;
end

if which_scan==1
    templatename = fullfile(subj_folder, 'Images_Pre/T1/Analyze/image_1.img');
    filename=fullfile(subj_folder, 'Images_Post/T1/Analyze/image_1.img');
else
    templatename = fullfile(subj_folder, 'Images_Post/T1/Analyze/image_1.img');
    filename=fullfile(subj_folder, 'Images_Pre/T1/Analyze/image_1.img');
end

spm_defaults
global defaults
%make highest quality, lowest speed
defaults.realign.estimate.quality = 1;
%post scan gets registered to pre scan, rather than mean of two images
defaults.realign.estimate.rtm=0;
%do not write mean image
defaults.realign.estimate.mean=0;
%smaller separation (more points sampled) yields more accurate results
defaults.realign.estimate.sep=2;
defaults.realign.estimate.which=1;
rmfield(defaults.realign.estimate, 'weight');
%do not flip image
defaults.analyze.flip=0;
%make voxels 1mm^3?
defaults.normalise.write.vox = [1 1 1];
defaults.normalise.estimate.wtsrc = 1;
%affine registration only?
defaults.normalise.estimate.nits=0;
%makes sure that all parts of the brain are weighted equally in the
%normalization
defaults.normalise.write.bb = [[-91 -125 -71];[89 91 109]];
%dnrm.estimate.nits = 0; %set to zero for linear only

%save the original nonrealigned image as oimage_1.img
[pathstr,name]=fileparts(filename);
% orig_filename=fullfile(pathstr,'oimage_1.img');
% copyfile(filename, orig_filename);
%make a copy of the .hdr file as well
orig_hdr=fullfile(pathstr,'oimage_1.hdr');
orig_img=fullfile(pathstr,'oimage_1.img');
hdr=fullfile(pathstr,'image_1.hdr');
img=fullfile(pathstr,'image_1.img');
copyfile(hdr,orig_hdr);
copyfile(img,orig_img);

%realigns post scan to pre scan
spm_realign({templatename; filename},defaults.realign.estimate)
spm_reslice({templatename; filename},defaults.realign.estimate) %when using a mask this step is
%not neccessary, but if you use a mask you get a "Dimensions do not match error"