diff --git a/AlignSubSpectra.m b/AlignSubSpectra.m index 5d30361..4600d07 100644 --- a/AlignSubSpectra.m +++ b/AlignSubSpectra.m @@ -70,7 +70,7 @@ else subSpecInd = [3 2 1 4]; end - case 'nifti' + case 'NIfTI' subSpecInd = [3 2 1 4]; case {'Philips','Philips_data','Philips_raw'} subSpecInd = [1 2 3 4]; @@ -89,7 +89,7 @@ end else switch MRS_struct.p.vendor - case {'GE','nifti'} + case {'GE','NIfTI'} subSpecInd = [3 2 1 4]; case {'Philips','Philips_data','Philips_raw'} subSpecInd = [1 4 3 2]; diff --git a/CoRegStandAlone/CoReg.m b/CoRegStandAlone/CoReg.m index 439a640..c136a3e 100755 --- a/CoRegStandAlone/CoReg.m +++ b/CoRegStandAlone/CoReg.m @@ -40,6 +40,17 @@ %Ultimately this switch will not be necessary... switch MRS_struct.p.vendor + + case 'GE' + [~,~,ext] = fileparts(struc{ii}); + if strcmp(ext,'.nii') + MRS_struct = GannetMask_GE_nii(MRS_struct.metabfile{ii}, struc{ii}, MRS_struct, ii, vox, kk); + else + MRS_struct = GannetMask_GE(MRS_struct.metabfile{ii}, struc{ii}, MRS_struct, ii, vox, kk); + end + + case 'NIfTI' + MRS_struct = GannetMask_NIfTI(MRS_struct.metabfile{ii}, struc{ii}, MRS_struct, ii, vox, kk); case 'Philips' sparname = [MRS_struct.metabfile{ii}(1:(end-4)) MRS_struct.p.spar_string]; @@ -73,20 +84,9 @@ case 'Siemens_rda' MRS_struct = GannetMask_SiemensRDA(MRS_struct.metabfile{ii}, struc{ii}, MRS_struct, ii, vox, kk); - case {'Siemens_twix', 'Siemens_dicom', 'dicom'} + case {'Siemens_twix', 'Siemens_dicom', 'DICOM'} MRS_struct = GannetMask_SiemensTWIX(MRS_struct.metabfile{ii}, struc{ii}, MRS_struct, ii, vox, kk); - case 'GE' - [~,~,ext] = fileparts(struc{ii}); - if strcmp(ext,'.nii') - MRS_struct = GannetMask_GE_nii(MRS_struct.metabfile{ii}, struc{ii}, MRS_struct, ii, vox, kk); - else - MRS_struct = GannetMask_GE(MRS_struct.metabfile{ii}, struc{ii}, MRS_struct, ii, vox, kk); - end - - case 'nifti' - error('NIfTI not yet supported.'); - end % Build output figure @@ -164,7 +164,7 @@ imagesc(squeeze(MRS_struct.mask.(vox{kk}).img{ii})); colormap('gray'); img = MRS_struct.mask.(vox{kk}).img{ii}(:); - caxis([0 mean(img(img > 0.01)) + 3*std(img(img > 0.01))]); + caxis([0 mean(img(img > 0.01)) + 3*std(img(img > 0.01))]); %#ok<*CAXIS> axis equal tight 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); @@ -172,8 +172,8 @@ title(t, 'FontName', 'Arial', 'FontSize', 15, 'Interpreter', 'none'); % Gannet logo - Gannet_logo = fullfile(fileparts(which('GannetLoad')), 'Gannet3_logo.png'); - I = imread(Gannet_logo,'png','BackgroundColor',[1 1 1]); + Gannet_logo = fullfile(fileparts(which('GannetLoad')), 'Gannet3_logo.jpg'); + I = imread(Gannet_logo); axes('Position', [0.825, 0.05, 0.125, 0.125]); imshow(I); text(0.9, 0, MRS_struct.version.Gannet, 'Units', 'normalized', 'FontName', 'Arial', 'FontSize', 14, 'FontWeight', 'bold', 'HorizontalAlignment', 'left'); diff --git a/CoRegStandAlone/CoRegStandAlone.m b/CoRegStandAlone/CoRegStandAlone.m index ee6ed4e..a45899b 100755 --- a/CoRegStandAlone/CoRegStandAlone.m +++ b/CoRegStandAlone/CoRegStandAlone.m @@ -91,13 +91,13 @@ switch MRS_struct.p.vendor - case 'dicom' + case 'DICOM' MRS_struct = DICOMRead(MRS_struct, metabfile{ii}); case 'GE' MRS_struct = GERead(MRS_struct, metabfile{ii}); - case 'nifti' + case 'NIfTI' MRS_struct = NIfTIMRSRead(MRS_struct, metabfile{ii}); case 'Philips' diff --git a/CoRegStandAlone/Seg.m b/CoRegStandAlone/Seg.m index dd945c0..3f3405e 100755 --- a/CoRegStandAlone/Seg.m +++ b/CoRegStandAlone/Seg.m @@ -248,8 +248,8 @@ title(t, 'FontName', 'Arial', 'FontSize', 15, 'Interpreter', 'none'); % Gannet logo - Gannet_logo = fullfile(fileparts(which('GannetLoad')), 'Gannet3_logo.png'); - I = imread(Gannet_logo,'png','BackgroundColor',[1 1 1]); + Gannet_logo = fullfile(fileparts(which('GannetLoad')), 'Gannet3_logo.jpg'); + I = imread(Gannet_logo); axes('Position', [0.825, 0.05, 0.125, 0.125]); imshow(I); text(0.9, 0, MRS_struct.version.Gannet, 'Units', 'normalized', 'FontName', 'Arial', 'FontSize', 14, 'FontWeight', 'bold', 'HorizontalAlignment', 'left'); diff --git a/DiscernDataType.m b/DiscernDataType.m index ee1e64d..9473956 100644 --- a/DiscernDataType.m +++ b/DiscernDataType.m @@ -5,6 +5,20 @@ switch lower(ext) case '.7' MRS_struct.p.vendor = 'GE'; + case '.dat' + MRS_struct.p.vendor = 'Siemens_twix'; + case '.data' + MRS_struct.p.vendor = 'Philips_data'; + case '.dcm' + MRS_struct.p.vendor = 'DICOM'; + case {'.gz','.nii'} + MRS_struct.p.vendor = 'NIfTI'; + case '.ima' + MRS_struct.p.vendor = 'Siemens_dicom'; + case '.raw' + MRS_struct.p.vendor = 'Philips_raw'; + case '.rda' + MRS_struct.p.vendor = 'Siemens_rda'; case '.sdat' MRS_struct.p.vendor = 'Philips'; if all(isstrprop(ext(end-3:end), 'upper')) @@ -12,22 +26,8 @@ else MRS_struct.p.spar_string = 'spar'; end - case '.data' - MRS_struct.p.vendor = 'Philips_data'; - case '.raw' - MRS_struct.p.vendor = 'Philips_raw'; - case '.rda' - MRS_struct.p.vendor = 'Siemens_rda'; - case '.dat' - MRS_struct.p.vendor = 'Siemens_twix'; - case '.ima' - MRS_struct.p.vendor = 'Siemens_dicom'; - case '.dcm' - MRS_struct.p.vendor = 'dicom'; - case {'.gz','.nii'} - MRS_struct.p.vendor = 'nifti'; otherwise - error('Unrecognized file type! Extension should be .7, .sdat, .data, .raw, .rda, .dat, .ima, .dcm, .gz, or .nii.'); + error('Unrecognized file type! Extension should be .7, .dat, .data, .dcm, .gz, .ima, .nii, .raw, .rda, or, .sdat.'); end end diff --git a/ExportToCSV.m b/ExportToCSV.m index 790a062..15ce617 100644 --- a/ExportToCSV.m +++ b/ExportToCSV.m @@ -9,7 +9,7 @@ function ExportToCSV(MRS_struct, vox, module) end out.MATLAB_ver = cellstr(repmat(version('-release'), n_rep)); out.Gannet_ver = cellstr(repmat(MRS_struct.version.Gannet, n_rep)); -out.date_of_analysis = cellstr(repmat(datestr(date, 'yyyy-mm-dd'), n_rep)); +out.date_of_analysis = cellstr(repmat(datestr(date, 'yyyy-mm-dd'), n_rep)); %#ok<*DATE,*DATST> %%% 1. Extract data from GannetFit %%% diff --git a/GERead.m b/GERead.m index c8b6404..9cac81c 100644 --- a/GERead.m +++ b/GERead.m @@ -215,7 +215,6 @@ MRS_struct.p.LarmorFreq(ii) = i_hdr_value(rdb_hdr_ps_mps_freq)/1e7; MRS_struct.p.sw(ii) = f_hdr_value(rdb_hdr_user0); - nechoes = hdr_value(rdb_hdr_nechoes); MRS_struct.p.GE.nechoes(ii) = nechoes; nex = hdr_value(rdb_hdr_navs); @@ -286,12 +285,12 @@ if nechoes == 1 if (dataframes + refframes) ~= nframes - mult = 1; + mult = 1; MRS_struct.p.GE.noadd(ii) = 1; - dataframes = dataframes * nex; - refframes = nframes - dataframes; + dataframes = dataframes * nex; + refframes = nframes - dataframes; else - mult = 1/nex; + mult = 1/nex; MRS_struct.p.GE.noadd(ii) = 0; end @@ -301,8 +300,8 @@ totalframes = totalframes - (refframes + 1); - MRS_struct.p.nrows(ii) = totalframes; - MRS_struct.p.Navg(ii) = dataframes * nex; + MRS_struct.p.nrows(ii) = totalframes; + MRS_struct.p.Navg(ii) = dataframes * nex; MRS_struct.p.Nwateravg(ii) = refframes * nex; else @@ -310,20 +309,14 @@ MRS_struct.p.Navg(ii) = dataframes * nex * nechoes; % RTN 2017 if (dataframes + refframes) ~= nframes - mult = nex/2; % RTN 2016 - multw = nex; % RTN 2016 - %mult = 1; % RTN 2017 - %multw = 1; % RTN 2017 + mult = nex/2; % RTN 2016 1; % RTN 2017 + multw = nex; % RTN 2016 1; % RTN 2017 MRS_struct.p.GE.noadd(ii) = 1; - dataframes = dataframes * nex; - refframes = nframes - dataframes; + dataframes = dataframes * nex; + refframes = nframes - dataframes; else - %mult = nex/2; % RTN 2016 - %multw = 1; % RTN 2016 - %mult = 1; % RTN 2017 - %multw = 1/nex; % RTN 2017 - mult = 1/nex; % MM 2020 - multw = 1; % MM 2020 + mult = 1/nex; % MM 2020 1; % RTN 2017 nex/2; % RTN 2016 + multw = 1; % MM 2020 1/nex; % RTN 2017 1; % RTN 2016 MRS_struct.p.GE.noadd(ii) = 0; end @@ -351,7 +344,7 @@ Y2 = 1 + refframes + (totalframes/nechoes) * (X2-1) + X1; MetabData = Y1 .* ShapeData(:,:,Y2,:) * mult; - totalframes = totalframes - (refframes + 1) * nechoes; % RTN 2017 + totalframes = totalframes - (refframes + 1) * nechoes; % RTN 2017 MRS_struct.p.nrows(ii) = totalframes; end @@ -361,7 +354,11 @@ WaterData = squeeze(complex(WaterData(1,:,:,:), WaterData(2,:,:,:))); WaterData = permute(WaterData, [3 1 2]); -% Generalized least squares method (An et al., JMRI, 2013, doi:10.1002/jmri.23941) +% Combine coils using generalized least squares method (An et al., JMRI, +% 2013, doi:10.1002/jmri.23941); the noise covariance matrix is more +% optionally estimated by using all averages as suggested by Rodgers & +% Robson (MRM, 2010, doi:10.1002/mrm.22230) + [nCh, nPts, nReps] = size(WaterData); noise_pts = false(1,nPts); noise_pts(ceil(0.75*nPts):end) = true; diff --git a/Gannet3_logo.jpg b/Gannet3_logo.jpg new file mode 100644 index 0000000..c6ba6ab Binary files /dev/null and b/Gannet3_logo.jpg differ diff --git a/Gannet3_logo.png b/Gannet3_logo.png deleted file mode 100644 index e660913..0000000 Binary files a/Gannet3_logo.png and /dev/null differ diff --git a/GannetCoRegister.m b/GannetCoRegister.m index edaf7dd..2b47a58 100644 --- a/GannetCoRegister.m +++ b/GannetCoRegister.m @@ -6,7 +6,7 @@ error('MATLAB:minrhs','Not enough input arguments.'); end -MRS_struct.version.coreg = '220923'; +MRS_struct.version.coreg = '221021'; warning('off'); % temporarily suppress warning messages @@ -36,7 +36,7 @@ run_count = 0; for ii = 1:MRS_struct.p.numScans - + [~,b,c] = fileparts(MRS_struct.metabfile{1,ii}); [~,e,f] = fileparts(struc{ii}); if strcmpi(f, '.gz') @@ -52,10 +52,10 @@ % Loop over voxels if PRIAM for kk = 1:length(vox) - + switch MRS_struct.p.vendor - case 'GE' + case 'GE' [~,~,ext] = fileparts(struc{ii}); if strcmp(ext,'.nii') MRS_struct = GannetMask_GE_nii(fname, struc{ii}, MRS_struct, ii, vox, kk); @@ -63,14 +63,13 @@ MRS_struct = GannetMask_GE(fname, struc{ii}, MRS_struct, ii, vox, kk); end - case 'nifti' - error('NIfTI not yet supported.'); -% MRS_struct = GannetMask_NIfTI(fname, struc{ii}, MRS_struct, ii, vox, kk); + case 'NIfTI' + MRS_struct = GannetMask_NIfTI(fname, struc{ii}, MRS_struct, ii, vox, kk); case 'Philips' sparname = [MRS_struct.metabfile{1,ii}(1:(end-4)) MRS_struct.p.spar_string]; MRS_struct = GannetMask_Philips(sparname, struc{ii}, MRS_struct, ii, vox, kk); - + case 'Philips_data' if exist(MRS_struct.metabfile_sdat,'file') MRS_struct.p.vendor = 'Philips'; @@ -100,7 +99,7 @@ fname = MRS_struct.metabfile{1,ii*2-1}; MRS_struct = GannetMask_SiemensRDA(fname, struc{ii}, MRS_struct, ii, vox, kk); - case {'Siemens_dicom', 'Siemens_twix', 'dicom'} + case {'Siemens_dicom', 'Siemens_twix', 'DICOM'} MRS_struct = GannetMask_SiemensTWIX(fname, struc{ii}, MRS_struct, ii, vox, kk); end @@ -122,10 +121,10 @@ set(h,'Color',[1 1 1]); figTitle = 'GannetCoRegister Output'; set(h,'Name',figTitle,'Tag',figTitle,'NumberTitle','off'); - + subplot(2,3,4:6); axis off; - + [~,tmp,tmp2] = fileparts(MRS_struct.mask.(vox{kk}).outfile{ii}); fname = [tmp tmp2]; if length(fname) > 30 @@ -133,22 +132,22 @@ end text(0.5, 0.75, 'Mask output: ', 'HorizontalAlignment' , 'right', 'FontName', 'Arial', 'FontSize', 13); text(0.5, 0.75, [' ' fname], 'FontName', 'Arial', 'FontSize', 13, 'Interpreter', 'none'); - + text(0.5, 0.63, 'Spatial parameters: ', 'HorizontalAlignment', 'right', 'FontName', 'Arial', 'FontSize', 13); text(0.5, 0.63, ' [LR, AP, FH]', 'FontName', 'Arial', 'FontSize', 13); - + tmp = [' ' num2str(MRS_struct.p.voxdim(ii,1)) ' \times ' num2str(MRS_struct.p.voxdim(ii,2)) ' \times ' num2str(MRS_struct.p.voxdim(ii,3)) ' mm^{3}']; text(0.5, 0.51, 'Dimensions: ', 'HorizontalAlignment', 'right', 'FontName', 'Arial', 'FontSize', 13); text(0.5, 0.51, tmp, 'FontName', 'Arial', 'FontSize', 13, 'Interpreter', 'tex'); - + tmp = [' ' num2str(prod(MRS_struct.p.voxdim(ii,:))/1e3) ' mL']; text(0.5, 0.39, 'Volume: ', 'HorizontalAlignment', 'right', 'FontName', 'Arial', 'FontSize', 13); text(0.5, 0.39, tmp, 'FontName', 'Arial', '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, 'Position: ', 'HorizontalAlignment', 'right', 'FontName', 'Arial', 'FontSize', 13); text(0.5, 0.27, tmp, 'FontName', 'Arial', 'FontSize', 13); - + if any(strcmp(MRS_struct.p.vendor,{'Philips','Philips_data'})) tmp = [' [' num2str(MRS_struct.p.voxang(ii,2), '%3.1f') ', ' num2str(MRS_struct.p.voxang(ii,1), '%3.1f') ', ' num2str(MRS_struct.p.voxang(ii,3), '%3.1f') '] deg']; else @@ -156,12 +155,12 @@ end text(0.5, 0.15, 'Angulation: ', 'HorizontalAlignment', 'right', 'FontName', 'Arial', 'FontSize', 13); text(0.5, 0.15, tmp, 'FontName', 'Arial', 'FontSize', 13); - + text(0.5, 0.03, 'CoRegVer: ', 'HorizontalAlignment', 'right', 'FontName', 'Arial', 'FontSize', 13); text(0.5, 0.03, [' ' MRS_struct.version.coreg], 'FontName', 'Arial', 'FontSize', 13); - + ha = subplot(2,3,1:3); - + if strcmp(MRS_struct.p.vendor,'Siemens_rda') [~,tmp,tmp2] = fileparts(MRS_struct.metabfile{1,ii*2-1}); else @@ -177,19 +176,19 @@ T1image = [T1image(1:12) '...' T1image(end-11:end)]; end t = ['Voxel from ' fname ' on ' T1image]; - + imagesc(MRS_struct.mask.(vox{kk}).img{ii}); axis equal tight 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); set(ha,'pos',[0 0.15 1 1]); title(t, 'FontName', 'Arial', 'FontSize', 15, 'Interpreter', 'none'); - + % Save output as PDF run_count = SavePDF(h, MRS_struct, ii, 1, kk, vox, mfilename, run_count); end - + end % Save MRS_struct as mat file diff --git a/GannetFit.m b/GannetFit.m index 4bf4010..fec1699 100644 --- a/GannetFit.m +++ b/GannetFit.m @@ -45,6 +45,7 @@ warning('off','stats:nlinfit:ModelConstantWRTParam'); warning('off','stats:nlinfit:IllConditionedJacobian'); +warning('off','stats:nlinfit:IterationLimitExceeded'); warning('off','MATLAB:rankDeficientMatrix'); % Loop over voxels if PRIAM @@ -1069,7 +1070,7 @@ size_max = size(MRS_struct.mask.img{ii},1); imagesc(MRS_struct.mask.img{ii}(:,size_max+(1:size_max))); colormap('gray'); - caxis([0 1]) + caxis([0 1]); %#ok<*CAXIS> axis equal tight off; subplot(2,2,4,'replace'); else @@ -1427,6 +1428,7 @@ warning('on','stats:nlinfit:ModelConstantWRTParam'); warning('on','stats:nlinfit:IllConditionedJacobian'); +warning('on','stats:nlinfit:IterationLimitExceeded'); warning('on','MATLAB:rankDeficientMatrix'); % Need to close hidden figures to show figures after Gannet is done running @@ -1435,6 +1437,8 @@ end + + %%%%%%%%%%%%%%%% GAUSS MODEL %%%%%%%%%%%%%%%% function F = GaussModel(x,freq) % Function for Gauss Model diff --git a/GannetLoad.m b/GannetLoad.m index 69b58f0..af86ada 100644 --- a/GannetLoad.m +++ b/GannetLoad.m @@ -17,8 +17,8 @@ error('MATLAB:minrhs','Not enough input arguments.'); end -MRS_struct.version.Gannet = '3.3.0-rc'; -MRS_struct.version.load = '221013'; +MRS_struct.version.Gannet = '3.3.0'; +MRS_struct.version.load = '221022'; VersionCheck(0, MRS_struct.version.Gannet); ToolboxCheck; @@ -212,11 +212,11 @@ else switch MRS_struct.p.vendor - case 'dicom' + case 'DICOM' loadFun = @DICOMRead; case 'GE' loadFun = @GERead; - case 'nifti' + case 'NIfTI' loadFun = @NIfTIMRSRead; case 'Philips' loadFun = @PhilipsRead; diff --git a/GannetMask_NIfTI.m b/GannetMask_NIfTI.m index ebdc96e..9d9512f 100644 --- a/GannetMask_NIfTI.m +++ b/GannetMask_NIfTI.m @@ -1,106 +1,55 @@ -function MRS_struct = GannetMask_NIfTI(fname, struc, MRS_struct, ii, vox, kk) +function MRS_struct = GannetMask_NIfTI(fname, nii_file, MRS_struct, ii, vox, kk) -% Co-register NIfTI-MRS files to structural images in NIfTI format. Code heavily -% based on Dr. Peter Van Schuerbbeek's (UZ Brussel) coreg_p code and Ralph -% Noeske's (GE Berlin) SV_MRI voxel co-registration code. +% Co-register NIfTI-MRS files to structural images in NIfTI format. Code +% heavily based on coreg_nifti.m from Osprey. -voxoff = MRS_struct.p.voxoff(ii,:); +% CREDITS: +% Chris Davies-Jenkins, Johns Hopkins University 2022. +% Xiangrui Li, Ph.D. for his helpful suggestions using nii_tool. -% Load in NIfTI file -V = spm_vol(struc); -[T1, XYZ] = spm_read_vols(V); -MRS_struct.mask.(vox{kk}).T1max(ii) = max(T1(:)); +nii_struc = nii_tool('load', nii_file); % load structural nifti +nii_mrsvox = nii_tool('load', fname); % load voxel nifti + +% nii_viewer(NiiStruct, NiiVox); % overlay voxel on structural + +% Assume MRS voxel and structural are in same space +nii_mrsvox.hdr.sform_code = nii_struc.hdr.sform_code; +nii_mrsvox.hdr.qform_code = nii_struc.hdr.qform_code; + +nii_mrsvox.img = 1; % overwrites image, so mask +nii_mrsvox.hdr.dim(4:end) = 1; % remove additional MRS dimensions from header -% Shift imaging voxel coordinates by half an imaging voxel so that the XYZ matrix -% tells us the x,y,z coordinates of the MIDDLE of that imaging voxel. -[~,voxdim] = spm_get_bbox(V,'fv'); -voxdim = abs(voxdim)'; -halfpixshift = -voxdim(1:3)/2; -halfpixshift(3) = -halfpixshift(3); -XYZ = XYZ + repmat(halfpixshift, [1 size(XYZ,2)]); - -dXYZ = sqrt((XYZ(1,:) - voxoff(1)).^2 + ... - (XYZ(2,:) - voxoff(2)).^2 + ... - (XYZ(3,:) - voxoff(3)).^2); -[~,refvox] = min(dXYZ); -[refvox_x, refvox_y, refvox_z] = ind2sub(V.dim, refvox(1)); - -e1_edge = voxoff + e1_SVS; -e2_edge = voxoff + e2_SVS; -e3_edge = voxoff + e3_SVS; - -de1_XYZ = sqrt((XYZ(1,:) - e1_edge(1)).^2 + ... - (XYZ(2,:) - e1_edge(2)).^2 + ... - (XYZ(3,:) - e1_edge(3)).^2); -[~,e1vox] = min(de1_XYZ); -[e1vox_x, e1vox_y, e1vox_z] = ind2sub(V.dim, e1vox(1)); - -de2_XYZ = sqrt((XYZ(1,:) - e2_edge(1)).^2 + ... - (XYZ(2,:) - e2_edge(2)).^2 + ... - (XYZ(3,:) - e2_edge(3)).^2); -[~,e2vox] = min(de2_XYZ); -[e2vox_x, e2vox_y, e2vox_z] = ind2sub(V.dim, e2vox(1)); - -de3_XYZ = sqrt((XYZ(1,:) - e3_edge(1)).^2 + ... - (XYZ(2,:) - e3_edge(2)).^2 + ... - (XYZ(3,:) - e3_edge(3)).^2); -[~,e3vox] = min(de3_XYZ); -[e3vox_x, e3vox_y, e3vox_z] = ind2sub(V.dim, e3vox(1)); - -% Create a mask with all voxels that are inside the voxel -mask = zeros(V.dim); - -nx = floor(sqrt((e1vox_x - refvox_x)^2 + (e1vox_y - refvox_y)^2 + (e1vox_z - refvox_z)^2)) * 2; -ny = floor(sqrt((e2vox_x - refvox_x)^2 + (e2vox_y - refvox_y)^2 + (e2vox_z - refvox_z)^2)) * 2; -nz = floor(sqrt((e3vox_x - refvox_x)^2 + (e3vox_y - refvox_y)^2 + (e3vox_z - refvox_z)^2)) * 2; - -stepx = ([e1vox_x, e1vox_y, e1vox_z] - [refvox_x, refvox_y, refvox_z]) / nx; -stepy = ([e2vox_x, e2vox_y, e2vox_z] - [refvox_x, refvox_y, refvox_z]) / ny; -stepz = ([e3vox_x, e3vox_y, e3vox_z] - [refvox_x, refvox_y, refvox_z]) / nz; - -mrs_box_ind = 1:(nx * ny * nz); -mrs_box_sub = zeros(3, nx * ny * nz); - -[mrs_box_sub_x, mrs_box_sub_y, mrs_box_sub_z] = ind2sub([nx, ny, nz], mrs_box_ind); - -mrs_box_sub(1,:) = mrs_box_sub_x; -mrs_box_sub(2,:) = mrs_box_sub_y; -mrs_box_sub(3,:) = mrs_box_sub_z; - -e1_stepx = repmat(stepx, [numel(mrs_box_sub(1,:)), 1])'; -e2_stepy = repmat(stepy, [numel(mrs_box_sub(1,:)), 1])'; -e3_stepz = repmat(stepz, [numel(mrs_box_sub(1,:)), 1])'; - -mrs_box_sub = repmat((mrs_box_sub(1,:) - 1), [3 1]) .* e1_stepx + ... - repmat((mrs_box_sub(2,:) - 1), [3 1]) .* e2_stepy + ... - repmat((mrs_box_sub(3,:) - 1), [3 1]) .* e3_stepz; -refvox_rep = repmat([refvox_x, refvox_y, refvox_z], [numel(mrs_box_sub(1,:)), 1])'; -mrs_box_sub = round(mrs_box_sub + refvox_rep); -mrs_box_ind = sub2ind(V.dim, mrs_box_sub(1,:), ... - mrs_box_sub(2,:), ... - mrs_box_sub(3,:)); -mask(mrs_box_ind) = 1; - -% Build output (code to make voxel mask yellow borrowed from SPM12) - -[a,b] = fileparts(fname); +[a,b,c] = fileparts(fname); if isempty(a) a = '.'; end -V_mask.fname = fullfile([a filesep b '_mask.nii']); -V_mask.descrip = 'MRS_voxel_mask'; -V_mask.dim = V.dim; +if strcmpi(c,'.gz') + b(end-3:end) = []; +end +mask_filename = fullfile([a filesep b '_mask.nii']); +% Transform voxel to image resolution and save under mask_filename +nii_xform(nii_mrsvox, nii_struc.hdr, mask_filename, 'linear', 0); + +% Load structural using SPM +V = spm_vol(nii_file); +T1 = spm_read_vols(V); +MRS_struct.mask.(vox{kk}).T1max(ii) = max(T1(:)); + +% Load mask using SPM to adapt some fields +V_mask = spm_vol(mask_filename); V_mask.dt = V.dt; -V_mask.mat = V.mat; -V_mask = spm_write_vol(V_mask, mask); +V_mask.descrip = 'MRS_voxel_mask'; +V_mask = spm_write_vol(V_mask, V_mask.private.dat(:,:,:)); % write mask MRS_struct.mask.(vox{kk}).outfile(ii,:) = cellstr(V_mask.fname); +% Not clear how to formulate the rotations for triple rotations (revisit) +MRS_struct.p.voxang(ii,:) = [NaN NaN NaN]; -% Transform structural image and co-registered voxel mask from voxel to -% world space for output -voxel_ctr = MRS_struct.p.voxoff(ii,:); -[img_t, img_c, img_s] = voxel2world_space(V, voxel_ctr); -[mask_t, mask_c, mask_s] = voxel2world_space(V_mask, voxel_ctr); +voxel_ctr = [nii_mrsvox.hdr.qoffset_x, nii_mrsvox.hdr.qoffset_y, nii_mrsvox.hdr.qoffset_z]; % CWDJ need to verify this + +MRS_struct.p.voxoff(ii,:) = voxel_ctr; +[img_t, img_c, img_s] = voxel2world_space(V, voxel_ctr); +[mask_t, mask_c, mask_s] = voxel2world_space(V_mask, voxel_ctr); w_t = zeros(size(img_t)); w_c = zeros(size(img_c)); @@ -163,6 +112,9 @@ three_plane_img(:,size_max*2+(1:size_max),:) = image_center(img_c, size_max); MRS_struct.mask.(vox{kk}).img{ii} = three_plane_img; -MRS_struct.mask.(vox{kk}).T1image(ii,:) = {struc}; +MRS_struct.mask.(vox{kk}).T1image(ii,:) = {nii_file}; + +end + + -end \ No newline at end of file diff --git a/GannetMask_SiemensRDA.m b/GannetMask_SiemensRDA.m index 07e2fc1..ddf39df 100644 --- a/GannetMask_SiemensRDA.m +++ b/GannetMask_SiemensRDA.m @@ -43,28 +43,27 @@ value = tline(occurence_of_colon+1:length(tline)); switch variable - case {'VOINormalSag' , 'VOINormalCor' , 'VOINormalTra' , 'VOIPositionSag', ... - 'VOIPositionCor', 'VOIPositionTra', 'VOIThickness','VOIReadoutFOV', ... - 'VOIPhaseFOV','VOIRotationInPlane'} + case {'VOINormalSag', 'VOINormalCor', 'VOINormalTra', 'VOIPositionSag', ... + 'VOIPositionCor', 'VOIPositionTra', 'VOIThickness', 'VOIReadoutFOV', ... + 'VOIPhaseFOV', 'VOIRotationInPlane'} rda.(variable) = str2double(value); end end end -MRS_struct.p.VoI_InPlaneRot(ii) = rda.VOIRotationInPlane; -MRS_struct.p.NormCor(ii) = rda.VOINormalCor; -MRS_struct.p.NormSag(ii) = rda.VOINormalSag; -MRS_struct.p.NormTra(ii) = rda.VOINormalTra; -MRS_struct.p.voxdim(ii,1) = rda.VOIPhaseFOV; -MRS_struct.p.voxdim(ii,2) = rda.VOIReadoutFOV; -MRS_struct.p.voxdim(ii,3) = rda.VOIThickness; -MRS_struct.p.voxoff(ii,1) = rda.VOIPositionSag; -MRS_struct.p.voxoff(ii,2) = rda.VOIPositionCor; -MRS_struct.p.voxoff(ii,3) = rda.VOIPositionTra; +MRS_struct.p.VoI_InPlaneRot(ii) = rda.VOIRotationInPlane; +MRS_struct.p.NormCor(ii) = rda.VOINormalCor; +MRS_struct.p.NormSag(ii) = rda.VOINormalSag; +MRS_struct.p.NormTra(ii) = rda.VOINormalTra; +MRS_struct.p.voxdim(ii,1) = rda.VOIPhaseFOV; +MRS_struct.p.voxdim(ii,2) = rda.VOIReadoutFOV; +MRS_struct.p.voxdim(ii,3) = rda.VOIThickness; +MRS_struct.p.voxoff(ii,1) = rda.VOIPositionSag; +MRS_struct.p.voxoff(ii,2) = rda.VOIPositionCor; +MRS_struct.p.voxoff(ii,3) = rda.VOIPositionTra; fclose(fid); - % Extract voxel position and rotation parameters from MRS_struct NormSag = MRS_struct.p.NormSag(ii); NormCor = MRS_struct.p.NormCor(ii); @@ -109,24 +108,24 @@ 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 = [MRS_struct.p.voxdim(ii,1) MRS_struct.p.voxdim(ii,2) MRS_struct.p.voxdim(ii,3)]; + 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 = [MRS_struct.p.voxdim(ii,1) MRS_struct.p.voxdim(ii,2) MRS_struct.p.voxdim(ii,3)]; 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 = [MRS_struct.p.voxdim(ii,1) MRS_struct.p.voxdim(ii,2) MRS_struct.p.voxdim(ii,3)]; + 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 = [MRS_struct.p.voxdim(ii,1) MRS_struct.p.voxdim(ii,2) MRS_struct.p.voxdim(ii,3)]; 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 = [MRS_struct.p.voxdim(ii,1) MRS_struct.p.voxdim(ii,2) MRS_struct.p.voxdim(ii,3)]; + 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 = [MRS_struct.p.voxdim(ii,1) MRS_struct.p.voxdim(ii,2) MRS_struct.p.voxdim(ii,3)]; end % The readout reference vector is the cross product of Norm and Phase @@ -137,18 +136,18 @@ M_R(1:3, 3) = Norm; % Define matrix for rotation around in-plane rotation angle -M3_Mu = [ cos(ROT) sin(ROT) 0 - -sin(ROT) cos(ROT) 0 - 0 0 1]; +M3_Mu = [ cos(ROT) sin(ROT) 0 + -sin(ROT) cos(ROT) 0 + 0 0 1]; -M3_R = M_R(1:3,1:3) * M3_Mu; -M_R(1:3,1:3) = M3_R; +M3_R = M_R(1:3,1:3) * M3_Mu; +M_R(1:3,1:3) = M3_R; % The MGH vox2ras matrix inverts the Readout column -M_R = M_R * [ 1 0 0 0 - 0 -1 0 0 - 0 0 1 0 - 0 0 0 1]; +M_R = M_R * [1 0 0 0 + 0 -1 0 0 + 0 0 1 0 + 0 0 0 1]; % Final rotation matrix rotmat = M_R(1:3,1:3); diff --git a/GannetPreInitialise.m b/GannetPreInitialise.m index 8a38910..a025987 100644 --- a/GannetPreInitialise.m +++ b/GannetPreInitialise.m @@ -3,7 +3,7 @@ % Acquisition parameters MRS_struct.p.target = {'GABAGlx'}; % Edited metabolite(s) of interest; permitted options are: % If MEGA-PRESS: - % {'GABAGlx'}, {'GSH'}, {'Lac'}, or {'EtOH'} + % {'GABA'}, {'GABAGlx'}, {'GSH'}, {'Lac'}, or {'EtOH'} % If HERMES: % {'GABAGlx','GSH'}, {'Lac','GSH'}, or {'EtOH','GABA','GSH'} % If HERCULES: @@ -38,8 +38,8 @@ 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 = 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.append = 0; % 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 end diff --git a/NIfTIMRSRead.m b/NIfTIMRSRead.m index d90b719..0b2db68 100644 --- a/NIfTIMRSRead.m +++ b/NIfTIMRSRead.m @@ -42,7 +42,7 @@ MRS_struct.p.LarmorFreq(ii) = hdr_ext.SpectrometerFrequency; MRS_struct.p.sw(ii) = 1/hdr.pixdim(5); -if strcmp(hdr_ext.Manufacturer,'GE') +if strcmp(hdr_ext.Manufacturer, 'GE') MRS_struct.p.TE(ii) = hdr_ext.EchoTime / 1e3; MRS_struct.p.TR(ii) = hdr_ext.RepetitionTime / 1e3; else @@ -80,7 +80,7 @@ fids_w = double(nii_w.img); MRS_struct.p.sw_water(ii) = 1/hdr_w.pixdim(5); - if strcmp(hdr_w_ext.Manufacturer,'GE') + if strcmp(hdr_w_ext.Manufacturer, 'GE') MRS_struct.p.TE_water(ii) = hdr_w_ext.EchoTime / 1e3; MRS_struct.p.TR_water(ii) = hdr_w_ext.RepetitionTime / 1e3; else @@ -111,36 +111,43 @@ if dims.coils > 0 + % Combine coils using generalized least squares method (An et al., + % JMRI, 2013, doi:10.1002/jmri.23941); the noise covariance matrix is + % more optionally estimated by using all averages as suggested by + % Rodgers & Robson (MRM, 2010, doi:10.1002/mrm.22230) + if nargin == 3 - [nCh, nPts, nReps] = size(fids_w); - noise_pts = false(1,nPts); + [nCh, nPts, nReps] = size(fids_w); + noise_pts = false(1,nPts); noise_pts(ceil(0.75*nPts):end) = true; - noise_pts = repmat(noise_pts, [1 nReps]); - tmp_fids_w = reshape(fids_w, [nCh nPts*nReps]); + noise_pts = repmat(noise_pts, [1 nReps]); + tmp_fids_w = reshape(fids_w, [nCh nPts*nReps]); - e = tmp_fids_w(:,noise_pts); - Psi = e*e'; + e = tmp_fids_w(:,noise_pts); + Psi = e*e'; fids_w_avg = mean(fids_w,3); - S = fids_w_avg(:,1); - w = (S'*(Psi\S))^-1 * S' / Psi; - fids_w = w.' .* fids_w; + S = fids_w_avg(:,1); + w = (S'*(Psi\S))^-1 * S' / Psi; + fids_w = w.' .* fids_w; + MRS_struct.fids.data_water = mean(squeeze(sum(fids_w,1)),2); end - [nCh, nPts, nReps] = size(fids); - noise_pts = false(1,nPts); + [nCh, nPts, nReps] = size(fids); + noise_pts = false(1,nPts); noise_pts(ceil(0.75*nPts):end) = true; - noise_pts = repmat(noise_pts, [1 nReps]); - tmp_fids = reshape(fids, [nCh nPts*nReps]); + noise_pts = repmat(noise_pts, [1 nReps]); + tmp_fids = reshape(fids, [nCh nPts*nReps]); - e = tmp_fids(:,noise_pts); - Psi = e*e'; + e = tmp_fids(:,noise_pts); + Psi = e*e'; if nargin == 2 fids_avg = mean(fids,3); - S = fids_avg(:,1); + S = fids_avg(:,1); end - w = (S'*(Psi\S))^-1 * S' / Psi; + w = (S'*(Psi\S))^-1 * S' / Psi; fids = w.' .* fids; + MRS_struct.fids.data = squeeze(sum(fids,1)); else @@ -168,7 +175,37 @@ end end - case 'Phillips' + case 'Philips' + % Undo phase cycling + corrph = repmat([-1 1], [1 size(MRS_struct.fids.data,2)/2]); + corrph = repmat(corrph, [size(MRS_struct.fids.data,1) 1]); + MRS_struct.fids.data = MRS_struct.fids.data .* corrph; + + % Re-introduce initial phase step + if MRS_struct.p.HERMES +% if strcmp(MRS_struct.p.ON_OFF_order,'offfirst') + phi = repelem(conj(MRS_struct.fids.data(1,2:2:end)) ./ abs(MRS_struct.fids.data(1,2:2:end)),2); +% elseif strcmp(MRS_struct.p.ON_OFF_order,'onfirst') +% ind1 = sort([1:4:size(MRS_struct.fids.data,2) 2:4:size(MRS_struct.fids.data,2)]); +% ind2 = sort([3:4:size(MRS_struct.fids.data,2) 4:4:size(MRS_struct.fids.data,2)]); +% phi(ind1) = repelem(conj(MRS_struct.fids.data(1,1:4:end)) ./ abs(MRS_struct.fids.data(1,1:4:end)),2); +% phi(ind2) = repelem(conj(MRS_struct.fids.data(1,4:4:end)) ./ abs(MRS_struct.fids.data(1,4:4:end)),2); +% end + MRS_struct.fids.data = MRS_struct.fids.data .* repmat(phi, [MRS_struct.p.npoints(ii) 1]); + else + if strcmp(MRS_struct.p.target{1}, 'GSH') + MRS_struct.fids.data = MRS_struct.fids.data .* ... + repmat(conj(mean(MRS_struct.fids.data(1,:))) ./ abs(mean(MRS_struct.fids.data(1,:))), size(MRS_struct.fids.data)); + else + MRS_struct.fids.data = MRS_struct.fids.data .* ... + repmat(conj(MRS_struct.fids.data(1,:)) ./ abs(MRS_struct.fids.data(1,:)), [MRS_struct.p.npoints(ii) 1]); + end + end + + if nargin == 3 + MRS_struct.fids.data_water = MRS_struct.fids.data_water .* ... + conj(MRS_struct.fids.data_water(1)) ./ abs(MRS_struct.fids.data_water(1)); + end end diff --git a/PhilipsRead_data.m b/PhilipsRead_data.m index 2c75e76..5300efb 100644 --- a/PhilipsRead_data.m +++ b/PhilipsRead_data.m @@ -23,10 +23,10 @@ ii = MRS_struct.ii; % Get the appropriate _act.spar -spar_files = dir('*_act.spar'); % it's automatically case-insensitive +spar_files = dir('*_act*.spar'); % it's automatically case-insensitive % If more than one .spar file is found, select one if isempty(spar_files) - spar_files = dir('*_act.SPAR'); % it's automatically case-insensitive + spar_files = dir('*_act*.SPAR'); % it's automatically case-insensitive end % If more than one .spar file is found, select one if length(spar_files) > 1 diff --git a/SavePDF.m b/SavePDF.m index a84af23..b04f585 100644 --- a/SavePDF.m +++ b/SavePDF.m @@ -1,8 +1,8 @@ function run_count = SavePDF(h, MRS_struct, ii, jj, kk, vox, module, run_count) % Gannet logo -Gannet_logo = fullfile(fileparts(which('GannetLoad')), 'Gannet3_logo.png'); -I = imread(Gannet_logo, 'png', 'BackgroundColor', [1 1 1]); +Gannet_logo = fullfile(fileparts(which('GannetLoad')), 'Gannet3_logo.jpg'); +I = imread(Gannet_logo); axes('Position', [0.825, 0.05, 0.125, 0.125]); imshow(I); text(0.9, 0, MRS_struct.version.Gannet, 'Units', 'normalized', 'FontName', 'Arial', 'FontSize', 14, 'FontWeight', 'bold', 'HorizontalAlignment', 'left'); @@ -21,7 +21,7 @@ d.h = 1-d.b; axes('Position', [d.l d.b d.w d.h]); text(0.0075, 0, ['Batch file: ' num2str(ii) ' of ' num2str(MRS_struct.p.numScans)], 'FontName', 'Arial', 'FontSize', 11, 'HorizontalAlignment', 'left'); -text(0.9925, 0, datestr(clock, 'dd-mmm-yyyy HH:MM:SS'), 'FontName', 'Arial', 'FontSize', 11, 'HorizontalAlignment', 'right'); +text(0.9925, 0, datestr(clock, 'dd-mmm-yyyy HH:MM:SS'), 'FontName', 'Arial', 'FontSize', 11, 'HorizontalAlignment', 'right'); %#ok<*CLOCK,*DATST> axis off; if any(strcmp(listfonts, 'Arial')) @@ -32,17 +32,17 @@ % If export_fig is installed, export PDF using it if MRS_struct.p.append && ~isempty(fileparts(which('export_fig'))) - + scr_sz = get(0,'ScreenSize'); fig_w = 11*72; fig_h = 8.5*72; set(gcf, 'Units', 'Pixels', 'Position', [(scr_sz(3)-fig_w)/2, (scr_sz(4)-fig_h)/2, fig_w, fig_h]); - + % Create output dir if ~exist(fullfile(pwd, 'Gannet_output'), 'dir') mkdir(fullfile(pwd, 'Gannet_output')); end - + pdf_name = fullfile(pwd, 'Gannet_output', [module '.pdf']); if exist(pdf_name, 'file') && (ii + jj) == 2 run_count = 1; @@ -58,18 +58,30 @@ elseif (ii + jj) > 2 && run_count > 0 pdf_name = fullfile(pwd, 'Gannet_output', [module '-' num2str(run_count) '.pdf']); end - - export_fig(pdf_name, '-pdf', '-painters', '-append', '-nocrop', '-nofontswap', '-silent', h); - + + try + export_fig(pdf_name, '-pdf', '-painters', '-append', '-nocrop', '-nofontswap', '-silent', h); + catch ME + switch ME.identifier + case 'MATLAB:UndefinedFunction' + error(['Cannot find the function ''export_fig.m''. ' ... + 'Please ensure that you have added the export_fig ', ... + 'folder in the main Gannet folder to your MATLAB ', ... + 'search path.']); + otherwise + rethrow(ME); + end + end + else - + set(h, 'PaperUnits', 'inches', 'PaperSize', [11 8.5], 'PaperPosition', [0 0 11 8.5]); - + % Create output folder if ~exist(fullfile(pwd, [module '_output']),'dir') mkdir(fullfile(pwd, [module '_output'])); end - + % For Philips .data if strcmp(MRS_struct.p.vendor, 'Philips_data') fullpath = MRS_struct.metabfile{1,ii}; @@ -77,20 +89,23 @@ fullpath = regexprep(fullpath, '\', '_'); fullpath = regexprep(fullpath, '/', '_'); end - + if strcmp(MRS_struct.p.vendor, 'Siemens_rda') [~,metabfile_nopath] = fileparts(MRS_struct.metabfile{1,ii*2-1}); else - [~,metabfile_nopath] = fileparts(MRS_struct.metabfile{1,ii}); + [~,metabfile_nopath,ext] = fileparts(MRS_struct.metabfile{1,ii}); + if strcmpi(ext,'.gz') + metabfile_nopath(end-3:end) = []; + end end - + module2 = lower(module); module2(1:6) = []; - + if strcmp(module2, 'fit') && MRS_struct.p.HERMES module2 = [MRS_struct.p.target{jj} '_' module2]; end - + if strcmp(MRS_struct.p.vendor, 'Philips_data') if isfield(MRS_struct.p, 'trimmed_avgs') pdf_name = fullfile(pwd, [module '_output'], [fullpath '_' vox{kk} '_' module2 '_' num2str(MRS_struct.p.Navg(ii)) '_avgs.pdf']); @@ -104,9 +119,9 @@ pdf_name = fullfile(pwd, [module '_output'], [metabfile_nopath '_' vox{kk} '_' module2 '.pdf']); end end - + saveas(h, pdf_name); - + end diff --git a/SiemensRead.m b/SiemensRead.m index 7699c21..c9a9c88 100644 --- a/SiemensRead.m +++ b/SiemensRead.m @@ -10,22 +10,22 @@ tline = fgets(fid); while isempty(strfind(tline, head_end_text)) - + tline = fgets(fid); - + if isempty(strfind(tline, head_start_text)) + isempty(strfind(tline, head_end_text)) == 2 - + % Store this data in the appropriate format occurence_of_colon = strfind(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; ']); + 'SeriesTime', 'SeriesDescription', 'ProtocolName', 'PatientPosition', 'ModelName', 'StationName', 'InstitutionName', ... + 'DeviceSerialNumber', 'InstanceDate', 'InstanceTime', 'InstanceComments', 'SequenceName', 'SequenceDescription', 'Nucleus', ... + 'TransmitCoil'} + rda.(variable) = value; case {'PatientSex'} % Sex converter! (int to M,F,U) switch value @@ -33,17 +33,16 @@ rda.sex = 'Unknown'; case 1 rda.sex = 'Male'; - case 2 + case 2 rda.sex = 'Female'; end - case {'SeriesNumber', 'InstanceNumber', 'AcquisitionNumber', 'NumOfPhaseEncodingSteps', 'NumberOfRows', 'NumberOfColumns', 'VectorSize'} %Integers - eval(['rda.', variable, ' = str2double(value); ']); + rda.(variable) = str2double(value); case {'PatientWeight', 'TR', 'TE', 'TM', 'DwellTime', 'NumberOfAverages', 'MRFrequency', 'MagneticFieldStrength', 'FlipAngle', ... 'SliceThickness', 'FoVHeight', 'FoVWidth', 'PercentOfRectFoV', 'PixelSpacingRow', 'PixelSpacingCol'} %Floats - eval(['rda.', variable, ' = str2double(value); ']); + rda.(variable) = str2double(value); case {'SoftwareVersion[0]'} rda.software_version = value; case {'CSIMatrixSize[0]'} @@ -70,16 +69,16 @@ rda.ColumnVector(2) = str2double(value); case {'ColumnVector[2]'} rda.ColumnVector(3) = str2double(value); - 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 + end %%RE 110726 Take the used bits of the header info @@ -143,22 +142,22 @@ tline = fgets(fid); while isempty(strfind(tline, head_end_text)) - + tline = fgets(fid); - - if isempty(strfind(tline, head_start_text)) + isempty(strfind(tline, head_end_text)) == 2 - - % Store this data in the appropriate format + + if isempty(strfind(tline, head_start_text)) + isempty(strfind(tline, head_end_text)) == 2 + + % Store this data in the appropriate format occurence_of_colon = strfind(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; ']); + 'SeriesTime', 'SeriesDescription', 'ProtocolName', 'PatientPosition', 'ModelName', 'StationName', 'InstitutionName', ... + 'DeviceSerialNumber', 'InstanceDate', 'InstanceTime', 'InstanceComments', 'SequenceName', 'SequenceDescription', 'Nucleus', ... + 'TransmitCoil'} + rda.(variable) = value; case {'PatientSex'} % Sex converter! (int to M,F,U) switch value @@ -166,17 +165,16 @@ rda.sex = 'Unknown'; case 1 rda.sex = 'Male'; - case 2 + case 2 rda.sex = 'Female'; end - case {'SeriesNumber', 'InstanceNumber', 'AcquisitionNumber', 'NumOfPhaseEncodingSteps', 'NumberOfRows', 'NumberOfColumns', 'VectorSize'} %Integers - eval(['rda.', variable, ' = str2double(value); ']); + rda.(variable) = str2double(value); case {'PatientWeight', 'TR', 'TE', 'TM', 'DwellTime', 'NumberOfAverages', 'MRFrequency', 'MagneticFieldStrength', 'FlipAngle', ... - 'SliceThickness', 'FoVHeight', 'FoVWidth', 'PercentOfRectFoV', 'PixelSpacingRow', 'PixelSpacingCol'} + 'SliceThickness', 'FoVHeight', 'FoVWidth', 'PercentOfRectFoV', 'PixelSpacingRow', 'PixelSpacingCol'} %Floats - eval(['rda.', variable, ' = str2double(value); ']); + rda.(variable) = str2double(value); case {'SoftwareVersion[0]'} rda.software_version = value; case {'CSIMatrixSize[0]'} @@ -203,17 +201,16 @@ rda.ColumnVector(2) = str2double(value); case {'ColumnVector[2]'} rda.ColumnVector(3) = str2double(value); - 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 %%RE 110726 Take the used bits of the header info @@ -243,32 +240,33 @@ MRS_struct.fids.ondata = hmm_complex; MRS_struct.fids.data = [MRS_struct.fids.ondata; MRS_struct.fids.offdata].'; -if nargin==4 +if nargin == 4 + %%%Now load in the Water data rda_filename=water_filename; %This is generally file3 fid = fopen(rda_filename); - + head_start_text = '>>> Begin of header <<<'; head_end_text = '>>> End of header <<<'; - + tline = fgets(fid); - + while (isempty(strfind(tline, head_end_text))) %#ok<*STREMP> - + tline = fgets(fid); - - if isempty(strfind(tline, head_start_text)) + isempty(strfind(tline, head_end_text)) == 2 + + if isempty(strfind(tline, head_start_text)) + isempty(strfind(tline, head_end_text)) == 2 % Store this data in the appropriate format occurence_of_colon = strfind(':',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; ']); + 'SeriesTime', 'SeriesDescription', 'ProtocolName', 'PatientPosition', 'ModelName', 'StationName', 'InstitutionName', ... + 'DeviceSerialNumber', 'InstanceDate', 'InstanceTime', 'InstanceComments', 'SequenceName', 'SequenceDescription', 'Nucleus', ... + 'TransmitCoil'} + rda.(variable) = str2double(value); case {'PatientSex'} % Sex converter! (int to M,F,U) switch value @@ -276,17 +274,16 @@ rda.sex = 'Unknown'; case 1 rda.sex = 'Male'; - case 2 + case 2 rda.sex = 'Female'; end - case {'SeriesNumber', 'InstanceNumber', 'AcquisitionNumber', 'NumOfPhaseEncodingSteps', 'NumberOfRows', 'NumberOfColumns', 'VectorSize'} %Integers - eval(['rda.', variable, ' = str2double(value); ']); + rda.(variable) = str2double(value); case {'PatientWeight', 'TR', 'TE', 'TM', 'DwellTime', 'NumberOfAverages', 'MRFrequency', 'MagneticFieldStrength', 'FlipAngle', ... - 'SliceThickness', 'FoVHeight', 'FoVWidth', 'PercentOfRectFoV', 'PixelSpacingRow', 'PixelSpacingCol'} + 'SliceThickness', 'FoVHeight', 'FoVWidth', 'PercentOfRectFoV', 'PixelSpacingRow', 'PixelSpacingCol'} %Floats - eval(['rda.', variable, ' = str2double(value); ']); + rda.(variable) = str2double(value); case {'SoftwareVersion[0]'} rda.software_version = value; case {'CSIMatrixSize[0]'} @@ -313,19 +310,17 @@ rda.ColumnVector(2) = str2double(value); case {'ColumnVector[2]'} rda.ColumnVector(3) = str2double(value); - 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 - + %%RE 110726 Take the used bits of the header info MRS_struct.p.LarmorFreq(MRS_struct.ii) = rda.MRFrequency; MRS_struct.p.npoints(MRS_struct.ii) = rda.VectorSize; @@ -334,24 +329,24 @@ % % Siemens documentation suggests that the data should be in a double complex format (8bytes for real, and 8 for imaginary?) % - + %bytes_per_point = 16; 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 hmm_complex = complex(hmm(1,:,:,:,:), hmm(2,:,:,:,:)); - + %RE 110726 This is the complex time domain data MRS_struct.fids.data_water = hmm_complex.'; - + end end diff --git a/SiemensTwixRead.m b/SiemensTwixRead.m index c8fc6ae..aaca686 100644 --- a/SiemensTwixRead.m +++ b/SiemensTwixRead.m @@ -41,7 +41,7 @@ % 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. +% 2022-10-20: Added support for XA30 sequence provided by JHU. ii = MRS_struct.ii; @@ -85,7 +85,7 @@ % If additional data points have been acquired before the echo starts, % remove these here. -MetabData = MetabData(:,(MRS_struct.p.pointsBeforeEcho+1):end,:); +MetabData = MetabData(:,(MRS_struct.p.pointsBeforeEcho+1):end,:); MRS_struct.p.npoints(ii) = MRS_struct.p.npoints(ii) - MRS_struct.p.pointsBeforeEcho; % Undo phase cycling @@ -100,6 +100,7 @@ % 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); MRS_struct.p.pointsBeforeEcho_water = WaterHeader.pointsBeforeEcho; MRS_struct.p.sw_water(ii) = 1/WaterHeader.dwellTime; @@ -118,7 +119,7 @@ % If additional data points have been acquired before the echo starts, % remove these here. - WaterData = WaterData(:,(MRS_struct.p.pointsBeforeEcho_water+1):end,:); + 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 @@ -130,45 +131,53 @@ 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; - + % Combine coils using generalized least squares method (An et al., + % JMRI, 2013, doi:10.1002/jmri.23941); the noise covariance matrix is + % more optionally estimated by using all averages as suggested by + % Rodgers & Robson (MRM, 2010, doi:10.1002/mrm.22230) + + [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]); + tmp_fids_w = reshape(WaterData, [nCh nPts*nReps]); + + e = tmp_fids_w(:,noise_pts); + Psi = e*e'; + fids_w_avg = mean(WaterData,3); + [~,ind] = max(abs(fids_w_avg),[],2); + ind = mode(ind); + S = fids_w_avg(:,ind); + w = (S'*(Psi\S))^-1 * S' / Psi; + WaterData = w.' .* WaterData; + MRS_struct.fids.data_water = double(mean(conj(squeeze(sum(WaterData,1))),2)); + end -% Coil combination of metabolite data -% If water data has been acquired, use the phase information from it to -% phase the metabolite data -if nargin == 3 - % 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). - [~,ind] = max(abs(mean(MetabData,3)),[],2); - ind = mode(ind); - firstpoint = mean(conj(MetabData(:,ind,:)),3); - channels_scale = squeeze(sqrt(sum(firstpoint .* conj(firstpoint)))); - firstpoint = repmat(firstpoint, [1 MRS_struct.p.npoints(ii) MRS_struct.p.nrows(ii)])/channels_scale; - MetabData = MetabData .* firstpoint; - MetabData = conj(squeeze(sum(MetabData,1))); - MRS_struct.fids.data = double(MetabData); +% Combine coils using generalized least squares method (An et al., JMRI, +% 2013, doi:10.1002/jmri.23941); the noise covariance matrix is more +% optionally estimated by using all averages as suggested by Rodgers & +% Robson (MRM, 2010, doi:10.1002/mrm.22230) + +[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]); +tmp_fids = reshape(MetabData, [nCh nPts*nReps]); + +e = tmp_fids(:,noise_pts); +Psi = e*e'; +if nargin == 2 + fids_avg = mean(MetabData,3); + [~,ind] = max(abs(fids_avg),[],2); + ind = mode(ind); + S = fids_avg(:,ind); end +w = (S'*(Psi\S))^-1 * S' / Psi; +MetabData = w.' .* MetabData; + +MRS_struct.fids.data = double(conj(squeeze(sum(MetabData,1)))); end @@ -336,7 +345,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)]); -elseif any(strcmp(TwixHeader.seqtype,{'MEGA-PRESS','MEGA-sLASER'})) % SH 20191213 +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. diff --git a/UpdateGannet.m b/UpdateGannet.m index 54d19c4..031dd72 100644 --- a/UpdateGannet.m +++ b/UpdateGannet.m @@ -48,7 +48,7 @@ zipURL = 'https://github.com/markmikkelsen/Gannet/archive/refs/heads/main.zip'; targetFolder = fullfile(pwd, ['tmp_' randsample(['A':'Z','0':'9'],5)]); mkdir(targetFolder); -targetFilename = fullfile(targetFolder, datestr(now,'yyyy-mm-dd.zip')); +targetFilename = fullfile(targetFolder, datestr(now,'yyyy-mm-dd.zip')); %#ok<*TNOW1,*DATST> websave(targetFilename, zipURL); newFilenames = unzip(targetFilename, targetFolder); diff --git a/VersionCheck.m b/VersionCheck.m index 9fa49d1..8d5d8ec 100644 --- a/VersionCheck.m +++ b/VersionCheck.m @@ -6,6 +6,7 @@ try java.net.InetAddress.getByName('www.google.com'); catch + warning('No internet connection. Skipping version check.'); return end @@ -29,7 +30,7 @@ end newVersionAvailable = 0; -if nargin < 2 || isempty(lastCheckTime) || etime(clock, lastCheckTime) > 86.4e3 +if nargin < 2 || isempty(lastCheckTime) || etime(clock, lastCheckTime) > 86.4e3 %#ok<*DETIM> url = 'https://raw.githubusercontent.com/markmikkelsen/Gannet/main/GannetLoad.m'; str = readURL(url); expression = '(?MRS_struct.version.Gannet = )''(?.*?)'''; @@ -46,7 +47,7 @@ fprintf(msg, latestVersion, currentVersion); end end - lastCheckTime = clock; + lastCheckTime = clock; %#ok<*CLOCK> end if nargout == 1 diff --git a/mapVBVD_Gannet.m b/mapVBVD_Gannet.m index d870541..ae84bd5 100644 --- a/mapVBVD_Gannet.m +++ b/mapVBVD_Gannet.m @@ -4,9 +4,9 @@ % requires twix_map_obj.m & read_twix_hdr.m % % -% Author: Philipp Ehses (philipp.ehses@tuebingen.mpg.de) -% -% +% Author: Philipp Ehses (philipp.ehses@dzne.de) +% +% % Philipp Ehses 11.02.07, original version % [..] % Philipp Ehses 22.03.11, port to VD11 @@ -33,7 +33,7 @@ % * For parsing, use an as-minimalistic-as-possible loop % to gather all mdhs in binary form. They are all stored % in one array "mdh_blob". -% * Translation of mdhs from binary to struct is vectorized +% * Translation of mdhs from binary to struct is vectorized % and almost instant. It's done in evalMDH(), replacing % evalMDHvb() and evalMDHvd(), no more freads inside! % * vectorized readMDH(), quasi-instant @@ -44,15 +44,15 @@ % * => Parsing speed improved by factor 3...7 or so % * Speed increase for reading data, esp. when slicing, % os-removal or reflected lines. Also for random acquisitions. -% Jonas Bause, 18.11.16 receiver phase for ramp-sampling fixed, now takes into account +% Jonas Bause, 18.11.16 receiver phase for ramp-sampling fixed, now takes into account % Chris Mirkes & PE offcenter shifts in readout direction -% -% +% +% % Input: -% +% % filename or simply meas. id, e.g. mapVBVD(122) (if file is in same path) % optional arguments (see below) -% +% % Output: twix_obj structure with elements (if available): % .image: image scan % .noise: for noise scan @@ -67,8 +67,8 @@ % .refscanPSRef1: phase stab. ref scan for reference data % .RTfeedback: realtime feedback data % .vop: vop rf data -% -% +% +% % The raw data can be obtained by calling e.g. twix_obj.image() or for % squeezed data twix_obj.image{''} (the '' are needed due to a limitation % of matlab's overloading capabilities). @@ -87,8 +87,8 @@ % with a semicolon. Thus, a call to e.g. twix_obj.image.NLin will produce % no output with or without semicolon termination. 'a = twix_obj.image.NLin' % will however produce the expected result. -% -% +% +% % Order of raw data: % 1) Columns % 2) Channels/Coils @@ -106,10 +106,10 @@ % 14) Idc % 15) Idd % 16) Ide -% -% +% +% % Optional parameters/flags: -% +% % removeOS: removes oversampling (factor 2) in read direction % doAverage: performs average (resulting avg-dim has thus size 1) % ignoreSeg: ignores segment mdh index (works basically the same as @@ -117,13 +117,13 @@ % rampSampRegrid optional on-the-fly regridding of ramp-sampled readout % doRawDataCorrect: enables raw data correction if used in the acquisition % (only works for VB atm) -% +% % These flags can also be set/unset later, e.g "twix_obj.image.flagRemoveOS = 1" -% -% +% +% % Examples: % twix_obj = mapVBVD(measID); -% +% % % return all image-data: % image_data = twix_obj.image(); % % return all image-data with all singular dimensions removed/squeezed: @@ -135,14 +135,14 @@ % % grouped into dim 5); but work with the squeezed data order % % => use '{}' instead of '()': % image_data = twix_obj.image{:,2:6,:,:,:}; -% +% % So basically it works like regular matlab array slicing (but 'end' is % not supported; note that there are still a few bugs with array slicing). -% +% % % NEW: unsorted raw data (in acq. order): % image_data = twix_obj.image.unsorted(); % no slicing supported atm -% -% +% +% % Suppress silly editor warnings in the entire file, barking about % unused variables: %#ok<*NASGU> @@ -157,7 +157,7 @@ else if ischar(filename) % assume that complete path is given - if ~strcmpi(filename(end-3:end),'.dat') + if ~strcmpi(filename(end-3:end),'.dat') filename=[filename '.dat']; %% adds filetype ending to file end else @@ -176,7 +176,7 @@ if filesfound == 0 error(['File with meas. id ' num2str(measID) ' not found.']); elseif filesfound > 1 - disp(['Multiple files with meas. id ' num2str(measID) ' found. Choosing first occurence.']); + %disp(['Multiple files with meas. id ' num2str(measID) ' found. Choosing first occurence.']); end end end @@ -288,8 +288,8 @@ % start of actual measurement data (sans header) fseek(fid,0,'bof'); -firstInt = fread(fid,1,'uint32'); -secondInt = fread(fid,1,'uint32'); +firstInt = fread(fid, 1, 'uint32'); +secondInt = fread(fid, 1, 'uint32'); % lazy software version check (VB or VD?) if and(firstInt < 10000, secondInt <= 64) @@ -304,7 +304,7 @@ measLength = cell(1, NScans); for k=1:NScans measOffset{k} = fread(fid,1,'uint64'); - measLength{k} = fread(fid,1,'uint64'); + measLength{k} = fread(fid,1,'uint64'); fseek(fid, 152 - 16, 'cof'); end else @@ -348,7 +348,7 @@ % end; if (mod(length(rawfactors),2) ~= 0) error('Error reading rawfactors'); - end; + end %note the transpose, this makes the vector %multiplication during the read easier arg.rawDataCorrectionFactors = rawfactors(1:2:end).' + 1i*rawfactors(2:2:end).'; @@ -357,7 +357,7 @@ end end end - %disp('Read raw data correction factors'); + %disp('Read raw data correction factors'); % MM (200610) end % data will be read in two steps (two while loops): @@ -374,8 +374,12 @@ % read header and calculate regridding (optional) rstraj = []; + isInplacePATref = false; if arg.bReadHeader [twix_obj{s}.hdr, rstraj] = read_twix_hdr_Gannet(fid); + if isfield(twix_obj{s}.hdr.MeasYaps, 'sPat') + isInplacePATref = str2double(twix_obj{s}.hdr.MeasYaps.sPat.ucRefScanMode(end))==2; + end end % declare data objects: @@ -392,10 +396,10 @@ twix_obj{s}.refscanPSRef1 = twix_map_obj_Gannet(arg,'refscan_phasestab_ref1',filename,version,rstraj); twix_obj{s}.RTfeedback = twix_map_obj_Gannet(arg,'rtfeedback',filename,version,rstraj); twix_obj{s}.vop = twix_map_obj_Gannet(arg,'vop',filename,version); % tx-array rf pulses - + % print reader version information if s==1 - %fprintf('Reader version: %d', twix_obj{s}.image.readerVersion); + %fprintf('Reader version: %d', twix_obj{s}.image.readerVersion); % MM (200610) isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0; if isOctave date_str = ctime(twix_obj{s}.image.readerVersion); @@ -403,15 +407,15 @@ else date_str = datetime(twix_obj{s}.image.readerVersion, 'ConvertFrom','posixtime'); end - %fprintf(' (UTC: %s)\n', char(date_str)); + %fprintf(' (UTC: %s)\n', char(date_str)); % MM (200610) end - + % jump to first mdh cPos = cPos + hdr_len; fseek( fid, cPos, 'bof' ); % find all mdhs and save them in binary form, first: - %fprintf('Scan %d/%d, read all mdhs:\n', s, NScans ) + %fprintf('Scan %d/%d, read all mdhs:\n', s, NScans ) % MM (200610) [mdh_blob, filePos, isEOF] = loop_mdh_read( fid, version, NScans, s, measOffset{s}, measLength{s}); % uint8; size: [ byteMDH Nmeas ] @@ -419,7 +423,7 @@ filePos( end ) = []; % get mdhs and masks for each scan, no matter if noise, image, RTfeedback etc: - [mdh, mask] = evalMDH( mdh_blob, version ); % this is quasi-instant (< 1s) :-) + [mdh, mask] = evalMDH( mdh_blob, version, isInplacePATref); % this is quasi-instant (< 1s) :-) clear mdh_blob % Assign mdhs to their respective scans and parse it in the correct twix objects. @@ -443,9 +447,9 @@ if arg.bReadRefScan clear tmpMdh isCurrScan = ( mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN )... - & ~( mask.MDH_PHASCOR | mask.MDH_PHASESTABSCAN | ... - mask.MDH_REFPHASESTABSCAN | ... - mask.MDH_RTFEEDBACK | mask.MDH_HPFEEDBACK); + & ~( mask.MDH_PHASCOR | mask.MDH_PHASESTABSCAN | ... + mask.MDH_REFPHASESTABSCAN | ... + mask.MDH_RTFEEDBACK | mask.MDH_HPFEEDBACK); for f = fieldnames( mdh ).' tmpMdh.(f{1}) = mdh.(f{1})( isCurrScan, : ); end @@ -468,14 +472,14 @@ end if arg.bReadPCScan % logic really correct? - + clear tmpMdh isCurrScan = mask.MDH_PHASCOR & ( ~mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); for f = fieldnames( mdh ).' tmpMdh.(f{1}) = mdh.(f{1})( isCurrScan, : ); end twix_obj{s}.phasecor.readMDH( tmpMdh, filePos(isCurrScan) ); - + clear tmpMdh isCurrScan = mask.MDH_PHASCOR & ( mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); for f = fieldnames( mdh ).' @@ -486,7 +490,7 @@ if arg.bReadPhaseStab clear tmpMdh isCurrScan = ( mask.MDH_PHASESTABSCAN & ~mask.MDH_REFPHASESTABSCAN ) ... - & (~mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); + & (~mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); for f = fieldnames( mdh ).' tmpMdh.(f{1}) = mdh.(f{1})( isCurrScan, : ); end @@ -494,7 +498,7 @@ clear tmpMdh isCurrScan = ( mask.MDH_PHASESTABSCAN & ~mask.MDH_REFPHASESTABSCAN ) ... - & ( mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); + & ( mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); for f = fieldnames( mdh ).' tmpMdh.(f{1}) = mdh.(f{1})( isCurrScan, : ); end @@ -502,7 +506,7 @@ clear tmpMdh isCurrScan = ( mask.MDH_REFPHASESTABSCAN & ~mask.MDH_PHASESTABSCAN ) ... - & (~mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); + & (~mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); for f = fieldnames( mdh ).' tmpMdh.(f{1}) = mdh.(f{1})( isCurrScan, : ); end @@ -510,15 +514,15 @@ clear tmpMdh isCurrScan = ( mask.MDH_REFPHASESTABSCAN & ~mask.MDH_PHASESTABSCAN ) ... - & ( mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); + & ( mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); for f = fieldnames( mdh ).' tmpMdh.(f{1}) = mdh.(f{1})( isCurrScan, : ); end twix_obj{s}.refscanPSRef0.readMDH( tmpMdh, filePos(isCurrScan) ); - + clear tmpMdh isCurrScan = ( mask.MDH_REFPHASESTABSCAN & mask.MDH_PHASESTABSCAN ) ... - & (~mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); + & (~mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); for f = fieldnames( mdh ).' tmpMdh.(f{1}) = mdh.(f{1})( isCurrScan, : ); end @@ -526,7 +530,7 @@ clear tmpMdh isCurrScan = ( mask.MDH_REFPHASESTABSCAN & mask.MDH_PHASESTABSCAN ) ... - & ( mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); + & ( mask.MDH_PATREFSCAN | mask.MDH_PATREFANDIMASCAN ); for f = fieldnames( mdh ).' tmpMdh.(f{1}) = mdh.(f{1})( isCurrScan, : ); end @@ -535,9 +539,9 @@ clear mdh tmpMdh filePos isCurrScan for scan = { 'image', 'noise', 'phasecor', 'phasestab', ... - 'phasestabRef0', 'phasestabRef1', 'refscan', ... - 'refscanPC', 'refscanPS', 'refscanPSRef0', ... - 'refscanPSRef1', 'RTfeedback', 'vop' } + 'phasestabRef0', 'phasestabRef1', 'refscan', ... + 'refscanPC', 'refscanPS', 'refscanPSRef0', ... + 'refscanPSRef1', 'RTfeedback', 'vop' } f = scan{1}; % remove unused fields @@ -555,7 +559,7 @@ end % NScans loop -if NScans == 1 +if strcmp(version, 'vb') twix_obj = twix_obj{1}; end @@ -582,173 +586,174 @@ % or threads using a java class (probably faster + no toolbox): % http://undocumentedmatlab.com/blog/explicit-multi-threading-in-matlab-part1 - switch version - case 'vb' - isVD = false; - byteMDH = 128; - case 'vd' - isVD = true; - byteMDH = 184; - szScanHeader = 192; % [bytes] - szChannelHeader = 32; % [bytes] - otherwise - % arbitrary assumptions: - isVD = false; - byteMDH = 128; - warning( [mfilename() ':UnknownVer'], 'Software version "%s" is not supported.', version ); - end +switch version + case 'vb' + isVD = false; + byteMDH = 128; + case 'vd' + isVD = true; + byteMDH = 184; + szScanHeader = 192; % [bytes] + szChannelHeader = 32; % [bytes] + otherwise + % arbitrary assumptions: + isVD = false; + byteMDH = 128; + warning( [mfilename() ':UnknownVer'], 'Software version "%s" is not supported.', version ); +end - cPos = ftell(fid); - n_acq = 0; - allocSize = 4096; - ulDMALength = byteMDH; - isEOF = false; - last_progress = 0; +cPos = ftell(fid); +n_acq = 0; +allocSize = 4096; +ulDMALength = byteMDH; +isEOF = false; +last_progress = 0; - mdh_blob = zeros( byteMDH, 0, 'uint8' ); - szBlob = size( mdh_blob, 2 ); - filePos = zeros(0, 1, class(cPos)); % avoid bug in Matlab 2013b: https://scivision.co/matlab-fseek-bug-with-uint64-offset/ +mdh_blob = zeros( byteMDH, 0, 'uint8' ); +szBlob = size( mdh_blob, 2 ); +filePos = zeros(0, 1, class(cPos)); % avoid bug in Matlab 2013b: https://scivision.co/matlab-fseek-bug-with-uint64-offset/ - fseek(fid,cPos,'bof'); +fseek(fid,cPos,'bof'); - % ====================================== - % constants and conditional variables - % ====================================== - bit_0 = uint8(2^0); - bit_5 = uint8(2^5); - mdhStart = 1-byteMDH; - - u8_000 = zeros( 3, 1, 'uint8'); % for comparison with data_u8(1:3) - - % 20 fill bytes in VD (21:40) - evIdx = uint8( 21 + 20*isVD); % 1st byte of evalInfoMask - dmaIdx = uint8((29:32) + 20*isVD); % to correct DMA length using NCol and NCha - if isVD - dmaOff = szScanHeader; - dmaSkip = szChannelHeader; - else - dmaOff = 0; - dmaSkip = byteMDH; - end - % ====================================== +% ====================================== +% constants and conditional variables +% ====================================== +bit_0 = uint8(2^0); +bit_5 = uint8(2^5); +mdhStart = 1-byteMDH; - isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0; - if isOctave % octave does not support a cancel button - h = waitbar(0,'','Name', sprintf('Reading Scan ID %d/%d', scan, Nscans)); - else - % MM (200610) - %h = waitbar(0,'','Name', sprintf('Reading Scan ID %d/%d', scan, Nscans),... - % 'CreateCancelBtn',... - % 'setappdata(gcbf,''canceling'',1)'); - %setappdata(h,'canceling',0) +u8_000 = zeros( 3, 1, 'uint8'); % for comparison with data_u8(1:3) + +% 20 fill bytes in VD (21:40) +evIdx = uint8( 21 + 20*isVD); % 1st byte of evalInfoMask +dmaIdx = uint8((29:32) + 20*isVD); % to correct DMA length using NCol and NCha +if isVD + dmaOff = szScanHeader; + dmaSkip = szChannelHeader; +else + dmaOff = 0; + dmaSkip = byteMDH; +end +% ====================================== + +% MM (200610) +%isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0; +%if isOctave % octave does not support a cancel button +% h = waitbar(0,'','Name', sprintf('Reading Scan ID %d/%d', scan, Nscans)); +%else +% +% h = waitbar(0,'','Name', sprintf('Reading Scan ID %d/%d', scan, Nscans),... +% 'CreateCancelBtn',... +% 'setappdata(gcbf,''canceling'',1)'); +% setappdata(h,'canceling',0) +%end + +t0 = tic; +while true + % Read mdh as binary (uint8) and evaluate as little as possible to know... + % ... where the next mdh is (ulDMALength / ushSamplesInScan & ushUsedChannels) + % ... whether it is only for sync (MDH_SYNCDATA) + % ... whether it is the last one (MDH_ACQEND) + % evalMDH() contains the correct and readable code for all mdh entries. + + try + % read everything and cut out the mdh + data_u8 = fread( fid, ulDMALength, 'uint8=>uint8' ); + data_u8 = data_u8( mdhStart+end : end ); + catch exc + warning( [mfilename() ':UnxpctdEOF'], ... + [ '\nAn unexpected read error occurred at this byte offset: %d (%g GiB)\n'... + 'Will stop reading now.\n' ... + '=== MATLABs error message ================\n' ... + exc.message ... + '\n=== end of error =========================\n' ... + ], cPos, cPos/1024^3 ) + isEOF = true; + break end - t0 = tic; - while true - % Read mdh as binary (uint8) and evaluate as little as possible to know... - % ... where the next mdh is (ulDMALength / ushSamplesInScan & ushUsedChannels) - % ... whether it is only for sync (MDH_SYNCDATA) - % ... whether it is the last one (MDH_ACQEND) - % evalMDH() contains the correct and readable code for all mdh entries. - - try - % read everything and cut out the mdh - data_u8 = fread( fid, ulDMALength, 'uint8=>uint8' ); - data_u8 = data_u8( mdhStart+end : end ); - catch exc - warning( [mfilename() ':UnxpctdEOF'], ... - [ '\nAn unexpected read error occurred at this byte offset: %d (%g GiB)\n'... - 'Will stop reading now.\n' ... - '=== MATLABs error message ================\n' ... - exc.message ... - '\n=== end of error =========================\n' ... - ], cPos, cPos/1024^3 ) - isEOF = true; - break - end + % MM (200610) + %if ~isOctave && getappdata(h,'canceling') + % break; + %end - % MM (200610) - %if ~isOctave && getappdata(h,'canceling') - % break; - %end - - bitMask = data_u8(evIdx); % the initial 8 bit from evalInfoMask are enough - - if isequal( data_u8(1:3), u8_000 ) ... % probably ulDMALength == 0 - || bitand(bitMask, bit_0) % MDH_ACQEND - - % ok, look closer if really all *4* bytes are 0: - data_u8(4)= bitget( data_u8(4),1); % ubit24: keep only 1 bit from the 4th byte - ulDMALength = double( typecast( data_u8(1:4), 'uint32' ) ); - - if ulDMALength == 0 || bitand(bitMask, bit_0) - cPos = cPos + ulDMALength; - % jump to next full 512 bytes - if mod(cPos,512) - cPos = cPos + 512 - mod(cPos,512); - end - break; - end - end - if bitand(bitMask, bit_5) % MDH_SYNCDATA - data_u8(4)= bitget( data_u8(4),1); % ubit24: keep only 1 bit from the 4th byte - ulDMALength = double( typecast( data_u8(1:4), 'uint32' ) ); - cPos = cPos + ulDMALength; - continue - end + bitMask = data_u8(evIdx); % the initial 8 bit from evalInfoMask are enough - % pehses: the pack bit indicates that multiple ADC are packed into one - % DMA, often in EPI scans (controlled by fRTSetReadoutPackaging in IDEA) - % since this code assumes one adc (x NCha) per DMA, we have to correct - % the "DMA length" - % if mdh.ulPackBit - % it seems that the packbit is not always set correctly - NCol_NCha = double( typecast( data_u8(dmaIdx), 'uint16' ) ); % [ushSamplesInScan ushUsedChannels] - ulDMALength = dmaOff + (8*NCol_NCha(1) + dmaSkip) * NCol_NCha(2); - - n_acq = n_acq + 1; - - % grow arrays in batches - if n_acq > szBlob - mdh_blob( :, end + allocSize ) = 0; - filePos( end + allocSize ) = 0; - szBlob = size( mdh_blob, 2 ); - end - mdh_blob(:,n_acq) = data_u8; - filePos( n_acq ) = cPos; - - % MM (200610) - %progress = (cPos-measOffset)/measLength; - % - %if progress > last_progress + 0.01 - % last_progress = progress; - % elapsed_time = toc(t0); - % time_left = elapsed_time * (1/progress-1); - % progress_str = sprintf('%3.0f %% read in %4.0f s;\nestimated time left: %4.0f s', round(100*progress), elapsed_time, time_left); - % waitbar(progress, h, progress_str); - %end + if isequal( data_u8(1:3), u8_000 ) ... % probably ulDMALength == 0 + || bitand(bitMask, bit_0) % MDH_ACQEND + % ok, look closer if really all *4* bytes are 0: + data_u8(4)= bitget( data_u8(4),1); % ubit24: keep only 1 bit from the 4th byte + ulDMALength = double( typecast( data_u8(1:4), 'uint32' ) ); + + if ulDMALength == 0 || bitand(bitMask, bit_0) + cPos = cPos + ulDMALength; + % jump to next full 512 bytes + if mod(cPos,512) + cPos = cPos + 512 - mod(cPos,512); + end + break; + end + end + if bitand(bitMask, bit_5) % MDH_SYNCDATA + data_u8(4)= bitget( data_u8(4),1); % ubit24: keep only 1 bit from the 4th byte + ulDMALength = double( typecast( data_u8(1:4), 'uint32' ) ); cPos = cPos + ulDMALength; - end % while true - %delete(h); % MM (200610) - - if isEOF - n_acq = n_acq-1; % ignore the last attempt + continue end - filePos( n_acq+1 ) = cPos; % save pointer to the next scan + % pehses: the pack bit indicates that multiple ADC are packed into one + % DMA, often in EPI scans (controlled by fRTSetReadoutPackaging in IDEA) + % since this code assumes one adc (x NCha) per DMA, we have to correct + % the "DMA length" + % if mdh.ulPackBit + % it seems that the packbit is not always set correctly + NCol_NCha = double( typecast( data_u8(dmaIdx), 'uint16' ) ); % [ushSamplesInScan ushUsedChannels] + ulDMALength = dmaOff + (8*NCol_NCha(1) + dmaSkip) * NCol_NCha(2); + + n_acq = n_acq + 1; + + % grow arrays in batches + if n_acq > szBlob + mdh_blob( :, end + allocSize ) = 0; + filePos( end + allocSize ) = 0; + szBlob = size( mdh_blob, 2 ); + end + mdh_blob(:,n_acq) = data_u8; + filePos( n_acq ) = cPos; + + % MM (200610) + %progress = (cPos-measOffset)/measLength; + + %if progress > last_progress + 0.01 + % last_progress = progress; + % elapsed_time = toc(t0); + % time_left = elapsed_time * (1/progress-1); + % progress_str = sprintf('%3.0f %% read in %4.0f s;\nestimated time left: %4.0f s', round(100*progress), elapsed_time, time_left); + % waitbar(progress, h, progress_str); + %end + + cPos = cPos + ulDMALength; +end % while true +%delete(h); % MM (200610) + +if isEOF + n_acq = n_acq-1; % ignore the last attempt +end + +filePos( n_acq+1 ) = cPos; % save pointer to the next scan - % discard overallocation: - mdh_blob = mdh_blob(:,1:n_acq); - filePos = reshape( filePos(1:n_acq+1), 1, [] ); % row vector +% discard overallocation: +mdh_blob = mdh_blob(:,1:n_acq); +filePos = reshape( filePos(1:n_acq+1), 1, [] ); % row vector - %fprintf('%8.1f MB read in %4.0f s\n', measLength/1024^2, round(toc(t0))); +%fprintf('%8.1f MB read in %4.0f s\n', measLength/1024^2, round(toc(t0))); % MM (200610) end % of loop_mdh_read() -function [mdh,mask] = evalMDH( mdh_blob, version ) +function [mdh,mask] = evalMDH( mdh_blob, version, isInplacePATref) % see pkg/MrServers/MrMeasSrv/SeqIF/MDH/mdh.h % and pkg/MrServers/MrMeasSrv/SeqIF/MDH/MdhProxy.h @@ -777,8 +782,8 @@ data_uint32 = reshape( data_uint32, [], Nmeas ).'; data_uint16 = reshape( data_uint16, [], Nmeas ).'; data_single = reshape( data_single, [], Nmeas ).'; - % byte pos. -%mdh.ulDMALength = data_uint32(:,1); % 1 : 4 +% byte pos. +mdh.ulDMALength = data_uint32(:,1); % 1 : 4 mdh.lMeasUID = data_uint32(:,2); % 5 : 8 mdh.ulScanCounter = data_uint32(:,3); % 9 : 12 mdh.ulTimeStamp = data_uint32(:,4); % 13 : 16 @@ -820,14 +825,20 @@ mask.MDH_PATREFANDIMASCAN = min(bitand(evalInfoMask1, 2^23),1); mask.MDH_REFLECT = min(bitand(evalInfoMask1, 2^24),1); mask.MDH_NOISEADJSCAN = min(bitand(evalInfoMask1, 2^25),1); -mask.MDH_VOP = min(bitand(mdh.aulEvalInfoMask(2), 2^(53-32)),1); % was 0 in VD +mask.MDH_VOP = min(bitand(mdh.aulEvalInfoMask(:,2), 2^(53-32)),1); % was 0 in VD mask.MDH_IMASCAN = ones( Nmeas, 1, 'uint32' ); noImaScan = ( mask.MDH_ACQEND | mask.MDH_RTFEEDBACK | mask.MDH_HPFEEDBACK ... - | mask.MDH_PHASCOR | mask.MDH_NOISEADJSCAN | mask.MDH_PHASESTABSCAN ... - | mask.MDH_REFPHASESTABSCAN | mask.MDH_SYNCDATA ... - | (mask.MDH_PATREFSCAN & ~mask.MDH_PATREFANDIMASCAN) ); + | mask.MDH_PHASCOR | mask.MDH_NOISEADJSCAN | mask.MDH_PHASESTABSCAN ... + | mask.MDH_REFPHASESTABSCAN | mask.MDH_SYNCDATA ); + +if ~isInplacePATref + noImaScan = (noImaScan | (mask.MDH_PATREFSCAN & ~mask.MDH_PATREFANDIMASCAN)); +end mask.MDH_IMASCAN( noImaScan ) = 0; end % of evalMDH() + + + diff --git a/precommit_unixtime b/precommit_unixtime index caaea4a..3ba6bfa 100644 --- a/precommit_unixtime +++ b/precommit_unixtime @@ -1 +1 @@ -1496740353 +1660732089 diff --git a/read_twix_hdr_Gannet.m b/read_twix_hdr_Gannet.m index f3c925a..1112f8e 100644 --- a/read_twix_hdr_Gannet.m +++ b/read_twix_hdr_Gannet.m @@ -1,9 +1,10 @@ function [prot,rstraj] = read_twix_hdr_Gannet(fid) -% Function to read raw data header information from siemens MRI scanners +% function to read raw data header information from siemens MRI scanners % (currently VB and VD software versions are supported and tested). % -% Author: Philipp Ehses (philipp.ehses@tuebingen.mpg.de), Mar/11/2014 +% Author: Philipp Ehses MPI Tuebingen, Mar/11/2014 +% email: philipp.ehses@dzne.de % Bug fix for data containing long strings: Alex Craven, Oct 2020 nbuffers = fread(fid, 1, 'uint32'); @@ -14,15 +15,15 @@ % ARC 20201020: previous version failed to find null termination on buffer names > 10 characters bl = 64; bufname = fread(fid, bl, 'uint8=>char').'; - + null_ix = find(bufname == 0,1); % ARC: index of null termination - + if isempty(null_ix) % This is not a simple buffer name as would be expected. % Since we don't have a good strategy for dealing with this, just % carry on reading until a null temrinator is found, then discard this - % buffer. - warning('read_twix_hdr: skipping extraordinarily long buffer.'); + % buffer. + warning('read_twix_hdr: skipping extraordinarily long buffer.'); % Keep reading until we find a null while isempty(null_ix) junk = fread(fid, bl, 'uint8=>char').'; @@ -30,14 +31,14 @@ error('read_twix_hdr: reached EOF before null. This is perplexing.'); end null_ix = find(junk == 0,1); - end - fseek(fid, null_ix-bl, 'cof'); + end + fseek(fid, null_ix-bl, 'cof'); continue % skip to the next buffer end - + bufname = bufname(1:null_ix-1); fseek(fid, numel(bufname) - (bl-1), 'cof'); - + buflen = fread(fid, 1,'uint32'); buffer = fread(fid, buflen, 'uint8=>char').'; buffer = regexprep(buffer,'\n\s*\n',''); % delete empty lines @@ -47,14 +48,15 @@ if nargout>1 rstraj = []; if isfield(prot,'Meas') && isfield(prot.Meas,'alRegridMode') && prot.Meas.alRegridMode(1) > 1 - ncol = prot.Meas.alRegridDestSamples(1); - dwelltime = prot.Meas.aflRegridADCDuration(1)/ncol; - gr_adc = zeros(1,ncol,'single'); - start = prot.Meas.alRegridRampupTime(1) - (prot.Meas.aflRegridADCDuration(1)-prot.Meas.alRegridFlattopTime(1))/2; - time_adc = start + dwelltime * (0.5:ncol); - ixUp = time_adc <= prot.Meas.alRegridRampupTime(1); - ixFlat = (time_adc <= prot.Meas.alRegridRampupTime(1)+prot.Meas.alRegridFlattopTime(1)) & ~ixUp; - ixDn = ~ixUp & ~ixFlat; + ncol = prot.Meas.alRegridDestSamples(1); + dwelltime = prot.Meas.aflRegridADCDuration(1)/ncol; + gr_adc = zeros(1,ncol,'single'); + % start = prot.Meas.alRegridRampupTime(1) - (prot.Meas.aflRegridADCDuration(1)-prot.Meas.alRegridFlattopTime(1))/2; + start = prot.Meas.alRegridDelaySamplesTime(1); + time_adc = start + dwelltime * (0.5:ncol); + ixUp = time_adc <= prot.Meas.alRegridRampupTime(1); + ixFlat = (time_adc <= prot.Meas.alRegridRampupTime(1)+prot.Meas.alRegridFlattopTime(1)) & ~ixUp; + ixDn = ~ixUp & ~ixFlat; gr_adc(ixFlat) = 1; if prot.Meas.alRegridMode(1) == 2 % trapezoidal gradient @@ -70,9 +72,9 @@ end % make sure that gr_adc is always positive (rstraj needs to be % strictly monotonic): - gr_adc = max(gr_adc,1e-4); - rstraj = (cumsum(gr_adc(:)) - ncol/2)/sum(gr_adc(:)); - rstraj = rstraj - rstraj(ncol/2+1); + gr_adc = max(gr_adc, 1e-4); + rstraj = (cumtrapz(gr_adc(:)) - ncol/2)/sum(gr_adc(:)); + rstraj = rstraj - mean(rstraj(ncol/2:ncol/2+1)); % scale rstraj by kmax (only works if all slices have same FoV!!!) kmax = prot.MeasYaps.sKSpace.lBaseResolution/... prot.MeasYaps.sSliceArray.asSlice{1}.dReadoutFOV; @@ -116,14 +118,14 @@ if (~isletter(name(1))) name = strcat('x', name); end - + value = char(strtrim(regexprep(tokens{m}(end), '("*)|( *<\w*> *[^\n]*)', ''))); value = regexprep(value, '\s*', ' '); - + try %#ok value = eval(['[' value ']']); % inlined str2num() end - + xprot.(name) = value; end end @@ -137,19 +139,19 @@ isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0; for var=vararray - + try value = eval(['[' var.value ']']); % inlined str2num() catch value = var.value; end - + % now split array name and index (if present) v = regexp(var.name,'(?\w*)\[(?[0-9]*)\]|(?\w*)','names'); - + cnt = 0; tmp = cell(2, numel(v)); - + breaked = false; for k=1:numel(v) if isOctave @@ -169,7 +171,7 @@ cnt = cnt+1; tmp{1,cnt} = '.'; tmp{2,cnt} = vk.name; - + if ~isempty(vk.ix) cnt = cnt+1; tmp{1,cnt} = '{}'; diff --git a/twix_map_obj_Gannet.m b/twix_map_obj_Gannet.m index 7195087..2196666 100644 --- a/twix_map_obj_Gannet.m +++ b/twix_map_obj_Gannet.m @@ -1,1131 +1,1178 @@ classdef twix_map_obj_Gannet < handle - -% class to hold information about raw data from siemens MRI scanners -% (currently VB and VD software versions are supported and tested). -% -% Author: Philipp Ehses (philipp.ehses@tuebingen.mpg.de), Aug/19/2011 -% -% -% Modified by Wolf Blecher (wolf.blecher@tuebingen.mpg.de), Apr/26/2012 -% Added reorder index to indicate which lines are reflected -% Added slice position for sorting, Mai/15/2012 -% -% Order of many mdh parameters are now stored (including the reflected ADC -% bit); PE, Jun/14/2012 -% -% data is now 'memory mapped' and not read until demanded; -% (see mapVBVD for a description) PE, Aug/02/2012 -% -% twix_obj.image.unsorted now returns the data in its acq. order -% [NCol,NCha,nsamples in acq. order], all average flags don't have an -% influence on the output, but 'flagRemoveOS' still works, PE, Sep/04/13 - - -properties(Dependent=true) - readerVersion - % flags - flagRemoveOS % removes oversampling in read (col) during read operation - flagDoAverage % averages over all avg during read operation - flagAverageReps % averages over all repetitions - flagAverageSets % averages over all sets - flagIgnoreSeg % sum over all segments during read operation - flagSkipToFirstLine % skips lines/partitions up to the first - % actually acquired line/partition - % (e.g. only the center k-space is acquired in - % refscans, we don't want all the leading zeros - % in our data) - % this is the default behaviour for everything - % but image scans (but can be changed manually) - flagRampSampRegrid % perform on-the-fly ramp sampling regridding - flagDoRawDataCorrect %SRY: apply raw data correction factors during read operation - - RawDataCorrectionFactors %SRY: allow the user to set/get the factors -end - -properties(GetAccess='public', SetAccess='public') - flagAverageDim % new: flags that determines whether certain dim. should be averaged/ignored - filename - softwareVersion - dataType - rampSampTrj -end - -properties(Dependent=true) - dataSize % this is the current output size, depends on fullSize + some flags - sqzSize - sqzDims -end - -properties(GetAccess='public', SetAccess='protected') - dataDims - - NCol % mdh information - NCha % mdh information - NLin % mdh information - NPar % mdh information - NSli % mdh information - NAve % mdh information - NPhs % mdh information - NEco % mdh information - NRep % mdh information - NSet % mdh information - NSeg % mdh information - NIda % mdh information - NIdb % mdh information - NIdc % mdh information - NIdd % mdh information - NIde % mdh information - NAcq % simple counter - - % mdh information - Lin - Par - Sli - Ave - Phs - Eco - Rep - Set - Seg - Ida - Idb - Idc - Idd - Ide - - centerCol - centerLin - centerPar - cutOff - coilSelect - ROoffcenter - timeSinceRF - IsReflected - IsRawDataCorrect %SRY: storage for MDH flag raw data correct - - slicePos - freeParam - iceParam - scancounter - timestamp - pmutime - - % memory position in file - memPos - - isBrokenFile % errors when parsing? -end - -properties(Hidden=true, SetAccess='protected') - arg % arguments - - fullSize % this is the full size of the data set according to the mdhs, i.e. flags - % like 'reduceOS' have no influence on it - - freadInfo - - skipLin - skipPar -end - -methods - % Constructor: - function this = twix_map_obj_Gannet(arg,dataType,fname,version,rstraj) - - if ~exist('dataType','var') - this.dataType = 'image'; - else - this.dataType = lower(dataType); - end - this.filename = fname; - this.softwareVersion = version; - - this.IsReflected = logical([]); - this.IsRawDataCorrect = logical([]); %SRY - this.NAcq = 0; - this.isBrokenFile = false; - - this.dataDims = {'Col','Cha','Lin','Par','Sli','Ave','Phs',... - 'Eco','Rep','Set','Seg','Ida','Idb','Idc','Idd','Ide'}; - - this.setDefaultFlags(); - if exist('arg','var') - % copy relevant arguments from mapVBVD argument list - names=fieldnames(arg); - for k=1:numel(names) - if isfield(this.arg,names{k}) - this.arg.(names{k}) = arg.(names{k}); - end - end - end + % class to hold information about raw data from siemens MRI scanners + % (currently VB and VD software versions are supported and tested). + % + % Author: Philipp Ehses, MPI Tuebingen, Aug/19/2011 + % email: philipp.ehses@dzne.de + % + % + % Modified by Wolf Blecher (wolf.blecher@tuebingen.mpg.de), Apr/26/2012 + % Added reorder index to indicate which lines are reflected + % Added slice position for sorting, Mai/15/2012 + % + % Order of many mdh parameters are now stored (including the reflected ADC + % bit); PE, Jun/14/2012 + % + % data is now 'memory mapped' and not read until demanded; + % (see mapVBVD for a description) PE, Aug/02/2012 + % + % twix_obj.image.unsorted now returns the data in its acq. order + % [NCol,NCha,nsamples in acq. order], all average flags don't have an + % influence on the output, but 'flagRemoveOS' still works, PE, Sep/04/13 + + + properties(Dependent=true) + readerVersion + % flags + flagRemoveOS % removes oversampling in read (col) during read operation + flagDoAverage % averages over all avg during read operation + flagAverageReps % averages over all repetitions + flagAverageSets % averages over all sets + flagIgnoreSeg % sum over all segments during read operation + flagSkipToFirstLine % skips lines/partitions up to the first + % actually acquired line/partition + % (e.g. only the center k-space is acquired in + % refscans, we don't want all the leading zeros + % in our data) + % this is the default behaviour for everything + % but image scans (but can be changed manually) + flagRampSampRegrid % perform on-the-fly ramp sampling regridding + flagDoRawDataCorrect %SRY: apply raw data correction factors during read operation + + RawDataCorrectionFactors %SRY: allow the user to set/get the factors + end - this.flagAverageDim(ismember(this.dataDims,'Ave')) = this.arg.doAverage; - this.flagAverageDim(ismember(this.dataDims,'Rep')) = this.arg.averageReps; - this.flagAverageDim(ismember(this.dataDims,'Set')) = this.arg.averageSets; - this.flagAverageDim(ismember(this.dataDims,'Seg')) = this.arg.ignoreSeg; - - switch this.softwareVersion - case 'vb' - % every channel has its own full mdh - this.freadInfo.szScanHeader = 0; % [bytes] - this.freadInfo.szChannelHeader = 128; % [bytes] - this.freadInfo.iceParamSz = 4; - case 'vd' - if ( this.arg.doRawDataCorrect ) - error('raw data correction for VD not supported/tested yet'); - end - this.freadInfo.szScanHeader = 192; % [bytes] - this.freadInfo.szChannelHeader = 32; % [bytes] - this.freadInfo.iceParamSz = 24; % vd version supports up to 24 ice params - otherwise - error('software version not supported'); - end + properties(GetAccess='public', SetAccess='public') + flagAverageDim % new: flags that determines whether certain dim. should be averaged/ignored + filename + softwareVersion + dataType + rampSampTrj + end - if exist('rstraj','var') - this.rampSampTrj = rstraj; - else - this.rampSampTrj = []; - this.arg.rampSampRegrid = false; - end + properties(Dependent=true) + dataSize % this is the current output size, depends on fullSize + some flags + sqzSize + sqzDims end - - % Copy function - replacement for matlab.mixin.Copyable.copy() to create object copies - % from http://undocumentedmatlab.com/blog/general-use-object-copy - function newObj = copy(this) - isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0; - if isOctave || verLessThan('matlab','7.11') - % R2010a or earlier - serialize via temp file (slower) - fname = [tempname '.mat']; - save(fname, 'obj'); - newObj = load(fname); - newObj = newObj.obj; - delete(fname); - else - % R2010b or newer - directly in memory (faster) - objByteArray = getByteStreamFromArray(this); - newObj = getArrayFromByteStream(objByteArray); - end + properties(GetAccess='public', SetAccess='protected') + dataDims + + NCol % mdh information + NCha % mdh information + NLin % mdh information + NPar % mdh information + NSli % mdh information + NAve % mdh information + NPhs % mdh information + NEco % mdh information + NRep % mdh information + NSet % mdh information + NSeg % mdh information + NIda % mdh information + NIdb % mdh information + NIdc % mdh information + NIdd % mdh information + NIde % mdh information + NAcq % simple counter + + % mdh information + Lin + Par + Sli + Ave + Phs + Eco + Rep + Set + Seg + Ida + Idb + Idc + Idd + Ide + + centerCol + centerLin + centerPar + cutOff + coilSelect + ROoffcenter + timeSinceRF + IsReflected + IsRawDataCorrect %SRY: storage for MDH flag raw data correct + + slicePos + freeParam + iceParam + evalInfoMask + scancounter + timestamp + pmutime + + % memory position in file + memPos + + isBrokenFile % errors when parsing? end - - function this = readMDH(this, mdh, filePos ) - % extract all values in all MDHs at once - % - % data types: - % Use double for everything non-logical, both ints and floats. Seems the - % most robust way to avoid unexpected cast-issues with very nasty side effects. - % Examples: eps(single(16777216)) == 2 - % uint32( 10 ) - uint32( 20 ) == 0 - % uint16(100) + 1e5 == 65535 - % size(1 : 10000 * uint16(1000)) == [1 65535] - % - % The 1st example always hits the timestamps. - - if ~isstruct( mdh ) || isempty( mdh ) - return - end + properties(Hidden=true, SetAccess='protected') + arg % arguments - this.NAcq = numel( filePos ); - sLC = double( mdh.sLC ) + 1; % +1: convert to matlab index style - evalInfoMask1 = double( mdh.aulEvalInfoMask(:,1) ).'; - - % save mdh information for each line - this.NCol = double( mdh.ushSamplesInScan ).'; - this.NCha = double( mdh.ushUsedChannels ).'; - this.Lin = sLC(:,1).' ; - this.Ave = sLC(:,2).' ; - this.Sli = sLC(:,3).' ; - this.Par = sLC(:,4).' ; - this.Eco = sLC(:,5).' ; - this.Phs = sLC(:,6).' ; - this.Rep = sLC(:,7).' ; - this.Set = sLC(:,8).' ; - this.Seg = sLC(:,9).' ; - this.Ida = sLC(:,10).'; - this.Idb = sLC(:,11).'; - this.Idc = sLC(:,12).'; - this.Idd = sLC(:,13).'; - this.Ide = sLC(:,14).'; - - this.centerCol = double( mdh.ushKSpaceCentreColumn ).' + 1; - this.centerLin = double( mdh.ushKSpaceCentreLineNo ).' + 1; - this.centerPar = double( mdh.ushKSpaceCentrePartitionNo ).' + 1; - this.cutOff = double( mdh.sCutOff ).'; - this.coilSelect = double( mdh.ushCoilSelect ).'; - this.ROoffcenter = double( mdh.fReadOutOffcentre ).'; - this.timeSinceRF = double( mdh.ulTimeSinceLastRF ).'; - this.IsReflected = logical(min(bitand(evalInfoMask1,2^24),1)); - this.scancounter = double( mdh.ulScanCounter ).'; - this.timestamp = double( mdh.ulTimeStamp ).'; - this.pmutime = double( mdh.ulPMUTimeStamp ).'; - this.IsRawDataCorrect = logical(min(bitand(evalInfoMask1,2^10),1)); %SRY - this.slicePos = double( mdh.SlicePos ).'; - this.iceParam = double( mdh.aushIceProgramPara ).'; - this.freeParam = double( mdh.aushFreePara ).'; - - this.memPos = filePos; - - end % of readMDH - - function this = tryAndFixLastMdh(this) - eofWarning = [mfilename() ':UnxpctdEOF']; % We have it inside this.readData() - warning( 'off', eofWarning ) % silence warnings for read... - warning( 'off', 'foo:bar' ) % ... and a stupid placeholder - - isLastAcqGood = false; - cnt = 0; - - while ~isLastAcqGood && this.NAcq > 0 && cnt < 100 - warning( 'foo:bar', 'baz') % make sure that lastwarn() does not return eofWarning - try - this.clean(); - this.unsorted(this.NAcq); - [~, warnid] = lastwarn(); - if strcmp( warnid, eofWarning ) - error( 'Make sure to go to the catch block.') - end - isLastAcqGood = true; - catch - this.isBrokenFile = true; - this.NAcq = this.NAcq-1; - end - cnt = cnt + 1; - end - -% if this.NAcq == 0 || cnt > 99 % everything is garbage -% warning( ) -% end + fullSize % this is the full size of the data set according to the mdhs, i.e. flags + % like 'reduceOS' have no influence on it - warning( 'on', eofWarning ) + freadInfo end - function this = clean(this) - - if this.NAcq == 0 - return; - end + properties(SetAccess='protected') + skipLin + skipPar + end + + methods + % Constructor: + function this = twix_map_obj_Gannet(arg,dataType,fname,version,rstraj) - % Cut mdh data to actual size. Maybe we rejected acquisitions at the end - % due to read errors. - fields = { 'NCol', 'NCha', ... - 'Lin', 'Par', 'Sli', 'Ave', 'Phs', 'Eco', 'Rep', ... - 'Set', 'Seg', 'Ida', 'Idb', 'Idc', 'Idd', 'Ide', ... - 'centerCol' , 'centerLin', 'centerPar', 'cutOff', ... - 'coilSelect' , 'ROoffcenter', 'timeSinceRF', 'IsReflected', ... - 'scancounter', 'timestamp', 'pmutime', 'IsRawDataCorrect', ... - 'slicePos' , 'iceParam', 'freeParam', 'memPos' }; - - nack = this.NAcq; - idx = 1:nack; - - for f = fields - f1 = f{1}; - if size(this.(f1),2) > nack % rarely - this.(f1) = this.(f1)(:,idx); % 1st dim: samples, 2nd dim acquisitions + if ~exist('dataType','var') + this.dataType = 'image'; + else + this.dataType = lower(dataType); end - end - this.NLin = max(this.Lin); - this.NPar = max(this.Par); - this.NSli = max(this.Sli); - this.NAve = max(this.Ave); - this.NPhs = max(this.Phs); - this.NEco = max(this.Eco); - this.NRep = max(this.Rep); - this.NSet = max(this.Set); - this.NSeg = max(this.Seg); - this.NIda = max(this.Ida); - this.NIdb = max(this.Idb); - this.NIdc = max(this.Idc); - this.NIdd = max(this.Idd); - this.NIde = max(this.Ide); - - % ok, let us assume for now that all NCol and NCha entries are - % the same for all mdhs: - this.NCol = this.NCol(1); - this.NCha = this.NCha(1); - - if strcmp(this.dataType,'refscan') - %pehses: check for lines with 'negative' line/partition numbers - %this can happen when the reference scan line/partition range - %exceeds the one of the actual imaging scan - if this.NLin>65500 %uint overflow check - this.Lin = mod(this.Lin + (65536 - min(this.Lin(this.Lin>65500))),65536)+1; - this.NLin = max(this.Lin); + this.filename = fname; + this.softwareVersion = version; + + this.IsReflected = logical([]); + this.IsRawDataCorrect = logical([]); %SRY + this.NAcq = 0; + this.isBrokenFile = false; + + this.dataDims = {'Col','Cha','Lin','Par','Sli','Ave','Phs',... + 'Eco','Rep','Set','Seg','Ida','Idb','Idc','Idd','Ide'}; + + this.setDefaultFlags(); + if exist('arg','var') + % copy relevant arguments from mapVBVD argument list + names=fieldnames(arg); + for k=1:numel(names) + if isfield(this.arg,names{k}) + this.arg.(names{k}) = arg.(names{k}); + end + end end - if this.NPar>65500 %uint overflow check - this.Par = mod(this.Par + (65536 - min(this.Par(this.Par>65500))),65536)+1; - this.NPar = max(this.Par); + + this.flagAverageDim(ismember(this.dataDims,'Ave')) = this.arg.doAverage; + this.flagAverageDim(ismember(this.dataDims,'Rep')) = this.arg.averageReps; + this.flagAverageDim(ismember(this.dataDims,'Set')) = this.arg.averageSets; + this.flagAverageDim(ismember(this.dataDims,'Seg')) = this.arg.ignoreSeg; + + switch this.softwareVersion + case 'vb' + % every channel has its own full mdh + this.freadInfo.szScanHeader = 0; % [bytes] + this.freadInfo.szChannelHeader = 128; % [bytes] + this.freadInfo.iceParamSz = 4; + case 'vd' + if ( this.arg.doRawDataCorrect ) + error('raw data correction for VD not supported/tested yet'); + end + this.freadInfo.szScanHeader = 192; % [bytes] + this.freadInfo.szChannelHeader = 32; % [bytes] + this.freadInfo.iceParamSz = 24; % vd version supports up to 24 ice params + otherwise + error('software version not supported'); end - end - % to reduce the matrix sizes of non-image scans, the size - % of the refscan_obj()-matrix is reduced to the area of the - % actually scanned acs lines (the outer part of k-space - % that is not scanned is not filled with zeros) - % this behaviour is controlled by flagSkipToFirstLine which is - % set to true by default for everything but image scans - if ~this.flagSkipToFirstLine - % the output matrix should include all leading zeros - this.skipLin = 0; - this.skipPar = 0; - else - % otherwise, cut the matrix size to the start of the - % first actually scanned line/partition (e.g. the acs/ - % phasecor data is only acquired in the k-space center) - this.skipLin = min(this.Lin)-1; - this.skipPar = min(this.Par)-1; - end - NLinAlloc = max(1, this.NLin - this.skipLin); - NParAlloc = max(1, this.NPar - this.skipPar); - - this.fullSize = [ this.NCol this.NCha NLinAlloc NParAlloc... - this.NSli this.NAve this.NPhs this.NEco... - this.NRep this.NSet this.NSeg this.NIda... - this.NIdb this.NIdc this.NIdd this.NIde ]; - - nByte = this.NCha*(this.freadInfo.szChannelHeader+8*this.NCol); - - % size for fread - this.freadInfo.sz = [2 nByte/8]; - % reshape size - this.freadInfo.shape = [this.NCol+this.freadInfo.szChannelHeader/8 ... - , this.NCha]; - % we need to cut MDHs from fread data - this.freadInfo.cut = this.freadInfo.szChannelHeader/8 + (1 : this.NCol); - - end % of clean - - - function varargout = subsref(this, S) - % this is where the magic happens - % Now seriously. Overloading of the subsref-method and working - % with a gazillion indices got really messy really fast. At - % some point, I should probably clean this code up a bit. But - % good news everyone: It seems to work. - switch S(1).type - case '.' - % We don't want to manage method/variable calls, so we'll - % simply call the built-in subsref-function in this case. - if nargout == 0 - varargout{1} = builtin('subsref', this, S); % CTR fix. - else - varargout = cell(1, nargout); - [varargout{:}] = builtin('subsref', this, S); - end - return; - case {'()','{}'} - otherwise - error('operator not supported'); + if exist('rstraj','var') + this.rampSampTrj = rstraj; + else + this.rampSampTrj = []; + this.arg.rampSampRegrid = false; + end end - [selRange,selRangeSz,outSize] = this.calcRange(S(1)); - - % calculate page table (virtual to physical addresses) - % this is now done every time, i.e. result is no longer saved in - % a property - slower but safer (and easier to keep track of updates) - ixToRaw = this.calcIndices; - - tmp = reshape(1:prod(double(this.fullSize(3:end))), this.fullSize(3:end)); - tmp = tmp(selRange{3:end}); - ixToRaw = ixToRaw(tmp); clear tmp; - ixToRaw = ixToRaw(:); - % delete all entries that point to zero (the "NULL"-pointer) - notAcquired = (ixToRaw == 0); - ixToRaw (notAcquired) = []; clear notAcquired; - - % calculate ixToTarg for possibly smaller, shifted + segmented - % target matrix: - cIx = ones(14, numel(ixToRaw), 'single'); - if ~this.flagAverageDim(3) - cIx( 1,:) = this.Lin(ixToRaw) - this.skipLin; - end - if ~this.flagAverageDim(4) - cIx( 2,:) = this.Par(ixToRaw) - this.skipPar; - end - if ~this.flagAverageDim(5) - cIx( 3,:) = this.Sli(ixToRaw); - end - if ~this.flagAverageDim(6) - cIx( 4,:) = this.Ave(ixToRaw); - end - if ~this.flagAverageDim(7) - cIx( 5,:) = this.Phs(ixToRaw); - end - if ~this.flagAverageDim(8) - cIx( 6,:) = this.Eco(ixToRaw); - end - if ~this.flagAverageDim(9) - cIx( 7,:) = this.Rep(ixToRaw); - end - if ~this.flagAverageDim(10) - cIx( 8,:) = this.Set(ixToRaw); - end - if ~this.flagAverageDim(11) - cIx( 9,:) = this.Seg(ixToRaw); - end - if ~this.flagAverageDim(12) - cIx(10,:) = this.Ida(ixToRaw); - end - if ~this.flagAverageDim(13) - cIx(11,:) = this.Idb(ixToRaw); - end - if ~this.flagAverageDim(14) - cIx(12,:) = this.Idc(ixToRaw); - end - if ~this.flagAverageDim(15) - cIx(13,:) = this.Idd(ixToRaw); - end - if ~this.flagAverageDim(16) - cIx(14,:) = this.Ide(ixToRaw); - end - % make sure that indices fit inside selection range - for k=3:numel(selRange) - tmp = cIx(k-2,:); - for l=1:numel(selRange{k}) - cIx(k-2,tmp==selRange{k}(l)) = l; + % Copy function - replacement for matlab.mixin.Copyable.copy() to create object copies + % from http://undocumentedmatlab.com/blog/general-use-object-copy + function newObj = copy(this) + isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0; + if isOctave || verLessThan('matlab','7.11') + % R2010a or earlier - serialize via temp file (slower) + fname = [tempname '.mat']; + save(fname, 'obj'); + newObj = load(fname); + newObj = newObj.obj; + delete(fname); + else + % R2010b or newer - directly in memory (faster) + objByteArray = getByteStreamFromArray(this); + newObj = getArrayFromByteStream(objByteArray); end end - sz = selRangeSz(3:end); % extra variable needed for octave compatibility - ixToTarg = this.sub2ind_double(sz, cIx(1,:),cIx(2,:),cIx(3,:),... - cIx(4,:),cIx(5,:),cIx(6,:),cIx(7,:),cIx(8,:),cIx(9,:),... - cIx(10,:),cIx(11,:),cIx(12,:),cIx(13,:),cIx(14,:)); - - mem = this.memPos(ixToRaw); - % sort mem for quicker access, sort cIxToTarg/Raw accordingly - [mem,ix] = sort(mem); - ixToTarg = ixToTarg(ix); - ixToRaw = ixToRaw(ix); - clear ix; - - % For a call of type data{:,:,1:3} matlab expects more than one - % output variable (three in this case) and will throw an error - % otherwise. This is a lazy way (and the only one I know of) to - % fix this. - varargout = cell(1, nargout); - varargout{1} = this.readData(mem,ixToTarg,ixToRaw,selRange,selRangeSz,outSize); - end % of subsref - - - function out = unsorted(this,ival) - % returns the unsorted data [NCol,NCha,#samples in acq. order] - if ~exist('ival','var') - mem = this.memPos; - else - mem = this.memPos(ival); - end - out = this.readData(mem); - end - - function out = readData(this,mem,cIxToTarg,cIxToRaw,selRange,selRangeSz,outSize) - - if ~exist('outSize','var') - selRange{1} = ':'; - selRange{2} = ':'; - outSize = [this.dataSize(1:2),numel(mem)]; - selRangeSz = outSize; - cIxToTarg = 1:selRangeSz(3); - cIxToRaw = cIxToTarg; - else - if isequal( selRange{1}(:), (1:this.dataSize(1)).' ) - selRange{1} = ':'; + function this = readMDH(this, mdh, filePos ) + % extract all values in all MDHs at once + % + % data types: + % Use double for everything non-logical, both ints and floats. Seems the + % most robust way to avoid unexpected cast-issues with very nasty side effects. + % Examples: eps(single(16777216)) == 2 + % uint32( 10 ) - uint32( 20 ) == 0 + % uint16(100) + 1e5 == 65535 + % size(1 : 10000 * uint16(1000)) == [1 65535] + % + % The 1st example always hits the timestamps. + + if ~isstruct( mdh ) || isempty( mdh ) + return end - if isequal( selRange{2}(:), (1:this.dataSize(2)).' ) - selRange{2} = ':'; + + this.NAcq = numel( filePos ); + sLC = double( mdh.sLC ) + 1; % +1: convert to matlab index style + + % save mdh information for each line + this.NCol = double( mdh.ushSamplesInScan ).'; + this.NCha = double( mdh.ushUsedChannels ).'; + this.Lin = sLC(:,1).' ; + this.Ave = sLC(:,2).' ; + this.Sli = sLC(:,3).' ; + this.Par = sLC(:,4).' ; + this.Eco = sLC(:,5).' ; + this.Phs = sLC(:,6).' ; + this.Rep = sLC(:,7).' ; + this.Set = sLC(:,8).' ; + this.Seg = sLC(:,9).' ; + this.Ida = sLC(:,10).'; + this.Idb = sLC(:,11).'; + this.Idc = sLC(:,12).'; + this.Idd = sLC(:,13).'; + this.Ide = sLC(:,14).'; + + this.evalInfoMask = double(mdh.aulEvalInfoMask).'; + evalInfoMask1 = this.evalInfoMask(1,:); + + this.centerCol = double( mdh.ushKSpaceCentreColumn ).' + 1; + this.centerLin = double( mdh.ushKSpaceCentreLineNo ).' + 1; + this.centerPar = double( mdh.ushKSpaceCentrePartitionNo ).' + 1; + this.cutOff = double( mdh.sCutOff ).'; + this.coilSelect = double( mdh.ushCoilSelect ).'; + this.ROoffcenter = double( mdh.fReadOutOffcentre ).'; + this.timeSinceRF = double( mdh.ulTimeSinceLastRF ).'; + this.IsReflected = logical(min(bitand(evalInfoMask1,2^24),1)); + this.scancounter = double( mdh.ulScanCounter ).'; + this.timestamp = double( mdh.ulTimeStamp ).'; + this.pmutime = double( mdh.ulPMUTimeStamp ).'; + this.IsRawDataCorrect = logical(min(bitand(evalInfoMask1,2^10),1)); %SRY + this.slicePos = double( mdh.SlicePos ).'; + this.iceParam = double( mdh.aushIceProgramPara ).'; + this.freeParam = double( mdh.aushFreePara ).'; + this.memPos = filePos; + + end % of readMDH + + function this = tryAndFixLastMdh(this) + eofWarning = [mfilename() ':UnxpctdEOF']; % We have it inside this.readData() + warning( 'off', eofWarning ) % silence warnings for read... + warning( 'off', 'foo:bar' ) % ... and a stupid placeholder + + isLastAcqGood = false; + cnt = 0; + + while ~isLastAcqGood && this.NAcq > 0 && cnt < 100 + warning( 'foo:bar', 'baz') % make sure that lastwarn() does not return eofWarning + try + this.clean(); + this.unsorted(this.NAcq); + [~, warnid] = lastwarn(); + if strcmp( warnid, eofWarning ) + error( 'Make sure to go to the catch block.') + end + isLastAcqGood = true; + catch + this.isBrokenFile = true; + this.NAcq = this.NAcq-1; + end + cnt = cnt + 1; end - end - out = complex(zeros(outSize,'single')); - out = reshape(out, selRangeSz(1), selRangeSz(2), []); - if isempty( mem ) - out = reshape(out,outSize); - return - end + % if this.NAcq == 0 || cnt > 99 % everything is garbage + % warning( ) + % end - cIxToTarg = this.cast2MinimalUint( cIxToTarg ); - - % subsref overloading makes this.that-calls slow, so we need to - % avoid them whenever possible - szScanHeader = this.freadInfo.szScanHeader; - readSize = this.freadInfo.sz; - readShape = this.freadInfo.shape; - readCut = this.freadInfo.cut; - keepOS = [1:this.NCol/4, 1+this.NCol*3/4:this.NCol]; - bRemoveOS = this.arg.removeOS; - bIsReflected = this.IsReflected(cIxToRaw); - bRegrid = this.flagRampSampRegrid && numel(this.rampSampTrj); - slicedata = this.slicePos(:, cIxToRaw); - %SRY store information about raw data correction - bDoRawDataCorrect = this.arg.doRawDataCorrect; - bIsRawDataCorrect = this.IsRawDataCorrect( cIxToRaw ); - isBrokenRead = false; - if (bDoRawDataCorrect) - rawDataCorrect = this.arg.rawDataCorrectionFactors; + warning( 'on', eofWarning ) end - % MiVö: Raw data are read line-by-line in portions of 2xNColxNCha float32 points (2 for complex). - % Computing and sorting(!) on these small portions is quite expensive, esp. when - % it employs non-sequential memory paths. Examples are non-linear k-space acquisition - % or reflected lines. - % This can be sped up if slightly larger blocks of raw data are collected, first. - % Whenever a block is full, we do all those operations and save it in the final "out" array. - % What's a good block size? Depends on data size and machine (probably L2/L3/L4 cache sizes). - % So...? Start with a small block, measure the time-per-line and double block size until - % a minimum is found. Seems sufficiently robust to end up in a close-to-optimal size for every - % machine and data. - blockSz = 2; % size of blocks; must be 2^n; will be increased - doLockblockSz = false; % whether blockSZ should be left untouched - tprev = inf; % previous time-per-line - blockCtr = 0; - blockInit = -inf(readShape(1), readShape(2), blockSz, 'single'); %init with garbage - blockInit = complex( blockInit ); - block = blockInit; - - if bRegrid - v1 = single(1:selRangeSz(2)); - v2 = single(1:blockSz); - rsTrj = {this.rampSampTrj,v1,v2}; - trgTrj = linspace(min(this.rampSampTrj),max(this.rampSampTrj),this.NCol); - trgTrj = {trgTrj,v1,v2}; - end - - % counter for proper scaling of averages/segments - count_ave = zeros([1 1 size(out,3)],'single'); - kMax = numel( mem ); % max loop index - - fid = this.fileopen(); - - for k = 1:kMax - % skip scan header - fseek(fid,mem(k) + szScanHeader,'bof'); - raw = fread(fid, readSize, 'float=>single').'; - - % MiVö: With incomplete files fread() returns less than readSize points. The subsequent reshape will therefore error out. - % We could check if numel(raw) == prod(readSize), but people recommend exception handling for performance - % reasons. Do it. - try - raw = reshape( complex(raw(:,1), raw(:,2)), readShape); - catch exc - offset_bytes = mem(k) + szScanHeader; - %remainingSz = readSize(2) - size(raw,1); - warning( [mfilename() ':UnxpctdEOF'], ... - [ '\nAn unexpected read error occurred at this byte offset: %d (%g GiB)\n'... - 'Actual read size is [%s], desired size was: [%s]\n' ... - 'Will ignore this line and stop reading.\n' ... - '=== MATLABs error message ================\n' ... - exc.message ... - '\n=== end of error =========================\n' ... - ], offset_bytes, offset_bytes/1024^3, num2str(size(raw)), num2str(readSize.') ) - - % Reject this data fragment. To do so, init with the values of blockInit - clear raw - raw( 1:prod(readShape) ) = blockInit(1); - raw = reshape( raw, readShape ); - isBrokenRead = true; % remember it and bail out later + function this = clean(this) + + if this.NAcq == 0 + return; + end + + % Cut mdh data to actual size. Maybe we rejected acquisitions at the end + % due to read errors. + fields = { 'NCol', 'NCha',... + 'Lin', 'Par', 'Sli', 'Ave', 'Phs', 'Eco', 'Rep',... + 'Set', 'Seg', 'Ida', 'Idb', 'Idc', 'Idd', 'Ide',... + 'centerCol', 'centerLin', 'centerPar', 'cutOff',... + 'coilSelect' , 'ROoffcenter', 'timeSinceRF',... + 'IsReflected', 'scancounter', 'timestamp', 'pmutime',... + 'IsRawDataCorrect', 'slicePos', 'iceParam',... + 'freeParam', 'evalInfoMask', 'memPos' }; + + nack = this.NAcq; + idx = 1:nack; + + for f = fields + f1 = f{1}; + if size(this.(f1),2) > nack % rarely + this.(f1) = this.(f1)(:,idx); % 1st dim: samples, 2nd dim acquisitions + end end - blockCtr = blockCtr + 1; - block(:,:,blockCtr) = raw; % fast serial storage in a cache array - - % Do expensive computations and reorderings on the gathered block. - % Unfortunately, a lot of code is necessary, but that is executed much less - % frequent, so its worthwhile for speed. - % TODO: Do *everything* block-by-block - if blockCtr == blockSz || k == kMax || (isBrokenRead && blockCtr > 1) - s = tic; % measure the time to process a block of data - - % remove MDH data from block: - block = block(readCut,:,:); - - if bRegrid - % correct for readout shifts - % the nco frequency is always scaled to the max. - % gradient amp and does account for ramp-sampling - ro_shift = this.calcOffcenterShiftRO(slicedata(:,k)); - deltak = max(abs(diff(rsTrj{1}))); - phase = (0:this.NCol-1).' * deltak * ro_shift; - phase_factor = exp(1j*2*pi*(phase - ro_shift * rsTrj{1})); - block = bsxfun(@times, block, phase_factor); - - % grid the data - F = griddedInterpolant(rsTrj, block); - block = F(trgTrj); + this.NLin = max(this.Lin); + this.NPar = max(this.Par); + this.NSli = max(this.Sli); + this.NAve = max(this.Ave); + this.NPhs = max(this.Phs); + this.NEco = max(this.Eco); + this.NRep = max(this.Rep); + this.NSet = max(this.Set); + this.NSeg = max(this.Seg); + this.NIda = max(this.Ida); + this.NIdb = max(this.Idb); + this.NIdc = max(this.Idc); + this.NIdd = max(this.Idd); + this.NIde = max(this.Ide); + + % ok, let us assume for now that all NCol and NCha entries are + % the same for all mdhs, or at least make sure to take the max() + this.NCol = max(this.NCol); + this.NCha = max(this.NCha); + + if strcmp(this.dataType,'refscan') + %pehses: check for lines with 'negative' line/partition numbers + %this can happen when the reference scan line/partition range + %exceeds the one of the actual imaging scan + if this.NLin>65500 %uint overflow check + this.Lin = mod(this.Lin + (65536 - min(this.Lin(this.Lin>65500))),65536)+1; + this.NLin = max(this.Lin); end - - ix = 1 + k - blockCtr : k; - if blockCtr ~= blockSz - block = block(:,:,1:blockCtr); + if this.NPar>65500 %uint overflow check + this.Par = mod(this.Par + (65536 - min(this.Par(this.Par>65500))),65536)+1; + this.NPar = max(this.Par); end - - if sum(isnan(block(:)))>0 - keyboard + end + + % to reduce the matrix sizes of non-image scans, the size + % of the refscan_obj()-matrix is reduced to the area of the + % actually scanned acs lines (the outer part of k-space + % that is not scanned is not filled with zeros) + % this behaviour is controlled by flagSkipToFirstLine which is + % set to true by default for everything but image scans + if ~this.flagSkipToFirstLine + % the output matrix should include all leading zeros + this.skipLin = 0; + this.skipPar = 0; + else + % otherwise, cut the matrix size to the start of the + % first actually scanned line/partition (e.g. the acs/ + % phasecor data is only acquired in the k-space center) + this.skipLin = min(this.Lin)-1; + this.skipPar = min(this.Par)-1; + end + NLinAlloc = max(1, this.NLin - this.skipLin); + NParAlloc = max(1, this.NPar - this.skipPar); + + this.fullSize = [ this.NCol this.NCha NLinAlloc NParAlloc... + this.NSli this.NAve this.NPhs this.NEco... + this.NRep this.NSet this.NSeg this.NIda... + this.NIdb this.NIdc this.NIdd this.NIde ]; + + nByte = this.NCha*(this.freadInfo.szChannelHeader+8*this.NCol); + + % size for fread + this.freadInfo.sz = [2 nByte/8]; + % reshape size + this.freadInfo.shape = [this.NCol+this.freadInfo.szChannelHeader/8 ... + , this.NCha]; + % we need to cut MDHs from fread data + this.freadInfo.cut = this.freadInfo.szChannelHeader/8 + (1 : this.NCol); + + end % of clean + + + function varargout = subsref(this, S) + % this is where the magic happens + % Now seriously. Overloading of the subsref-method and working + % with a gazillion indices got really messy really fast. At + % some point, I should probably clean this code up a bit. But + % good news everyone: It seems to work. + switch S(1).type + case '.' + % We don't want to manage method/variable calls, so we'll + % simply call the built-in subsref-function in this case. + if nargout == 0 + varargout{1} = builtin('subsref', this, S); % CTR fix. + else + varargout = cell(1, nargout); + [varargout{:}] = builtin('subsref', this, S); + end + return; + case {'()','{}'} + otherwise + error('operator not supported'); + end + + [selRange,selRangeSz,outSize] = this.calcRange(S(1)); + + % calculate page table (virtual to physical addresses) + % this is now done every time, i.e. result is no longer saved in + % a property - slower but safer (and easier to keep track of updates) + ixToRaw = this.calcIndices; + + tmp = reshape(1:prod(double(this.fullSize(3:end))), this.fullSize(3:end)); + tmp = tmp(selRange{3:end}); + ixToRaw = ixToRaw(tmp); clear tmp; + ixToRaw = ixToRaw(:); + % delete all entries that point to zero (the "NULL"-pointer) + notAcquired = (ixToRaw == 0); + ixToRaw (notAcquired) = []; clear notAcquired; + + % calculate ixToTarg for possibly smaller, shifted + segmented + % target matrix: + cIx = ones(14, numel(ixToRaw), 'single'); + if ~this.flagAverageDim(3) + cIx( 1,:) = this.Lin(ixToRaw) - this.skipLin; + end + if ~this.flagAverageDim(4) + cIx( 2,:) = this.Par(ixToRaw) - this.skipPar; + end + if ~this.flagAverageDim(5) + cIx( 3,:) = this.Sli(ixToRaw); + end + if ~this.flagAverageDim(6) + cIx( 4,:) = this.Ave(ixToRaw); + end + if ~this.flagAverageDim(7) + cIx( 5,:) = this.Phs(ixToRaw); + end + if ~this.flagAverageDim(8) + cIx( 6,:) = this.Eco(ixToRaw); + end + if ~this.flagAverageDim(9) + cIx( 7,:) = this.Rep(ixToRaw); + end + if ~this.flagAverageDim(10) + cIx( 8,:) = this.Set(ixToRaw); + end + if ~this.flagAverageDim(11) + cIx( 9,:) = this.Seg(ixToRaw); + end + if ~this.flagAverageDim(12) + cIx(10,:) = this.Ida(ixToRaw); + end + if ~this.flagAverageDim(13) + cIx(11,:) = this.Idb(ixToRaw); + end + if ~this.flagAverageDim(14) + cIx(12,:) = this.Idc(ixToRaw); + end + if ~this.flagAverageDim(15) + cIx(13,:) = this.Idd(ixToRaw); + end + if ~this.flagAverageDim(16) + cIx(14,:) = this.Ide(ixToRaw); + end + + % make sure that indices fit inside selection range + for k=3:numel(selRange) + tmp = cIx(k-2,:); + for l=1:numel(selRange{k}) + cIx(k-2,tmp==selRange{k}(l)) = l; end + end + + sz = selRangeSz(3:end); % extra variable needed for octave compatibility + ixToTarg = this.sub2ind_double(sz, cIx(1,:),cIx(2,:),cIx(3,:),... + cIx(4,:),cIx(5,:),cIx(6,:),cIx(7,:),cIx(8,:),cIx(9,:),... + cIx(10,:),cIx(11,:),cIx(12,:),cIx(13,:),cIx(14,:)); + + mem = this.memPos(ixToRaw); + % sort mem for quicker access, sort cIxToTarg/Raw accordingly + [mem,ix] = sort(mem); + ixToTarg = ixToTarg(ix); + ixToRaw = ixToRaw(ix); + clear ix; + + % For a call of type data{:,:,1:3} matlab expects more than one + % output variable (three in this case) and will throw an error + % otherwise. This is a lazy way (and the only one I know of) to + % fix this. + varargout = cell(1, nargout); + varargout{1} = this.readData(mem,ixToTarg,ixToRaw,selRange,selRangeSz,outSize); + + end % of subsref - if bRemoveOS % remove oversampling in read - block = ifft(block, [], 1); - block = fft(block(keepOS,:,:), [], 1); + + function out = unsorted(this,ival) + % returns the unsorted data [NCol,NCha,#samples in acq. order] + if ~exist('ival','var') + mem = this.memPos; + else + mem = this.memPos(ival); + end + out = this.readData(mem); + end + + + function out = readData(this,mem,cIxToTarg,cIxToRaw,selRange,selRangeSz,outSize) + + if ~exist('outSize','var') + selRange{1} = ':'; + selRange{2} = ':'; + outSize = [this.dataSize(1:2),numel(mem)]; + selRangeSz = outSize; + cIxToTarg = 1:selRangeSz(3); + cIxToRaw = cIxToTarg; + else + if isequal( selRange{1}(:), (1:this.dataSize(1)).' ) + selRange{1} = ':'; end - - if ( bDoRawDataCorrect && bIsRawDataCorrect(k) ) - %SRY apply raw data correction if necessary - block = bsxfun(@times, block, rawDataCorrect); + if isequal( selRange{2}(:), (1:this.dataSize(2)).' ) + selRange{2} = ':'; end - - isRefl = bIsReflected(ix); - block(:,:,isRefl) = block(end:-1:1,:,isRefl); - - if ~isequal(selRange{1},':') || ~isequal(selRange{2},':') - block = block( selRange{1}, selRange{2}, : ); % a bit slow + end + + try % much faster + out = zeros(outSize, 'like', complex(single(0))); + catch % legacy support + out = complex(zeros(outSize,'single')); + end + out = reshape(out, selRangeSz(1), selRangeSz(2), []); + + if isempty( mem ) + out = reshape(out,outSize); + return + end + + cIxToTarg = this.cast2MinimalUint( cIxToTarg ); + + % subsref overloading makes this.that-calls slow, so we need to + % avoid them whenever possible + szScanHeader = this.freadInfo.szScanHeader; + readSize = this.freadInfo.sz; + readShape = this.freadInfo.shape; + readCut = this.freadInfo.cut; + keepOS = [1:this.NCol/4, 1+this.NCol*3/4:this.NCol]; + bRemoveOS = this.arg.removeOS; + bIsReflected = this.IsReflected(cIxToRaw); + bRegrid = this.flagRampSampRegrid && numel(this.rampSampTrj); + ro_shift = this.ROoffcenter(cIxToRaw); + %SRY store information about raw data correction + bDoRawDataCorrect = this.arg.doRawDataCorrect; + bIsRawDataCorrect = this.IsRawDataCorrect( cIxToRaw ); + isBrokenRead = false; + if (bDoRawDataCorrect) + rawDataCorrect = this.arg.rawDataCorrectionFactors; + end + + % MiVö: Raw data are read line-by-line in portions of 2xNColxNCha float32 points (2 for complex). + % Computing and sorting(!) on these small portions is quite expensive, esp. when + % it employs non-sequential memory paths. Examples are non-linear k-space acquisition + % or reflected lines. + % This can be sped up if slightly larger blocks of raw data are collected, first. + % Whenever a block is full, we do all those operations and save it in the final "out" array. + % What's a good block size? Depends on data size and machine (probably L2/L3/L4 cache sizes). + % So...? Start with a small block, measure the time-per-line and double block size until + % a minimum is found. Seems sufficiently robust to end up in a close-to-optimal size for every + % machine and data. + blockSz = 2; % size of blocks; must be 2^n; will be increased + doLockblockSz = false; % whether blockSZ should be left untouched + tprev = inf; % previous time-per-line + blockCtr = 0; + blockInit = -inf(readShape(1), readShape(2), blockSz, 'single'); %init with garbage + blockInit = complex( blockInit ); + block = blockInit; + + if bRegrid + v1 = single(1:selRangeSz(2)*blockSz); + rsTrj = {this.rampSampTrj(:),v1(:)}; + trgTrj = linspace(min(this.rampSampTrj),max(this.rampSampTrj),this.NCol); + trgTrj = {trgTrj(:),v1(:)}; + end + + % counter for proper scaling of averages/segments + count_ave = zeros([1 1 size(out,3)],'single'); + kMax = numel( mem ); % max loop index + + fid = this.fileopen(); + + for k = 1:kMax + % skip scan header + fseek(fid,mem(k) + szScanHeader,'bof'); + raw = fread(fid, readSize, 'float=>single').'; + + if any(isnan(raw(:))) + %pehses: temporary fix for VE MPRAGE data + raw(isnan(raw(:))) = 0; end - - [sortIdx, I] = sort( cIxToTarg(ix), 'ascend' ); - block = block(:,:,I); % reorder according to sorted target indices - % Mark duplicate indices with 1; we'll have to treat them special for proper averaging - % Bonus: The very first storage can be made much faster, because it's in-place. - % Matlab urgently needs a "+=" operater, which makes "A(:,:,idx) = A(:,:,idx) + B" - % in-place and more readable. - isDupe = [ false, diff(sortIdx) == 0 ]; + % MiVö: With incomplete files fread() returns less than readSize points. The subsequent reshape will therefore error out. + % We could check if numel(raw) == prod(readSize), but people recommend exception handling for performance + % reasons. Do it. + try + raw = reshape(complex(raw(:,1), raw(:,2)), readShape); + catch exc + offset_bytes = mem(k) + szScanHeader; + %remainingSz = readSize(2) - size(raw,1); + warning( [mfilename() ':UnxpctdEOF'], ... + [ '\nAn unexpected read error occurred at this byte offset: %d (%g GiB)\n'... + 'Actual read size is [%s], desired size was: [%s]\n' ... + 'Will ignore this line and stop reading.\n' ... + '=== MATLABs error message ================\n' ... + exc.message ... + '\n=== end of error =========================\n' ... + ], offset_bytes, offset_bytes/1024^3, num2str(size(raw)), num2str(readSize.') ) + + % Reject this data fragment. To do so, init with the values of blockInit + clear raw + raw( 1:prod(readShape) ) = blockInit(1); + raw = reshape( raw, readShape ); + isBrokenRead = true; % remember it and bail out later + end - idx1 = sortIdx(~isDupe); % acquired once in this block - idxN = sortIdx( isDupe); % acquired multiple times + blockCtr = blockCtr + 1; + block(:,:,blockCtr) = raw; % fast serial storage in a cache array - count_ave(idx1) = count_ave(idx1) + 1; + % Do expensive computations and reorderings on the gathered block. + % Unfortunately, a lot of code is necessary, but that is executed much less + % frequent, so its worthwhile for speed. + % TODO: Do *everything* block-by-block + if blockCtr == blockSz || k == kMax || (isBrokenRead && blockCtr > 1) + s = tic; % measure the time to process a block of data - if isempty( idxN ) - % no duplicates - if all( count_ave(idx1) == 1 ) % first acquisition of this line - out(:,:,idx1) = block; % fast - else - out(:,:,idx1) = out(:,:,idx1) + block; % slow + % remove MDH data from block: + block = block(readCut,:,:); + + ix = 1 + k - blockCtr : k; + if blockCtr ~= blockSz + block = block(:,:,1:blockCtr); + end + + % reflect fids when necessary + isRefl = bIsReflected(ix); + block(:,:,isRefl) = block(end:-1:1,:,isRefl); + + if bRegrid + % correct for readout shifts + % the nco frequency is always scaled to the max. + % gradient amp and does account for ramp-sampling + deltak = max(abs(diff(rsTrj{1}))); + fovshift = reshape(ro_shift(ix),1,1,[]); + adcphase = deltak * bsxfun(@times, (0:this.NCol-1).', fovshift); + fovphase = bsxfun(@times, fovshift, rsTrj{1}); + phase_factor = exp(1j*2*pi*(adcphase - fovphase)); + block = bsxfun(@times, block, phase_factor); + + % grid the data + sz = size(block); + ninterp = prod(sz(2:end)); + rsTrj{2} = single(1:ninterp).'; + trgTrj{2} = rsTrj{2}; + if ninterp == 1 + blockdims = 1; + else + blockdims = [1, 2]; + end + F = griddedInterpolant(rsTrj(blockdims), block(:,:)); + block = reshape(F(trgTrj(blockdims)), sz); end - else - out(:,:,idx1) = out(:,:,idx1) + block(:,:,~isDupe); % slower - block = block(:,:,isDupe); - for n = 1:numel(idxN) - out(:,:,idxN(n)) = out(:,:,idxN(n)) + block(:,:,n); % snail :-) - count_ave(idxN(n)) = count_ave(idxN(n)) + 1; + if bRemoveOS % remove oversampling in read + block = ifft(block, [], 1); + block = fft(block(keepOS,:,:), [], 1); end - end - % At the first few iterations, evaluate the spent time-per-line and decide - % what to do with the block size. - if ~doLockblockSz - t = 1e6 * toc(s)/blockSz; % micro seconds + if ( bDoRawDataCorrect && bIsRawDataCorrect(k) ) + %SRY apply raw data correction if necessary + block = bsxfun(@times, block, rawDataCorrect); + end + + if this.flagAverageDim(1) + block = mean(block, 1); + elseif ~isequal(selRange{1},':') + block = block( selRange{1}, :, : ); % a bit slow + end - if t <= 1.1 * tprev % allow 10% inaccuracy. Usually bigger == better - % New block size was faster. Go a step further. - blockSz = blockSz * 2; - blockInit = cat(3, blockInit, blockInit); + if this.flagAverageDim(2) + block = mean(block, 2); else - % regression; reset size and lock it - blockSz = max( blockSz/2, 1 ); - blockInit = blockInit(:,:,1:blockSz); - doLockblockSz = true; + if ~isequal(selRange{2},':') + block = block( :, selRange{2}, : ); % a bit slow + end end - if bRegrid - rsTrj{3} = single(1:blockSz); - trgTrj{3} = rsTrj{3}; + + [sortIdx, I] = sort( cIxToTarg(ix), 'ascend' ); + block = block(:,:,I); % reorder according to sorted target indices + + % Mark duplicate indices with 1; we'll have to treat them special for proper averaging + % Bonus: The very first storage can be made much faster, because it's in-place. + % Matlab urgently needs a "+=" operater, which makes "A(:,:,idx) = A(:,:,idx) + B" + % in-place and more readable. + isDupe = [ false, diff(sortIdx) == 0 ]; + + idx1 = sortIdx(~isDupe); % acquired once in this block + idxN = sortIdx( isDupe); % acquired multiple times + + count_ave(idx1) = count_ave(idx1) + 1; + + if isempty( idxN ) + % no duplicates + if all( count_ave(idx1) == 1 ) % first acquisition of this line + out(:,:,idx1) = block; % fast + else + out(:,:,idx1) = out(:,:,idx1) + block; % slow + end + else + out(:,:,idx1) = out(:,:,idx1) + block(:,:,~isDupe); % slower + + block = block(:,:,isDupe); + for n = 1:numel(idxN) + out(:,:,idxN(n)) = out(:,:,idxN(n)) + block(:,:,n); % snail :-) + count_ave(idxN(n)) = count_ave(idxN(n)) + 1; + end + end + + % At the first few iterations, evaluate the spent time-per-line and decide + % what to do with the block size. + if ~doLockblockSz + t = 1e6 * toc(s)/blockSz; % micro seconds + + if t <= 1.1 * tprev % allow 10% inaccuracy. Usually bigger == better + % New block size was faster. Go a step further. + blockSz = blockSz * 2; + blockInit = cat(3, blockInit, blockInit); + else + % regression; reset size and lock it + blockSz = max( blockSz/2, 1 ); + blockInit = blockInit(:,:,1:blockSz); + doLockblockSz = true; + end + tprev = t; end - tprev = t; + + blockCtr = 0; + block = blockInit; % reset to garbage + end + + if isBrokenRead + this.isBrokenFile = true; + break end - - blockCtr = 0; - block = blockInit; % reset to garbage end - if isBrokenRead - this.isBrokenFile = true; - break + fclose(fid); + + % proper scaling (we don't want to sum our data but average it) + % For large "out" bsxfun(@rdivide,out,count_ave) is incredibly faster than + % bsxfun(@times,out,count_ave)! + % @rdivide is also running in parallel, while @times is not. :-/ + if any( reshape(count_ave,[],1) > 1 ) + clearvars -except out count_ave outSize + count_ave = max( 1, count_ave ); + out = bsxfun( @rdivide, out, count_ave); + end + + out = reshape(out,outSize); + end % of readData + + + function setDefaultFlags(this) + % method to set flags to default values + this.arg.removeOS = false; + this.arg.rampSampRegrid = false; + this.arg.doAverage = false; + this.arg.averageReps = false; + this.arg.averageSets = false; + this.arg.ignoreSeg = false; + this.arg.doRawDataCorrect = false; + this.flagAverageDim = false(1, 16); + + if strcmp(this.dataType,'image') || strcmp(this.dataType,'phasecor') || strcmp(this.dataType,'phasestab') + this.arg.skipToFirstLine = false; + else + this.arg.skipToFirstLine = true; + end + if ~isfield(this.arg,'rawDataCorrectionFactors') + this.arg.rawDataCorrectionFactors = []; end end - fclose(fid); + function dummy = resetFlags(this) + % method to reset flags to default values + this.flagRemoveOS = false; + this.flagRampSampRegrid = false; + this.flagDoRawDataCorrect = false; + this.flagAverageDim = false(1,16); - % proper scaling (we don't want to sum our data but average it) - % For large "out" bsxfun(@rdivide,out,count_ave) is incredibly faster than - % bsxfun(@times,out,count_ave)! - % @rdivide is also running in parallel, while @times is not. :-/ - if any( reshape(count_ave,[],1) > 1 ) - clearvars -except out count_ave outSize - count_ave = max( 1, count_ave ); - out = bsxfun( @rdivide, out, count_ave); + if strcmp(this.dataType,'image') || strcmp(this.dataType,'phasecor') || strcmp(this.dataType,'phasestab') + this.arg.skipToFirstLine = false; + else + this.arg.skipToFirstLine = true; + end + dummy = []; end - out = reshape(out,outSize); - end % of readData - - - function setDefaultFlags(this) - % method to set flags to default values - this.arg.removeOS = false; - this.arg.rampSampRegrid = false; - this.arg.doAverage = false; - this.arg.averageReps = false; - this.arg.averageSets = false; - this.arg.ignoreSeg = false; - this.arg.doRawDataCorrect = false; - this.flagAverageDim = false(1, 16); - - if strcmp(this.dataType,'image') || strcmp(this.dataType,'phasecor') || strcmp(this.dataType,'phasestab') - this.arg.skipToFirstLine = false; - else - this.arg.skipToFirstLine = true; + + function versiontime = get.readerVersion(~) + % returns utc-unixtime of last commit (from file precommit-unixtime) + p = fileparts(mfilename('fullpath')); + fid = fopen(fullfile(p, 'precommit_unixtime')); + if fid ~= -1 + versiontime = uint64(str2double(fgetl(fid))); + fclose(fid); + else + versiontime = 0; + end end - if ~isfield(this.arg,'rawDataCorrectionFactors') - this.arg.rawDataCorrectionFactors = []; + + function out = get.dataSize(this) + out = this.fullSize; + + if this.arg.removeOS + ix = ismember(this.dataDims, 'Col'); + out(ix) = this.NCol/2; + end + + out(this.flagAverageDim) = 1; end - end - function dummy = resetFlags(this) - % method to reset flags to default values - this.flagRemoveOS = false; - this.flagRampSampRegrid = false; - this.flagDoRawDataCorrect = false; - this.flagAverageDim = false(1,16); - - if strcmp(this.dataType,'image') || strcmp(this.dataType,'phasecor') || strcmp(this.dataType,'phasestab') - this.arg.skipToFirstLine = false; - else - this.arg.skipToFirstLine = true; + + function out = get.sqzDims(this) + out = this.dataDims(this.dataSize>1); end - dummy = []; - end - - - function versiontime = get.readerVersion(~) - % returns utc-unixtime of last commit (from file precommit-unixtime) - p = fileparts(mfilename('fullpath')); - fid = fopen(fullfile(p, 'precommit_unixtime')); - versiontime = uint64(str2double(fgetl(fid))); - fclose(fid); - end - - function out = get.dataSize(this) - out = this.fullSize; - - if this.arg.removeOS - ix = ismember(this.dataDims, 'Col'); - out(ix) = this.NCol/2; + + + function out = get.sqzSize(this) + out = this.dataSize(this.dataSize>1); end - - if this.flagAverageDim(1) || this.flagAverageDim(2) - warning('averaging in col and cha dim not supported, resetting flag'); - this.flagAverageDim(1:2) = false; + + + function set.flagRemoveOS(this,val) + % set method for removeOS + this.arg.removeOS = logical(val); end - - out(this.flagAverageDim) = 1; - end - - - function out = get.sqzDims(this) - out = this.dataDims(this.dataSize>1); - end - - - function out = get.sqzSize(this) - out = this.dataSize(this.dataSize>1); - end - - - function set.flagRemoveOS(this,val) - % set method for removeOS - this.arg.removeOS = logical(val); - end - - function out = get.flagRemoveOS(this) - out = this.arg.removeOS; - end + function out = get.flagRemoveOS(this) + out = this.arg.removeOS; + end - function set.flagDoAverage(this,val) - ix = ismember(this.dataDims, 'Ave'); - this.flagAverageDim(ix) = val; - end + function set.flagDoAverage(this,val) + ix = ismember(this.dataDims, 'Ave'); + this.flagAverageDim(ix) = val; + end - function out = get.flagDoAverage(this) - ix = ismember(this.dataDims, 'Ave'); - out = this.flagAverageDim(ix); - end - - function set.flagAverageReps(this,val) - ix = ismember(this.dataDims, 'Rep'); - this.flagAverageDim(ix) = val; - end + function out = get.flagDoAverage(this) + ix = ismember(this.dataDims, 'Ave'); + out = this.flagAverageDim(ix); + end - function out = get.flagAverageReps(this) - ix = ismember(this.dataDims, 'Rep'); - out = this.flagAverageDim(ix); - end + function set.flagAverageReps(this,val) + ix = ismember(this.dataDims, 'Rep'); + this.flagAverageDim(ix) = val; + end - function set.flagAverageSets(this,val) - ix = ismember(this.dataDims, 'Set'); - this.flagAverageDim(ix) = val; - end + function out = get.flagAverageReps(this) + ix = ismember(this.dataDims, 'Rep'); + out = this.flagAverageDim(ix); + end - function out = get.flagAverageSets(this) - ix = ismember(this.dataDims, 'Set'); - out = this.flagAverageDim(ix); - end + function set.flagAverageSets(this,val) + ix = ismember(this.dataDims, 'Set'); + this.flagAverageDim(ix) = val; + end - - function set.flagIgnoreSeg(this,val) - ix = ismember(this.dataDims, 'Seg'); - this.flagAverageDim(ix) = val; - end + function out = get.flagAverageSets(this) + ix = ismember(this.dataDims, 'Set'); + out = this.flagAverageDim(ix); + end - function out = get.flagIgnoreSeg(this) - ix = ismember(this.dataDims, 'Seg'); - out = this.flagAverageDim(ix); - end - - - function set.flagSkipToFirstLine(this,val) - val = logical(val); - if val ~= this.arg.skipToFirstLine - this.arg.skipToFirstLine = val; - - if this.arg.skipToFirstLine - this.skipLin = min(this.Lin)-1; - this.skipPar = min(this.Par)-1; - else - this.skipLin = 0; - this.skipPar = 0; - end - NLinAlloc = max(1, this.NLin - this.skipLin); - NParAlloc = max(1, this.NPar - this.skipPar); - this.fullSize(3:4) = [NLinAlloc NParAlloc]; + + function set.flagIgnoreSeg(this,val) + ix = ismember(this.dataDims, 'Seg'); + this.flagAverageDim(ix) = val; end - end - function out = get.flagSkipToFirstLine(this) - out = this.arg.skipToFirstLine; - end + function out = get.flagIgnoreSeg(this) + ix = ismember(this.dataDims, 'Seg'); + out = this.flagAverageDim(ix); + end - function out = get.flagRampSampRegrid(this) - out = this.arg.rampSampRegrid; - end - function set.flagRampSampRegrid(this, val) - val = logical(val); - if (val == true && isempty(this.rampSampTrj)) - error('No trajectory for regridding available'); + function set.flagSkipToFirstLine(this,val) + val = logical(val); + if val ~= this.arg.skipToFirstLine + this.arg.skipToFirstLine = val; + + if this.arg.skipToFirstLine + this.skipLin = min(this.Lin)-1; + this.skipPar = min(this.Par)-1; + else + this.skipLin = 0; + this.skipPar = 0; + end + NLinAlloc = max(1, this.NLin - this.skipLin); + NParAlloc = max(1, this.NPar - this.skipPar); + this.fullSize(3:4) = [NLinAlloc NParAlloc]; + end end - this.arg.rampSampRegrid = val; - end - %SRY: accessor methods for raw data correction - function out = get.flagDoRawDataCorrect(this) - out = this.arg.doRawDataCorrect; - end - function set.flagDoRawDataCorrect(this, val) - val = logical(val); - if (val == true && strcmp(this.softwareVersion, 'vd')) - error('raw data correction for VD not supported/tested yet'); + function out = get.flagSkipToFirstLine(this) + out = this.arg.skipToFirstLine; end - this.arg.doRawDataCorrect = val; - end + function out = get.flagRampSampRegrid(this) + out = this.arg.rampSampRegrid; + end - function out = get.RawDataCorrectionFactors(this) - out = this.arg.rawDataCorrectionFactors; - end + function set.flagRampSampRegrid(this, val) + val = logical(val); + if (val == true && isempty(this.rampSampTrj)) + error('No trajectory for regridding available'); + end + this.arg.rampSampRegrid = val; + end - function set.RawDataCorrectionFactors(this, val) - %this may not work if trying to set the factors before NCha has - %a meaningful value (ie before calling clean) - if (~isrow(val) || length(val) ~= this.NCha) - error('RawDataCorrectionFactors must be a 1xNCha row vector'); + %SRY: accessor methods for raw data correction + function out = get.flagDoRawDataCorrect(this) + out = this.arg.doRawDataCorrect; end - this.arg.rawDataCorrectionFactors = val; - end -end - -methods (Access='protected') - % helper functions - - function fid = fileopen(this) - % look out for unlikely event that someone is switching between - % windows and unix systems: - [path,name,ext] = fileparts(this.filename); - this.filename = fullfile(path,[name ext]); - - % test access - if numel(dir(this.filename))==0 - % update path when file of same name can be found in current - % working dir. -- otherwise throw error - [oldpath,name,ext] = fileparts(this.filename); - newloc = fullfile(pwd,[name ext]); - if numel(dir(newloc))==1 - fprintf('Warning: File location updated from "%s" to current working directory.\n',oldpath); - this.filename = newloc; - else - error(['File "' this.filename '" not found.']); + function set.flagDoRawDataCorrect(this, val) + val = logical(val); + if (val == true && strcmp(this.softwareVersion, 'vd')) + error('raw data correction for VD not supported/tested yet'); end + + this.arg.doRawDataCorrect = val; end - fid = fopen(this.filename); - end - - function [selRange,selRangeSz,outSize] = calcRange(this,S) + function out = get.RawDataCorrectionFactors(this) + out = this.arg.rawDataCorrectionFactors; + end - switch S.type - case '()' - bSqueeze = false; - case '{}' - bSqueeze = true; + function set.RawDataCorrectionFactors(this, val) + %this may not work if trying to set the factors before NCha has + %a meaningful value (ie before calling clean) + if (~isrow(val) || length(val) ~= this.NCha) + error('RawDataCorrectionFactors must be a 1xNCha row vector'); + end + this.arg.rawDataCorrectionFactors = val; end - selRange = num2cell(ones(1,numel(this.dataSize))); - outSize = ones(1,numel(this.dataSize)); + end - if ( isempty(S.subs) || strcmpi(S.subs(1),'') ) - % obj(): shortcut to select all data - % unfortunately, matlab does not allow the statement - % obj{}, so we can't use it... - % alternative: obj{''} (obj('') also works) - for k=1:numel(this.dataSize) - selRange{k} = 1:this.dataSize(k); + methods (Access='protected') + % helper functions + + function fid = fileopen(this) + % look out for unlikely event that someone is switching between + % windows and unix systems: + [path,name,ext] = fileparts(this.filename); + this.filename = fullfile(path,[name ext]); + + % test access + if numel(dir(this.filename))==0 + % update path when file of same name can be found in current + % working dir. -- otherwise throw error + [oldpath,name,ext] = fileparts(this.filename); + newloc = fullfile(pwd,[name ext]); + if numel(dir(newloc))==1 + fprintf('Warning: File location updated from "%s" to current working directory.\n',oldpath); + this.filename = newloc; + else + error(['File "' this.filename '" not found.']); + end end - if ~bSqueeze - outSize = this.dataSize; - else - outSize = this.sqzSize; + fid = fopen(this.filename); + end + + + function [selRange,selRangeSz,outSize] = calcRange(this,S) + + switch S.type + case '()' + bSqueeze = false; + case '{}' + bSqueeze = true; end - else - for k=1:numel(S.subs) + + selRange = num2cell(ones(1,numel(this.dataSize))); + outSize = ones(1,numel(this.dataSize)); + + if ( isempty(S.subs) || strcmpi(S.subs(1),'') ) + % obj(): shortcut to select all data + % unfortunately, matlab does not allow the statement + % obj{}, so we can't use it... + % alternative: obj{''} (obj('') also works) + for k=1:numel(this.dataSize) + selRange{k} = 1:this.dataSize(k); + end if ~bSqueeze - cDim = k; % nothing to do + outSize = this.dataSize; else - % we need to rearrange selRange from squeezed - % to original order - cDim = find(strcmp(this.dataDims,this.sqzDims{k}) == 1); + outSize = this.sqzSize; end - if strcmp(S.subs{k},':') - if k this.dataSize(k) + error('selection out of range'); end - elseif isnumeric(S.subs{k}) - selRange{cDim} = single(S.subs{k}); - else - error('unknown string in brackets (e.g. 1:end does not work here)'); end - outSize(k) = numel(selRange{cDim}); end + + selRangeSz = ones(1,numel(this.dataSize)); for k=1:numel(selRange) - if max(selRange{k}) > this.dataSize(k) - error('selection out of range'); - end + selRangeSz(k) = numel(selRange{k}); end - end - selRangeSz = ones(1,numel(this.dataSize)); - for k=1:numel(selRange) - selRangeSz(k) = numel(selRange{k}); - end - - % now select all indices for the dims that are averaged - for k = find(this.flagAverageDim) - selRange{k} = 1:this.fullSize(k); + % now select all indices for the dims that are averaged + for k = find(this.flagAverageDim) + selRange{k} = 1:this.fullSize(k); + end end - end - - function [ixToRaw, ixToTarget] = calcIndices(this) - % calculate indices to target & source(raw) - LinIx = this.Lin - this.skipLin; - ParIx = this.Par - this.skipPar; - sz = this.fullSize(3:end); % extra variable needed for octave compatibility - ixToTarget = this.sub2ind_double(sz,... - LinIx, ParIx, this.Sli, this.Ave, this.Phs, this.Eco,... - this.Rep, this.Set, this.Seg, this.Ida, this.Idb,... - this.Idc, this.Idd, this.Ide); - - % now calc. inverse index (page table: virtual to physical addresses) - % indices of lines that are not measured are zero - ixToRaw = zeros(1,prod(this.fullSize(3:end)),'double'); - - ixToRaw(ixToTarget) = 1:numel(ixToTarget); - end -end - -methods (Static) - % helper functions, accessible from outside via .function() - % without an instance of the class. - - function ro_shift = calcOffcenterShiftRO(slicedata) - % calculate ro offcenter shift from mdh's slicedata - - % slice position - pos = slicedata(1:3); - - %quaternion - a = slicedata(5); - b = slicedata(6); - c = slicedata(7); - d = slicedata(4); - - read_dir = zeros(3,1); - read_dir(1) = 2 * (a * b - c * d); - read_dir(2) = 1 - 2 * (a^2 + c^2); - read_dir(3) = 2 * (b * c + a * d); - - ro_shift = dot(pos, read_dir); - end - - function ndx = sub2ind_double(sz, varargin) - %SUB2IND_double Linear index from multiple subscripts. - % Works like sub2ind but always returns double - % also slightly faster, but no checks - %======================================== - sz = double(sz); - ndx = double(varargin{end}) - 1; - for i = length(sz)-1:-1:1 - ix = double(varargin{i}); - ndx = sz(i)*ndx + ix-1; + + function [ixToRaw, ixToTarget] = calcIndices(this) + % calculate indices to target & source(raw) + LinIx = this.Lin - this.skipLin; + ParIx = this.Par - this.skipPar; + sz = this.fullSize(3:end); % extra variable needed for octave compatibility + ixToTarget = this.sub2ind_double(sz,... + LinIx, ParIx, this.Sli, this.Ave, this.Phs, this.Eco,... + this.Rep, this.Set, this.Seg, this.Ida, this.Idb,... + this.Idc, this.Idd, this.Ide); + + % now calc. inverse index (page table: virtual to physical addresses) + % indices of lines that are not measured are zero + ixToRaw = zeros(1,prod(this.fullSize(3:end)),'double'); + + ixToRaw(ixToTarget) = 1:numel(ixToTarget); end - ndx = ndx + 1; - end % sub2ind_double - - function N = cast2MinimalUint( N ) - Nmax = max( reshape(N,[],1) ); - Nmin = min( reshape(N,[],1) ); - if Nmin < 0 || Nmax > intmax('uint64') - return + + function ix = getRawIx(this, varargin) + % ugly helper function + % index algebra seriously needs some reorganization + sel = ones(1,4); + sel(1:nargin-1) = varargin{:}; + sel = num2cell(sel); + [sli, eco, rep, set] = sel{:}; + + % selection order: sli, eco, rep, set + S.type = '()'; + S.subs = num2cell(ones(1,16)); + + S.subs{5} = sli; + S.subs{8} = eco; + S.subs{9} = rep; + S.subs{10} = set; + + backup_ave_flags = this.flagAverageDim; + this.flagAverageDim = true(1,16); + this.flagAverageDim(5) = false; + this.flagAverageDim(8:10) = false; + + [selRange, ~, ~] = this.calcRange(S); + ixToRaw = this.calcIndices; + tmp = reshape(1:prod(double(this.fullSize(3:end))), this.fullSize(3:end)); + tmp = tmp(selRange{3:end}); + ix = nonzeros(ixToRaw(tmp)); + ix = ix(1); + + this.flagAverageDim = backup_ave_flags; end - if Nmax > intmax('uint32') - idxClass = 'uint64'; - elseif Nmax > intmax('uint16') - idxClass = 'uint32'; - else - idxClass = 'uint16'; + function pos = getSlicePos(this, varargin) + pos = this.slicePos(:,this.getRawIx(varargin{:})); end + end + + methods (Static) + % helper functions, accessible from outside via .function() + % without an instance of the class. + + function ndx = sub2ind_double(sz, varargin) + %SUB2IND_double Linear index from multiple subscripts. + % Works like sub2ind but always returns double + % also slightly faster, but no checks + %======================================== + sz = double(sz); + ndx = double(varargin{end}) - 1; + for i = length(sz)-1:-1:1 + ix = double(varargin{i}); + ndx = sz(i)*ndx + ix-1; + end + ndx = ndx + 1; + end % sub2ind_double + + function N = cast2MinimalUint( N ) + Nmax = max( reshape(N,[],1) ); + Nmin = min( reshape(N,[],1) ); + if Nmin < 0 || Nmax > intmax('uint64') + return + end + + if Nmax > intmax('uint32') + idxClass = 'uint64'; + elseif Nmax > intmax('uint16') + idxClass = 'uint32'; + else + idxClass = 'uint16'; + end - N = cast( N, idxClass ); - end % cast2MinimalUint() + N = cast( N, idxClass ); + end % cast2MinimalUint() -end % of methods (Static) + end % of methods (Static) end % classdef + + +