hole_example1.m

实现PET/SPECT 幻影图像regression的matlab源代码 algorithms for Poisson emission tomography PET/SPECT/ Poisson

% hole_example1.m
%
% 2D example showing the "hole in sphere" effect for OSEM
%
% Copyright 2004-5-13, Jeff Fessler, The University of Michigan

if ~has_aspire, return, end

%
% generate data
%
if ~isvar('G0'), printm 'data'
	f.dir = test_dir;
	f.dsc0 = [f.dir 't0.dsc'];
	f.wtf0 = strrep(f.dsc0, 'dsc', 'wtf');
	if 1
		Fwhm0 = 3;
		os_run(['wt -chat 0 dsc 12 fwhm_detector 5 >! ' f.dsc0])
		os_run(['echo y | wt -chat 0 gen ' f.dsc0 ' row'])
%		os_run(['echo y | wt -chat 0 col2row ' f.wtf0 ' ' f.wtf0])
	end

	G0 = Gtomo2_wtmex(f.wtf0);
	ig = G0.arg.ig;
	sg = G0.arg.sg;
	nx = ig.nx; nb = sg.nb;
	ny = ig.ny; na = sg.na;
end

if ~isvar('Gb0'), printm 'Gb'
	f.nblock = 6;
	Gb0 = Gblock(G0, f.nblock, 0);
	Gb20 = Gblock(G0, 20, 0);
prompt
end


if ~isvar('yi'), printm 'data'
%	xtrue = [6.5 -0.5 5 5 0 100]; over = {'oversample', 4};
%	xtrue = ellipse_im(ig, xtrue, over{:});
	over = {'oversample', 1};
	xtrue1 = ellipse_im(ig, [18.5 -0.5 5 5 0 100], over{:});
	xtrue2 = ellipse_im(ig, [-3.5 -0.5 3 3 0 100], over{:});
	xtrue3 = ellipse_im(ig, [-18.5 -0.5 2 2 0 100], over{:});
	xtrue = xtrue1 + xtrue2 + xtrue3;
	im(xtrue, 'xtrue'), cbar

	ri = 1;
	yi = reshape(G0 * xtrue(ig.mask), nb, na) + ri;
%	yi = G0 * xtrue + ri;
	im(yi, 'yi'), cbar
prompt
end


% uniform initial image
xinit = ig.ones;

% FBP
if 0 & ~isvar('xfbp')
	tmp = fbp2(sg, ig);
	xfbp = fbp2(yi-ri, tmp, 'window', 'hann');
	im(xfbp, 'fbp'), cbar
prompt
end

if ~isvar('xe0'), printm 'EM'
	f.niter_em = 5;
	xe0 = eml_em(xinit(ig.mask), G0, yi(:), 1, ri, [], f.niter_em);
	xe0 = ig.embed(xe0);
	im clf, im(xe0, 'em'), cbar
prompt
end

xinit = xe0(:,:,end); % 5th EM iter

%
% OS-EM iterations
%
if ~isvar('xo0'), printm 'run os-em'
	f.niter = 10; % for testing
	if isempty(caller_name)
		f.niter = 200; % for full comparison
	end

	rri = ri * ones(nb,na);
	xo0 = eml_osem(xinit(ig.mask), Gb0, yi, [], rri, f.niter);
	xo0 = ig.embed(xo0);
	im clf, im(xo0(:,:,end), 'osem'), cbar
prompt
end

% regularization (edge preserving)
if ~isvar('R'), printm 'R'
	f.l2b = -9;
	f.l2b = -11; % for binary case
	f.potential = 'hyper3'; f.delta = 3;
%	f.potential = 'cauchy'; f.delta = 10; % this one used for grant
	f.potential = 'cauchy'; f.delta = 2;
	R = Robject(ig.mask, 'type_denom', 'matlab', ...
		'beta', 2^f.l2b, 'potential', f.potential, 'delta', f.delta);
	clear xp0
prompt
end

% Incremental EM-3, aka C-OS-3
if ~isvar('xp0')
%	xp0 = epl_os_emdp(xinit(ig.mask), Gb0, yi, ones(size(yi)), rri, R, ...
%		0+1*f.niter, 200, 1);
	xp0 = epl_inc(xinit(ig.mask), Gb20, yi, ones(size(yi)), rri, R, ...
		'niter', f.niter, 'hds', 3, 'os', 40, 'pixmax', 200, 'chat', 0);
	xp0 = ig.embed(xp0);
	im(xp0(:,:,end))
prompt
end

% regularization with perfect side info
if ~isvar('Rs'), printm 'Rs'
	ug = ugibb_form(xtrue > 50, 'threshold', 0.5, 'offsets', R.offsets);
	im clf, im(ug)
	f.l2b = -9;
	Rs = Robject(ig.mask, 'type_denom', 'matlab', ...
		'beta', 2^f.l2b, 'user_wt', ug);
	clear xs0
