diff --git a/README.md b/README.md index da98259..06c0fd5 100644 --- a/README.md +++ b/README.md @@ -6,9 +6,14 @@ A Cartesian 3D random Gaussian ESPReSSo subsampling is proposed which acquires a - precompiled 4D MR imaging sequence: CS_Retro (Siemens, VB20P) ## Reconstruction +retrospective motion gating and reconstruction of subsampled data - self-navigation signal extraction - gating procedure - for the CS reconstruction, please refer to https://github.com/thomaskuestner/CS_LAB +### Matlab +code included here +### Gadgetron +please refer to https://github.com/martinschwartz/CS_Gadgetron -------------------------------------------------------- Please read LICENSE file for licensing details. diff --git a/acquisition/CS_FLASH_retro.dll b/acquisition/CS_FLASH_retro.dll new file mode 100644 index 0000000..2863795 Binary files /dev/null and b/acquisition/CS_FLASH_retro.dll differ diff --git a/acquisition/CS_FLASH_retro.i86 b/acquisition/CS_FLASH_retro.i86 new file mode 100644 index 0000000..85ce932 Binary files /dev/null and b/acquisition/CS_FLASH_retro.i86 differ diff --git a/acquisition/CS_FLASH_retro.i87 b/acquisition/CS_FLASH_retro.i87 new file mode 100644 index 0000000..ad04c2c Binary files /dev/null and b/acquisition/CS_FLASH_retro.i87 differ diff --git a/reconstruction/fRetroGating.m b/reconstruction/fRetroGating.m new file mode 100644 index 0000000..fc79700 --- /dev/null +++ b/reconstruction/fRetroGating.m @@ -0,0 +1,499 @@ +function [dKSpace, dRespPhases] = fRetroGating(sFilename, iNPhases, dTolerance, dTime, sGatingMode, sGatingGather, sPrecision, iEcho, lDebug) +%FRETROGATING +% +% Copyright 2014-2016 Christian Wuerslin, University of Tuebingen, Germany +% christian.wuerslin@med.uni-tuebingen.de +% and Thomas Kuestner, University of Tuebingen, Germany +% thomas.kuestner@med.uni-tuebingen.de + +% inputs +% sFilename path to measurement file +% iNPhases number of respiratory gates +% dTolerance view sharing blending factor b of gates (1 (non-overlapping) <= b <= 2 (completely overlapping)) +% dTime crop scan time (0 := no cropping | scan time [s]) +% sGatingMode select gating procedure +% 'percentile' select gate centroids according to 5% and 95% percentile +% 'kmeans' perform k-means clustering to select gate centroids +% sGatingGather sample gathering mode: +% 'closest' take closest sample to gate centroid +% 'average' take all possible samples and weight them +% according to their distance to the gate centroid +% 'collect' gather all samples and let the +% reconstruction decide for the most favourable one +% 'none' no k-Space population, show simulations of the k-Space sampling density +% iEcho current to be processed echo +% lDebug debugging flag for additional plots +% +% outputs +% dKSpace gated kSpace +% dRespPhases respiratory phases (TR level) + +if nargin < 9, lDebug = false; end +if nargin < 8, iEcho = 1; end +if nargin < 7, sPrecision = 'single'; end +if nargin < 6, sGatingGather = 'closest'; end +if nargin < 5, sGatingMode = 'percentile'; end +if nargin < 4, dTime = 0; end +if nargin < 3, dTolerance = 1; end +if nargin < 2, iNPhases = 1; end + +% parse inputs +[sPath, sName, sExt] = fileparts(sFilename); +if isempty(sExt), sFilename = [sFilename, '.dat']; end +% sFilename = which(sFilename); +[sPath, sName] = fileparts(sFilename); +sFilename = [sPath, filesep, sName]; +iNPhasesLoop = 1:iNPhases; + +% ------------------------------------------------------------------------- +% Get some relevant data from the drecksMDH +drecksMDH = fMeas_mainLUTDrecks( sFilename ); +iNPartitions = drecksMDH.Geo.MatrixSize(3); +iNLines = drecksMDH.Geo.MatrixSize(2); +iBaseRes = drecksMDH.Geo.MatrixSize(1); +dTR = drecksMDH.Contrast.TR./1E6; +dNavPeriod = drecksMDH.Wip.NavPeriod; +if(isfield(drecksMDH.Geo,'EchoPosition')) + iEchoLine = drecksMDH.Geo.EchoPosition(1); + iEchoPartition = drecksMDH.Geo.EchoPosition(2); +else + iEchoLine = drecksMDH.Geo.MatrixSize(2)/2; + iEchoPartition = drecksMDH.Geo.MatrixSize(3)/2; +end +if(isfield(drecksMDH,'Wip')) + iNavPERes = drecksMDH.Wip.NavRes(1); + iNav3DRes = drecksMDH.Wip.NavRes(2); +else + iNavPERes = 8; + iNav3DRes = 1; +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Get loop counters from file and reduce to first channel and first echo +load([sFilename, '.mat'], 'iLC'); +load([sFilename, '.mat'], 'iSP'); +load([sFilename, '.mat'], 'dNavInt'); +iNChannels = max(unique(double(iLC(:,2)))); + +lMask = iLC(:, 1) == 2*iBaseRes & iLC(:, 11) < 65000 & iLC(:,7) == iEcho-1; % throw out dummy measurements (which are needed to capture surrogate signals and keep steady state) +iSP = iSP(lMask); +iLCReduced = iLC(lMask, :); % only "real" readouts +iSP = iSP(1:iNChannels:end,:); +iLCReduced = iLCReduced(1:iNChannels:end, :); + +fprintf(1, 'Scantime was %f3 s\n', length(iLCReduced)*dTR); + +% ------------------------------------------------------------------------- +% check cropping parameter dTime +dPRECROP = 5.0; % Crop first 5 s +fprintf(1, 'Cropped initial %.3f s for steady-state!\n', dPRECROP); +if(isscalar(dTime)) + dTimeStart = dPRECROP; + dTimeEnd = dTime; + dTime = num2str(dTime,'%03d'); +else + dTimeStart = dTime(:,1); + dTimeEnd = dTime(:,2); + dTime = dTime.'; % for backward compatibility + dTime = strrep(num2str(dTime(:).','%03d'),' ',''); +end + +if(dTimeStart(1) < dPRECROP) + dTimeStart(1) = dPRECROP; +end +if(any(dTimeEnd == 0)) + dTimeEnd(dTimeEnd == 0) = length(iLCReduced)*dTR; +end +if(dTimeEnd(end) > length(iLCReduced)*dTR) + dTimeEnd(end) = length(iLCReduced)*dTR; +end + +iCropRange = []; +for iI = 1:length(dTimeStart) + iCropRange = [iCropRange,round([dTimeStart(iI)/dTR:dTimeEnd(iI)/dTR])]; + fprintf(1, 'Scantime cutout from %.3f s to %.3f s\n', dTimeStart(iI),dTimeEnd(iI)); +end +iCropRange = unique(iCropRange); + +iLCReduced = iLCReduced(iCropRange,:); +iSP = iSP (iCropRange,:); +dNavInt = dNavInt (iCropRange); + +iLine = iLCReduced(:,3); +iPartition = iLCReduced(:,6); + +iADCLength = 128/4 + iBaseRes*4; % MDH + ADC + +if(usejava('jvm') && ~feature('ShowFigureWindows')) + % use text-based alternative + lflagDisp = false; +else + lflagDisp = true; +end + +% ------------------------------------------------------------------------- +% respiratory gating +% ------------------------------------------------------------------------- +if(all(iNPhasesLoop > 0) && exist('dNavInt', 'var')) + % Calculate the center positions of the phases + % --------------------------------------------------------------------- + % distinguish between inhale and exhale + iNPhasesOld = iNPhases; + + % gating method + dNavInt = dNavInt(1:length(iLine)); + dNavMin = min(dNavInt); + dNavMax = max(dNavInt); + dToleranceOld = dTolerance; + if(dNavMin == dNavMax) % BH acquisition + dGatePos = dNavMin; + iNPhases = 1; + else + switch sGatingMode + case 'percentile' + % 5% and 95% percentile gating + dX = linspace(dNavMin, dNavMax, 256); + iHist = hist(dNavInt, dX); + iSum = sum(iHist); + + iLowerLim = 1; + while sum(iHist(iLowerLim:end))/iSum > 0.9 + iLowerLim = iLowerLim + 1; + end + iUpperLim = 256; + while sum(iHist(1:iUpperLim))/iSum > 0.9 + iUpperLim = iUpperLim - 1; + end + + % equally spaced gate centroids + dMin = dX(iLowerLim); + dMax = dX(iUpperLim); + dGatePos = linspace(dMax, dMin, iNPhasesOld)'; + dTolerance = abs(diff(dGatePos)).*dTolerance./2; + if(isempty(dTolerance)) % breathhold + dTolerance = dToleranceOld; + end + dTolerance = repmat(dTolerance(1),[iNPhases, 1]); + + case 'kmeans' + % kmeans clustering of gate centroids + dGrad = gradient(dNavInt); + lExhale = dGrad >= 0; + lInhale = dGrad < 0; + iInd = find(lExhale,1,'first'):find(lInhale,1,'last'); % without outliers (due to dNav extraction) + lInd = false(1,length(dNavInt)); + lInd(iInd) = true; + [iClusterEx, dGatePosEx] = kmeans(dNavInt(lInd).',iNPhases); % without outliers (due to dNav extraction) + dGatePos = sort(dGatePosEx,'descend'); + + % get tolerance window width + dMinDist = zeros(length(dGatePos),2); % 1st col: upper dist, 2nd col: lower dist + for iI = 1:length(dGatePos) + lIndCl = lInd; + lIndCl(lInd) = iClusterEx == find(dGatePos(iI) == dGatePosEx); + dMinDist(iI,1) = max(dNavInt(lIndCl)); + dMinDist(iI,2) = min(dNavInt(lIndCl)); + end + dMinDist = abs(dMinDist - repmat(dGatePos, [1 2])); + dMaxDist = abs([[dMinDist(1,1); diff(dGatePos(1:iNPhasesOld))], [diff(dGatePos(1:iNPhasesOld)); dMinDist(iNPhasesOld,2)]]); + dTolerance = dMinDist + (dTolerance-1) .* dMaxDist; + end + end + % get respiratory phase mask + dRespPhases = zeros(length(dNavInt),1); + for iPh = iNPhasesLoop + if(strcmp(sGatingMode,'kmeans')) + iTol = [1 2]; + else + iTol = [1 1]; + end + lMask = dNavInt <= dGatePos(iPh) + dTolerance(iPh,iTol(1)) & dNavInt > dGatePos(iPh) - dTolerance(iPh,iTol(2)); + dRespPhases(lMask) = iPh; + end +else % ignore respiratory gating + dNavInt = ones(length(iLCReduced), 1)'; + iNPhases = 1; % needed in for loop + iNPhasesLoop = 1; + dGatePos = 1; + dTolerance = 1; + dRespPhases = ones(length(iLCReduced),1); +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Show the navigator signal, phase centers and tolerances +if lDebug + dColors = [0 1 0;0 0 1;1 0 0;0 0.965517241379310 1;1 0.482758620689655 0.862068965517241;1 0.862068965517241 0.482758620689655;0 0.551724137931035 1;0 0.413793103448276 0.0344827586206897;0.620689655172414 0.275862068965517 0.206896551724138;0.241379310344828 0 0.413793103448276;1 0.931034482758621 1]; + hfig_respSignal = figure('name',sprintf('Gating - Ph%02d, Tol%03.0f, t%s, Esp%s, InExh%s, %s, %s', iNPhases, 100*dToleranceOld, dTime, sGatingMode, sGatingGather)); + plot(dNavInt, 'Color', 'k'); + hold on; + dXMax = max(get(gca, 'XLim')); + for iI = 1:length(dGatePos) + if(strcmp(sGatingMode,'kmeans')) + line([0 dXMax], [dGatePos(iI) dGatePos(iI)], 'Color', dColors(iI,:), 'LineWidth', 1.5); + line([0 dXMax], [dGatePos(iI) + dTolerance(iI,1) dGatePos(iI) + dTolerance(iI,1)], 'Color', dColors(iI,:), 'LineStyle', '--'); + line([0 dXMax], [dGatePos(iI) - dTolerance(iI,2) dGatePos(iI) - dTolerance(iI,2)], 'Color', dColors(iI,:), 'LineStyle', '--'); + hold on; + lIndCl = lInd; + lIndCl(lInd) = iClusterEx == find(dGatePos(iI) == dGatePosEx); + if(iI < iNPhasesOld) + lIndVS = (iClusterEx == find(dGatePos(iI) == dGatePosEx) & iClusterEx == find(dGatePos(iI+1) == dGatePosEx)).'; % view sharing samples + else + lIndVS = false(1,length(lIndCl)); + end + lIndCl(lIndCl & lIndVS) = false; + + plot(find(lIndCl),dNavInt(lIndCl),'Color', dColors(iI,:), 'Marker', 'x', 'LineStyle', 'none'); + hold on; + plot(find(lIndVS),dNavInt(lIndVS),'Color', (dColors(iI,:)+dColors(iI+1,:))/2, 'Marker', 'x', 'LineStyle', 'none'); + hold on; + else % percentile + line([0 dXMax], [dGatePos(iI) dGatePos(iI)], 'Color', dColors(iI,:), 'LineWidth', 1.5); + line([0 dXMax], [dGatePos(iI) + dTolerance(iI) dGatePos(iI) + dTolerance(iI)], 'Color', dColors(iI,:), 'LineStyle', '--'); + line([0 dXMax], [dGatePos(iI) - dTolerance(iI) dGatePos(iI) - dTolerance(iI)], 'Color', dColors(iI,:), 'LineStyle', '--'); + hold on; + lIndCl = dNavInt <= dGatePos(iI) + dTolerance(iI) & dNavInt > dGatePos(iI) - dTolerance(iI); + if(iI == 1) + lIndVS = (dNavInt <= dGatePos(iI+1) + dTolerance(iI+1) & dNavInt > dGatePos(iI) - dTolerance(iI)); + elseif(iI == iNPhasesOld) + lIndVS = (dNavInt > dGatePos(iI-1) - dTolerance(iI-1) & dNavInt <= dGatePos(iI) + dTolerance(iI)); + else + lIndVS = (dNavInt <= dGatePos(iI+1) + dTolerance(iI+1) & dNavInt > dGatePos(iI) - dTolerance(iI)) | ... + (dNavInt > dGatePos(iI-1) - dTolerance(iI-1) & dNavInt <= dGatePos(iI) + dTolerance(iI)); + end + + lIndCl(lIndCl & lIndVS) = false; + + plot(find(lIndCl),dNavInt(lIndCl),'Color', dColors(iI,:), 'Marker', 'x', 'LineStyle', 'none'); + hold on; + plot(find(lIndVS),dNavInt(lIndVS),'Color', (dColors(iI,:)+dColors(iI+1,:))/2, 'Marker', 'x', 'LineStyle', 'none'); + hold on; + end + end + plot(dNavInt, 'Color', 'k'); +end +% ------------------------------------------------------------------------- + + +% ------------------------------------------------------------------------- +% Populate the k-Space +fprintf('*** k-Space statistics: ***\n'); +switch sGatingGather + + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + case 'closest' + if(lflagDisp), hF = waitbar(0, 'Populating k-Space'); end; + fid = fopen([sPath, filesep, sName, '.dat'], 'r'); + + % - - - - - - - - - - - - - - - - - - - - - - + % Pre-allocate k-Space memory + dKSpace = complex(zeros(iNLines, iBaseRes.*2, iNPartitions, length(iNPhasesLoop), iNChannels,sPrecision), ... + zeros(iNLines, iBaseRes.*2, iNPartitions, length(iNPhasesLoop), iNChannels,sPrecision)); + dKSpaceImg = zeros(iNLines, iNPartitions, length(iNPhasesLoop), 'uint8'); + % - - - - - - - - - - - - - - - - - - - - - - + + % - - - - - - - - - - - - - - - - - - - - - + % RESP PHASES loop + for iPh = iNPhasesLoop + lRespMask = dRespPhases == iPh; + dWeight = abs(dNavInt - dGatePos(iPh))'; + + % - - - - - - - - - - - - - - - + % LINES loop + for iL = 1:iNLines + lLineMask = iLine == (iL - 1); + + % - - - - - - - - - - - + % PARTITIONS loop + for iP = 1:iNPartitions + lMask = lRespMask & lLineMask & iPartition == (iP - 1); + dThisDist = dWeight(lMask); + if isempty(dThisDist), continue, end; + + [dVal, iInd] = min(dThisDist); + if(strcmp(sGatingMode,'kmeans')) + if(dNavInt(iInd) >= dGatePos(iPh)) + iTol = 1; + else + iTol = 2; + end; + else + iTol = 1; + end + if(dVal > dTolerance(iPh,iTol)) + continue; % just for backup + end + + iThisSP = iSP(lMask); + fseek(fid, double(iThisSP(iInd)), 'bof'); + dData = fread(fid, iADCLength.*iNChannels, 'float'); + dData = reshape(dData, [iADCLength, iNChannels]); + dData = dData(128/4 + 1:end,:); + dData = complex(dData(1:2:end,:), dData(2:2:end,:)); + dKSpace(iL, :, iP, iPh==iNPhasesLoop, :) = permute(dData, [3 1 4 2]); + dKSpaceImg(iL, iP, iPh==iNPhasesLoop) = 1; + end + % end of PARTITIONS loop + % - - - - - - - - - - - + if(lflagDisp), waitbar((double(iPh - 1).*double(iNLines) + double(iL))./(double(iNPhases).*double(iNLines)), hF); end; + + end + % end of LINES loop + % - - - - - - - - - - - - - - - + fprintf(1, ' Phase %u acceleration is %1.3f\n', iPh, double(numel(dKSpaceImg(:,:,iPh==iNPhasesLoop)))./double(sum(sum(dKSpaceImg(:,:,iPh==iNPhasesLoop) > 0)))); + + end + % end of PHASES loop + % - - - - - - - - - - - - - - - - - - - - - + fclose(fid); + if(lflagDisp), close(hF); end; + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + case 'collect' + if(lflagDisp), hF = waitbar(0, 'Populating k-Space'); end; + fid = fopen([sPath, filesep, sName, '.dat'], 'r'); + + % - - - - - - - - - - - - - - - - - - - - - - + % find maximal amount of samples per gate + dKSpaceImg = cell(iNLines, iNPartitions, iNPhases); + + % - - - - - - - - - - - - - - - - - - - - - - + % PHASES loop + for iPh = 1:iNPhases + if(strcmp(sGatingMode,'kmeans')) + iTol = [1 2]; + else + iTol = [1 1]; + end + dWeight = dNavInt <= dGatePos(iPh) + dTolerance(iPh,iTol(1)) & dNavInt > dGatePos(iPh) - dTolerance(iPh,iTol(2)); + dDist = abs(dNavInt - dGatePos(iPh))'; + + % - - - - - - - - - - - - - - - - + % LINES loop + for iL = 1:iNLines + lLineMask = iLine == (iL - 1); + + % - - - - - - - - - - - - + % PARTITIONS loop + for iP = 1:iNPartitions + lPartitionMask = iPartition == (iP - 1); + iInd = lLineMask & lPartitionMask & dWeight.'; + + if isempty(iInd), continue, end + [~,dThisDist] = sort(dDist(iInd)); + iInd = find(iInd); + iInd = iInd(dThisDist); + + dKSpaceImg{iL, iP, iPh} = iInd; + end + if(lflagDisp), waitbar((double(iPh - 1).*double(iNLines) + double(iL))./(double(iNPhases).*double(iNLines)), hF); end; + end + fprintf(1, ' Phase %u acceleration is %1.3f\n', iPh, (iNPartitions*iNLines)./(double(nnz(cellfun(@isempty,dKSpaceImg(:,:,iPh)))))); + end + if(lflagDisp), close(hF); end; + + % decide for the most occuring gate filling situations + % if any gate has more samples (k-space center), take the iNGate + % nearest samples + iSamples = cellfun(@length, dKSpaceImg); + n_elements = hist(iSamples(:),1:max(iSamples(:))); + iNGate = 1; + while(sum(n_elements(1:iNGate)) < 0.99*sum(n_elements(:))) + iNGate = iNGate + 1; + end + + lMask = iSamples <= iNGate & iSamples > 0; + lMask(iEchoLine-round(iNavPERes/2)+1:iEchoLine+round(iNavPERes/2), iEchoPartition-round(iNav3DRes/2)+1:iEchoPartition+round(iNav3DRes/2), :) = true; % keep navigator lines + lMask = find(~lMask); + for iMask = 1:length(lMask) + [cLine,cPar,cPha] = ind2sub(size(dKSpaceImg),lMask(iMask)); + dKSpaceImg{cLine,cPar,cPha} = []; + end + + % Pre-allocate k-Space memory + dKSpace = complex(zeros(iNLines, iBaseRes.*2, iNPartitions, iNPhases, iNGate, iNChannels, sPrecision), ... + zeros(iNLines, iBaseRes.*2, iNPartitions, iNPhases, iNGate, iNChannels, sPrecision)); + + if(lflagDisp), hF = waitbar(0, 'Populating k-Space'); end; + + % - - - - - - - - - - - - - - - - - - - - - - + % PHASES loop + for iPh = 1:iNPhases + + % - - - - - - - - - - - - - - - - + % LINES loop + for iL = 1:iNLines + + % - - - - - - - - - - - - + % PARTITIONS loop + for iP = 1:iNPartitions + if(isempty(dKSpaceImg{iL,iP,iPh})), continue; end; + + % - - - - - - - - - + % GATING loop + for iG = 1:length(dKSpaceImg{iL,iP,iPh}) + if(iG > iNGate), break, end; + iThisSP = iSP(dKSpaceImg{iL,iP,iPh}(iG)); + fseek(fid, double(iThisSP), 'bof'); + dData = fread(fid, iADCLength.*iNChannels, 'float'); + dData = reshape(dData, [iADCLength, iNChannels]); + dData = dData(128/4 + 1:end,:); + dData = complex(dData(1:2:end,:), dData(2:2:end,:)); + dKSpace(iL, :, iP, iPh, iG, :) = permute(dData, [3 1 4 2]); + end + % end of GATING loop + % - - - - - - - - - + end + % end of PARTITIONS loop + % - - - - - - - - - - - - + if(lflagDisp), waitbar((double(iPh - 1).*double(iNLines) + double(iL))./(double(iNPhases).*double(iNLines)), hF); end; + end + % end of LINES loop + % - - - - - - - - - - - - - - - - + end + % end of PHASES loop + % - - - - - - - - - - - - - - - - - - - - - - + fclose(fid); + if(lflagDisp), close(hF); end; + + + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + % No k-Space population, show simulations of the k-Space sampling + % density + case 'none' + dKSpaceImg = zeros(iNLines, iNPartitions, iNPhases, 'uint8'); + for iPh = 1:iNPhases + lMask = abs(dNavInt - dGatePos(iPh)) <= dTolerance; + + iThisLine = iLine(lMask); + iThisPartition = iPartition(lMask); + + for iI = 1:length(iThisLine) + dKSpaceImg(iThisLine(iI) + 1, iThisPartition(iI) + 1, iPh) = ... + dKSpaceImg(iThisLine(iI) + 1, iThisPartition(iI) + 1, iPh) + 1; + end + fprintf(1, ' Phase %u acceleration is %1.3f\n', iPh, double(numel(dKSpaceImg(:,:,iPh)))./double(sum(sum(dKSpaceImg(:,:,iPh) > 0)))); + end + dKSpace = dKSpaceImg; + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% show subsampling masks +if lDebug + dColMap = OptimalColor(256 + 20); + dColMap = dColMap(21:end, :); + dColMap(1,:) = [0, 0, 0]; + + for iPh = 1:iNPhases + figure, imshow(dKSpaceImg(:,:,iPh)', 'InitialMagnification', 400); + set(gcf, 'Name', sprintf('k-Space density for phase %u', iPh)); + colormap(dColMap); + set(gca, 'CLim', [0 255]); + end +end +% ------------------------------------------------------------------------- \ No newline at end of file diff --git a/reconstruction/fRetroGetNav.m b/reconstruction/fRetroGetNav.m new file mode 100644 index 0000000..ee700dd --- /dev/null +++ b/reconstruction/fRetroGetNav.m @@ -0,0 +1,153 @@ +function dNav = fRetroGetNav(sFilename, lDebug) +%FFLASHRETROGETNAV Extracts navigator data from CS_FLASH_retro sequence +% DNAV = FFLASHRETROGETNAV(SFILENAME) Extracts navigator data DNAV from +% the siemens meas-data file SFILENAME. Must have been parsed with +% FMeasCreateLUT before. +% +% See also: FMEASCREATELUT +% +% Copyright 2014-2016 Christian Wuerslin, University of Tuebingen, Germany +% christian.wuerslin@med.uni-tuebingen.de +% and Thomas Kuestner, University of Tuebingen, Germany +% thomas.kuestner@med.uni-tuebingen.de + +if(nargin < 1) + lDebug = false; +end +[sPath, sName, sExt] = fileparts(sFilename); + +% ------------------------------------------------------------------------- +% Get some relevant data from the drecksMDH +try + load(fullfile(sPath,[sName,'.mat']), 'SDrecksMDH'); + dBaseRes = SDrecksMDH.Geo.MatrixSize(1); + dTR = SDrecksMDH.Contrast.TR./1000; % in ms + if(isfield(SDrecksMDH.Geo,'EchoPosition')) + dEchoLine = SDrecksMDH.Geo.EchoPosition(1); + dEchoPartition = SDrecksMDH.Geo.EchoPosition(2); + else + dEchoLine = SDrecksMDH.Geo.MatrixSize(2)/2; + dEchoPartition = SDrecksMDH.Geo.MatrixSize(3)/2; + end + dNavPeriod = SDrecksMDH.Wip.NavPeriod; + dNavPERes = SDrecksMDH.Wip.NavRes(1); +catch + % Load the loop counters + load([sPath, filesep, sName, '.mat']); + dBaseRes = 256; + dEchoLine = 128; + dEchoPartition = 36; + dNavPeriod = 200; + dNavPERes = 1; +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Get ktSpace +[dKtSpace, iLC] = fMeasRead(sFilename, 'Set', [1 7], 'Eco', 0, 'Lin', dEchoLine, 'Par', dEchoPartition); +if(length(unique(iLC(:,10))) > 1) % in the case of switched on surrogate monitoring + lMask = true(size(iLC, 1), 1); + lMask = lMask & (iLC(:,10) == 1 | iLC(:,10) == 3 | iLC(:,10) == 5 | iLC(:,10) == 7); + iLC = iLC(lMask,:); + dKtSpace = dKtSpace(lMask,:); +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Determine data size +iNChannels = double(iLC(1, 2)); +iNSamples = size(dKtSpace, 2); +iNMeasurements = size(dKtSpace, 1)./(iNChannels); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Build up the kSpaces of all channels +dKtSpace = reshape(dKtSpace, [iNChannels, iNMeasurements, iNSamples]); +dKtSpace = permute(dKtSpace, [3, 2, 1]); % -> RO x t x CH +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Reconstruct the 1-D projections for all measuremetns and all channels +dImg = fftshift(ifft(ifftshift(dKtSpace(:,:,:)))); +dImg = flipdim(dImg, 1); % Invert the RO direction: 1-N -> H-F +dImg = dImg(iNSamples/4:iNSamples.*3/4 - 1, :, :, :); % RO x t x CH + +% Get x range of respiratory motion and sort out Channels with no relevant +% information in target area +dFreq = fft(dImg, [], 2); +dPower = dFreq.*conj(dFreq); +dIMGres = 1./(double(dNavPeriod)./1000.*double(iNMeasurements)); % The frequency resolution of dIMG in Hz +dPower = squeeze(sum(dPower(:, round(1./(5.*dIMGres)):round(1./(3.*dIMGres)), :), 2)); % RO x CH +dPowerInChan = dPower; % RO x CH +dPowerAcrossChannels = sum(dPowerInChan, 2); % RO x 1 +dPowerAcrossChannels(length(dPowerAcrossChannels).*3/4:end) = 0; % Prevent detection of regions in the abdomen +dPowerAcrossChannels = conv(dPowerAcrossChannels, fGaussianLP(20), 'same'); +[dMax, dX] = max(dPowerAcrossChannels); +dMax = max(dPowerInChan(dX - 20: dX + 20, :)); +lGoodChannels = dMax > 0.4.*max(dMax); + +dSOSImg = sqrt(sum(dImg(:,:,lGoodChannels).*conj(dImg(:,:,lGoodChannels)), 3)); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Get the navigator signal +iDisplacement = 80; % [mm] diaphragm displacement (max. +/- 80mm) +dDisplacement = round(iDisplacement/(SDrecksMDH.Geo.FOV(1)/SDrecksMDH.Geo.MatrixSize(1))); % [px] + +dRefImg = zeros(size(dSOSImg)); +dRefImg(:,end) = dSOSImg(:,end); +idx = 1:size(dRefImg,2); +dNav = zeros(1,size(dRefImg,2)); + +hF = waitbar(0, 'Getting Navigator'); +for i=size(dRefImg,2)-1:-1:2 + dRMSImg = zeros(size(dRefImg,1),2*dDisplacement+1); + tmp = dRefImg(:,idx(i+1:end)); + for iD = -dDisplacement:dDisplacement + dRMSImg(:,dDisplacement-iD+1) = sum((tmp - repmat(circshift(dSOSImg(:,idx(i)), iD),[1 size(tmp,2)])).^2,2); % RMS + waitbar((abs(i-size(dRefImg,2)-2)*(2*dDisplacement+1) + (iD + dDisplacement + 1))/((size(dRefImg,2)-2)*(2*dDisplacement+1)),hF); + end + [~, dNav(i)] = min(sum(dRMSImg(dX-round(dDisplacement/2):dX+round(dDisplacement/2),:))); % RMS + dNav(i) = dDisplacement + 1 - dNav(i); + dRefImg(:,idx(i)) = circshift(dSOSImg(:,idx(i)), dNav(i)); + waitbar(abs(i-size(dRefImg,2)-2)/(size(dRefImg,2)-2),hF); +end +close(hF); +dNav = -dNav; +dNav = conv(dNav, fGaussianLP(5), 'same'); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Visualize the result +if(lDebug) + figure('Name','Navi'), imagesc(dSOSImg(:, 3:end)); + colormap gray; + hold all; + plot(dNav(3:end)-mean(dNav)+dX); +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Interpolate the navigator data to TR intervals +load([sPath, filesep, sName, '.mat'], 'iLC'); +iLC = iLC(1:iNChannels:end, :); +iLC = iLC(iLC(:,7) == 0, :); % first echo +iLC = iLC(iLC(:, 1) == 2*dBaseRes, :); % only "real" readouts + +dNav = - dNav; +dNav = dNav - min(dNav); + +iNavInd = find((iLC(:,3) == dEchoLine) & (iLC(:,6) == dEchoPartition)); +if(length(iNavInd) ~= length(dNav)) % happens at the end when navi not fully acquired + iNavInd = iNavInd(1:length(dNav)); +end +dNavInt = interp1(iNavInd, dNav, 1:length(iLC), 'pchip'); +dNavInt_ms = interp1(iNavInd.*dTR, dNav, 1:length(iLC).*dTR, 'pchip'); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Append the navigator data to the LUT file +save([sPath, filesep, sName, '.mat'], 'dNavInt', '-append'); +save([sPath, filesep, sName, '.mat'], 'dNavInt_ms', '-append'); +save([sPath, filesep, sName, '.mat'], 'dSOSImg', '-append'); +% ------------------------------------------------------------------------- \ No newline at end of file diff --git a/reconstruction/fRetroGetNav2D.m b/reconstruction/fRetroGetNav2D.m new file mode 100644 index 0000000..c9ee372 --- /dev/null +++ b/reconstruction/fRetroGetNav2D.m @@ -0,0 +1,211 @@ +function [dNav, dKSpace] = fRetroGetNav2D(sFilename,lDebug) +%FFLASHRETROGETNAV2D Extracts navigator data from CS_FLASH_retro sequence +% DNAV = FFLASHRETROGETNAV(SFILENAME) Extracts navigator data DNAV from +% the siemens meas-data file SFILENAME. Must have been parsed with +% FMeasCreateLUT before. +% +% See also: FMEASCREATELUT +% +% Copyright 2014-2016 Christian Wuerslin, University of Tuebingen, Germany +% christian.wuerslin@med.uni-tuebingen.de +% and Thomas Kuestner, University of Tuebingen, Germany +% thomas.kuestner@med.uni-tuebingen.de + +if(nargin < 1) + lDebug = false; +end + +[sPath, sName, sExt] = fileparts(sFilename); + +% ------------------------------------------------------------------------- +% Get some relevant data from the drecksMDH +try + load(fullfile(sPath,[sName,'.mat']), 'SDrecksMDH'); + dBaseRes = SDrecksMDH.Geo.MatrixSize(1); + dTR = SDrecksMDH.Contrast.TR./1000; % in ms + if(isfield(SDrecksMDH.Geo,'EchoPosition')) + dEchoLine = SDrecksMDH.Geo.EchoPosition(1); + dEchoPartition = SDrecksMDH.Geo.EchoPosition(2); + else + dEchoLine = SDrecksMDH.Geo.MatrixSize(2)/2; + dEchoPartition = SDrecksMDH.Geo.MatrixSize(3)/2; + end + dNavPeriod = SDrecksMDH.Wip.NavPeriod; + dNavPERes = SDrecksMDH.Wip.NavRes(1); + +catch + % Load the loop counters + load([sPath, filesep, sName, '.mat']); + dBaseRes = 256; + dEchoLine = 128; + dEchoPartition = 36; + dNavPeriod = 200; + dNavPERes = 8; +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Read all projections (partition == echoPartition) AND (line == echoLine) +% Currently only uses the first echo, but could be extended to multi-echo +[dKSpace, iLC] = fMeasRead(sFilename, 'Set', [1 7], 'Seg', [0 60000], 'Eco', 0); +if(length(unique(iLC(:,10))) > 1) % in the case of switched on surrogate monitoring + lMask = true(size(iLC, 1), 1); + lMask = lMask & (iLC(:,10) == 1 | iLC(:,10) == 3 | iLC(:,10) == 5 | iLC(:,10) == 7); % binary encoded [ECG, Belt, Navi] -> cut out all cases in which Navi is on + iLC = iLC(lMask,:); + dKSpace = dKSpace(lMask,:); +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Determine data size +iNChannels = double(iLC(1, 2)); +iNSamples = size(dKSpace, 2); +dNavPERes = length(unique(iLC(:,3))); % correct number of phase-encoding lines +iNMeasurements = size(dKSpace, 1)./(iNChannels.*dNavPERes); +if(mod(iNMeasurements,1) ~= 0) % happens if a navi at the end was not completely acquired + dUniqueNavi = unique(iLC(:,3)); + iLastpos = find(iLC(:,3) == max(dUniqueNavi),1,'last'); + dKSpace = dKSpace(1:iLastpos,:); + iLC = iLC(1:iLastpos,:); + iNMeasurements = size(dKSpace, 1)./(iNChannels.*dNavPERes); +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Build up the kSpaces of all channels +dKSpace = reshape(dKSpace, [iNChannels, dNavPERes, iNMeasurements, iNSamples]); +dKSpace = permute(dKSpace, [4, 3, 2, 1]); % -> RO x t x PE x CH +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Reconstruct the 1-D projections for all measurements and all channels +dImg = fftshift(ifft(ifftshift(dKSpace))); +dImg = fftshift(ifft(ifftshift(dImg, 3), [], 3), 3); +dImg = flipdim(dImg, 1); % Invert the RO direction: 1-N -> H-F +dImg = dImg(iNSamples/4:iNSamples.*3/4 - 1, :, :, :); % RO x t x PE x CH + +dPower = sum(sum(sum(abs(dImg), 1), 2), 3); +dImg = dImg./repmat(dPower, [size(dImg, 1), size(dImg, 2), size(dImg, 3), 1]); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Get x range of respiratory motion +dFreq = fft(dImg, [], 2); +dPower = dFreq.*conj(dFreq); +dIMGres = 1./(double(dNavPeriod)./1000.*double(iNMeasurements)); % The frequency resolution of dIMG in Hz +dPower = squeeze(sum(dPower(:, round(1./(5.*dIMGres)):round(1./(3.*dIMGres)), :, :), 2)); % RO x PE x CH +if(ismatrix(dPower)) % happens for BH acquisition, i.e iNPHASES=1 + dPower = permute(dPower,[1 3 2]); +end +dPowerInChan = squeeze(sum(dPower, 2)); % RO x CH +dPowerAcrossChannels = sum(dPowerInChan, 2); % RO x 1 +dPowerAcrossChannels(length(dPowerAcrossChannels).*3./4:end) = 0; % Prevent detection of regions in the abdomen +dPowerAcrossChannels = conv(dPowerAcrossChannels, fGaussianLP(20), 'same'); +[dMax, dX] = max(dPowerAcrossChannels); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Sort out Channels with no relevant information in target area +dMax = max(dPowerInChan(dX - 20: dX + 20, :)); +lGoodChannels = dMax > 0.2.*max(dMax); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Get the best PE line +dPowerInPE = squeeze(sum(sum(dPower(dX - 20:dX + 20, :, lGoodChannels), 3))); +[dVal, dPos] = max(dPowerInPE); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% find best corresponding channels according to best phase encoding +% position +dFreq = fft(dImg(:,:,dPos,:),[],2); +dPower = dFreq.*conj(dFreq); +dIMGres = 1./(double(dNavPeriod)./1000.*double(iNMeasurements)); % The frequency resolution of dIMG in Hz +dPowerInChan = squeeze(sum(dPower(:, round(1./(5.*dIMGres)):round(1./(3.*dIMGres)), :), 2)); % RO x CH +dPowerAcrossChannels = sum(dPowerInChan, 2); % RO x 1 +dPowerAcrossChannels(length(dPowerAcrossChannels).*3./4:end) = 0; % Prevent detection of regions in the abdomen +dPowerAcrossChannels = conv(dPowerAcrossChannels, fGaussianLP(20), 'same'); +[dMax, dX] = max(dPowerAcrossChannels); +dMax = max(dPowerInChan(dX - 20: dX + 20, :)); +lGoodChannels = dMax > 0.2.*max(dMax); +if(lDebug) + [lMode, iChaSel, dPos] = fRetroNaviChoosePos( dImg, find(lGoodChannels), dPos, 1 ); + lGoodChannels = false(size(lGoodChannels)); + lGoodChannels(iChaSel) = true; +end +dRelevantImg = squeeze(dImg(:, :, dPos, lGoodChannels)); % In-plane res: 2 mm/px +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Sort out channels, which have distortions in the relevant area (bending +% artifact) +dRelevantImg = dRelevantImg.*conj(dRelevantImg); +dMax = max(max(dRelevantImg)); +dMax = repmat(dMax, [size(dRelevantImg, 1), size(dRelevantImg, 2), 1]); +dSOSImg = sqrt(sum(dRelevantImg./dMax, 3)); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Get the navigator signal +iDisplacement = 80; % [mm] diaphragm displacement (max. +/- 80mm) +dDisplacement = round(iDisplacement/(SDrecksMDH.Geo.FOV(1)/SDrecksMDH.Geo.MatrixSize(1))); % [px] + +dRefImg = zeros(size(dSOSImg)); +dRefImg(:,end) = dSOSImg(:,end); +idx = 1:size(dRefImg,2); +dNav = zeros(1,size(dRefImg,2)); + +hF = waitbar(0, 'Getting Navigator'); +for i=size(dRefImg,2)-1:-1:2 + dRMSImg = zeros(size(dRefImg,1),2*dDisplacement+1); + tmp = dRefImg(:,idx(i+1:end)); + for iD = -dDisplacement:dDisplacement + dRMSImg(:,dDisplacement-iD+1) = sum((tmp - repmat(circshift(dSOSImg(:,idx(i)), iD),[1 size(tmp,2)])).^2,2); % RMS + waitbar((abs(i-size(dRefImg,2)-2)*(2*dDisplacement+1) + (iD + dDisplacement + 1))/((size(dRefImg,2)-2)*(2*dDisplacement+1)),hF); + end + [~, dNav(i)] = min(sum(dRMSImg(dX-round(dDisplacement/2):dX+round(dDisplacement/2),:))); % RMS + %[~, dNav(i)] = min(sum(dRMSImg)); % RMS + dNav(i) = dDisplacement + 1 - dNav(i); + dRefImg(:,idx(i)) = circshift(dSOSImg(:,idx(i)), dNav(i)); + waitbar(abs(i-size(dRefImg,2)-2)/(size(dRefImg,2)-2),hF); +end +close(hF); +dNav = -dNav; +dNav = conv(dNav, fGaussianLP(5), 'same'); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Visualize the result +if(lDebug) + figure('Name','Navi'), imagesc(dSOSImg(:, 3:end)); + colormap gray; + hold all; + plot(dNav(3:end)-mean(dNav) + dX); +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Interpolate the navigator data to TR intervals +load([sPath, filesep, sName, '.mat'], 'iLC'); +iLC = iLC(1:iNChannels:end, :); +iLC = iLC(iLC(:,7) == 0, :); % first echo +iLC = iLC(iLC(:, 1) == 2*dBaseRes, :); % only "real" readouts + +dNav = -dNav; +dNav = dNav - min(dNav); + +iNavInd = find((iLC(:,3) == dEchoLine) & (iLC(:,6) == dEchoPartition)); +if(length(iNavInd) ~= length(dNav)) % happens at the end when navi not fully acquired + iNavInd = iNavInd(1:length(dNav)); +end +dNavInt = interp1(iNavInd, dNav, 1:length(iLC), 'pchip'); +dNavInt_ms = interp1(iNavInd.*dTR, dNav, 1:length(iLC).*dTR, 'pchip'); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Append the navigator data to the LUT file +save([sPath, filesep, sName, '.mat'], 'dNavInt', '-append'); +save([sPath, filesep, sName, '.mat'], 'dNavInt_ms', '-append'); +save([sPath, filesep, sName, '.mat'], 'dSOSImg', '-append'); +% ------------------------------------------------------------------------- \ No newline at end of file diff --git a/reconstruction/fRetroGetNav3D.m b/reconstruction/fRetroGetNav3D.m new file mode 100644 index 0000000..f8322e4 --- /dev/null +++ b/reconstruction/fRetroGetNav3D.m @@ -0,0 +1,241 @@ +function [dNav, dKSpace] = fRetroGetNav3D(sFilename,lDebug) +%FFLASHRETROGETNAV3D Extracts navigator data from CS_FLASH_retro sequence +% DNAV = FFLASHRETROGETNAV(SFILENAME) Extracts navigator data DNAV from +% the siemens meas-data file SFILENAME. Must have been parsed with +% FMeasCreateLUT before. +% +% See also: FMEASCREATELUT +% +% Copyright 2014-2016 Christian Wuerslin, University of Tuebingen, Germany +% christian.wuerslin@med.uni-tuebingen.de +% and Thomas Kuestner, University of Tuebingen, Germany +% thomas.kuestner@med.uni-tuebingen.de + +if(nargin < 1) + lDebug = false; +end + +[sPath, sName, sExt] = fileparts(sFilename); + +% ------------------------------------------------------------------------- +% Get some relevant data from the drecksMDH +try + load(fullfile(sPath,[sName,'.mat']), 'SDrecksMDH'); + dBaseRes = SDrecksMDH.Geo.MatrixSize(1); + dTR = SDrecksMDH.Contrast.TR./1000; % in ms + if(isfield(SDrecksMDH.Geo,'EchoPosition')) + dEchoLine = SDrecksMDH.Geo.EchoPosition(1); + dEchoPartition = SDrecksMDH.Geo.EchoPosition(2); + else + dEchoLine = SDrecksMDH.Geo.MatrixSize(2)/2; + dEchoPartition = SDrecksMDH.Geo.MatrixSize(3)/2; + end + dNavPeriod = SDrecksMDH.Wip.NavPeriod; + dNavPERes = SDrecksMDH.Wip.NavRes(1); + dNavSLRes = SDrecksMDH.Wip.NavRes(2); +catch + % Load the loop counters + load([sPath, filesep, sName, '.mat']); + dBaseRes = 256; + dEchoLine = 128; + dEchoPartition = 36; + dNavPeriod = 200; + dNavPERes = 8; + dNavSLRes = 8; +end +% ------------------------------------------------------------------------- + + +% ------------------------------------------------------------------------- +% Read all projections (partition == echoPartition) AND (line == echoLine) +% Currently only uses the first echo, but could be extended to multi-echo +[dKSpace, iLC] = fMeasRead(sFilename, 'Set', [1 7], 'Eco', 0); +if(length(unique(iLC(:,10))) > 1) % in the case of switched on surrogate monitoring + lMask = true(size(iLC, 1), 1); + lMask = lMask & (iLC(:,10) == 1 | iLC(:,10) == 3 | iLC(:,10) == 5 | iLC(:,10) == 7); + iLC = iLC(lMask,:); + dKSpace = dKSpace(lMask,:); +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Determine data size +iNChannels = double(iLC(1, 2)); +iNSamples = size(dKSpace, 2); +dNavPERes = length(unique(iLC(:,3))); % correct number of phase-encoding lines +dNavSLRes = length(unique(iLC(:,6))); +iNMeasurements = size(dKSpace, 1)./(iNChannels.*dNavPERes.*dNavSLRes); +if(mod(iNMeasurements,1) ~= 0) % happens if navi was acquired at the end and was not completed due to end of scan time + % crop away this remaining end + iPhaU = unique(iLC(:,3)); + iSlU = unique(iLC(:,6)); + iLastNav = find(iLC(:,3) == iPhaU(1) & iLC(:,6) == iSlU(1),1,'last'); + iLC = iLC(1:iLastNav-iNChannels,:); + dKSpace = dKSpace(1:iLastNav-iNChannels,:); + iNMeasurements = size(dKSpace, 1)./(iNChannels.*dNavPERes.*dNavSLRes); +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Build up the kSpaces of all channels +dKSpace = reshape(dKSpace, [iNChannels, dNavPERes, dNavSLRes, iNMeasurements, iNSamples]); +dKSpace = permute(dKSpace, [5, 4, 3, 2, 1]); % -> RO x t x SL x PE x CH +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Reconstruct the 1-D projections for all measurements and all channels +dImg = fftshift(ifft(ifftshift(dKSpace))); +dImg = fftshift(ifft(ifftshift(dImg, 3), [], 3), 3); +dImg = fftshift(ifft(ifftshift(dImg, 4), [], 4), 4); +dImg = flipdim(dImg, 1); % Invert the RO direction: 1-N -> H-F +dImg = dImg(iNSamples/4:iNSamples.*3/4 - 1, :, :, :, :); % RO x t x SL x PE x CH + +dPower = sum(sum(sum(sum(abs(dImg), 1), 2), 3), 4); +dImg = dImg./repmat(dPower, [size(dImg, 1), size(dImg, 2), size(dImg, 3), size(dImg, 4), 1]); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Get x range of respiratory motion +dFreq = fft(dImg, [], 2); +dPower = dFreq.*conj(dFreq); +dIMGres = 1./(double(dNavPeriod)./1000.*double(iNMeasurements)); % The frequency resolution of dIMG in Hz +dPower = squeeze(sum(dPower(:, round(1./(5.*dIMGres)):round(1./(3.*dIMGres)), :, :, :), 2)); % RO x SL x PE x CH +if(ismatrix(dPower)) % happens for BH acquisition, i.e iNPHASES=1 + dPower = permute(dPower,[1 3 2]); +end +dPowerInChan = squeeze(sum(sum(dPower, 2),3)); % RO x CH +dPowerAcrossChannels = sum(dPowerInChan, 2); % RO x 1 +dPowerAcrossChannels(length(dPowerAcrossChannels).*3./4:end) = 0; % Prevent detection of regions in the abdomen +dPowerAcrossChannels = conv(dPowerAcrossChannels, fGaussianLP(20), 'same'); +[dMax, dX] = max(dPowerAcrossChannels); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Sort out Channels with no relevant information in target area +dMax = max(dPowerInChan(dX - 20: dX + 20, :)); +lGoodChannels = dMax > 0.2.*max(dMax); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Get the best PE line +dPowerInPE = squeeze(sum(sum(sum(dPower(dX - 20:dX + 20, :, :, lGoodChannels), 2),4))); +[dVal, dPosPE] = max(dPowerInPE); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Get the best SL line +dPowerInSL = squeeze(sum(sum(dPower(dX - 20:dX + 20, :, dPosPE, lGoodChannels),4))); +[dVal, dPosSL] = max(dPowerInSL); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% find best corresponding channels according to best phase encoding +% position +dFreq = fft(dImg(:,:,dPosSL,dPosPE,:),[],2); +dPower = dFreq.*conj(dFreq); +dIMGres = 1./(double(dNavPeriod)./1000.*double(iNMeasurements)); % The frequency resolution of dIMG in Hz +dPowerInChan = squeeze(sum(dPower(:, round(1./(5.*dIMGres)):round(1./(3.*dIMGres)), :), 2)); % RO x CH +dPowerAcrossChannels = sum(dPowerInChan, 2); % RO x 1 +dPowerAcrossChannels(length(dPowerAcrossChannels).*3./4:end) = 0; % Prevent detection of regions in the abdomen +dPowerAcrossChannels = conv(dPowerAcrossChannels, fGaussianLP(20), 'same'); +[dMax, dX] = max(dPowerAcrossChannels); +dMax = max(dPowerInChan(dX - 20: dX + 20, :)); +lGoodChannels = dMax > 0.2.*max(dMax); +if(lDebug) + [lMode, iChaSel, dPosPE, dPosSL] = fRetroNaviChoosePos( dImg, find(lGoodChannels), dPosPE, dPosSL ); + lGoodChannels = false(size(lGoodChannels)); + lGoodChannels(iChaSel) = true; +end +dRelevantImg = squeeze(dImg(:, :, dPosSL, dPosPE, lGoodChannels)); % In-plane res: 2 mm/px +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Sort out channels, which have distortions in the relevant area (bending +% artifact) +dRelevantImg = dRelevantImg.*conj(dRelevantImg); +dMax = max(max(dRelevantImg)); +dMax = repmat(dMax, [size(dRelevantImg, 1), size(dRelevantImg, 2), 1]); +dSOSImg = sqrt(sum(dRelevantImg./dMax, 3)); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Get the navigator signal +% dNavImg1 = imfilter(dNavImg, [1 1]./2, 'replicate'); +% dFilter = fspecial('log', [11, 1]); +% for iI = 1:size(dNavImg, 2) +% dNavImgHP(:,iI) = conv(dSOSImg(:,iI), dFilter, 'same'); +% end +% dNavImgInt = interp1(dNavImgHP, linspace(1, size(dNavImgHP, 1), 256), 'pchip'); +% lMaxima = [zeros(1, size(dNavImgInt, 2)); diff(dNavImgInt,[], 1)] > 0 & [diff(dNavImgInt, [], 1); zeros(1, size(dNavImgInt, 2))] < 0; +% [dMax, dNav] = max(dNavImgInt); +% dMax = repmat(dMax, [size(dNavImgInt, 1), 1]); +% lMaxima(dNavImgInt < 0.5*dMax) = false; +% dNav = zeros(size(dNavImgInt, 2), 1); +% for iI = 1:length(dNav), dNav(iI) = find(lMaxima(:, iI), 1, 'first'); end +% dNav = dNav./256.*size(dNavImgHP, 1); +% dNav = conv(dNav, fGaussianLP(5), 'same'); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Get the navigator signal +iDisplacement = 80; % [mm] diaphragm displacement (max. +/- 80mm) +dDisplacement = round(iDisplacement/(SDrecksMDH.Geo.FOV(1)/SDrecksMDH.Geo.MatrixSize(1))); % [px] + +dRefImg = zeros(size(dSOSImg)); +dRefImg(:,end) = dSOSImg(:,end); +idx = 1:size(dRefImg,2); +dNav = zeros(1,size(dRefImg,2)); + +hF = waitbar(0, 'Getting Navigator'); +for i=size(dRefImg,2)-1:-1:2 + dRMSImg = zeros(size(dRefImg,1),2*dDisplacement+1); + tmp = dRefImg(:,idx(i+1:end)); + for iD = -dDisplacement:dDisplacement + dRMSImg(:,dDisplacement-iD+1) = sum((tmp - repmat(circshift(dSOSImg(:,idx(i)), iD),[1 size(tmp,2)])).^2,2); % RMS + waitbar((abs(i-size(dRefImg,2)-2)*(2*dDisplacement+1) + (iD + dDisplacement + 1))/((size(dRefImg,2)-2)*(2*dDisplacement+1)),hF); + end + [~, dNav(i)] = min(sum(dRMSImg(dX-round(dDisplacement/2):dX+round(dDisplacement/2),:))); % RMS + %[~, dNav(i)] = min(sum(dRMSImg)); % RMS + dNav(i) = dDisplacement + 1 - dNav(i); + dRefImg(:,idx(i)) = circshift(dSOSImg(:,idx(i)), dNav(i)); + waitbar(abs(i-size(dRefImg,2)-2)/(size(dRefImg,2)-2),hF); +end +close(hF); +dNav = -dNav; +dNav = conv(dNav, fGaussianLP(5), 'same'); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Visualize the result +if(lDebug) + figure('Name','Navi'), imagesc(dSOSImg(:, 3:end)); + colormap gray; + hold all; + plot(dNav(3:end)-mean(dNav) + dX); +end +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Interpolate the navigator data to TR intervals +load([sPath, filesep, sName, '.mat'], 'iLC'); +iLC = iLC(1:iNChannels:end, :); +iLC = iLC(iLC(:,7) == 0, :); % first echo +iLC = iLC(iLC(:, 1) == 2*dBaseRes, :); % only "real" readouts + +dNav = -dNav; +dNav = dNav - min(dNav); + +iNavInd = find((iLC(:,3) == dEchoLine) & (iLC(:,6) == dEchoPartition)); +if(length(iNavInd) ~= length(dNav)) % happens at the end when navi not fully acquired + iNavInd = iNavInd(1:length(dNav)); +end +dNavInt = interp1(iNavInd, dNav, 1:length(iLC), 'pchip'); +dNavInt_ms = interp1(iNavInd.*dTR, dNav, 1:length(iLC).*dTR, 'pchip'); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Append the navigator data to the LUT file +save([sPath, filesep, sName, '.mat'], 'dNavInt', '-append'); +save([sPath, filesep, sName, '.mat'], 'dNavInt_ms', '-append'); +save([sPath, filesep, sName, '.mat'], 'dSOSImg', '-append'); +% ------------------------------------------------------------------------- \ No newline at end of file diff --git a/reconstruction/fRetroRecon4D_main.m b/reconstruction/fRetroRecon4D_main.m new file mode 100644 index 0000000..eb6e2a5 --- /dev/null +++ b/reconstruction/fRetroRecon4D_main.m @@ -0,0 +1,200 @@ +function dImg = fRetroRecon4D_main(sFilename, iNPhases, dTolerance, dTime, sGating, sGatingGather, lDebug) +% +% main function for CS Retro reconstruction +% +% Copyright 2014-2016 Christian Wuerslin, University of Tuebingen, Germany +% christian.wuerslin@med.uni-tuebingen.de +% and Thomas Kuestner, University of Tuebingen, Germany +% thomas.kuestner@med.uni-tuebingen.de + +% inputs +% sFilename string of TWIX measurement file +% iNPhases number of gates +% dTolerance view sharing blending factor(s) b of +% respiratory gates: 1 (non-overlapping) <= b <= 2 (completely overlapping) +% dTime scalar/vector/array of crop scan time(s) (0 := no cropping | scan time [s]) +% sGating string to select gating procedure +% 'percentile' select gate centroids according to 5% and 95% percentile +% 'kmeans' perform k-means clustering to select gate centroids +% sGatingGather string of sample gathering mode: +% 'closest' take closest sample to gate centroid +% 'collect' gather all samples and let the +% reconstruction decide for the most +% favourable one; not working for cardiac +% 'none' no k-Space population, show simulations of +% the k-Space sampling density +% +% output +% dImg reconstructed image + +% ------------------------------------------------------------------------- +% parse inputs +if nargin < 7, lDebug = false; end +if nargin < 6, sGatingGather = 'closest'; end +if nargin < 5, sGating = 'percentile'; end +if nargin < 4, dTime = 0; end +if nargin < 3, dTolerance = 1; end +if nargin < 2, iNPhases = 4; end +if(nargin < 1 || isempty(sFilename)) + [sFile, sPath] = uigetfile({'*.dat', 'ADC measdata (*.dat)'} ,'Select ADC measdata file', 'MultiSelect', 'off', pwd); + if(~isempty(sFile)) + [~, sName] = fileparts(sFile); + else + return; + end + sFilename = [sPath,sFile]; +end + +if(~exist('CS_reconstruction','file')) + error('CS_reconstruction code not found! Please download it from "https://github.com/thomaskuestner/CS_LAB" or put it on Matlab path'); +end +currpath = fileparts(mfilename('fullpath')); +addpath(genpath([currpath,filesep,'utils'])); + +% calculation precision +sPrecision = 'single'; +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% load raw measurment file +[sPath, sName] = fileparts(sFilename); +[~, ~, ~, para.drecksMDH, para.iLCPositions] = fMeas_main([sPath, filesep, sName, '.dat'], false); +load([sPath, filesep, sName, '.mat'], 'iLC'); +para.iLC = iLC; +clear 'iLC'; +para.drecksMDH.Wip.NSamples = 0; +para.drecksMDH.Wip.Phases = iNPhases; +para.measPara.sequenceName = 'CS_Retro'; +para.measPara.dim = [para.drecksMDH.Geo.MatrixSize(2), 2*para.drecksMDH.Geo.MatrixSize(1), para.drecksMDH.Geo.MatrixSize(3), iNPhases, double(para.iLC(1,2))]; +para.measPara.LCall = ones(1,4); +dNavPERes = para.drecksMDH.Wip.NavRes(1); +dNavSLRes = para.drecksMDH.Wip.NavRes(2); +if(para.measPara.dim(3) == 1) + para.measPara.dimension = '2D'; +else + if(para.measPara.dim(4) == 1) + para.measPara.dimension = '3D'; + else + para.measPara.dimension = '4D'; + end +end +iNEchos = length(unique(para.iLC(:,7))); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% prepare output file +if(isscalar(dTime)) + dTShow = num2str(dTime,'%03d'); +else + dTShow = dTime.'; + dTShow = strrep(num2str(dTShow(:).','%03d'),' ',''); +end +sSaveFilename = [sPath, filesep, sName, sprintf('_Ph%1d_Tol%03.0f_t%s_%s_%s', iNPhases, 100*dTolerance, dTShow, sGating, sGatingGather)]; +% ------------------------------------------------------------------------- + +for iEcho = 1:iNEchos + if(iNEchos > 1) + sSaveFilename = [sSaveFilename,'_',num2str(iEcho,'E%02d')]; + end + [~, sFilenameShow, ~] = fileparts(sSaveFilename); + + + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + % extract self-navigation signal + if(~any(arrayfun(@(x) strcmp(x.name,'dNavInt'), whos('-file', [sPath, filesep, sName, '.mat'])))) + if(dNavPERes == 1 && dNavSLRes == 1) % 1D navigator signal + dNav = fRetroGetNav(sFilename,lDebug); + elseif(dNavPERes > 1 && dNavSLRes == 1) % 2D navigator signal + dNav = fRetroGetNav2D(sFilename,lDebug); + else % 3D navigator signal + dNav = fRetroGetNav3D(sFilename,lDebug); + end + end + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + % perform retrospective gating + if(~exist([sSaveFilename,'_kSpace.mat'],'file')) + kSpace = fRetroGating(sFilename, iNPhases, dTolerance, dTime, sGating, sGatingGather, sPrecision, iEcho, false); + save([sSaveFilename,'_kSpace.mat'], 'kSpace', '-v7.3'); + else + fprintf('Loading k-space: %s\n',sFilenameShow); + load([sSaveFilename,'_kSpace.mat']); + end + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + % Compressed Sensing reconstruction + if(strcmp(sGatingGather,'collect')) + kSpace = mat2cell(kSpace, size(kSpace, 1), size(kSpace, 2), size(kSpace, 3), size(kSpace, 4), size(kSpace,5), ones(1, size(kSpace, 6))); + kSpace = shiftdim(kSpace, 4); + else + kSpace = mat2cell(kSpace(:,:,:,:,:), size(kSpace, 1), size(kSpace, 2), size(kSpace, 3), size(kSpace, 4), ones(1, size(kSpace, 5))); + kSpace = shiftdim(kSpace, 3); + end + + nSize = size(kSpace{1}); + para = fPrepareRecon(para,nSize,sPath,sGatingGather); + [imgAll,imgCha] = CS_reconstruction(kSpace,para); + if(strcmp(sGatingGather,'collect') || ndims(imgAll) == 4) + dImg = imgAll; + else + dImg = zeros([size(imgAll,1), size(imgAll,2), nSize(3), nSize(4)]); + for i=1:iNPhases + dImg(:,:,:,i) = imgAll(:,:,(i-1)*nSize(3)+1:i*nSize(3)); + end + end + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + % save results + save([sSaveFilename,'_channel.mat'], 'imgCha', '-v7.3'); + save([sSaveFilename,'_recon.mat'], 'dImg', '-v7.3'); + % - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - +end +end + + +function para = fPrepareRecon(para,nSize,sPath,sGatingGather) +% set CS reconstruction parameters + nPha = nSize(1); + nFreq = nSize(2); + nZ = nSize(3); + nTime = nSize(4); + if(length(nSize) > 4) + nSmp = nSize(5); + else + nSmp = 1; + end + + % additional length of gating dimension + % drecksMDH v2.0 + para.drecksMDH.Wip.Phases = nTime; + if(strcmp(sGatingGather,'collect')) + para.drecksMDH.Wip.NSamples = nSmp; + end + + % CS recon parameter + if(strcmp(para.measPara.dimension,'3D')) + para.trafo.trafodim = [1 2 1 0; 1 0 1 0]; % 3D + elseif(strcmp(para.measPara.dimension,'4D')) + para.trafo.trafodim = [1 2 1 0; 0 0 0 1]; % 4D + end + para.cstype = 'FOCUSS'; + para.transformation = 'pca'; + para.trafo.scrambledim = [1 1 1 0]; + para.lambda = 0.01; + para.lambdaCalib = 0; % GRAPPA off + para.lambdaTV = 0; + para.lambdaESPReSSo = 0; + para.lambdaMC = 0; + para.sScaling = 'self'; + para.flagOversampling = logical([1 0 0]); % oversampling correction + para.prop.flagPlot = false; % Imagine off + para.prop.flagOptim = false; % no FFT optimization + para.savePath = sPath; + +end \ No newline at end of file diff --git a/reconstruction/utils/OptimalColor.m b/reconstruction/utils/OptimalColor.m new file mode 100644 index 0000000..2c1e469 --- /dev/null +++ b/reconstruction/utils/OptimalColor.m @@ -0,0 +1,37 @@ +function dColormap = OptimalColor(iNBins) +%OPTIMALCOLORS Example custom colormap for use with imagine +% DCOLORMAP = OPTIMALCOLOR(INBINS) returns a double colormap array of size +% (INBINS, 3). Use this template to implement you own custom colormaps. +% Imagine will interpret all m-files in this folder as potential colormap- +% generating functions an list them using the filename. + +% ------------------------------------------------------------------------- +% Process input +if ~nargin, iNBins = 256; end +iNBins = uint16(iNBins); +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Create look-up tables (pairs of x- and y-vectors) for the three colors +dYRed = [0; 0.471; 0.471; 0.518; 0.518; 1; 1]; +dXRed = [1; 31; 88; 100; 205; 254; 256]; + +dYGrn = [0; 0; 1; 1]; +dXGrn = [1; 30; 132; 256]; + +dYBlu = [0; 0; 1; 1]; +dXBlu = [1; 79; 181; 256]; +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% Interpolate and concatenate vectors to the final colormap +dRedInt = interp1(dXRed, dYRed, linspace(1, 255, iNBins)'); +dGrnInt = interp1(dXGrn, dYGrn, linspace(1, 255, iNBins)'); +dBluInt = interp1(dXBlu, dYBlu, linspace(1, 255, iNBins)'); + +dColormap = [dRedInt, dGrnInt, dBluInt]; +% ------------------------------------------------------------------------- + +% ========================================================================= +% *** END OF FUNCTION OptimalColor +% ========================================================================= \ No newline at end of file diff --git a/reconstruction/utils/fGaussianLP.m b/reconstruction/utils/fGaussianLP.m new file mode 100644 index 0000000..b24de5e --- /dev/null +++ b/reconstruction/utils/fGaussianLP.m @@ -0,0 +1,27 @@ +function dDataOut = fGaussianLP(dData, dW) +%GAUSSIANLP performs a 1D gaussian lowpass filter. +% +% HOUT = GAUSSIANLP(HIN, L, MODE) performs a gaussian lowpass of length L +% to the array in HIN. The gaussian lowpass is defined in the range +% [-floor(L/2), floor(L/2)] and the variance sigma of the Gaussinan is +% chosen such, that the definition boundaries of the Gaussian are equal to +% 3*sigma. MODE defines the method for padding the input array HIN and can +% be either 'circular', 'replicate' or 'symmetroc'. See also the +% documentation of 'padarray' for more information on the padding methods. + +if nargin == 1, dW = dData; end + +% ------------------------------------------------------------------------- +% Create Gaussian kernel. +dX = (-floor(dW/2):floor(dW/2))'; +dX = dX.*6./dW; +dG = exp(-0.5.*dX.^2); +dG = dG./sum(dG); +if nargin == 1 + dDataOut = dG; + return +end +% ------------------------------------------------------------------------- + +dDataOut = zeros(size(dData)); +for iI = 1:size(dData, 2), dDataOut(:,iI) = conv(dData(:,iI), dG, 'same'); end \ No newline at end of file