Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Update io_loadspec_twix #688

Merged
merged 1 commit into from
Jan 22, 2024
Merged

Update io_loadspec_twix #688

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
196 changes: 185 additions & 11 deletions libraries/FID-A/inputOutput/io_loadspec_twix.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,10 @@
%
% OUTPUTS:
% out = Input dataset in FID-A structure format.
% out_w = Input water reference data (only available for some
% sequences. This will be empty for others).

function out=io_loadspec_twix(filename);
function [out,out_w]=io_loadspec_twix(filename);


%read in the data using the new mapVBVD. This code has been adapted to
Expand Down Expand Up @@ -79,9 +81,10 @@
isWIP859=~isempty(strfind(sequence,'edit_859'));%Is this WIP 859 (MEGA-PRESS)?
isTLFrei=~isempty(strfind(sequence,'md_svs_edit')) ||... %Is Thomas Lange's MEGA-PRESS sequence
~isempty(strfind(sequence,'md_svs_slaser_edit')); %Is Thomas Lange's MEGA-s-LASER sequence
isMinn=~isempty(strfind(sequence,'eja_svs_')) ||... %Is this one of Eddie Auerbach's (CMRR, U Minnesota) sequences?
~isempty(strfind(sequence,'svs_slaser_dkd')) ||... % ... or Dinesh Deelchand's 2016 sLASER sequence?
~isempty(strfind(sequence,'svs_slaserVOI_dkd')); % ... or Dinesh Deelchand's 2022 'plug and play' sLASER sequence?
isMinn_eja=~isempty(strfind(sequence,'eja_svs_')); %Is this one of Eddie Auerbach's (CMRR, U Minnesota) sequences?
isMinn_dkd=~isempty(strfind(sequence,'svs_slaser_dkd')) ||... % ... or Dinesh Deelchand's 2016 sLASER sequence?
~isempty(strfind(sequence,'svs_slaserVOI_dkd'))||...% ... or Dinesh Deelchand's 2022 'plug and play' sLASER sequence?
~isempty(strfind(sequence,'svs_slaserVOI_dkd2'));%Is this Dinesh Deelchand's (CMRR, U Minnesota) sequence?
isSiemens=(~isempty(strfind(sequence,'svs_se')) ||... %Is this the Siemens PRESS seqeunce?
~isempty(strfind(sequence,'svs_st'))) && ... % or the Siemens STEAM sequence?
isempty(strfind(sequence,'eja_svs')); %And make sure it's not 'eja_svs_steam'.
Expand All @@ -105,7 +108,7 @@
end
end
end
elseif isMinn
elseif isMinn_eja || isMinn_dkd
if ~isempty(strfind(sequence,'mslaser'))
seq = 'MEGASLASER';
else
Expand All @@ -114,7 +117,7 @@
elseif ~isempty(strfind(sequence,'steam'))
seq = 'STEAM';
else
seq = 'MEGAPRESS';
seq = 'PRESS';
end
end
elseif isTLFrei && ~isempty(strfind(sequence,'md_svs_slaser_edit'))
Expand Down Expand Up @@ -225,11 +228,44 @@
Ncoils=twix_obj.hdr.Meas.iMaxNoOfRxChannels;

%Find the TE:
TE = twix_obj.hdr.MeasYaps.alTE{1}; %Franck Lamberton
if ~isMinn_dkd
TE = twix_obj.hdr.MeasYaps.alTE{1}; %Franck Lamberton
else
TE = twix_obj.hdr.MeasYaps.alTE{1} + twix_obj.hdr.MeasYaps.alTE{2} + twix_obj.hdr.MeasYaps.alTE{3};
end

%Find the TR:
TR = twix_obj.hdr.MeasYaps.alTR{1}; %Franck Lamberton

wRefs=false; %Flag to identify if there are automatic water reference scans acquired. There are none by default.

%Noticed that in Dinesh Deelchand's CMRR sLASER seuqence
%(svs_slaserVOI_dkd2) contains some reference scans, which are acquired at
%the very end. We will save these and store them as a separate water
%reference structure (out_w).
if isMinn_dkd
%If nRefs is non-zero, then there are automatic water reference scans.
%These will be saved as "outw".
nRefs=twix_obj.hdr.MeasYaps.sSpecPara.lAutoRefScanNo;
if nRefs==1
nRefs=0;%Becuase Nrefs seems to have a default value of 1 even if there are no reference scans.
end
if nRefs>0
wRefs=true;
end

if ndims(fids)==4
fids_w=fids(:,:,:,end-(nRefs-1):end);
fids=fids(:,:,:,end-(nRefs+Naverages-1):end-nRefs);
elseif ndims(fids)==3
fids_w=fids(:,:,end-(nRefs-1):end);
fids=fids(:,:,end-(nRefs+Naverages-1):end-nRefs);
elseif ndims(fids)==2
fids_w=fids(:,end-(nRefs-1):end);
fids=fids(:,end-(nRefs+Naverages-1):end-nRefs);
end
end

