 function C = Cfd2(nx, ny)
%function C = Cfd2(nx, ny)
%|
%| Create 2D finite-difference object (fatrix2) for regularization.
%| Simple version for just horizontal and vertical differences
%| with periodic boundary conditions.
%| 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)]
%|
%| 2013-09-10 Jeff Fessler, University of Michigan

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

% 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 - circshift(x, [1 0]), ...
		x - circshift(x, [0 1]));

% adjoint operation
% input y is [nx ny 2] and output is [nx ny]
back = @(arg,y) ...
		(y(:,:,1) - circshift(y(:,:,1), [-1 0])) + ...
		(y(:,:,2) - 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);


function Cfd_test
nx = 8;
ny = 6;
C = Cfd2(nx, ny);
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);

if 0 % test vs Cdiff1
	Ch = Cdiff1([nx ny], 'type_diff', 'circshift', 'offset', [1 0]);
	Cv = Cdiff1([nx ny], 'type_diff', 'circshift', 'offset', [0 1]);
	Cd = [Ch; Cv];
	jf_equal(full(Cd), full(C))
end
