artifacts appeared using simple fdk

[JMLiu-SEU]

Expected Behavior

I made an ideal spherical shell phantom (3D tomography) and using Ax to generate projections and use FDK to reconstruct. the geometry are set to fixed ratio of default geometry and the geo.accuracy is set to 0.5. the num of projs are 400, I think is sufficient for fdk.

Actual Behavior

The reconstructed images are added two paralleal streak artifacts compared with phantom I made, do not know why. just copy the code below into a .m file into tigre root folder will reproduce the result

Code to reproduce the problem (If applicable)

%code of generate phantom: size are 512*512*512
N = 512;
R1 = 200;
R2 = 180;
[X,Y,Z] = meshgrid(-N/2+0.5:N/2-0.5);

D = sqrt(X.^2+Y.^2+Z.^2);
D = R1+sqrt(3)/2-D;
D(D<0) = 0;
D(D>sqrt(3)) = sqrt(3);
D = D/sqrt(3);
volumedata = D;

D = sqrt(X.^2+Y.^2+Z.^2);
D = R2+sqrt(3)/2-D;
D(D<0) = 0;
D(D>sqrt(3)) = sqrt(3);
D = D/sqrt(3);

volumedata = volumedata - D;

m_min = min(volumedata,'','all');
m_max = max(volumedata,'','all');
mkdir('ball');
for i = 1:N
    filename = sprintf('%.3d.dcm',i);
    R = uint16((volumedata(:,:,i) - m_min) / (m_max - m_min) * 16383 * 0.4);
    meta.Modality = 'Micro CT';
    meta.SeriesInstanceUID = '1.2.0.20190424.115306.20190424.115251-sio2-9';
    meta.InstanceNumber = i-1;
    meta.PatientID = sprintf('%.4d',i);
    dicomwrite(R, ['ball\',filename], meta);
end

% code of making projections
geo.DSD = 8000*4*3.375e-3;                         % Distance Source Detector      (mm)
geo.DSO = geo.DSD/2;                             % Distance Source Origin        (mm)
% Detector parameters
geo.nDetector=[512; 512];                        % number of pixels              (px)
geo.dDetector=[4*3.375e-3; 4*3.375e-3];          % size of each pixel            (mm)
geo.sDetector=geo.nDetector.*geo.dDetector;      % total size of the detector    (mm)
% Image parameters
geo.nVoxel=[512;512;512];                   % number of voxels              (vx)
geo.sVoxel=geo.nVoxel*4*1.6875e-3;          % total size of the image       (mm)
geo.dVoxel=geo.sVoxel./geo.nVoxel;          % size of each voxel            (mm)
% Offsets
geo.offOrigin =[0;0;0];                     % Offset of image from origin   (mm)              
% geo.offOrigin =[0;10*4*1.6875e-3;0];                     % Offset of image from origin   (mm)              
% geo.offOrigin =[10*4*1.6875e-3;0;0];                     % Offset of image from origin   (mm)              
% geo.offOrigin =[0;0;0];                     % Offset of image from origin   (mm)              
geo.offDetector=[0;0];                     % Offset of Detector            (mm)
                                            % These two can be also defined
                                            % per angle
% Auxiliary 
geo.accuracy=0.5;                           % Variable to define accuracy of
                                            % 'interpolated' projection
                                            % It defines the amoutn of
                                            % samples per voxel.
                                            % Recommended <=0.5             (vx/sample)
% Optional Parameters
% There is no need to define these unless you actually need them in your
% reconstruction                   
geo.COR=0;                                  % y direction displacement for 
                                            % centre of rotation
                                            % correction                   (mm)
                                            % This can also be defined per
                                            % angle
geo.rotDetector=[0;0;0];                    % Rotation of the detector, by 
                                            % X,Y and Z axis respectively. (rad)
                                            % This can also be defined per
                                            % angle       
geo.mode='cone';                            % Or 'parallel'. Geometry type.  

%% Define angles of projection and load phatom image
% define projection angles (in radians)
angles = linspace(0,2*pi,400);
file = '.\ball\';
filefolder = dir([file,'*.dcm']);
imagehei = 512;
imagewid = 512;
imagesize = imagehei * imagewid;
imagenum = length(filefolder);
volumedata = zeros(imagehei,imagewid,imagenum);
for i = 1:imagenum
    filename = filefolder(i).name;
    volumedata(:,:,i) = dicomread([file filename]);
end
volumedata = single(volumedata);

projections = Ax(volumedata,geo,angles,'interpolated');
%% Save projections as DICOM file
savedcm = projections;
file =  ['.\ball_proj_',sprintf('%.2f',geo.offOrigin(1)*1e3),'_y_',sprintf('%.2f',geo.offOrigin(2)*1e3),'_z_',sprintf('%.2f',geo.offOrigin(3)*1e3),'\'];
mkdir(file);
imagenum = size(savedcm,3);
m_min = min(savedcm,'','all');
m_max = max(savedcm,'','all');
for i=1:imagenum
    filename = sprintf('%.4d.dcm',i);
    R = uint16((savedcm(:,:,i) - m_min) / (m_max - m_min) * 16383);
    meta.Modality = 'Micro CT';
    meta.InstanceNumber = i-1;
    meta.PatientID = sprintf('%.4d',i);
    dicomwrite(R, [file,filename], meta);
end

%code of reconstruct
file = '.\ball_proj_0.00_y_0.00_z_0.00\';
filefolder = dir([file,'*.dcm']);
imagehei = 512;
imagewid = 512;
imagesize = imagehei * imagewid;
imagenum = length(filefolder);
ProjData = zeros(imagehei,imagewid,imagenum);
for i = 1:imagenum
    filename = filefolder(i).name;
    ProjData(:,:,i) = dicomread([file filename]);
end
ProjData = single(ProjData);

%% Define angles of projection
angles=linspace(0,2*pi,400);
%% Reconstruct image using FDK
imgFDK=FDK(ProjData,geo,angles);
imshow(imgFDK(:,:,200),[]);

Specifications

• MATLAB version: r2020a
• OS: windows10
• CUDA version: cuda10.0

[JMLiu-SEU]

the left side and right side are two parallel artifacts more explicit if adjust the window level

[AnderBiguri]

anlges=linspace(0,2*pi,400) repeats 1 angle (0==2pi) so you got extra artefacts in that direction. If you just make `angles=linspace(0,2pi,401);angles(end)=[]or the equivalentangles=linspace(0,2pi-2pi/400,400);` they go away.

[JMLiu-SEU]

thankyou so much