| |
Regularization in medical imaging
Example
Consider the backprojection code for forward
and inverse CT scanning shown
here.
We develop code to implement an inverse solution
for backprojection for the following regularization
schemes:
- Tikhonov Regularization
- One advanced regularization technique, such as:
MAP, ML, etc.
- Code to imlement one technique to estimate
the hyperparameter value: eg. GCV
- One edge preserving regularized inverse
Forward problem
Unfiltered backprojection may be formulated as using
the transpose of the sensitivity matrix as the inverse.
This Matlab/octave code will implement that function,
using the function makeproj.m.
(Increase spc to use less memory or run faster)
spc=.05; rlim=1;
[x,y]= meshgrid(-rlim:spc:rlim,-rlim:spc:rlim);
plen= size(x,1);
img = (x.^2 + y.^2) > rlim;
img( x >.45 & x<.65 & y>-.05 & y<.45) =1;
img( x >-.55 & x<-.25 & y>.45 & y<.65) =1;
angls= [0:15:179.9];
% create sensitivity matrix prm
prm= [];
for ang= angls
prma= makeproj(ang*(pi/180),x,y);
prm= [prm;prma];
end
prm= sparse(prm);
% calculate measurements
proj = prm* img(:);
Original Image
imlim= ( x.^2 + y.^2 > rlim );
imwrite('img1.png',img*255);
|
Backprojection (unregularized)
imbp= reshape( prm' * proj, plen, plen);
imbp = imbp -min(imbp(:)); imbp= imbp/max(imbp(:));
imbp= imbp.*(~imlim) + imlim;
imwrite('img2.png',imbp*255);
| 
|
|
Tikhonov Regularization (λ=1.0)
H=prm;
I=eye(size(H'*H));
imt=(H'*H+I)\H'*proj;
imt= reshape( imt , plen, plen);
imt= imt-min(imt(:)); imt= imt/max(imt(:));
imt = imt.*(~imlim) + imlim;
imwrite('img3.png',imt*255);
| 
|
Tikhonov Regularization (λ=0.01)
imt=(H'*H+.01*I)\H'*proj;
imt= reshape( imt , plen, plen);
imt= imt-min(imt(:)); imt= imt/max(imt(:));
imt = imt.*(~imlim) + imlim;
imwrite('img4.png',imt*255);
| 
|
Tikhonov Regularization (λ=10×−6)
imt=(H'*H+1e-6*I)\H'*proj;
imt= reshape( imt , plen, plen);
imt= imt-min(imt(:)); imt= imt/max(imt(:));
imt = imt.*(~imlim) + imlim;
imwrite('img5.png',imt*255);
|
| |
Last Updated:
$Date: 2004-06-29 07:57:32 -0400 (Tue, 29 Jun 2004) $
|
Engineering