  function ob = GFD(varargin)
%|function ob = GFD([args])
%|
%| Construct object that computes image gradient for regularization based on sparsity such as in compressed sensing.
%|
%| option
%|	'mask'	logical [nx ny]	image-domain mask, often: true(nx,ny)
%|
%| out:
%|	ob	[np]	Fatrix object, where np = sum(mask(:))
%|
%| Sathish Ramani, University of Michigan, 2011-04-14

if nargin < 1, help(mfilename), error(mfilename), end
if nargin == 1 && streq(varargin{1}, 'test'), GFD_test, return, end

% options
arg.mask = [];

% finite-difference filters in 4 directions in 2D
fz = zeros( 3 );
fz( 2, 2 ) = -1;
fdfilters = zeros( numel(fz), 4 );
for ii = 6 : 9
    fi = fz; fi = fi(:); fi(ii) = 1; 
    arg.fdfilters( :, :, ii - 5 ) = reshape( fi, [ 3 3 ] );
end
clear fz fi ii;
arg.nf = 4; % FD in 4 directions

arg = vararg_pair(arg, varargin);

[nx ny] = size( arg.mask );
arg.odim = [nx ny arg.nf];

if isempty(arg.mask), fail 'must provide a mask', end
if ~islogical(arg.mask), error 'mask must be logical', end

% build fatrix2 object
ob = fatrix2('imask', arg.mask, 'odim', arg.odim, 'arg', arg, ...
	'forw', @GFD_forw, 'back', @GFD_back);


% GFD_forw(): y = A * x 
% image gradient
function y = GFD_forw( arg, x )
y = zeros( [ size( arg.mask ) arg.nf ] );
for ii = 1 : arg.nf
    y( :, :, ii ) = imfilter( x, arg.fdfilters( :, :, ii ), 'circular' );
end


% GFD_back(): x = A' * y
% adjoint of image gradient
function x = GFD_back( arg, y )
x = 0;
for ii = 1 : arg.nf
    x = x + imfilter( y(:, :, ii ), arg.fdfilters( end : -1 : 1, end : -1 : 1, ii ), 'circular' );
end


% Gwave2_test()
function GFD_test

mask = true(2^5,2^4); mask(1) = false; % stress
C = GFD('mask', mask);

x = zeros(size(mask)); x(end/2+1,floor(end/2+1)) = 1;
tmp = C * x;
im(tmp), cbar

% todo: test complex case too

if 1
	fatrix2_tests(C)
	test_adjoint(C, 'big', 1)
end