% ir_et_ml_em_demo.m
% illustrate emission tomography with E-ML-EM algorithm

if ~isvar('sino')
	ig = image_geom('nx', 60, 'ny', 64, 'dx', 1);
	[xtrue xell] = ellipse_im(ig, 'shepplogan-emis', 'oversample', 3);

	im plc 2 4, clim = [0 8];
	im(1, xtrue, clim, 'xtrue'), cbar

	sg = sino_geom('par', 'nb', 60, 'na', 40, 'dr', ig.dx, ...
		'strip_width', 'd', 'orbit', 180);
	sino = ellipse_sino(sg, xell, 'oversample', 2);

	im(5, sg.s, sg.ad, sino)
	xtick([-1 0 1] * sg.rfov), ytick([0 90])
end

if ~isvar('yi'), printm 'yi'
	rng(0)
	scale = 1e5 / sum(sino(:)); % 1e5 counts
	yi = poisson(scale * sino) / scale;
	im(6, sg.s, sg.ad, yi)
	xtick([-1 0 1] * sg.rfov), ytick([0 90])
end

if ~isvar('xfbp'), printm 'fbp'
	tmp = fbp2(sg, ig);
	xfbp = fbp2(yi, tmp);
	im(2, xfbp, 'fbp', clim), cbar
prompt
end

if ~isvar('A'), printm 'A'
	A = Gtomo2_strip(sg, ig); % tomography system matrix

	if 0
		tmp = A * xtrue;
		im(tmp), cbar
		max_percent_diff(tmp, sino)
	end
end

const = sum(col(yi - yi .* log(yi + eps)), 'double'); % trick for plot
kl = @(y, v) (sum(col(v - yi .* log(v)), 'double') - const) / 1e5 * 2;

reg = 0; % replace with "10" (e.g.) for exponential regularized case

if 1 % ML-EM
	niter = 30;
	nrmse = nan(niter,1);
	nlogl = nan(niter,1); % negative log likelihood
	x = ig.ones;
	a = A' * sg.ones;
	for iter=1:niter
		x = x ./ (a + reg) .* (A' * div0(yi, A*x)); % E-ML-EM
		nrmse(iter) = 100 * nrms(x(:), xtrue(:));
		Ax = A * x;
		nlogl(iter) = kl(yi, Ax);
		if 1
			im(3, x, clim), cbar
			im(7, sg.s, sg.ad, Ax)
			xtick([-1 0 1] * sg.rfov), ytick([0 90])
			im subplot 4
			plot(1:niter, nrmse, '-o'), axis([0 niter 0 70])
			titlef 'NRMSE %%'
			im subplot 8
			plot(1:niter, nlogl, '-o'), axis([0 niter 0 0.6])
			xlabelf 'iter', titlef 'Neg. Log. Like'
			drawnow
		end
	end
end

%	x = eql_sps_os(x, A, yi, 1, 0, R, ...) % regularized version
