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
Get the source code here. 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!
Thank you very much for providing the source code. It is amazing.
ReplyDeleteHappy to help, have a nice day
DeleteThanks Stefano. I wonder if I can actually use it for the 3d volume reconstruction of the 2D projections obtained in Transmission Electron Microscopy on a tilt axis +- 60?
DeleteHi,
DeleteIt may be adaptable since the basic principle is more or less the same: get a lot of pictures of an object from various angles then try to reconstruct what it looked like. It doesn't matter if you get the images by bombarding a huge object on a table with X-rays with a CT scanner or from a tiny particle under the microscope.
But the resulting image is always the 2D "top" view since the images themselves are all taken on the same plane from a side perspective. What you usually do then, is slice the object into tiny segments, reconstruct each one singularly, and finally overlay them all to create the 3D representation.
Pretty much like those 3D sculpture puzzles for example (http://www.hasbro.com/common/instruct/3DSculptures-StarWarsJarJarBinks.PDF)
As for the tilt axis, the main issue is usually blurring if the range is not wide enough but these reconstruction techniques (see also http://groglogs.blogspot.it/2011/07/matlab-em-expectation-maximization.html) were created exactly for this scenario, were you have a limited range of projections to work with and they work amazingly well I'd say.
If I remember correctly, we tested these algorithms with ranges as low as 45 degrees (90 total) with surprisingly good results.
Additionally, there are a lot of public papers describing reconstruction techniques to apply when the viewing range is limited, and many of those are adaptation of the "basic" SIRT, you may find them useful.
Have a nice day
Thanks for the response Stefano. However I must admit that i am new to these techniques and still feel myself confused in some of the basic things. I wonder if you can write me an algorithm to implement the steps for the 3d volume reconstruction.
Deletein my understanding this code will help me to reconstruct again a 2D image as you mentioned also. now i dont know how to slice these projections and then reconstruct again.
Have a nice day !!
As i understood the code:
Deleteinstead of the phantom image, i will use my own central image (at 0 degree) and then instead of sinograms in the image i will use my tilt series but again if i use the tilt series i will eventually get a 2d reconstructed image.
I have never done TEM reconstructions and I wouldn't know where to begin. Isn't there a way you can scan the object to get the projection, then move up/down a fixed axis then scan again to get another projection and so on?
ReplyDeleteThe 2D image you get is only one of the many 2D images (slices) you need to recreate the 3D object. The picture 8-4 on page 359 from "Scanning Transmission Electron Microscopy: Imaging and Analysis" by Stephen J. Pennycook and Peter D. Nellist (http://books.google.it/books?id=N_98Hef0ZTYC&lpg=PA359&ots=JNQ9otS1Tm&pg=PA359#v=onepage&q&f=false) might give you a better view of the whole process.
What remains to check is how adaptable the SIRT technique I described is towards your goal. SIRT is a great technique and I'm confident enough that there's someone who explored an adaptation for your particular case, but I do not have any reference to give you unfortunately. I would suggest searching Google Scholar (http://scholar.google.it/) which is a great source for academic material
Yes in TEM you have the choice to scan the object on both axis or only single axis.
ReplyDeleteIn my case the images that i have are taken only on single axis. and the range of the axis is -61 to +61 degrees. hence total 121 images.
hi Stefano,
ReplyDeletei tried to run your source code with an .Bmp image with 2027x2027 dimension, but it's not working with me :(
Hi, if you look closely at the code, you'll see we tested it on a phantom image created by MatLab (http://it.mathworks.com/help/images/ref/phantom.html - this is what we'd like to reconstruct) which we then "scanned" using the Radon function (http://it.mathworks.com/help/images/ref/radon.html - this is the sinogram, something not readily comprehensible tot he human eye).
DeleteI don't know how you modified your code and what does that bitmap represent but I'd suggest you to check if the image is being read correctly and that you didn't forget to remove or edit some lines from my code. Eg: if your image represents something, switch F = phantom with F = your image; if your image is a sinogram, switch R = radon instead. You can also use imshow (http://it.mathworks.com/help/images/ref/imshow.html) to visually check what's happening in your program.
i tried to run your code with a Bitmap image 2026x2026 but with no result ;(
ReplyDeletei tried to change the value from double to uint8 but i get nothing ;(
%Author: Stefano "grog" Ghio, Michele Bozzano, Bruno Mazzarello
%Released under CreativeCommons and GPLv2 licenses
%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
%clear matlab environment
clear all;
close all;
theta = 0:179; %projection angle 180 degrees
F = imread('akrem.bmp'); %create new test phantom 128x128 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
%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
%At applied to g aka noisy alien
At = iradon(Noised,theta,'linear', 'none', 1,2026); %reconstruct noisy alien
figure, imshow(At,[]),title('At G'); %show noisy alien
S2 = (sum(sum(At))); %calculate pixels sum on At
At = (((At)/S2)*S1); %normalize At so that pixel counts match
n = 10;%iterations
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,2026);% reconstruct alien using Fk unfiltered sinogram
%normalize
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
%delete values <0 aka not real!
index = Fk<0;
Fk(index)=0;
%delete values >1 aka not real!
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
Arrerr(k) = sum(sum((F - Fk).^2));
%stop when there is no more improvement
if((k>2) &&(Arrerr(k)>Arrerr(k-1)))
break;
end
end
figure, plot(Arrerr),title('Arrerr');%show error improvement over all iterations
Mmm.. do
Deletefigure, imshow(F,[]),title('Alien vs Predator'); %show alien
and
figure, imshow(R,[]),title('Sinoalien');%show sinogram
show something? The first one should produce your image, the second one should produce the sinogram from which you'll reconstruct the image.
Also, I don't know if this is going to work on non grayscale images.
This comment has been removed by the author.
ReplyDeleteI forgot something so i removed it and posted it again under this post.
DeleteWhy I tried your code on the Lena 256*256 image but the reconstruction results are not acceptable? And especially after Figure 4, it just gives white and black areas that divided the whole image. And the error seems does not converge but oscillate. Could you possibly tell me anything why it does not work well on this Lena 256*256 image? I have changed the parameters because of the change of the input image size. Anyway really appreciate your work, it works so well on the phantom image!
ReplyDeleteHello Chen,
DeleteI don't know if this is going to work on non-grayscale images. Could you try by editing Lena's image before processing it?
Regards
Hi Stefano, the Lena 256*256 image I have is a grayscale image...
DeleteWell it seems worked somehow. I adjusted the normalization and deletion of <0 and >1 pixels part
DeleteGlad to hear that. I'm sorry I didn't look at the code before, but did you still have the lines to add the noise?
DeleteNoised=imnoise(R,'poisson');%add Poissonian noise
The normalization and subsequent value deletion part was tied to this. Though the algorithm should work just as well!
If I find the time I'll try with my own greyscale Lenna to check what was was the issue.
In the meantime thank you for your answers,
Have a nice day
@chen mu may i know what you did or adjust? Cause i got the same problem
DeleteYou can adjust the normalization to different scales, like 0 to 255 and some normalizations maybe not needed. Try change the normalization parts for colored images.
DeleteThank you very much. I am doing my thesis and your reply was such a big big help. I didnt expect that you will respond. Thank you so much
DeleteBut i always get error in Arrerr. Do you know how to solve that?
DeleteHi I have another problem about you code. You used 'linear' and 'none' command in your iradon function, which corresponds to unfiltered backprojection (See http://www.mathworks.com/help/images/ref/iradon.html). And I tried to remove the 'linear' and 'none' command, so iradon will result in filtered backprojection. But the result of the error will not even converge. Why the filtered backprojection work even worse than the unfiltered backprojection? Thanks!!
ReplyDeleteI verified that it has sth to do with the poisson noise you add to R. Actually why you add the poisson noise? And I tried if you do not do that and remove the 'linear' and 'none' command in your iradon function, which means you use the filtered backprojection, you get even smaller error, as small as less than 5, compared with the unfiltered backprojection, which has error around 100.
DeleteDear Chen, I'm glad you managed to solve your issue and I thank you for taking the time of reporting it back to me.
DeleteI will have to review the code again since apparently I need to better comment all the sections, sorry for the confusion.
Obviously it's nonsense to add noise or skip on the enhancing features the iradon function offers in a real life scenario, why would one intentionally harden his job?
In this example we did that to show just how good this algorithm is. The idea was to show the algorithm at work, not to provide a direct generic implementation, which can nevertheless be derived from the posted code, as you did. Plus, not using the filtered back projection was by design, since if you have it, you would be able to apply other algorithms other than SIRT which particularly shines under this adverse conditions; also, it's not always possible to obtain a decent filtered back projection, eg if your angle is narrow or the overall scan quality is poor.
If you look at all the imshow we make and firstly check the noised image and then the reconstruction you get to fully appreciate it. Had I simply shown it without the added noise, it would have been harder to spot the differences between this approach and other reconstruction algorithms, especially if you consider that this delivers high quality results even if you apply it under such harsh circumstances.
Thank you again for your feedback
Have a nice day
How is SIRT different from ART?
ReplyDeleteHey Jorge,
Deletenot much. The basic idea is the same and SIRT is just an improved ART in the sense that ART will calculate and immediately apply the correction factor to single cells during iteration, while SIRT will average those corrections before applying them.
SIRT is then slower but the image quality is better because it attenuates the impact of neighbour corrections.
If you'd like to see the actual math behind them, there are some whitepapers you can Google eg: http://www.iaea.org/inis/collection/NCLCollectionStore/_Public/43/048/43048827.pdf
Unfortunately sites such as Wikipedia are a bit poor on this topic.
Cheers
Thank you very much! :) I really like how the paper describes the per cell "patch"/correction before applying to all arrays as you explained. I understand the main concept a little bit better. It sounds like SIRT is acting like a soft thresholding and ART is not. Now, I am trying to translate that into your code. Basically, I am trying to figure out how to modify your code to create a ART reconstruction.
Deleteit seems that the important steps are:
St = sum(sum(t)); % here is where you add all elements of the re
t = (t/St)*S1; % where is where the averaging happens
Can I just get rid of the St? Why defined S1 as sum(sum(F))? in the ideal world you won't know your image before hand. Any help will be highly appreciated. Thanks for your previous and rapid answer.
Yes, the most information you can find online is in papers, I just picked a random one, maybe there are better ones out there.
DeleteHowever, why do you want ART? I don't remember the whole math but I think even back when I was studying, it was already stated that one does not use ART in practice since you have better alternatives in terms of iterative algorithms or just go with filtered backprojection whenever you have a good enough angle.
Now, unfortunately I did this more than 5 years ago so it takes a while to remember stuff and I dug my university notes and Matlab references just for you :)
As a result I would say (more of a fuzzy best guess honestly) you want to change this part:
St = sum(sum(t)); --gives you the average of correction using all rays at once
t = (t/St)*S1; --applies average correction to image and normalize
you can check the references here: https://ch.mathworks.com/help/matlab/ref/sum.html, https://ch.mathworks.com/help/images/ref/radon.html and https://ch.mathworks.com/help/images/ref/iradon.html.
So the two lines above mean that we get ALL rays at once to determine how to correct our guess, but since in ART we just focus on a single ray at a time instead, I think in Matlab this means something along the lines of (not actual code, can't remember how to write it):
for i=t_start:t_end --for all matrix rows
St[i] = sum(t)[i]; --get the correction for a single ray
t = t/St[i]; --apply correction from single ray
end
t = t*S1; --normalize
Also the normalization is not mandatory or even to be done exactly like we did. This is also why the algorithm works if you do not have the starting image already, there would be no point in even trying to reconstruct then anyway. Note that we also use it here:
Arrerr(k) = sum(sum((F - Fk).^2)); --F is the actual image we are reconstructing
to determine how good are we correcting, but in practice you would use the current (Fk) and previous (Fk - 1) reconstruction steps to understand how good you're doing to determine when to stop. Basically the reconstruction will improve (reducing the approximation error) before degrading again (noise amplification error), we want to stop as soon as we determine the noise amplification error is becoming greater than the approximation reduction.
Now that's the best I can do I'm sorry.
Also is almost midnight on a Friday, let's say my brain is just not in the game anymore.
Cheers!
thanks again, I ended up changing the code to:
Deletefor column = 1:size(t,2)
St = sum(sum(t(:,column)));
t(:,column) = (t(:,column)/St);
end
t = lamda*t/(max(t(:)));
where lamda is just a weighting factor, .22 for my project, but you can tune it.
I am interesting in ART since I am the TA for a medical imaging class and I am building a CT reconstruction Demo for the students. I just wanted to introduce them to the most popular techniques and ART is "historically" important. My research field is MRI, therefore this demo was a little out of my league but I really appreciate your help. As expected, I will reference and cite your functions. Thanks!
I understand, glad it helped :)
Deletehi can i have these codes ?
DeleteHello, sure the code is available at https://digilander.libero.it/otacoconvention/source/bi/SIRT.m
DeleteCheers
HELLO I WANT TO MAKE A 2D AND 3D RECONSTRUCTION ON MATLAB BY X TOMOGRAPHY USING ALGEBRIC METHODS
ReplyDeleteHOW CAN HELP ME
Hey Fatma,
DeleteI think you could start by running one of these examples. The one in this post shows how to reconstruct a test image after adding some noise to it.
Cheers
Why do i get error in using "-" in F - Fk ?
ReplyDeleteHello,
Deleteit is a bit difficult with little information like this, I would suggest you do some debug first.
Maybe F or Fk is null? It might depend on how they are initialized. Is the image loaded there existing?
Also when you copy the code, check if some strange formatting is added, it might happen sometimes when copying from browsers; maybe that minus sign looks like "-" but is some strange dash instead.
Cheers
Hi. Its my first time analyzing codes in matlab and i am not really familliar with sirt and other algorithms. I dont what i am going to adjust or edit when i am using my own grayscale picture. Can you help me please?
ReplyDeleteCan I still able request about this ?
ReplyDeleteHello, yes we can try
Delete