CERN/TIGRE

Sinogram/Projections are zero for angle zero

roflmaostc opened this issue · 11 comments

Thanks for this package!

For the angles equal to 0, 2pi, 4pi the sinogram seems to be zero.

Expected Behavior

Should be non zero

Actual Behavior

Did I assume the orientation correctly? [Z, Y, X]?

For angle 0, the sinogram should show the same shape. I shifted them to make it more clearly.
Figure_1

Code to reproduce the problem (If applicable)

import numpy as np
import tigre
import tigre.algorithms as algs
from tigre.utilities import sample_loader
from tigre.utilities.Measure_Quality import Measure_Quality
import tigre.utilities.gpu as gpu
import matplotlib.pyplot as plt

listGpuNames = gpu.getGpuNames()
if len(listGpuNames) == 0:
    print("Error: No gpu found")
else:
    for id in range(len(listGpuNames)):
        print("{}: {}".format(id, listGpuNames[id]))

gpuids = gpu.getGpuIds(listGpuNames[0])
print(gpuids)


sz = (1, 200, 200)
geo = tigre.geometry(mode="parallel", nVoxel=np.array(sz), default=True)

sample = np.zeros(sz, dtype=np.float32)
x = np.linspace(-sz[1] // 2, sz[1] // 2, sz[1])
X, Y = np.meshgrid(x + 40,x)#
sample[0, :, :] += ((X**2 + Y**2) <= 10**2)


angles = np.linspace(0, np.pi, 4)
proj = tigre.Ax(sample, geo, angles, gpuids=gpuids)


plt.figure()
for i in range(4):
    plt.plot(i + proj[i, 0, :], label=angles[i])
    plt.legend()

plt.show()

Specifications

  • MATLAB/python version: 3.10
  • OS: Ubuntu
  • CUDA version: 12.1

Hum, worrying indeed.
I will try to fix ASAP. Can you do a small test for me in the meantime?

Can you make the detector double in number of pixels nDetector and half in size of pixels dDetector and try again?
Also, can you try to give a small epsilon to the angles (e.g. 0.00000000001)?

With an angle offset of 2e-6 it still occurs (2e-6 is still 40 times larger than single float machine epsilon)

import numpy as np
import tigre
import tigre.algorithms as algs
from tigre.utilities import sample_loader
from tigre.utilities.Measure_Quality import Measure_Quality
import tigre.utilities.gpu as gpu
import matplotlib.pyplot as plt

listGpuNames = gpu.getGpuNames()
if len(listGpuNames) == 0:
    print("Error: No gpu found")
else:
    for id in range(len(listGpuNames)):
        print("{}: {}".format(id, listGpuNames[id]))

gpuids = gpu.getGpuIds(listGpuNames[0])
print(gpuids)


sz = (1, 200, 200)
geo = tigre.geometry(mode="parallel", nVoxel=np.array(sz), default=True)

sample = np.zeros(sz, dtype=np.float32)
x = np.linspace(-sz[1] // 2, sz[1] // 2, sz[1])
X, Y = np.meshgrid(x + 40,x)#
sample[0, :, :] += ((X**2 + Y**2) <= 10**2)


angles = np.linspace(0, np.pi, 4) + 0.000002 # 2e-6
proj = tigre.Ax(sample, geo, angles, gpuids=gpuids)


plt.figure()
for i in range(4):
    plt.plot(i + proj[i, 0, :], label=angles[i])
    plt.legend()

plt.show()

Not sure how I change the nDetector since the geometry class throws an error.

     21 geox = copy.deepcopy(geo)
---> 22 geox.check_geo(angles)
     23 """
     24 Here we cast all values in geo to single point precision float. This way we know what behaviour
     25 to expect from pytigre to Cuda and can change single parameters accordingly.
     26 """
     27 geox.cast_to_single()

AttributeError: 'Geometry' object has no attribute 'check_geo'
``

Interesting. I will investigate.

I think the error is because you changed nDetector, but not dDetector (for the same sDetector, both need to change, they are redundant variables, sDetector=nDetector*dDetector), but don't worry, I think I need to dig this myself.

sounds more like a function (check_geo) is missing in the template?

Yes, you are right, I didn't read the message well.
Huh, another weird thing, the class Goemetry should have the method check_geo() https://github.com/CERN/TIGRE/blob/master/Python/tigre/utilities/geometry.py#L25

I was copying the geometry template from the web page.

Maybe there's something missing?

@roflmaostc ah, on second read, that webpage is seriously outdated.... I suggest looking at the demos in the demos folder, those should be up to date and running.

1e-5 works, so it definetly is a numerical precission stuff.

But I need to dig the actual core to see what causes it....

@roflmaostc Fixed! :)

Thanks for quick fix!