CERN / TIGRE

TIGRE: Tomographic Iterative GPU-based Reconstruction Toolbox
BSD 3-Clause "New" or "Revised" License
529 stars 180 forks source link

Sinogram/Projections are zero for angle zero #453

Closed roflmaostc closed 12 months ago

roflmaostc commented 1 year ago

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

AnderBiguri commented 1 year ago

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)?

roflmaostc commented 1 year ago

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'
``
AnderBiguri commented 1 year ago

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.

roflmaostc commented 1 year ago

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

AnderBiguri commented 1 year ago

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

roflmaostc commented 1 year ago

I was copying the geometry template from the web page.

Maybe there's something missing?

AnderBiguri commented 1 year ago

@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.

AnderBiguri commented 1 year ago

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

AnderBiguri commented 1 year ago

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

AnderBiguri commented 12 months ago

@roflmaostc Fixed! :)

roflmaostc commented 12 months ago

Thanks for quick fix!