Python fails at 2D tomography

[AnderBiguri]

The following code will return projection data that are all zeros.

from __future__ import division, print_function
from __future__ import division

import tigre
import tigre.algorithms as algs
import numpy as np
from tigre.demos.Test_data import data_loader
from tigre.utilities.Measure_Quality import Measure_Quality
from matplotlib import pyplot as plt

#geo1 = tigre.geometry(mode='cone', high_quality=False, default=True)
geo = tigre.geometry(mode='cone')
geo.DSD = 1535
geo.DSO = 1000
geo.nVoxel = np.array([1, 256, 256])
geo.sVoxel = geo.nVoxel * 1.0
geo.dVoxel = geo.nVoxel / geo.sVoxel
geo.dDetector = np.array([0.8, 0.8])*2               # size of each pixel            (mm)
geo.nDetector = np.array([1, 256])     # has to be true 3D for python version
geo.sDetector = geo.dDetector * geo.nDetector
geo.offOrigin = np.array([0, 0, 0])
geo.offDetector = np.array([0.0, 0.0]) # an known bug for v1.5
geo.rotDetector = np.array([0, 0, 0])
geo.accuracy=0.5

niter = 20
nangles = 100
angles = np.linspace(0, 2 * np.pi, nangles, dtype=np.float32)
head = data_loader.load_head_phantom(geo.nVoxel)
# test_data = scipy.io.loadmat('head.mat')
# data = test_data['img'].copy()
#head = data[:,:,51:53].transpose(2,1,0).copy()
#head = np.expand_dims(data[100,:,:].copy(),axis=0)

proj = tigre.Ax(head,geo,angles)

plt.figure()
plt.imshow(proj[50,:,:])
plt.show()
exit()

Debugging test made:

• Changing geo.nVoxel to geo.nVoxel = np.array([4, 256, 256]) will make the projections be non-zero.

• I was worried that offsets may be wrong, but changing geo.nDetector = np.array([10,1]) shows that the image is well aligned to the center of the detector.

[genusn]

@AnderBiguri , Something is happening in siddon?

proj = tigre.Ax(head,geo,angles, krylov="interpolated")
valMin = np.min(proj)
valMax = np.max(proj)
print(valMin, valMax)

showed

0.0 80.21488

The comment https://github.com/CERN/TIGRE/pull/189#issuecomment-725458261 might be related to this?

[AnderBiguri]

@genusn hmmm I am not sure if the comment is related. In that comment, I found out that the problem with the base was that it did not take into account half a pixel in locating the detector in space, thus when projecting a cube, a symmetry was missing. But your PR instead had that symmetry good (i.e. half top flipped was the same as half-bottom of detector) so if anything, the PR should have "solved" the problem we see here.

