20/07/2011

[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

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!

30 comments:

  1. Thank you very much for providing the source code. It is amazing.

    ReplyDelete
    Replies
    1. Happy to help, have a nice day

      Delete
    2. Thanks 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?

      Delete
    3. Hi,

      It 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

      Delete
    4. 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.
      in 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 !!


      Delete
    5. As i understood the code:
      instead 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.

      Delete
  2. 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?

    The 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

    ReplyDelete
  3. Yes in TEM you have the choice to scan the object on both axis or only single axis.
    In 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.




    ReplyDelete
  4. hi Stefano,
    i tried to run your source code with an .Bmp image with 2027x2027 dimension, but it's not working with me :(

    ReplyDelete
    Replies
    1. 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).

      I 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.

      Delete
  5. i tried to run your code with a Bitmap image 2026x2026 but with no result ;(
    i 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



















    ReplyDelete
    Replies
    1. Mmm.. do

      figure, 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.

      Delete
  6. This comment has been removed by the author.

    ReplyDelete
    Replies
    1. I forgot something so i removed it and posted it again under this post.

      Delete
  7. Why 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!

    ReplyDelete
    Replies
    1. Hello Chen,

      I 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

      Delete
    2. Hi Stefano, the Lena 256*256 image I have is a grayscale image...

      Delete
    3. Well it seems worked somehow. I adjusted the normalization and deletion of <0 and >1 pixels part

      Delete
    4. Glad 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?

      Noised=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

      Delete
  8. Hi 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!!

    ReplyDelete
    Replies
    1. I 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.

      Delete
    2. Dear Chen, I'm glad you managed to solve your issue and I thank you for taking the time of reporting it back to me.

      I 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

      Delete
  9. How is SIRT different from ART?

    ReplyDelete
    Replies
    1. Hey Jorge,

      not 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

      Delete
    2. 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.
      it 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.

      Delete
    3. Yes, the most information you can find online is in papers, I just picked a random one, maybe there are better ones out there.

      However, 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!

      Delete
    4. thanks again, I ended up changing the code to:
      for 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!

      Delete
    5. I understand, glad it helped :)

      Delete
  10. HELLO I WANT TO MAKE A 2D AND 3D RECONSTRUCTION ON MATLAB BY X TOMOGRAPHY USING ALGEBRIC METHODS
    HOW CAN HELP ME

    ReplyDelete
    Replies
    1. Hey Fatma,

      I 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

      Delete

With great power comes great responsibility.

Da grandi poteri derivano grandi responsabilità.