% Extract voxel dimensions
if (strcmp(version,'vd') || strcmp(version,'vb') || strcmp(version,'XA30'))
TwixHeader.VoI_RoFOV = twix_obj.hdr.Config.VoI_RoFOV; % Voxel size in readout direction [mm]
Expand Down Expand Up @@ -309,7 +345,7 @@

%Now index the dimension of the averages
if strcmp(version,'vd') || strcmp(version,'ve') || strcmp(version,'XA30')
if isMinn || isConnectom
if isMinn_eja || isMinn_dkd || isConnectom
dims.averages=find(strcmp(sqzDims,'Set'));
else
dims.averages=find(strcmp(sqzDims,'Ave'));
Expand Down Expand Up @@ -354,7 +390,7 @@
else
dims.subSpecs=find(strcmp(sqzDims,'Ida'));
end
elseif isWIP529 || isMinn
elseif isWIP529 || isMinn_eja || isMinn_dkd
dims.subSpecs=find(strcmp(sqzDims,'Eco'));
elseif ismodWIP
if strcmp(version,'vd') || strcmp(version,'ve') || strcmp(version,'vb')
Expand Down Expand Up @@ -403,45 +439,81 @@
if length(sqzDims)==5
fids=permute(fids,[dims.t dims.coils dims.averages dims.subSpecs dims.extras]);
dims.t=1;dims.coils=2;dims.averages=3;dims.subSpecs=4;dims.extras=5;
if wRefs
fids_w=permute(fids_w,[dims.t dims.coils dims.averages dims.subSpecs dims.extras]);
end
elseif length(sqzDims)==4
if dims.extras==0
fids=permute(fids,[dims.t dims.coils dims.averages dims.subSpecs]);
dims.t=1;dims.coils=2;dims.averages=3;dims.subSpecs=4;dims.extras=0;
if wRefs
fids_w=permute(fids_w,[dims.t dims.coils dims.averages dims.subSpecs]);
end
elseif dims.subSpecs==0
fids=permute(fids,[dims.t dims.coils dims.averages dims.extras]);
dims.t=1;dims.coils=2;dims.averages=3;dims.subSpecs=0;dims.extras=4;
if wRefs
fids_w=permute(fids_w,[dims.t dims.coils dims.averages dims.extras]);
end
elseif dims.averages==0
fids=permute(fids,[dims.t dims.coils dims.subSpecs dims.extras]);
dims.t=1;dims.coils=2;dims.averages=0;dims.subSpecs=3;dims.extras=4;
if wRefs
fids_w=permute(fids_w,[dims.t dims.coils dims.subSpecs dims.extras]);
end
elseif dims.coils==0
fids=permute(fids,[dims.t dims.averages dims.subSpecs dims.extras]);
dims.t=1;dims.coils=0;dims.averages=2;dims.subSpecs=3;dims.extras=4;
if wRefs
fids_w=permute(fids_w,[dims.t dims.averages dims.subSpecs dims.extras]);
end
end
elseif length(sqzDims)==3
if dims.extras==0 && dims.subSpecs==0
fids=permute(fids,[dims.t dims.coils dims.averages]);
dims.t=1;dims.coils=2;dims.averages=3;dims.subSpecs=0;dims.extras=0;
if wRefs
fids_w=permute(fids_w,[dims.t dims.coils dims.averages]);
end
elseif dims.extras==0 && dims.averages==0
fids=permute(fids,[dims.t dims.coils dims.subSpecs]);
dims.t=1;dims.coils=2;dims.averages=0;dims.subSpecs=3;dims.extras=0;
if wRefs
fids_w=permute(fids_w,[dims.t dims.coils dims.subSpecs]);
end
elseif dims.extras==0 && dims.coils==0
fids=permute(fids,[dims.t dims.averages dims.subSpecs]);
dims.t=1;dims.coils=0;dims.averages=2;dims.subSpecs=3;dims.extras=0;
if wRefs
fids_w=permute(fids_w,[dims.t dims.averages dims.subSpecs]);
end
end
elseif length(sqzDims)==2
if dims.extras==0 && dims.subSpecs==0 && dims.averages==0
fids=permute(fids,[dims.t dims.coils]);
dims.t=1;dims.coils=2;dims.averages=0;dims.subSpecs=0;dims.extras=0;
if wRefs
fids_w=permute(fids_w,[dims.t dims.coils]);
end
elseif dims.extras==0 && dims.subSpecs==0 && dims.coils==0
fids=permute(fids,[dims.t dims.averages]);
dims.t=1;dims.coils=0;dims.averages=2;dims.subSpecs=0;dims.extras=0;
if wRefs
fids_w=permute(fids_w,[dims.t dims.averages]);
end
elseif dims.extras==0 && dims.averages==0 && dims.coils==0
fids=permute(fids,[dims.t dims.subSpecs]);
dims.t=1;dims.coils=0;dims.averages=0;dims.subSpecs=2;dims.extras=0;
if wRefs
fids_w=permute(fids_w,[dims.t dims.subSpecs]);
end
end
elseif length(sqzDims)==1
fids=permute(fids,[dims.t]);
dims.t=1;dims.coils=0;dims.averages=0;dims.subSpecs=0;dims.extras=0;
if wRefs
fids_w=permute(fids_w,[dims.t]);
end
end

