Contents

% cone_beam_ct_example.m
% Illustrate cone-beam X-ray CT image reconstruction via FDK and iterative.
% Uses a tiny system size because cone-beam projection can be slow.
% Copyright 2005-1-21, Jeff Fessler, University of Michigan

First run the FDK example. It generates the true image xtrue

and noiseless projection views "proj" and noisy data "yi" and generates (noisy) FDK recon "xfdk" for comparison / initialization.

if ~isvar('xfdk')
	bi = 1e6; % 1M photons / ray
	ri = 0; % no scatter etc. for now
	dfs = inf; % flat!
	feldkamp_example
prompt
end

feldkamp_example 6: cg: cone-beam CT geometry
feldkamp_example 17: fov rmax=255.318
feldkamp_example 21: ig: image geometry
feldkamp_example 31: ell: ellipsoid object
feldkamp_example 41: xtrue: true image volume
feldkamp_example 49: proj: analytical ellipsoid projection views
feldkamp_example 57: li_hat: projection view log data
feldkamp_example 71: fdk
fdk backprojection time: 0.085593
feldkamp_example 143: xtrue: min=0 max=0.02 
FDK error min=-0.00927172 max=0.00884152 
feldkamp_example 145: nrmse = 17.2666%

generate system matrix needed for iterative reconstruction

if ~isvar('A'), printm 'A'
	f.nz_pad = 18; % add this many slices to each side.  todo: how many?
	ig_pad = ig.expand_nz(f.nz_pad);
	A = Gcone(cg, ig_pad, 'type', 'sf2', 'class', 'Fatrix');
%	f.sys_type = aspire_pair(cg, ig_pad, 'system', '3l'); % old way
%	A = Gtomo3(f.sys_type, ig_pad.mask, ig.nx, ig.ny, ig_pad.nz, ...
%		'chat', 0, 'permute213', true, 'checkmask', im&0);
end


% block object for ordered-subsets iterations
if ~isvar('Ab'), printm 'Ab'
	f.nblock = 8;
	Ab = Gblock(A, f.nblock);
end


if 1 % padding functions to help with long object problem
	f.rep = @(x) cat(3, repmat(x(:,:,1), [1 1 f.nz_pad]), x, ...
		repmat(x(:,:,end), [1 1 f.nz_pad]));
	f.pad = @(x) cat(3, zeros(ig.nx, ig.ny, f.nz_pad), x, ...
		zeros(ig.nx, ig.ny, f.nz_pad));
end


% check discrete vs analytical projection (note long object problem)
if 0, printm 'proj check'
	for ii=1:1
		cpu etic
		pp = Ab * f.rep(xtrue);
		cpu etoc 'proj time'
	end
	nrms(pp, proj)
	im clf, im_toggle(proj(:,:,1:12:end), pp(:,:,1:12:end), [0 4.4])
prompt
end

cone_beam_ct_example 19: A
cone_beam_ct_example 30: Ab

regularization object

if ~isvar('R'), printm 'regularizer'
	f.l2b = 2^4.5;
	f.delta = 100/1000;
	R = Reg1(ig_pad.mask, 'type_denom', 'matlab', ...
                'pot_arg', {'hyper3', f.delta}, 'beta', 2^f.l2b);
	W = Gdiag(yi);
	if 1 % check spatial resolution (away from edges)
		psf = qpwls_psf(A, R, 1, ig_pad.mask, W, 'fwhmtype', 'profile');
	prompt
	end
end


if ~isvar('si'), printm 'log sinogram'
	si = log(bi ./ yi); % log sinogram
	xinit = ig_pad.maskit(f.rep(xfdk)); % initalize iterations
	im(ig_pad.embed(xinit))
end

cone_beam_ct_example 58: regularizer
Warn: form_ctc_strum 311: delta=0.1 but cgrad_delta=0.01.  BAD!
qpwls_psf(143): reale: imaginary fraction of  is 0.16445
qpwls_psf 149: gtg, negative FT: -1.14316%, n=78538
qpwls_psf 159: ctc, negative FT: -5.22115e-06%, n=1
qpwls_psf 188: approximate condition number 458.875
qpwls_psf_fwhm 441: fwhm = 1.34187 1.53845 1.094 [pixel]
cone_beam_ct_example 71: log sinogram

OS-SQS iterations for PWLS

if ~isvar('xos'), printm 'start os-sqs pwls iterations'
	f.niter_os = 10;
	xos = pwls_sqs_os(xinit, Ab, reshaper(si, '2d'), R, ...
			'wi', reshaper(W.arg.diag, '2d'), 'niter', f.niter_os);
	xos = ig_pad.embed(xos);
	im(xos), cbar
end

cone_beam_ct_example 79: start os-sqs pwls iterations
Warn: iter_saver 13: using default isave "last"
pwls_sqs_os: 2 of 10 0.1
pwls_sqs_os: 4 of 10 0.2
pwls_sqs_os: 6 of 10 0.2
pwls_sqs_os: 8 of 10 0.2
pwls_sqs_os: 10 of 10 0.3

CG iterations

if ~isvar('xpcg1'), printm 'PWLS reconstruction with CG (no preconditioner)'
	f.niter1 = 20;
	[xpcg1 cost_pcg1] = pwls_pcg1(xinit, A, W, si(:), R, ...
		'userfun', @userfun_cost1, ...
		'niter', f.niter1, 'stop_threshold', 1e-3);
	xpcg1 = ig_pad.embed(xpcg1);
	im(xpcg1), cbar
