Closed AnderBiguri closed 2 months ago
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....
Similar bug to #335 in projection, I wonder?
@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.
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?
@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.
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.
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)