Artefacts in SART and OS_SART reconstruction for DBT Scanner

[pathros123]

Expected Behavior

I want to reconstruct a DBT data using 16 projection data into 48 images

Actual Behavior

I tried FDK, OS_SART and SART but some king of duplicating structures are present in SART and OS_SART reconstructed images. Surprisingly these structures are considerably less in FDK but the contrast of output image is very poor. I have added all three output images below for your reference.

I used FDK to initialise SART and OS_SART.
Number of iterations : 8
Block size : 8

Code to reproduce the problem (If applicable)

By examining the SART and OS_SART code, i found that for each iteration, a duplicate of these structures are added to the image. The line is given below

The input res and res after one iteration is shown in the image

  if nesterov
            % The nesterov update is quite similar to the normal update, it
            % just uses this update, plus part of the last one.
            ynesterov=res+ bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,index_angles(:,jj)).*(proj(:,:,index_angles(:,jj))-Ax(res,geo,angles_reorder(:,jj),'gpuids',gpuids)),geo,angles_reorder(:,jj),'gpuids',gpuids));
            res=(1-gamma)*ynesterov+gamma*ynesterov_prev;
        else
            res=res+lambda* bsxfun(@times,1./V(:,:,jj),Atb(W(:,:,index_angles(:,jj)).*(proj(:,:,index_angles(:,jj))-Ax(res,geo,angles_reorder(:,jj),'gpuids',gpuids)),geo,angles_reorder(:,jj),'gpuids',gpuids));
        end

Specifications

• MATLAB version: 2021b
• OS: windows 10
• CUDA version: 10.0

[AnderBiguri]

Thanks! Can you also test SIRT, or test SART and OS-SART with the angle reordering turned off?

Also,are you displaying the same slices? As you know, DBT has very limited depth resolution and artifacts look similar to what you show for SART when visualizing slices out of the focal point, for all algorithms

[pathros123]

Thankyou for the quick response! We tried SIRT and OS_SART in order strategies 'ordered' and 'random'. But we got the same artefacts!

Yes. I am displaying the same slices.

Regarding your last comment,

1. Among the 48 reconstructed output images, only slices 19, 20, 21 shows output without these artefacts.
2. If this is due to visualising slices out of the focal point, then why FDK algorithm doesn't show any artefacts?
3. I tried reducing the number of slices for reconstruction from 48 to 40. Still the issue persists.

[AnderBiguri]

