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

@@ -74,7 +74,7 @@
 %
 %   Alexandros-Stavros Iliopoulos       [email protected]
 %
-% VERSION
+% RELEASE
 %
 %   1.0.0 - October 31, 2018
 %

+dvf/coorDisplace.m

@@ -124,7 +124,7 @@
 %
 %   Alexandros-Stavros Iliopoulos       [email protected]
 %
-% VERSION
+% RELEASE
 %
 %   1.0.0 - October 31, 2018
 %

+dvf/eigJacobian.m

@@ -102,7 +102,7 @@
 %
 %   Alexandros-Stavros Iliopoulos       [email protected]
 %
-% VERSION
+% RELEASE
 %
 %   1.0.0 - October 31, 2018
 %

+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                      [email protected]
 %   Xiaobai Sun                         [email protected]
 %
-% 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
 %
 % ------------------------------------------------------------

+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       [email protected]
 %
-% VERSION
+% RELEASE
 %
 %   1.0.0 - October 31, 2018
 %

+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       [email protected]
 %
-% 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
 %
 % ------------------------------------------------------------

+dvf/interpnVf.m

@@ -113,7 +113,7 @@
 %
 %   Alexandros-Stavros Iliopoulos       [email protected]
 %
-% VERSION
+% RELEASE
 %
 %   1.0.0 - October 31, 2018
 %