Issue with Helical Reconstruction Geometry Parameters in Openrecon

Hello,

I'm working on helical reconstruction using Openrecon and trying to match the reconstruction geometry parameters to those provided in Leap. Here is my code:

SO = 110;          % The distance from source to object (cm)
DO = 30;           % The distance from detector to object center(cm)
YL = 512;          % Detector cell number along the horizontal direction of detector array
YLC= (YL-1)*0.5;   % Detector center along the horizontal direction of detector array
ZL = 128;          % Detector cell number along the vertical direction of detector array
ZLC= (ZL-1)*0.5;   % Detector center along the vertical direction of detector array
DecAngle = 0.2388; % Detector array fan beam angle along the horizontal direction (rad)0.2388
DecHeigh = 8.32;   % Detector array height along the vertical direction (cm)
N_Turn   = 10;     % The number of turns for the whole helical scan 
N_2pi    = 720;    % The projections/views number for each turn of scan
h        = 1;      % Helical pitch related to detector height
%The following line 25-28 are used to define the reconstruction paramters
WtInd    = 3;      %1:Redundancy weighting; 2:2D weigthing; 3: 3D weighting 
k1       = 5;      % The order to define the 3D weighting function
delta    = 60;     % The range to define smoothness of 2D weigthing function 
HSCoef   = 0.6;    % This is used to define the half-scan range              

ObjR  = 13.7;      % Diameter of imaginh object or radius of FOV
SSO   = 512 ;     % Define the size of the reconstructed image
XN    = SSO;
YN    = SSO;
ZN    = SSO;
XC    = (XN-1)*0.5;
YC    = (YN-1)*0.5;
ZC    = (ZN-1)*0.5;
BetaS = -N_Turn*pi;
BetaE =  N_Turn*pi;
ViewN =  N_Turn*N_2pi;

DecWidth = tan(DecAngle*0.5)*(SO+DO)*2;
dYL   =  DecWidth/YL;
dZL   =  DecHeigh/ZL;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ScanGeom = struct( 'ScDet',  [SO DO YL ZL DecAngle DecHeigh YLC ZLC h],  ...
                   'Proj',   [BetaS BetaE ViewN N_2pi], ...
                   'Obj',    [ObjR XN YN ZN XC YC ZC], ...
                   'Rec',    [delta HSCoef k1]);            
% Generating projections or load the projections
disp('Generating the cone beam projections');
dataStruct = load('LEAP_PROJECTION.mat');% B: Spectrum of x-ray source
disp(fieldnames(dataStruct));
images = dataStruct.nifti_data;
swapped_data_array = permute(images, [1, 3, 2]);
Proj=swapped_data_array;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rebinning conebeam to parallel beam 
disp('Paralllel rebinning from the conebeam projections');
Proj = Parallel_Rebinning_CB_Curve(Proj,ScanGeom);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Perform Ramp filtering
disp('Filtering the cone beam projection');
n = -YL:YL;
hsl=-2./((pi^2)*(4*n.*n-1))/dYL;
%plot(n,hsl); 
NN = 2^ceil(log(YL*3)/log(2));
HS = zeros(1,NN);
HS(1:YL+1)= hsl(YL+1:2*YL+1);
HS(end-YL+1:end)=hsl(1:YL);
FtS = fft(HS);
fpwd = zeros(size(Proj));
for i=1:ViewN
    for k=1:ZL
     FtP = fft(Proj(i,:,k),NN);
     tep = ifft(FtS.*FtP);
     fpwd(i,:,k) = real(tep(1:YL)); 
    end
end
save leap_sinogram_ramp_zeropad.mat fpwd;

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Backproject the filtered data into the 3D space
ScanGeom.ScDet(5) = DecWidth;
disp('BackProjection the filtered projection');
ctimage = zeros(XN,YN,ZN);
switch(WtInd)
    case 1 %1/2 redundancy weigthing
        Parallel_FDK_Helical_NOWeighting(ScanGeom,fpwd,ctimage);
        save RecResNoW ctimage;
    case 2 % 2D weighting
        Parallel_FDK_Helical_2DWeighting(ScanGeom,fpwd,ctimage);
        save RecRes2DW ctimage;
    case 3 %3D weighting
        Parallel_FDK_Helical_3DWeighting(ScanGeom,fpwd,ctimage);
        save RecRes3DW ctimage;
    otherwise %1/2 redundancy weigthing
        Parallel_FDK_Helical_NOWeighting(ScanGeom,fpwd,ctimage);
        save RecResNoW ctimage;
end

volumeViewer(ctimage); % or volshow(ctimage);

% Calculate and display elapsed time
time = cputime - time;
disp(['Elapsed time: ', num2str(time), ' seconds']);

% Define the output folder
output_folder = 'ct_images_leap1';
if ~exist(output_folder, 'dir')
    mkdir(output_folder);
end

processed_images = zeros(512, 512);

% Save each slice as a PNG image
for z = 1:ZN
    slice = ctimage(:, :, z);
    max_ = max(slice(:));
    min_ = min(slice(:));
    slice = (slice - min_) / (max_ - min_);

    filename = fullfile(output_folder, sprintf('slice_%03d.png', z));
    imwrite(slice, filename);
end

I used the geometry from Leap and set the values from this script However, my reconstruction results do not match my expectations. Here is what I get:

Slice 001: slice_001

Slice 082: slice_082

Slice 227: slice_227

I was expecting a more accurate reconstruction. Could you please help me identify what might be going wrong?Any insights or suggestions would be greatly appreciated. Thank you!

Best regards,

xcist / main

Issue with Helical Reconstruction Geometry Parameters in Openrecon #78