%Now reorder the fids for the Universal MEGA implementation
Expand Down Expand Up @@ -482,10 +554,15 @@

%Now get the size of the data array:
sz=size(fids);
if wRefs
sz_w=size(fids_w);
end

%Now take fft of time domain to get fid:
specs=fftshift(fft(fids,[],dims.t),dims.t);

if wRefs
specs_w=fftshift(ifft(fids_w,[],dims.t),dims.t);
end

%Now get relevant scan parameters:*****************************

Expand Down Expand Up @@ -524,17 +601,33 @@
if dims.averages~=0
averages=sz(dims.averages)*sz(dims.subSpecs);
rawAverages=averages;
if wRefs
averages_w=sz_w(dims.averages)*sz_w(dims.subSpecs);
rawAverages_w=averages_w;
end
else
averages=sz(dims.subSpecs);
rawAverages=1;
if wRefs
averages_w=sz_w(dims.subSpecs);
rawAverages_w=1;
end
end
else
if dims.averages~=0
averages=sz(dims.averages);
rawAverages=averages;
if wRefs
averages_w=sz_w(dims.averages);
rawAverages_w=averages_w;
end
else
averages=1;
rawAverages=1;
if wRefs
averages_w=1;
rawAverages_w=1;
end
end
end

Expand Down Expand Up @@ -568,7 +661,7 @@
else
leftshift = twix_obj.image.iceParam(5,1);
end
elseif isMinn || isConnectom || isDondersMRSfMRI
elseif isMinn_eja || isMinn_dkd || isConnectom || isDondersMRSfMRI
try
leftshift = twix_obj.image.iceParam(5,1);
catch
Expand Down Expand Up @@ -668,4 +761,85 @@
out.flags.isHERCULES = 1;
end

if wRefs
%FILLING IN DATA STRUCTURE
out_w.fids=fids_w;
out_w.specs=specs_w;
out_w.sz=sz_w;
out_w.ppm=ppm;
out_w.t=t;
out_w.spectralwidth=spectralwidth;
out_w.dwelltime=dwelltime;
out_w.txfrq=txfrq;
out_w.date=date;
out_w.dims=dims;
out_w.Bo=Bo;
out_w.averages=averages_w;
out_w.rawAverages=rawAverages_w;
out_w.subspecs=subspecs;
out_w.rawSubspecs=rawSubspecs;
out_w.seq=seq;
out_w.te=TE/1000;
out_w.tr=TR/1000;
out_w.pointsToLeftshift=leftshift;
out_w.centerFreq = centerFreq;
out_w.geometry = geometry;
if isfield(twix_obj.hdr.Config,'Nucleus')
out_w.nucleus = twix_obj.hdr.Config.Nucleus ;
end
if isfield(twix_obj.hdr.Dicom,'SoftwareVersions')
out_w.software = [version ' ' twix_obj.hdr.Dicom.SoftwareVersions];
else
out_w.software = version;
end


%FILLING IN THE FLAGS
out_w.flags.writtentostruct=1;
out_w.flags.gotparams=1;
out_w.flags.leftshifted=0;
out_w.flags.filtered=0;
out_w.flags.zeropadded=0;
out_w.flags.freqcorrected=0;
out_w.flags.phasecorrected=0;
out_w.flags.averaged=0;
out_w.flags.addedrcvrs=0;
out_w.flags.subtracted=0;
out_w.flags.writtentotext=0;
out_w.flags.downsampled=0;
if out_w.dims.subSpecs==0
out_w.flags.isFourSteps=0;
else
out_w.flags.isFourSteps=(out_w.sz(out_wop_pl.dims.subSpecs)==4);
end
% Add info for niiwrite
out_w.PatientPosition = twix_obj.hdr.Config.PatientPosition;
out_w.Manufacturer = 'Siemens';
[~,filename,ext] = fileparts(filename);
out_w.OriginalFile = [filename ext];
% Sequence flags
out_w.flags.isUnEdited = 0;
out_w.flags.isMEGA = 0;
out_w.flags.isHERMES = 0;
out_w.flags.isHERCULES = 0;
out_w.flags.isPRIAM = 0;
out_w.flags.isMRSI = 0;
if strcmp(seq,'PRESS') || strcmp(seq,'STEAM') || strcmp(seq,'SLASER')
out_w.flags.isUnEdited = 1;
end
if contains(seq,'MEGA')
out_w.flags.isMEGA = 1;
end
if strcmp(seq,'HERMES')
out_w.flags.isHERMES = 1;
end
if strcmp(seq,'HERCULES')
out_w.flags.isHERCULES = 1;
end
else
%No water reference data found. Returning empty struct for out_w:
disp('No water reference data found. Returning empty field');
out_w=[];
end

%DONE