 function C = Cfd2_stack(nx, ny)
%function C = Cfd2_stack(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)]
%|
%| This version probably has *slightly* more overhead than Cfd2.m
%| because it makes a "meta" fatrix object from two fatrix objects,
%| but the computation time difference is less than 10% on m laptop.
%|
%| 2013-09-10 Jeff Fessler, University of Michigan

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

% object for computing horizontal finite differences
forw = @(arg,x) x - circshift(x, [1 0]);
back = @(arg,x) x - circshift(x, [-1 0]);
C1 = fatrix2('forw', forw, 'back', back, 'idim', [nx ny], 'odim', [nx ny]);

% object for computing vertical finite differences
forw = @(arg,x) x - circshift(x, [0 1]);
back = @(arg,x) x - circshift(x, [0 -1]);
C2 = fatrix2('forw', forw, 'back', back, 'idim', [nx ny], 'odim', [nx ny]);

% stack them up:
C = [C1; C2]; % this works because vertcat() is overloaded


function Cfd2_stack_test
nx = 8;
ny = 6;
C = Cfd2_stack(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);

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

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
