|
Dear Marko, Ged, All
I have successfully used the fieldmap undistortion in SPM5, with a few
tweaks to the fieldmap toolbox (updated Fieldmap.m attached). Firstly,
"spm_get"'s are replaced with "spm_select"'s, and in spm_vol structures, the
4th dimension of the dim field is replaced with the dt field.
I've also attached the file I use to apply the unwarping using the toolbox
- cspm_FieldMap_ngui.m, which is a hack of the FieldMap_ngui.m in the toolbox.
It assumes Siemens fieldmaps (phase/magnitude), which needs to be set in the
pm_defaults.m:
pm_def.INPUT_DATA_FORMAT = 'PM'; (default is 'RI')
I've put some help (ful?) comments at the start that should get you started.
I have made a further change to allow for the fieldmaps not quite covering
the whole brain:
line 1453 of Fieldmap.m
IP.uepiP.dat=uepi.*msk;
replaced with
IP.uepiP.dat=uepi;
This change came about because we missed a few mm in the posterior region of
a few fieldmaps, which led to the area being masked out in the unwarped images.
Hope this helps.
Best wishes,
Paul
On Mon, 6 Nov 2006 16:14:59 +0100, Marko Wilke
<[log in to unmask]> wrote:
>Dear Ged, All,
>
...
>
>I still wonder, therefore, if anyone has managed to create and use a
>fieldmap (voxel displacement map, to be exact) in spm5 ?
>
>Thanks (again) for all hints,
>best,
>Marko
>
function varargout=FieldMap(varargin)
%
% FieldMap is an SPM2 Toolbox for creating field maps and unwarping EPI.
% A full description of the toolbox and a usage manual can be found in
% FieldMap.man. This can launched by the toolbox help button or using
% spm_help FieldMap.man.The theoretical and practical principles behind
% the toolbox are described in principles.man.
%
% FORMAT FieldMap
%
% FieldMap launches the gui-based toolbox. Help is available via the
% help button (which calls spm_help FieldMap.man). FieldMap is a multi
% function function so that the toolbox routines can also be accessed
% without using the gui. A description of how to do this can be found
% in FieldMap_ngui.m
%
% Input parameters and the mode in which the toolbox works can be
% customised using the defaults file called pm_defaults.m.
%
% Main data structure:
%
% IP.P : 4x1 cell array containing real part short TE,
% imaginary part short TE, real part long TE and
% imaginary part long TE.
% IP.pP : Cell containing pre-calculated phase map. N.B.
% IP.P and IP.pP are mutually exclusive.
% IP.epiP : Cell containing EPI image used to demonstrate
% effects of unwarping.
% IP.uepiP : Cell containing unwarped EPI image.
% IP.nwarp : Cell containing non-distorted image.
% IP.vdmP : Cell containing the voxel displacement map (VDM)
% IP.et : 2x1 Cell array with short and long echo-time (ms).
% IP.epifm : Flag indicating EPI based field map (1) or not (0).
% IP.blipdir : Direction of phase-encode blips for k-space traversal
% (1 = positive or -1 = negative)
% IP.ajm : Flag indicating if Jacobian modulation should be applied
% (1) or not (0).
% IP.tert : Total epi readout time (ms).
%
% IP.uflags : Struct containing parameters guiding the unwrapping.
% Further explanations of these parameters are in
% FieldMap.man and pm_make_fieldmap.m
% .iformat : 'RI' or 'PM'
% .method : 'Huttonish', 'Mark3D' or 'Mark2D'
% .fwhm : FWHM (mm) of Gaussian filter for field map smoothing
% .pad : Size (in-plane voxels) of padding kernel.
% .etd : Echo time difference (ms).
%
% IP.fm : Struct containing field map information
% IP.fm.upm : Phase-unwrapped field map (Hz).
% IP.fm.mask : Binary mask excluding the noise in the phase map.
% IP.fm.opm : "Raw" field map (Hz) (not unwrapped).
% IP.fm.fpm : Phase-unwrapped, regularised field map (Hz).
% IP.fm.jac : Partial derivative of field map in y-direction.
%
% IP.vdm : Struct with voxel displacement map information
% IP.vdm.vdm : Voxel displacement map (scaled version of IP.fm.fpm).
% IP.vdm.jac : Jacobian-1 of forward transform.
% IP.vdm.ivdm : Inverse transform of voxel displacement
% (used to unwarp EPI image if field map is EPI based)
% (used to warp flash image prior to coregistration when
% field map is flash based (or other T2 weighting).
% IP.vdm.ijac : Jacobian-1 of inverse transform.
% IP.jim : Jacobian sampled in space of EPI.
%
% IP.cflags : Struct containing flags for coregistration
% (these are the default SPM coregistration flags -
% defaults.coreg).
% .cost_fun
% .sep
% .tol
% .fwhm
%
%_______________________________________________________________________
% Refs and Background reading:
%
% Jezzard P & Balaban RS. 1995. Correction for geometric distortion in
% echo planar images from Bo field variations. MRM 34:65-73.
%
% Hutton C et al. 2002. Image Distortion Correction in fMRI: A Quantitative
% Evaluation, NeuroImage 16:217-240.
%
% Cusack R & Papadakis N. 2002. New robust 3-D phase unwrapping
% algorithms: Application to magnetic field mapping and
% undistorting echoplanar images. NeuroImage 16:754-764.
%
% Jenkinson M. 2003. Fast, automated, N-dimensional phase-
% unwrapping algorithm. MRM 49:193-197.
%
%_______________________________________________________________________
% Acknowledegments
%
%_______________________________________________________________________
% FieldMap.m Jesper Andersson and Chloe Hutton 23/03/04
persistent PF FS WS PM % GUI related constants
persistent ID % Image display
persistent IP % Input and results
persistent DGW % Delete Graphics Window (if we created it)
global st % Global for spm_orthviews
global defaults
spm_defaults
if nargin == 0
Action = 'welcome';
else
Action = varargin{1};
end
%
% We have tried to divide the Actions into 3 categories:
% 1) Functions that can be called directly from the GUI are at the
% beginning.
% 2) Next come other GUI-related 'utility' functions - ie those that
% care of the GUI behind the scenes.
% 3) Finally, call-back functions that are not GUI dependent and can
% be called from a script are situated towards the end.
% See FieldMap_ngui.m for details on how to use these.
%
switch lower(Action)
%=======================================================================
%
% Create and initialise GUI for FieldMap toolbox
%
%=======================================================================
case 'welcome'
DGW = 0;
% Get all default values (these may effect GUI
IP = FieldMap('Initialise');
%
% Since we are using persistent variables we better make sure
% there is no-one else out there.
%
if ~isempty(PM)
figure(PM);
return
end
S = get(0,'ScreenSize');
WS = spm('WinScale');
FS = spm('FontSizes');
PF = spm_platform('fonts');
rect = [100 100 510 520].*WS;
PM = figure('IntegerHandle','off',...
'Name',sprintf('%s%s','FieldMap Toolbox, ver. beta 1.0',...
spm('GetUser',' (%s)')),...
'NumberTitle','off',...
'Tag','FieldMap',...
'Position',rect,...
'Resize','off',...
'Pointer','Arrow',...
'Color',[1 1 1]*.8,...
'MenuBar','none',...
'DefaultTextFontName',PF.helvetica,...
'DefaultTextFontSize',FS(12),...
'DefaultUicontrolFontName',PF.helvetica,...
'DefaultUicontrolFontSize',FS(12),...
'HandleVisibility','on',...
'Visible','on',...
'DeleteFcn','FieldMap(''Quit'');');
%-Create phase map controls
%-Frames and text
uicontrol(PM,'Style','Frame','BackgroundColor',spm('Colour'),...
'Position',[10 270 490 240].*WS);
uicontrol(PM,'Style','Text','String','Create field map in Hz',...
'Position',[100 480 300 020].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times,'FontAngle','Italic')
uicontrol(PM,'Style','Text','String','Short TE',...
'Position',[25 457 60 20].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
uicontrol(PM,'Style','Text','String','ms',...
'Position',[78 435 30 20].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
uicontrol(PM,'Style','Text','String','Long TE',...
'Position',[25 397 60 20].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
uicontrol(PM,'Style','Text','String','ms',...
'Position',[78 375 30 20].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
uicontrol(PM,'Style','Text','String',{'Precalculated','field map:'},...
'Position',[17 305 100 40].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
uicontrol(PM,'Style','Text','String',{'Field map value:'},...
'Position',[240 280 100 20].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
uicontrol(PM,'Style','Text','String',{'Hz'},...
'Position',[403 280 20 20].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
%-Objects with Callbacks
uicontrol(PM,'Style','Edit',...
'Position',[30 435 50 024].*WS,...
'ToolTipString','Give shortest echo-time in ms',...
'CallBack','FieldMap(''EchoTime'');',...
'UserData',1,...
'Tag','ShortTime');
uicontrol(PM,'Style','Edit',...
'Position',[30 375 50 024].*WS,...
'ToolTipString','Give longest echo-time in ms',...
'CallBack','FieldMap(''EchoTime'');',...
'UserData',2,...
'Tag','LongTime');
% Load Real and Imaginary or Phase and Magnitude Buttons
if IP.uflags.iformat=='PM'
uicontrol(PM,'String','Load Phase',...
'Position',[125 450 90 022].*WS,...
'ToolTipString','Load Analyze image containing first phase image',...
'CallBack','FieldMap(''LoadFilePM_Gui'');',...
'UserData',1,...
'Tag','LoadPhase1');
uicontrol(PM,'String','Load Mag.',...
'Position',[125 425 90 022].*WS,...
'ToolTipString','Load Analyze image containing first magnitude image',...
'CallBack','FieldMap(''LoadFilePM_Gui'');',...
'UserData',2,...
'Tag','LoadMag1');
uicontrol(PM,'String','Load Phase',...
'Position',[125 390 90 022].*WS,...
'ToolTipString','Load Analyze image containing second phase image',...
'CallBack','FieldMap(''LoadFilePM_Gui'');',...
'UserData',3,...
'Tag','LoadPhase2');
uicontrol(PM,'String','Load Mag.',...
'Position',[125 365 90 022].*WS,...
'ToolTipString','Load Analyze image containing second magnitudeimage',...
'CallBack','FieldMap(''LoadFilePM_Gui'');',...
'UserData',4,...
'Tag','LoadMag2');
else
uicontrol(PM,'String','Load Real',...
'Position',[125 450 90 022].*WS,...
'ToolTipString','Load Analyze image containing real part of short echo-time image',...
'CallBack','FieldMap(''LoadFile_Gui'');',...
'UserData',1,...
'Tag','LoadShortReal');
uicontrol(PM,'String','Load Imag',...
'Position',[125 425 90 022].*WS,...
'ToolTipString','Load Analyze image containing imaginary part of short echo-time image',...
'CallBack','FieldMap(''LoadFile_Gui'');',...
'UserData',2,...
'Tag','LoadShortImag');
uicontrol(PM,'String','Load Real',...
'Position',[125 390 90 022].*WS,...
'ToolTipString','Load Analyze image containing real part of long echo-time image',...
'CallBack','FieldMap(''LoadFile_Gui'');',...
'UserData',3,...
'Tag','LoadLongReal');
uicontrol(PM,'String','Load Imag',...
'Position',[125 365 90 022].*WS,...
'ToolTipString','Load Analyze image containing imaginary part of long echo-time image',...
'CallBack','FieldMap(''LoadFile_Gui'');',...
'UserData',4,...
'Tag','LoadLongImag');
end
uicontrol(PM,'String','Calculate',...
'Position',[265 337 90 022].*WS,...
'Enable','Off',...
'ToolTipString','Go ahead and create unwrapped field map (in Hz)',...
'CallBack','FieldMap(''CreateFieldmap_Gui'');',...
'Tag','CreateFieldMap');
uicontrol(PM,'String','Write',...
'Position',[370 337 90 022].*WS,...
'Enable','Off',...
'ToolTipString','Write Analyze image containing fininshed field map (in Hz)',...
'CallBack','FieldMap(''WriteFieldMap_Gui'');',...
'Tag','WriteFieldMap');
uicontrol(PM,'String','Load',...
'Position',[125 312 90 022].*WS,...
'ToolTipString','Load Analyze image containing fininshed field map (in Hz)',...
'CallBack','FieldMap(''LoadFieldMap_Gui'');',...
'Tag','LoadFieldMap');
% Empty string written to Fieldmap value box
uicontrol(PM,'Style','Text',...
'Position',[340 280 50 020].*WS,...
'HorizontalAlignment','left',...
'String','');
%-Create voxel-displacement-map controls
%-Frames and text
uicontrol(PM,'Style','Frame','BackgroundColor',spm('Colour'),...
'Position',[10 10 490 250].*WS);
uicontrol(PM,'Style','Text',...
'String','Create voxel displacement map (VDM) and unwarp EPI',...
'Position',[70 230 350 020].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times,'FontAngle','Italic')
uicontrol(PM,'Style','Text',...
'String','EPI-based field map',...
'Position',[50 200 200 020].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
uicontrol(PM,'Style','Text',...
'String','Polarity of phase-encode blips',...
'Position',[50 170 200 020].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
uicontrol(PM,'Style','Text',...
'String','Apply Jacobian modulation',...
'Position',[50 140 200 020].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
uicontrol(PM,'Style','Text',...
'String','Total EPI readout time',...
'Position',[50 110 200 020].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
uicontrol(PM,'Style','Text',...
'String','ms',...
'Position',[370 110 30 020].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times)
%-Objects with Callbacks
uicontrol(PM,'style','RadioButton',...
'String','Yes',...
'Position',[300 200 60 022].*WS,...
'ToolTipString','Field map is based on EPI data and needs inverting before use',...
'CallBack','FieldMap(''RadioButton'');',...
'Tag','EPI',...
'UserData',1);
uicontrol(PM,'style','RadioButton',...
'String','No',...
'Position',[370 200 60 022].*WS,...
'ToolTipString','Phase-map is based on non-EPI (Flash) data and can be used directly',...
'CallBack','FieldMap(''RadioButton'');',...
'Tag','EPI',...
'UserData',2);
uicontrol(PM,'style','RadioButton',...
'String','+ve',...
'Position',[300 170 60 022].*WS,...
'ToolTipString','K-space is traversed using positive phase-encode blips',...
'CallBack','FieldMap(''RadioButton'');',...
'Tag','BlipDir',...
'UserData',1);
uicontrol(PM,'style','RadioButton',...
'String','-ve',...
'Position',[370 170 60 022].*WS,...
'ToolTipString','K-space is traversed using negative phase-encode blips',...
'CallBack','FieldMap(''RadioButton'');',...
'Tag','BlipDir',...
'UserData',2);
uicontrol(PM,'style','RadioButton',...
'String','Yes',...
'Position',[300 140 60 022].*WS,...
'ToolTipString','Do Jacobian intensity modulation to compensate for compression/stretching',...
'CallBack','FieldMap(''RadioButton'');',...
'Tag','Jacobian',...
'UserData',1);
uicontrol(PM,'style','RadioButton',...
'String','No',...
'Position',[370 140 60 022].*WS,...
'ToolTipString','Don''t do Jacobian intensity modulation to compensate for compression/stretching',...
'CallBack','FieldMap(''RadioButton'');',...
'Tag','Jacobian',...
'UserData',2);
uicontrol(PM,'Style','Edit',...
'Position',[300 110 70 024].*WS,...
'ToolTipString','Give total time for readout of EPI echo train ms',...
'CallBack','FieldMap(''ReadOutTime'');',...
'Tag','ReadTime');
uicontrol(PM,'String','Load EPI image',...
'Position',[30 70 120 022].*WS,...
'ToolTipString','Load sample modulus EPI Analyze image',...
'CallBack','FieldMap(''LoadEpi_Gui'');',...
'Tag','LoadEpi');
uicontrol(PM,'String','Match VDM to EPI',...
'Position',[165 70 120 022].*WS,...
'ToolTipString','Match vdm to EPI (but only if you want to...)',...
'CallBack','FieldMap(''MatchVDM_gui'');',...
'Tag','MatchVDM');
uicontrol(PM,'String','Write unwarped',...
'Position',[300 70 120 022].*WS,...
'ToolTipString','Write unwarped EPI Analyze image',...
'CallBack','FieldMap(''WriteUnwarped_Gui'');',...
'Tag','WriteUnwarped');
uicontrol(PM,'String','Load structural',...
'Position',[30 30 120 022].*WS,...
'ToolTipString','Load structural image (but only if you want to...)',...
'CallBack','FieldMap(''LoadStructural_Gui'');',...
'Tag','LoadStructural');
uicontrol(PM,'String','Match structural',...
'Position',[164 30 120 022].*WS,...
'ToolTipString','Match structural image to EPI (but only if you want to...)',...
'CallBack','FieldMap(''MatchStructural_Gui'');',...
'Tag','MatchStructural');
uicontrol(PM,'String','Help',...
'Position',[320 30 80 022].*WS,...
'ToolTipString','Help',...
'CallBack','FieldMap(''Help_Gui'');',...
'ForegroundColor','g',...
'Tag','Help');
uicontrol(PM,'String','Quit',...
'Position',[420 30 60 022].*WS,...
'Enable','On',...
'ToolTipString','Exit toolbox',...
'CallBack','FieldMap(''Quit'');',...
'Tag','Quit');
%-Apply defaults
if ~isempty(IP.et{1})
h = findobj(get(PM,'Children'),'Tag','ShortTime');
set(h,'String',sprintf('%2.2f',IP.et{1}));
end
if ~isempty(IP.et{2})
h = findobj(get(PM,'Children'),'Tag','LongTime');
set(h,'String',sprintf('%2.2f',IP.et{2}));
end
if ~isempty(IP.epifm)
if IP.epifm
h = findobj(get(PM,'Children'),'Tag','EPI','UserData',1);
maxv = get(h,'Max');
set(h,'Value',maxv);
IP.epifm = 1;
else
h = findobj(get(PM,'Children'),'Tag','EPI','UserData',2);
maxv = get(h,'Max');
set(h,'Value',maxv);
IP.epifm = 0;
end
else % Default to yes
h = findobj(get(PM,'Children'),'Tag','EPI','UserData',1);
maxv = get(h,'Max');
set(h,'Value',maxv);
IP.epifm = 1;
end
if ~isempty(IP.ajm)
if IP.ajm
h = findobj(get(PM,'Children'),'Tag','Jacobian','UserData',1);
maxv = get(h,'Max');
set(h,'Value',maxv);
IP.ajm = 1;
else
h = findobj(get(PM,'Children'),'Tag','Jacobian','UserData',2);
maxv = get(h,'Max');
set(h,'Value',maxv);
IP.ajm = 0;
end
else % Default to no
h = findobj(get(PM,'Children'),'Tag','Jacobian','UserData',2);
maxv = get(h,'Max');
set(h,'Value',maxv);
IP.ajm = 0;
end
if ~isempty(IP.blipdir)
if IP.blipdir
h = findobj(get(PM,'Children'),'Tag','BlipDir','UserData',1);
maxv = get(h,'Max');
set(h,'Value',maxv);
IP.blipdir = 1;
else
h = findobj(get(PM,'Children'),'Tag','BlipDir','UserData',2);
maxv = get(h,'Max');
set(h,'Value',maxv);
IP.blipdir = -1;
end
else % Default to 1 = positive
h = findobj(get(PM,'Children'),'Tag','BlipDir','UserData',1);
maxv = get(h,'Max');
set(h,'Value',maxv);
IP.blipdir = 1;
end
if ~isempty(IP.tert)
h = findobj(get(PM,'Children'),'Tag','ReadTime');
set(h,'String',sprintf('%2.1f',IP.tert));
end
%
% Disable necessary buttons and parameters until phase-map is finished.
%
FieldMap('Reset_Gui');
%=======================================================================
%
% Load real or imaginary part of dual echo-time data - gui version.
%
%=======================================================================
case 'loadfile_gui'
FieldMap('Reset_Gui');
%
% File selection using gui
%
index = get(gcbo,'UserData');
FieldMap('LoadFile',index);
ypos = [445 420 385 360];
uicontrol(PM,'Style','Text','String',spm_str_manip(IP.P{index}.fname,['l',num2str(50)]),...
'Position',[220 ypos(index) 260 20].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times,...
'FontSize',FS(8));
FieldMap('Assert');
%=======================================================================
%
% Load phase and magnitude images - gui version.
%
%=======================================================================
case 'loadfilepm_gui'
FieldMap('Reset_Gui');
%
% File selection using gui
%
index = get(gcbo,'UserData');
FieldMap('LoadFilePM',index);
ypos = [445 420 385 360];
uicontrol(PM,'Style','Text','String',spm_str_manip(IP.P{index}.fname,['l',num2str(50)]),...
'Position',[220 ypos(index) 260 20].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times,...
'FontSize',FS(8));
FieldMap('Assert');
%=======================================================================
%
% Create unwrapped phase-map - gui version
%
%=======================================================================
case 'createfieldmap_gui'
%
% Create fieldmap from complex images...
%
status=FieldMap('CreateFieldMap',IP);
if ~isempty(status)
FieldMap('DisplayImage',FieldMap('MakedP'),[.05 .75 .95 .2],1);
% Toggle relevant buttons
FieldMap('ToggleGUI','Off','CreateFieldMap');
FieldMap('ToggleGUI','On',str2mat('LoadFieldMap','WriteFieldMap'));
%
% Check that have correct parameters ready before allowing
% an image to be loaded for unwarping and hence conversion of
% fieldmap to voxel shift map.
%
if (FieldMap('GatherVDMData'))
FieldMap('FM2VDM',IP);
FieldMap('ToggleGUI','On',str2mat('LoadEpi','EPI','BlipDir',...
'Jacobian','ReadTime'));
end
end
%=======================================================================
%
% Load already prepared fieldmap (in Hz).
%
%=======================================================================
case 'loadfieldmap_gui'
%
% Reset all previously loaded field maps and images to unwarp
%
FieldMap('Reset_gui');
%
% Select a precalculated phase map
%
FieldMap('LoadFieldMap');
uicontrol(PM,'Style','Text','String',spm_str_manip(IP.pP.fname,['l',num2str(50)]),...
'Position',[220 307 260 20].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'),...
'FontName',PF.times,...
'FontSize',FS(8));
FieldMap('DisplayImage',FieldMap('MakedP'),[.05 .75 .95 .2],1);
FieldMap('ToggleGUI','Off',str2mat('CreateFieldMap','WriteFieldMap',...
'MatchVDM','WriteUnwarped',...
'LoadStructural', 'MatchStructural'));
% Check that have correct parameters ready before allowing
% an image to be loaded for unwarping and hence conversion of
% fieldmap to voxel shift map.
%
if (FieldMap('GatherVDMData'))
FieldMap('FM2VDM',IP);
FieldMap('ToggleGUI','On',str2mat('LoadEpi','EPI','BlipDir',...
'Jacobian','ReadTime'));
end
%=======================================================================
%
% Write out field map in Hz - using GUI
%
%=======================================================================
case 'writefieldmap_gui'
FieldMap('Write',IP.P{1},IP.fm.fpm,'fpm_',64,'Fitted phase map');
%=======================================================================
%
% Load sample EPI image using gui, unwarp and display result
%
%=======================================================================
case 'loadepi_gui'
%
% Select and display EPI image
%
FieldMap('LoadEPI');
FieldMap('DisplayImage',IP.epiP,[.05 .5 .95 .2],2);
% The fm2vdm parameters may have been changed so must calculate
% the voxel shift map with current parameters.
if (FieldMap('GatherVDMData'))
FieldMap('FM2VDM',IP);
FieldMap('UnwarpEPI',IP);
FieldMap('DisplayImage',IP.uepiP,[.05 .25 .95 .2],3);
% Once EPI has been unwarped can enable unwarp checks and structural stuff
FieldMap('ToggleGUI','On',str2mat('MatchVDM','WriteUnwarped','LoadStructural'));
end
%
% Redisplay phasemap in case .mat file associated with
% EPI image different from images that filadmap was
% estimated from.
%
FieldMap('DisplayImage',FieldMap('MakedP'),[.05 .75 .95 .2],1);
%=======================================================================
%
% Coregister fieldmap magnitude image to EPI to unwarp using GUI
%
%=======================================================================
case 'matchvdm_gui'
if (FieldMap('GatherVDMData'))
FieldMap('MatchVDM',IP);
FieldMap('DisplayImage',FieldMap('MakedP'),[.05 .75 .95 .2],1)
FieldMap('UnwarpEPI',IP);
FieldMap('DisplayImage',IP.epiP,[.05 .5 .95 .2],2);
FieldMap('DisplayImage',IP.uepiP,[.05 .25 .95 .2],3);
% Once EPI has been unwarped can enable unwarp checks
FieldMap('ToggleGUI','On',str2mat('MatchVDM','WriteUnwarped','LoadStructural'));
end
%=======================================================================
%
% Load structural image
%
%=======================================================================
case 'loadstructural_gui'
FieldMap('LoadStructural');
FieldMap('DisplayImage',IP.nwarp,[.05 0.0 .95 .2],4);
%
% Redisplay the other images to allow for recalculation
% of size of display.
%
FieldMap('DisplayImage',FieldMap('MakedP'),[.05 .75 .95 .2],1);
FieldMap('DisplayImage',IP.epiP,[.05 .5 .95 .2],2);
FieldMap('DisplayImage',IP.uepiP,[.05 .25 .95 .2],3);
FieldMap('ToggleGUI','On','MatchStructural');
%=======================================================================
%
% Match structural image to unwarped EPI
%
%=======================================================================
case 'matchstructural_gui'
FieldMap('MatchStructural',IP);
FieldMap('DisplayImage',IP.nwarp,[.05 0.0 .95 .2],4);
%=======================================================================
%
% Write out unwarped EPI
%
%=======================================================================
case 'writeunwarped_gui'
IP.uepiP = FieldMap('Write',IP.epiP,IP.uepiP.dat,'u',IP.epiP.dt(1),'Unwarped image');
FieldMap('ToggleGUI','On', 'LoadStructural');
FieldMap('ToggleGUI','Off', 'WriteUnwarped');
%=======================================================================
%
% Write out voxel shift map
%
%=======================================================================
case 'writevdm_gui'
if ~isempty(IP.pP)
if IP.epifm == 1
FieldMap('Write',IP.pP,IP.vdm.ivdm,'vdm_',4,'Voxel shift map');
else
FieldMap('Write',IP.pP,IP.vdm.vdm,'vdm_',4,'Voxel shift map');
end
else
if IP.epifm == 1
FieldMap('Write',IP.P,IP.vdm.ivdm,'vdm_',4,'Voxel shift map');
else
FieldMap('Write',IP.P,IP.vdm.vdm,'vdm_',4,'Voxel shift map');
end
end
%=======================================================================
%
% Help using spm_help
%
%=======================================================================
case 'help_gui'
FieldMap('help');
%=======================================================================
%
% Quit Toolbox
%
%=======================================================================
case 'quit'
Fgraph=spm_figure('FindWin','Graphics');
if ~isempty(Fgraph)
if DGW
delete(Fgraph);
Fgraph=[];
DGW = 0;
else
spm_figure('Clear','Graphics');
end
end
PM=spm_figure('FindWin','FieldMap');
if ~isempty(PM)
delete(PM);
PM=[];
end
%=======================================================================
%
% The following functions are GUI related functions that go on behind
% the scenes.
%=======================================================================
%
%=======================================================================
%
% Reset gui inputs when a new field map is to be calculated
%
%=======================================================================
case 'reset_gui'
% Clear any precalculated fieldmap name string
uicontrol(PM,'Style','Text','String','',...
'Position',[230 307 260 20].*WS,...
'ForegroundColor','k','BackgroundColor',spm('Colour'));
% Clear Fieldmap values
uicontrol(PM,'Style','Text',...
'Position',[350 280 50 020].*WS,...
'HorizontalAlignment','left',...
'String','');
% Disable input routes to make sure
% a new fieldmap is calculated.
%
FieldMap('ToggleGUI','Off',str2mat('CreateFieldMap','WriteFieldMap',...
'LoadEpi','WriteUnwarped','MatchVDM',...
'LoadStructural','MatchStructural',...
'EPI','BlipDir',...
'Jacobian','ReadTime'));
%
% A new input image implies all processed data is void.
%
IP.fm = [];
IP.vdm = [];
IP.jim = [];
IP.pP = [];
IP.epiP = [];
IP.uepiP = [];
IP.vdmP = [];
ID = cell(4,1);
spm_orthviews('Reset');
spm_figure('Clear','Graphics');
%=======================================================================
%
% Echo time was changed or entered
%
%=======================================================================
case 'echotime'
FieldMap('Assert');
%=======================================================================
%
% Check if inputs are correct for calculating new phase map
%
%=======================================================================
case 'assert'
if ~isempty(IP.pP) % If we're using precalculated fieldmap.
go = 0;
else
go = 1;
for i=1:2
if isempty(IP.P{i}) go = 0; end
end
for i=3:4
if (isempty(IP.P{i}) & IP.uflags.iformat=='RI') go = 0; end
end
h = findobj(get(PM,'Children'),'Tag','ShortTime');
IP.et{1} = str2num(get(h,'String'));
h = findobj(get(PM,'Children'),'Tag','LongTime');
IP.et{2} = str2num(get(h,'String'));
if isempty(IP.et{1}) | isempty(IP.et{2}) | IP.et{2} < IP.et{1}
go = 0;
end
end
if go
FieldMap('ToggleGui','On','CreateFieldMap');
else % Switch off fieldmap creation
FieldMap('ToggleGui','Off','CreateFieldMap');
end
%=======================================================================
%
% Enable or disable specified buttons or parameters in GUI
%
%=======================================================================
case 'togglegui'
h = get(PM,'Children');
tags=varargin{3};
for n=1:size(varargin{3},1)
set(findobj(h,'Tag',deblank(tags(n,:))),'Enable',varargin{2});
end
%=======================================================================
%
% A radiobutton was pressed.
%
%=======================================================================
case 'radiobutton'
%
% Enforce radio behaviour.
%
index = get(gcbo,'UserData');
if index==1
partner=2;
else
partner=1;
end
tag = get(gcbo,'Tag');
value = get(gcbo,'Value');
maxv = get(gcbo,'Max');
h = findobj(get(PM,'Children'),'Tag',tag,'UserData',partner);
if value == maxv
set(h,'Value',get(h,'Min'));
else
set(h,'Value',get(h,'Max'));
end
% Update Hz to voxel displacement map if the input parameters are ok
if (FieldMap('GatherVDMData'))
FieldMap('FM2VDM',IP);
%
% Update unwarped EPI if one is loaded
%
if ~isempty(IP.epiP)
IP.uepiP = FieldMap('UnwarpEPI',IP);
FieldMap('DisplayImage',IP.uepiP,[.05 .25 .95 .2],3);
FieldMap('ToggleGUI','On',str2mat('WriteUnwarped'));
end
end
%=======================================================================
%
% Total readout-time was changed or entered
%
%=======================================================================
case 'readouttime'
%
% Update unwarped EPI
%
% This means the transformation phase-map to voxel displacement-map
% has changed, which means the inversion will have to be redone,
% which means (in the case of flash-based field map) coregistration
% between field-map and sample EPI image should be redone. Phew!
%
if (FieldMap('GatherVDMData'))
FieldMap('FM2VDM',IP);
% Update unwarped EPI if one is loaded
if ~isempty(IP.epiP)
IP.uepiP = FieldMap('UnwarpEPI',IP);
FieldMap('DisplayImage',IP.uepiP,[.05 .25 .95 .2],3);
FieldMap('ToggleGUI','On',str2mat('WriteUnwarped',...
'LoadStructural'));
end
FieldMap('ToggleGUI','On',str2mat('LoadEpi','EPI','BlipDir',...
'Jacobian','ReadTime'));
else
% If readtime is missing switch everything off...
FieldMap('ToggleGUI','Off',str2mat('LoadEpi','EPI',...
'BlipDir','Jacobian',...
'WriteUnwarped','LoadStructural',...
'MatchStructural', 'MatchVDM'));
end
%=======================================================================
%
% Sample UIControls pertaining to information needed for transformation
% phase-map -> voxel displacement-map
%
%=======================================================================
case 'gathervdmdata'
h = findobj(get(PM,'Children'),'Tag','EPI','UserData',1);
v = get(h,{'Value','Max'});
if v{1} == 1
IP.epifm = 1;
else
IP.epifm = 0;
end
h = findobj(get(PM,'Children'),'Tag','BlipDir','UserData',1);
v = get(h,{'Value','Max'});
if v{1} == 1 %% CHLOE: Changed
IP.blipdir = 1;
else
IP.blipdir = -1;
end
h = findobj(get(PM,'Children'),'Tag','Jacobian','UserData',1);
v = get(h,{'Value','Max'});
if v{1} == 1 %% CHLOE: Changed
IP.ajm = 1;
else
IP.ajm = 0;
end
h = findobj(get(PM,'Children'),'Tag','ReadTime');
IP.tert = str2num(get(h,'String'));
if isempty(IP.tert) | isempty(IP.fm) | isempty(IP.fm.fpm)
varargout{1} = 0;
FieldMap('ToggleGui','On','ReadTime');
else
varargout{1} = 1;
end
%=======================================================================
%
% Draw transversal and sagittal image.
%
%=======================================================================
case 'redrawimages'
global curpos;
for i=1:length(ID)
if ~isempty(ID{i}) & ~isempty(st.vols{i})
set(st.vols{i}.ax{2}.ax,'Visible','Off'); % Disable event delivery in Coronal
set(st.vols{i}.ax{2}.d,'Visible','Off'); % Make Coronal invisible
set(st.vols{i}.ax{1}.ax,'Position',ID{i}.tra_pos);
set(st.vols{i}.ax{1}.ax,'ButtonDownFcn',['FieldMap(''Orthviews'');']);
set(get(st.vols{i}.ax{1}.ax,'YLabel'),'String',ID{i}.label);
set(st.vols{i}.ax{3}.ax,'Position',ID{i}.sag_pos);
set(st.vols{i}.ax{3}.ax,'ButtonDownFcn',['FieldMap(''Orthviews'');']);
end
end
%=======================================================================
%
% Callback for orthviews
%
%=======================================================================
case 'orthviews'
if strcmp(get(gcf,'SelectionType'),'normal')
spm_orthviews('Reposition');,...
%spm_orthviews('set_pos2cm');,...
elseif strcmp(get(gcf,'SelectionType'),'extend')
spm_orthviews('Reposition');,...
%spm_orthviews('set_pos2cm');,...
spm_orthviews('context_menu','ts',1);
end;
curpos = spm_orthviews('pos',1);
set(st.in, 'String',sprintf('%3.3f',spm_sample_vol(st.vols{1},curpos(1),curpos(2),curpos(3),st.hld)));
%=======================================================================
%
% Display image.
%
%=======================================================================
case 'displayimage'
Fgraph=spm_figure('FindWin','Graphics');
if isempty(Fgraph)
st.fig=spm_figure('Create','Graphics','Graphics');
DGW = 1;
end
if ~isempty(ID{varargin{4}})
spm_orthviews('Delete',ID{varargin{4}}.h);
ID{varargin{4}} = [];
end
h = spm_orthviews('Image',varargin{2},[.01 .01 .01 .01]);
% Set bounding box to allow display of all images
spm_orthviews('MaxBB');
% Put everything in space of original EPI image.
if varargin{4} == 2
% spm_orthviews('Space',h);
end
%
% Get best possible positioning and scaling for
% transversal and sagittal views.
%
tra_pos = get(st.vols{varargin{4}}.ax{1}.ax,'Position');
sag_pos = get(st.vols{varargin{4}}.ax{3}.ax,'Position');
field_pos = varargin{3};
x_scale = field_pos(3) / (tra_pos(3) + sag_pos(3));
height = max([tra_pos(4) sag_pos(4)]);
y_scale = field_pos(4) / height;
if x_scale > y_scale % Height limiting
scale = y_scale;
dx = (field_pos(3) - scale*(tra_pos(3) + sag_pos(3))) / 3;
dy = 0;
else % Width limiting
scale = x_scale;
dx = 0;
dy = (field_pos(4) - scale*height) / 2;
end
tra_pos = [field_pos(1)+dx field_pos(2)+dy scale*tra_pos([3 4])];
sag_pos = [field_pos(1)+tra_pos(3)+2*dx field_pos(2)+dy scale*sag_pos([3 4])];
label = {'Fieldmap in Hz',...
'Warped EPI',...
'Unwarped EPI',...
'Structural'};
ID{varargin{4}} = struct('tra_pos', tra_pos,...
'sag_pos', sag_pos,...
'h', h,...
'label', label{varargin{4}});
FieldMap('RedrawImages');
st.in = uicontrol(PM,'Style','Text',...
'Position',[340 280 50 020].*WS,...
'HorizontalAlignment','left',...
'String','');
%=======================================================================
%=======================================================================
%
% The functions below are called by the various gui buttons but are
% not dependent on the gui to work. These functions can therefore also
% be called from a script bypassing the gui and can return updated
% variables.
%
%=======================================================================
%=======================================================================
%
%=======================================================================
%
% Create and initialise parameters for FieldMap toolbox
%
%=======================================================================
case 'initialise'
%
% Initialise parameters
%
IP.P = cell(1,4);
IP.pP = [];
IP.epiP = [];
IP.uepiP = [];
IP.nwarp = [];
IP.fm = [];
IP.vdm = [];
IP.et = cell(1,2);
IP.epifm = [];
IP.blipdir = [];
IP.ajm = [];
IP.tert = [];
IP.uflags = struct('iformat','','method','','fwhm',[],'pad',[],'etd',[],'ws',[]);
% Get default parameters
pm_defaults;
% Define parameters for fieldmap creation
IP.et{1} = pm_def.SHORT_ECHO_TIME;
IP.et{2} = pm_def.LONG_ECHO_TIME;
% Set parameters for unwrapping
IP.uflags.iformat = pm_def.INPUT_DATA_FORMAT;
IP.uflags.method = pm_def.UNWRAPPING_METHOD;
IP.uflags.fwhm = pm_def.FWHM;
IP.uflags.pad = pm_def.PAD;
IP.uflags.ws = pm_def.WS;
IP.uflags.etd = pm_def.LONG_ECHO_TIME - pm_def.SHORT_ECHO_TIME;
% Set parameters for unwarping
IP.ajm = pm_def.DO_JACOBIAN_MODULATION;
IP.blipdir = pm_def.K_SPACE_TRAVERSAL_BLIP_DIR;
IP.tert = pm_def.TOTAL_EPI_READOUT_TIME;
IP.epifm = pm_def.EPI_BASED_FIELDMAPS;
varargout{1}= IP;
%=======================================================================
%
% Load real or imaginary part of dual echo-time data - NO gui.
%
%=======================================================================
case 'loadfile'
index=varargin{2};
prompt_string = {'Select short echo-time real',...
'Select short echo-time imaginary',...
'Select long echo-time real',...
'Select long echo-time imaginary'};
IP.P{index} = spm_vol(spm_get(1,'*.img',prompt_string{index}));
if isfield(IP,'pP') & ~isempty(IP.pP)
IP.pP = [];
end
varargout{1} = IP.P{index};
%=======================================================================
%
% Load phase or magnitude part of dual (possibly) echo-time data - NO gui.
%
%=======================================================================
case 'loadfilepm'
index=varargin{2};
prompt_string = {'Select phase image',...
'Select magnitude image',...
'Select second phase image or press done',...
'Select magnitude image or press done'};
if index<3
IP.P{index} = spm_vol(spm_get(1,'*.img',prompt_string{index}));
else
IP.P{index} = spm_vol(spm_get([0 1],'*.img',prompt_string{index}));
end
if isfield(IP,'pP') & ~isempty(IP.pP)
IP.pP = [];
end
varargout{1} = IP.P{index};
%=======================================================================
%
% Load already prepared fieldmap (in Hz) - no gui
%
%=======================================================================
case 'loadfieldmap'
%
% Select field map
%
% 24/03/04 - Chloe change below to spm_get(1,'fpm_*.img')...
IP.pP = spm_vol(spm_get(1,'fpm_*.img','Select field map'));
IP.fm.fpm = spm_read_vols(IP.pP);
IP.fm.jac = pm_diff(IP.fm.fpm,2);
if isfield(IP,'P') & ~isempty(IP.P{1})
IP.P = cell(1,4);
end
varargout{1} = IP.fm;
varargout{2} = IP.pP;
%=======================================================================
%
% Create unwrapped phase-map - NO gui
%
%=======================================================================
case 'createfieldmap'
IP=varargin{2};
% First check that images are in same space.
if size([IP.P{1} IP.P{2} IP.P{3} IP.P{4}],2)==4
ip_dim=cat(1,[IP.P{1}.dim' IP.P{2}.dim' IP.P{3}.dim' IP.P{4}.dim']');
ip_mat=cat(2,[IP.P{1}.mat(:) IP.P{2}.mat(:) IP.P{3}.mat(:) IP.P{4}.mat(:)]');
else
ip_dim=cat(1,[IP.P{1}.dim' IP.P{2}.dim']');
ip_mat=cat(2,[IP.P{1}.mat(:) IP.P{2}.mat(:)]');
end
if any(any(diff(ip_dim,1,1),1)&[1,1,1])
errordlg({'Images don''t all have same dimensions'});
drawnow;
varargout{1}=[];
elseif any(any(diff(ip_mat,1,1),1))
errordlg({'Images don''t all have same orientation & voxel size'});
drawnow;
varargout{1}=[];
else
% Update flags for unwarping (in case TEs have been adjusted
IP.uflags.etd = IP.et{2}-IP.et{1};
IP.fm = pm_make_fieldmap([IP.P{1} IP.P{2} IP.P{3} IP.P{4}],IP.uflags);
varargout{1} = IP.fm;
end
%=======================================================================
%
% Convert field map to voxel displacement map and
% do necessary inversions of displacement fields.
%
%=======================================================================
case 'fm2vdm'
IP=varargin{2};
%
% If we already have memory mapped image pointer to voxel
% displacement map let's reuse it (so not to void possible
% realignment). If field-map is non-EPI based the previous
% realignment (based on different parameters) will be non-
% optimal and we should advice user to redo it.
%
% If no pointer already exist we'll make one.
%
if isfield(IP,'vdmP') & ~isempty(IP.vdmP)
if IP.epifm ~= 1
msgbox({'Changing this parameter means that if previously',...
'you matched VDM to EPI, this result may no longer',...
'be optimal. In this case we recommend you redo the',...
'Match VDM to EPI.'},...
'Coregistration notice','warn','modal');
end
else
if ~isempty(IP.pP)
IP.vdmP = struct('dim', IP.pP.dim,...
'dt',[4 spm_platform('bigend')],...
'mat', IP.pP.mat);
IP.vdmP.fname=fullfile(spm_str_manip(IP.pP.fname, 'h'),['vdm_' deblank(spm_str_manip(IP.pP.fname,'t'))]);
else
IP.vdmP = struct('dim', IP.P{1}.dim,...
'dt',[4 spm_platform('bigend')],...
'mat', IP.P{1}.mat);
IP.vdmP.fname=fullfile(spm_str_manip(IP.P{1}.fname, 'h'),['vdm_' deblank(spm_str_manip(IP.P{1}.fname,'t'))]);
end
end
% Scale field map and jacobian by total EPI readout time
IP.vdm.vdm = IP.blipdir*IP.tert*1e-3*IP.fm.fpm;
IP.vdm.jac = IP.blipdir*IP.tert*1e-3*IP.fm.jac;
if IP.epifm ==1
spm_progress_bar('Init',3,'Inverting displacement map','');
spm_progress_bar('Set',1);
% Invert voxel shift map and multiply by -1...
IP.vdm.ivdm = -1*pm_invert_phasemap(IP.vdm.vdm);
IP.vdm.ijac = -1*pm_diff(IP.vdm.ivdm,2);
spm_progress_bar('Set',2);
spm_progress_bar('Clear');
FieldMap('Write',IP.vdmP,IP.vdm.ivdm,'',IP.vdmP.dt(1),'Voxel shift map');
else
FieldMap('Write',IP.vdmP,IP.vdm.vdm,'',IP.vdmP.dt(1),'Voxel shift map');
end
varargout{1} = IP.vdm;
varargout{2} = IP.vdmP;
%=======================================================================
%
% Write out image
%
%=======================================================================
case 'write'
V=varargin{2};
vol=varargin{3};
name=varargin{4};
% Write out image
img=struct('dim', V.dim,...
'dt',[varargin{5} spm_platform('bigend')],...
'mat', V.mat);
img.fname=fullfile(spm_str_manip(V.fname, 'h'),[name deblank(spm_str_manip(V.fname,'t'))]);
img.descrip=varargin{6};
%img = spm_create_vol(img);
img=spm_write_vol(img,vol);
varargout{1}=img;
%=======================================================================
%
% Create fieldmap (Hz) struct for use when displaying image.
%
%=======================================================================
case 'makedp'
if isfield(IP,'vdmP') & ~isempty(IP.vdmP);
dP = struct('dim', IP.vdmP.dim,...
'dt',[64 spm_platform('bigend')],...
'pinfo', [1 0]',...
'mat', IP.vdmP.mat,...
'dat', reshape(IP.fm.fpm,IP.vdmP.dim(1:3)),...
'fname', 'display_image');
else
if isfield(IP,'P') & ~isempty(IP.P{1})
dP = struct('dim', IP.P{1}.dim,...
'dt',[64 spm_platform('bigend')],...
'pinfo', [1 0]',...
'mat', IP.P{1}.mat,...
'dat', reshape(IP.fm.fpm,IP.P{1}.dim(1:3)),...
'fname', 'display_image');
elseif isfield(IP,'pP') & ~isempty(IP.pP)
dP = struct('dim', IP.pP.dim,...
'dt',[64 spm_platform('bigend')],...
'pinfo', [1 0]',...
'mat', IP.pP.mat,...
'dat', reshape(IP.fm.fpm,IP.pP.dim(1:3)),...
'fname', 'display_image');
end
end
varargout{1} = dP;
%=======================================================================
%
% Load sample EPI image - NO gui
%
%=======================================================================
case 'loadepi'
%
% Select EPI
%
IP.epiP = spm_vol(spm_get(1,'*.img','Select sample EPI image'));
varargout{1} = IP.epiP;
%=======================================================================
%
% Create unwarped epi - NO gui
%
%=======================================================================
case 'unwarpepi'
%
% Update unwarped EPI
%
IP=varargin{2};
IP.uepiP = struct('fname', 'Image in memory',...
'dim', IP.epiP.dim,...
'dt',[64 spm_platform('bigend')],...
'pinfo', IP.epiP.pinfo(1:2),...
'mat', IP.epiP.mat);
% Need to sample EPI and voxel shift map in space of EPI...
[x,y,z] = ndgrid(1:IP.epiP.dim(1),1:IP.epiP.dim(2),1:IP.epiP.dim(3));
xyz = [x(:) y(:) z(:)];
% Space of EPI is IP.epiP{1}.mat and space of
% voxel shift map is IP.vdmP{1}.mat
tM = inv(IP.epiP.mat\IP.vdmP.mat);
x2 = tM(1,1)*x + tM(1,2)*y + tM(1,3)*z + tM(1,4);
y2 = tM(2,1)*x + tM(2,2)*y + tM(2,3)*z + tM(2,4);
z2 = tM(3,1)*x + tM(3,2)*y + tM(3,3)*z + tM(3,4);
xyz2 = [x2(:) y2(:) z2(:)];
%
% Make mask since it is only meaningful to calculate undistorted
% image in areas where we have information about distortions.
%
msk = reshape(double(xyz2(:,1)>=1 & xyz2(:,1)<=IP.vdmP.dim(1) &...
xyz2(:,2)>=1 & xyz2(:,2)<=IP.vdmP.dim(2) &...
xyz2(:,3)>=1 & xyz2(:,3)<=IP.vdmP.dim(3)),IP.epiP.dim(1:3));
% Read in voxel displacement map in correct space
tvdm = reshape(spm_sample_vol(spm_vol(IP.vdmP.fname),xyz2(:,1),...
xyz2(:,2),xyz2(:,3),1),IP.epiP.dim(1:3));
% Voxel shift map must be subtracted from the y-coordinates
uepi = reshape(spm_sample_vol(IP.epiP,xyz(:,1),...
xyz(:,2)-tvdm(:),xyz(:,3),1),IP.epiP.dim(1:3));
% Sample Jacobian in correct space and apply if required
if IP.ajm==1
if IP.epifm==1 % If EPI, use inverted jacobian
IP.jim = reshape(spm_sample_vol(IP.vdm.ijac,xyz2(:,1),...
xyz2(:,2),xyz2(:,3),1),IP.epiP.dim(1:3));
else
IP.jim = reshape(spm_sample_vol(IP.vdm.jac,xyz2(:,1),...
xyz2(:,2),xyz2(:,3),1),IP.epiP.dim(1:3));
end
uepi = uepi.*(1+IP.jim);
end
%IP.uepiP.dat=uepi.*msk; PMM
IP.uepiP.dat=uepi;
varargout{1}=IP.uepiP;
%=======================================================================
%
% Coregister fieldmap magnitude image to EPI to unwarp
%
%=======================================================================
case 'matchvdm'
IP=varargin{2};
%
% Need a fieldmap magnitude image
%
if isempty(IP.pP) & ~isempty(IP.P{1})
IP.fmagP=struct('dim', IP.P{1}.dim,...
'dt',IP.P{1}.dt,...
'pinfo', IP.P{1}.pinfo,...
'mat', IP.P{1}.mat);
IP.fmagP.fname=fullfile(spm_str_manip(IP.P{1}.fname, 'h'),['mag_' deblank(spm_str_manip(IP.P{1}.fname,'t'))]);
% If using real and imaginary data, calculate using sqrt(i1.^2 + i2.^2).
% If using phase and magnitude, use magnitude image.
if IP.uflags.iformat=='RI'
IP.fmagP = spm_imcalc(spm_vol([IP.P{1}.fname;IP.P{2}.fname]),IP.fmagP,'sqrt(i1.^2 + i2.^2)');
else
IP.fmagP = IP.P{2};
end
else
IP.fmagP = spm_vol(spm_get(1,'*.img','Select field map magnitude image'));
end
% Now we have field map magnitude image, we want to coregister it to the
% EPI to be unwarped.
% If using an EPI field map:
% 1) Coregister magnitude image to EPI.
% 2) Apply resulting transformation matrix to voxel shift map
% If using a non-EPI field map:
% 1) Forward warp magnitude image
% 2) Coregister warped magnitude image to EPI.
% 3) Apply resulting transformation matrix to voxel shift map
if IP.epifm==1
[mi,M] = FieldMap('Coregister',IP.epiP,IP.fmagP);
MM = IP.fmagP.mat;
else
% Need to sample magnitude image in space of EPI to be unwarped...
[x,y,z] = ndgrid(1:IP.epiP.dim(1),1:IP.epiP.dim(2),1:IP.epiP.dim(3));
xyz = [x(:) y(:) z(:)];
% Space of EPI is IP.epiP{1}.mat and space of fmagP is IP.fmagP.mat
tM = inv(IP.epiP.mat\IP.fmagP.mat);
x2 = tM(1,1)*x + tM(1,2)*y + tM(1,3)*z + tM(1,4);
y2 = tM(2,1)*x + tM(2,2)*y + tM(2,3)*z + tM(2,4);
z2 = tM(3,1)*x + tM(3,2)*y + tM(3,3)*z + tM(3,4);
xyz2 = [x2(:) y2(:) z2(:)];
wfmag = reshape(spm_sample_vol(IP.fmagP,xyz2(:,1),...
xyz2(:,2),xyz2(:,3),1),IP.epiP.dim(1:3));
% Invert voxel shift map since it is a forward warp.
vdm = pm_invert_phasemap(IP.vdm.vdm);
% Need to sample voxel shift map in space of EPI to be unwarped
tvdm = reshape(spm_sample_vol(vdm,xyz2(:,1),...
xyz2(:,2),xyz2(:,3),0),IP.epiP.dim(1:3));
% Now apply warps to resampled forward warped magnitude image...
wfmag = reshape(spm_sample_vol(wfmag,xyz(:,1),xyz(:,2)-tvdm(:),...
xyz(:,3),1),IP.epiP.dim(1:3));
% Write out forward warped magnitude image
IP.wfmagP = struct('dim', IP.epiP.dim,...
'dt',[64 spm_platform('bigend')],...
'pinfo', IP.epiP.pinfo,...
'mat', IP.epiP.mat);
IP.wfmagP = FieldMap('Write',IP.epiP,wfmag,'wfmag_',4,'Voxel shift map');
% Now coregister warped magnitude field map to EPI
[mi,M] = FieldMap('Coregister',IP.epiP,IP.wfmagP);
% Update the .mat file of the forward warped mag image
spm_get_space(deblank(IP.wfmagP.fname),M*IP.wfmagP.mat);
% Get the original space of the fmap magnitude
MM = IP.fmagP.mat;
end
% Update .mat file for voxel displacement map
IP.vdmP.mat=M*MM;
spm_get_space(deblank(IP.vdmP.fname),M*MM);
varargout{1} = IP.vdmP;
%=======================================================================
%
% Load structural image
%
%=======================================================================
case 'loadstructural'
IP.nwarp = spm_vol(spm_get(1,'*.img','Select structural image'));
varargout{1} = IP.nwarp;
%=======================================================================
%
% Coregister structural image to unwarped EPI
%
%=======================================================================
case 'matchstructural'
IP = varargin{2};
[mi,M] = FieldMap('Coregister',IP.uepiP,IP.nwarp);
MM = IP.nwarp.mat;
IP.nwarp.mat=M*MM;
spm_get_space(deblank(IP.nwarp.fname),IP.nwarp.mat);
varargout{1}=mi;
%=======================================================================
%
% Coregister two images - NO gui
%
%=======================================================================
case 'coregister'
%
% Coregisters image VF to image VG (use VF and VG like in SPM code)
%
VG = varargin{2};
VF = varargin{3};
% Define flags for coregistration...
IP.cflags.cost_fun = defaults.coreg.estimate.cost_fun;
IP.cflags.sep = defaults.coreg.estimate.sep;
IP.cflags.tol = defaults.coreg.estimate.tol;
IP.cflags.fwhm = defaults.coreg.estimate.fwhm;
IP.cflags.graphics = 0;
% Voxel sizes (mm)
vxg = sqrt(sum(VG.mat(1:3,1:3).^2));
vxf = sqrt(sum(VF.mat(1:3,1:3).^2));
% Smoothnesses
fwhmg = sqrt(max([1 1 1]*IP.cflags.sep(end)^2 - vxg.^2, [0 0 0]))./vxg;
fwhmf = sqrt(max([1 1 1]*IP.cflags.sep(end)^2 - vxf.^2, [0 0 0]))./vxf;
% Need to load data smoothed in unit8 format (as in spm_coreg_ui)
% if isfield(VG,'dat')
% VG=rmfield(VG,'dat');
% end
if ~isfield(VG,'uint8')
VG.uint8 = loaduint8(VG);
VG = smooth_uint8(VG,fwhmg); % Note side effects
end
% if isfield(VF,'dat')
% VF=rmfield(VF,'dat');
% end
if ~isfield(VF,'uint8')
VF.uint8 = loaduint8(VF);
VF = smooth_uint8(VF,fwhmf); % Note side effects
end
x=spm_coreg(VG,VF,IP.cflags);
M = inv(spm_matrix(x));
MM = spm_get_space(deblank(VF.fname));
varargout{1}=spm_coreg(x,VG,VF,2,IP.cflags.cost_fun,IP.cflags.fwhm);
varargout{2}=M;
%=======================================================================
%
% Use spm_help to display help for FieldMap toolbox
%
%=======================================================================
case 'help'
spm_help('FieldMap.man');
return
end
%=======================================================================
%
% Smoothing functions required for spm_coreg
%
%=======================================================================
function V = smooth_uint8(V,fwhm)
% Convolve the volume in memory (fwhm in voxels).
lim = ceil(2*fwhm);
s = fwhm/sqrt(8*log(2));
x = [-lim(1):lim(1)]; x = smoothing_kernel(fwhm(1),x); x = x/sum(x);
y = [-lim(2):lim(2)]; y = smoothing_kernel(fwhm(2),y); y = y/sum(y);
z = [-lim(3):lim(3)]; z = smoothing_kernel(fwhm(3),z); z = z/sum(z);
i = (length(x) - 1)/2;
j = (length(y) - 1)/2;
k = (length(z) - 1)/2;
spm_conv_vol(V.uint8,V.uint8,x,y,z,-[i j k]);
return;
%_______________________________________________________________________
%_______________________________________________________________________
function krn = smoothing_kernel(fwhm,x)
% Variance from FWHM
s = (fwhm/sqrt(8*log(2)))^2+eps;
% The simple way to do it. Not good for small FWHM
% krn = (1/sqrt(2*pi*s))*exp(-(x.^2)/(2*s));
% For smoothing images, one should really convolve a Gaussian
% with a sinc function. For smoothing histograms, the
% kernel should be a Gaussian convolved with the histogram
% basis function used. This function returns a Gaussian
% convolved with a triangular (1st degree B-spline) basis
% function.
% Gaussian convolved with 0th degree B-spline
% int(exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= -0.5..0.5)
% w1 = 1/sqrt(2*s);
% krn = 0.5*(erf(w1*(x+0.5))-erf(w1*(x-0.5)));
% Gaussian convolved with 1st degree B-spline
% int((1-t)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t= 0..1)
% +int((t+1)*exp(-((x+t))^2/(2*s))/sqrt(2*pi*s),t=-1..0)
w1 = 0.5*sqrt(2/s);
w2 = -0.5/s;
w3 = sqrt(s/2/pi);
krn = 0.5*(erf(w1*(x+1)).*(x+1) + erf(w1*(x-1)).*(x-1) - 2*erf(w1*x ).* x)...
+w3*(exp(w2*(x+1).^2) + exp(w2*(x-1).^2) - 2*exp(w2*x.^2));
krn(krn<0) = 0;
return;
%=======================================================================
%
% Load data function required for spm_coreg
%
%=======================================================================
function udat = loaduint8(V)
% Load data from file indicated by V into an array of unsigned bytes.
if size(V.pinfo,2)==1 & V.pinfo(1) == 2,
mx = 255*V.pinfo(1) + V.pinfo(2);
mn = V.pinfo(2);
else,
spm_progress_bar('Init',V.dim(3),...
['Computing max/min of ' spm_str_manip(V.fname,'t')],...
'Planes complete');
mx = -Inf; mn = Inf;
for p=1:V.dim(3),
img = spm_slice_vol(V,spm_matrix([0 0 p]),V.dim(1:2),1);
mx = max([max(img(:))+paccuracy(V,p) mx]);
mn = min([min(img(:)) mn]);
spm_progress_bar('Set',p);
end;
end;
spm_progress_bar('Init',V.dim(3),...
['Loading ' spm_str_manip(V.fname,'t')],...
'Planes loaded');
udat = uint8(0);
udat(V.dim(1),V.dim(2),V.dim(3))=0;
rand('state',100);
for p=1:V.dim(3),
img = spm_slice_vol(V,spm_matrix([0 0 p]),V.dim(1:2),1);
acc = paccuracy(V,p);
if acc==0,
udat(:,:,p) = uint8(round((img-mn)*(255/(mx-mn))));
else,
% Add random numbers before rounding to reduce aliasing artifact
r = rand(size(img))*acc;
udat(:,:,p) = uint8(round((img+r-mn)*(255/(mx-mn))));
end;
spm_progress_bar('Set',p);
end;
spm_progress_bar('Clear');
return;
function acc = paccuracy(V,p)
if ~spm_type(V.dt(1),'intt'),
acc = 0;
else,
if size(V.pinfo,2)==1,
acc = abs(V.pinfo(1,1));
else,
acc = abs(V.pinfo(1,p));
end;
end;
%_______________________________________________________________________
function Pout = cspm_FieldMap_ngui( 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
%
% 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
Pin = spm_select(Inf,'image','Select phase image');
end
if nargin < 3
Pin = spm_select(Inf,'image','Select magnitude image');
end
if nargin < 4
TE = spm_input('Enter two echo-times (TE''s in ms)',1,'r','0 0',2);
end
if nargin < 5
type = 'fmri';
elseif strcmpi(type,'dti')
check_diffusion_toolbox
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
check_fieldmap_toolbox
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 check_fieldmap_toolbox
% Check that fieldmap toolbox is on path
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
% --------------------------------------------------------------------
function check_diffusion_toolbox
% Check that Diffusion toolbox is on path
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
|
