 function C = Cfd2_abs(nx, ny)
%function C = Cfd2_abs(nx, ny)
%|
%| Create 2D finite-difference object (fatrix2) for regularization.
%| Simple version for just horizontal and vertical differences
%| with periodic boundary conditions.
%| This version supports "abs" needed for separable quadratic surrogates.
%| See Cdiffs and Reg1 for more general regularizers. 
%|
%| in
%|	nx,ny		size of input image
%| out
%|	C		fatrix2 of size [nx ny] by [nx ny 2]
%|
%| if g = C * f then
%|	g[n,m,1] = f[n,m] - f[mod(n-1,N), m]
%|	g[n,m,2] = f[n,m] - f[n, mod(m-1,M)]
%|
%| if g = abs(C) * f then
%|	g[n,m,1] = f[n,m] + f[mod(n-1,N), m]
%|	g[n,m,2] = f[n,m] + f[n, mod(m-1,M)]
%|
%| 2013-11-03 Jeff Fessler, University of Michigan

if nargin == 1 && streq(nx, 'test'), Cfd2_abs_test, return, end
if nargin < 2, help(mfilename), error(mfilename), end

arg.coef = -1; % gets replaced by +1 for abs(C)

% function for computing horiz and vert finite differences
% input x is [nx ny] and output is [nx ny 2]
forw = @(arg,x) cat(3, ...
		x + arg.coef * circshift(x, [1 0]), ...
		x + arg.coef * circshift(x, [0 1]));

% adjoint operation
% input y is [nx ny 2] and output is [nx ny]
back = @(arg,y) ...
		(y(:,:,1) + arg.coef * circshift(y(:,:,1), [-1 0])) + ...
		(y(:,:,2) + arg.coef * circshift(y(:,:,2), [0 -1]));

idim = [nx ny];
odim = [nx ny 2]; % trick: output is 2 finite difference images
C = fatrix2('forw', forw, 'back', back, 'idim', [nx ny], 'odim', odim, ...
	'arg', arg, 'abs', @Cfd2_abs_fun);


function ob = Cfd2_abs_fun(ob)
ob.arg.coef = 1;


function Cfd2_abs_test
nx = 8;
ny = 6;
C = Cfd2_abs(nx, ny);
%fatrix2_tests(C)

C = abs(C);
x = zeros(C.idim); x(nx/2,ny/2) = 1;
y = C * x;
b = C' * y;

im plc 1 3
im(1, x)
im(2, y)
im(3, b)

test_adjoint(C);
