The EM (Expectation Maximization) is the other major method used to produce images in emission tomography.
To understand this demo, you should first see the demo demonstrating FBP (Filtered Back Projection).
To start it directly, use: playshow bbdemo_tomo2
For comparison, we shall use the same phantom and scanner model as the one we used in the demo of Filtered Back Projection.
I_dim=[150,150]; img=bbphantom(I_dim,'restest',5,[],.1); theta=(0:5:175); [A,S_dim]=bbradon(I_dim,theta); S=A*img(:); clim=[0,2.5]; bbimagein(img,clim,'l'); axis image off; title('Test phantom');
As with the FBP algorithm, we only record a noisy version of the sinogram. While the noise characteristics is mostly ignored by the FBP-algorithm, the EM-algorithm improves on this by modelling the statistical nature of the noise.
S=bbrand('poisson',S); dist=1-2*(0:S_dim(1)-1)/(S_dim(1)-1); S_max=max(S(:)); imagesc(theta,dist,reshape(S,S_dim),[-.5,S_max+.5]); colormap(gray(S_max+1)); colorbar; title('Sinogram with Poisson noise'); xlabel('Projection angle'); ylabel('Distance from center');
The EM (Expectation Maximization) algorithm is not actually "an" algorithm, but a general method which can be used to estimate a maximum likelihood solution for a statistical problem.
In the current setting, we are faced with a problem where the measurements are drawn from a Poisson distribution. BBTools implements this case in the routine BBEMPOISSON, which was derived by Lange and Carson.
Unlike FBP, the EM-algorithm is non-linear and can not be described in terms of a black-box operator. It is an iterative solver and must be called explicitly.
The EM-algorithm only assumes that the system operator and measurements are non-negative, and the measurements are Poisson distributed with a mean given by applying the operator to a solution. This makes it useful for many problems outside tomography.
img_em1=bbempoisson(A,S); bbimagein(reshape(img_em1,I_dim),clim,'l'); axis image off; title('EM-reconstruction');
Just as the previous case, we see the the reconstructed image is very noisy. One popular way to handle this is to post-filter the image, i.e. smoothing the image after it has been reconstructed.
An easy way to do this in BBTools is to use BBCONVN.
K=bbkernel([20,20],'gauss',3/20); F=bbconvn(I_dim,K,'same'); img_em2=F*img_em1; bbimagein(reshape(img_em2,I_dim),clim,'l'); axis image off; title('EM-reconstruction with post-filter');
Another popular way to avoid excessive noise, is early stopping. This instructs the EM-algorithm to stop before it converges to the maximum-likelihood solution.
While this is a popular method, there is not much theory to support it in the case of the EM-algorithm. However, it appears to work well in practice, although it may be difficult to choose the correct number of iterations.
img_em3=bbempoisson(A,S,'maxi',8); bbimagein(reshape(img_em3,I_dim),clim,'l'); axis image off; title('EM-reconstruction with early stopping');
A sieve is a restriction of solutions produced by the EM algorithm. It can be non-linear, e.g. a restriction of the largest concentration in a pixel. In BBTools, a sieve is implemented by processing the image after completing each iteration.
The most common method, however, is a convolutional sieve. This is just a filter applied at each iteration.
K=bbkernel([9,9],'gauss',1/9); F=bbconvn(I_dim,K,'same'); img_em4=bbempoisson(A,S,'sieve',F); bbimagein(reshape(img_em4,I_dim),clim,'l'); axis image off; title('EM-reconstruction with Gaussian sieve');
