From 711fafc9878f9bc16d60ca06f796d5e70f59b002 Mon Sep 17 00:00:00 2001 From: "ryan.topfer@polymtl.ca" Date: Wed, 19 Aug 2020 15:41:48 -0400 Subject: [PATCH 1/4] BUG: Construct MaRdI object based on discovered image volumes, and throw warning when different from dcm header (series repetition number, i.e. number of measurements) --- Img/MaRdI.m | 65 +++++++++++++++++++++++++++-------------------------- 1 file changed, 33 insertions(+), 32 deletions(-) diff --git a/Img/MaRdI.m b/Img/MaRdI.m index 6c723263..497c5914 100644 --- a/Img/MaRdI.m +++ b/Img/MaRdI.m @@ -32,10 +32,6 @@ % .Hdrs % Cell array of (truncated) DICOM headers courtesy of dicominfo(). % (One entry for every image) -% -% ========================================================================= -% Author::ryan.topfer@polymtl.ca -% ========================================================================= %% ======================================================================== % @@ -118,12 +114,6 @@ nRows = SpecialHdr.Height ; nColumns = SpecialHdr.Width ; - if ~myisfield( SpecialHdr.MrProt, 'lRepetitions' ) - nVolumes = 1 ; - else - nVolumes = SpecialHdr.MrProt.lRepetitions + 1 ; - end - % containers for a few terms varying between dicoms: sliceLocations = zeros( nImages, 1 ) ; % [units: mm] echoTimes = zeros( nImages, 1 ) ; % [units: ms] @@ -151,8 +141,34 @@ echoTimes = sort( unique( echoTimes ), 1, 'ascend' ) ; nEchoes = length( echoTimes ) ; + if ~myisfield( SpecialHdr.MrProt, 'lRepetitions' ) + nVolumesAcquired = 1 ; + else + nVolumesAcquired = SpecialHdr.MrProt.lRepetitions + 1 ; + end + + if nImages == nSlices*nEchoes*nVolumesAcquired + nVolumes = nVolumesAcquired; + else + acquisitionTimes = zeros(nImages, 1); + + for iImg = 1 : nImages + acquisitionTimes(iImg) = str2num(RawHdrs{iImg}.AcquisitionTime); + end + + nVolumes = numel(unique(acquisitionTimes))/(nSlices*nEchoes); + msg = sprintf( [ ... + 'Missing images or invalid DICOM header.' ... + '\n*****\n' ... + 'Of the ' num2str(nVolumesAcquired) ' image repetitions (volume measurements) ' ... + 'claimed by the header, only ' num2str(nVolumes) ' were found.\n' ... + 'This scenario is untested and may cause errors.' ... + '\n*****\n' ] ); + warning(msg); + end + Img.img = zeros( nRows, nColumns, nSlices, nEchoes, nVolumes ) ; - + % copy of RawHdrs to be reorganized: Img.Hdrs = cell( nSlices, nEchoes, nVolumes ) ; @@ -597,16 +613,9 @@ % sampling, beginning at the k_min periphery, has been considered in the % current implementation! -nSlices = Img.getnumberofslices() ; -nEchoes = length( Img.getechotime() ) ; - -if myisfield( Img.Hdr.MrProt, 'lRepetitions' ) - nMeasurements = ( Img.Hdr.MrProt.lRepetitions + 1) ; -else - nMeasurements = 1; -end - -assert( nMeasurements == size( Img.img, 5 ), 'Invalid .Hdr' ) ; +nSlices = Img.getnumberofslices() ; +nEchoes = length( Img.getechotime() ) ; +nMeasurements = Img.getnumberofmeasurements(); tAcq = Img.getacquisitiontime() ; @@ -746,17 +755,9 @@ % For EPI MOSAIC (which has a single AcquisitionTime value for each volume), % t( iSlice ) = AcquisitionTime (first slice) + iSlice*(volumeTR/nSlices) -nSlices = Img.getnumberofslices() ; -nEchoes = length( Img.getechotime() ) ; - -if myisfield( Img.Hdr.MrProt, 'lRepetitions' ) - nMeasurements = ( Img.Hdr.MrProt.lRepetitions + 1) ; -else - nMeasurements = 1; -end - -%assert( nMeasurements == size( Img.img, 5 ), 'Invalid .Hdr' ) ; % EAO: -%commented this out so I could save a time series +nSlices = Img.getnumberofslices() ; +nEchoes = length( Img.getechotime() ) ; +nMeasurements = Img.getnumberofmeasurements(); t = zeros( nSlices, nEchoes, nMeasurements ) ; From 01f7ac31a8cbdc6d1d5d45a42f9306a2f2e34fd9 Mon Sep 17 00:00:00 2001 From: "ryan.topfer@polymtl.ca" Date: Wed, 19 Aug 2020 18:08:23 -0400 Subject: [PATCH 2/4] BUG: MaRdI.resliceimg(), fix erroneous evaluation of grid dimensions. Method now returns the correct array size but the extrapolation is still off (likely due to the ROI/mask assignment on which it is based) --- Img/MaRdI.m | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/Img/MaRdI.m b/Img/MaRdI.m index 497c5914..c87feb4c 100644 --- a/Img/MaRdI.m +++ b/Img/MaRdI.m @@ -1250,16 +1250,25 @@ error('Dimensions of input Img.img must be >= 2, and <= 5') ; end -if ndims( X_Ep ) == 2 % interpolating down to 2d single-slice - gridSizeEp = [ size(X_Ep) 1 ] ; -elseif ndims( X_Ep ) == 3 - gridSizeEp = size(X_Ep) ; -else - error('Expected 2d or 3d target interpolation grid') +% ----------------- +%% Define variables + +switch ndims(X_Ep) + case 2 + % interpolating down to a single-slice + gridSizeEp = [ size(X_Ep) 1 ] ; + case 3 + % interpolating to a 3d volume + gridSizeEp = size(X_Ep) ; + otherwise + % sanity check + error('Expected 2d or 3d target interpolation grid') end -%% ----- -% Define variables +% sanity check +assert( [all(gridSizeEp>=1)] & [nnz(gridSizeEp==1)<=1], ... + 'Unexpected result for size of grid coordinates:' ); + gridSizeIp = size( X_Ip ) ; if myisfieldfilled( Img.Hdr, 'MaskingImage' ) @@ -1268,11 +1277,6 @@ maskIp = true( gridSizeIp ) ; end -gridSizeEp = size( X_Ep ) ; -if ndims(gridSizeEp) == 2 - gridSizeEp = [gridSizeEp 1] ; -end - if ~exist('maskEp') if ~isUsingScatteredInterpolant warning('No logical mask provided: For faster results, restrict the target/output voxels to those of interest by providing this mask!') ; From 0bb44b6d60b46e5e1eb23605ed1c40fcc17a0762 Mon Sep 17 00:00:00 2001 From: "ryan.topfer@polymtl.ca" Date: Wed, 19 Aug 2020 18:47:08 -0400 Subject: [PATCH 3/4] =?UTF-8?q?BUG:=20Fix=20overwrite=20of=20gridSize=20va?= =?UTF-8?q?riable=E2=80=94what=20should=20be=20a=20const=20(only=20matlab?= =?UTF-8?q?=20does=20not=20support=20them...)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Img/MaRdI.m | 34 +++++++++++++--------------------- 1 file changed, 13 insertions(+), 21 deletions(-) diff --git a/Img/MaRdI.m b/Img/MaRdI.m index c87feb4c..cc0f4a9b 100644 --- a/Img/MaRdI.m +++ b/Img/MaRdI.m @@ -1235,33 +1235,26 @@ end end -isUsingScatteredInterpolant = [] ; +assert( [ndims(Img.img) >= 2] & [ndims(Img.img) <= 5],... + 'Dimensions of input Img.img must be >= 2, and <= 5') ; -if (ndims(Img.img) > 1) && (ndims(Img.img) <= 5) +% ----------------- +%% Define variables +gridSizeIp = Img.getgridsize(); - gridSizeIp = Img.getgridsize() ; - - if gridSizeIp(3) > 1 - isUsingScatteredInterpolant = true ; - else - isUsingScatteredInterpolant = false ; - end +if gridSizeIp(3) > 1 + isUsingScatteredInterpolant = true ; else - error('Dimensions of input Img.img must be >= 2, and <= 5') ; + isUsingScatteredInterpolant = false ; end -% ----------------- -%% Define variables switch ndims(X_Ep) - case 2 - % interpolating down to a single-slice + case 2 % interpolating down to a single-slice gridSizeEp = [ size(X_Ep) 1 ] ; - case 3 - % interpolating to a 3d volume + case 3 % interpolating to a 3d volume gridSizeEp = size(X_Ep) ; - otherwise - % sanity check + otherwise % sanity check error('Expected 2d or 3d target interpolation grid') end @@ -1269,8 +1262,6 @@ assert( [all(gridSizeEp>=1)] & [nnz(gridSizeEp==1)<=1], ... 'Unexpected result for size of grid coordinates:' ); -gridSizeIp = size( X_Ip ) ; - if myisfieldfilled( Img.Hdr, 'MaskingImage' ) maskIp = logical( sum( sum( Img.Hdr.MaskingImage, 5 ), 4 ) ) ; else @@ -1279,7 +1270,8 @@ if ~exist('maskEp') if ~isUsingScatteredInterpolant - warning('No logical mask provided: For faster results, restrict the target/output voxels to those of interest by providing this mask!') ; + warning(['No logical mask provided: For faster results, ' ... + 'restrict the target/output voxels to those of interest by providing this mask!']) ; end maskEp = true( gridSizeEp ) ; end From 206a2442de56323fb7ff87568062913642883f12 Mon Sep 17 00:00:00 2001 From: "ryan.topfer@polymtl.ca" Date: Wed, 19 Aug 2020 19:08:24 -0400 Subject: [PATCH 4/4] MaRdI.resliceimg(): Add check on array dimensions of ROI extrapolation mask (must match dims of extrapolation coordinate arrays) --- Img/MaRdI.m | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/Img/MaRdI.m b/Img/MaRdI.m index cc0f4a9b..f9a3c3af 100644 --- a/Img/MaRdI.m +++ b/Img/MaRdI.m @@ -1274,6 +1274,15 @@ 'restrict the target/output voxels to those of interest by providing this mask!']) ; end maskEp = true( gridSizeEp ) ; +else + if gridSizeEp(3) > 1 + gridSizeMaskEp = size(maskEp) + else + gridSizeMaskEp = [size(maskEp) 1] + end + assert( isequal(gridSizeMaskEp, gridSizeEp), ... + [ 'Dimensions of the arrays of extrapolation coordinate ' ... + 'and the assigned extrapolation ROI mask must be identical' ] ); end iEp = find( maskEp(:) ) ; % indices of target voxels