[MatLab] EM – Expectation Maximization reconstruction technique implementation

EM – Expectation Maximization is an iterative algorithm used in tomographic images (as in CT) reconstruction, very useful when the FBP – Filtered Back Projection technique is not applicable.

Basic formula is:

f_k+1 = (f_k / alpha) (At (g / (A f_k)))

where:

• f_k is solution (the resulting reconstructed image) at k-th iteration, at first iteration it is our guess
• g is image sinogram (what we get from the scanning)
• A f_k is Radon transform of f_k
• At is inverse Radon
• alpha is inverse Radon of a sinogram with all values 1 which represents our scanning machine

%EM - Expectation Maximization
%Formula is:
%f_k+1 = (f_k / alpha) (At (g / (A f_k)))
%f_k is solution at k-th iteration, at first iteration it is our guess
%g is image sinogram
%A f_k is Radon transform of f_k
%At is inverse Radon of its argument
%alpha is inverse Radon of a sinogram with all values 1 which represents our scanning machine

%clear matlab environment
clear all;
close all;
theta=[0:89]; %limited projection angle 90 degrees instead of 180
F=phantom(64); %create test phantom 64x64 pixels aka alien
figure, imshow(F,[]),title('Alien vs Predator'); %show alien
S1 = sum(sum(F));%calculate pixels sum on F
R = radon(F,theta); %apply Radon aka scan the alien to create a sinogram
figure, imshow(R,[]),title('Sinoalien');%show sinogram
%set all values <0 to 0
index=R<0;
R(index)=0;
n = 2000;%iterations
Fk=ones(64);%our initial guess
%create alpha aka pixels sum of the projection matrix (our scanning machine)
sinogramma1=ones(95,90);
alpha=iradon(sinogramma1, theta,'linear', 'none', 1,64);
%calculate relative error
relerrfact=sum(sum(F.^2));
for  k=1:n
     Afk = radon(Fk,theta);%create sinogram from current step solution
     %calculate g / A f_k=Noised./(A f_k+eps); aka initial sinogram / current step sinogram
     %eps is matlab thing to prevent 0 division
     GsuAFK=R./(Afk+eps);
     retro=iradon(GsuAFK, theta, 'linear', 'none', 1,64);%At (g / (A f_k))
     %multiply current step with previous step result and divide for alpha updating f_k
     ratio=Fk./alpha;
     Fk=ratio.*retro;
     %normalize
     St = sum(sum(Fk));
     Fk = (Fk/St)*S1;
	 %calculate step improvement
     Arrerr(k) = sum(sum((F - Fk).^2))/relerrfact;
	 %stop when there is no more improvement
     if((k>2) &&(Arrerr(k)>Arrerr(k-1)))
        break;
     end     
end
figure, imshow(Fk,[]),title('Fk');%show reconstructed alien
figure, plot(Arrerr),title('Arrerr');%show error improvement over all iterations

%compare our result with the one we would have had using the FBP - Filtered Back Projection
easy=iradon(R,theta, 'Shepp-Logan',1,64);
figure, imshow(easy,[]),title('FBP');

%calculate error between EM and FBP - with limited image size and projection degree FBP is bad!
FBPerr=sum(sum((F - easy).^2))/relerrfact;

[MatLab] SIRT - Simultaneous Iterative Reconstruction Technique implementation

SIRT – Simultaneous (algebraic) Reconstruction Technique is an iterative algorithm used in tomographic images (as in CT) reconstruction, very useful when the FBP – Filtered Back Projection technique is not applicable.

Basic formula is:

f_(k+1) = f_k + At (g - A f_K)

where:

• f_k is solution (the resulting reconstructed image) at k-th iteration, at first iteration it is our guess
• g is image sinogram (what we get from the scanning)
• A f_k is Radon transform of f_k
• At is inverse Radon

Note: this sample code is to showcase the algorithm. It does not use FBP for the iradon, it adds noise, it normalizes the image values as per OUR needs. Given these difficult conditions, it still produces an amazing result.

Please check the comments carefully and tailor it to your use case!

    

%SIRT - Simultaneous Iterative Reconstruction Technique
%Formula is:
%f_(k+1) = f_k + At (g - A f_K)
%f_k is solution at k-th iteration, at first iteration it is our guess
%g is image sinogram
%A f_k is Radon transform of f_k
%At is inverse Radon of its argument

%This example shows how to scan an image with a limited angle, obtain the sinogram, add noise to make things worse, and reconstruct the initial image with great accuracy
%You may want to remove/edit some sections if you intend to apply this to your use case:
%- skip/correct normalization
%- skip noise
%- use FBP for the iradon

