 function C = Cfd2_nx2(nx, ny)
%function C = Cfd2_nx2(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 [2*nx ny]
%|
%| 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)]
%|
%| This version uses one fatrix2 but output dimension (odim) differs
%| from that of Cfd2.m in a way that seems somewhat less convenient.
%| I recommend using Cfd2 or Cfd2_stack instead.
%|
%| 2013-09-16 Jeff Fessler, University of Michigan

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

% object for computing both horizontal and vertical finite differences
forw = @(arg,x) [ ...
	x - circshift(x, [1 0]); % this semicolon concatates along x
	x - circshift(x, [0 1])];
back = @(arg,y) ...
	y(1:end/2,:) - circshift(y(1:end/2,:), [-1 0]) + ...
	y(end/2+1:end,:) - circshift(y(end/2+1:end,:), [0 -1]);
C = fatrix2('forw', forw, 'back', back, 'idim', [nx ny], 'odim', [2*nx ny]);


function Cfd2_nx2_test
nx = 8;
ny = 6;
C = Cfd2_nx2(nx, ny);
x = zeros(C.idim); x(nx/2,ny/2) = 1;
y = C * x;
b = C' * y;

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

test_adjoint(C);
prompt

if 1 % 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];
	Cd = full(Cd);
	C = full(C);
	% the next 3 lines are inconvenient permutations
	C = reshape(C, [nx 2 ny nx*ny]);
	C = permute(C, [1 3 2 4]);
	C = reshape(C, [nx*2*ny nx*ny]);
	im plc 1 2, im(1, C'), im(2, Cd')
	jf_equal(full(Cd), full(C))
end

if 0 % timing tests
	n = 2^12;
	C0 = Cfd2(n,n);
	Cs = Cfd2_stack(n,n);
	tmp = ones(n,n);
	tic, C0 * tmp; toc
	tic, Cs * tmp; toc
end
