CERN / TIGRE

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

Fix python BP type variable name and default #305

Closed tsadakane closed 3 years ago

tsadakane commented 3 years ago

Summary

TODO

tsadakane commented 3 years ago

TEST

Code

  1. Run the following test code on the base workspace.
  2. Copy three npy files to the PR workspace.
  3. Modify the test code (False to True) to use new variable names
  4. Run the modified test code using the PR code.
    
    import numpy as np
    import tigre
    from tigre.utilities import sample_loader
    import tigre.utilities.gpu as gpu
    # from tigre.utilities import CTnoise

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)

exit()

Geometry

geo = tigre.geometry(mode="cone", nVoxel=np.array([256, 256, 256]), default=True) geo.dDetector = np.array([0.8, 0.8]) 2 # size of each pixel (mm) geo.sDetector = geo.dDetector geo.nDetector

print(geo)

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

Prepare projection data

head = sample_loader.load_head_phantom(geo.nVoxel) proj = tigre.Ax(head, geo, angles, gpuids=gpuids)

noise_proj = CTnoise.add(proj, Poisson=1e5, Gaussian=np.array([0, 10]))

if False: # SWITCH PR and BASE

PR

backproj_default = tigre.Atb(proj, geo, angles, gpuids=gpuids)
backproj_FDK     = tigre.Atb(proj, geo, angles, gpuids=gpuids, backprojection_type="FDK")
backproj_matched = tigre.Atb(proj, geo, angles, gpuids=gpuids, backprojection_type="matched")
backproj_default_base = np.load("arr_backproj_default.npy")
backproj_FDK_base     = np.load("arr_backproj_FDK.npy"    )    
backproj_matched_base = np.load("arr_backproj_matched.npy")

diff_base    = np.abs(backproj_default - backproj_default_base)
diff_FDK     = np.abs(backproj_FDK - backproj_FDK_base)
diff_matched = np.abs(backproj_matched - backproj_matched_base)
diff_FDK_default = np.abs(backproj_FDK - backproj_default)
print("max diff_base    = {}".format(np.max(diff_base)))
print("max diff_FDK     = {}".format(np.max(diff_FDK)))
print("max diff_matched = {}".format(np.max(diff_matched)))
print("Check if new default is FDK")
print("max backproj_FDK-backproj_default = {}".format(np.max(diff_FDK_default)))

else:

base

backproj_default = tigre.Atb(proj, geo, angles, gpuids=gpuids)
backproj_FDK     = tigre.Atb(proj, geo, angles, gpuids=gpuids, krylov="FDK")
backproj_matched = tigre.Atb(proj, geo, angles, gpuids=gpuids, krylov="matched")
np.save("arr_backproj_default.npy", backproj_default)
np.save("arr_backproj_FDK.npy"    , backproj_FDK)
np.save("arr_backproj_matched.npy", backproj_matched)
## Result

max diff_base = 208.908935546875 max diff_FDK = 0.0 max diff_matched = 0.0 Check if new default is FDK max backproj_FDK-backproj_default = 0.0


## Environment
* Windows 10
*  Python 3.8.5 (anaconda)
*  Visual studio 2019 (16.9.6)
AnderBiguri commented 3 years ago

Thanks! Will test and merge :+1:

AnderBiguri commented 3 years ago

The CGLS algorithm (https://github.com/CERN/TIGRE/blob/master/Python/tigre/algorithms/krylov_subspace_algorithms.py) uses not FDK, which is wrong. The 3 Atb() calls in this file need to be matched

AnderBiguri commented 3 years ago

This is the same for statistical algorithms, I believe.

tsadakane commented 3 years ago

@AnderBiguri Thank you for your review. I didn't see around there.

This is the same for statistical algorithms, I believe.

MLEM of MATLAB also uses 'FDK'. Shall I change both or leave them alone?

AnderBiguri commented 3 years ago

@tsadakane I think lets change both. I need to do some numerical tests about this, as it worked well with FDK, but I have been recently advised in a different project that MLEM/OSEM need adjoint backprojectors. Do change them I will investigate the inpact on image quality.