%clear matlab environment
clear all;
close all;
theta = [0:179]; %projection angle 180 degrees
F = phantom(128); %create new test phantom 128x128 pixels aka alien. You can use your own image here if you want. It was not tested on non-grayscale images
figure, imshow(F,[]),title('Alien vs Predator'); %show alien
S1 = sum(sum(F));%calculate pixels sum on F
R = radon(F,theta);%apply Radon aka scan the alien to create a sinogram
figure, imshow(R,[]),title('Sinoalien');%show sinogram
%add image noise to the sinogram. You can skip this part, it's just to showcase the algorithm, it's not obviously needed
maximum=max(max(R));
minimum=min(min(R));
R=(R-minimum)/(maximum-minimum);
%normalize between 0 and 1. You can tune this as per your needs
fact=1001;%set the number of X-rays, the higher the better (and deadlier)
R=(fact/10^12)*R;
Noised=imnoise(R,'poisson');%add Poissonian noise
figure, imshow(Noised,[]),title('Dirty alien');%show noisy sinogram
%This part below you actually need it! Edit accordingly to your case
%At applied to g aka noisy alien
%Check the parameters here if you edited the code above! You might also want to tune the parameters for the iradon. Here we're showing that this works even when the FBP is not available!
At = iradon(Noised,theta,'linear', 'none', 1,128); %reconstruct noisy alien.
figure, imshow(At,[]),title('At G'); %show noisy alien
%algorithm starts crunching here
S2 = sum(sum(At)); %calculate pixels sum on At
At = (At/S2)*S1; %normalize At so that pixel counts match. Edit as per your needs
n = 100;%iterations. Might be more, there's always a limit over which it doesn't make sense to keep iterating though!
Fk = At;%Matrix Fk is our solution at the k-th step, now it is our initial guess
for  k=1:n
    t = iradon(radon(Fk,theta),theta, 'linear', 'none', 1,128);% reconstruct alien using Fk unfiltered sinogram. Maybe use FBP if you want here
    %normalize. Again, as per your needs
    St = sum(sum(t));
    t = (t/St)*S1;
    %update using (At g - At A f_k) 
    %new Fk = Fk + difference between reconstructed aliens At_starting - t_previuous_step
    Fk = Fk + At - t;
	%remember that our alien is normalized between 0 and 1. Might not be your case!
    %delete values <0 aka not real! Might not be your case!
    index = Fk<0;
    Fk(index)=0;
    %delete values >1 aka not real! Might not be your case!
    index = find(Fk>1);
    Fk(index)=1;
    %show reconstruction progress every 10 steps
    if(mod(k,10)== 0)
    figure,imshow(Fk,[]),title('Fk');
    end
    %calculate step improvement between steps. Tune as per your needs
    Arrerr(k) = sum(sum((F - Fk).^2));
    %stop when there is no more improvement. Tune as per your needs
    if((k>2) &&(Arrerr(k)>Arrerr(k-1)))
       break;
    end
end
figure, plot(Arrerr),title('Arrerr');%show error improvement over all iterations

[MatLab] Filter a tomographic image in Fourier’s space

MatLab has built-in functions to simulate acquisition and elaboration of tomographic (as in CT) images. When it comes to filtering the image prior to back-projecting it, it is possible to do it yourself without relying on the (good as in Shepp-Logan) filters MatLab has.

 

We will filter our image in Fourier’s space. For each column of the sinogram, we:

• apply the Fourier transform
• filter it
• apply the Fourier antitransform

When we reconstruct our image WITHOUT having MatLab apply any filter, we’ll see the image, filtered with our filter, as a result.

%filter an image in Fourier's space

%clear matlab environment
clear all;
close all;
theta = [0:179]; %projection angle 180 degrees
F = phantom(256); %create new test phantom 256x256 pixels aka alien
figure, imshow(F,[]),title('Alien vs Predator'); %show alien
R = radon(F,theta); %apply Radon aka scan the alien to create a sinogram
figure, imshow(R,[]),title('Sinoalien'); %show sinogram
%get Shepp-Logan filter, you can use any filter you want
Fi = phantom(128);
Ri = radon(Fi,theta);
[I, Filter]=iradon(Ri,theta, 'Shepp-Logan');
%add image noise to the sinogram
maximum=max(max(R));
minimum=min(min(R));
R=(R-minimum)/(maximum-minimum);
%normalize between 0 and 1
fact=1001;%set the number of X-rays, the higher the better (and deadlier)
R=(fact/10^12)*R;
Noised=imnoise(R,'poisson');%add Poissonian noise
figure, imshow(Noised,[]),title('Dirty alien');%show noisy sinogram
%add zero-padding
Padded(1:512, 1:180)=0;
Padded(1:367, :) = Noised;
figure, imshow(Padded,[]),title('Padded alien');%show padded sinogram
%filter the noise with our filter in Fourier's space for each column of the sinogram
for k=1:180
	%FPadded = filtered and padded
    FPadded(:, k)=fft(Padded(:, k));%Fourier transform
    FPadded(:, k)=FPadded(:, k).*Filter;%Filter in Fourier
    Padded(:, k)=real(ifft(FPadded(:, k)));%Fourier antitransform of the filtered and still padded sinogram column
end
%remove padding
R=Padded(1:367,:);
figure, imshow(R,[]),title('Unpadded filtered alien');%show final sinogram
I=iradon(R, theta, 'None'); %back project without filtering to show final result
figure, imshow(I,[]),title('Final alien');