prompt
end


% circulant preconditioner based on middle voxel and quadratic regularizer
if ~isvar('pre2'), printm 'circulant preconditioner'
	Rq = Reg1(ig_pad.mask, 'beta', 2^f.l2b); % quadratic regularizer
	pre2 = qpwls_precon('circ0', {A, W}, Rq.C, ig_pad.mask); % preconditioner
end

cone_beam_ct_example 89: PWLS reconstruction with CG (no preconditioner)
pwls_pcg1: 3 of 20 0.1
pwls_pcg1: 5 of 20 0.1
pwls_pcg1: 7 of 20 0.1
pwls_pcg1: 9 of 20 0.1
pwls_pcg1: 11 of 20 0.1
pwls_pcg1: 13 of 20 0.2
pwls_pcg1: 15 of 20 0.2
pwls_pcg1: 17 of 20 0.2
pwls_pcg1: 19 of 20 0.2
cone_beam_ct_example 101: circulant preconditioner
qpwls_precon_circ_Mdft(169): reale: imaginary fraction of  is 0.16445
qpwls_precon_circ_Mdft 171: setting -1.14316% to zero

PCG

if 0 | ~isvar('xpcg2'), printm 'PWLS with PCG and circulant preconditioner'
	[xpcg2 cost_pcg2] = pwls_pcg1(xinit, A, W, si(:), R, ...
		'userfun', @userfun_cost1, ...
		'niter', f.niter1, 'stop_threshold', 1e-3, 'precon', pre2);
	xpcg2 = ig_pad.embed(xpcg2);
	im(xpcg2), cbar
prompt
end


if 0 % compare "convergence rate" (of cost) without and with preconditioning
	if im
		clf
		plot(1:f.niter1, cost_pcg1, '-o', 1:f.niter1, cost_pcg2, '-x')
		xlabel 'iteration', ylabel 'cost'
		legend('CG (no precon)', 'PCG2 (circulant precon)')
	end
return
end


% reshape data to be "2d arrays" for OS iterations (subset over last dim)
if ~isvar('os_data'), printm 'os_data'
	if isscalar(bi) && isscalar(ri)
		os_data = reshaper(yi, '2d');
		os_data = {os_data, ...
			bi * ones(size(os_data)), ri * ones(size(os_data))};
	else
		os_data = {reshaper(yi, '2d'), reshaper(bi, '2d'), ...
			reshaper(ri, '2d')}; % all data as 2d arrays
	end
end

cone_beam_ct_example 108: PWLS with PCG and circulant preconditioner
pwls_pcg1: 3 of 20 0.1
pwls_pcg1: 5 of 20 0.1
pwls_pcg1: 7 of 20 0.1
pwls_pcg1: 9 of 20 0.1
pwls_pcg1: 11 of 20 0.2
pwls_pcg1: 13 of 20 0.2
pwls_pcg1: 15 of 20 0.2
pwls_pcg1: 17 of 20 0.2
pwls_pcg1: 19 of 20 0.3
cone_beam_ct_example 130: os_data

OS-SPS iterations for transmission penalized likelihood

if ~isvar('xpl'), printm 'start tpl os-sqs iterations'
	f.niter = 20;
	xs = tpl_os_sps(xinit, Ab, os_data{:}, R, 1+f.niter);
	xs = ig_pad.embed(xs);
	xpl = xs(:,:,:,end);
	im(xpl)
end

cone_beam_ct_example 143: start tpl os-sqs iterations
tpl_os_sps: 3 of 21 0.1
tpl_os_sps: 4 of 21 0.1
tpl_os_sps: 5 of 21 0.1
tpl_os_sps: 6 of 21 0.1
tpl_os_sps: 7 of 21 0.1
tpl_os_sps: 8 of 21 0.1
tpl_os_sps: 9 of 21 0.2
tpl_os_sps: 10 of 21 0.2
tpl_os_sps: 11 of 21 0.2
tpl_os_sps: 12 of 21 0.2
tpl_os_sps: 13 of 21 0.2
tpl_os_sps: 14 of 21 0.3
tpl_os_sps: 15 of 21 0.3
tpl_os_sps: 16 of 21 0.3
tpl_os_sps: 17 of 21 0.3
tpl_os_sps: 18 of 21 0.3
tpl_os_sps: 19 of 21 0.4
tpl_os_sps: 20 of 21 0.4
tpl_os_sps: 21 of 21 0.4

compare FDK vs iterative results

if 1
	iz_good = f.nz_pad + [1:ig.nz]; % original slices
	xpcg1_good = xpcg1(:,:,iz_good);
	xpl_good = xpl(:,:,iz_good);
	xos_good = xos(:,:,iz_good);
	show = @(x) xlabelf('rms = %.1f HU', rms(col(xtrue-x))/0.01*1000);
	clim = [0 0.02];
	im plc 2 3
	im(1, xtrue, 'true', clim), cbar
	im(2, xfdk, 'FDK', clim), cbar
	show(xfdk)
	im(4, xpcg1_good, 'PWLS-PCG', clim), cbar
	show(xpcg1_good)
	im(5, xpl_good, 'PL-OS', clim), cbar
	show(xpl_good)
	im(6, xos_good, 'PWLS-OS', clim), cbar
	show(xos_good)
end

Published with MATLAB® 7.12