The linear solver in BBTools was developed to invert linear systems in the presence of noise.
This demo shows how BBTools can be utilized for the reconstruction. This demo may need substantially longer to complete than the other demos, as it use techniques that are more general than the other demos.
To start it directly, use: playshow bbdemo_tomo3
For this demo we shall use a phantom widely used in tomography: the Shepp-Logan phantom. We shall use a variant by Peter Toft, which is substantially easier to handle.
This is available as PHANTOM in the image processing toolbox, but we shall use BBPHANTOM which renders it substantially more accurately. We normalize the intensity to measure 20000 counts (on average).
The phantom have an intense outer ring, a large interior area, ventricles of intensity zero, and several other features that are significantly more difficult to see.
theta=(0:5:175); I_dim=[200,200]; img=bbphantom(I_dim,'toft'); img=img*(20000/(sum(img(:))*length(theta))); m=max(img(:)); clim=[0,1.25]*m; bbimagein(img,clim,'l'); axis image off; colorbar; title('Test phantom');
As with the FBP algorithm, we create the Radon operator and compute the sinogram. Also, we only have access to a noisy measurement of the sinogram, so we create the sinogram by drawing random samples from a Poisson distribution.
[A,S_dim]=bbradon(I_dim,theta); S=A*img(:); 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(sprintf('Sinogram (%i counts)',sum(S(:)))); xlabel('Projection angle'); ylabel('Distance from center');
We start off with a FBP reconstruction for comparison. In most cases, this would be the first step, at very least as a preview.
bb_fbp=bbfbp(S_dim,I_dim,'cut',.23,'hamming'); img_fbp=reshape(bb_fbp*S,I_dim); bbimagein(img_fbp,clim,'l'); axis image off; colorbar; title('FBP reconstruction');
Starting modestly, we will apply the BBTools solver in its pure form. That is, we compute a simple, regularized, solution without doing anything fancy.
The selection of the parameter for the regularization will not be discussed in this demo. Please refer to the image deblurring demos for ways to handle this.
opt=struct('reg',true,'rmaxi',10); pre1=bbpresolve(A,S(:),opt); img1=bbregsolve(pre1,'tikhonov',3.5); img1=reshape(img1,I_dim); bbimagein(img1,clim,'l'); axis image off; colorbar; title('Pure Tikhonov');
Iteration: 3, residual: 4.9437e-01 Iteration: 6, residual: 4.8848e-01 Iteration: 9, residual: 4.8428e-01 Halting: max. number of iterations reached. (Increase rmaxi or rkeep to improve the solution.) Residual norm: 4.8947e-01 (resolved triplets only) Residual norm: 9.1864e-02 (unregularized solution)
The image above is clearly very noisy, but it does capture many important aspects of the image. We may therefore apply a post-filter.
Note that the overall appearance is similar to the FBP reconstrucion. This is very important in tomographic imaging. Doctors spend years training interpreting such images, and we should exploit their skills.
K1=bbkernel([21,21],'gauss',6/21); F1=bbconvn(I_dim,K1,'same'); img1f=reshape(F1*img1(:),I_dim); bbimagein(img1f,clim,'l'); axis image off; colorbar; title('Tikhonov with post-filter');
The parameters for the two reconstruction methods were chosen to provide the same resolution and overall appearance. One may rightly ask why we should use a method that spends a great deal more effort to achieve a similar result. Most of this answer will be addressed in the following demo.
For now, we shall take a look at the so-called "point spread functions" (PSF) for each method. A PSF is the reconstructed image of a single point without noise.
Any image can be viewed as a collection of points, so a linear reconstruction method can be fully characterized by a set of PSFs.
NOTE: this does not hold for a non-linear method such as the EM-algorithm. In fact, the existence of a PSF can be used to define linearity. The lack of a PSF in the EM-algorithm can be an obstacle for interpreting an image.
We start by looking at a collection of PSFs for FBP.
img_psf=bbphantom(I_dim,'points',20,1,'hard'); S_psf=A*img_psf(:); psf_fbp=reshape(bb_fbp*S_psf,I_dim); bbimagein(psf_fbp,'l');axis image off; colorbar; title('PSF for FBP');
Then we do the same for the post-filtered Tikhonov solution.
Even if the match is not perfect, they are at least comparable as claimed. This shows that we are not cheating.
However, when we look at the PSFs we also notice another thing: the PSF of the post-filtered Tikhonov solution is more round than the equivalent FBP image, particularly near the edges of the image.
This may appear to be a subtle difference, but many of the most widely used statistical programs used to analyze the brain relies on this assumption.
pre2=bbpresolve(A,S_psf(:),opt); img2=bbregsolve(pre2,'tikhonov',2.0); img2=F1*img2; img2=reshape(img2,I_dim); bbimagein(img2,'l'); axis image off; colorbar; title('PSF for Tikhonov with post-filter');
Iteration: 3, residual: 6.9610e-01 Iteration: 6, residual: 6.9597e-01 Iteration: 9, residual: 6.7793e-01 Halting: max. number of iterations reached. (Increase rmaxi or rkeep to improve the solution.) Residual norm: 6.7676e-01 (resolved triplets only) Residual norm: 6.8985e-04 (unregularized solution)
We should also compare with the EM-algorithm. Unfortunately, we can no longer compare the PSF to check if we have a fair comparison. As explained previously, the EM-reconstruction is non-linear and reconstructing a point have no direct bearings on a point that forms part of an image.
The most fair comparison we can do is to use post-filtering with the same filter we used for the Tikhonov solution.
img_em=bbempoisson(A,S); img_emf=reshape(F1*img_em,I_dim); bbimagein(img_emf,clim,'l'); axis image off; colorbar; title('EM Reconstruction with post-filter');
At first sight, the EM-image may appear "nicer" than the other images. However, it is questionable whether it actually contains more information. There is also a tendency to create false features that look natural.
Perhaps the one quality that makes the EM-image look better is that it does not contain negative counts. We have consistently shown the negative values also, but is is certainly easy enough to set these to zero.
Below we show the post-filtered Tikhonov clipped at zero, but we would get practically the same result if we used the FBP image. It is now much more questionable whether the EM reconstruction offers a significant advantage.
bbimagein(img1f.*(img1f>0),clim,'l'); axis image off; colorbar; title('Clipped image');
