Hi
Someone pointed out a bug in the SPM5 routine I emailed last month for using
the FieldMap undistortion toolbox on multiple files - I've attached a
corrected version; save and run "cspm_FieldMap" at the command prompt.
You'll need the FieldMap.m file I posted in Nov 2006 (for SPM5), the
remainder of the FieldMap toolbox (see SPM extensions), and the fieldmap
(phase magnitude images coregistered to the EPI/volumes to unwarp, and two
TE's used to derive the fieldmap).
Best wishes,
Paul
function Pout = cspm_FieldMap( Pin,FM_phase,FM_mag,TE,type,overwrite)
% CSPM_FIELDMAP_NGUI - Cobled together from Fieldmap_ngui
% Applies fieldmap to images using the fieldmap toolbox. Works for SPM5
% using the modifications to "Fieldmap.m" from Paul Macey (posted to spm
% mailing list Nov 2006). This version only works with phase/mag, but could
% probably he hacked.
%
% CSPM_FIELDMAP_NGUI will prompt for files, TE's, type
%
% CSPM_FIELDMAP_NGUI( P, FM_phase, FM_mag, TE )
% P = images to correct
% FM_phase = phase image filename
% FM_mag = magnitude image filename
% TE = 2 x 1 vector [TE_first_echo TE_second_echo]
%
% Phase image will be scaled to between -pi and pi
%
% CSPM_FIELDMAP_NGUI( P, FM_phase, FM_mag, TE, type, overwrite )
% type - if "dti" the DTI parameters will be applied to the corrected images
% using the Diffusion toolbox
% overwrite - 0 / 1
%
% Pout = CSPM_FIELDMAP_NGUI( P, FM_phase, FM_mag, TE, type, overwrite )
% Returns corrected filenames
%
% This assumes pm_defaults.m set to phase magnitude
% pm_def.INPUT_DATA_FORMAT = 'PM'; % 'RI' = load two real and
% % imaginary image pairs
% % 'PM' = load one or two
% % phase and magnitude image
% % pairs.
if nargin < 1
Pin = spm_select(Inf,'image','Select files to unwarp');
end
if nargin < 2
FM_phase = spm_select(Inf,'image','Select phase image');
disp('If the fieldmaps were not collected at the same slice locations as the ')
disp('images to unwarp, the phase and magnitude images must be coregistered ')
disp('with the images to unwarp (before running the current routine). ')
disp('Use the magnitude image to coregister, and apply the parameters to the phase image.')
end
if nargin < 3
FM_mag = spm_select(Inf,'image','Select magnitude image');
end
if nargin < 4
disp(['If you do not have the TE''s but know the TE difference, you can enter',...
'any numbers that give the same difference. For example, if you know the',...
'TE difference is 3.9 ms you can enter "0 3.9". (Philips give the TE ',...
' difference in the PAR file.)'])
TE = spm_input('Enter two echo-times (TE''s in ms)',1,'r','0 0',2);
end
if nargin < 5
type = spm_input('Which type of scans?',2,'b',{'fMRI (and other)' 'DTI'},{'fmri' 'dti'},1);
end
if strcmpi(type,'dti')
if ~check_diffusion_toolbox
error('Cannot apply DTI parameters to unwarped files without diffusion toolbox.')
end
end
if nargin < 6
overwrite = 0;
end
if ~iscell(Pin)
for i = 1:size(Pin,1)
P{i} = deblank(Pin(i,:)); %#ok<AGROW>
end
Pin = P;
clear P
end
if ~check_fieldmap_toolbox
disp('Cannot continue without FieldMap toolbox; see http://www.fil.ion.ucl.ac.uk/spm/ext/.')
disp('For SPM5, also need modifications to "Fieldmap.m" from Paul Macey ')
disp('posted to spm mailing list Nov 2006.')
error('Cancelled.')
end
FM_phase = cspm_fieldmap_scalephase( FM_phase, overwrite );
% Determine unwarped filenames
for i = 1:length(Pin)
[pth,nm,ext] = fileparts(Pin{i});
Pout{i} = [pth,filesep,'u',nm,ext]; %#ok<AGROW>
end
if ~overwrite && cspm_filesallexist(Pout)
return
end
% Template script for FieldMap toolbox
% _____________________________________________________________________
%
% This script gives an example of how to run the FieldMap toolbox
% FieldMap.m without using the GUI. It can be expanded using standard
% matlab code to for example create multiple field maps or use a single
% fieldmap to unwarp multiple images.
%
% As it stands the script uses routines from the FieldMap toolbox to
% create a single field map which is matched to an EPI and then used
% to unwarp it. A structural image is loaded and matched to the unwarped
% EPI.
%
% For details about the FieldMap toolbox, see FieldMap.man. For a
% description of the components of the structure IP, see FieldMap.m.
% For an introduction to the theoretcial and practical principles behind
% the toolbox, see principles.man.
%_________________________________________________________________
% FieldMap_ngui.m Chloe Hutton & Jesper Andersson 27/01/04
%
%----------------------------------------------------------------------
% Set up default parameters and structures
%----------------------------------------------------------------------
%spm_defaults
IP = FieldMap('Initialise'); % Gets default params from pm_defaults
IP.uflags.iformat = 'PM';
IP.et{1} = min(TE);
IP.et{2} = max(TE);
%----------------------------------------------------------------------
% Load measured field map data - phase and magnitude or real and imaginary
%----------------------------------------------------------------------
IP.P{1} = spm_vol(FM_phase);
IP.P{2} = spm_vol(FM_mag);
% if IP.uflags.iformat=='PM'
% for index=1:4
% IP.P{index} = FieldMap('LoadFilePM',index);
% end
% else
% for index=1:4
% IP.P{index} = FieldMap('LoadFile',index);
% end
% end
%----------------------------------------------------------------------
% Or you may want to load a precalculated Hz phase map instead...
%----------------------------------------------------------------------
% [IP.fm, IP.pP] = FieldMap('LoadFieldMap');
%----------------------------------------------------------------------
% Create field map (in Hz) - this routine calls the unwrapping
%----------------------------------------------------------------------
IP.fm = FieldMap('CreateFieldMap',IP);
%----------------------------------------------------------------------
% Write out field map
% Outputs -> fpm_NAME-OF-FIRST-INPUT-IMAGE.img
%----------------------------------------------------------------------
FieldMap('Write',IP.P{1},IP.fm.fpm,'fpm_',64,'Smoothed phase map');
%----------------------------------------------------------------------
% Convert Hz to voxels and write voxel displacement map
% Outputs -> vdm_NAME-OF-FIRST-INPUT-IMAGE.img
%----------------------------------------------------------------------
[IP.vdm, IP.vdmP]=FieldMap('FM2VDM',IP);
%----------------------------------------------------------------------
% Match voxel displacement map to image
% Outputs -> mag_NAME-OF-FIRST-INPUT-IMAGE.img
%----------------------------------------------------------------------
IP.epiP = spm_vol(Pin{1});%FieldMap('LoadEPI');
IP.vdmP = FieldMap('MatchVDM',IP);
for i = 1:length(Pin)
%----------------------------------------------------------------------
% Select an EPI to unwarp
%----------------------------------------------------------------------
IP.epiP = spm_vol(Pin{i});%FieldMap('LoadEPI');
%----------------------------------------------------------------------
% Unwarp EPI
%----------------------------------------------------------------------
IP.uepiP = FieldMap('UnwarpEPI',IP);
%----------------------------------------------------------------------
% Write unwarped EPI
% Outputs -> uNAME-OF-EPI.img
%----------------------------------------------------------------------
IP.uepiP = FieldMap('Write',IP.epiP,IP.uepiP.dat,'u',IP.epiP.dt(1),'Unwarped image');
end
% If DTI, apply direction/gradient parameters
if strcmpi(type,'dti')
% Apply gradient correction
jobs{1}.tools{1}.vgtbx_Diffusion{1}.dti_reorient_gradient.srcimgs = Pin;
jobs{1}.tools{1}.vgtbx_Diffusion{1}.dti_reorient_gradient.tgtimgs = Pout;
jobs{1}.tools{1}.vgtbx_Diffusion{1}.dti_reorient_gradient.snparams = {''};
spm_jobman('run',jobs)
clear jobs
end
if nargout == 0
clear Pout
end
% ------------------------------------------------------------------------------
function Pout = cspm_fieldmap_scalephase( P, overwrite )
% Script to scale a phase map so that max = pi and min =-pi radians.
% Writes out scaled image prepended with 'sc'.
if nargin < 2
overwrite = 1;
end
cellout = 0;
if nargin == 0
P = spm_select(1, 'IMAGE','Select fieldmap phase image to scale');
elseif iscell(P)
cellout = 1;
P = char(P);
end
V=spm_vol(P);
vol=spm_read_vols(V);
mn=min(vol(:));
mx=max(vol(:));
svol=-pi+(vol-mn)*2*pi/(mx-mn);
oV=V;
name='sc';
oV.fname=[spm_str_manip(V.fname, 'h') ['/' name] deblank(spm_str_manip(V.fname,'t'))];
if overwrite || ~exist(oV.fname,'file')
oV.dt(1)=spm_type('float32');
oV = spm_create_vol(oV);
spm_write_vol(oV,svol);
end
if nargout > 0
Pout = oV.fname;
if cellout
Pout = {Pout};
end
end
% ------------------------------------------------------------------------------
function ok = cspm_filesallexist( P, prefix )
if iscell(P)
for i = 1:length(P)
if nargin < 2
fn = P{i};
else
[pth,nm,xt] = fileparts(P{i});
fn = [pth,filesep,prefix,nm,xt];
end
if ~exist(fn,'file')
ok = 0;
return
end
end
ok = 1;
end
% ------------------------------------------------------------------------------
function ok = check_fieldmap_toolbox
% Check that fieldmap toolbox is on path
ok = 0;
if ~exist(which('FieldMap'),'file')
tb = spm('TBs');
I = strmatch('FieldMap',{tb(:).name});
if isempty(I) || I == 0
msgbox('FieldMap toolbox not installed!! Go to SPM website.',...
'MRI Analysis','error')
return
end
spm('TBlaunch',tb,I)
end
ok = 1;
% --------------------------------------------------------------------
function ok = check_diffusion_toolbox
% Check that Diffusion toolbox is on path
ok = 0;
if ~exist(which('dti_set_dtidata'),'file')
tb = spm('TBs');
I = strmatch('Diffusion',{tb(:).name});
if isempty(I) || I == 0
msgbox('Diffusion toolbox not installed!! Go to SPM website.',...
'MRI Analysis','error')
return
end
spm('TBlaunch',tb,I)
end
ok = 1;
|
