LAVI-USP / DBT-Reconstruction

LAVI open-source reconstruction toolbox for digital breast tomosynthesis (DBT)
GNU General Public License v3.0
54 stars 19 forks source link

Help needed in setting geometry parameters #5

Open roshtha opened 4 years ago

roshtha commented 4 years ago

Hello,

I run FBP algorithm for reconstructing images from DBT breast phantom raw data. The reconstructed images look blurry and lack details when compared with image outputs from actual DBT system. The actual number of slices is 171. Here I am reconstructing only 16 slices and parameter.dz was modified accordingly.
I believe I have some serious issues with geometry parameter settings. Could you please have a look at the following geometry definition and comment on any flaws. I attach a pictorial representation for the geometry of DBT system.

detector

The geometry definition follows.

% Breast voxels density
parameter.nx = 2816;    % number of voxels (columns)
parameter.ny = 3584;    % number of voxels (rows)
parameter.nz = 16;     % actual number of slices=171   number of voxels (slices)

% Detector panel pixel density
parameter.nu = 2816;    % number of pixels (columns) 
parameter.nv = 3584;    % number of pixels (rows) 

% Single voxel real size (mm)
parameter.dx = 0.25;    % on the X axis (mm)
parameter.dy = 0.25;    % on the Y axis (mm)
parameter.dz = 2.5;      % 0.25;  on the Z axis (mm)

% Single detector real size (mm)
parameter.du = 0.25;    % on the X axis (mm)
parameter.dv = 0.25;   % on the Y axis (mm)

% X-ray source and detector distances
parameter.DSD = 650;         % Distance from source to detector (mm)
parameter.DSO = 593;         % Distance from source to the top of object (mm)
parameter.DDR = 17;          % Distance from detector to pivot (mm)
parameter.DSR = parameter.DSD - parameter.DDR;  % Distance from source to pivot (mm)
parameter.DAG = 15;                             % Distance of Air Gap (mm)

% Number of Projections
parameter.nProj = 16;  

% Angle settings (Degrees)
parameter.tubeAngle = 15;   % Tube Angle
parameter.tubeDeg = linspace(-parameter.tubeAngle/2,parameter.tubeAngle/2,parameter.nProj);

parameter.detAngle = 0;   % Detector Angle
parameter.detectorDeg = linspace(-parameter.detAngle/2,parameter.detAngle/2,parameter.nProj);

Raw data matrix shared at https://drive.google.com/drive/folders/10er-5LklkIFZc3BR0Z_AKhpyFiaivTIC?usp=sharing

I also attach reconstructed image outputs from code and machine.

code_output.zip machine_output.zip

It would be very helpful if you can provide some inputs on this.

Thanks.

rodrigovimieiro commented 4 years ago

Hi, thanks for posting this issue.

The problem was occurring with your data. There were some zero values in the edges of the projections and were causing inf values after log. Make sure you don't have inf or NaN after log. This is very common with the detectors in the edges. I cropped around 10 pixels from both the beginning and the end of the projections rows:

dataProj = dataProj(10:end-10,:,:);

Also, as you sent the data through .mat, I manually set the PixelDepth to 16 in order to the pre-processing code works.

parameter.bitDepth = 16;

In the Parameter settings, some line codes were missing so I added them for you, and I also correctly set the number of slices to:

parameter.nz = 171;

If you set it to 16, it would only reconstruct the first 16 slices. If you want to save only a certain slice, you might want to set parameter.sliceRange correctly.

I also corrected some geometry problems. Please, check carefully the values that I set:

% Breast voxels density
parameter.nx = 1317;%2816;    % number of voxels (columns)
parameter.ny = 3584;    % number of voxels (rows)
parameter.nz = 171;     % actual number of slices=171   number of voxels (slices)

% Detector panel pixel density
parameter.nu = 2816;    % number of pixels (columns) 
parameter.nv = 3584;    % number of pixels (rows) 

% Single voxel real size (mm)
parameter.dx = 0.083;    % on the X axis (mm)
parameter.dy = 0.083;    % on the Y axis (mm)
parameter.dz = 0.25;      % 0.25;  on the Z axis (mm)

% Single detector real size (mm)
parameter.du = 0.083;    % on the X axis (mm)
parameter.dv = 0.083;   % on the Y axis (mm)

% X-ray source and detector distances
parameter.DSD = 650;         % Distance from source to detector (mm)
parameter.DSO = 593;         % Distance from source to the top of object (mm)
parameter.DDR = 88;          % Distance from detector to pivot (mm)
parameter.DSR = parameter.DSD - parameter.DDR;  % Distance from source to pivot (mm)
parameter.DAG = 17;                             % Distance of Air Gap (mm)

% Detector and object full real sizes (mm)
parameter.sx = parameter.nx.*parameter.dx;  
parameter.sy = parameter.ny.*parameter.dy;  
parameter.sz = (parameter.nz.*parameter.dz)+parameter.DAG;  
parameter.su = parameter.nu.*parameter.du;  
parameter.sv = parameter.nv.*parameter.dv;  

% Detector and object Volume grid settings 
parameter.xs = (parameter.nx-1:-1:0)*parameter.dx;
parameter.ys = (-(parameter.ny-1)/2:1:(parameter.ny-1)/2)*parameter.dy;
parameter.zs = (0:1:parameter.nz-1)*parameter.dz + parameter.DAG;
parameter.us = (parameter.nu-1:-1:0)*parameter.du;
parameter.vs = (-(parameter.nv-1)/2:1:(parameter.nv-1)/2)*parameter.dv;

% Number of Projections
parameter.nProj = 16;  

% Angle settings (Degrees)
parameter.tubeAngle = 15;   % Tube Angle
parameter.tubeDeg = linspace(-parameter.tubeAngle/2,parameter.tubeAngle/2,parameter.nProj);

parameter.detAngle = 0;   % Detector Angle
parameter.detectorDeg = linspace(-parameter.detAngle/2,parameter.detAngle/2,parameter.nProj);

% Slice range to be saved
parameter.sliceRange = 1:parameter.nz; 
% Region of interest (ROI) to store
parameter.iROI = 1:parameter.ny;     
parameter.jROI = 1:parameter.nx;  

% Bit number quatization
parameter.bitDepth = 16;     % Load from dicom header

Screen Shot 2020-04-29 at 15 19 16

There are still errors in the reconstruction because de MC is not getting on focus anywhere. I suspect they are from the angles. The projections angles are linearly spaced? Double check your geometry and let me know.

I am sending you the results that I had through this link.

Let me know how it goes.

Best.

roshtha commented 4 years ago

Thanks for the code and reply. I made modifications and run the code. Got images similar to yours. The projection angles are linearly spaced spanning from -7.5 to +7.5. The tube motion is continuous. The FBP reconstruction takes around 40 minutes for 171 slices. I could not find cuda version code for projection and backprojection scripts. I think backprojectionDDb_mex_CUDA and projectionDDb_mex_CUDA are for rotating detectors. Ours is a static detector.

rodrigovimieiro commented 4 years ago

Hi,

Great. However, I still think there are errors in the geometry. You might want to use animation = uint8(1); to see your geometry.

The source codes from the CUDA acceleration are in the Functions/Sources folder. Please, follow the instructions to compile those codes. If you are under Windows, or Linux.

Please, let me know if you face any errors. I can provide some binaries for Linux (Ubuntu) if you can't compile them.

Best.