I suspect it may have to do with odd detector voxels, or something like that (possibly related #208). Still not sure, I need to find time to dig into this (too many things on my plate!!)

[genusn]

so if anything, the PR should have "solved" the problem we see here.

I see.

How about proj = tigre.Ax(head,geo,angles, krylov="interpolated")? The result was non zero.

[AnderBiguri]

I don't have much time to test now, but its a big iffy to test this issue with interpolated. Interpolated will also interpolate outside the image (up to 0.5 pixels outside the boundary of the border pixel), so if the problem is a minor offset in compute_deltas, its harder to see with interpolated. That said, it could also just bee that interpolated is correct.

[genusn]

Thank you for your comment.

Back to "ray-voxel", modifying

geo.nVoxel = np.array([2, 256, 256])   ## Sorry. I've fogotten this line at the first edition.
geo.nDetector = np.array([2, 256]) 

and printout the range of value of the projections

for idxY in range(geo.nDetector[0]):
    valMin = np.min(proj[:,idxY, :])
    valMax = np.max(proj[:,idxY, :])
    print(idxY, valMin, valMax)

the result was

0 0.0 80.25215
1 0.0 0.0

the values on v=1 are all zero.

[genusn]

Sorry, my comment above seems to be a different issue. In that condition, the result of

for idxZ in range(geo.nVoxel[0]):
    valMin = np.min(head[idxZ, :,:])
    valMax = np.max(head[idxZ, :,:])
    print(idxZ, valMin, valMax)

is

0 0.0 0.83519334
1 0.0 0.0

, which means head is zero in the plane of z=1 and the result shown in my last comment is the expected one.

[genusn]

I wanted to know the relative position of the center of the image and the x-ray source, or where the rotation axis (RA) is in the image, or where the RA is projected on. I prepared an image like a wire parallel to the RA, and bumped into this issue.

It seems to occur only when nVoxelZ and nDetectorV are both odd.

[AnderBiguri]

@genusn I thought all this could be caused by integer divisions, but I have not found such thing in the source.

That thing you show its quite strange. What are the angles? If you are making a vertical wire, its possible to get different intensities as we are dealing with voxels, and therefore ray-voxel intersections will have different lengths with rotation angle. However that "zero" is weird....

[genusn]

Here is the code to reproduce the above image and the result from other view:

import numpy as np
from matplotlib import pyplot as plt
import tigre

test_case = 3
if test_case == 0:
    # Both even
    nVoxelX=10
    nDetectorU=10
elif test_case == 1:
    # nVoxel Even, nDetec Odd
    nVoxelX=10
    nDetectorU=11
elif test_case == 2:
    # nVoxel Odd, nDetec Even
    nVoxelX=11
    nDetectorU=10
else:
    # Both Odd
    # Central plane all zero if projection_type=="Siddon"
    nVoxelX=11
    nDetectorU=11

# isotropic
nVoxelY=nVoxelX
nVoxelZ=nVoxelX # 2 # 3 # 5 #
nDetectorV=11#nDetectorU

## FOV is ~10mm cube.
## Detector size: ~20mm square
geo=tigre.geometry(mode='cone', default=True)
geo.DSD = 2000
geo.DSO = 1000
geo.nVoxel = np.array([nVoxelZ,nVoxelY,nVoxelX])
geo.dVoxel = np.array([1.0, 1.0, 1.0])
geo.sVoxel = geo.nVoxel * geo.dVoxel
geo.dDetector = np.array([2,2])               # size of each pixel            (mm)
geo.nDetector = np.array([nDetectorV, nDetectorU])
geo.sDetector = geo.dDetector * geo.nDetector
geo.offOrigin = np.array([0, 0, 0])
geo.offDetector = np.array([0.0, 0.0])
geo.rotDetector = np.array([0, 0, 0])
geo.accuracy=0.01

nangles=10
angles = np.linspace(0, 2 * np.pi, nangles, endpoint=False, dtype=np.float32)

wire = np.zeros(geo.nVoxel, dtype = np.float32)
wire[:, np.int32(geo.nVoxel[1]//2), np.int32(geo.nVoxel[2]//2)] = 1

# proj = tigre.Ax(wire,geo,angles, "interpolated")
proj = tigre.Ax(wire,geo,angles, "Siddon")
for iv in range(geo.nDetector[0]):
    plt.plot(proj[:, iv, geo.nDetector[1]//2], label='v={}'.format(iv))
plt.legend()
plt.savefig("plot_case{}.png".format(test_case))
plt.show()
tigre.plotProj(proj, savegif="proj_case{}.gif".format(test_case))
# fig, axes=plt.subplots(10,1)
# for it  in range(proj.shape[0]):
#     axes[it].imshow(proj[it])
# plt.show()

[genusn]

This also happens on MATLAB.

%% Init 
clear;
close all;
%% Test parameters
test_case = 3

if test_case == 0
    % Both even
    nVoxelZ=10;
    nDetectorV=10;
elseif test_case == 1
    % nVoxel Even, nDetec Odd
    nVoxelZ=10;
    nDetectorV=11;
elseif test_case == 2
    % nVoxel Odd, nDetec Even
    nVoxelZ=11;
    nDetectorV=10;
else
    % nVoxel Odd, nDetec Odd
    nVoxelZ=11;
    nDetectorV=11;
end

nVoxelX = nVoxelZ;
nVoxelY = nVoxelZ;
nDetectorU = nDetectorV;

nangles=10;

mode='cone'  % 'parallel'
ptype= 'Siddon' % 'interpolated' % 'ray-voxel'
%%
geo=defaultGeometry('mode',mode);
geo.DSD = 2000;
geo.DSO = 1000;
geo.nVoxel = [nVoxelX;nVoxelY;nVoxelZ];
geo.dVoxel = [1.0; 1.0; 1.0];
geo.sVoxel = geo.nVoxel .* geo.dVoxel;
geo.dDetector = [2;2] ;              % size of each pixel            (mm)
geo.nDetector = [nDetectorU; nDetectorV];     % has to be true 3D for python version
geo.sDetector = geo.dDetector .* geo.nDetector;
geo.offOrigin = [0;0;0];
geo.offDetector = [0.0; 0.0] ;
geo.rotDetector = [0; 0; 0];
geo.accuracy=0.01;

%%
angles = linspace(0, 2 * pi - 2*pi/nangles, nangles);

wire = zeros(geo.nVoxel', 'single');
wire(int32((geo.nVoxel(1))/2):int32((geo.nVoxel(1)+1)/2), int32((geo.nVoxel(2))/2):int32((geo.nVoxel(2)+1)/2), :) = 1;
%%
proj = Ax(wire,geo,angles, ptype);
%% Plot
cellLeg = {};
for iv =1:nDetectorV
    plot(squeeze(proj(iv, int32(nDetectorU/2), :)));
    cellLeg{iv} = sprintf('V=%d', iv);
    hold on
end
legend(cellLeg);
hold off
plotProj(proj, angles, 'savegif', sprintf('proj_%s_%s.gif', mode, ptype))

[AnderBiguri]

@genusn good to have confirmation, but expected. This has to do with the source code, I need to get to the bottom of why its happening. But unfortunately it will need to wait, as I don't have time now :(

[yliu88au]

What I found is that, instead of detector 0.5 offset on v, if we offset 0.5 on u direction, then the projection will work. This shows that, perhaps internally, python version has u and v exchanged. Hope this helps. @AnderBiguri @genusn

[genusn]

@yliu88au Thank you for the info. This happens both on Matlab and Python, not only 2d, but also 3d.

What I am looking at is a variable 'ray' in the kernel function. One or two components can be zero. If zero, the integral of the ray becomes zero(strange!) and if nonzero, such as 1e-8, the value is as expected. The "non-generic ray" does not seem to be treated proprely, so far to me. I'll look into further. What do you think @AnderBiguri @yliu88au ?

[AnderBiguri]

For uneven detector, I more or less understand what is happening here with the center pixel then, with your info.

The ray variable indeed suddenly starts having all zeroes due to the following code, as ray.z is definetly zero, while ray.x and ray.y can be zero.

https://github.com/CERN/TIGRE/blob/09030ea1a79b8d62250687e7d8ac4cef088abe9b/MATLAB/Source/Siddon_projection.cu#L160-L165

Maybe we need to introduce some other line like this:

https://github.com/CERN/TIGRE/blob/09030ea1a79b8d62250687e7d8ac4cef088abe9b/MATLAB/Source/Siddon_projection.cu#L672

But I am thinking that the best solution should be to add an epsilon to those divisions. A bit worried about their impact on the projection, I guess it needs to be tested.

That said, not sure if I understand what is going on with the other pixels in the image...maybe its just too few angles, running the same with 360 gives a smooth transition of intensity there? Should me something around [1, sqrt(2)]*pval

[AnderBiguri]

@genusn disregard the previous comment, eps in there does not seem to fix the issue.... Confused.

[AnderBiguri]

Here is the issue in the kernels. The following epsilons do not help us get rid of the division by zero, they only make the problem worse!

https://github.com/CERN/TIGRE/blob/09030ea1a79b8d62250687e7d8ac4cef088abe9b/MATLAB/Source/Siddon_projection.cu#L204-L206

In the horizontal line on uneven detectors, ray.z==0, which means that there is a division by zero. This means that there will be zero intersections in the z direction in the image, which makes sense.

The epsilon helps us avoid a numerical issue, but introduces another one. az becomes something like -5000000000000, thus aminc also does, the initial intersection with the voxel array.

https://github.com/CERN/TIGRE/blob/09030ea1a79b8d62250687e7d8ac4cef088abe9b/MATLAB/Source/Siddon_projection.cu#L212

The in the loop, we will start then from this value to read image values (L247)

https://github.com/CERN/TIGRE/blob/09030ea1a79b8d62250687e7d8ac4cef088abe9b/MATLAB/Source/Siddon_projection.cu#L236-L254

azu is correct, however, the "unit" change between intersections in z. Its Inf.

Thus, after the first image read with this wrong az value, the next intersection becomes Inf, so none of the if conditions apply anymore.

https://github.com/CERN/TIGRE/blob/09030ea1a79b8d62250687e7d8ac4cef088abe9b/MATLAB/Source/Siddon_projection.cu#L253

As the other parameters are correct ( most importantly, Np number of intersections), this iterates Np times without accessing any of the conditions, thus returns zero.

This was an example for that horizontal line we see in this demo, but obviously the same happens with x and y.

---

I need to carefully think how these special "trivial" lines are/should be handled by the kernel. One could assume that there is something there to do with azu and az, but adding an epsilon may not work. Maybe we need to catch Infs, but these things tend to be computationally costly..... hum....

[AnderBiguri]

Maybe a temporary solution is to test properly interpolated projector for these cases and enforce it, particularly in 2D.

[genusn]

@AnderBiguri Another solution is to introduce If-trivial-clause. If the trivial cases are rare, branch divergence may not be a problem.... I am not sure yet.

Or, add tiny offset of the angle (to avoid vertical blackout) and V (to avoid horizontal blackout)?

[AnderBiguri]

@genusn yes, I will test the trivial case solution, maybe worth considering.

I will first see if I can find a solution without touching the geo (even if we already to it quite often...), but that also can be a solution. Looking into this now.

[AnderBiguri]

I think I have solved the axis-parallel ray issue in #218

[genusn]

Thank you very much. Could you please apply 24d7101 to python?

[AnderBiguri]

@genusn forgot to remove the eps, but otherwise should be there. Give me a sec.

[yliu88au]

PR218 It seems fixed the 2D problem as well as 0.5 offset. Great!

On the other hand, algorithms in Python are still buggy, particularly for those "POCS" family. I am trying to fix a few...

[AnderBiguri]

@yliu88au ah that's great! One bug at a time. Can you open an issue with what you see? Which algorithms aside from Pocs you see go bad?