% Returns tv(U) and a subgradient G at U.
function [tv, G] = tv_l2_subgrad(U)
    [m, n] = size(U);
    % The intensity differences across rows.
    row_diff = U(1:m-1, 2:n) - U(1:m-1, 1:n-1);
    % The intensity differences across columns.
    col_diff = U(2:m, 1:n-1) - U(1:m-1, 1:n-1);
    % Compute the total variation.
    grad_norms = norms([row_diff(:) col_diff(:)], 2, 2);
    tv = sum(grad_norms)/((m-1)*(n-1));
    % Find a subgradient.
    G = zeros(m, n);
    % When non-differentiable, set to 0.
    grad_mat = reshape(grad_norms, m-1, n-1);
    grad_mat(grad_mat == 0) = inf;
    G(1:m-1, 1:n-1) = G(1:m-1, 1:n-1) - row_diff./grad_mat;
    G(1:m-1, 2:n) = G(1:m-1, 2:n) + row_diff./grad_mat;
    G(1:m-1, 1:n-1) = G(1:m-1, 1:n-1) - col_diff./grad_mat;
    G(2:m, 1:n-1) = G(2:m, 1:n-1) + col_diff./grad_mat;
end