Function discrete_gradient_into 

fn discrete_gradient_into(
    x: &[f64],
    out: &mut [f64],
    grid: &Grid,
    axis: usize,
    transpose: bool,
)

Backward difference along axis, scaled by voxel size, zero on the leading face.

transpose gives the operator upstream pairs with it: a forward difference, zero on the trailing face. The two agree with a true adjoint pair except on the two boundary faces, which carry no TV weight (the brain mask is zero there), so the discrepancy never reaches the solution.