function varargout = spm_shoot3di(v0,prm,args)
% Geodesic shooting
% FORMAT [theta,Jtheta,v1,phi,Jphi] = spm_shoot3di(v0,prm,args)
% v0     - Initial velocity field n1*n2*n3*3 (single prec. float)
% prm  - 8 settings
%        - [1][2][3] Voxel sizes
%        - [4][5][6][7][8] Regularisation settings.
%          Regularisation uses the sum of
%          - [4] - absolute displacements
%          - [5] - laplacian
%          - [6] - bending energy
%          - [7] - linear elasticity mu
%          - [8] - linear elasticity lambda
% args   - Integration parameters
%          - [1] Num time steps
%
% theta  - Inverse deformation field n1*n2*n3*3 (single prec. float)
% Jtheta - Inverse Jacobian tensors n1*n2*n3 (single prec. float)
% v1     - Final velocity field n1*n2*n3*3 (single prec. float)
% phi    - Forward deformation field n1*n2*n3*3 (single prec. float)
% Jphi   - Forward Jacobian tensors n1*n2*n3 (single prec. float)
%
% This code generates deformations and their Jacobian determinans from
% initial velocity fields by gedesic shooting.  See the work of Miller,
% Younes and others.
%
% LDDMM (Beg et al) uses the following evolution equation:
%     d\phi/dt = v_t(\phi_t)
% where a variational procedure is used to find the stationary solution
% for the time varying velocity field.
% In principle though, once the initial velocity is known, then the
% velocity at subsequent time points can be computed.  This requires
% initial momentum (m_0), computed (using differential operator L) by:
%     m_0 = L v_0
% Then (Ad_{\phi_t})^* m_0 is computed:
%     m_t = |d \phi_t| (d\phi_t)^T m_0(\phi_t)
% The velocity field at this time point is then obtained by using
% multigrid to solve:
%     v_t = L^{-1} m_t
%
% These equations can be found in:
% Younes (2007). "Jacobi fields in groups of diffeomorphisms and
% applications". Quarterly of Applied Mathematics, vol LXV,
% number 1, pages 113-134 (2007).
%
% Multigrid is currently used to obtain v_t = L^{-1} m_t, but
% this could also be done by convolution with the Greens function
% N = L^{-1} (see e.g. Bro-Nielson).
%
%________________________________________________________
% (c) Wellcome Trust Centre for NeuroImaging (2009)

% John Ashburner
% $Id: spm_shoot3di.m 7461 2018-10-29 15:59:58Z john $

args0 = [8 4 4];
if nargin<3
    args = args0;
else
    if numel(args)<numel(args0)
        args = [args args0((numel(args)+1):end)];
    end
end
verb     = false;
N        = args(1);   % # Time steps
d        = size(v0);
d        = d(1:3);
vt       = v0;

if ~isfinite(N)
    % Number of time steps from an educated guess about how far to move
    N = double(floor(sqrt(max(max(max(v0(:,:,:,1).^2+v0(:,:,:,2).^2+v0(:,:,:,3).^2)))))+1);
end

F = spm_shoot_greens('kernel',d,prm);

if verb, fprintf('N=%g:', N); end

m0    = spm_diffeo('vel2mom',v0,prm); % Initial momentum (m_0 = L v_0)
if verb, fprintf('\t%g', 0.5*sum(v0(:).*m0(:))/prod(d)); end

% Compute initial small deformation and Jacobian matrices from the velocity.
% The overall deformation and its Jacobian tensor field is generated by
% composing a series of small deformations.
[theta,Jtheta] = spm_diffeo('smalldef', vt,-1/N);

% If required, also compute the forward and possibly also its Jacobian
% tensor field. Note that the order of the compositions will be the
% reverse of those for building the forward deformation.
if nargout>=5
    [ phi, Jphi]  = spm_diffeo('smalldef', vt,1/N);
elseif nargout>=4
    phi           = spm_diffeo('smalldef', vt,1/N);
end

for t=2:abs(N)
    mt        = pullg(m0,theta,Jtheta);
    %vt       = spm_diffeo('mom2vel', mt, [prm 2 2],vt); Solve via V-cycles
    vt        = spm_shoot_greens(mt,F,prm);
    if verb, fprintf('\t%g', 0.5*sum(vt(:).*mt(:))/prod(d)); end

    [   dp,    dJ] = spm_diffeo('smalldef',  vt,-1/N);    % Small deformation
    [theta,Jtheta] = spm_diffeo('comp', theta,dp, Jtheta,dJ); % Build up large def. from small defs
    clear dp dJ

    % If required, build up forward deformation and its Jacobian tensor field from
    % small deformations
    if nargout>=5
        [ dp,dJ]    = spm_diffeo('smalldef',  vt,1/N);    % Small deformation
        [phi, Jphi] = spm_diffeo('comp', dp, phi, dJ, Jphi); % Build up large def. from small defs
        clear dp dJ
    elseif nargout>=4
        dp          = spm_diffeo('smalldef',  vt,1/N);    % Small deformation
        phi         = spm_diffeo('comp', dp, phi);        % Build up large def. from small defs
        clear dp
    end
    drawnow
end
if verb, fprintf('\n'); end


varargout{1} = theta;
varargout{2} = Jtheta;
if nargout>=3
    mt           = pullg(m0,theta,Jtheta);
    varargout{3} = spm_shoot_greens(mt,F,prm);
end
if nargout>=4, varargout{4} = phi;  end
if nargout>=5, varargout{5} = Jphi; end
%__________________________________________________________________________________

%__________________________________________________________________________________
function mt = pullg(m0,phi,J)
% The Ad* operation
mt = zeros(size(m0),'single');
for i=1:size(m0,3)
    phii        = phi(:,:,i,:);
    mr1         = spm_diffeo('pullc',m0(:,:,:,1), phii);
    mr2         = spm_diffeo('pullc',m0(:,:,:,2), phii);
    mr3         = spm_diffeo('pullc',m0(:,:,:,3), phii);

    dj          = J(:,:,i,1,1).*(J(:,:,i,2,2).*J(:,:,i,3,3)-J(:,:,i,2,3).*J(:,:,i,3,2))...
                 +J(:,:,i,2,1).*(J(:,:,i,1,3).*J(:,:,i,3,2)-J(:,:,i,1,2).*J(:,:,i,3,3))...
                 +J(:,:,i,3,1).*(J(:,:,i,1,2).*J(:,:,i,2,3)-J(:,:,i,1,3).*J(:,:,i,2,2));

    mt(:,:,i,1) = (J(:,:,i,1,1).*mr1 + J(:,:,i,2,1).*mr2 + J(:,:,i,3,1).*mr3).*dj;
    mt(:,:,i,2) = (J(:,:,i,1,2).*mr1 + J(:,:,i,2,2).*mr2 + J(:,:,i,3,2).*mr3).*dj;
    mt(:,:,i,3) = (J(:,:,i,1,3).*mr1 + J(:,:,i,2,3).*mr2 + J(:,:,i,3,3).*mr3).*dj;
end
%__________________________________________________________________________________

%__________________________________________________________________________________