prompt
end

if ~isvar('xs0'), printm 'recon with side info'
	xs0 = epl_inc(xinit(ig.mask), Gb20, yi, ones(size(yi)), rri, Rs, ...
		'niter', f.niter, 'hds', 3, 'os', 40, 'pixmax', 200, 'chat', 0);
	xs0 = ig.embed(xs0);
	im(xs0(:,:,end))
prompt
end


if 1 & im % examine "hole in sphere"
	iy = ny/2;
	t = squeeze(xo0(:,iy,:));
	t = squeeze(xp0(:,iy,:));
	im(t)
	subplot(211)
	ix = 1:nx;
	ix = 10:60;
	plot(ix, xtrue(ix,iy), 'y-o', 1:nx, t);

	plot(	ix, xtrue(ix,iy), 'y-x', ...
		ix, xo0(ix,iy,end), 'c-o', ...
		ix, xp0(ix,iy,end), 'g.-', ...
		ix, xs0(ix,iy,end), 'r.-')
		legend('True', 'ML: OS-EM', 'PL: C-OS-3')
	xlabel 'Horizontal Pixel Index'
	ylabel 'Activity [arbitrary units]'
	title 'Profile'
	axisy([0 140])
end

if 0 & im % total activity per iteration
	t = xtrue3 > 0;
	sum_true = sum(xtrue(t));
	t1 = reshape(xo0, [], f.niter+1); t1 = sum(t1(t,:)) / sum_true;
	t2 = reshape(xp0, [], f.niter+1); t2 = sum(t2(t,:)) / sum_true;
	t3 = reshape(xs0, [], f.niter+1); t3 = sum(t3(t,:)) / sum_true;
	ii = 0:f.niter;
	subplot(212)
	plot( ...
		ii(1), t3(1), 'y-s', ...
		ii(1), t2(1), 'g.-', ...
		ii(1), t1(1), 'c--o', ...
		ii(1:8:end), t3(1:8:end), 'ys', ...
		ii(1:8:end), t2(1:8:end), 'g.', ...
		ii(1:8:end), t1(1:8:end), 'co', ...
		ii, t3, 'y-', ...
		ii, t2, 'g-', ...
		ii, t1, 'c--', ...
		ii, 1, 'y:')
	legend('PL side', 'PL: C-OS-3', 'ML: OS-EM', 4)
	legend('PL: C-OS-3 with side information', 'ML: OS-EM', 4)
	xlabel 'Iteration', ylabel 'Activity Recovery'
	axis([0 120 0.8 1.02])
	axis([0 120 0.6 1.02])
%	savefig hole_example_recovery1 c
%	set(gcf, 'InvertHardCopy', 'off'); print -djpeg90 hole_example_recovery2
return
end

if 1 & im % examine "hole in sphere"
	t = xtrue1(:) > 99;
	t1 = reshape(xo0, [nx*ny f.niter+1]);
	t1 = t1(t, :);
	t2 = reshape(xp0, [nx*ny f.niter+1]);
	t2 = t2(t, :);
	ii = 0:f.niter;
	subplot(212)
	plot(	...
		ii(1), max(t1(:,1), [], 1), 'c-o', ...
		ii(1), min(t1(:,1), [], 1), 'c--x', ...
		ii(1), max(t2(:,1), [], 1), 'g-', ...
		ii(1), min(t2(:,1), [], 1), 'g--', ...
		ii(1:8:end), max(t1(:,1:8:end), [], 1), 'co', ...
		ii(1:8:end), min(t1(:,1:8:end), [], 1), 'cx', ...
...%		ii(1:8:end), max(t2(:,1:8:end), [], 1), 'g-', ...
...%		ii(1:8:end), min(t2(:,1:8:end), [], 1), 'g--', ...
		ii, max(t1, [], 1), 'c-', ...
		ii, min(t1, [], 1), 'c--', ...
		ii, max(t2, [], 1), 'g-', ...
		ii, min(t2, [], 1), 'g--', ...
		ii, 100, 'y:')
	legend('ML: OSEM max', 'ML: OSEM min', 'PL: C-OS-3 max', 'PL: C-OS-3 min')
	xlabel 'Iteration'
	ylabel 'Activity [arbitrary units]'
	title 'Maximum and Minimum Values within ROI'
	axisy([0 165])
return
end

if 1 % figure showing reconstructions for Ken
	ix = 1:nx; iy=1:ny;
%	ix = 7+[nx/4:(3*nx/4)]; iy=ny/4:(3*ny/4);
	t = [xtrue(ix,iy); xo0(ix,iy,end); xp0(ix,iy,end)];
	t(t>120) = 120; t = uint8(255 * t/120);
	im clf, im(t), cbar
%	imwrite(t', 'hole.tif')
end

% savefig fig_hole_example1_new