% psf_mismatch_example1.m%% 1D example showing the effect of PSF mismatch on ML-EM algorithm%% Copyright 2001-8-24, Jeff Fessler, The University of Michigan%% generate data%if ~isvar('yi'), printm 'data'	n.x = 128;	n.y = 1;	xtrue = zeros(n.x, n.y) + 1;	xtrue(n.x/2+1+[-5:5]) = 100;	ghist = 'normcdf(x+1/2, 0, fwhm/sqrt(log(256))) - normcdf(x-1/2, 0, fwhm/sqrt(log(256)))';	ghist = inline(ghist, 'x', 'fwhm');	nk = 21;	xk = [-(nk-1)/2 : (nk-1)/2]';	f.fwhm0 = 5;	f.fwhm1 = 7;	f.fwhm2 = 2;	kern0 = ghist(xk, f.fwhm0);	printf('kern0 discard %g', 1-sum(kern0))	kern0 = kern0 / sum(kern0);	G0 = filtmat('1d', kern0, n.x);	kern1 = ghist(xk, f.fwhm1);	printf('kern1 discard %g', 1-sum(kern1))	kern1 = kern1 / sum(kern1);	G1 = filtmat('1d', kern1, n.x);	kern2 = ghist(xk, f.fwhm2);	printf('kern2 discard %g', 1-sum(kern1))	kern2 = kern2 / sum(kern2);	G2 = filtmat('1d', kern2, n.x);	ri = 5;	yi = G0 * xtrue + ri;	if im		plot(1:n.x, yi, '-o'), xlabel i, ylabel 'y_i', title Data	endpromptend% uniform initial imagexinit = ones(size(xtrue));%% ML-EM iterations%if ~isvar('x2'), printm 'ML-EM for various psf models'	f.niter = 81;	x0 = eml_em(xinit, G0, yi, 1, ri, [], f.niter);	x1 = eml_em(xinit, G1, yi, 1, ri, [], f.niter);	x2 = eml_em(xinit, G2, yi, 1, ri, [], f.niter);	im clf, im pl 2 2	im(1, x0, 'x0')	im(2, x1, 'x1')	im(3, x2, 'x2')promptend	fw_true = fwhm1(xtrue, 'imid', n.x/2+1);	fw0 = fwhm1(x0(:,2:end), 'imid', n.x/2+1);	fw1 = fwhm1(x1(:,2:end), 'imid', n.x/2+1);	fw2 = fwhm1(x2(:,2:end), 'imid', n.x/2+1);if 1 & im	ii = 1:(f.niter-1);	clf, plot(ii(1:4:end), fw0(1:4:end), 'yo', ...		ii(1:4:end), fw1(1:4:end), 'cx', ...		ii(1:4:end), fw2(1:4:end), 'g+', ...		[min(ii) max(ii)], fw_true * [1 1], 'r--', ...		ii, fw0, 'y-', ...		ii, fw1, 'c-', ...		ii, fw2, 'g-')	legend(	sprintf('True system PSF, FWHM=%g', f.fwhm0), ...		sprintf('Model PSF too big, FWHM=%g', f.fwhm1), ...		sprintf('Model PSF too small, FWHM=%g', f.fwhm2), ...		'ideal FWHM')	xlabel 'Iteration of ML-EM algorithm'	ylabel 'FWHM of reconstructed signal'end