diff --git a/.DS_Store b/.DS_Store index dc62a7f..3fc6363 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/CoRegStandAlone/CoReg.m b/CoRegStandAlone/CoReg.m new file mode 100755 index 0000000..dfbcd87 --- /dev/null +++ b/CoRegStandAlone/CoReg.m @@ -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 + + + + diff --git a/CoRegStandAlone/CoRegStandAlone.m b/CoRegStandAlone/CoRegStandAlone.m new file mode 100755 index 0000000..e1f6dfe --- /dev/null +++ b/CoRegStandAlone/CoRegStandAlone.m @@ -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 + + + diff --git a/CoRegStandAlone/Seg.m b/CoRegStandAlone/Seg.m new file mode 100644 index 0000000..1cf1d28 --- /dev/null +++ b/CoRegStandAlone/Seg.m @@ -0,0 +1,218 @@ +function MRS_struct = Seg(MRS_struct) + +% Relies on SPM being installed +% +% Runs segmentation script if segmented images not present according to +% file convention of c1, c2 and c3 as prefixes on the anatomical image name +% for the GM, WM and CSF segmentations. If these files are present, they +% are loaded and used for the voxel segmentation +% +% This script does not require any information from the fit. +% This is useful if only the tissue segmentation information is supposed to +% be obtained. + +MRS_struct.version.segment = '180403'; +vox = {MRS_struct.p.Vox{1}}; + +% Set up SPM for batch processing +spm('defaults','fmri'); +spm_jobman('initcfg'); +kk = 1; +for ii = 1:length(MRS_struct.metabfile) + + % 1 - Take nifti from GannetCoRegister and segment it in SPM + + [T1dir, T1name, T1ext] = fileparts(MRS_struct.mask.(vox{kk}).T1image{ii}); + anatimage = MRS_struct.mask.(vox{kk}).T1image{ii}; + + % Check to see if segmentation already done - if not, do it + % Check which SPM version is installed and segment accordingly + tmp = [T1dir '/c1' T1name T1ext]; + if ~exist(tmp,'file') + spmversion = fileparts(which('spm')); + if strcmpi(spmversion(end-4:end),'spm12') + CallSPM12segmentation(anatimage); + else + CallSPM8segmentation(anatimage); + end + end + + % 2 - Determine GM, WM and CSF fractions for each voxel + + if strcmp(T1dir,'') + T1dir='.'; + end + + GM = [T1dir '/c1' T1name T1ext]; + WM = [T1dir '/c2' T1name T1ext]; + CSF = [T1dir '/c3' T1name T1ext]; + + GMvol = spm_vol(GM); + WMvol = spm_vol(WM); + CSFvol = spm_vol(CSF); + + % Loop over voxels if PRIAM + for kk = 1:length(vox) + + voxmaskvol = spm_vol(cell2mat(MRS_struct.mask.(vox{kk}).outfile(ii))); + + % GM + O_GMvox.fname = [T1dir '/c1' T1name '_GM_mask.nii']; + O_GMvox.descrip = 'GMmasked_MRS_Voxel_Mask'; + O_GMvox.dim = voxmaskvol.dim; + O_GMvox.dt = voxmaskvol.dt; + O_GMvox.mat = voxmaskvol.mat; + GM_voxmask_vol = GMvol.private.dat(:,:,:) .* voxmaskvol.private.dat(:,:,:); + O_GMvox = spm_write_vol(O_GMvox, GM_voxmask_vol); + + % WM + O_WMvox.fname = [T1dir '/c2' T1name '_WM_mask.nii']; + O_WMvox.descrip = 'WMmasked_MRS_Voxel_Mask'; + O_WMvox.dim = voxmaskvol.dim; + O_WMvox.dt = voxmaskvol.dt; + O_WMvox.mat = voxmaskvol.mat; + WM_voxmask_vol = WMvol.private.dat(:,:,:) .* voxmaskvol.private.dat(:,:,:); + O_WMvox = spm_write_vol(O_WMvox, WM_voxmask_vol); + + % CSF + O_CSFvox.fname = [T1dir '/c3' T1name '_CSF_mask.nii']; + O_CSFvox.descrip = 'CSFmasked_MRS_Voxel_Mask'; + O_CSFvox.dim = voxmaskvol.dim; + O_CSFvox.dt = voxmaskvol.dt; + O_CSFvox.mat = voxmaskvol.mat; + CSF_voxmask_vol = CSFvol.private.dat(:,:,:) .* voxmaskvol.private.dat(:,:,:); + O_CSFvox = spm_write_vol(O_CSFvox, CSF_voxmask_vol); + + % ***NB*** + % For a subject with multiple voxels, the segmented voxels will + % get overwritten + + % 3 - Calculate an adjusted gabaiu and output it to the structure + + gmsum = sum(sum(sum(O_GMvox.private.dat(:,:,:)))); + wmsum = sum(sum(sum(O_WMvox.private.dat(:,:,:)))); + csfsum = sum(sum(sum(O_CSFvox.private.dat(:,:,:)))); + + gmfra = gmsum/(gmsum+wmsum+csfsum); + wmfra = wmsum/(gmsum+wmsum+csfsum); + csffra = csfsum/(gmsum+wmsum+csfsum); + + tissuefra = gmfra+wmfra; + + MRS_struct.out.(vox{kk}).tissue.GMfra(ii) = gmfra; + MRS_struct.out.(vox{kk}).tissue.WMfra(ii) = wmfra; + MRS_struct.out.(vox{kk}).tissue.CSFfra(ii) = csffra; + + % 4 - Build output + if ishandle(104) + clf(104); % MM (170831) + end + h = figure(104); + % MM (170831): 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 = 'GannetSegment Output'; + set(gcf,'Name',figTitle,'Tag',figTitle,'NumberTitle','off'); + + % Voxel co-registration + subplot(2,2,1); + size_max = size(MRS_struct.mask.(vox{kk}).img{ii},1); + imagesc(MRS_struct.mask.(vox{kk}).img{ii}(:,size_max+(1:size_max))); + colormap('gray'); + caxis([0 0.5]); + axis equal; + axis tight; + axis off; + + % 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]; + title(t, 'Interpreter', 'none'); + + % Output results + subplot(2,2,2); + axis off; + + text_pos = 1; + + tmp1 = 'GM voxel fraction'; + tmp2 = sprintf(': %.3g', MRS_struct.out.(vox{kk}).tissue.GMfra(ii)); + text(0, text_pos-0.1, tmp1, 'FontName', 'Helvetica', 'FontSize', 10); + text(0.5, text_pos-0.1, tmp2, 'FontName', 'Helvetica', 'FontSize', 10); + + tmp1 = 'WM voxel fraction'; + tmp2 = sprintf(': %.3g', MRS_struct.out.(vox{kk}).tissue.WMfra(ii)); + text(0, text_pos-0.2, tmp1, 'FontName', 'Helvetica', 'FontSize', 10); + text(0.5, text_pos-0.2, tmp2, 'FontName', 'Helvetica', 'FontSize', 10); + + tmp1 = 'CSF voxel fraction'; + tmp2 = sprintf(': %.3g', MRS_struct.out.(vox{kk}).tissue.CSFfra(ii)); + text(0, text_pos-0.3, tmp1, 'FontName', 'Helvetica', 'FontSize', 10); + text(0.5, text_pos-0.3, tmp2, 'FontName', 'Helvetica', 'FontSize', 10); + + % MM (180112) + tmp1 = 'Filename'; + if strcmp(MRS_struct.p.vendor,'Siemens_rda') + [~,tmp2,tmp3] = fileparts(MRS_struct.metabfile{ii*2-1}); + else + [~,tmp2,tmp3] = fileparts(MRS_struct.metabfile{ii}); + end + text(0, text_pos-0.4, tmp1, 'FontName', 'Helvetica', 'FontSize', 10); + text(0.5, text_pos-0.4, [': ' tmp2 tmp3], 'FontName', 'Helvetica', 'FontSize', 10, 'Interpreter', 'none'); + + tmp1 = 'Anatomical image'; + [~,tmp2,tmp3] = fileparts(MRS_struct.mask.(vox{kk}).T1image{ii}); % MM (180112) + text(0, text_pos-0.5, tmp1, 'FontName', 'Helvetica', 'FontSize', 10); + text(0.5, text_pos-0.5, [': ' tmp2 tmp3], 'FontName', 'Helvetica', 'FontSize', 10, 'Interpreter', 'none'); + + tmp1 = 'SegmentVer'; + tmp2 = [': ' MRS_struct.version.segment]; + text(0, text_pos-0.6, tmp1, 'FontName', 'Helvetica', 'FontSize', 10); + text(0.5, text_pos-0.6, tmp2, 'FontName', 'Helvetica', 'FontSize', 10); + + % Gannet logo + subplot(2,2,4); + axis off; + script_path = which('GannetFit'); + Gannet_logo = [script_path(1:(end-12)) '/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 '_segment.pdf']); % MM (180112) + else + pdfname = fullfile('CoRegStandAlone_output', [metabfile_nopath '_segment.pdf']); % MM (180112) + end + saveas(gcf, pdfname); + + + end + +end + +end + + + diff --git a/DICOMRead.m b/DICOMRead.m index 8985e31..01a1a44 100644 --- a/DICOMRead.m +++ b/DICOMRead.m @@ -1,19 +1,17 @@ -function MRS_struct = DICOMRead(MRS_struct,folder,waterfolder) -%% MRS_struct = DICOMRead(MRS_struct,folder,waterfolder) +function MRS_struct = DICOMRead(MRS_struct,metabfile,waterfile) +%% MRS_struct = DICOMRead(MRS_struct,metabfile,waterfile) % This function is designed to load edited MR spectroscopy data in the % general form of DICOM data into a Gannet file structure. Files usually % have the extension '.DCM' or '.dcm', and contain exactly 1 FID per -% file, i.e. an acquisition of 320 averages will yield 320 IMA files. +% file, i.e. an acquisition of 320 averages will yield 320 DCM files. % -% The user must specify the folder containing all of these averages. It -% is assumed that they are ordered in the order of acquisition. +% It is assumed that they are ordered in the order of acquisition. % Water-suppressed and water-unsuppressed files need to be stored in -% separate folders, which need to be defined accordingly: +% separate folders, e.g. '/user/data/subject01/dcm_gaba/' and +% '/user/data/subject01/dcm_water/', respectively. % % Example: -% folder = '/user/data/subject01/dcm_gaba/'; -% waterfolder = '/user/data/subject01/dcm_water/'; (optional) -% MRS_struct = DICOMRead(MRS_struct,folder,waterfolder); +% MRS_struct = DICOMRead(MRS_struct,'/user/data/subject01/dcm_gaba/metab.dcm','/user/data/subject01/dcm_water/water.dcm'); % % Author: % Dr. Georg Oeltzschner (Johns Hopkins University, 2016-11-10) @@ -32,7 +30,9 @@ % of Minnesota) (2017-11-20). Thanks to Jim Lagopoulos. % 0.95: Fills missing voxel geometry parameters in DICOM header with zero % values. Thanks to Alen Tersakyan. -% +% 0.96: Fixed to accomodate batch processing of coregister/segmentation. +% (2018-09-19) +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -46,6 +46,7 @@ % no overlap. % dcm_file_list = [dir([folder,'/*.DCM']); dir([folder,'/*.dcm'])]; % may % cause problems on win/unix systems, take out for now % GO 11/16/2016 +[folder,~,~] = fileparts(metabfile); % GO 11/01/2016 dcm_file_list = dir([folder,'/*.dcm']); % GO 11/16/2016 fprintf('%d water-suppressed DCM files detected in %s.\n',length(dcm_file_list),folder); @@ -60,7 +61,7 @@ %%% /PREPARATION %%% %%% HEADER INFO PARSING %%% -DicomHeader = read_dcm_header(dcm_file_names{1}); +DicomHeader = read_dcm_header(metabfile); MRS_struct.p.seq = DicomHeader.sequenceFileName; MRS_struct.p.TR(ii) = DicomHeader.TR; MRS_struct.p.TE(ii) = DicomHeader.TE; @@ -100,12 +101,14 @@ % It appears that IMA stores the transients weirdly, 1-n/2 are all ONs, and % n/2-n are all OFFS. Shuffle them below. -a = MRS_struct.fids.data(:,1:end/2); -b = MRS_struct.fids.data(:,1+end/2:end); -c = zeros(size(MRS_struct.fids.data)); -c(:,1:2:end) = a; -c(:,2:2:end) = b; -MRS_struct.fids.data = c; +if size(MRS_struct.fids.data,2) > 1 + a = MRS_struct.fids.data(:,1:end/2); + b = MRS_struct.fids.data(:,1+end/2:end); + c = zeros(size(MRS_struct.fids.data)); + c(:,1:2:end) = a; + c(:,2:2:end) = b; + MRS_struct.fids.data = c; +end %%% /DATA LOADING %%% @@ -117,6 +120,7 @@ % Set up the file name array. if nargin == 3 + [waterfolder,~,~] = fileparts(waterfile); water_file_list = dir([waterfolder,'/*.DCM']); fprintf('%d water-unsuppressed DCM files detected in %s.\n',length(water_file_list),waterfolder); disp('Reading water-unsuppressed files...') diff --git a/GannetLoad.m b/GannetLoad.m index 7c9ab23..c36bb3e 100755 --- a/GannetLoad.m +++ b/GannetLoad.m @@ -47,7 +47,7 @@ % 1. Pre-initialise %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -MRS_struct.version.load = '180326'; % set to date when final updates have been made +MRS_struct.version.load = '180919'; % set to date when final updates have been made MRS_struct.ii = 0; MRS_struct.metabfile = metabfile; MRS_struct = GannetPreInitialise(MRS_struct); @@ -136,17 +136,10 @@ case 'Siemens_dicom' % GO 11/01/2016 if exist('waterfile','var') - % Load the data. GannetLoad specifies a file in the first - % place, so take apart that filename, and feed the containing - % folder into the SiemensDICOMRead function. % GO 11/01/2016 - [imafolder,~,~] = fileparts(metabfile{ii}); % GO 02/05/2017 - [waterfolder,~,~] = fileparts(waterfile{ii}); % GO 02/05/2017 - MRS_struct = SiemensDICOMRead(MRS_struct,imafolder,waterfolder); % GO 02/05/2017 + MRS_struct = SiemensDICOMRead(MRS_struct,metabfile{ii},waterfile{ii}); % GO 02/05/2017 WaterData = MRS_struct.fids.data_water; else - % Same as above, but without parsing the waterfolder. % GO 02/05/2017 - [imafolder,~,~] = fileparts(metabfile{ii}); % GO 11/01/2016 - MRS_struct = SiemensDICOMRead(MRS_struct,imafolder); % GO 11/01/2016 + MRS_struct = SiemensDICOMRead(MRS_struct,metabfile{ii}); % GO 11/01/2016 end FullData = MRS_struct.fids.data; % Determine order of ON and OFF acquisitions @@ -154,17 +147,10 @@ case 'dicom' % GO 11/30/2016 if exist('waterfile','var') - % Load the data. GannetLoad specifies a file in the first - % place, so take apart that filename, and feed the containing - % folder into the DICOMRead function. % GO 11/01/2016 - [dcmfolder,~,~] = fileparts(metabfile{ii}); % GO 02/05/2017 - [waterfolder,~,~] = fileparts(waterfile{ii}); % GO 02/05/2017 - MRS_struct = DICOMRead(MRS_struct,dcmfolder,waterfolder); % GO 02/05/2017 + MRS_struct = DICOMRead(MRS_struct,metabfile{ii},waterfile{ii}); % GO 02/05/2017 WaterData = MRS_struct.fids.data_water; else - % Same as above, but without parsing the waterfolder. % GO 02/05/2017 - [dcmfolder,~,~] = fileparts(metabfile{ii}); % GO 11/01/2016 - MRS_struct = DICOMRead(MRS_struct,dcmfolder); % GO 11/01/2016 + MRS_struct = DICOMRead(MRS_struct,metabfile{ii}); % GO 11/01/2016 end FullData = MRS_struct.fids.data; diff --git a/SiemensDICOMRead.m b/SiemensDICOMRead.m index 67a986a..95e46e8 100644 --- a/SiemensDICOMRead.m +++ b/SiemensDICOMRead.m @@ -1,19 +1,17 @@ -function MRS_struct = SiemensDICOMRead(MRS_struct,folder,waterfolder) -%% MRS_struct = SiemensDICOMRead(MRS_struct,folder,waterfolder) +function MRS_struct = SiemensDICOMRead(MRS_struct,metabfile,waterfile) +%% MRS_struct = SiemensDICOMRead(MRS_struct,metabfile,waterfile) % This function is designed to load edited MR spectroscopy data in the % Siemens flavour of DICOM data into a Gannet file structure. Files usually % have the extension '.IMA' or '.ima', and contain exactly 1 FID per % file, i.e. an acquisition of 320 averages will yield 320 IMA files. % -% The user must specify the folder containing all of these averages. It -% is assumed that they are ordered in the order of acquisition. +% It is assumed that they are ordered in the order of acquisition. % Water-suppressed and water-unsuppressed files need to be stored in -% separate folders, which need to be defined accordingly: +% separate folders, e.g. '/user/data/subject01/ima_gaba/' and +% '/user/data/subject01/ima_water/', respectively. % % Example: -% folder = '/user/data/subject01/ima_gaba/'; -% waterfolder = '/user/data/subject01/ima_water/'; (optional) -% MRS_struct = SiemensDICOMRead(MRS_struct,folder,waterfolder); +% MRS_struct = SiemensDICOMRead(MRS_struct,'/user/data/subject01/ima_gaba/metab.ima','/user/data/subject01/ima_water/water.ima'); % % Author: % Dr. Georg Oeltzschner (Johns Hopkins University, 2016-11-10) @@ -32,6 +30,8 @@ % of Minnesota) (2017-11-20). Thanks to Jim Lagopoulos. % 0.95: Fills missing voxel geometry parameters in DICOM header with zero % values. Thanks to Alen Tersakyan. +% 0.96: Fixed to accomodate batch processing of coregister/segmentation. +% (2018-09-19) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -46,6 +46,7 @@ % no overlap. % ima_file_list = [dir([folder,'/*.IMA']); dir([folder,'/*.ima'])]; % may % cause problems on win/unix systems, take out for now % GO 11/16/2016 +[folder,~,~] = fileparts(metabfile); % GO 11/01/2016 ima_file_list = dir([folder,'/*.IMA']); % GO 11/16/2016 fprintf('%d water-suppressed IMA files detected in %s.\n',length(ima_file_list),folder); @@ -60,7 +61,7 @@ %%% /PREPARATION %%% %%% HEADER INFO PARSING %%% -DicomHeader = read_dcm_header(ima_file_names{1}); +DicomHeader = read_dcm_header(metabfile); MRS_struct.p.seq = DicomHeader.sequenceFileName; MRS_struct.p.TR(ii) = DicomHeader.TR; MRS_struct.p.TE(ii) = DicomHeader.TE; @@ -106,12 +107,14 @@ % It appears that IMA stores the transients weirdly, 1-n/2 are all ONs, and % n/2-n are all OFFS. Shuffle them below. -a = MRS_struct.fids.data(:,1:end/2); -b = MRS_struct.fids.data(:,1+end/2:end); -c = zeros(size(MRS_struct.fids.data)); -c(:,1:2:end) = a; -c(:,2:2:end) = b; -MRS_struct.fids.data = c; +if size(MRS_struct.fids.data,2) > 1 + a = MRS_struct.fids.data(:,1:end/2); + b = MRS_struct.fids.data(:,1+end/2:end); + c = zeros(size(MRS_struct.fids.data)); + c(:,1:2:end) = a; + c(:,2:2:end) = b; + MRS_struct.fids.data = c; +end %%% /DATA LOADING %%% @@ -123,6 +126,7 @@ % Set up the file name array. if nargin == 3 + [waterfolder,~,~] = fileparts(waterfile); water_file_list = dir([waterfolder,'/*.IMA']); fprintf('%d water-unsuppressed IMA files detected in %s.\n',length(water_file_list),waterfolder); disp('Reading water-unsuppressed files...')