Code covered by the BSD License  

Highlights from
3D Cone beam CT (CBCT) projection backprojection FDK MLEM reconstruction MATLAB codes for students

4.66667

4.7 | 9 ratings Rate this file 106 Downloads (last 30 days) File Size: 232 KB File ID: #35548
image thumbnail

3D Cone beam CT (CBCT) projection backprojection FDK MLEM reconstruction MATLAB codes for students

by

 

10 Mar 2012 (Updated )

3D Cone beam CT (CBCT) projection backprojection MLEM FDK reconstruction MATLAB source codes

| Watch this File

File Information
Description

Hello,
My name is Kyungsang Kim.
This program is about 3D cone-beam CT.

This is made for students who learn the medical imaging.

Please let me know if the program has problems.
(E-mail: kssigari(at)gmail)

Thank you

Required Products MATLAB
MATLAB release MATLAB 7.8 (R2009a)
MATLAB Search Path
/
/bin
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (18)
12 Jun 2014 mat  
29 May 2014 Hamdi Mahdjoub  
28 May 2014 Bella

Would 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!

24 May 2014 Rama Teja A S

Please include analytical equation for better understanding. Or please explain the algorithm in projection.m. If possible please include analytical equation references

11 Nov 2013 Kyung Sang Kim

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

11 Nov 2013 erick

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

28 Oct 2013 Kyung Sang Kim

Hi 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,

28 Oct 2013 erick

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

27 Oct 2013 Kyungsang

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

27 Oct 2013 erick

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

27 Oct 2013 erick

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

04 Apr 2013 Kyung Sang Kim

Hi, Masoud Hashemi,
This program is for flat detectors,
I think you should modify for the arc detector geometry,
Thanks for using my software.

30 Mar 2013 Masoud Hashemi

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

02 Mar 2013 Tai Chieh

Great!

16 Dec 2012 ljpzgx

The 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!

13 Mar 2012 Kyung Sang Kim

This 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".

12 Mar 2012 Kyung Sang Kim

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

12 Mar 2012 Kyung Sang Kim

Please 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, :)

Updates
20 Mar 2012

I optimized the source codes

16 Nov 2012

Fixed a bug

13 Dec 2012

I changed title, tag

14 Jan 2013

Change description

24 Apr 2013

I modified codes simply.

17 Dec 2013

Update code and debug

28 Feb 2014

Debug and add some descriptions

25 Oct 2014

Fixed filtering function

Contact us