# gradientNeumannAdj

Computes the discretized adjoint of the gradient with Neumann boundary conditions.

## Contents

- Syntax
- Description
- Input Arguments
- Output Arguments
- More About
- References
- See Also
- Code

## Syntax

    res = gradientNeumannAdj(yg,seti)

## Description

`res = gradientNeumannAdj(yg,seti)` computes the adjoint of the gradient using standard finite differences with Neumann boundary conditions.

This function corresponds to gradientNeumann.html.

## Input Arguments

- `yg` is in 2D: 2*dim x nInv x nInv, see gradientNeumann.html.
- `yg` is in 3D: 2*dim x nInv x nInv x nInv, see gradientNeumann.html.
- seti : structural array.
- seti.T, seti.S : identify complex with real, see setIdImagReal.html.
- seti.hInv, seti.nInv : see setGridScale.html.
- seti.iG : write matrix as vector, see setReshapeVecMat.html.

## Output Arguments

- res : resulting adjoint of gradient

## More About

Analytical

Analytical we have the relation

,

i.e. , see [1, Sec. 6.1]$.

Discretization

For discretization we use explicitly in case of 2D:

- (case i = 1):
- (case 1 < i < n):
- (case i = n):

The case is analog.

An extension to 3D is straightforward.

## References

- [1] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1):120-145, 2011. URL: https://doi.org/10.1007/s10851-010-0251-1.

## See Also

- gradientNeumann.html

## Code

    function res = gradientNeumannAdj(yg,seti)

    ygz = seti.T(yg); % ygz: dim x nInv x nInv (x nInv): so GD and complex
    % if seti.dim == 2
    %     x = divergence(squeeze(ygz(1,:,:)),squeeze(ygz(2,:,:)));
    % elseif seti.dim == 3
    %     x = divergence(squeeze(ygz(1,:,:,:)),squeeze(ygz(2,:,:,:)),squeeze(ygz(3,:,:,:)));
    % end
    % % x is a complex matrix of dimension 2 or 3
    % res = seti.iG(x); % write matrix as vector
    % res = seti.S(res); % result: matrix in R^2
    % %res is xnRVD

    h = seti.hInv;
    N = seti.nInv;
    % G = seti.G (vector -> matrix)

    gAdj = gradAdjCalc(ygz,h,N);
    res = seti.iG(gAdj); % write matrix as vector
    res = seti.S(res); % result: matrix in R^2: res is xnRVD

    end

Code: subfunction: guAdjCalc (fast)

    function res = gradAdjCalc(y,h,N)
    % N = nInv
    dim = size(y,1); % y is ygz, so dim x nInv x nInv (x nInv)

    if dim == 2
        gradAdjComp = zeros(dim,N,N);

        % gradAdjComp(1,:,:) = gradAdj(y)_{i,j}^1
        % = (case i = 1):     -y(1,j)/h
        %   (case 1 < i < n): -(y(i,j)-y(i-1,j))/h
        %   (case i = n):     y(n-1,j)/h
        yRowShift = circshift(y,[0 +1 0]);
        gradAdjComp(1,:,:) = squeeze(-(y(1,:,:)-yRowShift(1,:,:))/h); % case 1 < i < n
        gradAdjComp(1,1,:) = squeeze(-y(1,1,:)/h); % case i = 1
        gradAdjComp(1,end,:) = squeeze(y(1,end-1,:)/h); % case i = n

        % gradAdjComp(2,:,:) = gradAdj(y)_{i,j}^2
        yColShift = circshift(y,[0 0 +1]);
        gradAdjComp(2,:,:) = squeeze(-(y(2,:,:)-yColShift(2,:,:))/h); % case 1 < j < n
        gradAdjComp(2,:,1) = squeeze(-y(2,:,1)/h); % case j = 1
        gradAdjComp(2,:,end) = squeeze(y(2,:,end-1)/h); % case j = n

        % gradAdj = gradAdj(y)_{i,j}^1 + gradAdj(y)_{i,j}^2
        gradAdj = gradAdjComp(1,:,:) + gradAdjComp(2,:,:); % size N x N
        res = gradAdj;
    elseif dim == 3
        gradAdjComp = zeros(dim,N,N,N);

        % gradAdjComp(1,:,:,:) = gradAdj(y)_{i,j,k}^1
        y1Shift = circshift(y,[0 +1 0 0]);
        gradAdjComp(1,:,:,:) = squeeze(-(y(1,:,:,:)-y1Shift(1,:,:,:))/h); % case 1 < i < n
        gradAdjComp(1,1,:,:) = squeeze(-y(1,1,:,:)/h); % case i = 1
        gradAdjComp(1,end,:,:) = squeeze(y(1,end-1,:,:)/h); % case i = n

        % gradAdjComp(2,:,:) = gradAdj(y)_{i,j,k}^2
        y2Shift = circshift(y,[0 0 +1 0]);
        gradAdjComp(2,:,:,:) = squeeze(-(y(2,:,:,:)-y2Shift(2,:,:,:))/h); % case 1 < j < n
        gradAdjComp(2,:,1,:) = squeeze(-y(2,:,1,:)/h); % case j = 1
        gradAdjComp(2,:,end,:) = squeeze(y(2,:,end-1,:)/h); % case j = n

        % gradAdjComp(3,:,:) = gradAdj(y)_{i,j,k}^3
        y3Shift = circshift(y,[0 0 0 +1]);
        gradAdjComp(3,:,:,:) = squeeze(-(y(3,:,:,:)-y3Shift(3,:,:,:))/h); % case 1 < k < n
        gradAdjComp(3,:,:,1) = squeeze(-y(3,:,:,1)/h); % case k = 1
        gradAdjComp(3,:,:,end) = squeeze(y(3,:,:,end-1)/h); % case k = n

        % gradAdj = gradAdj(y)_{i,j,k}^1 + gradAdj(y)_{i,j,k}^2 + gradAdj(y)_{i,j,k}^3
        gradAdj = gradAdjComp(1,:,:) + gradAdjComp(2,:,:) + gradAdjComp(3,:,:); % size N x N x N
        res = gradAdj;

    else
        error('gradientNeumannAdj.m: dimension must be 2D or 3D.')
    end
    end