diff --git a/libraries/FID-A/inputOutput/io_loadspec_rda.m b/libraries/FID-A/inputOutput/io_loadspec_rda.m index a0643c1c..1556b304 100644 --- a/libraries/FID-A/inputOutput/io_loadspec_rda.m +++ b/libraries/FID-A/inputOutput/io_loadspec_rda.m @@ -1,422 +1,422 @@ -%io_loadspec_rda.m -% -%Jay Hennessy, McGill University 2016. -%Heiner N. Raum, UKM Muenster 2024 (adaption of geometry header fields from XA30 to VE). -% -% USAGE: -% [out]=io_loadspec_rda(rda_filename); -% -% DESCRIPTION: -% Reads in siemens rda data (.rda file). -% -% op_loadspec_rda outputs the data in structure format, with fields corresponding to time -% scale, fids, frequency scale, spectra, and header fields containing -% information about the acquisition. The resulting matlab structure can be -% operated on by the other functions in this MRS toolbox. -% -% INPUTS: -% pathname = pathname of Siemens rda data to load. Can be folder or file, -% which will autmatically be distinguished. -% -% OUTPUTS: -% out = Input dataset in FID-A structure format. - -function [out] = io_loadspec_rda(pathname) - - if isfile(pathname) - filesInFolder{1} = pathname; - else - % Create list of complete filenames (incl. path) in the folder - dirFolder = dir(pathname); - filesInFolder = dirFolder(~[dirFolder.isdir]); - filesInFolder = filesInFolder(~ismember({filesInFolder.name}, {'.','..','.DS_Store'})); - %hidden = logical(ones(1,length(filesInFolder))); - hidden = true(length(filesInFolder)); - for jj = 1:length(filesInFolder) - if strcmp(filesInFolder(jj).name(1),'.') - hidden(jj) = 0; - end - end - filesInFolder = filesInFolder(hidden);%delete hidden files - filesInFolder = fullfile(pathname, {filesInFolder.name}); - end - - fid = fopen(filesInFolder{1}); - - head_start_text = '>>> Begin of header <<<'; - head_end_text = '>>> End of header <<<'; - - tline = fgets(fid); - - while ~contains(tline,head_end_text) - %while (isempty(strfind(tline , head_end_text))) - - tline = fgets(fid); - - if ~contains(tline,head_start_text) && ~contains (tline,head_end_text) - %if ( isempty(strfind (tline , head_start_text)) + isempty(strfind (tline , head_end_text )) == 2) - - % Store this data in the appropriate format - - occurence_of_colon = findstr(':',tline); - variable = tline(1:occurence_of_colon-1) ; - value = tline(occurence_of_colon+1 : length(tline)) ; - - switch variable - case { 'PatientID' , 'PatientName' , 'StudyDescription' , 'PatientBirthDate' , 'StudyDate' , 'StudyTime' , 'PatientAge' , 'SeriesDate' , ... - 'SeriesTime' , 'SeriesDescription' , 'ProtocolName' , 'PatientPosition' , 'ModelName' , 'StationName' , 'InstitutionName' , ... - 'DeviceSerialNumber', 'InstanceDate' , 'InstanceTime' , 'InstanceComments' , 'SequenceName' , 'SequenceDescription' , 'Nucleus' ,... - 'TransmitCoil' } - eval(['rda.' , variable , ' = value; ']); - case { 'PatientSex' } - % Sex converter! (int to M,F,U) - switch value - case 0 - rda.sex = 'Unknown'; - case 1 - rda.sex = 'Male'; - case 2 - - rda.sex = 'Female'; - end - - case { 'SeriesNumber' , 'InstanceNumber' , 'AcquisitionNumber' , 'NumOfPhaseEncodingSteps' , 'NumberOfRows' , 'NumberOfColumns' , 'VectorSize' } - %Integers - eval(['rda.' , variable , ' = str2num(value); ']); - case { 'PatientWeight' , 'TR' , 'TE' , 'TM' , 'DwellTime' , 'NumberOfAverages' , 'MRFrequency' , 'MagneticFieldStrength' , 'FlipAngle' , ... - 'SliceThickness' , 'FoVHeight' , 'FoVWidth' , 'PercentOfRectFoV' , 'PixelSpacingRow' , 'PixelSpacingCol', ... - 'VOIRotationInPlane', 'VOINormalTra', 'VOINormalCor', 'VOINormalSag', 'VOIReadoutFOV', 'VOIPhaseFOV', 'VOIThickness', ... - 'VOIPositionTra', 'VOIPositionCor', 'VOIPositionSag'} - %Floats - eval(['rda.' , variable , ' = str2num(value); ']); - case {'SoftwareVersion[0]' } - rda.software_version = value; - case {'CSIMatrixSize[0]' } - rda.CSIMatrix_Size(1) = str2num(value); - case {'CSIMatrixSize[1]' } - rda.CSIMatrix_Size(2) = str2num(value); - case {'CSIMatrixSize[2]' } - rda.CSIMatrix_Size(3) = str2num(value); - case {'PositionVector[0]' } - rda.PositionVector(1) = str2num(value); - case {'PositionVector[1]' } - rda.PositionVector(2) = str2num(value); - case {'PositionVector[2]' } - rda.PositionVector(3) = str2num(value); - case {'RowVector[0]' } - rda.RowVector(1) = str2num(value); - case {'RowVector[1]' } - rda.RowVector(2) = str2num(value); - case {'RowVector[2]' } - rda.RowVector(3) = str2num(value); - case {'ColumnVector[0]' } - rda.ColumnVector(1) = str2num(value); - case {'ColumnVector[1]' } - rda.ColumnVector(2) = str2num(value); - case {'ColumnVector[2]' } - rda.ColumnVector(3) = str2num(value); - %Modified HR% - case {'SlabThickness[0]' } - rda.SlabThickness(1) = str2num(value); - case {'SlabThickness[1]' } - rda.SlabThickness(2) = str2num(value); - case {'SlabThickness[2]' } - rda.SlabThickness(3) = str2num(value); - case {'SlabOrientation[0,0]' } - rda.SlabOrientation(1,1) = str2num(value); - case {'SlabOrientation[0,1]' } - rda.SlabOrientation(1,2) = str2num(value); - case {'SlabOrientation[0,2]' } - rda.SlabOrientation(1,3) = str2num(value); - case {'SlabOrientation[1,0]' } - rda.SlabOrientation(2,1) = str2num(value); - case {'SlabOrientation[1,1]' } - rda.SlabOrientation(2,2) = str2num(value); - case {'SlabOrientation[1,2]' } - rda.SlabOrientation(2,3) = str2num(value); - case {'SlabOrientation[2,0]' } - rda.SlabOrientation(3,1) = str2num(value); - case {'SlabOrientation[2,1]' } - rda.SlabOrientation(3,2) = str2num(value); - case {'SlabOrientation[2,2]' } - rda.SlabOrientation(3,3) = str2num(value); - %mod. end% - otherwise - % We don't know what this variable is. Report this just to keep things clear - %disp(['Unrecognised variable ' , variable ]); - end - - else - % Don't bother storing this bit of the output - end - - end - %% - - % Get the header of the first file to make some decisions. - seqtype = rda.SeriesDescription; - - - % Prepare voxel geometry information - % If a parameter is set to zero (e.g. if no voxel rotation is - % 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 = {'VOIRotationInPlane', 'VOINormalTra', 'VOINormalCor', 'VOINormalSag', 'VOIReadoutFOV', 'VOIPhaseFOV', 'VOIThickness', ... - 'VOIPositionTra', 'VOIPositionCor', 'VOIPositionSag'}; - for pp = 1:length(VoI_Params) - if ~isfield(rda, VoI_Params{pp}) - rda.(VoI_Params{pp}) = realmin('double'); - %Modified HR% - flag.(VoI_Params{pp}) = false; - % - else - flag.(VoI_Params{pp}) = true; - end - end - - % Write them - %Modified HR% - if contains(rda.software_version,'XA30','IgnoreCase',true) && ~any(structfun(@(x) x == true,flag)) - geometry.size.VoI_RoFOV = rda.SlabThickness(3); % Voxel size in readout direction [mm] - geometry.size.VoI_PeFOV = rda.SlabThickness(2); % Voxel size in phase encoding direction [mm] - geometry.size.VoIThickness = rda.SlabThickness(1); % Voxel size in slice selection direction [mm] - geometry.pos.PosSag = rda.PositionVector(1); % Sagittal coordinate of voxel [mm] - geometry.pos.PosCor = rda.PositionVector(2); % Coronal coordinate of voxel [mm] - geometry.pos.PosTra = rda.PositionVector(3); % Transversal coordinate of voxel [mm] - geometry.rot.NormSag = rda.SlabOrientation(1,1); % Sagittal component of normal vector of voxel - geometry.rot.NormCor = rda.SlabOrientation(1,2); % Coronal component of normal vector of voxel - geometry.rot.NormTra = rda.SlabOrientation(1,3); % Transversal component of normal vector of voxel - - %Calculate InPlaneRot (HR), adapted from coreg_siemens.m; still in - %DICOM coordinates - Norm = rda.SlabOrientation(1,:)'; - [~, maxdir] = max([abs(Norm(1)) abs(Norm(2)) abs(Norm(3))]); - switch maxdir - case 1 - vox_orient = 's'; % 't' = transversal, 's' = sagittal', 'c' = coronal; - case 2 - vox_orient = 'c'; % 't' = transversal, 's' = sagittal', 'c' = coronal; - case 3 - vox_orient = 't'; % 't' = transversal, 's' = sagittal', 'c' = coronal; - end - Phase = zeros(3, 1); - switch vox_orient - case 't' - % For transversal voxel orientation, the phase reference vector lies in - % the sagittal plane - Phase(1) = 0; - Phase(2) = Norm(3)*sqrt(1/(Norm(2)*Norm(2)+Norm(3)*Norm(3))); - Phase(3) = -Norm(2)*sqrt(1/(Norm(2)*Norm(2)+Norm(3)*Norm(3))); - VoxDims = [geometry.size.VoI_PeFOV geometry.size.VoI_RoFOV geometry.size.VoIThickness]; - case 'c' - % For coronal voxel orientation, the phase reference vector lies in - % the transversal plane - Phase(1) = Norm(2)*sqrt(1/(Norm(1)*Norm(1)+Norm(2)*Norm(2))); - Phase(2) = -Norm(1)*sqrt(1/(Norm(1)*Norm(1)+Norm(2)*Norm(2))); - Phase(3) = 0; - VoxDims = [geometry.size.VoI_PeFOV geometry.size.VoI_RoFOV geometry.size.VoIThickness]; - case 's' - % For sagittal voxel orientation, the phase reference vector lies in - % the transversal plane - Phase(1) = -Norm(2)*sqrt(1/(Norm(1)*Norm(1)+Norm(2)*Norm(2))); - Phase(2) = Norm(1)*sqrt(1/(Norm(1)*Norm(1)+Norm(2)*Norm(2))); - Phase(3) = 0; - VoxDims = [geometry.size.VoI_PeFOV geometry.size.VoI_RoFOV geometry.size.VoIThickness]; - end - PhaseRot = rda.SlabOrientation(2,:)'; - % swapped signs because VoI_InPlaneRot rotates counterclockwise - % around Norm - if dot(cross(Phase,PhaseRot),Norm) <= 0 - geometry.rot.VoI_InPlaneRot = acos(dot(Phase,PhaseRot)); - else - geometry.rot.VoI_InPlaneRot = -acos(dot(Phase,PhaseRot)); - end - %End of calculating InPlaneRot - else - geometry.size.VoI_RoFOV = rda.VOIReadoutFOV; % Voxel size in readout direction [mm] - geometry.size.VoI_PeFOV = rda.VOIPhaseFOV; % Voxel size in phase encoding direction [mm] - geometry.size.VoIThickness = rda.VOIThickness; % Voxel size in slice selection direction [mm] - geometry.pos.PosCor = rda.VOIPositionCor; % Coronal coordinate of voxel [mm] - geometry.pos.PosSag = rda.VOIPositionSag; % Sagittal coordinate of voxel [mm] - geometry.pos.PosTra = rda.VOIPositionTra; % Transversal coordinate of voxel [mm] - geometry.rot.VoI_InPlaneRot = rda.VOIRotationInPlane; % Voxel rotation in plane - geometry.rot.NormCor = rda.VOINormalCor; % Coronal component of normal vector of voxel - geometry.rot.NormSag = rda.VOINormalSag; % Sagittal component of normal vector of voxel - geometry.rot.NormTra = rda.VOINormalTra; % Transversal component of normal vector of voxel - end - %mod. end% - out.geometry = geometry; - - - % - % So now we should have got to the point after the header text - % - % Siemens documentation suggests that the data should be in a double complex format (8bytes for real, and 8 for imaginary?) - % - fclose(fid); - bytes_per_point = 16; - - % Preallocate array in which the FIDs are to be extracted. - fids = zeros(length(filesInFolder),rda.VectorSize); - % Collect all FIDs and sort them into fids array - for kk = 1:length(filesInFolder) - fid = fopen(filesInFolder{kk}); - tline = fgets(fid); - %while (isempty(strfind(tline , head_end_text))) - while ~contains(tline , head_end_text) - tline = fgets(fid); - end - complex_data = fread(fid , rda.CSIMatrix_Size(1) * rda.CSIMatrix_Size(1) *rda.CSIMatrix_Size(1) *rda.VectorSize * 2 , 'double'); - %fread(fid , 1, 'double'); %This was a check to confirm that we had read all the data (it passed!) - fclose(fid); - - % Now convert this data into something meaningful - - %Reshape so that we can get the real and imaginary separated - hmm = reshape(complex_data, 2 , rda.VectorSize , rda.CSIMatrix_Size(1) , rda.CSIMatrix_Size(2) , rda.CSIMatrix_Size(3) ); - - %Combine the real and imaginary into the complex matrix - fids(kk,:) = complex(hmm(1,:,:,:,:),hmm(2,:,:,:,:)); - end - fids = fids'; - %Remove the redundant first element in the array - % Time_domain_data = reshape(hmm_complex, rda.VectorSize , rda.CSIMatrix_Size(1) , rda.CSIMatrix_Size(2) , rda.CSIMatrix_Size(3)); - % Time_domain_data = fft(fft(fftshift(fftshift(padarray(fftshift(fftshift(ifft(ifft(Time_domain_data,[],2),[],3),2),3),[0 16 16],'both'),2),3),[],2),[],3); - - % rewind first point phase - corrph = conj(fids(1,:))./abs(fids(1,:)); - corrph = repmat(corrph, [size(fids, 1), 1]); - fids = fids .* corrph; - - % get the spectrum from the fid - specs = fftshift(fft(fids,[],1),1); - - % make calculations for the output mrs structure - sz = size(fids); - dwelltime = rda.DwellTime/1000000; - spectralwidth=1/dwelltime; - txfrq = rda.MRFrequency*1000000; - dims.t = 1; - dims.subSpecs = 0; - dims.coils = 0; - dims.extras = 0; - dims.averages = 2; - % Bo = rda.MagneticFieldStrength; - Bo = rda.MRFrequency/42.577; - rawAverages = rda.NumberOfAverages; - if length(filesInFolder) >= rawAverages - if contains(seqtype, 'edit') - rawAverages = length(filesInFolder); - averages = rda.NumberOfAverages; - subspecs = fix(rawAverages/averages); - rawSubspecs = rda.NumberOfAverages; - else - rawAverages = length(filesInFolder); - averages = rawAverages; - subspecs = 1; - rawSubspecs = 1; - end - elseif length(filesInFolder) == 1 - averages = 1; - subspecs =1; - rawSubspecs = 'na'; - else - if contains(seqtype,'edit') || contains(seqtype,'mpress') - subspecs = 2; - rawSubspecs = 2; - % Distinguish whether MEGA data have been averaged on-scanner or - % not - if length(filesInFolder) == 2 - dims.subSpecs = 3; - dims.averages = 2; - averages = 1; - rawAverages = 1; - fids = reshape(fids,[sz(1),averages,2 ]); - specs = reshape(specs,[sz(1),averages,2 ]); - else - dims.subSpecs = 3; - dims.averages = 2; - averages = rawAverages/2; - fids = reshape(fids,[sz(1),averages,2 ]); - specs = reshape(specs,[sz(1),averages,2 ]); - end - sz = size(fids); - end - end - - date = rda.StudyDate; - seq = rda.SequenceDescription; - TE = rda.TE; - TR = rda.TR; - pointsToLeftShift = 0; - - %Calculate t and ppm arrays using the calculated parameters: - f=[(-spectralwidth/2)+(spectralwidth/(2*sz(1))):spectralwidth/(sz(1)):(spectralwidth/2)-(spectralwidth/(2*sz(1)))]; - ppm=f/(Bo*42.577); - % Siemens data assumes the center frequency to be 4.65 ppm: - centerFreq = 4.65; - ppm=ppm + centerFreq; - - t=[0:dwelltime:(sz(1)-1)*dwelltime]; - - %FILLING IN DATA STRUCTURE - out.fids=fids; - out.specs=specs; - out.sz=sz; - out.ppm=ppm; - out.t=t; - out.spectralwidth=spectralwidth; - out.dwelltime=dwelltime; - out.txfrq=txfrq; - out.date=date; - out.dims=dims; - out.Bo=Bo; - out.averages=averages; - out.rawAverages=rawAverages; - out.subspecs=subspecs; - out.rawSubspecs=rawSubspecs; - out.seq=seq; - out.te=TE; - out.tr=TR; - out.pointsToLeftshift=pointsToLeftShift; - out.centerFreq = centerFreq; - out.geometry = geometry; - out.software = rda.software_version; - - %FILLING IN THE FLAGS - out.flags.writtentostruct=1; - out.flags.gotparams=1; - out.flags.leftshifted=0; - out.flags.filtered=0; - out.flags.zeropadded=0; - out.flags.freqcorrected=0; - out.flags.phasecorrected=0; - if out.averages == 1 - out.flags.averaged=1; - else - out.flags.averaged=0; - end - out.flags.addedrcvrs=1; - out.flags.subtracted=0; - out.flags.writtentotext=0; - out.flags.downsampled=0; - if out.dims.subSpecs==0 - out.flags.isFourSteps=0; - else - out.flags.isFourSteps=(out.sz(out.dims.subSpecs)==4); - end - - % Add info for niiwrite - out.PatientPosition = rda.PatientPosition; - out.Manufacturer = 'Siemens'; - [~,filename,ext] = fileparts(filesInFolder{1}); - out.OriginalFile = [filename ext]; - - % Sequence flags - out.flags.isUnEdited = 0; - out.flags.isMEGA = 0; - out.flags.isHERMES = 0; - out.flags.isHERCULES = 0; - out.flags.isPRIAM = 0; - out.flags.isMRSI = 0; -end +%io_loadspec_rda.m +% +%Jay Hennessy, McGill University 2016. +%Heiner N. Raum, UKM Muenster 2024 (adaption of geometry header fields from XA30 to VE). +% +% USAGE: +% [out]=io_loadspec_rda(rda_filename); +% +% DESCRIPTION: +% Reads in siemens rda data (.rda file). +% +% op_loadspec_rda outputs the data in structure format, with fields corresponding to time +% scale, fids, frequency scale, spectra, and header fields containing +% information about the acquisition. The resulting matlab structure can be +% operated on by the other functions in this MRS toolbox. +% +% INPUTS: +% pathname = pathname of Siemens rda data to load. Can be folder or file, +% which will autmatically be distinguished. +% +% OUTPUTS: +% out = Input dataset in FID-A structure format. + +function [out] = io_loadspec_rda(pathname) + + if isfile(pathname) + filesInFolder{1} = pathname; + else + % Create list of complete filenames (incl. path) in the folder + dirFolder = dir(pathname); + filesInFolder = dirFolder(~[dirFolder.isdir]); + filesInFolder = filesInFolder(~ismember({filesInFolder.name}, {'.','..','.DS_Store'})); + %hidden = logical(ones(1,length(filesInFolder))); + hidden = true(length(filesInFolder)); + for jj = 1:length(filesInFolder) + if strcmp(filesInFolder(jj).name(1),'.') + hidden(jj) = 0; + end + end + filesInFolder = filesInFolder(hidden);%delete hidden files + filesInFolder = fullfile(pathname, {filesInFolder.name}); + end + + fid = fopen(filesInFolder{1}); + + head_start_text = '>>> Begin of header <<<'; + head_end_text = '>>> End of header <<<'; + + tline = fgets(fid); + + while ~contains(tline,head_end_text) + %while (isempty(strfind(tline , head_end_text))) + + tline = fgets(fid); + + if ~contains(tline,head_start_text) && ~contains (tline,head_end_text) + %if ( isempty(strfind (tline , head_start_text)) + isempty(strfind (tline , head_end_text )) == 2) + + % Store this data in the appropriate format + + occurence_of_colon = findstr(':',tline); + variable = tline(1:occurence_of_colon-1) ; + value = tline(occurence_of_colon+1 : length(tline)) ; + + switch variable + case { 'PatientID' , 'PatientName' , 'StudyDescription' , 'PatientBirthDate' , 'StudyDate' , 'StudyTime' , 'PatientAge' , 'SeriesDate' , ... + 'SeriesTime' , 'SeriesDescription' , 'ProtocolName' , 'PatientPosition' , 'ModelName' , 'StationName' , 'InstitutionName' , ... + 'DeviceSerialNumber', 'InstanceDate' , 'InstanceTime' , 'InstanceComments' , 'SequenceName' , 'SequenceDescription' , 'Nucleus' ,... + 'TransmitCoil' } + eval(['rda.' , variable , ' = value; ']); + case { 'PatientSex' } + % Sex converter! (int to M,F,U) + switch value + case 0 + rda.sex = 'Unknown'; + case 1 + rda.sex = 'Male'; + case 2 + + rda.sex = 'Female'; + end + + case { 'SeriesNumber' , 'InstanceNumber' , 'AcquisitionNumber' , 'NumOfPhaseEncodingSteps' , 'NumberOfRows' , 'NumberOfColumns' , 'VectorSize' } + %Integers + eval(['rda.' , variable , ' = str2num(value); ']); + case { 'PatientWeight' , 'TR' , 'TE' , 'TM' , 'DwellTime' , 'NumberOfAverages' , 'MRFrequency' , 'MagneticFieldStrength' , 'FlipAngle' , ... + 'SliceThickness' , 'FoVHeight' , 'FoVWidth' , 'PercentOfRectFoV' , 'PixelSpacingRow' , 'PixelSpacingCol', ... + 'VOIRotationInPlane', 'VOINormalTra', 'VOINormalCor', 'VOINormalSag', 'VOIReadoutFOV', 'VOIPhaseFOV', 'VOIThickness', ... + 'VOIPositionTra', 'VOIPositionCor', 'VOIPositionSag'} + %Floats + eval(['rda.' , variable , ' = str2num(value); ']); + case {'SoftwareVersion[0]' } + rda.software_version = value; + case {'CSIMatrixSize[0]' } + rda.CSIMatrix_Size(1) = str2num(value); + case {'CSIMatrixSize[1]' } + rda.CSIMatrix_Size(2) = str2num(value); + case {'CSIMatrixSize[2]' } + rda.CSIMatrix_Size(3) = str2num(value); + case {'PositionVector[0]' } + rda.PositionVector(1) = str2num(value); + case {'PositionVector[1]' } + rda.PositionVector(2) = str2num(value); + case {'PositionVector[2]' } + rda.PositionVector(3) = str2num(value); + case {'RowVector[0]' } + rda.RowVector(1) = str2num(value); + case {'RowVector[1]' } + rda.RowVector(2) = str2num(value); + case {'RowVector[2]' } + rda.RowVector(3) = str2num(value); + case {'ColumnVector[0]' } + rda.ColumnVector(1) = str2num(value); + case {'ColumnVector[1]' } + rda.ColumnVector(2) = str2num(value); + case {'ColumnVector[2]' } + rda.ColumnVector(3) = str2num(value); + %Modified HR% + case {'SlabThickness[0]' } + rda.SlabThickness(1) = str2num(value); + case {'SlabThickness[1]' } + rda.SlabThickness(2) = str2num(value); + case {'SlabThickness[2]' } + rda.SlabThickness(3) = str2num(value); + case {'SlabOrientation[0,0]' } + rda.SlabOrientation(1,1) = str2num(value); + case {'SlabOrientation[0,1]' } + rda.SlabOrientation(1,2) = str2num(value); + case {'SlabOrientation[0,2]' } + rda.SlabOrientation(1,3) = str2num(value); + case {'SlabOrientation[1,0]' } + rda.SlabOrientation(2,1) = str2num(value); + case {'SlabOrientation[1,1]' } + rda.SlabOrientation(2,2) = str2num(value); + case {'SlabOrientation[1,2]' } + rda.SlabOrientation(2,3) = str2num(value); + case {'SlabOrientation[2,0]' } + rda.SlabOrientation(3,1) = str2num(value); + case {'SlabOrientation[2,1]' } + rda.SlabOrientation(3,2) = str2num(value); + case {'SlabOrientation[2,2]' } + rda.SlabOrientation(3,3) = str2num(value); + %mod. end% + otherwise + % We don't know what this variable is. Report this just to keep things clear + %disp(['Unrecognised variable ' , variable ]); + end + + else + % Don't bother storing this bit of the output + end + + end + %% + + % Get the header of the first file to make some decisions. + seqtype = rda.SeriesDescription; + + + % Prepare voxel geometry information + % If a parameter is set to zero (e.g. if no voxel rotation is + % 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 = {'VOIRotationInPlane', 'VOINormalTra', 'VOINormalCor', 'VOINormalSag', 'VOIReadoutFOV', 'VOIPhaseFOV', 'VOIThickness', ... + 'VOIPositionTra', 'VOIPositionCor', 'VOIPositionSag'}; + for pp = 1:length(VoI_Params) + if ~isfield(rda, VoI_Params{pp}) + rda.(VoI_Params{pp}) = realmin('double'); + %Modified HR% + flag.(VoI_Params{pp}) = false; + % + else + flag.(VoI_Params{pp}) = true; + end + end + + % Write them + %Modified HR% + if contains(rda.software_version,{'XA30','XA60'},'IgnoreCase',true) && ~any(structfun(@(x) x == true,flag)) + geometry.size.VoI_RoFOV = rda.SlabThickness(3); % Voxel size in readout direction [mm] + geometry.size.VoI_PeFOV = rda.SlabThickness(2); % Voxel size in phase encoding direction [mm] + geometry.size.VoIThickness = rda.SlabThickness(1); % Voxel size in slice selection direction [mm] + geometry.pos.PosSag = rda.PositionVector(1); % Sagittal coordinate of voxel [mm] + geometry.pos.PosCor = rda.PositionVector(2); % Coronal coordinate of voxel [mm] + geometry.pos.PosTra = rda.PositionVector(3); % Transversal coordinate of voxel [mm] + geometry.rot.NormSag = rda.SlabOrientation(1,1); % Sagittal component of normal vector of voxel + geometry.rot.NormCor = rda.SlabOrientation(1,2); % Coronal component of normal vector of voxel + geometry.rot.NormTra = rda.SlabOrientation(1,3); % Transversal component of normal vector of voxel + + %Calculate InPlaneRot (HR), adapted from coreg_siemens.m; still in + %DICOM coordinates + Norm = rda.SlabOrientation(1,:)'; + [~, maxdir] = max([abs(Norm(1)) abs(Norm(2)) abs(Norm(3))]); + switch maxdir + case 1 + vox_orient = 's'; % 't' = transversal, 's' = sagittal', 'c' = coronal; + case 2 + vox_orient = 'c'; % 't' = transversal, 's' = sagittal', 'c' = coronal; + case 3 + vox_orient = 't'; % 't' = transversal, 's' = sagittal', 'c' = coronal; + end + Phase = zeros(3, 1); + switch vox_orient + case 't' + % For transversal voxel orientation, the phase reference vector lies in + % the sagittal plane + Phase(1) = 0; + Phase(2) = Norm(3)*sqrt(1/(Norm(2)*Norm(2)+Norm(3)*Norm(3))); + Phase(3) = -Norm(2)*sqrt(1/(Norm(2)*Norm(2)+Norm(3)*Norm(3))); + VoxDims = [geometry.size.VoI_PeFOV geometry.size.VoI_RoFOV geometry.size.VoIThickness]; + case 'c' + % For coronal voxel orientation, the phase reference vector lies in + % the transversal plane + Phase(1) = Norm(2)*sqrt(1/(Norm(1)*Norm(1)+Norm(2)*Norm(2))); + Phase(2) = -Norm(1)*sqrt(1/(Norm(1)*Norm(1)+Norm(2)*Norm(2))); + Phase(3) = 0; + VoxDims = [geometry.size.VoI_PeFOV geometry.size.VoI_RoFOV geometry.size.VoIThickness]; + case 's' + % For sagittal voxel orientation, the phase reference vector lies in + % the transversal plane + Phase(1) = -Norm(2)*sqrt(1/(Norm(1)*Norm(1)+Norm(2)*Norm(2))); + Phase(2) = Norm(1)*sqrt(1/(Norm(1)*Norm(1)+Norm(2)*Norm(2))); + Phase(3) = 0; + VoxDims = [geometry.size.VoI_PeFOV geometry.size.VoI_RoFOV geometry.size.VoIThickness]; + end + PhaseRot = rda.SlabOrientation(2,:)'; + % swapped signs because VoI_InPlaneRot rotates counterclockwise + % around Norm + if dot(cross(Phase,PhaseRot),Norm) <= 0 + geometry.rot.VoI_InPlaneRot = abs(acos(dot(Phase,PhaseRot))); + else + geometry.rot.VoI_InPlaneRot = -abs(acos(dot(Phase,PhaseRot))); + end + %End of calculating InPlaneRot + else + geometry.size.VoI_RoFOV = rda.VOIReadoutFOV; % Voxel size in readout direction [mm] + geometry.size.VoI_PeFOV = rda.VOIPhaseFOV; % Voxel size in phase encoding direction [mm] + geometry.size.VoIThickness = rda.VOIThickness; % Voxel size in slice selection direction [mm] + geometry.pos.PosCor = rda.VOIPositionCor; % Coronal coordinate of voxel [mm] + geometry.pos.PosSag = rda.VOIPositionSag; % Sagittal coordinate of voxel [mm] + geometry.pos.PosTra = rda.VOIPositionTra; % Transversal coordinate of voxel [mm] + geometry.rot.VoI_InPlaneRot = rda.VOIRotationInPlane; % Voxel rotation in plane + geometry.rot.NormCor = rda.VOINormalCor; % Coronal component of normal vector of voxel + geometry.rot.NormSag = rda.VOINormalSag; % Sagittal component of normal vector of voxel + geometry.rot.NormTra = rda.VOINormalTra; % Transversal component of normal vector of voxel + end + %mod. end% + out.geometry = geometry; + + + % + % So now we should have got to the point after the header text + % + % Siemens documentation suggests that the data should be in a double complex format (8bytes for real, and 8 for imaginary?) + % + fclose(fid); + bytes_per_point = 16; + + % Preallocate array in which the FIDs are to be extracted. + fids = zeros(length(filesInFolder),rda.VectorSize); + % Collect all FIDs and sort them into fids array + for kk = 1:length(filesInFolder) + fid = fopen(filesInFolder{kk}); + tline = fgets(fid); + %while (isempty(strfind(tline , head_end_text))) + while ~contains(tline , head_end_text) + tline = fgets(fid); + end + complex_data = fread(fid , rda.CSIMatrix_Size(1) * rda.CSIMatrix_Size(1) *rda.CSIMatrix_Size(1) *rda.VectorSize * 2 , 'double'); + %fread(fid , 1, 'double'); %This was a check to confirm that we had read all the data (it passed!) + fclose(fid); + + % Now convert this data into something meaningful + + %Reshape so that we can get the real and imaginary separated + hmm = reshape(complex_data, 2 , rda.VectorSize , rda.CSIMatrix_Size(1) , rda.CSIMatrix_Size(2) , rda.CSIMatrix_Size(3) ); + + %Combine the real and imaginary into the complex matrix + fids(kk,:) = complex(hmm(1,:,:,:,:),hmm(2,:,:,:,:)); + end + fids = fids'; + %Remove the redundant first element in the array + % Time_domain_data = reshape(hmm_complex, rda.VectorSize , rda.CSIMatrix_Size(1) , rda.CSIMatrix_Size(2) , rda.CSIMatrix_Size(3)); + % Time_domain_data = fft(fft(fftshift(fftshift(padarray(fftshift(fftshift(ifft(ifft(Time_domain_data,[],2),[],3),2),3),[0 16 16],'both'),2),3),[],2),[],3); + + % rewind first point phase + corrph = conj(fids(1,:))./abs(fids(1,:)); + corrph = repmat(corrph, [size(fids, 1), 1]); + fids = fids .* corrph; + + % get the spectrum from the fid + specs = fftshift(fft(fids,[],1),1); + + % make calculations for the output mrs structure + sz = size(fids); + dwelltime = rda.DwellTime/1000000; + spectralwidth=1/dwelltime; + txfrq = rda.MRFrequency*1000000; + dims.t = 1; + dims.subSpecs = 0; + dims.coils = 0; + dims.extras = 0; + dims.averages = 2; + % Bo = rda.MagneticFieldStrength; + Bo = rda.MRFrequency/42.577; + rawAverages = rda.NumberOfAverages; + if length(filesInFolder) >= rawAverages + if contains(seqtype, 'edit') + rawAverages = length(filesInFolder); + averages = rda.NumberOfAverages; + subspecs = fix(rawAverages/averages); + rawSubspecs = rda.NumberOfAverages; + else + rawAverages = length(filesInFolder); + averages = rawAverages; + subspecs = 1; + rawSubspecs = 1; + end + elseif length(filesInFolder) == 1 + averages = 1; + subspecs =1; + rawSubspecs = 'na'; + else + if contains(seqtype,'edit') || contains(seqtype,'mpress') + subspecs = 2; + rawSubspecs = 2; + % Distinguish whether MEGA data have been averaged on-scanner or + % not + if length(filesInFolder) == 2 + dims.subSpecs = 3; + dims.averages = 2; + averages = 1; + rawAverages = 1; + fids = reshape(fids,[sz(1),averages,2 ]); + specs = reshape(specs,[sz(1),averages,2 ]); + else + dims.subSpecs = 3; + dims.averages = 2; + averages = rawAverages/2; + fids = reshape(fids,[sz(1),averages,2 ]); + specs = reshape(specs,[sz(1),averages,2 ]); + end + sz = size(fids); + end + end + + date = rda.StudyDate; + seq = rda.SequenceDescription; + TE = rda.TE; + TR = rda.TR; + pointsToLeftShift = 0; + + %Calculate t and ppm arrays using the calculated parameters: + f=[(-spectralwidth/2)+(spectralwidth/(2*sz(1))):spectralwidth/(sz(1)):(spectralwidth/2)-(spectralwidth/(2*sz(1)))]; + ppm=f/(Bo*42.577); + % Siemens data assumes the center frequency to be 4.65 ppm: + centerFreq = 4.65; + ppm=ppm + centerFreq; + + t=[0:dwelltime:(sz(1)-1)*dwelltime]; + + %FILLING IN DATA STRUCTURE + out.fids=fids; + out.specs=specs; + out.sz=sz; + out.ppm=ppm; + out.t=t; + out.spectralwidth=spectralwidth; + out.dwelltime=dwelltime; + out.txfrq=txfrq; + out.date=date; + out.dims=dims; + out.Bo=Bo; + out.averages=averages; + out.rawAverages=rawAverages; + out.subspecs=subspecs; + out.rawSubspecs=rawSubspecs; + out.seq=seq; + out.te=TE; + out.tr=TR; + out.pointsToLeftshift=pointsToLeftShift; + out.centerFreq = centerFreq; + out.geometry = geometry; + out.software = rda.software_version; + + %FILLING IN THE FLAGS + out.flags.writtentostruct=1; + out.flags.gotparams=1; + out.flags.leftshifted=0; + out.flags.filtered=0; + out.flags.zeropadded=0; + out.flags.freqcorrected=0; + out.flags.phasecorrected=0; + if out.averages == 1 + out.flags.averaged=1; + else + out.flags.averaged=0; + end + out.flags.addedrcvrs=1; + out.flags.subtracted=0; + out.flags.writtentotext=0; + out.flags.downsampled=0; + if out.dims.subSpecs==0 + out.flags.isFourSteps=0; + else + out.flags.isFourSteps=(out.sz(out.dims.subSpecs)==4); + end + + % Add info for niiwrite + out.PatientPosition = rda.PatientPosition; + out.Manufacturer = 'Siemens'; + [~,filename,ext] = fileparts(filesInFolder{1}); + out.OriginalFile = [filename ext]; + + % Sequence flags + out.flags.isUnEdited = 0; + out.flags.isMEGA = 0; + out.flags.isHERMES = 0; + out.flags.isHERCULES = 0; + out.flags.isPRIAM = 0; + out.flags.isMRSI = 0; +end