File Exchange

image thumbnail

Cone beam CT simulation

version 1.0 (51.5 KB) by

To simulate a CBCT system by computing the 2D X-ray projection images.

7 Ratings



View License

Cone beam CT simulation

By: Deshan Yang, PhD
Department of radiation oncology, school of medicine, Washington University in Saint Louis

compute_projections.m - the main program, to compute all projection and save the 2D projection images into a folder
compute_one_projection.m - Called by the main program to compute only one projection image
straight_line_integral.m - The core function to compute a line integral through a 3D matrix
A lot of lines are commented out in this file because the file was originally written
in MATLAB only, and laterly rewritten into c for speed improvement.
straight_line_integral_inner.c - This is what rewritten in C from the straight_line_integral.m
sort_ts1_ts2.c - A utility function to combine two sorted array into one sorted array.
It is written in c to improve the speed

create_CSV_file_for_OSCaR.m - To create the CSV file which is used in OSCaR CBCT reconstruction
However, OSCaR does not support the saved mat data files without modification.
I have the modified OSCaR functions.

load_all_projection_images.m - Loading all projection images from a folder into a 3D volume
make_sure_positive.m - A small utility function

Binary files:
*.mexw32, *.mexw64 The compiled binary MEX files for 32-bit and 64-bit windows

Other files: - For different ways to compute projections. They are not in use anymore because
they are either much slower or less accurate

How to setup the 3D data:
You need to prepare an image 3D volume as the digital phantom (or patient, or object, whatever you call it).
For this main function: compute_projections(xs,ys,zs,data3d,mode,output_folder)
xs: the x coordinates of the image pixels in an array
ys: the y coordinates of the image pixels in an array
zs: the z coordinates of the image pixels in an array

xs, ys and zs should be centered at [0,0,0] which is the CBCT isocenter.

About computation performance:
The straight_line_integral function is partially written in c. It is very accurate and fairly fast.
The computation of projection implemented here is both faster and more accurate than the radon transformation code in MATLAB.
For a 3D matrix (512x512x176), it takes about 35 seconds to compute one 2D projection in resolution of 512x384 on a DELL with 2.33 GHz CPU.
Computation does not need a lot of RAM from my understanding.

Comments and Ratings (23)


akram (view profile)

hi Yang,
i have a 3d matrix of 3D sheep-logan, i want to apply the radon transform on it, so i wonder where can i find this function to call it!!
thank you for advance

mary noor

mary noor

please whate meaning the vailable data3d if its meaning image in 3d how can i call it pleas give me example


Ben (view profile)

The projection image looks wired. There are many empty lines that look like artifacts, as Samuel reported before. Anyone successfully went through the example code (posted by Deshan Yang at 09 May 2011)?


Rajesh (view profile)

Can I get forward Projections using this code. Please suggest.


Rajesh (view profile)

I have a volume data. I want forward Projections. Can I get it from this algorithm?

Yu Tang


LI (view profile)

and I wnat to try"% Method 2, using my own data". But I didn't see file named'my_3d_data.mat'.

Easley Hu

I found the Coordinate system and the unit uesed in this code are different from that used in OSCaR. Simply flip left/right also does't work. So many people have the problem how to use the
simulated data to reconstruct with OSCaR. Therefore, can you be kindly to give us a example to show how to solve this probem? Thanks very much!

zahid ansari

I have 360 dicom imges taken from each angle and i want to know how to make .csv file so that i can open it from OSCaR

Beyung Doo Jo

Sorry, I've just found your e-mail address.

Deshan Yang

Deshan Yang (view profile)

Beyung Doo Jo, as I stated previously, OSCAR won't work without modification. Please contact me in private if you want me to help. However, you really shall figure it out by yourself, by just fitting the data into OSCAR data format.

BTW, the projection images created by this code need to be flipped left/right to be compatible with OSCAR.

Beyung Doo Jo

Deshan Yang, i have a question for you, By using your create_CSV_file_for_OSCaR.m, i got CSV file. but when i used this file for OSCAr, it didn't work. How can i use this SCV file with OSCAR? Could you help me?


Samuel (view profile)


thank you very much for this work. :-))

