Simple code for iterative image reconstruction (see

[@arturjosefreitas8260]

Loved the video, it had a lot of helpful information to start my project. Thanks a lot, keep the great content!!

[@AndrewJReader]

Many thanks for the feedback, really appreciated. Glad it was helpful for starting your project!

[@be419]

Thank you so much for the informative video! I'm new to using MatLab-- where did you import the data of the sinograph for your code to work on it?

[@AndrewJReader]

Thanks for the feedback! In this code I used the phantom that is available in Matlab, and then used the radon function to create the sinogram data from that phantom. Hence I did not need to import any data for the sinogram.

[@charalmposkor7949]

In ratio command inside the loop it says to me that it's an error because 'arrays have incompatible sizes for this operation'. How do i solve that?

[@AndrewJReader]

Hi, have you checked the sizes of the arrays? If so, what size are they? Are you using "./" in Matlab to ensure element-wise division? Hope this helps you to find the solution?

[@charalmposkor7949]

@@AndrewJReader I have to say that this video was very useful and no other helped me so much. Yes i used "./" but how i check the arrays? Furthermore the function sens is not recognised.

[@AndrewJReader]

@@charalmposkor7949 I was suggesting that you could print to screen the size of the arrays to try and understand what is happening (e.g. using the function "size" in Matlab). In the video, I think I used "sens" as an image array, based on backprojecting a sinogram filled with ones. So it is not a Matlab function. Hope that helps?

[@charalmposkor7949]

​​@@AndrewJReaderwhy sens is not recognised?Also it says that ineed positive integers. How do i change that?

[@AndrewJReader]

@@charalmposkor7949 sens is defined later in the video, nearer the top of the code. Sorry that the text is not crystal clear

[@marianiemusarudin1575]

Hi Prof. Do you have any video explaining on OSEM reconstruction using matlab?

[@AndrewJReader]

Hi, thanks for the comment. I have not done a video on this yet, but it is not a difficult extension. It would look something like the following (the below is not optimised, nor tested, but hopefully clear enough to give the idea, feedback welcome): % Use ordered subsets (OS), just update based on a subset of the measured data number_subsets = 1; for ii = 1:num_iterations for ss=1:number_subsets subset_of_angles_to_use = azi_angles(ss:number_subsets:end); sensitivity_image = iradon(sinogram_ones(:,subset_of_angles_to_use),subset_of_angles_to_use,'none',xd); % Forward model the current reconstruction estimate fp_rec = radon(mlem_recon, subset_of_angles_to_use); % Ratio ratio = noisy_sino_gt_psf(:,subset_of_angles_to_use) ./ (fp_rec + 0.001); % BP the ratio bp_ratio = iradon(ratio,subset_of_angles_to_use,'none',xd); % Divide by SIGMA_i a_ij (the sensitivity image) correction_image = bp_ratio ./ (sensitivity_image + 0.001); mlem_recon = mlem_recon .* correction_image; end % Subset loop end % Iteration loop

[@arturjosefreitas8260]

Hello. I'd like to know which articles were used to perform this algorithm, thank you

[@AndrewJReader]

Hi, this algorithm is based on the 1982 paper by Shepp and Vardi, published in IEEE Transactions on Medical Imaging (for maximum likelihood image reconstruction). Most people refer to the method as ML-EM, or EM-ML. A faster version is called ordered subsets EM (OSEM).

[@arturjosefreitas8260]

@@AndrewJReader Once again, thank you :))

[@anasssofti9271]

Thank you for your wonderful simulation and implementation, I have a simple question about why the BP operator (A.T) is actually the Transpose of the forward projection (A), I couldn't get the relation if you can please guide me in this matter.

[@AndrewJReader]

Many thanks for the feedback! When deriving MLEM, we find the need for A^T. This means we need to swap the rows and columns of A. The matrix A contains rows which are lines, each line being used in scalar product with the input vector. When considering A^T, it therefore has columns which are lines. To apply A^T to an input vector means that the input vector contains weighting factors for each column (i.e. a weigthed set of lines, which is backprojection). For more info please see this video: kzbin.info/www/bejne/d4OqeX-NoMtor7c

