artifacts appeared using simple fdk
JMLiu-SEU opened this issue · 3 comments
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
anlges=linspace(0,2*pi,400)
repeats 1 angle (0==2*pi) so you got extra artefacts in that direction.
If you just make angles=linspace(0,2*pi,401);angles(end)=[]
or the equivalent angles=linspace(0,2*pi-2*pi/400,400);
they go away.
thankyou so much