I have installed it on my Linux Ubuntu 8.4 LTS x86_64 + Octave 3.2.
Of course I had to rebuild the mex files for linux and need to download the phantom3d.m file to have a testing phantom.

All this is done and you code seems to work. So I have a large set of projections.

The point now is that it does not seems to be the projections I expected. There is a lot of empty lines and columns on the projections, as if there were an aliasing issue when computing lines integrals or something like that.

I use your sample code in your "example post" up here, but I update the Shepp-Logan phantom/matrix to 128x128x128 voxels of 2x2x2 mm.

hum, well now I want to attach a screen capt of the picture and a copy of one projection_x_y.mat ... but I don't found the "attach file" button on this web page :-((

An other question,
X-Ray attenuation is a I=I0*exp(-mu*x) decrease function, with 'mu' specific for each biological tissu and each X-Ray intensity we use.
So my question is : do someone know what mean the values of the Shepp-Logan phantom in "phantom3d.m" ? values rank from -5.5511e-17 to 1.0, witch seems to be a bug ... as we cannot have negative 'mu' coef ! But we can have negative Hounsfiel units ... but it usually ranks from -1000 to +3000 UH ...

Thank you for any help.

Deshan Yang

Deshan Yang (view profile)

One more note here:
The sinogram images calculated by the simulation might need to be flipped left/right before being used into reconstruction, for instance, using OSCaR. This of course depends on your CBCT reconstruction algorithm implementation.

Easley Hu

It's perfect! Thank you very much!

judy love

thank you very much!

Deshan Yang

Deshan Yang (view profile)

Many people have questions about how to setup the phantom. Here are some examples:

% CBCT_simulation_demo.m
output_folder = 'c:\CBCT_output_folder'; % Set the output folder

% Loading 3D phantom
% Method 1, using the Shepp-Logan phantom
data3d = phantom3d(256); % data3d will be 256x256x256, assuming the pixel size is 1mmx1mmx1mm
xs = 1:256; % 1 mm pixel size
xs = xs-mean(xs); % Center at the [0,0,0]
ys = xs; % 1 mm pixel size
zs = xs; % 1 mm pixel size
compute_projections(xs,ys,zs,data3d,1,output_folder); % Compute the projections

% Method 2, using my own data
data3d = load('my_3d_data.mat'); % Load saved 3D matrix
dx = 1; dy = 1; dz = 1; % The pixel size, is 1mm x 1mm x 1mm, giving in mm, not cm
% Image pixel size should be modified according to the real physical pixel dimension
dim = size(data3d);
xs = ((1:dim(2))-1)*dx; % Set xs
xs = xs - mean(xs); % Center to x=0
ys = ((1:dim(1))-1)*dy; % Set ys
ys = ys - mean(ys); % Center to y = 0
zs = ((1:dim(3))-1)*dz; % Set zs
zs = zs - mean(zs); % Center to z = 0

compute_projections(xs,ys,zs,data3d,1,output_folder); % Compute the projections


Tian (view profile)

Thank you very much for your great job.
I have a misunderstanding when running the program.
The result have a lot of NaNs. And I dont know if it is the problem I set the parameters.
Could you please upload the your sample data and how do you set the parameter?
Thank you very much for your help.

Easley Hu

Thank you for your share!
I have a problem when running these codes with my own 3D data.I found the d=[dx dy dz] have the value of [1 0 0], that means p1 and p2 can be –Inf or Inf. So there is an error on "ny = ceil(p1(2)):floor(p2(2));".
In order to run this toolbox correctly would you please give me a simple example of the 3D data ?

Deshan Yang

Deshan Yang (view profile)


Good catch. The "SDD" is the detector-to-isocenter distance. It should have been called something like IDD.


Great job! But I have one misunderstanding. Infile compute_projections.m I read:
SDD = 500; % Detector to axis distance = 50 cm
SAD = 1000; % X-ray source to axis distance = 100 cm,
but how SAD (source-axis distance) may be greater than SDD (source-detector distance) if axis is between the source and the detector?


MATLAB Release
MATLAB 7.6 (R2008a)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video

Win prizes and improve your MATLAB skills

Play today