[@anasssofti9271]

@@AndrewJReader Thanks for the clear response + video recommendation, It's really clear in my mind, your videos are like a knowledge reference to me.😄

[@biamorais2914]

Hey!! Thank you for the video. I am having some troubles. I have the code done till line 8 but when I run it it only shows the figure 1 when nothing in it. I cant see the sheep Logan phantom in it. Idk what's wrong

[@AndrewJReader]

Can you copy and paste the precise lines of code that you have typed in? Perhaps there might be an error that has crept in?

[@danish07delhi]

Thanks a lot, its very helpful!

[@AndrewJReader]

Thanks for the feedback, glad it helped!

[@TheFreeSpiritKID]

How to do that to cbct data?

[@AndrewJReader]

You would need a function for the forward and backprojection model for the cone-beam CT geometry. In this video, just parallel projections are used.

[@witnesschirume333]

Thanks for a very informative video, but i have a question on what you are referring to as measured data. How do you get this data is it experimental data? and if so in what format will it be and how do you import this data to matlab?

[@AndrewJReader]

Thanks for the feedback. In this video the "measured data" is simulated, to represent what we would obtain from a scan. Raw scan data is often available in sinogram format, or else in list-mode format (which can be binned into sinograms). Perhaps this video can help: kzbin.info/www/bejne/aXOmYZWkjNSVgtU or indeed this one for actual reconstruction of measured data: kzbin.info/www/bejne/fGi1YoN5bLugqpI Hope that helps?

[@kaoutharchh3167]

Can you help me I look for Reconstruction 3d images brain tumor in matlab

[@AndrewJReader]

For 3D there are many possible scanner geometries and of course different modalities (PET, SPECT, CT,...). Hence more specific software is often needed. Which modality and which scanner are you interested in?

[@runekotzerhzn5511]

Hey Andrew, thanks for the video! I am using this algorithm for my bachelors thesis. It works well, but I ran into a problem. I want to add gaussian noise to the measured sinogram, but whenever this is done my reconstruction is barely recognizable. Maybe you can give me a tipp on where to change my code for this to work , thanks : )

[@AndrewJReader]

Thanks for your comments. The simplest fixes would be to use the algorithm for a smaller number of iterations, or to smooth the image after finishing a chosen number of iterations. Otherwise you need a MAPEM method rather than MLEM. I do have lectures on MAPEM on this channel. Or you can put a CNN into the method, .... there are many possible noise compensation strategies (regularisation strategies).

[@runekotzerhzn5511]

@@AndrewJReader Thanks, I'll check out your video on the MAPEM method. So if I understood correctly there is unfortunately no way of adapting my MLEM algorithm to work with a noisy sinogram.

[@AndrewJReader]

@@runekotzerhzn5511 MLEM in this video is intended for sinograms with Poisson noise. Even then, images which are reconstructed from noisy data, based only on data fidelity, tend to be noisy. You can change the spatial basis functions, to avoid insisting on pixel perfect reconstructions. So you could use Gaussian spatial basis functions, for less noisy (and less sharp) images.

[@runekotzerhzn5511]

@@AndrewJReader I fixed the issue with MLEM not working at all (reconstruction is barely recognizable), I referred to in my original comment. I think the problem was that my addition of gaussian noise caused negative sinogram values, which should not be the case! Now I only add positive noise values to the sinogram and the reconstruction improves for a few iterations, then quickly gets worse - as you pointed out in your first answer (small number of iterations). I also tried out the SART algo. from Skimage, it also can't handle noise in the way I am currently doing it. Shame because I would like to show that iterative algorithms can perform well with noise. I was not able to find your video on the MAPEM method, can give me the video title? Thanks

[@AndrewJReader]

@@runekotzerhzn5511 Right, indeed can only have positions values with MLEM. Please see about 38 minutes into this video kzbin.info/www/bejne/i33biYqOn9GYopY

[@amianifineug1353]

Please code back projection I want to use it thank you

