Skip to content

Commit

Permalink
New function CoRegStandAlone.m to co-register and segment MRS data (n…
Browse files Browse the repository at this point in the history 
…ot necessarily edited) in all supported format. Small changes to SiemensDICOMRead.m and DICOMRead.m that fixed flawed post-processing of co-registering and segmenting DICOM MRS data.
  • Loading branch information
@schorschinho
schorschinho committed Sep 19, 2018
1 parent 156813f commit 147f5b5
Show file tree
Hide file tree
Showing 7 changed files with 605 additions and 51 deletions.
Binary file modified .DS_Store
View file
Binary file not shown.
214 changes: 214 additions & 0 deletions CoRegStandAlone/CoReg.m
View file
@@ -0,0 +1,214 @@
function MRS_struct = CoReg(MRS_struct, nii_name)

%Coregistration of MRS voxel volumes to imaging datasets, based on headers.

MRS_struct.version.coreg = '180718';

% First check if SPM12 is installed and on the search path
spmversion = fileparts(which('spm'));
if isempty(spmversion)
error('SPM not found! Please install SPM12 (https://www.fil.ion.ucl.ac.uk/spm/software/spm12) and make sure it is on your search path.');
elseif strcmpi(spmversion(end-3:end),'spm8')
error(['SPM8 detected! Gannet 3.0 no longer supports SPM8. ' ...
'Please install SPM12 (https://www.fil.ion.ucl.ac.uk/spm/software/spm12) and make sure it is on your search path.']);
end

if MRS_struct.ii ~= length(nii_name)
error('The number of nifti files does not match the number of MRS files processed by CoRegStandAlone.');
end

numscans = numel(MRS_struct.metabfile);
vox = {MRS_struct.p.Vox{1}};

for ii = 1:numscans
% Loop over voxels if PRIAM
for kk = 1:length(vox)

%Ultimately this switch will not be necessary...
switch MRS_struct.p.vendor

case 'Philips'
fname = MRS_struct.metabfile{ii};
sparname = [fname(1:(end-4)) MRS_struct.p.spar_string];
MRS_struct = GannetMask_Philips(sparname, nii_name{ii}, MRS_struct, ii, vox, kk);

case 'Philips_data'
if exist(MRS_struct.metabfile_sdat,'file') % MM (170720)
MRS_struct.p.vendor = 'Philips';
MRS_struct.metabfile_data = MRS_struct.metabfile;
MRS_struct.metabfile = MRS_struct.metabfile_sdat;
MRS_struct = GannetCoRegister(MRS_struct,nii_name);
MRS_struct.metabfile = MRS_struct.metabfile_data;
MRS_struct.p.vendor = 'Philips_data';
else
error([MRS_struct.p.vendor ' format does not include voxel location information in the header. See notes in GannetCoRegister.']);
%If this comes up, once GannetLoad has been read:
%1. Switch vendor to Philips
% MRS_struct.p.vendor = 'Philips';
%2. Copy .data filenames.
% MRS_struct.metabfile_data = MRS_struct.metabfile;
%3. Replace the list with the corrsponding SDAT files (in correct order)
% MRS_struct.metabfile = {'SDATfile1.sdat' 'SDATfile2.SDAT'};
%4. Rerun GannetCoRegister
%
%5. Copy .sdat filenames and replace .data ones. Tidy up.
% MRS_struct.metabfile_sdat = MRS_struct.metabfile;
% MRS_struct.metabfile = MRS_struct.metabfile_data;
% MRS_struct.p.vendor = 'Philips_data'
end

case 'Siemens_rda'
fname = MRS_struct.metabfile{ii};
MRS_struct = GannetMask_SiemensRDA(fname, nii_name{ii}, MRS_struct, ii, vox, kk);

case {'Siemens_twix', 'Siemens_dicom'}
fname = MRS_struct.metabfile{ii};
MRS_struct = GannetMask_SiemensTWIX(fname, nii_name{ii}, MRS_struct, ii, vox, kk);

case 'GE'
fname = MRS_struct.metabfile{ii};
MRS_struct = GannetMask_GE(fname, nii_name{ii}, MRS_struct, ii, vox, kk);

end

% Build output figure
if ishandle(103)
clf(103); % MM (170720)
end
h = figure(103);
% MM (170629): Open figure in center of screen
scr_sz = get(0, 'ScreenSize');
fig_w = 1000;
fig_h = 707;
set(h, 'Position', [(scr_sz(3)-fig_w)/2, (scr_sz(4)-fig_h)/2, fig_w, fig_h]);
set(h,'Color',[1 1 1]);
figTitle = 'GannetCoRegister Output';
set(gcf,'Name',figTitle,'Tag',figTitle,'NumberTitle','off');

