Tiny symmetry mismatch, unsure where it comes from.

[AnderBiguri]

Expected Behavior

Projections of a symetric phantom need to produce a symetric projection (for certain geometries)

Actual Behavior

A tiny tiny tiny numerical offset exists. Much smaller than half of a pixel.

Code to reproduce the problem (If applicable)

     geo=defaultGeometry('nDetector',[512,512],'nVoxel',[512,512,512]);
    geo.dDetector=[0.280;0.280]*4;
    geo.sDetector=geo.dDetector.*geo.nDetector;
    geo.dVoxel=[1,1,1]';
    geo.sVoxel=geo.nVoxel.*geo.dVoxel;
    
    geo.DSO=600;
    geo.DSD=900;


        sphere=@(x,y,z,r,x0,y0,z0)((x-x0).^2+(y-y0).^2+(z-z0).^2<=r.^2);
        n_spheres_per_plane=6;
        
        % define locations of spheres and radious (in pixel units)
        r_sphere=5;
        l=100*geo.dVoxel(1); %for later
        z0=[100+geo.nVoxel(3)/2 -100+geo.nVoxel(3)/2+1]; %positive/negative
        
        t = linspace(0,2*pi,n_spheres_per_plane+1); t(end)=[];
        x0 = r.*cos(t)+geo.nVoxel(1)/2+0.5;
        y0 = r.*sin(t)+geo.nVoxel(2)/2+0.5;
        
        % get coordinates of all pixels
        [x,y,z]=meshgrid(1:geo.nVoxel(1),1:geo.nVoxel(2),1:geo.nVoxel(3));
        
        % crate all spheres
        sample=zeros(geo.nVoxel','single');
        for ii=1:n_spheres_per_plane
            sample=sample + sphere(x,y,z,r_sphere,x0(ii),y0(ii),z0(1))*(1+single(ii==1));
            sample=sample + sphere(x,y,z,r_sphere,x0(ii),y0(ii),z0(2))*(1+single(ii==1));
        end
    end


%%
        geo.offDetector=[30;0].*geo.dDetector;
        angles=0;
        proj=Ax(sample,geo,angles,'interpolated'); % happens for both projectors

       % symmetry check:
       imshow([proj(1:end/2,:,1)-flipud(proj(end/2+1:end,:,1))],[]);colorbar
       % notice the error is a shift upwards of the sample. its about 0.5% of error in pixel value. 
       % for this geometry, this number is the minimum error one I found:
       geo.offDetector=[0;-0.001236]; % this is the offset that causes the best value. 
       %Still can be seen in the projection, but at this stage we can attribute it to numerical error. However, I have no idea where that number comes from and why its consistent between projectors. I don't see a relationship of the number with the geometry.

[AnderBiguri]

A better code to reproduce it:

%% Initialize
clear;
close all;
%% Geometry
geo=defaultGeometry();
%% Define angles of projection and load phatom image

% define projection angles (in radians)
angles=0;



%% left right symetry
head=headPhantom(geo.nVoxel);
head=permute(head,[2,1,3]);
head(:,1:end/2,:)=fliplr(head(:,end/2+1:end,:));

proj_inter=Ax(head,geo,angles,'interpolated');
proj_siddon=Ax(head,geo,angles,'Siddon');

figure(1)
test=[proj_inter(:,1:end/2,1)-fliplr(proj_inter(:,end/2+1:end,1))]./proj_inter(:,1:end/2,1);
test(~isfinite(test))=0;
imshow(test,[]);colorbar;colormap(redblue)
title("Interpolated projection left-right symetry break in %")
caxis([-0.01,0.01])
figure(2)
test=[proj_siddon(:,1:end/2,1)-fliplr(proj_siddon(:,end/2+1:end,1))]./proj_siddon(:,1:end/2,1);
test(isnan(test))=0;
imshow(test,[]);colorbar;colormap(redblue)
title("Siddon projection left-right symetry break in %")
caxis([-0.01,0.01])

%% top/bottom symetry:
head=headPhantom(geo.nVoxel);
head=permute(head,[2,1,3]);
head(:,:,end/2+1:end)=flip(head(:,:,1:end/2),3);

proj_inter=Ax(head,geo,angles,'interpolated');
proj_siddon=Ax(head,geo,angles,'Siddon');

figure(1)
test=[proj_inter(1:end/2,:,1)-flipud(proj_inter(end/2+1:end,:,1))]./proj_inter(1:end/2,:,1);
test(~isfinite(test))=0;
imshow(test,[]);colorbar;colormap(redblue)
title("Interpolated projection top-bottom symetry break in %")
caxis([-0.01,0.01])
figure(2)
test=[proj_siddon(1:end/2,:,1)-flipud(proj_siddon(end/2+1:end,:,1))]./proj_siddon(1:end/2,:,1);
test(isnan(test))=0;
imshow(test,[]);colorbar;colormap(redblue)
title("Siddon projection top-bottom symetry break in %")
caxis([-0.01,0.01])

It could be numerical error, but it kinda looks like a numerical error shifted in one direction. Still unsure where this comes from....

[yliu88au]

Similar bug to #335 in projection, I wonder?

[AnderBiguri]

@yliu88au no, not at all unfortunately. #335 is just a weight, while this mismatch is much smaller, yet it happens in the forward proejction, so its not weighted. I think its numerical errors, but not sure.

[yliu88au]

For python version, there are some compilation warning for projection.cpp:
"projection.cpp
../Common/CUDA/projection.cpp(15): warning C4244: 'return': conversion from 'double' to 'float', possible loss of data
../Common/CUDA/projection.cpp(26): warning C4244: '=': conversion from 'double' to 'float', possible loss of data
../Common/CUDA/projection.cpp(30): warning C4244: '=': conversion from 'double' to 'float', possible loss of data
../Common/CUDA/projection.cpp(34): warning C4244: '=': conversion from 'double' to 'float', possible loss of data"

Could these be the cause?

[AnderBiguri]

@yliu88au no still not. In theory all the code should work with floats. But perhaps that makes all create this numerical error. I guess a test that I should carry is to recompile the entire projector with doubles, and see if the effect is still there. That would, in theory, discriminate if the cause is numerical errors or an actual bug.

[AnderBiguri]

This seems to come from floating point arithmetics. Someone proposed a way to solve it:

In the geometric transform code of the GPU, do not translate directional vectors, only rotate them. It seems to mitigate most of the numerical issues.