[@AndrewJReader]

Backprojection can be done just with iradon and no filter in this video. To code it from first principles, you can see what I did in Python in this video: kzbin.info/www/bejne/eIm7faKMrqaHZrM (that creates the matrix for forward projection, and so one can simply do the transpose of that to get backprojection). Better though to adapt the forward projection code directly to backprojection code.

[@amelieteache3852]

Hello, first of all, thank you for awesome video. It helped a lot. Could I just ask one question ? I have sinogram of phantom from SPECT. If i want to use it as an input, should I make inverse RT, then RT to obtain sinogram or should I use directly my measured data? I reconstructed FBP with Ram-Lak filter and with no filter. If I use the reconstructed image with Ram-Lak as measured data and not filtred image as first approximation it works perfectly. Once I put my measured sinogram without any modification, I recieve just black figure. Thank you in advance!

[@AndrewJReader]

Hi, thanks for the question. For this kind of iterative reconstruction, MLEM, you should use the raw measured data. Of course, for SPECT, there will however be many other effects that you should ideally model (attenuation, for example) or else precorrect for (precorrection is not ideal though).(For an example in PET, see: kzbin.info/www/bejne/fGi1YoN5bLugqpI)A raw measured SPECT sinogram, in general, will not be consistent with the Radon transform, unless all the other effects are first compensated for by modelling or else corrected. How did you do your FBP with Ram-Lak filter reconstruction? Did you use your raw measured sinogram data, and then use iradon in Matlab? Perhaps you used other software, that does the precorrections?

[@amelieteache3852]

@@AndrewJReader Yes, exactly, I used my raw data and used iradon in Matlab. Also I was thinking the problem could be cause by the size of matrix? When I use my raw measured data, I have 128x120, but when I make radon transform of FBP as first approximation I recieve 131x120 so I forced the transform to use output_size=max(size(raw_data). If I make radon transform of my FBP with Ram-Lak, it creates matrix of the same size as non filtred image I used as first approximation. Could this cause some kind of problem?

[@AndrewJReader]

@@amelieteache3852 Have you tried removing the last argument (xd) that I put into my iradon? That forces the output of iradon to be that size in x and y. Not sure if you used an argument like that in your FBP? It is certainly better to use MLEM on the raw measured data, and not on a forward projection of a FBP reconstructed image. Indeed, just need to get the sizes of the sinogram and images to be correct. Also, please be sure to have no negative values in your measured sinogram data from the start, as MLEM is a non-negative algorithm.

[@amelieteache3852]

@@AndrewJReader I have following code: raw_data=abs(S); %measured sinogram output_size=max(size(raw_data)); FBP=abs(Ram_Lak); %reconstructed image with Ram-Lak filter angles = 1:1:179; figure; for i = 1:20 FP=radon(FBP,angles,output_size); comp= raw_data./(FP+0.0001); BP=iradon(comp, angles, 'none'); FBP=FBP.*BP; subplot(4,5,i); imagesc(FBP);colormap(gca,gray); end I also compared measured sinogram and front projection of reconstruction. The measured sinogram really does not correspond very well with the FP and when I divide it, it already shows whole black image during the first iteration.

[@AndrewJReader]

@@amelieteache3852 If I were you I would display every single image and sinogram at every single stage of your code, step by step. As you know, I use imshow with min and max specified, but if you know and trust imagesc and colormap, then ok. So, display the raw_data after abs(S). Also, report its dimensions, check its min and max. Likewise for FBP. Then, in your loop, display FP, and also report its dimensions (e.g. print them to screen). Presumably the dimensions match for raw_data and FP? Next, display comp. Probably there is an issue with the division? So you can display the min and max of comp, using min(comp(:)) and max(comp(:)). Do the values look ok, or are there Inf or Nan results? Then display BP and inspect its size, and the range of values. I think you can see my approach to solving bugs - display everything (images, dimensions, min and max), and check everything - it will soon be apparent where the results do not match what we expect. I could probably debug if I had the raw data, but that is probably going way beyond a comment / chat on KZbin!