version 1.14 (172 KB) by
Kyungsang Kim

3D Cone beam CT (CBCT) projection backprojection FDK, iterative reconstruction Matlab examples

This program is about 3D cone-beam CT for flat detector.

This is made for students who start to learn the CT medical imaging.

GPU option:

- I provide the projection, back-projection codes using built-in GPU functions (semi-GPU codes)

- Above version Matlab 2013b with parallel computing toolbox are needed.

- you can switch on/off the usage of GPU.

- The computational time of this code is several times faster than CPU-based code. Speed depends on your GPU.

Examples:

- FDK ( known as filtered backprojection)

- MLEM (Maximum Likelihood Expectation Maximization),

- SART (Simultaneous Algebraic Reconstruction Technique)

- SQS (Separable Quadratic Surrogates).

You can run:

0. Please check "ParamSetting.m" carefully.

1. Make a measurement: MeasurementGen.m: F5

1.1 Poisson Noise can be imposed.

2. Reconstruction examples

2.1 FDK

2.2 MLEM

2.3 SART

2.4 SQS

Please see this: https://onedrive.live.com/redir?resid=4599730E4EF6090E!4249&authkey=!AJPJM8_qtvrhuGc&ithint=file%2cpdf

Thank you

Kyungsang Kim (2021). 3D Cone beam CT (CBCT) projection backprojection FDK, iterative reconstruction Matlab examples (https://www.mathworks.com/matlabcentral/fileexchange/35548-3d-cone-beam-ct-cbct-projection-backprojection-fdk-iterative-reconstruction-matlab-examples), MATLAB Central File Exchange. Retrieved .

Created with
R2013b

Compatible with any release

**Inspired:**
CT Basic reconstruction algorithms, TIGRE: Tomographic Iterative GPU-based Reconstruction Toolbox, LAVI DBT-Reconstruction

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Meysam Tavakolizhiyuan quxuanhaowawawa,thank you very much! It help me a lot!

Shu MengRuk LeeThanks a lot for the code.

Can you please explain the algorithm in projection.m file?

Misha HilarioHi Kyunsang where is the img28.mat file?

Heng XuAnyway know how to deal with Gantry tilted case?

Tengiz LobzhanidzeLuo ronghuiromdhane hamidahow can i visualize the volume 3D after the reconstruction by MLEM?

Thinh NguyenScott MorochIf I have a set of 2D projections, how do I convert these images into the "img128.m" data file?

John HoffmanSuper helpful code. Was able to run it and get a reconstruction on the first try with my own data. Saved me from having to implement FDK from scratch! As with all reconstruction code: read the documentation *first* before attempting to do a reconstruction so that you understand how the code is configured. And *understand your geometry* and CT geometry in general.

Nadim Abrarhow can i run this code ?? My CUDA driver or my GPU device is not proper.

Amy BeckerWhere is the input into the Matlab code for the location of the central ray from the x-ray source onto the detector? This would be a location on the detector in pixels [vertical and horizontal location]

liang xiaokunfkreuzHi, how can I make the 'img' from my own 2D images. Thank you

lei denghitesh tekchandaniWhich Body Part it is?

121Christopher RockTo use this package,

1. Extract the ZIP file to a folder and add that folder to your path.

2. Run MeasurementGen(). This creates the proj.mat file which is missing when you try to run the other files.

Hany KashaniHi all I dont seem to find the main Demo.m code , could someone please advise?

thanks

Hany

Carla RomanoI put in input a matrix 138*157*101. The only way to run MeasurementGen.m is to impose in param setting x=157 and y=138. Is it normal? After that the code back projection gave me this error

Error in backprojection (line 25)

vol(:,:,iz) = (Ratio.*interp2(proj',pu,pv,param.interptype));

Could someone help me?

Ka MirulHi Aubrey, you should convert your projection image to double or single data format. When you have image input yo can use imread and then stack it into 3D data.

Aubrey Kelley-CogdellHello All,

What format should i insert 'img' as?

I have the series of 2D image projections as TIF

Matheus SilvaHello, Kyungsang Kim!

Could you refer me to the paper where you studied to implement the MLEM?

Thanks in advance.

hyt 1220Who knows which paper this file refered to？

yiyuchentangthanks you for the codes!

it`s helpful for me!

dectnice program. just wonder the following line in filtering.m

proj(:,:,i) = fproj(end/2-param.nu/2+1:end/2+param.nu/2,:)/2/param.du*(2*pi/param.nProj)/2*(param.DSD/param.DSO);

can you explain why we need the last several terms instead of just fproj(end/2-param.nu/2+1:end/2+param.nu/2,:)?

thanks,

sigridkatamaranThank you for the code.

Is there beam hardening correction algorithm?

bo zhouThe projection x-ray is mono-energitic or poly-energitic spectrum and what with kVp and mAs?

Thank you so much

Huang YanYangood job!besides I wonder if you can tell me which paper did your work refer to ?

Mu Sungood code.

Ander BiguriNote that the projection implementation has an asymptotically increasing error when it reaches to vertical edges at 45 degrees. This is due to the fact that the way to project and back project is to "stretch" the image to transform the cone to a parallel, thus making the closest single voxels (at 45) incredibly wide. You can test this by projecting a single image with a full white "cube image" at 45 degrees. Strange artifacts arise. This is less relevant the closer to the centre of the image one is. The toolbox is still great to learn about medical imaging, but this error makes it (unfortunately) not a good toolbox for research purposes or real image reconstruction.

Nick TsuiYou need to run MeasurementGen.m first to generate the proj.mat. Could have been a lot less confusing if proj.mat is pre-generated in the folder. Thanks though.

Nick TsuiWhere is your proj.mat file? cannot find it?

Farid ZDear Kyungsang, may I ask two questions?

1- In comments you mentioned page 96 Jiang Hsieh 2nd edition; however, I cannot follow the code in projection.m, i.e., how it is representing the formula 3.59. May I ask for your hints on it?

2- How can we reach from a CT images of a patient to the sample file which you have provided?

Thank you.

GERARDO RAMIREZyw628Sorry my mistake, I found proj.mat

yw628Hello,

Where is proj.mat file?

Diego GarciaGreat contribution, thank you very much for the code.

Ali MeghoufelAli MeghoufelGreat contribution for others!

Thank you Kyungsang Kim!

AntonReally cool! Thank you for this code! Especially I like to see a real phantom and animated reconstructions.

SetaThank you Kyung :-)

Kyungsang KimHi Seta,

Please do as following:

1) Open "ParamSetting.m"

2) Change to "param.gpu = 0"

Then, please try again.

This option needs "Parallel Computing toolbox" and "GPU (Hardware)",

If you have a GPU, please check memory of GPU and number of cores.

Also, it would be useful above 2013b version.

SetaThank you for the code. Is there a description file to explain what each file does and how to run the code please?

When I run Recon1_FDK or MeasurmentGen, they both come up with errors saying some file or function is missing or undefined. More specifically, when running Recon1_FDK error suggests proj.mat does not exist and when I run MeasurementGen the error suggest gpuArray (called in Projection.m) is undefined.

Amny thanks for assistance in advance.

matHamdi MahdjoubBellaWould it be possible to use this program to forward project a 3D image where dx=dy, but dx=dy!=dz? Right now, my parameters are:

%% Parameter setting %%

param.nx = 400;

param.ny = 400;

param.nz = 66;

%The real detector panel pixel density (number of pixels)

param.nu = 400;

param.nv = 400;

% Detector setting, according to Varian Trilogy OBI (real size)

param.su = 350; % mm

param.sv = 350; % mm

% X-ray source and detector setting

param.DSD = 1500; % Distance source to detector

param.DSO = 1000; % X-ray source to object axis distance

% This confusing stuff

res=(param.DSO/param.DSD)*(param.su/param.nu);

% (1000/1500)*(350/400)

param.sx = param.nx*res; % mm

param.sy = param.ny*res; % mm

param.sz = param.nz*res; % mm

but this doesn't create the correct forward projections, so far as I have been able to check. Am I trying to do something that is impossible in image processing?

Thank you!

Rama Teja A SPlease include analytical equation for better understanding. Or please explain the algorithm in projection.m. If possible please include analytical equation references

Kyungsang KimHi, Erick,

That's a good point.

This factor is used in the weighted backprojection for fan-beam and cone-beam geometries.

(see Jiang Hsieh 2nd edition, Figs 3.35 (reason) and 3.37(geometry))

Result is similar compared to the previous one, In real geometry, this factor is not critical factor for the image, because CT rotates 360 deg (average factor ~ 1) as well as itself is close to 1.

However, I think I omitted it for fast calculation, soon I will fix it for accuracy.

Thanks a lot.

ps) you can directly send e-mail to me, it is sometimes difficult to check comments here.

erickHi, I compared you work with the FDK equation in some papers, it seems that your FDK algorithm omited a factor. I added this factor in the "backprojection.m" of your algorithm, the code is shown as follows.

Ratio=param.DSO.^2./(param.DSO-ry).^2;

vol(:,:,iz) = Ratio.*interp2(uu,vv ,proj',pu,pv,'linear');

Do you think it is necessary to add this factor? I am waiting for your reply, Thank you.

Kyungsang KimHi Erick,

First Question:

If you have real projections ("RealProj"), then first calculate:

1. I0 = max(RealProj(:)); % if you don't know I0

2. Proj = -log(RealProj/I0);

Here, "I0" is the blank image (intensity of source) as mentioned in the text book.

If you think that the final image has a problem, check 1) "RealProj" has 0 values or 2) "I0" is too high than background values.

Solution of cases

1) step 3. Proj(isinf(Proj)) = 0;

2) step 1. Define I0 as a background (high intensity) average.

step 2. Proj = -log(min(RealProj/I0,1));

Then, you can execute FDK or MLEM as in Demo.m.

Second Question:

Parameter setting is very important. Image resolution is defined by detector size and distances of DSO, DSO. In your case, DSD is 700 and DSO is 560. so ratio is 560/700 = 0.8. Then detector size of a pixel is 50/1024, so maximum resolution of a voxel is 50/1024*0.8.

Now, voxel resolution "dx=dy=dz" is defined. then define the number of voxels that covers whole images (nx, ny, nz).

Thus, sx = nx*dx, sy=ny*dy, sz=nz*dz;

Last Question:

I modified just for fast calculation. But if you can see Jiang Hsieh's book (second edition,page 96), it is easy to understand.

Thanks,

erickThank you for your answers. I have a few questions left.

If I have already got real projections from x-ray imaging system, for example, I got the projections of a mouse, do I need to do the "CTprojection"? Under this circumstance, if I only need to do "CTbackprojection", how to set following params about the reconstructed object:nx, ny, nz, sx, sy, and sz.

The params of the detector are setted as follows:

param.su=50;%mm

param.sv=50;%mm

param.nu=1024;

param.nv=1024;

param.DSD=700;

param.DSO=560;

Last question, would you please give some reference about FDK algorithm not using iso-center(DSO) domain.

Thank you.

KyungsangHi Erick,

1. w = cos(theta), to calculate w, my code considers distances in detectors, as your comment, many books calculate in iso-center (DSO) domain. In this case, uu and vv should be differently calculated by considering distance ratio from source to center and source to detector. Please think about angle between iso-center line and extension line source-voxel-detector.

2. In my code, pi/2 rotate shift is for calculating exactly your comment. "interp2" function results in 90 degree rotated image. So I did "angle" -> "angle-pi/2", and I think DSO+ry is right.

3. If you want to get as real projection data, first set blank scan intensity "b", this is related to "dose", then measurement M = b*exp(-projection(obj)).

Thanks your comment.

erickI am wondering if I have real projected images, do I need to use the function "CTprojection"?

erickGood work. I have some questions.

Firstly, in the file of "backprojection.m", you set "w = (param.DSD)./sqrt((param.DSD)^2+uu.^2 + vv.^2);", but in many books, they said that w should be setted as "w = (param.DSO)./sqrt((param.DSO)^2+uu.^2 + vv.^2);". Can you give me some reason?

Secondly, also in the "backprojection.m", I think the description of rx, ry, pu and pv should be:

rx=xx.*cos(angle)+yy.*sin(angle_rad);

ry=-xx.*sin(angle)+yy.*cos(angle_rad);

pu=rx.*(param.DSO)./(param.DSO-ry);

pv=param.zs(iz)*(param.DSO)./(param.DSO-ry);

Would you please give me some suggestion, thank you.

Kyungsang KimHi, Masoud Hashemi,

This program is for flat detectors,

I think you should modify for the arc detector geometry,

Thanks for using my software.

Masoud HashemiI am wondering if the projection function is for Arc detectors (3rd Generation) or Flat detectors?

Tai ChiehGreat!

ljpzgxThe codes can run but the result is not very good,can you upload the referrences you refer for the codes,so tha we can understand clearly. Thank you！

Kyungsang KimThis program is similar to other programs but

much faster.

To accelerate speed,

(1) define the reconstruction object and projection images as a "single" precision,

and then (2) use "matlabpool".

Kyungsang KimYou can use "matlabpool",

From

for i=1:n (angle)

projection or backprojection

end

To

parfor i=1:n (angle)

projection or backprojection

!!please delete figure plot functions

end

This increases the performance more than twice in quad cores.

Kyungsang KimPlease change the function code "projection_pixel.m":

From

dist = (SAD+SDD)./sqrt((SAD+SDD)^2 + uu.^2 + vv.^2)*abs(ys(2)-ys(1));

To

dist = sqrt((SAD+SDD)^2 + uu.^2 + vv.^2)./(SAD+SDD)*abs(ys(2)-ys(1));

Thank you, :)