From 1ba4049cd23722a4bd4d7216741342db41083688 Mon Sep 17 00:00:00 2001 From: markmikkelsen Date: Thu, 20 Oct 2022 15:08:44 -0400 Subject: [PATCH] Several minor changes - Updated README.md - Coil combination of TWIX data now performed using GLS - Added support for new JHU universal sequence for Siemens data - Turned off warnings of iteration limits in SpectralRegistrationHERMES.m - Some cosmetic changes --- GERead.m | 60 +++++------ GannetPreInitialise.m | 2 +- README.md | 2 +- SiemensTwixRead.m | 192 +++++++++++++---------------------- SpectralRegistrationHERMES.m | 4 + 5 files changed, 109 insertions(+), 151 deletions(-) diff --git a/GERead.m b/GERead.m index 2bae64a..c8b6404 100644 --- a/GERead.m +++ b/GERead.m @@ -300,7 +300,7 @@ MetabData = ShapeData(:,:,refframes+2:end,:) * mult / 2; totalframes = totalframes - (refframes + 1); - + MRS_struct.p.nrows(ii) = totalframes; MRS_struct.p.Navg(ii) = dataframes * nex; MRS_struct.p.Nwateravg(ii) = refframes * nex; @@ -335,20 +335,20 @@ ShapeData = reshape(raw_data, [2 MRS_struct.p.npoints(ii) totalframes nreceivers]); - [X1,X2] = ndgrid(1:refframes, 1:nechoes); - X1 = X1'; X1 = X1(:); - X2 = X2'; X2 = X2(:); - Y1 = (-1).^(MRS_struct.p.GE.noadd(ii) * (X1-1)); - Y1 = permute(repmat(Y1, [1 MRS_struct.p.npoints(ii) 2 nreceivers]), [3 2 1 4]); - Y2 = 1 + (totalframes/nechoes) * (X2-1) + X1; + [X1,X2] = ndgrid(1:refframes, 1:nechoes); + X1 = X1'; X1 = X1(:); + X2 = X2'; X2 = X2(:); + Y1 = (-1).^(MRS_struct.p.GE.noadd(ii) * (X1-1)); + Y1 = permute(repmat(Y1, [1 MRS_struct.p.npoints(ii) 2 nreceivers]), [3 2 1 4]); + Y2 = 1 + (totalframes/nechoes) * (X2-1) + X1; WaterData = Y1 .* ShapeData(:,:,Y2,:) * multw; - [X1,X2] = ndgrid(1:dataframes, 1:nechoes); - X1 = X1'; X1 = X1(:); - X2 = X2'; X2 = X2(:); - Y1 = (-1).^(MRS_struct.p.GE.noadd(ii) * (X1-1)); - Y1 = permute(repmat(Y1, [1 MRS_struct.p.npoints(ii) 2 nreceivers]), [3 2 1 4]); - Y2 = 1 + refframes + (totalframes/nechoes) * (X2-1) + X1; + [X1,X2] = ndgrid(1:dataframes, 1:nechoes); + X1 = X1'; X1 = X1(:); + X2 = X2'; X2 = X2(:); + Y1 = (-1).^(MRS_struct.p.GE.noadd(ii) * (X1-1)); + Y1 = permute(repmat(Y1, [1 MRS_struct.p.npoints(ii) 2 nreceivers]), [3 2 1 4]); + Y2 = 1 + refframes + (totalframes/nechoes) * (X2-1) + X1; MetabData = Y1 .* ShapeData(:,:,Y2,:) * mult; totalframes = totalframes - (refframes + 1) * nechoes; % RTN 2017 @@ -362,30 +362,32 @@ WaterData = permute(WaterData, [3 1 2]); % Generalized least squares method (An et al., JMRI, 2013, doi:10.1002/jmri.23941) -[nCh, nPts, nReps] = size(WaterData); -noise_pts = false(1,nPts); +[nCh, nPts, nReps] = size(WaterData); +noise_pts = false(1,nPts); noise_pts(ceil(0.75*nPts):end) = true; -noise_pts = repmat(noise_pts, [1 nReps]); -tmpWaterData = reshape(WaterData, [nCh nPts*nReps]); +noise_pts = repmat(noise_pts, [1 nReps]); +tmpWaterData = reshape(WaterData, [nCh nPts*nReps]); -e = tmpWaterData(:,noise_pts); -Psi = e*e'; +e = tmpWaterData(:,noise_pts); +Psi = e*e'; WaterData_avg = mean(WaterData,3); -S = WaterData_avg(:,1); -w = (S'*(Psi\S))^-1 * S' / Psi; -WaterData = w.' .* WaterData; +S = WaterData_avg(:,1); +w = (S'*(Psi\S))^-1 * S' / Psi; +WaterData = w.' .* WaterData; + MRS_struct.fids.data_water = mean(squeeze(sum(WaterData,1)),2); -[nCh, nPts, nReps] = size(MetabData); -noise_pts = false(1,nPts); +[nCh, nPts, nReps] = size(MetabData); +noise_pts = false(1,nPts); noise_pts(ceil(0.75*nPts):end) = true; -noise_pts = repmat(noise_pts, [1 nReps]); -tmpMetabData = reshape(MetabData, [nCh nPts*nReps]); +noise_pts = repmat(noise_pts, [1 nReps]); +tmpMetabData = reshape(MetabData, [nCh nPts*nReps]); -e = tmpMetabData(:,noise_pts); -Psi = e*e'; -w = (S'*(Psi\S))^-1 * S' / Psi; +e = tmpMetabData(:,noise_pts); +Psi = e*e'; +w = (S'*(Psi\S))^-1 * S' / Psi; MetabData = w.' .* MetabData; + MRS_struct.fids.data = squeeze(sum(MetabData,1)); end diff --git a/GannetPreInitialise.m b/GannetPreInitialise.m index f3e733f..8a38910 100644 --- a/GannetPreInitialise.m +++ b/GannetPreInitialise.m @@ -38,7 +38,7 @@ MRS_struct.p.mat = 0; % Save MRS_struct as a .mat file MRS_struct.p.csv = 0; % Extract useful data from MRS_struct and export them to a .csv file (applies to GannetFit, % GannetSegment and GannetQuantify) - MRS_struct.p.append = 0; % Append PDF outputs into one PDF (separately for each module) (requires export_fig in the Gannet + MRS_struct.p.append = 1; % Append PDF outputs into one PDF (separately for each module) (requires export_fig in the Gannet % folder to be added to the search path and GhostScript to be installed) MRS_struct.p.hide = 0; % Do not display output figures diff --git a/README.md b/README.md index cff511d..bf0b799 100644 --- a/README.md +++ b/README.md @@ -60,7 +60,7 @@ __It is highly recommended that you only add the main SPM12 folder (_spm12_) to ## Compatibility -Gannet is currently being developed in MATLAB R2021a in macOS 11.4 Big Sur. While reasonable effort is made to ensure legacy and cross-OS compatibility, an error-free user experience is not guaranteed. +Gannet is currently being developed in MATLAB R2022b in macOS 12.6 Monterey. While reasonable effort is made to ensure legacy and cross-OS compatibility, an error-free user experience is not guaranteed. ## Supported file formats diff --git a/SiemensTwixRead.m b/SiemensTwixRead.m index 494eb3a..c8fc6ae 100644 --- a/SiemensTwixRead.m +++ b/SiemensTwixRead.m @@ -39,37 +39,39 @@ % combination (currently under dev.) % 2021-03-30: Added support of CMRR PRESS sequence and Michael % Dacko's MEGA-PRESS sequence. +% 2022-10-13: Coil combination now performed using generalized least +% squares. +% 2022-10-20: Added support for new JHU universal sequence. ii = MRS_struct.ii; % Get the raw data and header info from the MEGA-PRESS files. -[MetabData, MetabHeader] = GetTwixData(fname); -% Populate MRS_struct with relevant info. -MRS_struct.p.pointsBeforeEcho = MetabHeader.pointsBeforeEcho; -MRS_struct.p.sw(ii) = 1/MetabHeader.dwellTime; -MRS_struct.p.LarmorFreq(ii) = MetabHeader.tx_freq; -MRS_struct.p.TR(ii) = MetabHeader.TR; -MRS_struct.p.TE(ii) = MetabHeader.TE; -MRS_struct.p.npoints(ii) = size(MetabData,2); -MRS_struct.p.nrows(ii) = size(MetabData,3); -MRS_struct.p.Navg(ii) = size(MetabData,3); -MRS_struct.p.VoI_InPlaneRot(ii) = MetabHeader.VoI_InPlaneRot; -MRS_struct.p.NormCor(ii) = MetabHeader.NormCor; -MRS_struct.p.NormSag(ii) = MetabHeader.NormSag; -MRS_struct.p.NormTra(ii) = MetabHeader.NormTra; -MRS_struct.p.voxdim(ii,1) = MetabHeader.VoI_PeFOV; -MRS_struct.p.voxdim(ii,2) = MetabHeader.VoI_RoFOV; -MRS_struct.p.voxdim(ii,3) = MetabHeader.VoIThickness; -MRS_struct.p.voxoff(ii,1) = MetabHeader.PosSag; -MRS_struct.p.voxoff(ii,2) = MetabHeader.PosCor; -MRS_struct.p.voxoff(ii,3) = MetabHeader.PosTra; -MRS_struct.p.TablePosition(ii,1) = MetabHeader.TablePosSag; -MRS_struct.p.TablePosition(ii,2) = MetabHeader.TablePosCor; -MRS_struct.p.TablePosition(ii,3) = MetabHeader.TablePosTra; -MRS_struct.p.seqorig = MetabHeader.seqorig; +[MetabData, MetabHeader] = GetTwixData(fname); +MRS_struct.p.pointsBeforeEcho = MetabHeader.pointsBeforeEcho; +MRS_struct.p.sw(ii) = 1/MetabHeader.dwellTime; +MRS_struct.p.LarmorFreq(ii) = MetabHeader.tx_freq; +MRS_struct.p.TR(ii) = MetabHeader.TR; +MRS_struct.p.TE(ii) = MetabHeader.TE; +MRS_struct.p.npoints(ii) = size(MetabData,2); +MRS_struct.p.nrows(ii) = size(MetabData,3); +MRS_struct.p.Navg(ii) = size(MetabData,3); +MRS_struct.p.VoI_InPlaneRot(ii) = MetabHeader.VoI_InPlaneRot; +MRS_struct.p.NormCor(ii) = MetabHeader.NormCor; +MRS_struct.p.NormSag(ii) = MetabHeader.NormSag; +MRS_struct.p.NormTra(ii) = MetabHeader.NormTra; +MRS_struct.p.voxdim(ii,1) = MetabHeader.VoI_PeFOV; +MRS_struct.p.voxdim(ii,2) = MetabHeader.VoI_RoFOV; +MRS_struct.p.voxdim(ii,3) = MetabHeader.VoIThickness; +MRS_struct.p.voxoff(ii,1) = MetabHeader.PosSag; +MRS_struct.p.voxoff(ii,2) = MetabHeader.PosCor; +MRS_struct.p.voxoff(ii,3) = MetabHeader.PosTra; +MRS_struct.p.TablePosition(ii,1) = MetabHeader.TablePosSag; +MRS_struct.p.TablePosition(ii,2) = MetabHeader.TablePosCor; +MRS_struct.p.TablePosition(ii,3) = MetabHeader.TablePosTra; +MRS_struct.p.seqorig = MetabHeader.seqorig; if isfield(MetabHeader,'deltaFreq') - MRS_struct.p.Siemens.deltaFreq.metab(ii) = MetabHeader.deltaFreq; + MRS_struct.p.Siemens.deltaFreq.metab(ii) = MetabHeader.deltaFreq; end if isfield(MetabHeader,'editRF') @@ -89,16 +91,16 @@ % Undo phase cycling % Seems to be needed for some of Jamie Near's sequences if strcmp(MRS_struct.p.seqorig,'JN') - corrph = repmat([-1 1], [size(MetabData,2) size(MetabData,3)/2]); - corrph = repmat(corrph, [size(MetabData,1) 1 1]); - corrph = reshape(corrph, [size(MetabData,1) size(MetabData,2) size(MetabData,3)]); + corrph = repmat([-1 1], [size(MetabData,2) size(MetabData,3)/2]); + corrph = repmat(corrph, [size(MetabData,1) 1 1]); + corrph = reshape(corrph, [size(MetabData,1) size(MetabData,2) size(MetabData,3)]); MetabData = MetabData .* corrph; end % If water reference is provided, load this one as well, and populate % MRS_struct with water reference specific information. if nargin == 3 - [WaterData, WaterHeader] = GetTwixData(fname_water); + [WaterData, WaterHeader] = GetTwixData(fname_water); MRS_struct.p.pointsBeforeEcho_water = WaterHeader.pointsBeforeEcho; MRS_struct.p.sw_water(ii) = 1/WaterHeader.dwellTime; MRS_struct.p.TR_water(ii) = WaterHeader.TR; @@ -113,101 +115,48 @@ MRS_struct.p.Siemens = reorderstructure(MRS_struct.p.Siemens, 'editRF', 'deltaFreq'); end end - + % If additional data points have been acquired before the echo starts, % remove these here. WaterData = WaterData(:,(MRS_struct.p.pointsBeforeEcho_water+1):end,:); MRS_struct.p.npoints_water(ii) = MRS_struct.p.npoints_water(ii) - MRS_struct.p.pointsBeforeEcho_water; - + % Undo phase cycling % Seems to be needed for some of Jamie Near's sequences if strcmp(MRS_struct.p.seqorig,'JN') - corrph = repmat([-1 1], [size(WaterData,2) size(WaterData,3)/2]); - corrph = repmat(corrph, [size(WaterData,1) 1 1]); - corrph = reshape(corrph, [size(WaterData,1) size(WaterData,2) size(WaterData,3)]); + corrph = repmat([-1 1], [size(WaterData,2) size(WaterData,3)/2]); + corrph = repmat(corrph, [size(WaterData,1) 1 1]); + corrph = reshape(corrph, [size(WaterData,1) size(WaterData,2) size(WaterData,3)]); WaterData = WaterData .* corrph; end + + % Coil combination + % Generalized least squares method (An et al., JMRI, 2013, doi:10.1002/jmri.23941) + WaterData_avg = mean(WaterData,3); + e = WaterData_avg(:,ceil(0.75*size(WaterData_avg,2)):end); + Psi = e*e'; + [~,ind] = max(abs(WaterData_avg),[],2); + ind = mode(ind); + S = WaterData_avg(:,ind); + w = (S'*(Psi\S))^-1 * S' / Psi; + w = repmat(w.', [1 size(WaterData,2) size(WaterData,3)]); + WaterData = w .* WaterData; - % Coil combination and prephasing - -% firstpoint_water = conj(WaterData(:,1,:)); -% channels_scale = squeeze(sqrt(sum(firstpoint_water .* conj(firstpoint_water),1))); -% channels_scale = repmat(channels_scale, [1 size(WaterData,1) MRS_struct.p.npoints_water(ii)]); -% channels_scale = permute(channels_scale, [2 3 1]); -% firstpoint_water = repmat(firstpoint_water, [1 MRS_struct.p.npoints_water(ii) 1])./channels_scale; -% WaterData = WaterData .* firstpoint_water; -% WaterData = conj(squeeze(sum(WaterData,1))); -% WaterData = squeeze(mean(WaterData,2)); - - % Calculate mean of water FID for all channels - WaterMean = mean(WaterData,3); - % Determine signal strength for each channel - [~,ind] = max(abs(WaterMean),[],2); - ind = mode(ind); - sig = abs(WaterMean(:,ind)); - % Normalize so the sum is 1 - sig = sig/norm(sig); - % Determine phase of each channel - phi = angle(WaterMean(:,ind)); - % Apply changes - WaterData = WaterData .* repmat(exp(-1i*phi) .* sig, [1 size(WaterData,2) size(WaterData,3)]); - % Sum over coils - WaterData = squeeze(sum(conj(WaterData),1)); - % Average across averages - WaterData = squeeze(mean(WaterData,2)); - MRS_struct.fids.data_water = double(WaterData); - -% % Generalized least squares approach (MM: under dev.) -% % Align water signal over each avg. for each coil element -% lsqnonlinopts = optimoptions(@lsqnonlin); -% lsqnonlinopts = optimoptions(lsqnonlinopts,'Algorithm','levenberg-marquardt','Display','off'); -% -% t1 = 0:(1/MRS_struct.p.sw(ii)):(size(WaterData,2)-1)*(1/MRS_struct.p.sw(ii)); -% tMax = find(t1 <= 0.2,1,'last'); -% t2 = 0:(1/MRS_struct.p.sw(ii)):(tMax-1)*(1/MRS_struct.p.sw(ii)); -% -% for jj = 1:size(WaterData,1) -% flatdata(:,1,:) = real(WaterData(jj,1:tMax,:)); -% flatdata(:,2,:) = imag(WaterData(jj,1:tMax,:)); -% target = squeeze(flatdata(:,:,1)); -% target = target(:); -% for kk = 2:size(WaterData,3) -% transient = squeeze(flatdata(:,:,kk)); -% fun = @(x) SpecReg(double(transient(:)), double(target), t2, x); -% p = lsqnonlin(fun, [0 0], [], [], lsqnonlinopts); -% WaterData(jj,:,kk) = WaterData(jj,:,kk) .* ... -% exp(1i*p(1)*2*pi*t1) * exp(1i*pi/180*p(2)); -% end -% end -% -% WaterData_avg = mean(WaterData,3); -% e = WaterData_avg(:,ceil(0.75*size(WaterData_avg,2)):end); -% Psi = e*e'; -% [~,ind] = max(abs(WaterData_avg),[],2); -% ind = mode(ind); -% S = WaterData_avg(:,ind); -% w = (S'*(Psi\S))^-1 * S' / Psi; -% w = repmat(w.', [1 size(WaterData,2) size(WaterData,3)]); -% WaterData = w .* WaterData; -% MRS_struct.fids.data_water = double(mean(conj(squeeze(sum(WaterData,1))),2)); + MRS_struct.fids.data_water = double(mean(conj(squeeze(sum(WaterData,1))),2)); end -% Coil combination and prephasing of metabolite data. +% Coil combination of metabolite data % If water data has been acquired, use the phase information from it to -% phase the metabolite data. +% phase the metabolite data if nargin == 3 - MetabData = MetabData .* repmat(exp(-1i*phi) .* sig, [1 size(MetabData,2) size(MetabData,3)]); - MetabData = squeeze(sum(conj(MetabData),1)); - MRS_struct.fids.data = double(MetabData); - -% % Generalized least squares approach (MM: under dev.) -% MetabData_avg = mean(MetabData,3); -% e = MetabData_avg(:,ceil(0.75*size(MetabData_avg,2)):end); -% Psi = e*e'; -% w = (S'*(Psi\S))^-1 * S' / Psi; -% w = repmat(w.', [1 size(MetabData,2) size(MetabData,3)]); -% MetabData = w .* MetabData; -% MRS_struct.fids.data = double(conj(squeeze(sum(MetabData,1)))); + % Generalized least squares method (An et al., JMRI, 2013, doi:10.1002/jmri.23941) + MetabData_avg = mean(MetabData,3); + e = MetabData_avg(:,ceil(0.75*size(MetabData_avg,2)):end); + Psi = e*e'; + w = (S'*(Psi\S))^-1 * S' / Psi; + w = repmat(w.', [1 size(MetabData,2) size(MetabData,3)]); + MetabData = w .* MetabData; + MRS_struct.fids.data = double(conj(squeeze(sum(MetabData,1)))); else % If no water data provided, combine data based upon max point of % metabolite data (average all transients). @@ -254,7 +203,7 @@ % Read information from .hdr part of the TWIX object TwixHeader.readoutOSFactor = twix_obj.hdr.Config.ReadoutOSFactor; % Data are oversampled by this factor compared to exam card setting TwixHeader.removeOS = twix_obj.hdr.Config.RemoveOversampling; % Is the oversampling removed in the RDA files? -TwixHeader.TR = twix_obj.hdr.Config.TR * 1e-3; % TR [ms] +TwixHeader.TR = twix_obj.hdr.Config.TR(1) * 1e-3; % TR [ms] TwixHeader.vectorSize = twix_obj.hdr.Config.VectorSize; % Data points specified on exam card TwixHeader.VoI_InPlaneRot = twix_obj.hdr.Config.VoI_InPlaneRotAngle; % Voxel rotation in plane TwixHeader.VoI_RoFOV = twix_obj.hdr.Config.VoI_RoFOV; % Voxel size in readout direction [mm] @@ -274,7 +223,7 @@ % performed), the respective field is left empty in the TWIX file. This % case needs to be intercepted. Setting to the minimum possible value. VoI_Params = {'VoI_InPlaneRot','VoI_RoFOV','VoI_PeFOV','VoIThickness','NormCor','NormSag','NormTra', ... - 'PosCor','PosSag','PosTra','TablePosSag','TablePosCor','TablePosTra'}; + 'PosCor','PosSag','PosTra','TablePosSag','TablePosCor','TablePosTra'}; for pp = 1:length(VoI_Params) if isempty(TwixHeader.(VoI_Params{pp})) TwixHeader.(VoI_Params{pp}) = realmin('double'); @@ -354,6 +303,9 @@ elseif strfind(TwixHeader.sequenceFileName,'eja_svs_press') TwixHeader.seqtype = 'PRESS'; TwixHeader.seqorig = 'CMRR'; +elseif strfind(TwixHeader.sequenceFileName,'smm_svs_herc') + TwixHeader.seqtype = 'MEGA-PRESS'; + TwixHeader.seqorig = 'Universal'; else TwixHeader.seqorig = TwixHeader.sequenceString; error(['Unknown sequence: ' TwixHeader.seqorig '. Please contact the Gannet team for support.']) @@ -362,7 +314,7 @@ % Now reorder the FID data array according to software version and sequence % origin and sequence type. if strcmp(TwixHeader.seqtype,'PRESS') - + % For PRESS data, the first dimension of the 4D data array contains the % time-domain FID datapoints. The second dimension contains the number % of the coils. The third dimension contains the number of averages. @@ -376,16 +328,16 @@ if ndims(TwixData) == 4 TwixData = TwixData(:,:,:,2); end - + % For the standard Siemens svs_se sequence, the number of points % acquired before the echo maximum are stored here: TwixHeader.pointsBeforeEcho = twix_obj.image.freeParam(1); - + TwixData = permute(TwixData,[dims.coils dims.points dims.dyn dims.averages]); TwixData = reshape(TwixData,[size(TwixData,1) size(TwixData,2) size(TwixData,3)*size(TwixData,4)]); - + elseif any(strcmp(TwixHeader.seqtype,{'MEGA-PRESS','MEGA-sLASER'})) % SH 20191213 - + % For all known MEGA-PRESS implementations, the first dimension of the 4D % data array contains the time-domain FID datapoints. dims.points = 1; @@ -425,7 +377,7 @@ dims.dyn = find(strcmp(TwixHeader.sqzDims,'Ide')); end end - + % It looks like newer CMRR implementations may have another (5th) % dimension of the FID array: if strcmp(TwixHeader.seqorig,'CMRR') && length(TwixHeader.sqzDims) > 4 @@ -436,7 +388,7 @@ TwixData = permute(TwixData,[dims.coils dims.points dims.dyn dims.averages]); TwixData = reshape(TwixData,[size(TwixData,1), size(TwixData,2), size(TwixData,3) * size(TwixData,4)]); end - + % MEGA-PRESS sequences store the number of points acquired before the % echo maximum in different fields, depending on the origin of the % sequence: @@ -454,7 +406,7 @@ else TwixHeader.pointsBeforeEcho = twix_obj.image.freeParam(1); end - + end end diff --git a/SpectralRegistrationHERMES.m b/SpectralRegistrationHERMES.m index 54f136d..f09652a 100644 --- a/SpectralRegistrationHERMES.m +++ b/SpectralRegistrationHERMES.m @@ -2,6 +2,8 @@ % Align using a multistep variant of spectral registration (Mikkelsen et % al. Magn Reson Med. 2018;80(1):21-28. doi:10.1002/mrm.27027) +warning('off','stats:nlinfit:IterationLimitExceeded'); % temporarily suppress warning messages about iteration limit + showPlots = 0; % Looping parameters @@ -279,6 +281,8 @@ end fprintf('\n'); +warning('on','stats:nlinfit:IterationLimitExceeded'); % turn warning about about iteration limit back on + end