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.
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.
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,
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.
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.
Comment only