function Gtv = make_TV_fatrix(varargin)

arg.class = 'fatrix2';
arg.im_dims = [];
arg.tv_dims = [];

arg = vararg_pair(arg, varargin);

ndim = length(arg.im_dims);
if isempty(arg.tv_dims), arg.tv_dims = 1:min(3,ndim); end
arg.tv_dims = unique(arg.tv_dims,'sorted');

arg.tv_ffts = arrayfun(@(x) fft([1,-1],x),arg.im_dims(arg.tv_dims),'UniformOutput',false);

arg.idim = [prod(arg.im_dims),1];
arg.odim = [prod(arg.im_dims)*length(arg.tv_dims),1];

Gtv = fatrix2('arg',arg,'forw',@TV_forw,'back',@TV_back);

end

function y = TV_forw(arg,x)

x = reshape(x,arg.im_dims);

dims1 = arrayfun(@(x) 1:x,arg.im_dims,'UniformOutput',false);
y = zeros([prod(arg.im_dims),length(arg.tv_dims)],class(x));
for id = 1:length(arg.tv_dims)
    dims2 = dims1;
    dims2{arg.tv_dims(id)} = circshift(dims2{arg.tv_dims(id)},[0,-1]);
    y(:,id) = col(x(dims2{:}) - x(dims1{:}));
end
y = col(y);

end

function x = TV_back(arg,y)

y = reshape(y,[prod(arg.im_dims),length(arg.tv_dims)]);

x = 0;

dims1 = arrayfun(@(x) 1:x,arg.im_dims,'UniformOutput',false);
for id = 1:length(arg.tv_dims)
    yval = reshape(y(:,id),arg.im_dims);
    dims2 = dims1;
    dims2{arg.tv_dims(id)} = circshift(dims2{arg.tv_dims(id)},[0,1]);
    x = x + (yval(dims2{:}) - yval(dims1{:}));
end
x = col(x);

end