subplot(2,3,4:6)
axis off;

tmp = 'Mask output: ';
text(0.5, 0.75, tmp,'HorizontalAlignment','right', ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);

[~,tmp,tmp2] = fileparts(MRS_struct.mask.(vox{kk}).outfile{ii});
text(0.5, 0.75, [' ' tmp tmp2], ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13,'Interpreter','none');

tmp = 'Spatial parameters: ';
text(0.5, 0.63, tmp, 'HorizontalAlignment', 'right', ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);
tmp = ' [LR, AP, FH]';
text(0.5, 0.63, tmp, ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);

tmp = 'Dimension: ';
text(0.5, 0.51, tmp, 'HorizontalAlignment', 'right', ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);
tmp = [' [' num2str(MRS_struct.p.voxdim(ii,1)) ', ' num2str(MRS_struct.p.voxdim(ii,2)) ', ' num2str(MRS_struct.p.voxdim(ii,3)) '] mm'];
text(0.5, 0.51, tmp, ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);

tmp = 'Volume: ';
text(0.5, 0.39, tmp, 'HorizontalAlignment', 'right', ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);
vol = MRS_struct.p.voxdim(ii,1)*MRS_struct.p.voxdim(ii,2)*MRS_struct.p.voxdim(ii,3)*.001;
tmp = [' ' num2str(vol) ' mL'];
text(0.5, 0.39, tmp, ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);

tmp = 'Position: ';
text(0.5, 0.27, tmp, 'HorizontalAlignment', 'right', ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);
tmp = [' [' num2str(MRS_struct.p.voxoff(ii,1), '%3.1f') ', ' num2str(MRS_struct.p.voxoff(ii,2), '%3.1f') ', ' num2str(MRS_struct.p.voxoff(ii,3), '%3.1f') '] mm'];
text(0.5, 0.27, tmp, ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);

tmp = 'Angulation: ';
text(0.5, 0.15, tmp, 'HorizontalAlignment', 'right', ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);
tmp = [' [' num2str(MRS_struct.p.voxang(ii,1), '%3.1f') ', ' num2str(MRS_struct.p.voxang(ii,2), '%3.1f') ', ' num2str(MRS_struct.p.voxang(ii,3), '%3.1f') '] deg'];
text(0.5, 0.15, tmp, ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);

tmp = 'CoRegVer: ';
text(0.5, 0.03, tmp, 'HorizontalAlignment', 'right', ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);
tmp = [' ' MRS_struct.version.coreg];
text(0.5, 0.03, tmp, ...
'VerticalAlignment', 'top', ...
'FontName', 'Helvetica','FontSize',13);

h = subplot(2,3,1:3);

% MM (180112)
[~,tmp,tmp2] = fileparts(MRS_struct.metabfile{ii});
[~,tmp3,tmp4] = fileparts(MRS_struct.mask.(vox{kk}).T1image{ii});
t = ['Voxel from ' tmp tmp2 ' on ' tmp3 tmp4];

imagesc(squeeze(MRS_struct.mask.(vox{kk}).img{ii}));
colormap('gray');
img = MRS_struct.mask.(vox{kk}).img{ii};
img = img(:);
caxis([0 mean(img(img>0.01)) + 3*std(img(img>0.01))]); % MM (180807)
axis equal;
axis tight;
axis off;
text(10,size(MRS_struct.mask.(vox{kk}).img{ii},1)/2,'L','Color',[1 1 1],'FontSize',20);
text(size(MRS_struct.mask.(vox{kk}).img{ii},2)-20,size(MRS_struct.mask.(vox{kk}).img{ii},1)/2,'R','Color',[1 1 1],'FontSize',20);
get(h,'pos');
set(h,'pos',[0.0 0.15 1 1]);
title(t, 'FontName', 'Helvetica', 'FontSize', 15, 'Interpreter', 'none');

script_path = which('GannetLoad');
Gannet_logo = [script_path(1:(end-13)) '/Gannet3_logo.png'];
A2 = imread(Gannet_logo,'png','BackgroundColor',[1 1 1]);
axes('Position',[0.80, 0.05, 0.15, 0.15]);
image(A2);
axis off;
axis square;

