From 64d8ff25ccb81c043bf9e73ca236eca2d2537e7e Mon Sep 17 00:00:00 2001 From: Alexandros-Stavros Iliopoulos <1577182+ailiop@users.noreply.github.com> Date: Mon, 11 Feb 2019 18:21:52 -0500 Subject: [PATCH] bug fixes (local control window & 2D optimal control values); domain prepadding; 3D deformation demo script and more informative demo visualizations; documentation update and extension per JOSS review comments --- +dvf/changeUnit.m | 2 +- +dvf/coorDisplace.m | 2 +- +dvf/eigJacobian.m | 2 +- +dvf/feedbackControlVal.m | 45 ++- +dvf/icResidual.m | 8 +- +dvf/imdeform.m | 43 +-- +dvf/interpnVf.m | 2 +- +dvf/inversion.m | 264 +++++++++++------ +dvf/issingularJacobian.m | 4 +- +dvf/jacobian.m | 2 +- +dvf/maskDomain.m | 12 +- +dvf/ntdcMeasures.m | 17 +- +dvf/rescale.m | 2 +- +dvf/sizeVf.m | 2 +- +util/{imageGrid.m => imageGridSmooth.m} | 63 ++-- +util/ndgridGpu.m | 2 +- +util/parseOptArgs.m | 2 +- +vis/imageSlices.m | 165 +++++++++++ +vis/mfigure.m | 100 +++++++ CODE_OF_CONDUCT.md | 87 ++++++ CONTRIBUTING.md | 134 +++++++++ README.md | 305 +++++++++++++------ codemeta.json | 4 +- demo_inversion_2d.m | 171 ++++++++--- demo_inversion_3d.m | 245 ---------------- demo_inversion_3d_z0.m | 359 +++++++++++++++++++++++ demo_inversion_3d_zsin.m | 353 ++++++++++++++++++++++ paper.md | 6 +- references.bib | 4 +- 29 files changed, 1866 insertions(+), 541 deletions(-) rename +util/{imageGrid.m => imageGridSmooth.m} (63%) create mode 100644 +vis/imageSlices.m create mode 100644 +vis/mfigure.m create mode 100644 CODE_OF_CONDUCT.md create mode 100644 CONTRIBUTING.md delete mode 100644 demo_inversion_3d.m create mode 100644 demo_inversion_3d_z0.m create mode 100644 demo_inversion_3d_zsin.m diff --git a/+dvf/changeUnit.m b/+dvf/changeUnit.m index dc7d74c..b46f684 100644 --- a/+dvf/changeUnit.m +++ b/+dvf/changeUnit.m @@ -74,7 +74,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+dvf/coorDisplace.m b/+dvf/coorDisplace.m index 6dca98c..b05606f 100644 --- a/+dvf/coorDisplace.m +++ b/+dvf/coorDisplace.m @@ -124,7 +124,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+dvf/eigJacobian.m b/+dvf/eigJacobian.m index 60f786a..b02e9e0 100644 --- a/+dvf/eigJacobian.m +++ b/+dvf/eigJacobian.m @@ -102,7 +102,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+dvf/feedbackControlVal.m b/+dvf/feedbackControlVal.m index 6bd70c5..a60d387 100644 --- a/+dvf/feedbackControlVal.m +++ b/+dvf/feedbackControlVal.m @@ -116,9 +116,9 @@ % MU(x) = 1 - 2 / (Lmax(x) + Lmin(x)) [case R] % % where Lc and Lr refer to the complex and real eigenvalues, -% respectively, at some point. Case R means all eigenvalues at x are -% real, while case C means that two eigenvalues at x form a complex -% conjugate pair. +% respectively, at some point x. Case R means all eigenvalues at x +% are real, while case C means that two eigenvalues at x form a +% complex conjugate pair. % % * 'alternating' % @@ -147,11 +147,13 @@ % % [1] A. Dubey*, A.-S. Iliopoulos*, X. Sun, F.-F. Yin, and L. Ren, % "Iterative inversion of deformation vector fields with feedback -% control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, May 2018. -% DOI: 10.1002/mp.12962. +% control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, 2018. +% - DOI: 10.1002/mp.12962 +% - arXiv: 1610.08589 [cs.CV] % % [2] A. Dubey, "Symmetric completion of deformable registration via -% bi-residual inversion," PhD thesis, Duke University, Durham, NC, USA. +% bi-residual inversion," PhD thesis, Duke University, Durham, NC, USA, +% 2018. % % % See also dvf.jacobian, dvf.eigJacobian, dvf.inversion @@ -322,26 +324,34 @@ MaskCtrl = reshape( MaskCtrl, [prod(szDom), 1] ); Mu = zeros( [prod(szDom), 1], 'like', real(Lambda) ); + % real & complex eigenvalue masks (within controllable region) + MaskReal = (max( abs(imag(Lambda)), [], 2 ) < tauZero) & MaskCtrl; + MaskComplex = ~MaskReal & MaskCtrl; + % 2D/3D DVF switch dim case 2 % ======== 2D (real/conjugate pair) - Mu(MaskCtrl) = 1 - 2 ./ sum( real(Lambda(MaskCtrl,:)), 2 ); + % ----- REAL EIGENVALUES (CASE R) + + Mu(MaskReal) = 1 - 2 ./ sum( real(Lambda(MaskReal,:)), 2 ); + + % ----- COMPLEX EIGENVALUES (CASE C) + + % (complex eigenvalues always appear as conjugate pairs in 2D) + Mu(MaskComplex) = 1 - (real(Lambda(MaskComplex,1)) ./ ... + abs(Lambda(MaskComplex,1)).^2); case 3 % ======== 3D (real & complex lambda) % ----- REAL EIGENVALUES (CASE R) - MaskReal = (max( abs(imag(Lambda)), [], 2 ) < tauZero) & MaskCtrl; Mu(MaskReal) = 1 - 2 ./ (max( real(Lambda(MaskReal,:)), [], 2 ) + ... min( real(Lambda(MaskReal,:)), [], 2 )); % ----- COMPLEX EIGENVALUES (CASE C) - % complex-eigenvalues domain mask - MaskComplex = ~MaskReal & MaskCtrl; - % row-column indices of real and (upper halfspace) complex eigenvalues [~, jR] = min( abs(imag(Lambda(MaskComplex,:))), [], 2 );% real col-idx [~, jC] = max( imag(Lambda(MaskComplex,:)) , [], 2 );% complex col-idx @@ -388,8 +398,17 @@ % Abhishek Dubey abhisdub@cs.duke.edu % Xiaobai Sun xiaobai@cs.duke.edu % -% VERSION +% RELEASE +% +% 1.0.2 - February 11, 2019 +% +% CHANGELOG +% +% 1.0.2 (Feb 11, 2019) - Alexandros +% ! fixed bug with pointwise optimal control scheme ('optimal') +% whereby the complex-eigenvalues case was not handled for 2D DVFs % -% 1.0.0 - October 31, 2018 +% 1.0.0 (Oct 31, 2018) +% . initial release % % ------------------------------------------------------------ diff --git a/+dvf/icResidual.m b/+dvf/icResidual.m index 74f9f29..394a2ed 100644 --- a/+dvf/icResidual.m +++ b/+dvf/icResidual.m @@ -78,10 +78,12 @@ % [1] A. Dubey*, A.-S. Iliopoulos*, X. Sun, F.-F. Yin, and L. Ren, % "Iterative inversion of deformation vector fields with feedback % control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, 2018. -% [DOI:10.1002/mp.12962] +% - DOI: 10.1002/mp.12962 +% - arXiv: 1610.08589 [cs.CV] % % [2] A. Dubey, "Symmetric completion of deformable registration via -% bi-residual inversion," PhD thesis, Duke University, Durham, NC, USA. +% bi-residual inversion," PhD thesis, Duke University, Durham, NC, USA, +% 2018. % % [3] G. E. Christensen and H. J. Johnson, "Consistent image % registration," IEEE Transactions on Medical Imaging, vol. 20, no. 7, @@ -148,7 +150,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+dvf/imdeform.m b/+dvf/imdeform.m index 1a9b976..bce4136 100644 --- a/+dvf/imdeform.m +++ b/+dvf/imdeform.m @@ -1,4 +1,4 @@ -function Idfm = imdeform (Iref, DVF, mapping, interpMethod, extrapVal) +function IDfm = imdeform (IRef, DVF, mapping, interpMethod, extrapVal) % % IMDEFORM - 2D/3D image deformation % @@ -30,7 +30,7 @@ % % INPUT % -% IREF Reference (to-be-deformed) 3D image [NX x NY x NZ] +% IREF Reference (to-be-deformed) image [NX x NY x NZ] % DVF Deformation vector field [NX x NY x NZ x 3] % (displacement measured in pixels) % MAPPING Deformation mapping [string] @@ -45,9 +45,9 @@ % % OUTPUT % -% IDFM Deformed 3D image [NX x NY x NZ] +% IDFM Deformed image [NX x NY x NZ] % -% NOTES +% NOTE % % The image/vector field array sizes in the documentation refer to the % 3D case. In the 2D case, the sizes are amended accordingly: @@ -60,19 +60,20 @@ % IDFM = IMDEFORM(IREF,DVF,'bwd') deforms the input reference 3D image, % such that % -% IDFM(x) = IREF(x + DVF(x)) +% IDFM(x) = IREF(x + DVF(x)) , % -% for all x in the image domain. The above model is referred to as the -% backward-mapping deformation model. +% for all x in the deformed-image domain. The above model is referred to +% as the backward-mapping deformation model. % % IDFM = IMDEFORM(IREF,DVF,'fwd') applies the forward-mapping deformation % model, instead: % -% IDFM(x + DVF(x)) = IREF(x) . +% IDFM(x + DVF(x)) = IREF(x) , % -% This tends to be slower, as it is based on scattered data -% interpolation, rather than gridded data interpolation (which is the -% case for backward-mapping). +% for all x in the reference-image domain. This is much slower, +% especially for large images, as it is based on scattered data +% interpolation. Instead, the backward-mapping model is based on gridded +% data interpolation. % % IDFM = IMDEFORM(IREF,DVF,...,MINTERP) also specifies the interpolation % method to be used when regridding the deformed image onto the regular @@ -116,7 +117,7 @@ % spatial domain size and dimensionality [szDom, dim] = dvf.sizeVf( DVF ); - % ensure 2D/3D case + % ensure 2D/3D case (the code actually supports ND image deformation) if ~( ((dim == 3) && (length(szDom) == 3)) || ... ((dim == 2) && (length(szDom) == 2)) ) error( [mfilename ':InvalidDimensionality'], ... @@ -142,12 +143,12 @@ case {'bwd','backward'} % ----- BACKWARD MAPPING - Idfm = interpn( Iref, gridDfm{:}, interpMethod, NaN ); + IDfm = interpn( IRef, gridDfm{:}, interpMethod, NaN ); case {'fwd','forward'} % ----- FORWARD MAPPING gridDfm = cellfun( @(x) x(:), gridDfm, 'UniformOutput', false ); - Idfm = griddata( gridDfm{:}, Iref(:), gridRef{:}, interpMethod ); + IDfm = griddata( gridDfm{:}, IRef(:), gridRef{:}, interpMethod ); otherwise % ----- ERROR (UNKNOWN MAPPING) @@ -157,7 +158,7 @@ end % switch (deformation mapping) % set any out-of-domain (NaN) values to 0 - Idfm( isnan(Idfm) ) = extrapVal; + IDfm( isnan(IDfm) ) = extrapVal; end @@ -171,8 +172,16 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % -% 1.0.0 - October 31, 2018 +% 1.0.2 - Ferbuary 11, 2019 +% +% CHANGELOG +% +% 1.0.2 (Feb 11, 2019) - Alexandros +% . documentation clean-up +% +% 1.0.0 (Oct 31, 2018) - Alexandros +% . initial release % % ------------------------------------------------------------ diff --git a/+dvf/interpnVf.m b/+dvf/interpnVf.m index 9b613ce..88b87bf 100644 --- a/+dvf/interpnVf.m +++ b/+dvf/interpnVf.m @@ -113,7 +113,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+dvf/inversion.m b/+dvf/inversion.m index 4f15e9d..2038dd0 100644 --- a/+dvf/inversion.m +++ b/+dvf/inversion.m @@ -58,6 +58,9 @@ % * K is the total number of iteration steps % * PRCT_RG_K and PRCT_RF_K contain the 50th, 90th, 95th, 98th, % and 100th percentiles of IC residual magnitudes per step +% * Computation of the step-wise residual percentiles will +% increase the execution time of INVERSION (especially in the +% case of multi-scale iteration steps) % % OPTIONS % @@ -70,7 +73,7 @@ % control value scheme [1]: % % * 'midrange': local midrange values; -% * 'optimal': pointwise optimal values (sensitive to initial guess); +% * 'optimal': locally optimal values; % * 'alternating': alternating values (midrange percentiles); or % * 'global-midrange': global midrange value. % @@ -140,8 +143,7 @@ % stable, linear-rate convergence towards the inverse. % % Once the error has been made sufficiently small, ||E_k(x)||_oo <= 1, -% the iteration proceeds to the second phase to accelerate convergence -% [2]: +% the iteration proceeds to phase two, accelerating convergence [2]: % % G_k+1(x) = G_k(x) - RF_k(x + G_k(x)) , % @@ -196,10 +198,12 @@ % [1] A. Dubey*, A.-S. Iliopoulos*, X. Sun, F.-F. Yin, and L. Ren, % "Iterative inversion of deformation vector fields with feedback % control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, 2018. -% DOI: 10.1002/mp.12962. +% - DOI: 10.1002/mp.12962 +% - arXiv: 1610.08589 [cs.CV] % % [2] A. Dubey, "Symmetric completion of deformable registration via -% bi-residual inversion," PhD thesis, Duke University, Durham, NC, USA. +% bi-residual inversion," PhD thesis, Duke University, Durham, NC, USA, +% 2018. % % % See also dvf.feedbackControlVal, dvf.icResidual @@ -227,12 +231,13 @@ opt.Mask = dvf.maskDomain( F ); % static parameters - accAsync = true; % asunchronous acceleration stage? + accAsync = true; % asynchronous acceleration stage? interp = 'linear'; % interpolation method + flagPrepad = true; % prepad spatial domain? prcIcMeasure = [50 90 95 98 100]; % ICR percentile ranks icrFullResolution = true; % calculate ICR in full resolution? - windowControl = 0.05 * min(szDom); % max local-control window size - fallbackControl = 0.8; % mu value at incontrollable points + windowControl = []; % max local-neighborhood size + fallbackControl = 0.8; % mu value at uncontrollable points prcMidAlternating = [98 50]; % alternating control percentiles prcErrorBound = 99; % Jacobian norm percentile rank tauSingularJ = 1e-6; % Jacobian singularity threshold @@ -243,13 +248,13 @@ opt = util.parseOptArgs( opt, varargin{:} ); % optional computation flags - flag.prctR = (nargout > 3); + flagOut.prctR = (nargout > 3); % disable multi-scale iteration with more than 2 scales if length( opt.scale ) > 2 error( [mfilename ':TooManyScales'], ... - ['Multi-scale iteration with more than 2 scales is currently ' ... - 'disabled.'] ); + ['Multi-scale iteration with more than 2 scales' ... + ' is currently disabled.'] ); end @@ -280,11 +285,31 @@ %% PRE-PROCESSING - [Mu, bndNormJF, szDomScl, MaskCtrl] = ... - preprocessing( F, opt.control, opt.scale, opt.Mask, windowControl, ... + % differential and spectral pre-computations + % *** must precede pre-padding: the MATLAB GRADIENT function uses a + % single-sided difference model at boundary points (vs central + % differences at interior points), hence pre-padding would affect + % boundary values, which in turn could influence the (quasi-)global + % bound on the deformation Jacobian norm which controls transition + % to the implicit-Newton acceleration iteration + [Mu, bndNormJF, MaskCtrl] = ... + preprocessing( F, opt.control, opt.Mask, windowControl, ... fallbackControl, prcErrorBound, prcMidAlternating, ... opt.numIteration, tauSingularJ, flagDetJ, dim ); + % pre-pad spatial domain (to avoid missing boundary values when down- and + % up-sampling with IMRESIZE in the case of multi-scale computation) + if flagPrepad + [F, G, opt.Mask, MaskCtrl, Mu, npad] = ... + pad( F, G, opt.Mask, MaskCtrl, Mu, min(opt.scale) ); + szDom = szDom + 2*npad; % (padded domain size at full resolution) + end + + % turn scale factors to scaled-domain sizes to avoid floor/ceiling issues + % when scaling down and then up with IMRESIZE (e.g. 11->6->12) + szDomScl = arrayfun( @(s) ceil( szDom * s ), opt.scale, ... + 'UniformOutput', false ); + %% DVF INVERSION ITERATION @@ -309,10 +334,10 @@ opt.stopThreshold(s), opt.accThreshold, ... bndNormJF, opt.numIteration(s), interp, ... valExtrap, opt.scale(s)*windowControl, ... - fallbackControl, prcIcMeasure, szDom, accAsync, ... - flag.prctR, icrFullResolution ); + fallbackControl, prcIcMeasure, szDom, ... + accAsync, flagOut.prctR, icrFullResolution ); - % scale back to original domain size + % scale back to original (padded) domain size G = dvf.rescale( G, szDom ); % update "total" iteration counter @@ -327,13 +352,18 @@ ks = ks + 1; % final IC residual measures - if flag.prctR + if flagOut.prctR [~, ~, ~, ~, prctRGIter(:,ks), prctRFIter(:,ks)] = ... icResidualCalc( F, G, interp, valExtrap, {}, opt.Mask, ... - prcIcMeasure, flag.prctR ); + prcIcMeasure, flagOut.prctR ); end % if (final IC residual measures) - + % de-pad spatial domain + if flagPrepad + [G, MaskCtrl, Mu] = depad( G, MaskCtrl, Mu, npad ); + end + + end @@ -341,13 +371,12 @@ %% ================================================== % PRE-PROCESSING -function [Mu, bndJFNorm, szDomScl, MaskCtrl, JF, LambdaF] = ... - preprocessing (F, control, scales, Mask, wndCtrl, valCtrlFallback, ... +function [Mu, bndJFNorm, MaskCtrl] = ... + preprocessing (F, control, Mask, wndCtrl, valCtrlFallback, ... prcErrBnd, prcMidAlt, numIter, tauJ, flagDetJ, dim) % IN F forward DVF [NX x NY x NZ x 3] % control feedback control param. [ | string] % tauAcc acceleration threshold [scalar] -% scales iteration resol. scales [S-vector] % Mask domain of interest mask [NX x NY x NZ] (logical) % wndCtrl feedback control n/hood [scalar | WX x WY x WZ] % valCtrlFallback control fallback value [scalar] @@ -359,11 +388,9 @@ % dim domain dimensionality [2 | 3] % OUT Mu control parameter values[K-vector | NX x NY x NZ] % bndJFNorm ||JF(x)||_oo bound [scalar] -% szDomScl scaled domain sizes [S-cell(3-vector)] % MaskCtrl controllable domain mask[NX x NY x NZ] (logical) -% JF deformation Jacobian [NX x NY x NZ x 3 x 3] -% LambdaF Jacobian eigenvalues [NX x NY x NZ x 3] - +% npad + % forward deformation Jacobian JF = dvf.jacobian( F ); @@ -381,7 +408,7 @@ % local control neighborhood (if adaptive control neighborhoods are used, % set the trivial neighborhood---local calculations will be done % during the iteration rather than preprocessing) - if isscalar( wndCtrl ) + if isscalar( wndCtrl ) || isempty( wndCtrl ) wndCtrl = 1; end @@ -402,34 +429,28 @@ end % if (adaptive/prefixed control?) - % non-invertible & incontrollable domain warning - if any( MaskNonInv, 'all' ) + % non-invertible & uncontrollable domain warning + if any( MaskNonInv & Mask, 'all' ) warning( [mfilename ':NonInvertibleRegion'], ... 'Input DVF is non-invertible at %d/%d points', ... nnz(MaskNonInv), numel(MaskNonInv) ); end - if any( ~MaskCtrl, 'all' ) - warning( [mfilename ':IncontrollableRegion'], ... - 'The inversion error is incontrollable at %d/%d points', ... + if any( ~MaskCtrl & Mask, 'all' ) + warning( [mfilename ':UncontrollableRegion'], ... + 'The inversion error is uncontrollable at %d/%d points', ... nnz(~MaskCtrl), numel(MaskCtrl) ); end MaskCtrl = MaskCtrl & ~MaskNonInv; % in case of constant/alternating control values, replicate for convenience if isvector( Mu ) - Mu = repmat( Mu(:), [ceil( max(numIter) * length(Mu) ), 1] ); + Mu = repmat( Mu(:), [ceil( max(numIter) / length(Mu) ), 1] ); end % "soft" upper bound of Jacobian oo-norm over controllable/masked domain bndJFNorm = max( sum( abs(JF), dim+2 ), [], dim+1 ); bndJFNorm = prctile( bndJFNorm(MaskCtrl), prcErrBnd ); - % turn scale factors to scaled-domain sizes to avoid floor/ceiling issues - % when scaling down and then up with IMRESIZE (e.g. 11->6->12) - [szDom, ~] = dvf.sizeVf( F ); - szDomScl = arrayfun( @(s) ceil( szDom * s ), scales, ... - 'UniformOutput', false ); - end @@ -519,7 +540,7 @@ end % local neighborhood radius - if isscalar( radiusMuMax ) + if isscalar( radiusMuMax ) || isempty( radiusMuMax ) radiusMu = min( [max( RFNorm(MaskCtrl) ), ... max( RGNorm(MaskCtrl) ) * bndNormJF, ... radiusMuMax] ); @@ -644,11 +665,12 @@ % IC residual: RF(y) = F(y) + G(y + F(y)) [y = x + G(x)] % = (F(x+G(x)) + G(x)) + G(x+G(x)+F(x+G(x))) - G(x) % = RG(x) + G(x + RG(x)) - G(x) - % *** RF(x + G(x)) can be numerically unstable (in part because it - % involves non-small deformations when RF is used as feedback in - % the implicit Newton iteration), plus it is inefficient compared - % to the formulation below (involving an additional interpolation - % step) + % *** direct calculation of RF(x + G(x)) is relatively unstable + % numerically (in part because it involves non-small deformations + % when RF is used as feedback in the implicit Newton iteration); it + % is also inefficient and very resolution-limited, involving an + % additional interpolation step compared to the equivalent + % formulation in terms of RG(x) gridRG = cell( size(gridG) ); [gridRG{:}] = dvf.coorDisplace( RG ); RF = RG + (dvf.interpnVf( G, gridRG, interp, extrap ) - G); @@ -669,6 +691,55 @@ +%% ================================================== +% IC RESIDUAL MAGNITUDE PERCENTILES + +function prctR = icrMgnPercentiles (R, Mask, prcPop, flag) +% IN R IC residual [NX x NY x NZ x 3] +% Mask spatial domain mask [NX x NY x NZ] (logical) +% prcPop ICR percentile ranks [P-vector] +% flag compute ICR percentiles?[boolean] +% OUT prctR ICR magnitude prctiles [P-vector] + + if flag % ---- compute ICR percentiles + dim = ndims( Mask ); + RMgn = sqrt( sum( R.^2, dim+1 ) ); + [~, RMgn] = sanitizeIcResidual( R, RMgn ); + prctR = prctile( RMgn(Mask), prcPop, 'all' ); + else % ---- do not compute anything + prctR = NaN; + end % if + +end + + + +%% ================================================== +% SANITIZE NAN/INF INTERPOLATION VALUES IN IC RESIDUAL + +function [R, RMgn] = sanitizeIcResidual (R, RMgn) +% R IC residual vectorfield [NX x NY x NZ x 3] +% RMgn IC residual magnitudes [NX x NY x NZ] + + % domain dimensionality + dim = ndims( RMgn ); + + % spatial mask of "bad" interpolation values (NaN/Inf) + MaskBad = isnan( RMgn ) | isinf( RMgn ); + + % set invalid IC residual magnitudes equal to max IC residual (in + % order to not skew statistics) + RMgn(MaskBad) = max( RMgn(~MaskBad) ); + + % set invalid IC residual vectors equal to 0 (no feedback control) + R = reshape( R, [numel(MaskBad), dim] ); + R(MaskBad,:) = 0; + R = reshape( R, [size(MaskBad), dim] ); + +end + + + %% ================================================== % FIELD & MASK RESCALING PRIOR TO ITERATION CORE @@ -762,49 +833,68 @@ %% ================================================== -% SANITIZE NAN/INF INTERPOLATION VALUES IN IC RESIDUAL +% PAD SPATIAL DOMAIN -function [R, RMgn] = sanitizeIcResidual (R, RMgn) -% R IC residual vectorfield [NX x NY x NZ x 3] -% RMgn IC residual magnitudes [NX x NY x NZ] - - % domain dimensionality - dim = ndims( RMgn ); - - % spatial mask of "bad" interpolation values (NaN/Inf) - MaskBad = isnan( RMgn ) | isinf( RMgn ); - - % set invalid IC residual magnitudes equal to max IC residual (in - % order to not skew statistics) - RMgn(MaskBad) = max( RMgn(~MaskBad) ); - - % set invalid IC residual vectors equal to 0 (no feedback control) - R = reshape( R, [numel(MaskBad), dim] ); - R(MaskBad,:) = 0; - R = reshape( R, [size(MaskBad), dim] ); +function [FPad, GPad, MaskDomPad, MaskCtrlPad, MuPad, npad] = ... + pad (F, G, MaskDom, MaskCtrl, Mu, minScale) +% IN F forward DVF [NX x NY x NZ x 3] +% G inverse DVF estimate [NX x NY x NZ x 3] +% MaskDom valid domain mask [NX x NY x NZ] (logical) +% MaskCtrl controllable domain mask[NX x NY x NZ] (logical) +% Mu feedback control values [K-vector | NX x NY x NZ] +% minScale min domain scaling [scalar] +% OUT FPad padded forward DVF [PX x PY x PZ x 3] +% GPad padded inverse DVF [PX x PY x PZ x 3] +% MaskDomPad padded valid domain [PX x PY x PZ] (logical) +% MaskCtrlPad padded controllable dom.[PX x PY x PZ] (logical) +% MuPad padded control values [K-vector | PX x PY x PZ] +% npad dimension-wise pad size [(PX-NX) (PY-NY) (PZ-NZ)]/2 + + % scaling-adpative pad size (to avoid missing boundary values when down- + % and up-sampling with IMRESIZE in the case of multi-scale computation) + [~, dim] = dvf.sizeVf( F ); + npad = repmat( ceil( 1 / minScale ), [1 dim] ); + + % pad spatial domain in relevant vector- and scalar-field arrays + FPad = padarray( F, [npad 0], 0, 'both' ); + GPad = padarray( G, [npad 0], 0, 'both' ); + MaskDomPad = padarray( MaskDom , npad, false, 'both' ); + MaskCtrlPad = padarray( MaskCtrl, npad, false, 'both' ); + if ~isvector( Mu ) + MuPad = padarray( Mu, npad, 0, 'both' ); + else + MuPad = Mu; + end end %% ================================================== -% IC RESIDUAL MAGNITUDE PERCENTILES - -function prctR = icrMgnPercentiles (R, Mask, prcPop, flag) -% IN R IC residual [NX x NY x NZ x 3] -% Mask spatial domain mask [NX x NY x NZ] (logical) -% prcPop ICR percentile ranks [P-vector] -% flag compute ICR percentiles?[boolean] -% OUT prctR ICR magnitude prctiles [P-vector] - - if flag % ---- compute ICR percentiles - dim = ndims( Mask ); - RMgn = sqrt( sum( R.^2, dim+1 ) ); - [~, RMgn] = sanitizeIcResidual( R, RMgn ); - prctR = prctile( RMgn(Mask), prcPop, 'all' ); - else % ---- do not compute anything - prctR = NaN; - end % if +% DE-PAD SPATIAL DOMAIN + +function [GDepad, MaskDepad, MuDepad] = depad (G, Mask, Mu, npad) +% IN G inverse DVF estimate [PX x PY x PZ x 3] +% Mask domain mask [PX x PY x PZ] (logical) +% Mu feedback control values [K-vector | PX x PY x PZ] +% npad dimension-wise pad size [(PX-NX) (PY-NY) (PZ-NZ)]/2 +% OUT GDepad de-padded inverse DVF [NX x NY x NZ x 3] +% MaskDepad de-padded domain mask [NX x NY x NZ] (logical) +% MuDepad de-padded control values[K-vector | NX x NY x NZ] + + % dimension-wise internal (de-padded) spatial domain array indices + [szDomPad, ~] = dvf.sizeVf( G ); + iDepad = arrayfun( @(p,s) (1+p : s-p), npad, szDomPad, ... + 'UniformOutput', false ); + + % de-pad spatial domain in relevant vector- and scalar-field arrays + GDepad = G(iDepad{:},:); + MaskDepad = Mask(iDepad{:}); + if ~isvector( Mu ) + MuDepad = Mu(iDepad{:}); + else + MuDepad = Mu; + end end @@ -816,16 +906,22 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % -% 1.0.1 - November 07, 2018 +% 1.0.2 - Ferbuary 11, 2019 % % CHANGELOG % +% 1.0.2 (Feb 11, 2019) - Alexandros +% ! removed upper limit on neighborhood size for local +% control-pararameter value search +% + added scaling-factor-dependent pre-padding of spatial domain +% to obviate boundary interpolation artifacts due to resizing +% % 1.0.1 (Nov 07, 2018) - Alexandros % ! fixed error in non-invertible region calculation with 2D DVFs % % 1.0.0 (Oct 31, 2018) - Alexandros -% . initial release +% . initial version % % ------------------------------------------------------------ diff --git a/+dvf/issingularJacobian.m b/+dvf/issingularJacobian.m index 577f5bb..b0492b5 100644 --- a/+dvf/issingularJacobian.m +++ b/+dvf/issingularJacobian.m @@ -115,7 +115,7 @@ % 2D/3D DETERMINANT (EXPLICIT CALCULATION) function detJ = detJacobian (J) -% IN J Jacobian matrices [N x D x D] +% IN J Jacobian matrices [N x D x D] (D = 2,3) % OUT detJ Jacobian determinants [N x 1] switch size(J,3) @@ -151,7 +151,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.1 - November 7, 2018 % diff --git a/+dvf/jacobian.m b/+dvf/jacobian.m index db138fb..7f1c51b 100644 --- a/+dvf/jacobian.m +++ b/+dvf/jacobian.m @@ -108,7 +108,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+dvf/maskDomain.m b/+dvf/maskDomain.m index 23f8f99..7dc4d1b 100644 --- a/+dvf/maskDomain.m +++ b/+dvf/maskDomain.m @@ -55,8 +55,9 @@ % % [1] A. Dubey*, A.-S. Iliopoulos*, X. Sun, F.-F. Yin, and L. Ren, % "Iterative inversion of deformation vector fields with feedback -% control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, May 2018. -% DOI: 10.1002/mp.12962. +% control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, 2018. +% - DOI: 10.1002/mp.12962 +% - arXiv: 1610.08589 [cs.CV] % % % See also dvf.inversion, dvf.coorDisplace, imclose @@ -76,6 +77,9 @@ % max displacement in each dimension [1 x D] maxDisp = ceil( reshape( max( max( abs(F), [], 1:dim ), 0 ), [1 dim] ) ); + if isa( F, 'gpuArray' ) + maxDisp = gather( maxDisp ); + end %% DOMAIN MASK @@ -92,7 +96,7 @@ % deformed domain (equivalent to nearest-neighbor scattered interpolation % from deformed coordinates onto padded reference grid) - MaskDfm = false( szDom + 2*maxDisp + 1 ); + MaskDfm = false( szDom + 2*maxDisp + 1 ); % *padded domain range) gridDfm = cellfun( @(x,d) round(x) + d, gridDfm, num2cell(maxDisp), ... 'UniformOutput', false ); idxDfm = sub2ind( size(MaskDfm), gridDfm{:} ); @@ -119,7 +123,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+dvf/ntdcMeasures.m b/+dvf/ntdcMeasures.m index f5f8ab2..48fca3d 100644 --- a/+dvf/ntdcMeasures.m +++ b/+dvf/ntdcMeasures.m @@ -49,10 +49,12 @@ % deformation Jacobian and its eigevalues using the input DVF, and then % calculates the three spectral maps. % -% *NOTE* The distinction between the two is made solely by whether or -% not the input array is of 'real' or 'complex' type. If the input is -% to be eigenvalues that happen to be strictly real, it should be -% passed as complex(LAMBDA). +% NOTES +% +% The distinction between the two cases (eigenvalue vs DVF input) is made +% solely by whether or not the input array is of 'real' or 'complex' +% type. If the input array contains eigenvalues that happen to be +% strictly real, it should be passed as complex(LAMBDA). % % MEASURES % @@ -98,8 +100,9 @@ % % [1] A. Dubey*, A.-S. Iliopoulos*, X. Sun, F.-F. Yin, and L. Ren, % "Iterative inversion of deformation vector fields with feedback -% control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, May 2018. -% DOI: 10.1002/mp.12962. +% control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, 2018. +% - DOI: 10.1002/mp.12962 +% - arXiv: 1610.08589 [cs.CV] % % % See also dvf.jacobian, dvf.eigJacobian, dvf.feedbackControlVal @@ -153,7 +156,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+dvf/rescale.m b/+dvf/rescale.m index d87bfef..b5d2f91 100644 --- a/+dvf/rescale.m +++ b/+dvf/rescale.m @@ -174,7 +174,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+dvf/sizeVf.m b/+dvf/sizeVf.m index cd995e7..dab3738 100644 --- a/+dvf/sizeVf.m +++ b/+dvf/sizeVf.m @@ -65,7 +65,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+util/imageGrid.m b/+util/imageGridSmooth.m similarity index 63% rename from +util/imageGrid.m rename to +util/imageGridSmooth.m index 6c92ca8..bc7e47a 100644 --- a/+util/imageGrid.m +++ b/+util/imageGridSmooth.m @@ -1,6 +1,6 @@ -function I = imageGrid (sizeI, width, bufferBnd) +function I = imageGridSmooth (sizeI, width, bufferBnd) % -% IMAGEGRID - Grid-like 2D image generation +% IMAGEGRIDSMOOTH - Smooth grid-like image generation % % ---------------------------------------------------------------------- % @@ -23,30 +23,30 @@ % % SYNTAX % -% I = IMAGEGRID( SIZE ) -% I = IMAGEGRID( SIZE, WIDTH ) -% I = IMAGEGRID( SIZE, WIDTH, BUFFER ) +% I = IMAGEGRIDSMOOTH( SIZE ) +% I = IMAGEGRIDSMOOTH( SIZE, WIDTH ) +% I = IMAGEGRIDSMOOTH( SIZE, WIDTH, BUFFER ) % % INPUT % % SIZE Size of output image [D-vector] % [D1, D2, ...] -% WIDTH Checkered-pattern strips width [D-vector|scalar] -% {default: SIZE/21} +% WIDTH Grid spacing (sine half-period) [D-vector|scalar] +% {default: SIZE/10} % BUFFER Dimension-wise buffer from [D-vector|scalar] -% checkered-pattern to boundary +% image to boundary % {default: 0} % % OUTPUT % -% I Checkered-pattern image [D1-by-D2-by-...] +% I Smooth-grid image [D1-by-D2-by-...] % % DESCRIPTION % -% I = IMAGEGRID(SIZE,WIDTH,BUFFER) generates an image with a -% D-dimensional checkered pattern as per dimension-wise specified strip -% widths. The image intensities are normalized such that they are K/D, -% where K is the number of coinciding strips (for each pixel). +% I = IMAGEGRIDSMOOTH(SIZE,WIDTH,BUFFER) generates an image with a +% D-dimensional smooth grid-like pattern. The image is a Kronecker +% product of 1-dimensional sinusoids with specified half-periods. The +% image intensities are normalized to [0,1] range. % % If BUFFER is non-zero for the d-th dimension, then BUFFER(d) points are % masked off on both ends of the d-th axis. @@ -59,7 +59,7 @@ %% DEFAULTS if ~exist( 'width', 'var' ) || isempty( width ) - width = sizeI / 21; + width = sizeI / 10; end if~exist( 'bufferBnd', 'var' ) || isempty( bufferBnd ) bufferBnd = 0; @@ -71,8 +71,7 @@ % number of dimensions nDim = length( sizeI ); - % if single strip/buffer width was input, use the same for all - % dimensions + % if single width/buffer width was input, use the same for all dimensions if isscalar( width ) width = repmat( width, [nDim, 1] ); end @@ -86,10 +85,10 @@ bufferBnd = reshape( bufferBnd, [1, nDim] ); - %% DIMENSION-WISE (1D) CHECKERED PATTERNS AND BUFFER-MASKS + %% DIMENSION-WISE (1D) IMAGE PROFILE AND BUFFER-MASKS - chk1d = arrayfun( @(w,d,sz) ... - reshape( (mod( 0:(sz-1), 2*w ) < w), ... + img1d = arrayfun( @(w,d,sz) ... + reshape( sin( pi*w * linspace(0,1,sz) ), ... [ones(1,d-1), sz, 1] ), ... width, (1:nDim), sizeI, ... 'UniformOutput', false ); @@ -101,15 +100,15 @@ 'UniformOutput', false ); - %% CHECKERED-PATTERN IMAGE (WITH BUFFER-ZONES) + %% MULTI-DIMENSIONAL IMAGE (WITH BUFFER-ZONES) - % initialize checkered-pattern image and buffer-zone mask - I = chk1d{1}; + % initialize image and buffer-zone mask + I = img1d{1}; M = mask1d{1}; - % Kronecker summation of all remaining 1D patterns + % Kronecker product of all remaining 1D profiles for d = 2 : nDim - I = I + chk1d{d}; + I = I .* img1d{d}; M = M .* mask1d{d}; end @@ -127,9 +126,21 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % -% 1.0.0 - October 31, 2018 +% 1.0.2 - February 11, 2019 +% +% CHANGELOG +% +% 1.0.2 (Feb 11, 2019) - Alexandros +% . changed 1D image profile from binary-valued (sharp edges) to +% sinusoidal (smooth) +% . multiplicative composition of 1D image profiles into ND image +% (used to be additive) +% . renamed function: imageGrid --> imageGridSmooth +% +% 1.0.0 (Oct 31, 2018) - Alexandros +% . initial release % % ------------------------------------------------------------ diff --git a/+util/ndgridGpu.m b/+util/ndgridGpu.m index 72205da..70ad256 100644 --- a/+util/ndgridGpu.m +++ b/+util/ndgridGpu.m @@ -104,7 +104,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+util/parseOptArgs.m b/+util/parseOptArgs.m index f52ad9d..e0001a5 100644 --- a/+util/parseOptArgs.m +++ b/+util/parseOptArgs.m @@ -106,7 +106,7 @@ % % Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu % -% VERSION +% RELEASE % % 1.0.0 - October 31, 2018 % diff --git a/+vis/imageSlices.m b/+vis/imageSlices.m new file mode 100644 index 0000000..2e241cb --- /dev/null +++ b/+vis/imageSlices.m @@ -0,0 +1,165 @@ +function imageSlices (I, varargin) +% +% IMAGESLICES - Slice visualization of 3D image in subplots +% +% ---------------------------------------------------------------------- +% +% Copyright (C) 2018, Department of Computer Science, Duke University +% +% This program is free software: you can redistribute it and/or modify it +% under the terms of the GNU General Public License as published by the +% Free Software Foundation, either version 3 of the License, or (at your +% option) any later version. +% +% This program is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +% Public License for more details. +% +% You should have received a copy of the GNU General Public License along +% with this program. If not, see . +% +% ---------------------------------------------------------------------- +% +% SYNTAX +% +% IMAGESLICES( I ) +% IMAGESLICES( I, 'Name', Value, ... ) +% +% INPUT +% +% I 3D image [NX x NY x NZ] +% +% OPTIONS +% +% 'Slices' [3-vector] {ceil( [NX NY NZ] / 2 )} +% +% Image slice indices (along X, Y, and Z axes). +% +% 'Colormap' [C-by-3] {gray} +% +% IMAGESC colormap. +% +% 'CAxis' [2-vector | 'auto'] {'auto'} +% +% Color-axis limits (input to CAXIS) function. +% +% 'Rows' [scalar] {1} +% +% Number of subplot rows. +% +% 'Columns' [scalar] {3} +% +% Number of subplot columns. +% +% 'StartIndex' [scalar] {1} +% +% First subplot index. The slices are visualized in Rows x Columns +% subplots, from index StartIndex (XY/axial slice), to indices +% StartIndex+1 (XZ/coronal slice) and StartIndex+2 (YZ/sagittal slice). +% +% 'Title' [string] {''} +% +% Title string to be shared among slice subplots. +% +% 'PColor' [boolean] {false} +% +% If TRUE, then IMAGESLICES visualizes the image slices using PCOLOR +% instead of IMAGESC. This places a little more strain on the +% renderer, but can be useful for distinguishing NaN values. +% +% DESCRIPTION +% +% IMAGESLICES is a convenience function for visualizing 3 dimension-wise +% slices of a 3D image in subplots. +% +% DEPENDENCIES +% +% util.parseOptArgs +% +% +% See also imagesc, pcolor, slice, volumeViewer +% + + + %% PARAMETERS + + % options and default values + opt.slices = ceil( size(I) / 2 ); + opt.colormap = gray; + opt.caxis = 'auto'; + opt.rows = 1; + opt.columns = 3; + opt.startIndex = 1; + opt.title = ''; + opt.pcolor = false; + + % parse optional arguments + opt = util.parseOptArgs( opt, varargin{:} ); + + + %% IMAGESC vs PCOLOR + + if opt.pcolor + hPlot = @pcolor; + else + hPlot = @imagesc; + end + + + %% VISUALIZATION + + % XY slice + subplot( opt.rows, opt.columns, opt.startIndex ) + hPlot( squeeze( I(:,:,opt.slices(3)) ).' ) + axis image + shading flat + caxis( opt.caxis ) + colormap( opt.colormap ) + colorbar + title( {opt.title, sprintf( 'XY slice (z = %d)', opt.slices(3) )} ) + drawnow + + % XZ slice + subplot( opt.rows, opt.columns, opt.startIndex+1 ) + hPlot( squeeze( I(:,opt.slices(2),:) ).' ) + axis image + shading flat + caxis( opt.caxis ) + colormap( opt.colormap ) + colorbar + title( {opt.title, sprintf( 'XZ slice (y = %d)', opt.slices(2) )} ) + drawnow + + % YZ slice + subplot( opt.rows, opt.columns, opt.startIndex+2 ) + hPlot( squeeze( I(opt.slices(1),:,:) ).' ) + axis image + shading flat + caxis( opt.caxis ) + colormap( opt.colormap ) + colorbar + title( {opt.title, sprintf( 'YZ slice (x = %d)', opt.slices(1) )} ) + drawnow + + +end + + + +%%------------------------------------------------------------ +% +% AUTHORS +% +% Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu +% +% RELEASE +% +% 1.0.2 - February 11, 2019 +% +% CHANGELOG +% +% 1.0.2 (Feb 11, 2019) - Alexandros +% . initial version +% +% ------------------------------------------------------------ diff --git a/+vis/mfigure.m b/+vis/mfigure.m new file mode 100644 index 0000000..e542f8d --- /dev/null +++ b/+vis/mfigure.m @@ -0,0 +1,100 @@ +function [varargout] = mfigure (varargin) +% +% MFIGURE - Figure object with screen-adaptive size +% +% SYNTAX +% +% MFIGURE +% MFIGURE( SCALE ) +% H = MFIGURE( ... ) +% +% INPUT +% +% SCALE Figure-size to screen-size ratio [scalar] +% {default: 0.6} +% +% OUTPUT +% +% H Figure object handle [matlab.ui.figure] +% +% DESCRIPTION +% +% MFIGURE is a wrapper around the FIGURE function. It creates a figure +% object whose size is set relative to the screen size. This can be +% useful when using large, high-resolution screens. +% +% DEPENDENCIES +% +% +% +% +% See also figure +% + + + % default value for optional scaling factor + switch nargin + case 0 + figureScale = 0.6; + case 1 + figureScale = varargin{1}; + end + + % Get screen parameters + graphicalRoot = groot; % Get groot handle + screenSizeGroot = graphicalRoot.ScreenSize(3:4); % Get screen size + % width and height + + % Compute screen offset + figSizeEdgeOffset = (1 - figureScale) * graphicalRoot.ScreenSize(3:4) / 2; + + % Open figure + hf = figure; % Initialize figure + + % Get/set units + figUnits = hf.Units; % Get current figure units (users may change defaults) + hf.Units = graphicalRoot.Units; % Force units the same + + % Compute figure size in terms of width and height + figSize = screenSizeGroot - figSizeEdgeOffset*2; % width, height offsets + + % Position and resize figure + hf.Position = [figSizeEdgeOffset(1) figSizeEdgeOffset(2) ... + figSize(1) figSize(2)]; % left bottom width height + + % Set figure units back + hf.Units = figUnits; + + % optional output: figure handle + if nargout > 0 + varargout{1} = hf; + end + + +end + + + +%%------------------------------------------------------------ +% +% AUTHORS +% +% Kevin Mattheus Moerman +% Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu +% +% RELEASE +% +% 1.0.2 - February 11, 2019 +% +% CHANGELOG +% +% 1.0.2 (February 11, 2019) - Alexandros +% + rudimentary documentation +% +% 0.0.1 (Dec 10, 2018) - Kevin +% . initial version +% [https://github.com/ailiop/idvf/issues/4#issue-389261860] +% +% ------------------------------------------------------------ + + diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md new file mode 100644 index 0000000..20ae749 --- /dev/null +++ b/CODE_OF_CONDUCT.md @@ -0,0 +1,87 @@ +## Contributor Covenant Code of Conduct + + +### Our Pledge + +In the interest of fostering an open and welcoming environment, we as +contributors and maintainers pledge to making participation in our project +and our community a harassment-free experience for everyone, regardless of +age, body size, disability, ethnicity, sex characteristics, gender identity +and expression, level of experience, education, socio-economic status, +nationality, personal appearance, race, religion, or sexual identity and +orientation. + + +### Our Standards + +Examples of behavior that contributes to creating a positive environment +include: + +- Using welcoming and inclusive language +- Being respectful of differing viewpoints and experiences +- Gracefully accepting constructive criticism +- Focusing on what is best for the community +- Showing empathy towards other community members + +Examples of unacceptable behavior by participants include: + +- The use of sexualized language or imagery and unwelcome sexual attention + or advances +- Trolling, insulting/derogatory comments, and personal or political + attacks +- Public or private harassment +- Publishing others' private information, such as a physical or electronic + address, without explicit permission +- Other conduct which could reasonably be considered inappropriate in a + professional setting + + +### Our Responsibilities + +Project maintainers are responsible for clarifying the standards of +acceptable behavior and are expected to take appropriate and fair +corrective action in response to any instances of unacceptable behavior. + +Project maintainers have the right and responsibility to remove, edit, or +reject comments, commits, code, wiki edits, issues, and other contributions +that are not aligned to this Code of Conduct, or to ban temporarily or +permanently any contributor for other behaviors that they deem +inappropriate, threatening, offensive, or harmful. + + +### Scope + +This Code of Conduct applies both within project spaces and in public +spaces when an individual is representing the project or its community. +Examples of representing a project or community include using an official +project e-mail address, posting via an official social media account, or +acting as an appointed representative at an online or offline event. +Representation of a project may be further defined and clarified by project +maintainers. + + +### Enforcement + +Instances of abusive, harassing, or otherwise unacceptable behavior may be +reported by contacting the project team at idvf@cs.duke.edu. All +complaints will be reviewed and investigated and will result in a response +that is deemed necessary and appropriate to the circumstances. The project +team is obligated to maintain confidentiality with regard to the reporter +of an incident. Further details of specific enforcement policies may be +posted separately. + +Project maintainers who do not follow or enforce the Code of Conduct in +good faith may face temporary or permanent repercussions as determined by +other members of the project's leadership. + + +### Attribution + +This Code of Conduct is adapted from the [Contributor Covenant][homepage], +version 1.4, available at +https://www.contributor-covenant.org/version/1/4/code-of-conduct.html + +[homepage]: https://www.contributor-covenant.org + +For answers to common questions about this code of conduct, see +https://www.contributor-covenant.org/faq diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..7fcbfc0 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,134 @@ +## Contributing to idvf + + +We welcome contributions that help improve the `idvf` project! Please read +through the guidelines in this document before reporting an issue or +submitting a request. We will do our best to respond to all issues and +requests, but please bear in mind that it may take us a while. + + + + + +### Contents + + +- [Types of contributions](#contrib-types) +- [Bug reports](#bug-reports) +- [Pull requests](#pull-requests) +- [Code and documentation style](#style) +- [Feature requests](#feature-requests) + + + + + +### Types of contributions + + +- *Bug reports.* +- *Compatibility patches:* support for earlier MATLAB versions; + alternatives to MATLAB Toolbox functions. +- *Minor patches:* documentation clarification and typos; + code/documentation formatting; naming convention amendments. +- *Functionality updates.* +- *Testing:* testing or demo scripts for existing functionality. +- *Re-implementations:* performance improvements; system/language support + extensions. +- Anything else, as long as its utility and functionality is described. + + + + + +### Bug reports + + +Please [open a new issue][github-new-issue] for each unreported bug. +Specify "[BUG] *1-sentence-description-of-bug*" as the issue title, and +list the following information in the issue body: + +- Brief summary and background. +- Bug description: what should happen, and what happens instead. +- Version of MATLAB and relevant toolboxes; operating system; GPU + hardware and drivers (if the bug affects GPU execution). +- Code for a concise MATLAB script that reproduces and illustrates the + bug. +- Any other relevant notes (e.g., what you think causes the bug, any + steps you may have taken to identify or resolve it, etc). + + +[github-new-issue]: https://help.github.com/articles/creating-an-issue/ + + + + + +### Pull requests + + +Please [submit a GitHub pull request][github-pull-request] for each code or +documentation contribution to `idvf`. When submitting a pull request, +please adhere to the following. + +- Clearly identify the [type of your contribution](#contrib-types) in the + title and body of your pull request. + - If your contributions span multiple types, please separate them + into individual pull requests. Minor patches should be lumped into + a single pull request. +- Include a brief description of the rationale, functionality, and + implementation of your contribution. +- Include a testing or demo MATLAB script (named `test_xxx.m` or + `demo_xxx.m`) that can be used to illustrate and validate your + contribution. Include a brief description of the script in the pull + request body. +- [Squash partial commits][github-squash-commit]. +- If applicable, draft some relevant text to be added to or amended in + the README. Please include the text in the pull request body, *not* as + part of the commit. + +We encourage you to open a new issue to discuss any intended contributions +prior to developing or submitting a pull request. + + +[github-pull-request]: https://help.github.com/articles/about-pull-requests/ + +[github-squash-commit]: https://help.github.com/articles/about-pull-request-merges/ + + + + + +### Code and documentation style + + +Please try to follow the style conventions in the `idvf` repository when +submitting pull requests. We recommend that you use the source code in +`+dvf/inversion.m` as a style template for functions, and +`demo_inversion_2d.m` for scripts. We generally try to observe the +following rules: + +- The code should be clear, stable, and efficient. Clarity and stability + take precedence over efficiency and performance. The code should be + self-documented if possible (avoid referring to descriptions in + existing issues or pull requests). +- Function documentation should be comprehensive and follow the format of + existing functions (e.g., `+dvf/inversion.m`). +- Function and variable names are in `camelCase`; script names are in + `snake_case`. Generally, matrix/array names start with an uppercase + letter, while scalar/vector/function names start with a lowercase + letter. +- All code blocks should be briefly documented. +- We prefer 4-space indentation (no tabs), operator/operand alignment + across multiple lines, and 80-column line width. + + + + + +### Feature requests + + +It is unlikely that we will do much development for new features, unless +they are essential or supported by theoretical advances. We do, however, +encourage the development and submission of new features via pull requests. diff --git a/README.md b/README.md index 9e51f2a..f0601ce 100644 --- a/README.md +++ b/README.md @@ -1,46 +1,57 @@ -[![DOI](https://zenodo.org/badge/153499710.svg)](https://zenodo.org/badge/latestdoi/153499710) -[![GitHub release](https://img.shields.io/github/release/ailiop/idvf.svg)](https://GitHub.com/ailiop/idvf/releases/) -[![GitHub license](https://img.shields.io/github/license/ailiop/idvf.svg)](https://github.com/ailiop/idvf/blob/master/LICENSE) -[![Github all releases](https://img.shields.io/github/downloads/ailiop/idvf/total.svg)](https://GitHub.com/ailiop/idvf/releases/) -[![GitHub issues](https://img.shields.io/github/issues/ailiop/idvf.svg)](https://GitHub.com/ailiop/idvf/issues/) +# idvf: _Iterative inversion of deformation vector field with adaptive bi-residual feedback control_ + +[![Zenodo DOI](https://zenodo.org/badge/153499710.svg)](https://zenodo.org/badge/latestdoi/153499710) +[![GitHub release](https://img.shields.io/github/release/ailiop/idvf.svg)](https://github.com/ailiop/idvf/releases/) +[![GitHub license](https://img.shields.io/github/license/ailiop/idvf.svg)](https://github.com/ailiop/idvf/blob/master/LICENSE) +[![Github all releases](https://img.shields.io/github/downloads/ailiop/idvf/total.svg)](https://github.com/ailiop/idvf/releases/) +[![GitHub issues](https://img.shields.io/github/issues/ailiop/idvf.svg)](https://github.com/ailiop/idvf/issues/) -# idvf: _Iterative inversion of deformation vector field
with adaptive bi-residual feedback control_ - + ## Contents + - [Algorithm overview](#algorithm) - [Two inverse consistency residuals](#ic-residuals) - [DVF inversion iteration with bi-residual feedback control](#inversion-iteration) - [References](#references) +- [Getting started](#getting-started) + - [Installation](#installation) + - [Testing](#testing) - [Software description](#software) - [DVF inversion](#inversion-function) - [Preprocessing & DVF characterization](#preprocessing) - [Post-evaluation](#post-evaluation) - [GPU utilization](#gpu) + - [Eigenvalue calculations](#eigenvalues) - [System environment](#system-reqs) +- [License and community guidelines](#license-contrib-reports) - [Contributors](#contributors) - + + ## Algorithm overview + This repository contains a set of MATLAB functions for fast and accurate inversion of 2D/3D deformation vector fields (DVFs), which are also known -as flow or dense motion fields. The inversion algorithm framework and a -unified analysis are described in References [1,2]. We give a brief -overview of the DVF inversion algorithm which underlies the repository -functions. +as flow or dense motion fields. The inversion algorithm framework and a +unified analysis are described in References +[[1](#medphys2018),[2](#phddubey2018)]. To our knowledge, this is the +first DVF inversion framework that is supported by theoretical guarantees +under mild and verifiable conditions. We give a brief overview of the DVF +inversion algorithm which underlies the repository functions. + We assume that we are provided with a deformation vector field (DVF) -denoted by _u_. The DVF describes a (non-linear) mapping from a reference +denoted by _u_. The DVF describes a (non-linear) mapping from a reference image _Ir_ onto a study image _Is_ via point-wise displacement, i.e., @@ -48,7 +59,7 @@ displacement, i.e., Ir(y) = Is(y + u(y)),

-at all locations _y_ in the reference domain. We aim at computing the +at all locations _y_ in the reference domain. We aim at computing the inverse DVF _v_ such that

@@ -59,12 +70,13 @@ at all _x_ in the study domain. - + ### Two inverse consistency residuals + If two DVFs _u_ and _v_ are inverse of each other, they must meet the -*inverse consistency* (IC) condition. We define the study IC residual, +*inverse consistency* (IC) condition. We define the study IC residual,

s(x) = v(x) + u(x + v(x)), @@ -76,28 +88,32 @@ and the reference IC residual, r(y) = u(y) + v(y + u(y)).

-The DVFs _u_ and _v_ are inverse-consistent if both residuals are zero. The -residuals are related to the inversion error, _e(x) = v(x) − v*(x)_, -where _v*_ is the true inverse DVF. The IC residuals are -computationally available, while the inversion error is unknown. We use the -IC residuals as feedback in iterative DVF inversion, and exercise adaptive -control over the feedback for global convergence and local acceleration. +The DVFs _u_ and _v_ are inverse-consistent if both residuals are zero. +The residuals are related to the inversion error, _e(x) = v(x) − +v*(x)_, where _v*_ is the true inverse DVF. The IC residuals are +computationally available, while the inversion error is unknown. We use +the IC residuals as feedback in iterative DVF inversion, and exercise +adaptive control over the feedback for global convergence and local +acceleration. - + ### DVF inversion iteration with bi-residual feedback control + Due to the non-linear nature of the DVF mapping, one resorts to iterative -procedures for DVF inversion. Our procedure has two stages: preprocessing +procedures for DVF inversion. Our procedure has two stages: preprocessing and adaptive inversion iteration. + In the preprocessing stage, we compute the DVF Jacobians and their -eigenvalues over the spatial domain. We use this spectral information to +eigenvalues over the spatial domain. We use this spectral information to determine data-adaptive parameters for feedback control. -The inversion iteration proceeds in two phases. During phase one, at step + +The inversion iteration proceeds in two phases. During phase one, at step _k+1_, we update the current inverse estimate _vk_ using the study IC residual _sk_ as feedback: @@ -106,14 +122,13 @@ study IC residual _sk_ as feedback: vk(x) − (1 − μk(x)) sk(x),

-where _μ_ is the feedback control parameter. With the adaptive control +where _μ_ is the feedback control parameter. With the adaptive control schemes provided in this repository, the iteration is guaranteed to -converge globally, under certain mild conditions [1]. +converge globally, under certain mild conditions [[1](#medphys2018)]. + -Inversion errors can be estimated using the IC residuals and the DVF -spectral information. Once the errors are made sufficiently small, the -iteration is switched to phase two: +Once the errors are made sufficiently small, the iteration is switched to +phase two:

vk+1(x) = @@ -121,138 +136,246 @@ iteration is switched to phase two:

where the reference IC residual _rk_ at displaced locations is -used as feedback. The local convergence rate is quadratic during phase -two. We may refer to this iteration as implicit Newton iteration. However, -phase-two iteration steps do not entail explicit formation and inversion of -the Newton matrix; the explicit Newton step suffers from several numerical -problems. +used as feedback. Inversion errors can be estimated using the IC residuals +and the DVF spectral information. The local convergence rate is quadratic +during phase two. We may refer to this iteration as implicit Newton +iteration. However, phase-two iteration steps do not entail explicit +formation and inversion of the Newton matrix; the explicit Newton step +suffers from several numerical problems. -Phase transition and integration in the inversion iteration is enabled by a -multi-resolution scheme, elaborated in Reference [2]. +Phase transition and integration in the inversion iteration is facilitated +by a multi-resolution scheme, elaborated in Reference [[2](#phddubey2018)]. - + + ### References -[1] A. Dubey*, A.S. Iliopoulos*, X. Sun, -F.F. Yin, and L. Ren, "Iterative inversion of -deformation vector fields with feedback control", *Medical Physics*, -45(7):3147-3160, 2018. -[2] A. Dubey, *Symmetric Completion of Deformable -Registration via Bi-residual Inversion*, PhD dissertation, Duke University, -Durham, NC, USA, 2018. + + +[1] A. Dubey*, A.S. Iliopoulos*, X. Sun, F.F. Yin, and L. Ren, +["Iterative inversion of deformation vector fields with feedback +control"][medphys2018-doi], *Medical Physics*, 45(7):3147-3160, 2018. +[(arXiv)][medphys2018-arxiv] + +[medphys2018-doi]: https://doi.org/10.1002/mp.12962 + +[medphys2018-arxiv]: https://arxiv.org/abs/1610.08589 + + + + +[2] A. Dubey, *Symmetric completion of deformable registration via +bi-residual inversion*, PhD thesis, Duke University, Durham, NC, USA, 2018. + + + + + + +## Getting started + + + + + +### Installation + + +To use idvf, simply add its top-level directory to the [MATLAB +path][matlab-path]. All functions are organized in +[packages][matlab-packages]. + +[matlab-path]: https://www.mathworks.com/help/matlab/matlab_env/what-is-the-matlab-search-path.html + +[matlab-packages]: https://www.mathworks.com/help/matlab/matlab_oop/scoping-classes-with-packages.html + + + + +### Testing + +To test idvf, run the included demo scripts (`demo_inversion_2d`, +`demo_inversion_3d_z0`, `demo_inversion_3d_zsin`). Each script +demonstrates inversion of a sample forward DVF using 8 different feedback +control settings, and produces the following visualizations: +- Deformation of a synthetic grid-like image by the forward DVF. +- Spectral measure maps of the forward DVF + ([preprocessing](#preprocessing)). +- IC residual magnitudes ([post-evaluation](#post-evaluation)) with each + control setting: percentiles at each step of the iteration, and spatial + maps at the end of the iteration. +- Image-space difference maps in reference image recovery with each + control setting. - + +The adaptive control schemes are guaranteed to converge; that is, the +corresponding IC residuals tend towards zero, at a rate that depends on the +control setting. The spectral measure maps provide useful information +regarding the invertibility and controllability of the forward DVF over its +spatial domain. + + + + + ## Software description - + ### DVF inversion -The main function is `dvf.inversion`. It takes a forward DVF _u_ and -returns the inverse DVF _v_. A 3D DVF is represented by a 4D array of size + +The main function is `dvf.inversion`. It takes a forward DVF _u_ and +returns the inverse DVF _v_. A 3D DVF is represented by a 4D array of size _Nx_ x _Ny_ x _Nz_ x _3_, in `single` or -`double` precision. Specifically, `U(i,j,k,1:3)` is the 3D forward -displacement vector _u(x)_ at voxel _x = (i,j,k)_. Two-dimensional DVFs are -represented similarly. +`double` precision. Specifically, `U(i,j,k,1:3)` is the 3D forward +displacement vector _u(x)_ at voxel _x = (i,j,k)_. Two-dimensional DVFs +are represented analogously. + We provide a number of options in function `dvf.inversion`, which are input -as name-value pairs. These options pertain to the choice of feedback -control schemes and iteration parameters. The function documentation -details the provided options. +as name-value pairs. These options pertain to the choice of feedback +control schemes and iteration parameters. The function documentation +details the supported options. + The memory requirement of `dvf.inversion` is approximately 5 times the amount of memory for the input DVF data array. -We provide two scripts, `demo_inversion_2d` and `demo_inversion_3d`, to -demonstrate how to use the repository functions. The scripts include -post-evaluation of the inversion results in terms of IC residual magnitude -percentiles throughout the iteration, per several control schemes. - - + ### Preprocessing & DVF characterization + Preprocessing functions are invoked automatically by the main function when -adaptive feedback control is chosen. Functions `dvf.jacobian` and +adaptive feedback control is chosen. Functions `dvf.jacobian` and `dvf.eigJacobian` compute the DVF Jacobians and their eigenvalues, -respectively, at all pixels/voxels. The function `dvf.feedbackControlVal` -takes the eigenvalues and returns adaptive values for the feedback control -parameter _μ_ per the chosen scheme. +respectively, at all pixels/voxels. Function `dvf.feedbackControlVal` +takes the eigenvalues and returns adaptive feedback control parameter +values per the chosen scheme. + -The function `dvf.ntdcMeasures` calculates three characteristic spectral -measures of a DVF. Specifically, it returns the algebraic control index +Function `dvf.ntdcMeasures` calculates three characteristic spectral +measures of a DVF. Specifically, it returns the algebraic control index map, the non-translational component spectral radius map, and the determinant map. - + ### Post-evaluation -The `dvf.inversion` function returns post-evaluation measures as optional -output. The measures are IC residual magnitude percentiles at each -iteration step. + +The two IC residuals enable post-evaluation of inverse DVF estimates in the +absence of the true inverse DVF. Function `dvf.inversion` returns IC +residual magnitude percentiles at each iteration step as optional output. +Function `dvf.icResidual` returns an IC residual field given a DVF pair +(reference or study IC residual, depending on the order of the input DVFs). + Certain regions of the study domain are invalid for IC residual -calculations. These regions cover points that are either mapped outside the -domain by the forward DVF, or lie outside the deformed reference -domain. Function `dvf.maskDomain` calculates the valid domain. +calculations. These regions cover points that are either mapped outside +the domain by the forward DVF, or lie outside the deformed reference +domain. Function `dvf.maskDomain` calculates the valid domain. - + ### GPU utilization + All `dvf.*` functions are GPU-enabled, supported by the MATLAB Parallel -Computing Toolbox. When a GPU is available, GPU computation is invoked by -casting the input data arrays into `gpuArray` objects. For example: +Computing Toolbox. When a compatible GPU is available, GPU computation is +invoked by casting the input data arrays into `gpuArray` objects; for +example: - V = dvf.inversion(gpuArray(U)) + V = dvf.inversion(gpuArray(U)); The output data arrays are returned as `gpuArray` objects. - + + +### Eigenvalue calculations + + +Eigenvalue calculations are currently slow with large DVFs. They are +implemented as a loop that calls the MATLAB `eig` function for each +Jacobian over the spatial domain. This approach is inefficient but has the +benefit of numerical stability. In our experience with respiratory DVFs +over thoracic and abdominal CT domains, we have not encountered any major +issues with efficient but numerically unstable alternatives such as +closed-form root calculation. Nevertheless, we do not include such +alternatives, in the interest of stability and simplicity over efficiency. + + + + + ## System environment -The repository code was developed and tested on MATLAB R2018b. The + +The repository code was developed and tested on MATLAB R2018b. The following functions in MATLAB toolboxes are used: - `imresize`, `imresize3`, `imdilate`, `imerode`, and `imclose` (Image - Processing Toolbox, version 10.3); -- `prctile` (Statistics and Machine Learning Toolbox, version 11.4); and -- `gpuArray` (Parallel Computing Toolbox, version 6.13). + Processing Toolbox, tested with version 10.3); +- `prctile` (Statistics and Machine Learning Toolbox, tested with version + 11.4); and +- `gpuArray` (Parallel Computing Toolbox, tested with version 6.13). + -GPU computation was tested on an NVIDIA Quantum TXR113-1000R GPU and CUDA +GPU computations were tested on an NVIDIA Quantum TXR113-1000R and CUDA 10.0 drivers. - + + + +## License and community guidelines + + +The idvf code is licensed under the [GNU general public license +v3.0][license]. If you wish to contribute to idvf or report any +bugs/issues, please see our [contribution guidelines][contrib] and [code of +conduct][conduct]. + +[license]: https://github.com/ailiop/idvf/blob/master/LICENSE + +[contrib]: https://github.com/ailiop/idvf/blob/master/CONTRIBUTING.md + +[conduct]: https://github.com/ailiop/idvf/blob/master/CODE_OF_CONDUCT.md + + + + ## Contributors -- *Design, development, and testing:* + +- *Design, development, testing:* Alexandros-Stavros Iliopoulos, Abhishek Dubey, and Xiaobai Sun Department of Computer Science, Duke University - *Consulting on medical CT image analysis and applications:* Lei Ren and Fang-Fang Yin Department of Radiation Oncology, Duke University Medical School + +- *Review, suggestions, bug reports:* + @Kevin-Mattheus-Moerman diff --git a/codemeta.json b/codemeta.json index 5ebf67a..4e0f822 100644 --- a/codemeta.json +++ b/codemeta.json @@ -25,11 +25,11 @@ "identifier": "https://doi.org/10.5281/zenodo.1476601", "codeRepository": "https://github.com/ailiop/idvf", "datePublished": "2018-11-01", - "dateModified": "2018-11-01", + "dateModified": "2019-02-11", "dateCreated": "2018-11-01", "description": "Iterative inversion of deformation vector field with adaptive bi-residual feedback control", "keywords": "deformation vector field, optical flow, deformable registration, fixed-point iteration, inverse consistency, adaptive feedback control, medical image analysis, computer vision, MATLAB", "license": "GPL v3.0", "title": "idvf", - "version": "v1.0.0" + "version": "v1.0.2" } diff --git a/demo_inversion_2d.m b/demo_inversion_2d.m index 8cc74c2..22e1a90 100644 --- a/demo_inversion_2d.m +++ b/demo_inversion_2d.m @@ -26,10 +26,12 @@ % [1] A. Dubey*, A.-S. Iliopoulos*, X. Sun, F.-F. Yin, and L. Ren, % "Iterative inversion of deformation vector fields with feedback % control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, 2018. -% DOI: 10.1002/mp.12962. +% - DOI: 10.1002/mp.12962 +% - arXiv: 1610.08589 [cs.CV] % % [2] A. Dubey, "Symmetric completion of deformable registration via -% bi-residual inversion," PhD thesis, Duke University, Durham, NC, USA. +% bi-residual inversion," PhD thesis, Duke University, Durham, NC, USA, +% 2018. % @@ -46,7 +48,6 @@ pathData = 'data/c-deformation.mat'; ioData = matfile( pathData ); F = ioData.F; -MaskEval = dvf.maskDomain( F ); [szDom, ~] = dvf.sizeVf( F ); % inversion parameters @@ -56,14 +57,14 @@ opt.stopThreshold = 1e-6; opt.accThreshold = 0; opt.InitialValue = []; -opt.Mask = MaskEval; +opt.Mask = dvf.maskDomain( F ); -% create an array of inversion parameters to compare multiple results -opt = repmat( opt, [7 1] ); +% create an array of inversion parameters to compare multiple schemes +opt = repmat( opt, [8 1] ); opt(1).control = 0.0; % constant mu = 0 opt(2).control = 0.5; % constant mu = 0.5 opt(3).control = 'alternating'; % alternating -opt(4).control = 'optimal'; % locally optimal +opt(4).control = 'optimal'; % pointwise optimal opt(5).control = 'midrange'; % local midrange opt(6).control = 'midrange'; % local midrange & acc. opt(6).accThreshold = 1; % . @@ -71,18 +72,29 @@ opt(7).accThreshold = 1; % . opt(7).scale = [0.5, 1.0]; % . opt(7).numIteration = [4 16]; % . +opt(8).control = 'optimal'; % pointwise optimal & acc. & multiscale +opt(8).accThreshold = 1; % . +opt(8).scale = [0.5, 1.0]; % . +opt(8).numIteration = [4 16]; % . % legend labels for each control setting lgnd = {'\mu = 0', '\mu = 0.5', '\mu_{oe}', '\mu_*(x)', '\mu_m(x)', ... - '\mu_m(x) & acc', '\mu_m(x) & acc & MS'}; + '\mu_m(x) & acc', '\mu_m(x) & acc & MS', '\mu_*(x) & acc & MS'}; -% error & IC residual curve visualization options +% IC residual percentile curve visualization options prcIcMeasure = [50 90 95 98 100]; argVisIcCurve = {'-o', 'LineWidth', 2}; +clrIcCurve = [0.8650 0.8110 0.4330; % yellow + 0.9718 0.5553 0.7741; % pink + 0.6859 0.4035 0.2412; % brown + 1.0000 0.5482 0.1000; % orange + 0.6365 0.3753 0.6753; % purple + 0.3718 0.7176 0.3612; % green + 0.2941 0.5447 0.7494; % blue + 0.9047 0.1918 0.1988]; % red -% error-map visualization options -climError = [0 30]; -cmapError = [1 1 1; jet(30)]; +% IC residual maps visualization options +climIcMgn = [0 1]; % color-axis limits %% ==================== (BEGIN) @@ -94,10 +106,11 @@ fprintf( '...visualizing synthetic reference and deformed images...\n' ); -IRef = util.imageGrid( szDom ); -IDfm = dvf.imdeform( IRef, F ); +IRef = util.imageGridSmooth( szDom, [15 15] ); +IStd = dvf.imdeform( IRef, F ); -figure +hFig = vis.mfigure; +hFig.Name = 'reference and study (deformed) images'; % reference image subplot( 1, 2, 1 ) @@ -107,9 +120,9 @@ title( 'reference image' ) drawnow -% deformed (study) image +% study (deformed) image subplot( 1, 2, 2 ) -imagesc( IDfm.' ) +imagesc( IStd.' ) axis image colormap( gray ) title( 'deformed image' ) @@ -120,9 +133,10 @@ fprintf( '...calculating & visualizing spectral NTDC measures...\n' ); -[CtrlIdx, Rho, Det] = dvf.ntdcMeasures( F ); +[CtrlIdx, Rho, Det, Lambda] = dvf.ntdcMeasures( F ); -figure +hFig = vis.mfigure; +hFig.Name = 'spectral NTDC measures of forward DVF'; % control index subplot( 1, 3, 1 ) @@ -160,16 +174,20 @@ Mu = cell( size(opt) ); prctRG = cell( size(opt) ); prctRF = cell( size(opt) ); +RG = cell( size(opt) ); +RF = cell( size(opt) ); for i = 1 : length(opt) fprintf( ' - scheme #%d (%s)...\n', i, lgnd{i} ); [G{i}, Mu{i}, ~, prctRG{i}, prctRF{i}] = dvf.inversion( F, opt(i) ); + RG{i} = dvf.icResidual( G{i}, F ); + RF{i} = dvf.icResidual( F, G{i} ); end %% ==================== VISUALIZATION: IC RESIDUALS PER ITERATION -fprintf( '...visualizing IC residual & error percentile curves...\n' ); +fprintf( '...visualizing step-wise IC residual percentile curves...\n' ); % resolution-normalized iteration steps kstep = cell( size(opt) ); @@ -180,11 +198,13 @@ kstep{i} = cumsum( horzcat( 0, kstep{i}{:} ) ); end -figure -numPrct = length( prcIcMeasure ); +hFig = vis.mfigure; +hFig.Name = 'step-wise IC residual percentile curves'; +numPrct = length( prcIcMeasure ); % legend -subplot( 5, numPrct, (1:numPrct) ) +subplot( 5, numPrct, (1:numPrct), ... + 'NextPlot', 'ReplaceChildren', 'ColorOrder', clrIcCurve ) plot( nan( 2, length(opt) ), argVisIcCurve{:} ); title( 'step-wise IC residual percentile curves' ) legend( lgnd, 'Orientation', 'horizontal', 'Location', 'south' ) @@ -193,26 +213,28 @@ for p = 1 : numPrct % ---- each percentile - % IC residual (RG) + % study IC residuals subplot( 5, numPrct, p + (1:2)*numPrct ) hold on for i = 1 : length(opt) - plot( kstep{i}, prctRG{i}(p,:), argVisIcCurve{:} ); + plot( kstep{i}, prctRG{i}(p,:), argVisIcCurve{:}, ... + 'Color', clrIcCurve(i,:) ); end - ylabel( sprintf( 'r_G[%d%%]', prcIcMeasure(p) ) ); - xlabel( 'amortizerd iteration step' ) + ylabel( sprintf( 's_k[%d%%]', prcIcMeasure(p) ) ); + xlabel( 'amortizerd iteration step (k)' ) set( gca, 'YScale', 'log' ) grid on drawnow - % IC residual (RF) + % reference IC residuals subplot( 5, numPrct, p + (3:4)*numPrct ) hold on for i = 1 : length(opt) - plot( kstep{i}, prctRF{i}(p,:), argVisIcCurve{:} ); + plot( kstep{i}, prctRF{i}(p,:), argVisIcCurve{:}, ... + 'Color', clrIcCurve(i,:) ); end - ylabel( sprintf( 'r_F[%d%%]', prcIcMeasure(p) ) ); - xlabel( 'amortized iteration step' ) + ylabel( sprintf( 'r_k[%d%%]', prcIcMeasure(p) ) ); + xlabel( 'amortized iteration step (k)' ) set( gca, 'YScale', 'log' ) grid on drawnow @@ -220,6 +242,81 @@ end % for (p) +%% ==================== VISUALIZATION: IC RESIDUAL MAPS + +fprintf( '...visualizing final IC residual maps...\n' ); +fprintf( ' - white pixels indicate NaN (out-of-bounds) values\n' ); + +% calculate point-wise IC residual magnitudes +RGMgn = cellfun( @(r) sqrt( sum( r.^2, 3 ) ), RG, 'UniformOutput', false ); +RFMgn = cellfun( @(r) sqrt( sum( r.^2, 3 ) ), RF, 'UniformOutput', false ); + +% study IC residuals +hFig = vis.mfigure; +hFig.Name = 'study IC residual maps'; +for i = 1 : length(opt) + subplot( 2, 4, i ) + pcolor( RGMgn{i}.' ) + axis image + shading flat + colormap jet + caxis( climIcMgn ) + colorbar + title( {lgnd{i}, 'study IC residuals s(x)'} ) + drawnow +end + +% reference IC residuals +hFig = vis.mfigure; +hFig.Name = 'reference IC residual maps'; +for i = 1 : length(opt) + subplot( 2, 4, i ) + pcolor( RFMgn{i}.' ) + axis image + shading flat + colormap jet + caxis( climIcMgn ) + colorbar + title( {lgnd{i}, 'reference IC residuals r(x)'} ) + drawnow +end + + +%% ==================== VISUALIZATION: RECOVERED REFERENCE IMAGE + +fprintf( '...reference image recovery...\n' ); +fprintf( ' - I_rec(x) := I_std(x + v(x))\t' ); +fprintf( ' [I_std(y) = I_ref(y + u(y))]\n' ); +fprintf( ' = I_ref(x + v(x) + u(x + v(x)))\n' ); +fprintf( ' = I_ref(x + s(x))\n' ); + +% ---------- calculation +fprintf( ' - calculation...\n' ); + +IRefRec = cell( size(opt) ); +for i = 1 : length(opt) + IRefRec{i} = dvf.imdeform( IRef, RG{i} ); +end + +% ---------- visualization +fprintf( ' - difference image visualization (I_ref - I_rec)...\n' ); + +hFig = vis.mfigure; +hFig.Name = 'image-space errors in reference image recovery'; +for i = 1 : length(opt) + subplot( 2, 4, i ) + pcolor( (IRef - IRefRec{i}).' ); + axis image + shading flat + caxis( [-1, +1] ) + colormap parula + colorbar + title( {lgnd{i}, 'I_{ref}(x) - I_{ref}(x+s(x))'} ) + drawnow +end + + + %% ==================== (END) fprintf( '\n***** END (%s) *****\n\n', mfilename ); @@ -234,11 +331,19 @@ % % VERSION % -% 0.1 - October 08, 2018 +% 0.2 - December 21, 2018 % % CHANGELOG % +% 0.2 (Dec 21, 2018) - Alexandros +% + added IC residual magnitude maps +% + added image-space error maps (reference image recovery) +% + added control scheme: pointwise optimal control values with +% local search, local acceleration, and two-scale iteration +% . changed synthetic reference image to smooth version +% . parameter clean-up and explicit visualization options +% % 0.1 (Oct 08, 2018) - Alexandros -% * initial implementation +% . initial implementation % % ------------------------------------------------------------ diff --git a/demo_inversion_3d.m b/demo_inversion_3d.m deleted file mode 100644 index 03cd4a1..0000000 --- a/demo_inversion_3d.m +++ /dev/null @@ -1,245 +0,0 @@ -% demo_inversion_3d.m -% -% 3D DVF inversion [1,2] demo script (3D-expanded 2D DVF). -% -% ---------------------------------------------------------------------- -% -% Copyright (C) 2018, Department of Computer Science, Duke University -% -% This program is free software: you can redistribute it and/or modify it -% under the terms of the GNU General Public License as published by the -% Free Software Foundation, either version 3 of the License, or (at your -% option) any later version. -% -% This program is distributed in the hope that it will be useful, but -% WITHOUT ANY WARRANTY; without even the implied warranty of -% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General -% Public License for more details. -% -% You should have received a copy of the GNU General Public License along -% with this program. If not, see . -% -% ---------------------------------------------------------------------- -% -% REFERENCES -% -% [1] A. Dubey*, A.-S. Iliopoulos*, X. Sun, F.-F. Yin, and L. Ren, -% "Iterative inversion of deformation vector fields with feedback -% control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, 2018. -% DOI: 10.1002/mp.12962. -% -% [2] A. Dubey, "Symmetric completion of deformable registration via -% bi-residual inversion," PhD thesis, Duke University, Durham, NC, USA. -% - - - -%% ==================== CLEAN UP - -clear variables -close all - - -%% ==================== PARAMETERS - -% example DVF data -pathData = 'data/c-deformation.mat'; -ioData = matfile( pathData ); -F = ioData.F; -[szDom, ~] = dvf.sizeVf( F ); - -% 2D-to-3D DVF embedding (illustration with zero Z-axis displacement) -zdim = 15; -F = repmat( F, [1 1 1 zdim] ); % (NX x NY x 2 x NZ) -F = permute( F, [1 2 4 3] ); % (NX x NY x NZ x 2) -F = cat( 4, F, zeros( [szDom zdim] ) ); % (NX x NY x NZ x 3) - -% inversion parameters -opt.control = 'midrange'; -opt.scale = 1; -opt.numIteration = 20; -opt.stopThreshold = 1e-6; -opt.accThreshold = 0; -opt.InitialValue = []; -opt.Mask = dvf.maskDomain( F ); - -% create an array of inversion parameters to compare multiple results -opt = repmat( opt, [7 1] ); -opt(1).control = 0.0; % constant mu = 0 -opt(2).control = 0.5; % constant mu = 0.5 -opt(3).control = 'alternating'; % alternating -opt(4).control = 'optimal'; % locally optimal -opt(5).control = 'midrange'; % local midrange -opt(6).control = 'midrange'; % local midrange & acc. -opt(6).accThreshold = 1; % . -opt(7).control = 'midrange'; % local midrange & acc. & multiscale -opt(7).accThreshold = 1; % . -opt(7).scale = [0.5, 1.0]; % . -opt(7).numIteration = [4 16]; % . - -% legend labels for each control setting -lgnd = {'\mu = 0', '\mu = 0.5', '\mu_{oe}', '\mu_*(x)', '\mu_m(x)', ... - '\mu_m(x) & acc', '\mu_m(x) & acc & MS'}; - -% IC residual percentile curve visualization options -prcIcMeasure = [50 90 95 98 100]; -argVisIcCurve = {'-o', 'LineWidth', 2}; - - -%% ==================== (BEGIN) - -fprintf( '\n***** BEGIN (%s) *****\n\n', mfilename ); - - -%% ==================== VISUALIZATION: REFERENCE & DEFORMED IMAGES - -fprintf( '...visualizing synthetic reference and deformed images...\n' ); - -IRef = util.imageGrid( [szDom zdim] ); -IDfm = dvf.imdeform( IRef, F ); - -figure - -% reference image -subplot( 1, 2, 1 ) -imagesc( IRef(:,:,(zdim+1)/2).' ) -axis image -colormap( gray ) -title( {'reference image', '(central axial slice)'} ) -drawnow - -% deformed (study) image -subplot( 1, 2, 2 ) -imagesc( IDfm(:,:,(zdim+1)/2).' ) -axis image -colormap( gray ) -title( {'deformed image', '(central axial slice)'} ) -drawnow - - -%% ==================== VISUALIZATION: SPECTRAL NTDC MEASURES - -fprintf( '...calculating & visualizing spectral NTDC measures...\n' ); - -[CtrlIdx, Rho, Det, Lambda] = dvf.ntdcMeasures( F ); - -figure - -% control index -subplot( 1, 3, 1 ) -imagesc( CtrlIdx(:,:,(zdim+1)/2).' ) -axis image -colormap( jet ) -colorbar -title( {'algebraic control index', '(central axial slice)'} ); -drawnow - -% NTDC spectral radius -subplot( 1, 3, 2 ) -imagesc( Rho(:,:,(zdim+1)/2).' ) -axis image -colormap( jet ) -colorbar -title( {'NTDC spectral radius', '(central axial slice)'} ); -drawnow - -% determinant map -subplot( 1, 3, 3 ) -imagesc( Det(:,:,(zdim+1)/2).' ) -axis image -colormap( jet ) -colorbar -title( {'determinant map', '(central axial slice)'} ); -drawnow - - -%% ==================== INVERSION - -fprintf( '...computing inverse with %d iteration schemes...\n', length(opt) ); - -G = cell( size(opt) ); -prctRG = cell( size(opt) ); -prctRF = cell( size(opt) ); - -for i = 1 : length(opt) - fprintf( ' - scheme #%d (%s)...\n', i, lgnd{i} ); - opt(i).control = dvf.feedbackControlVal( Lambda, opt(i).control ); - [G{i}, ~, ~, prctRG{i}, prctRF{i}] = dvf.inversion( F, opt(i) ); -end - - -%% ==================== VISUALIZATION: IC RESIDUALS PER ITERATION - -fprintf( '...visualizing IC residual & error percentile curves...\n' ); - -% resolution-normalized iteration steps -kstep = cell( size(opt) ); -for i = 1 : length(opt) - kstep{i} = arrayfun( @(s,k) repmat( s^2 + (s~=1)/k, [1 k] ), ... - opt(i).scale, opt(i).numIteration, ... - 'UniformOutput', false ); - kstep{i} = cumsum( horzcat( 0, kstep{i}{:} ) ); -end - -figure -numPrct = length( prcIcMeasure ); - -% legend -subplot( 5, numPrct, (1:numPrct) ) -plot( nan( 2, length(opt) ), argVisIcCurve{:} ); -title( 'step-wise IC residual percentile curves' ) -legend( lgnd, 'Orientation', 'horizontal', 'Location', 'south' ) -axis( gca, 'off' ) -drawnow - -for p = 1 : numPrct % ---- each percentile - - % IC residual (RG) - subplot( 5, numPrct, p + (1:2)*numPrct ) - hold on - for i = 1 : length(opt) - plot( kstep{i}, prctRG{i}(p,:), argVisIcCurve{:} ); - end - ylabel( sprintf( 'r_G[%d%%]', prcIcMeasure(p) ) ); - xlabel( 'amortizerd iteration step' ) - set( gca, 'YScale', 'log' ) - grid on - drawnow - - % IC residual (RF) - subplot( 5, numPrct, p + (3:4)*numPrct ) - hold on - for i = 1 : length(opt) - plot( kstep{i}, prctRF{i}(p,:), argVisIcCurve{:} ); - end - ylabel( sprintf( 'r_F[%d%%]', prcIcMeasure(p) ) ); - xlabel( 'amortized iteration step' ) - set( gca, 'YScale', 'log' ) - grid on - drawnow - -end % for (p) - - -%% ==================== (END) - -fprintf( '\n***** END (%s) *****\n\n', mfilename ); - - - -%%------------------------------------------------------------ -% -% AUTHORS -% -% Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu -% -% VERSION -% -% 0.1 - October 08, 2018 -% -% CHANGELOG -% -% 0.1 (Oct 08, 2018) - Alexandros -% * initial implementation -% -% ------------------------------------------------------------ diff --git a/demo_inversion_3d_z0.m b/demo_inversion_3d_z0.m new file mode 100644 index 0000000..0de5ad9 --- /dev/null +++ b/demo_inversion_3d_z0.m @@ -0,0 +1,359 @@ +% demo_inversion_3d_z0.m +% +% 3D DVF inversion [1,2] demo script (3D-expanded 2D DVF; no displacement +% along Z-axis). +% +% ---------------------------------------------------------------------- +% +% Copyright (C) 2018, Department of Computer Science, Duke University +% +% This program is free software: you can redistribute it and/or modify it +% under the terms of the GNU General Public License as published by the +% Free Software Foundation, either version 3 of the License, or (at your +% option) any later version. +% +% This program is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +% Public License for more details. +% +% You should have received a copy of the GNU General Public License along +% with this program. If not, see . +% +% ---------------------------------------------------------------------- +% +% REFERENCES +% +% [1] A. Dubey*, A.-S. Iliopoulos*, X. Sun, F.-F. Yin, and L. Ren, +% "Iterative inversion of deformation vector fields with feedback +% control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, 2018. +% - DOI: 10.1002/mp.12962 +% - arXiv: 1610.08589 [cs.CV] +% +% [2] A. Dubey, "Symmetric completion of deformable registration via +% bi-residual inversion," PhD thesis, Duke University, Durham, NC, USA, +% 2018. +% + + + +%% ==================== CLEAN UP + +clear variables +close all + + +%% ==================== PARAMETERS + +% example DVF data +pathData = 'data/c-deformation.mat'; +ioData = matfile( pathData ); +F = ioData.F; +[szDom, ~] = dvf.sizeVf( F ); + +% 2D-to-3D DVF embedding (illustration with zero Z-axis displacement) +zdim = 15; +F = repmat( F, [1 1 1 zdim] ); % (NX x NY x 2 x NZ) +F = permute( F, [1 2 4 3] ); % (NX x NY x NZ x 2) +F = cat( 4, F, zeros( [szDom zdim] ) ); % (NX x NY x NZ x 3) +szDom = [szDom zdim]; + +% inversion parameters +opt.control = 'midrange'; +opt.scale = 1; +opt.numIteration = 20; +opt.stopThreshold = 1e-6; +opt.accThreshold = 0; +opt.InitialValue = []; +opt.Mask = dvf.maskDomain( F ); + +% create an array of inversion parameters to compare multiple schemes +opt = repmat( opt, [8 1] ); +opt(1).control = 0.0; % constant mu = 0 +opt(2).control = 0.5; % constant mu = 0.5 +opt(3).control = 'alternating'; % alternating +opt(4).control = 'optimal'; % pointwise optimal +opt(5).control = 'midrange'; % local midrange +opt(6).control = 'midrange'; % local midrange & acc. +opt(6).accThreshold = 1; % . +opt(7).control = 'midrange'; % local midrange & acc. & multiscale +opt(7).accThreshold = 1; % . +opt(7).scale = [0.5, 1.0]; % . +opt(7).numIteration = [4 16]; % . +opt(8).control = 'optimal'; % pointwise optimal & acc. & multiscale +opt(8).accThreshold = 1; % . +opt(8).scale = [0.5, 1.0]; % . +opt(8).numIteration = [4 16]; % . + +% legend labels for each control setting +lgnd = {'\mu = 0', '\mu = 0.5', '\mu_{oe}', '\mu_*(x)', '\mu_m(x)', ... + '\mu_m(x) & acc', '\mu_m(x) & acc & MS', '\mu_*(x) & acc & MS'}; + +% IC residual percentile curve visualization options +prcIcMeasure = [50 90 95 98 100]; +argVisIcCurve = {'-o', 'LineWidth', 2}; +clrIcCurve = [0.8650 0.8110 0.4330; % yellow + 0.9718 0.5553 0.7741; % pink + 0.6859 0.4035 0.2412; % brown + 1.0000 0.5482 0.1000; % orange + 0.6365 0.3753 0.6753; % purple + 0.3718 0.7176 0.3612; % green + 0.2941 0.5447 0.7494; % blue + 0.9047 0.1918 0.1988]; % red + +% IC residual maps visualization options +zIdx = (zdim + 1) / 2; % axial slice index +climIcMgn = [0 1]; % color-axis limits + + +%% ==================== (BEGIN) + +fprintf( '\n***** BEGIN (%s) *****\n\n', mfilename ); + + +%% ==================== VISUALIZATION: REFERENCE & STUDY IMAGES + +fprintf( '...visualizing synthetic reference and study images...\n' ); + +IRef = util.imageGridSmooth( [szDom zdim], [15 15 5] ); +IStd = dvf.imdeform( IRef, F ); + +hFig = vis.mfigure; +hFig.Name = 'reference and study (deformed) images'; + +% reference image +subplot( 1, 2, 1 ) +imagesc( IRef(:,:,zIdx).' ) +axis image +colormap( gray ) +title( {'reference image', sprintf( '(axial slice z = %d)', zIdx )} ) +drawnow + +% study (deformed) image +subplot( 1, 2, 2 ) +imagesc( IStd(:,:,zIdx).' ) +axis image +colormap( gray ) +title( {'study image', sprintf( '(axial slice z = %d)', zIdx )} ) +drawnow + + +%% ==================== VISUALIZATION: SPECTRAL NTDC MEASURES + +fprintf( '...calculating & visualizing spectral NTDC measures...\n' ); + +[CtrlIdx, Rho, Det, Lambda] = dvf.ntdcMeasures( F ); + +hFig = vis.mfigure; +hFig.Name = 'spectral NTDC measures of forward DVF'; + +% control index +subplot( 1, 3, 1 ) +imagesc( CtrlIdx(:,:,zIdx).' ) +axis image +colormap( jet ) +colorbar +title( {'algebraic control index', sprintf( '(axial slice z = %d)', zIdx )} ); +drawnow + +% NTDC spectral radius +subplot( 1, 3, 2 ) +imagesc( Rho(:,:,zIdx).' ) +axis image +colormap( jet ) +colorbar +title( {'NTDC spectral radius', sprintf( '(axial slice z = %d)', zIdx )} ); +drawnow + +% determinant map +subplot( 1, 3, 3 ) +imagesc( Det(:,:,zIdx).' ) +axis image +colormap( jet ) +colorbar +title( {'determinant map', sprintf( '(axial slice z = %d)', zIdx )} ); +drawnow + + +%% ==================== INVERSION + +fprintf( '...computing inverse with %d iteration schemes...\n', length(opt) ); + +G = cell( size(opt) ); +prctRG = cell( size(opt) ); +prctRF = cell( size(opt) ); +RG = cell( size(opt) ); +RF = cell( size(opt) ); + +for i = 1 : length(opt) + fprintf( ' - scheme #%d (%s)...\n', i, lgnd{i} ); + opt(i).control = dvf.feedbackControlVal( Lambda, opt(i).control ); + [G{i}, ~, ~, prctRG{i}, prctRF{i}] = dvf.inversion( F, opt(i) ); + RG{i} = dvf.icResidual( G{i}, F ); + RF{i} = dvf.icResidual( F, G{i} ); +end + + +%% ==================== VISUALIZATION: IC RESIDUALS PER ITERATION + +fprintf( '...visualizing step-wise IC residual percentile curves...\n' ); + +% resolution-normalized iteration steps +kstep = cell( size(opt) ); +for i = 1 : length(opt) + kstep{i} = arrayfun( @(s,k) repmat( s^2 + (s~=1)/k, [1 k] ), ... + opt(i).scale, opt(i).numIteration, ... + 'UniformOutput', false ); + kstep{i} = cumsum( horzcat( 0, kstep{i}{:} ) ); +end + +hFig = vis.mfigure; +hFig.Name = 'step-wise IC residual percentile curves'; +numPrct = length( prcIcMeasure ); + +% legend +subplot( 5, numPrct, (1:numPrct), ... + 'NextPlot', 'ReplaceChildren', 'ColorOrder', clrIcCurve ) +plot( nan( 2, length(opt) ), argVisIcCurve{:} ); +title( 'step-wise IC residual percentile curves' ) +legend( lgnd, 'Orientation', 'horizontal', 'Location', 'south' ) +axis( gca, 'off' ) +drawnow + +for p = 1 : numPrct % ---- each percentile + + % study IC residuals + subplot( 5, numPrct, p + (1:2)*numPrct ) + hold on + for i = 1 : length(opt) + plot( kstep{i}, prctRG{i}(p,:), argVisIcCurve{:}, ... + 'Color', clrIcCurve(i,:) ); + end + ylabel( sprintf( 's_k[%d%%]', prcIcMeasure(p) ) ); + xlabel( 'amortizerd iteration step (k)' ) + set( gca, 'YScale', 'log' ) + grid on + drawnow + + % reference IC residuals + subplot( 5, numPrct, p + (3:4)*numPrct ) + hold on + for i = 1 : length(opt) + plot( kstep{i}, prctRF{i}(p,:), argVisIcCurve{:}, ... + 'Color', clrIcCurve(i,:) ); + end + ylabel( sprintf( 'r_k[%d%%]', prcIcMeasure(p) ) ); + xlabel( 'amortized iteration step (k)' ) + set( gca, 'YScale', 'log' ) + grid on + drawnow + +end % for (p) + + +%% ==================== VISUALIZATION: IC RESIDUAL MAPS + +fprintf( '...visualizing final IC residual maps...\n' ); +fprintf( ' - white pixels indicate NaN (out-of-bounds) values\n' ); + +% calculate point-wise IC residual magnitudes +RGMgn = cellfun( @(r) sqrt( sum( r.^2, 4 ) ), RG, 'UniformOutput', false ); +RFMgn = cellfun( @(r) sqrt( sum( r.^2, 4 ) ), RF, 'UniformOutput', false ); + +% study IC residuals +hFig = vis.mfigure; +hFig.Name = 'study IC residual maps'; +for i = 1 : length(opt) + subplot( 2, 4, i ) + pcolor( RGMgn{i}(:,:,zIdx).' ) + axis image + shading flat + colormap jet + caxis( climIcMgn ) + colorbar + title( {sprintf( 'study IC residuals (%s)', lgnd{i} ), ... + sprintf( '(axial slice z = %d)', zIdx )} ) + drawnow +end + +% reference IC residuals +hFig = vis.mfigure; +hFig.Name = 'reference IC residual maps'; +for i = 1 : length(opt) + subplot( 2, 4, i ) + pcolor( RFMgn{i}(:,:,zIdx).' ) + axis image + shading flat + colormap jet + caxis( climIcMgn ) + colorbar + title( {sprintf( 'reference IC residuals (%s)', lgnd{i} ), ... + sprintf( '(axial slice z = %d)', zIdx )} ) + drawnow +end + + +%% ==================== VISUALIZATION: RECOVERED REFERENCE IMAGE + +fprintf( '...reference image recovery...\n' ); +fprintf( ' - I_rec(x) := I_std(x + v(x))\t' ); +fprintf( ' [I_std(y) = I_ref(y + u(y))]\n' ); +fprintf( ' = I_ref(x + v(x) + u(x + v(x)))\n' ); +fprintf( ' = I_ref(x + s(x))\n' ); + +% ---------- calculation +fprintf( ' - calculation...\n' ); + +IRefRec = cell( size(opt) ); +for i = 1 : length(opt) + IRefRec{i} = dvf.imdeform( IRef, RG{i} ); +end + +% ---------- visualization +fprintf( ' - difference image visualization (I_ref - I_rec)...\n' ); + +hFig = vis.mfigure; +hFig.Name = 'image-space errors in reference image recovery'; +for i = 1 : length(opt) + subplot( 2, 4, i ) + pcolor( (IRef(:,:,zIdx) - IRefRec{i}(:,:,zIdx)).' ); + axis image + shading flat + caxis( [-1, +1] ) + colormap parula + colorbar + title( {lgnd{i}, 'I_{ref}(x) - I_{ref}(x+s(x))'} ) + drawnow +end + + +%% ==================== (END) + +fprintf( '\n***** END (%s) *****\n\n', mfilename ); + + + +%%------------------------------------------------------------ +% +% AUTHORS +% +% Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu +% +% VERSION +% +% 0.2 - December 21, 2018 +% +% CHANGELOG +% +% 0.2 (Dec 21, 2018) - Alexandros +% + added IC residual magnitude maps +% + added image-space error maps (reference image recovery) +% + added control scheme: pointwise optimal control values with +% local search, local acceleration, and two-scale iteration +% . changed synthetic reference image to smooth version +% . parameter clean-up and explicit visualization options +% +% 0.1 (Oct 08, 2018) - Alexandros +% . initial implementation +% +% ------------------------------------------------------------ diff --git a/demo_inversion_3d_zsin.m b/demo_inversion_3d_zsin.m new file mode 100644 index 0000000..ddae491 --- /dev/null +++ b/demo_inversion_3d_zsin.m @@ -0,0 +1,353 @@ +% demo_inversion_3d_zsin.m +% +% 3D respiratory DVF inversion [1,2] demo script (3D-expanded 2D DVF; +% Z-axis deformation with sinusoidal profile). +% +% ---------------------------------------------------------------------- +% +% Copyright (C) 2018, Department of Computer Science, Duke University +% +% This program is free software: you can redistribute it and/or modify it +% under the terms of the GNU General Public License as published by the +% Free Software Foundation, either version 3 of the License, or (at your +% option) any later version. +% +% This program is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +% Public License for more details. +% +% You should have received a copy of the GNU General Public License along +% with this program. If not, see . +% +% ---------------------------------------------------------------------- +% +% REFERENCES +% +% [1] A. Dubey*, A.-S. Iliopoulos*, X. Sun, F.-F. Yin, and L. Ren, +% "Iterative inversion of deformation vector fields with feedback +% control," Medical Physics, vol. 45, no. 7, pp. 3147-3160, 2018. +% - DOI: 10.1002/mp.12962 +% - arXiv: 1610.08589 [cs.CV] +% +% [2] A. Dubey, "Symmetric completion of deformable registration via +% bi-residual inversion," PhD thesis, Duke University, Durham, NC, USA, +% 2018. +% + + + +%% ==================== CLEAN UP + +clear variables +close all + + +%% ==================== PARAMETERS + +% example CT and DVF data +pathData = 'data/c-deformation.mat'; +ioData = matfile( pathData ); +F = ioData.F; +[szDom, ~] = dvf.sizeVf( F ); + +% 2D-to-3D DVF embedding +zdim = 15; +F = repmat( F, [1 1 1 zdim] ); % (NX x NY x 2 x NZ) +F = permute( F, [1 2 4 3] ); % (NX x NY x NZ x 2) +F = cat( 4, F, 3*ones( [szDom zdim] ) ); % (NX x NY x NZ x 3) +szDom = [szDom zdim]; + +% sinusoidal displacement scaling along Z-axis +F(:,:,:,3) = F(:,:,:,3) .* ... + reshape( sin( linspace( 0 , pi,size(F,2)) ), 1, [], 1 ) .* ... + reshape(1+cos( linspace(-pi,3*pi,size(F,1)) ), [], 1, 1 ) .* ... + reshape( -cos( linspace( 0 , pi,size(F,3)) ), 1, 1, [] ); + +% inversion parameters +opt.control = 'midrange'; +opt.scale = 1; +opt.numIteration = 20; +opt.stopThreshold = 1e-6; +opt.accThreshold = 0; +opt.InitialValue = []; +opt.Mask = dvf.maskDomain( F ); + +% create an array of inversion parameters to compare multiple results +opt = repmat( opt, [8 1] ); +opt(1).control = 0.0; % constant mu = 0 +opt(2).control = 0.5; % constant mu = 0.5 +opt(3).control = 'alternating'; % alternating +opt(4).control = 'optimal'; % pointwise optimal +opt(5).control = 'midrange'; % local midrange +opt(6).control = 'midrange'; % local midrange & acc. +opt(6).accThreshold = 1; % . +opt(7).control = 'midrange'; % local midrange & acc. & multiscale +opt(7).accThreshold = 1; % . +opt(7).scale = [0.5, 1.0]; % . +opt(7).numIteration = [4 16]; % . +opt(8).control = 'optimal'; % pointwise optimal & acc. & multiscale +opt(8).accThreshold = 1; % . +opt(8).scale = [0.5, 1.0]; % . +opt(8).numIteration = [4 16]; % . + +% legend labels for each control setting +lgnd = {'\mu = 0', '\mu = 0.5', '\mu_{oe}', '\mu_*(x)', '\mu_m(x)', ... + '\mu_m(x) & acc', '\mu_m(x) & acc & MS', '\mu_*(x) & acc & MS'}; + +% slices for image and spectral measure maps visualization +slices = [108 51 10]; %14]; + +% IC residual percentile curve visualization options +prcIcMeasure = [50 90 95 98 100]; +argVisIcCurve = {'-o', 'LineWidth', 2}; +clrIcCurve = [0.8650 0.8110 0.4330; % yellow + 0.9718 0.5553 0.7741; % pink + 0.6859 0.4035 0.2412; % brown + 1.0000 0.5482 0.1000; % orange + 0.6365 0.3753 0.6753; % purple + 0.3718 0.7176 0.3612; % green + 0.2941 0.5447 0.7494; % blue + 0.9047 0.1918 0.1988]; % red + +% IC residual maps visualization options +zIdx = slices(3); +climIcMgn = [0 1]; % color-axis limits + +% image-space visualization options +climImg = [-1 1]; % color-axis limits + + +%% ==================== (BEGIN) + +fprintf( '\n***** BEGIN (%s) *****\n\n', mfilename ); + + +%% ==================== VISUALIZATION: REFERENCE & STUDY IMAGES + +fprintf( '...visualizing synthetic reference and study images...\n' ); + +IRef = util.imageGridSmooth( szDom, [15 15 5] ); +IStd = dvf.imdeform( IRef, F ); + +hFig = vis.mfigure; +hFig.Name = 'reference and study (deformed) images'; + +% reference image +vis.imageSlices( IRef, 'slices', slices, 'colormap', gray, ... + 'title', 'reference image', 'caxis', climImg, ... + 'rows', 2, 'columns', 3, 'startIndex', 1 ); + +% study (deformed) image +vis.imageSlices( IStd, 'slices', slices, 'colormap', gray, ... + 'title', 'study image', 'caxis', climImg, ... + 'rows', 2, 'columns', 3, 'startIndex', 4 ); + + +%% ==================== VISUALIZATION: SPECTRAL NTDC MEASURES + +fprintf( '...calculating & visualizing spectral NTDC measures...\n' ); + +[CtrlIdx, Rho, Det, Lambda] = dvf.ntdcMeasures( F ); + +hFig = vis.mfigure; +hFig.Name = 'spectral NTDC measures of forward DVF'; + +% control index +vis.imageSlices( CtrlIdx, 'slices', slices, 'colormap', jet, ... + 'title', 'algebraic control index', ... + 'rows', 3, 'columns', 3, 'startIndex', 1 ); + +% NTDC spectral radius +vis.imageSlices( Rho, 'slices', slices, 'colormap', jet, ... + 'title', 'NTDC spectral radius', ... + 'rows', 3, 'columns', 3, 'startIndex', 4 ); + +% determinant map +vis.imageSlices( Det, 'slices', slices, 'colormap', jet, ... + 'title', 'determinant map', ... + 'rows', 3, 'columns', 3, 'startIndex', 7 ); + + +%% ==================== INVERSION + +fprintf( '...computing inverse with %d iteration schemes...\n', length(opt) ); + +G = cell( size(opt) ); +prctRG = cell( size(opt) ); +prctRF = cell( size(opt) ); +RG = cell( size(opt) ); +RF = cell( size(opt) ); + +for i = 1 : length(opt) + fprintf( ' - scheme #%d (%s)...\n', i, lgnd{i} ); + opt(i).control = dvf.feedbackControlVal( Lambda, opt(i).control ); + [G{i}, ~, ~, prctRG{i}, prctRF{i}] = dvf.inversion( F, opt(i) ); + RG{i} = dvf.icResidual( G{i}, F ); + RF{i} = dvf.icResidual( F, G{i} ); +end + + +%% ==================== VISUALIZATION: IC RESIDUALS PER ITERATION + +fprintf( '...visualizing step-wise IC residual percentile curves...\n' ); + +% resolution-normalized iteration steps +kstep = cell( size(opt) ); +for i = 1 : length(opt) + kstep{i} = arrayfun( @(s,k) repmat( s^2 + (s~=1)/k, [1 k] ), ... + opt(i).scale, opt(i).numIteration, ... + 'UniformOutput', false ); + kstep{i} = cumsum( horzcat( 0, kstep{i}{:} ) ); +end + +hFig = vis.mfigure; +hFig.Name = 'step-wise IC residual percentile curves'; +numPrct = length( prcIcMeasure ); + +% legend +subplot( 5, numPrct, (1:numPrct), ... + 'NextPlot', 'ReplaceChildren', 'ColorOrder', clrIcCurve ) +plot( nan( 2, length(opt) ), argVisIcCurve{:} ); +title( 'step-wise IC residual percentile curves' ) +legend( lgnd, 'Orientation', 'horizontal', 'Location', 'south' ) +axis( gca, 'off' ) +drawnow + +for p = 1 : numPrct % ---- each percentile + + % srudy IC residuals + subplot( 5, numPrct, p + (1:2)*numPrct ) + hold on + for i = 1 : length(opt) + plot( kstep{i}, prctRG{i}(p,:), argVisIcCurve{:}, ... + 'Color', clrIcCurve(i,:) ); + end + ylabel( sprintf( 's_k[%d%%]', prcIcMeasure(p) ) ); + xlabel( 'amortizerd iteration step (k)' ) + set( gca, 'YScale', 'log' ) + grid on + drawnow + + % reference IC residuals + subplot( 5, numPrct, p + (3:4)*numPrct ) + hold on + for i = 1 : length(opt) + plot( kstep{i}, prctRF{i}(p,:), argVisIcCurve{:}, ... + 'Color', clrIcCurve(i,:) ); + end + ylabel( sprintf( 'r_k[%d%%]', prcIcMeasure(p) ) ); + xlabel( 'amortized iteration step (k)' ) + set( gca, 'YScale', 'log' ) + grid on + drawnow + +end % for (p) + + +%% ==================== VISUALIZATION: IC RESIDUAL MAPS + +fprintf( '...visualizing final IC residual maps...\n' ); +fprintf( ' - white pixels indicate NaN (out-of-bounds) values\n' ); + +% calculate point-wise IC residual magnitudes +RGMgn = cellfun( @(r) sqrt( sum( r.^2, 4 ) ), RG, 'UniformOutput', false ); +RFMgn = cellfun( @(r) sqrt( sum( r.^2, 4 ) ), RF, 'UniformOutput', false ); + +% study IC residuals +hFig = vis.mfigure; +hFig.Name = 'study IC residual maps'; +for i = 1 : length(opt) + subplot( 2, 4, i ) + pcolor( RGMgn{i}(:,:,zIdx).' ) + axis image + shading flat + colormap jet + caxis( climIcMgn ) + colorbar + title( {sprintf( 'study IC residuals (%s)', lgnd{i} ), ... + sprintf( '(axial slice z = %d)', zIdx )} ) + drawnow +end + +% reference IC residuals +hFig = vis.mfigure; +hFig.Name = 'reference IC residual maps'; +for i = 1 : length(opt) + subplot( 2, 4, i ) + pcolor( RFMgn{i}(:,:,zIdx).' ) + axis image + shading flat + colormap jet + caxis( climIcMgn ) + colorbar + title( {sprintf( 'reference IC residuals (%s)', lgnd{i} ), ... + sprintf( '(axial slice z = %d)', zIdx )} ) + drawnow +end + + +%% ==================== VISUALIZATION: RECOVERED REFERENCE IMAGE + +fprintf( '...reference image recovery...\n' ); +fprintf( ' - I_rec(x) := I_std(x + v(x))\t' ); +fprintf( ' [I_std(y) = I_ref(y + u(y))]\n' ); +fprintf( ' = I_ref(x + v(x) + u(x + v(x)))\n' ); +fprintf( ' = I_ref(x + s(x))\n' ); + +% ---------- calculation +fprintf( ' - calculation...\n' ); + +IRefRec = cell( size(opt) ); +for i = 1 : length(opt) + IRefRec{i} = dvf.imdeform( IRef, RG{i} ); +end + +% ---------- visualization +fprintf( ' - difference image visualization (I_ref - I_rec)...\n' ); + +hFig = vis.mfigure; +hFig.Name = 'image-space errors in reference image recovery'; +for i = 1 : length(opt) + subplot( 2, 4, i ) + pcolor( (IRef(:,:,zIdx) - IRefRec{i}(:,:,zIdx)).' ); + axis image + shading flat + caxis( climImg ) + colormap parula + colorbar + title( {lgnd{i}, 'I_{ref}(x) - I_{ref}(x+s(x))'} ) + drawnow +end + + +%% ==================== (END) + +fprintf( '\n***** END (%s) *****\n\n', mfilename ); + + + +%%------------------------------------------------------------ +% +% AUTHORS +% +% Alexandros-Stavros Iliopoulos ailiop@cs.duke.edu +% +% VERSION +% +% 0.2 - December 21, 2018 +% +% CHANGELOG +% +% 0.2 (Dec 21, 2018) - Alexandros +% + added IC residual magnitude maps +% + added image-space error maps (reference image recovery) +% + added control scheme: pointwise optimal control values with +% local search, local acceleration, and two-scale iteration +% . changed synthetic reference image to smooth version +% . parameter clean-up and explicit visualization options +% +% 0.1 (Oct 08, 2018) - Alexandros +% . initial implementation +% +% ------------------------------------------------------------ diff --git a/paper.md b/paper.md index 4b7bbe7..041f0d8 100644 --- a/paper.md +++ b/paper.md @@ -43,15 +43,15 @@ for DVF inversion with guaranteed convergence. We use an iterative inversion procedure and employ adaptive bi-residual feedback control to achieve global convergence and local acceleration -[@medphys2018; @phddubey2018]. The inverse DVF estimate and the forward DVF -must meet the inverse consistency (IC) condition. We define two IC +[@medphys2018; @phddubey2018]. The inverse DVF estimate and the forward +DVF must meet the inverse consistency (IC) condition. We define two IC residuals, which measure the inconsistency between two DVFs in the study and reference domains. The residuals, which are computationally available, are related to the unknown inversion error. We use them as feedback in a two-phase iteration. In the first phase, we modulate the study-domain IC residual with adaptive feedback control for guaranteed global convergence, under certain mild conditions [@medphys2018]. Once the error is made -sufficiently small we switch to phase two, where we use the +sufficiently small, we switch to phase two where we use the reference-domain IC residual and achieve locally quadratic convergence rate. Phase transition and integration is enabled by a multi-resolution scheme [@phddubey2018]. diff --git a/references.bib b/references.bib index 14e2678..f1ac028 100644 --- a/references.bib +++ b/references.bib @@ -14,8 +14,8 @@ @Article{medphys2018 @PhdThesis{phddubey2018, author = {Dubey, Abhishek}, - title = {Symmetric Completion of Deformable Registration via - Bi-Residual Inversion}, + title = {Symmetric completion of deformable registration via + bi-residual inversion}, school = {Duke University}, year = 2018, address = {Durham, NC, USA}