1- Indeed, the focal point, likely (depends on the geometry, but you'd need to share that).
2- Why does FDK not show them?? Surprised to be fair, it should. Admitedly sometimes iterative algorithms go bad when you overiterate in extremely bad adquisition settings, like DBT
3- Its not a thing you can change by changing the reconstruction slices, its caused by the acquisition geometry.

In order for me to verify this is not a bug, can you provide a compete demo of this, with fake data generation in TIGRE, such that I Can copy paste and investigate? If its an issue of TIGRE, and not the data/acquisition, then you should be able to reproduce it without the real data, using e.g. shepp-logan. Can you try that?

[pathros123]

The reason for FDK not showing this error maybe due to the poor contrast of output images.

We will try this using shepp-logan and will update you. Thankyou

[rupikaraj]

Were you able to reach any conclusion?

[AnderBiguri]

Hi @rupikaraj. I am quite busy these days so I have limited time to investigate. It would help if you would post a minimal example reproducing your issue that does not use your data, but tigres geometry. That is much easier for me (and others) to help with.

[pathros123]

@AnderBiguri I tried shepp-logan phantom on the d19_DBT_example original code using the default parameters and got the same artefacts. I have added some of the reconstructed images below.

Please note that i have not changed any other parameters from the default code.

[AnderBiguri]

Thanks! Can you also post the code to generate that? It's not good security advice to download unknown person drive files, sorry.

[pathros123]

I used the d19_DBTexample.m code in MATLAB/Demos in this repository. Didn't change any other parameters. Just changed the head phantom to shepp-logan.

`shepp = phantom3dAniso(geo.nvoxel) % test data creation function

shepp = single(shepp) `

[AnderBiguri]

Hi @pathros123 . This is supsicius, in the sense that I think its not showing the same effect. I think the shepp logan example is showing the errors of the discretization of the phantom, not errors in algorithms. In general, if you change the data and the effect is not there anymore/appears suddenly, its the data issue, not the algorithm issue.

I will investigate, but still unsure if this is a bug on the code, or just you have some issue with the geometry/angles in your real case, as I am not 100% sure this demo is showing the same effect as your real data.

[pathros123]

@AnderBiguri I did some research but couldn't find any similar issues in DBT reconstruction. In my data and in shepp-logan phantom the error found in the image seemed to be similar for me. I am pretty new to this field. Please check and give your response if it is an issue with the data.

[AnderBiguri]

@pathros123 can you show me slices in the source detector axis for both cases?

[pathros123]

I'm not sure what 'slices in source detector axis' means. Could you please clarify your question?

[AnderBiguri]

@pathros123 These are 3D volumes, can you show me cuts of the volumes in different axis?

[pathros123]

Here are the source to detector axis images of last slices in both shepp-logan and my data.

[AnderBiguri]

hi @pathros123 , apologies, can you show me images of the 3 axis? Not sure how you made these, but they should be the TIGRE XY axis, i.e. Z should not appear in the plot. Can you show a slice of all 3 cuts?

[pathros123]

@AnderBiguri I have added the code below for shepp-logan phantom.
`clc;clear;close all

% Distances
geo.DSD = 660; % Distance Source Detector (mm)
geo.DSO = 620; % Distance Source Origin (mm)
% Detector parameters
geo.nDetector=[3064;2396]/2; % number of pixels (px)
geo.dDetector=[0.1; 0.1]; % size of each pixel (mm)
geo.sDetector=geo.nDetector.*geo.dDetector; % total size of the detector (mm)
% Image parameters
geo.nVoxel=[128;2054;784]/2; % number of voxels (vx)
geo.sVoxel=[63.3;205.3;78.4]/2; % total size of the image (mm)
geo.dVoxel=geo.sVoxel./geo.nVoxel; % size of each voxel (mm)
% Offsets
Airgap = 22; % DBT airgap (mm)
geo.offOrigin =[((geo.sVoxel(1)/2)-...
(geo.DSD-geo.DSO)+Airgap);0;geo.sVoxel(3)/2]; % Offset of image from origin (mm)
geo.offDetector=[0; geo.sDetector(2)/2]; % Offset of Detector (mm)

geo.accuracy=0.5; % Variable to define accuracy of

geo.COR=0; % y direction displacement for

nprojs = 9; % Number of projections
tubeangle = 25; % Angle range
angles=deg2rad(linspace(-tubeangle/2,tubeangle/2,nprojs)); % Angles in rad
geo.rotDetector=[0;0;0]; % Rotation of the detector, by

geo.mode='cone'; % Or 'parallel'. Geometry type.

%% Adapt CT geo to DBT
geo=staticDetectorGeo(geo,angles);
%% Adapt DBT projections to TIGRE CT projections

% head=headPhantom(geo.nVoxel);

shepp=phantom3dAniso(geo.nVoxel);
shepp=single(shepp);

projections=Ax(shepp,geo,angles,'interpolated');

%% Simple BP
recFDK = FDK( projections,geo,angles);

%% SART
recSART = SART(projections,geo,angles,2,'OrderStrategy','ordered');

%% plot

plotImg(recSART,'dim','z','clims',[0 1], 'slice', 100)

`

The plots for all three cuts are added below.

I tried to resize the shepp-logan to 128x128x128 size for obtaining a better view. Here are the images i got.

In the gammex phantom data that I am using, geo.nVoxel = 48,1150,903. Hence the XY axis image is a very narrow rectangle.

I hope this helps.

[AnderBiguri]

Sorry I am quite busy, I'll have a look soon :)