This demo explores the same scenario as the previous demo. However, we use less noisy measurements to see what is happening, and use more advanced methods.
To start it directly, use: playshow bbdemo_tomo4
We start with the same phantom as the previous demo. However, this time we allow for 75000 counts, so we can filter much less severely.
theta=(0:5:175); I_dim=[200,200]; img=bbphantom(I_dim,'toft'); img=img*(75000/(sum(img(:))*length(theta))); [A,S_dim]=bbradon(I_dim,theta); S=bbrand('poisson',A*img(:)); 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 with a FBP reconstruction.
This time we can clearly see the ventricles in the center, and we can even start to see that the region in between them is elevated compared to the background.
bb_fbp=bbfbp(S_dim,I_dim,'cut',.45,'hamming'); img_fbp=reshape(bb_fbp*S,I_dim); m=max(img(:)); clim=[0,1.2]*m; bbimagein(img_fbp,clim,'l'); axis image off; colorbar; title('FBP reconstruction');
As before, we may imitate FBP using BBTools. This is important, but hardly exciting given we spend to much time to compute the solution.
To simulate the FBP-image, we regularized by following and followed up using a post-filter. The regularization constant is approximately set using the L-curve heuristic. Please see the deblurring demos for further explanation.
opt=struct('reg',true,'rmaxi',10); pre1=bbpresolve(A,S(:),opt); img1=bbregsolve(pre1,'tikhonov',1.5); K1=bbkernel([21,21],'gauss',2.5/21); F1=bbconvn(I_dim,K1,'same'); img1=reshape(F1*img1(:),I_dim); bbimagein(img1,clim,'l'); axis image off; colorbar; title('Tikhonov with post-filter');
Iteration: 3, residual: 3.2612e-01 Iteration: 6, residual: 3.0982e-01 Iteration: 9, residual: 2.9733e-01 Halting: max. number of iterations reached. (Increase rmaxi or rkeep to improve the solution.) Residual norm: 3.0619e-01 (resolved triplets only) Residual norm: 5.4214e-02 (unregularized solution)
We will now put more of the filtering in the Tikhonov regularization, and just apply a small amount of post-filtering.
Clearly, the result is a drastic improvement over the FBP image. Also note that the solution is very fast to compute, since the expensive presolving step only needs to be done once.
Thus, we may choose the filter parameters after the computation. This could be done by real-time updates as the results are shown on-screen.
img2=bbregsolve(pre1,'tikhonov',15); K2=bbkernel([21,21],'gauss',1.5/21); F2=bbconvn(I_dim,K2,'same'); img2=reshape(F2*img2(:),I_dim); bbimagein(img2,clim,'l'); axis image off; colorbar; title('Tikhonov with post-filter');
The most obvious anatomical information is simply to include information about the location of the head. In BBTools, this can be accomplished using a mask. (This is a part of clinical practice on many sites; it is typically performed to correct for attenuation.) Please see the deblurring demos for details.
However, other anatomical information may be available using information obtained from other sources. This might for example be CT or MRI images scanned separately.
For this example, we will assume we have identified the interior fluid-filled chambers (the so-called ventricles), and most of the large constant region (corresponding to so-called white matter).
We also assume it is known that these regions have constant concentrations. This is typically the kind of information that is known prior to the experiment.
In practice, it is necessary to coregistrate the images, i.e. moving and rotating them to represent the same locations. This can only be done with some amount of certainty, so the regions are drawn with a border around them.
Tools for region drawing are too specialized to be included in BBTools. We cheat and call a separate program to generate them from the original phantom.
mask=bbphantom(I_dim,'ellip',I_dim.*[.48,.36],'crisp'); bbdemo_tomoreg; tmp=reg+mask-1; bbimagein(tmp+img.*(tmp==0),[-1,0,3,3],'l'); axis image; colorbar title('Mask and anatomical regions');
We now put it all into the pot and stir. It is tempting just to reconstruct the image, but this would just look like the regions were averaged on one of the previously reconstructed images.
Tikhonov regularization works by favoring images with constant concentration, provided the data allows it. Therefore, we do not gain much if restricting nearly constant regions to be exactly constant.
What we can do, however, is to use a very small amount of regularization. The other areas of the image may be very bad, but if the constant regions are sufficiently large, then they will be accurate. That is, we use the reconstruction as an advanced form of partial volume correction.
inx=find(mask); Mi=bbindex(inx,numel(img)); Mr=bbmean(reg,'same'); A3=A*Mr*Mi'; pre3=bbpresolve(A3,S(:),opt); img3=Mi'*bbregsolve(pre3,'tikhonov',2.0); img3=reshape(img3,I_dim); bbimagein(img3,clim,'l'); axis image off; colorbar title('Solution with anatomical constraints');
Iteration: 5, residual: 2.9370e-01 Iteration: 10, residual: 2.8357e-01 Halting: max. number of iterations reached. (Increase rmaxi or rkeep to improve the solution.) Residual norm: 2.8357e-01 (resolved triplets only) Residual norm: 9.5406e-02 (unregularized solution)
Now that we have determined some parts of the image with high accuracy, we may return to the remaining image and reconstruct that separately. That is, we subtract the known part from the right hand side, and disclude it from further consideration.
We now start to see several details in the orginal phantom. For instance, we start to discern the two largest regions in the bottom of the phantom. These would never be distinguishable without a solid amount of anatomical information.
img4r=img3.*(reg>0); S4=S-A*img4r(:); Z4=bbdiag(reg(:)==0); A4=A*Z4*Mi'; pre4=bbpresolve(A4,S4(:),opt); img4=Mi'*bbregsolve(pre4,'tikhonov',15); img4=reshape(F2*(img4+img4r(:)),I_dim); bbimagein(img4,clim,'l'); axis image off; colorbar title('Reconstruction using anatomy and regularization');
Iteration: 5, residual: 4.3906e-01 Iteration: 9, residual: 4.3820e-01 Halting: max. number of iterations reached. (Increase rmaxi or rkeep to improve the solution.) Residual norm: 4.2113e-01 (resolved triplets only) Residual norm: 1.3234e-01 (unregularized solution)