% For Philips .data
if strcmpi(MRS_struct.p.vendor,'Philips_data')
fullpath = MRS_struct.metabfile{ii};
fullpath = regexprep(fullpath, '.data', '_data');
fullpath = regexprep(fullpath, '\', '_');
fullpath = regexprep(fullpath, '/', '_');
end

% MM (180112)
[~,metabfile_nopath] = fileparts(MRS_struct.metabfile{ii});

% Save PDF output
set(gcf,'PaperUnits','inches');
set(gcf,'PaperSize',[11 8.5]);
set(gcf,'PaperPosition',[0 0 11 8.5]);
if strcmpi(MRS_struct.p.vendor,'Philips_data')
pdfname = fullfile('CoRegStandAlone_output', [fullpath '_' vox{kk} '_coreg.pdf']); % MM (180112)
else
pdfname = fullfile('CoRegStandAlone_output', [metabfile_nopath '_' vox{kk} '_coreg.pdf']); % MM (180112)
end
saveas(gcf, pdfname);


end

end




128 changes: 128 additions & 0 deletions CoRegStandAlone/CoRegStandAlone.m
View file
@@ -0,0 +1,128 @@
function MRS_struct = CoRegStandAlone(metabfile,niifile)
% CoRegStandAlone(metabfile,niifile)
% Reads the relevant header information from MRS files, performs
% SPM-based co-registration to a 3D anatomical image in nifti (*.nii)
% format, and returns the tissue class fractions of gray matter, white
% matter, and cerebrospinal fluid in the MRS voxel.
%
% Requires:
% - SPM 12 (https://www.fil.ion.ucl.ac.uk/spm/software/spm12/)
%
% Input:
% metabfile - cell containing the path to an MRS data file in one of
% the following formats:
% Philips SDAT/SPAR, Siemens TWIX (*.dat), Siemens
% RDA (*.rda), GE (.7), DICOM (*.ima, *.dcm)
% niifile - cell containing the path to a *.nii file
%
% It is possible to batch process data. In this case, the number of
% elements of the metabfile cell needs to match the number of
% elements of the niifile cell, with matching order.
%
% Example:
% MRS_struct = CoRegStandAlone({'MRS1.dat, 'MRS2.dat'}, {'nii.dat', 'nii.dat'})
%
% This example will co-register and segment two Siemens MRS voxels to
% the same structural.
%
% Author:
% Dr. Georg Oeltzschner (Johns Hopkins University, 2018-09-19)
% goeltzs1@jhmi.edu
%
% History:
% 2018-09-19: First version of the code.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Pre-initialise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

MRS_struct.version.load = '180326'; % set to date when final updates have been made
MRS_struct.ii = 0;
MRS_struct.metabfile = metabfile;

% Flags
MRS_struct.p.mat = 1; % Save results in *.mat output structure? (0 = NO, 1 = YES (default)).
MRS_struct.p.Vox = {'vox1'}; % Name of the voxel

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. Determine data parameters from header
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Discern input data format
MRS_struct = GannetDiscernDatatype(metabfile{1}, MRS_struct);

% Create output folder
if ~exist('CoRegStandAlone_output','dir')
mkdir CoRegStandAlone_output;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3. Load data from files
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for ii = 1:length(metabfile) % Loop over all files in the batch (from metabfile)

MRS_struct.ii = ii;

switch MRS_struct.p.vendor

case 'GE'
MRS_struct = GERead(MRS_struct, metabfile{ii});

case 'Siemens_twix'
MRS_struct = SiemensTwixRead(MRS_struct, metabfile{ii});

case 'Siemens_dicom' % GO 11/01/2016
MRS_struct = SiemensDICOMRead(MRS_struct,metabfile{ii}); % GO 11/01/2016

case 'dicom' % GO 11/30/2016
MRS_struct = DICOMRead(MRS_struct,metabfile{ii}); % GO 11/01/2016

case 'Siemens_rda'
MRS_struct = SiemensRead(MRS_struct, metabfile{ii}, metabfile{ii});

case 'Philips'
MRS_struct = PhilipsRead(MRS_struct, metabfile{ii});

case 'Philips_data'
MRS_struct = PhilipsRead_data(MRS_struct, metabfile{ii});

case 'Philips_raw'
MRS_struct = PhilipsRawLoad(MRS_struct,metabfile{ii},3,0); % GO 11/02/2016

end % end of vendor switch loop for data load

end % end of load-and-processing loop over datasets

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 4. Call coregister function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for ii = 1:length(metabfile) % Loop over all files in the batch (from metabfile)
MRS_struct = CoReg(MRS_struct, niifile);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 5. Call segment function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for ii = 1:length(metabfile) % Loop over all files in the batch (from metabfile)
MRS_struct = Seg(MRS_struct);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 6. Clean up, save data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Save MRS_struct as mat file
MRS_struct = rmfield(MRS_struct,'fids');
if ii == length(metabfile) && MRS_struct.p.mat
% Set up filename
mat_name = ['CoRegStandAlone_output/MRS_struct_CoRegStandAlone.mat'];
save(mat_name,'MRS_struct');
end

end



0 comments on commit 147f5b5

Please sign in to comment.