llnl/leap

comprehensive library of 3D transmission Computed Tomography (CT) algorithms with Python and C++ APIs and fully integrated with PyTorch

LivermorE AI Projector for Computed Tomography (LEAP)

This is a C++/CUDA library (Linux, Windows, and Mac*) of 3D tomographic algorithms (pre-processing algorithms, projectors, and analytic (FBP) and iterative reconstruction algorithms) with a Python interface. The projectors (forward and back projection) are implemented for both multi-GPU and multi-core CPU and we provide bindings to PyTorch to achieve differentiable forward and backward projectors for AI/ML-driven Computed Tomography (CT) applications.

There are a lot of CT reconstruction packages out there, so why choose LEAP? In short, LEAP has more accurate projectors and FBP algorithms, more features, and most algorithms run as fast or faster than other popular CT reconstruction packages, but here is a more detailed list:

1. Seamless integration with PyTorch using torch.nn.Module and torch.autograd.Function to enable differentiable forward and backward projectors for AI/ML-driven Computed Tomography (CT) applications.
2. Quantitatively accurate, matched (forward and back) projector pairs that model the finite size of the voxel and detector pixel; very similar to the Separable Footprint method [Long, Fessler, and Balter, TMI, 2010]. These matched projectors ensure convergence and provide accurate, smooth results. Unmatch projectors or those projectors that do not model the finite size of the voxel or detector pixel may produce artifacts when used over enough iterations [DeMan and Basu, PMB, 2004].
3. Multi-GPU and multi-core CPU implementations of all algorithms that are as fast or faster than other popular CT reconstruction packages.
4. Algorithms not limited by the amount of GPU memory.
5. Flexible 3D CT geometry specification that allows users to specify arbitrary shifts of the source and detector positions, non-uniform angular spacing, and more.
6. Flexible 3D CT volume specification.
7. Quantitatively accurate and flexible analytic reconstruction algorithms, i.e., Filtered Backprojection (FBP).
8. Can avoid costly CPU-to-GPU data transfers by performing operations on data already on a GPU.
9. Special-case FBP algorithms that are rarely included in other packages, such as helical, truncated projections, offset detector scan, and Attenuated Radon Transform.
10. Special-case models such as the Attenuated Radon Transform (SPECT and VAM applications) and reconstruction of cylindrically-symmetric objects (flash x-ray applications).
11. Iterative reconstruction algorithms: OSEM, OS-SART, ASD-POCS, RWLS, RDLS, ML-TR, IFBP (RWLS-SARR)
12. Fast multi-GPU 3D densoing methods.
13. Pre-processing algorithms: outlier correction, detector deblur, ring removal, scatter correction, metal artifact reduction (MAR), multi-material beam hardening correction (BHC), dual energy decomposition, and SIRZ
14. Easy-to-use, simple API.
15. Easy-to-build executable because the only dependency is CUDA. Python API can be run with or without PyTorch (of course the neural network stuff requires PyTorch).
16. Permissible license.

Physics-based modeling and correction algorithms (e.g., scatter correction, beam hardening correction (BHC), dual energy decomposition, and SIRZ) can be applied when used with the XrayPhysics package.

*Mac version does not have GPU support and some featurings are missing.

Installation and Usage

Documentation is available here

Installation and usage information is posted on the wiki page

Demo scripts for most functionality in the demo_leapctype directory

Demo scripts for AI/ML/DL applications in the demo_leaptorch directory

Example Results

As a simple demonstration of the accuracy of our projectors we show below the results of FDK reconstructions using ASTRA and LEAP of the walnut CT data. The LEAP reconstruction has 1.7 times higher SNR and reconstructed this data 7.5 times faster than ASTRA.

Future Releases

For the next releases, we are working on the following:

1. Fixes of bugs reported by our users
2. Feature requests from our users
3. More noise reduction filters
4. AMD GPU Support
5. cone-parallel geometry support
6. PyQt GUI
7. beam hardening correction algorithms that account for variable takeoff angle and graded collimator/ bowtie filter

Authors

Kyle Champley ([email protected])

Hyojin Kim ([email protected])

License

LEAP is distributed under the terms of the MIT license. All new contributions must be made under this license. See LICENSE in this directory for the terms of the license. See LICENSE for more details.
SPDX-License-Identifier: MIT
LLNL-CODE-848657

Please cite our work by referencing this github page and citing our article:

Hyojin Kim and Kyle Champley, "Differentiable Forward Projector for X-ray Computed Tomography”, ICML, 2023

leap's Issues

errors when running demo

I have installed the leapct, but the errors encounter when I run the demo as below:

Can you help me with this issue?

encountered an error while running the demo

Merry Christmas !
After installation, I encountered an error while running the code /LEAP/demo_leaptorch/test_fproject_and_FBP.py. How should I handle this?

Traceback (most recent call last):
File "/data4/liqiaoxin/.pycharm_helpers/pydev/pydevconsole.py", line 364, in runcode
coro = func()
File "", line 1, in
File "/data4/liqiaoxin/.pycharm_helpers/pydev/_pydev_bundle/pydev_umd.py", line 198, in runfile
pydev_imports.execfile(filename, global_vars, local_vars) # execute the script
File "/data4/liqiaoxin/.pycharm_helpers/pydev/_pydev_imps/_pydev_execfile.py", line 18, in execfile
exec(compile(contents+"\n", file, 'exec'), glob, loc)
File "/data4/liqiaoxin/code/LEAP/demo_leaptorch/test_fproject_and_FBP.py", line 25, in
from leaptorch import Projector
File "/home/liqiaoxin/anaconda3/envs/pytorch/lib/python3.9/site-packages/leaptorch.py", line 13, in
lct = tomographicModels()
File "/home/liqiaoxin/anaconda3/envs/pytorch/lib/python3.9/site-packages/leapctype.py", line 78, in init
self.libprojectors = cdll.LoadLibrary(os.path.join(current_dir, "../build/lib/libleap.so"))
File "/home/liqiaoxin/anaconda3/envs/pytorch/lib/python3.9/ctypes/init.py", line 460, in LoadLibrary
return self._dlltype(name)
File "/home/liqiaoxin/anaconda3/envs/pytorch/lib/python3.9/ctypes/init.py", line 382, in init
self._handle = _dlopen(self._name, mode)
OSError: /home/liqiaoxin/anaconda3/envs/pytorch/lib/python3.9/site-packages/../build/lib/libleap.so: cannot open shared object file: No such file or directory

powell method optimizer function implementation by gpu code

@kylechampley Do you have any plan to implement current python's scipy.optimize.minimize into your own gpu c++ code ?
That could bring faster auto geo-calibration result.

Discrepancy (in Cone-beam geometry) between applying filterProjections then backproject and directly using FBP

Hi,

First of all, thank you for the amazing work!

I'm encountering some discrepancy issue for conebeam geometry. Basically, I tried to filter the projection first and then backproject the filtered projection. However, the results are different from directly calling FBP function (the latter seems to give me correct attenuation coefficient value in reconstruction, but the former was about ~0.1x the correct value). I tried to find what FBP does uniquely, and thought that maybe weightbackprojection was the cause, but even applying weightedbackprojection to filtered projection, the issue was still there. I would really appreciate some help on this.

Here is some part of the code that I use for filter and then backproject:

  projector = leaptorch.Projector(forward_project=False,use_static=False,use_gpu=True,gpu_device=torch.device(device),batch_size=1)
  projector.set_conebeam(numAngles, 
                              numRows, 
                              numCols, 
                              pixelHeight, 
                              pixelWidth, 
                              centerRow, 
                              centerCol, 
                              phis, 
                              sod, 
                              sdd)
  projector.set_volume(numX, numY, numZ, voxelWidth, voxelHeight)
  projector.allocate_batch_data()
  projF = projector.leapct.filterProjections(proj.clone())
  vol = projector(projF)

when trying to use the weightedbackprojector, I replaced the last line of code with this
vol = projector.leapct.weightedBackproject(projF[0], projector.vol_data[0])
But it did not help

I want specific descriptions about some sub-functions on the tomographicModels

Hi!

I'm very thrilled using the LEAP CT tool, adopting things i used on MATLAB.

Especially, the FDK runtime reconstructing my own real-data from C-arm felt brilliant.
Such fast and correct reconstruction was applied.

However, I have two problems using this tool on my research project.

First, for using the function " leapct.load_param(param_fn) " I can't guess what form of data would be appropriate for the "param_fn".
I want easily allocating the CT geometry parameters by using this function.

Second, I an really looking ahead for the projection tool on LEAP.

The example code shows that syntax like " f = leapct.load_volume(img_fn) " can help loading volumes for forward projection

but i also can't find the data form of img_fn . I tried really hard searching for guides, and leaked for the source code, but i have no idea about it. Please give me a guide and an example of the img_fn data.

Thank you.

CUDA related error

When I use ASD-POCS using my data. Some cuda-related error happen.
The forward projection and FBP works. I'm using LEAP v1.13
In Linux the error is

  File "/home/exouser/stefen/code/DDS/physics/multileap4GPU.py", line 44, in IR
    leapct.ASDPOCS(proj,rec,iter_num,10,10,filters)
  File "/home/exouser/.local/lib/python3.10/site-packages/leapctype.py", line 2971, in ASDPOCS
    P1[P1==0.0] = 1.0
RuntimeError: CUDA error: an illegal memory access was encountered
CUDA kernel errors might be asynchronously reported at some other API call, so the stacktrace below might be incorrect.
For debugging consider passing CUDA_LAUNCH_BLOCKING=1.
Compile with `TORCH_USE_CUDA_DSA` to enable device-side assertions.

In Windows, the error is

  File "D:\code\improved diffusion\DDS\physics\multileap4GPU.py", line 44, in IR
    leapct.ASDPOCS(proj,rec,iter_num,10,10,filters)
  File "C:\Users\Shuo_Han\AppData\Local\anaconda3\envs\diff\lib\site-packages\leapctype.py", line 2987, in ASDPOCS
    Pf_minus_g = self.allocateData(g)
  File "C:\Users\Shuo_Han\AppData\Local\anaconda3\envs\diff\lib\site-packages\leapctype.py", line 1429, in allocateData
    y = torch.zeros([x.shape[0], x.shape[1], x.shape[2]], dtype=torch.float32, device=torch.device('cuda:'+str(self.get_gpu())))
RuntimeError: CUDA error: an illegal memory access was encountered
CUDA kernel errors might be asynchronously reported at some other API call,so the stacktrace below might be incorrect.
For debugging consider passing CUDA_LAUNCH_BLOCKING=1.

Calculate the HU value

The Hounsfield scale of CT is set around water measuring 0 HU. By altering the attenuation level (mid-HU value) and range (extent of gray scale) at which the image is viewed, the tissue to be displayed can be determined.
Will LEAP add a function to calculate the HU value automatically?

Support AAPM low-dose CT helical reconstruction

Although Low Dose CT Grand Challenge is 2016 dataset. Nowadays still many paper use it. But it seems not a clear code to guide how to reconstruct to CT image from the raw projection data.
If LEAP can do that, I can be use to the down stream task(sparse-view & limited angle).

curved detector

Hey! Does LEAP handle curved detectors? I am working with medical data and they all have curved detectors so I need software for this.

Oblique CT supporting

Is there any demo example of oblique CBCT(Laminography CT) at LEAP packet?

Angular range problem

Thanks for your contribution! In the following code does 360 means 2pi?
leapct.set_conebeam(numAngles, numRows, numCols, pixelSize, pixelSize, 0.5*(numRows-1), 0.5*(numCols-1), leapct.setAngleArray(numAngles, 360.0), 1100, 1400)
I'm wondering what if I want to use the sparse-view or limited-angle setting how do I change the code. And how to adjust the start_angle, cause in TIGRE it's easy to do that as following.
start_angle = np.deg2rad(140)
end_angle = start_angle + 2 * np.pi
angles = np.linspace(start_angle, end_angle, 984, dtype=np.float32)
And for the iterative algorithm, I saw you have SART reconstruction and TV denoising, is it easy to combine them together to do the sparse-view reconstruction? Like usually when we use SART we will add one regularization term.

Compile error No cuda toolset found

Hi, thanks for your work. Based on the instruction I already installed cmake 3.29 and visual studio 2019, CUDA version is as follows, but I ran into the compile issue, could you please instruct me how to fix this? And how to use the precompile file? Just put the two file into cmake file folder? libleap.dll
libleap.so
(leap) PS D:\ssh_down\test_code\recons_code\LEAP> nvcc -V
nvcc: NVIDIA (R) Cuda compiler driver
Copyright (c) 2005-2023 NVIDIA Corporation
Built on Tue_Aug_15_22:09:35_Pacific_Daylight_Time_2023
Cuda compilation tools, release 12.2, V12

(leap) PS D:\ssh_down\test_code\recons_code\LEAP> pip install -v .
Using pip 23.3.1 from C:\Users\Shuo_Han\AppData\Local\anaconda3\envs\leap\lib\site-packages\pip (python 3.10)
Processing d:\ssh_down\test_code\recons_code\leap
Running command python setup.py egg_info

D:\ssh_down\test_code\recons_code\LEAP>rd /s /q win_build

D:\ssh_down\test_code\recons_code\LEAP>mkdir win_build

D:\ssh_down\test_code\recons_code\LEAP>cd win_build

D:\ssh_down\test_code\recons_code\LEAP\win_build>cmake .. -G "Visual Studio 16" -DDEV_MODE=1
CMake Deprecation Warning at CMakeLists.txt:1 (cmake_minimum_required):
Compatibility with CMake < 3.5 will be removed from a future version of
CMake.

Update the VERSION argument <min> value or use a ...<max> suffix to tell
CMake that the project does not need compatibility with older versions.

-- The C compiler identification is MSVC 19.29.30154.0
-- The CXX compiler identification is MSVC 19.29.30154.0
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working C compiler: C:/Program Files (x86)/Microsoft Visual Studio/2019/Enterprise/VC/Tools/MSVC/14.29.30133/bin/Hostx64/x64/cl.exe - skipped
-- Detecting C compile features
-- Detecting C compile features - done
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: C:/Program Files (x86)/Microsoft Visual Studio/2019/Enterprise/VC/Tools/MSVC/14.29.30133/bin/Hostx64/x64/cl.exe - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
CMake Error at D:/code/ct recons/cmake-3.29.0-windows-x86_64/share/cmake-3.29/Modules/CMakeDetermineCompilerId.cmake:563 (message):
No CUDA toolset found.
Call Stack (most recent call first):
D:/code/ct recons/cmake-3.29.0-windows-x86_64/share/cmake-3.29/Modules/CMakeDetermineCompilerId.cmake:8 (CMAKE_DETERMINE_COMPILER_ID_BUILD)
D:/code/ct recons/cmake-3.29.0-windows-x86_64/share/cmake-3.29/Modules/CMakeDetermineCompilerId.cmake:53 (__determine_compiler_id_test)
D:/code/ct recons/cmake-3.29.0-windows-x86_64/share/cmake-3.29/Modules/CMakeDetermineCUDACompiler.cmake:131 (CMAKE_DETERMINE_COMPILER_ID)
src/CMakeLists.txt:9 (project)

-- Configuring incomplete, errors occurred!

D:\ssh_down\test_code\recons_code\LEAP\win_build>cmake --build . --config Release
Microsoft (R) Build Engine version 16.11.2+f32259642 for .NET Framework
Copyright (C) Microsoft Corporation. All rights reserved.

MSBUILD : error MSB1009: Project file does not exist.
Switch: ALL_BUILD.vcxproj
Failed to compile!
Preparing metadata (setup.py) ... done
ERROR: No .egg-info directory found in C:\Users\Shuo_Han\AppData\Local\Temp\pip-pip-egg-info-8d_p5qgr

cbct post processing

What post-processing can be used to convert the right image into the left image？thanks

bug in axially-aligned modular-beam forward projector

There is a bug in axially-aligned modular-beam forward projector causing incorrect values. It gets worse as the detector rotation around the optical axis increases. If there is no rotation, everything is fine. The backprojector is not affected.

I've identified the bug and am working on resolving it. I intend on fixing this and releasing a new version in the next few days.

Total Variation functions do not work for "large" arrays

The Total Variation (TV) functions (cost, gradient quadratic form) do not work for volumes larger than about 1K^3. We are currently working on implementing these algorithms so that they split up the volume into smaller chunks that can be processed on a GPU. We will also make these function compute using multiple GPUs while we are at it.

erros when runing demos

I have successfully installed the leapct "pip install ." as below.

But when I run the demo "test_recon_NN.py", the error enconters:

Can you help me with this issue?

remove cone-beam artifacts

For 128-slice(or more) ct , cone-beam artifacts are severe,Whether LEAP supports removal of cone-beam artifacts?

Helical FBP does not work properly

The helical FBP reconstruction does not work properly and we are currently working on fixing this.

Volume coordination

I have a question regarding the orientation of the resulting volume. I believe the row direction of the volume corresponds to the left-right (LR) direction, and the column direction corresponds to the anterior-posterior (AP) direction. However, as shown in the image above, the row and column directions seem to be swapped in the resulting volume. Could you clarify the correct orientation of the rows and columns in the output?

Build error

When I try to install with "pip install ." it fails to build the wheel returning the following error. There seem to be several cases named parallelBeam rather than parallelbeam.

src/main_projector.cpp:606:15: error: ‘class tomographicModels’ as no member named ‘set_parallelBeam’; did you mean ‘set_parallelbeam’?
606 | tempModel.set_parallelBeam(numAngles, numRows, numCols, pixelHeight, pixelWidth, centerRow, centerCol, phis);
| ^~~~~~~~~~~~~~~~
| set_parallelbeam
error: command '/usr/bin/gcc' failed with exit code 1

any advice on projection angle correction

Thank you for the recent update on geometric correction, which has been extremely useful for solving real-world problems. In my CBCT system, I am facing an issue with an unreliable angle sensor, so I am forced to estimate the angles from 0 to 360 degrees, which unfortunately degrades the quality of the reconstructions. I have tried to maximize the gradient/sharpness of volume, but such loss may maximize the noise. I believe that adjusting the 'sweep_parameter' and 'consistency_cost' might help address this issue. I would appreciate any advice you can offer on this matter. Thank you once again.

geometric calibration can find centerRow ?

offsetScan not working?

issue: d12_geometric_calibration.py

The new geometric_calibration.py in leap v1.13 is very useful. However, it would be even more helpful if it could also work with offset scans.

Additionally, I have two issues:
First issue:
Is leapct.convert_to_modular() correct at line 214 of d12_geometric_calibration.py? I suspect it should be convert_to_modularbeam().

Second issue:
In the example of d12_geometric_calibration.py, it seems to find up to three variables correctly. However, it doesn't seem to find four geometric variables (centerRow, centerCol, tau, and detector tilt). I tested it by modifying the code as follows. Is the code I changed incorrect, or is it that optimization doesn't work well with four variables?

    # First we perturb the true geometry specification
    Delta_centerRow = 10.0
    Delta_centerCol = 20.0
    Delta_tau = 3.0

    leapct.set_centerRow(leapct.get_centerRow() + Delta_centerRow)
    leapct.set_centerCol(leapct.get_centerCol() + Delta_centerCol)
    leapct.set_tau(Delta_tau)

    # Specify the number of times the minimize function is called
    # This may help estimate a more accurate estimate
    numIter = 2

    from scipy.optimize import minimize

    # Start by trying to find the center column
    leapct.find_centerCol(g)

    numIter = 10
    costFcn_rct = lambda x: leapct.consistency_cost(g, x[0], x[1], x[2], x[3])
    x0 = (0.0, 0.0, 0.0, 0.0)
    bounds = ((-50, 50), (-50, 50), (-5, 5), (-2, 2))
    for n in range(numIter):
        res = minimize(costFcn_rct, x0, method='Powell', bounds=bounds)
        x0 = res.x

    # Print out estimated parameters and various cost values
    print(res.message)
    print('estimated perturbations: ' + str(res.x))
    print('optimized cost = ' + str(costFcn_rct(res.x)))

    # A significant detector rotation was detected, so let's assume this is correct
    # In this case, we have to switch to modular-beam geometry
    leapct.set_centerRow(leapct.get_centerRow() + res.x[0])
    leapct.set_centerCol(leapct.get_centerCol() + res.x[1])
    leapct.set_tau(leapct.get_tau() + res.x[2])
    leapct.convert_to_modularbeam()
    leapct.rotate_detector(res.x[3])

    print('Current estimate of detector center: ' + str(leapct.get_centerRow()) + ', ' + str(leapct.get_centerCol()))

    print('Elapsed time: ' + str(time.time() - startTime) + ' seconds')

Can't use libleapct.dll wheninstalled manually on Windows

When installed manually with CUDA Toolkit 12.4, there is a dependency missing.

FileNotFoundError: Could not find module 'C:\Users\user\.conda\envs\LEAP\Lib\site-packages\libleapct.dll' (or one of its dependencies). Try using the full path with constructor syntax.

The library libleapct.dll needs cufft64_10.dll, but CUDA Toolkit's file is named cufft64_11.dll.
The issues can be resolved by copying cufft64_11.dll in cuda toolkit installation folder and renaming the copy to cufft64_10.dll.

Different results on CPU and GPU

Hello! I have tested the code and observed that I get different results when forward projecting on the CPU vs the GPU. See the image below. Is this expected behavior?

Code to reproduce:

import torch
from leaptorch import Projector

def main():
    device = torch.device("cuda:0")

    proj = Projector(use_gpu=True, gpu_device=device)

    dimx = 128
    dimy = 128
    dimz = 64
    offsetx = 0
    offsety = 0
    offsetz = 0
    nangles = 180
    nrows = dimz
    ncols = max(dimx, dimy)
    sdd = 1350
    sod = 930
    pheight = sdd/sod
    pwidth = sdd/sod
    width = 1 
    height = 1 
    arange = 180
    crow = 0.5*float(nrows - 1)
    ccol = 0.5*float(ncols - 1)
    phis = torch.arange(0, arange, arange//nangles).float()

    proj.set_volume(dimx=dimx,
                    dimy=dimy,
                    dimz=dimz,
                    width=width,
                    height=height,
                    offsetx=offsetx,
                    offsety=offsety,
                    offsetz=offsetz)

    proj.set_cone_beam(nangles=nangles,
                    nrows=nrows,
                    ncols=ncols,
                    pheight=pheight,
                    pwidth=pwidth,
                    crow=crow,
                    ccol=ccol,
                    arange=arange,
                    phis=phis,
                    sod=sod,
                    sdd=sdd)

    volume = torch.ones((1,64,128,128)).float().to(device)
    volume[..., 40:80, 40:80] = 3.0  # Create a cuboid

    sinogram_gpu = proj(volume)

    device = torch.device("cpu")
    proj_cpu = Projector(use_gpu=False)
    proj_cpu.set_volume(dimx=dimx,
                    dimy=dimy,
                    dimz=dimz,
                    width=width,
                    height=height,
                    offsetx=offsetx,
                    offsety=offsety,
                    offsetz=offsetz)

    proj_cpu.set_cone_beam(nangles=nangles,
                    nrows=nrows,
                    ncols=ncols,
                    pheight=pheight,
                    pwidth=pwidth,
                    crow=crow,
                    ccol=ccol,
                    arange=arange,
                    phis=phis,
                    sod=sod,
                    sdd=sdd)

    sinogram_cpu = proj_cpu(volume.to(device))
    import matplotlib.pyplot as plt

    fig, axs = plt.subplots(1,2)
    axs[0].imshow(sinogram_gpu[0,:,32,:].detach().cpu().numpy())
    axs[0].set_title("GPU")
    axs[1].imshow(sinogram_cpu[0,:,32,:].numpy())
    axs[1].set_title("CPU")
    fig.show()
    fig.savefig("test_diff.jpg")
if __name__ == "__main__":
    main()

How to define an arbitrarily oriented geometry?

import leaptorch error

leapctype.py", line 149, in create_new_model
self.libprojectors.create_new_model.restype = ctypes.c_int
AttributeError: 'NoneType' object has no attribute 'create_new_model'

Help me: Offsetscan

Thank you for sharing the excellent toolkit.

I am currently interested in offset scan and am trying it with the data I have. However, I seem to have made a mistake, and the reconstruction isn’t working well.

Could you possibly help me if I share the geometry information and data with you?

Thank you very much.

Pytorch optimization through FBP

Hi:

Firstly, thank you very much for the amazing tool.

I have a question for the optimization through the FBP algorithm. For example, if I use the forward projection, the optimization is working (no error).

Set up the projector

proj = Projector(forward_project=True, use_static=True, use_gpu=True, gpu_device=0, batch_size=1)
proj.set_volume(512, 512, 4, 1, 1, 0, 0, 0)
proj.set_conebeam(360, 8, 512, 1, 1, 256, 256, proj.leapct.setAngleArray(360, 360.0), 1000, 1500, tau=0.0, helicalPitch=0.0)

Opitmize the g_truth through the forward projection

g_truth= torch.zeros([1,4,512,512])
g_truth.requires_grad = True
f_estimated = self.projector(g_truth)
loss = self.loss_func(f_estimated.cpu().float(), torch.ones([1,360, 8, 512]).cpu().float())
loss.backward()

Opitmize the f_truth through the FBP reconstruction

However, if I use the FBP, the pytorch optimize is unable to find any gradient:
f_truth= torch.zeros([1,360, 8, 512])
f_truth.requires_grad = True
g_estimated = self.projector.fbp(f_truth)
loss = self.loss_func(g_estimated.cpu().float(), torch.ones([1,4,512,512]).cpu().float())
loss.backward()

gives "RuntimeError: element 0 of tensors does not require grad and does not have a grad_fn".

Could you help me to see whether the FBP algorithm is not supporting the optimization?

Thank you so much.

Modular geometry

Hello,
Is it possible to use the modular geometry configuration to compute projections pixel-by-pixel independently (Line projection) and combine batches of pixels from different viewing angles? I could not find enough documentation on how to use this option in the library. Could you please share some documentation or an example on this topic?

Thank you in advance.

Best regards,

Simulated fan-beam limited-angle reconstruction

Hi, I want to use LEAP to simulate the fan-beam limited-angle reconstruction, the following is my geomerty.

leapct.set_fanbeam(180, numRows, numCols, pixelSize, pixelSize, 0.5*(numRows-1), 0.5*(numCols-1), leapct.setAngleArray(numAngles, 90.0), sod, sdd)

The complete scanning process encompasses 720 views that are evenly distributed across 360 degrees. I want simulate the projection data with 180 views over a 90-degree span. But it will rasie the error

setParkerWeights: Not enough data!

Previouslt in astra toolbox it is easy to simulate like

limited_projGeom180 = astra.create_proj_geom('fanflat', 2.0, 768, np.linspace(0, np.pi/2, 180, endpoint=False), 512, 512)

Is there any similar way to do in LEAP? Thanks!

Helical reconstruction

Hello,

Thank you for the excellent library you have developed.

I have reviewed the Python demo code you wrote for helical (d03_helical.py), and it seems to only perform an FBP for reconstruction. However, I noticed that there are other functions written for helical in the C++ source.

Could you please provide a C++ source for helical reconstruction that I can follow? Additionally, I read the LEAP.pdf file but did not fully understand which equation was used for implementing helical reconstruction with cone and curve detector geometry.

Thank you for your assistance.

helical fbp runs slow

Hey! LEAP is soooo impressive! So many features. And it produces great results! Thanks for sharing this with us!! One little problem I have noticed is the helical FBP runs much slower than your other FBP algorithms. I mean the helical projectors are fast but the helical FBP is slow. Is there a reason for this? Can it be faster?

Is the FP and FBP part of the algorithm in LEAP based on ray-driven, pixel-driven, or distance-driven?

RWLS not working

Hey! RWLS reconstruction returns all zeros. Could you look into this?

Is there an official Docker image for LEAP?

Hi,

Thank you for your great CT geometry tools! I am really using well in window system.

So my question is, is there a Docker image to install and run on Linux system?

Although there is a way of installing LEAP on Linux using .sh file, I want to utilize LEAP by using docker image.

Is there any way to make a Docker image for LEAP?

Thanks:)

Error when building on Windows@AMD

So, I tried to build the project without success on Windows with AMD EPYC processor and RTX4070 GPU.
I renamed the setup_AMD.py to setup.py and run pip install -v . and failed on building.
Here is the error that caused the failure:

total_variation.cu
  tmpxft_00001eb8_00000000-7_total_variation.cudafe1.cpp
  "C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Tools\MSVC\14.39.33519\bin\HostX86\x64\link.exe" /nologo /INCREMENTAL:NO /LTCG /DLL /MANIFEST:EMBED,ID=2 /MANIFESTUAC:NO /LIBPATH:C:\Users\Reko2\.conda\envs\LEAP\Lib\site-packages\torch\lib "/LIBPATH:C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v12.4\lib\x64" /LIBPATH:C:\Users\Reko2\.conda\envs\LEAP\libs /LIBPATH:C:\Users\Reko2\.conda\envs\LEAP /LIBPATH:C:\Users\Reko2\.conda\envs\LEAP\PCbuild\amd64 "/LIBPATH:C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Tools\MSVC\14.39.33519\ATLMFC\lib\x64" "/LIBPATH:C:\Program Files\Microsoft Visual Studio\2022\Community\VC\Tools\MSVC\14.39.33519\lib\x64" "/LIBPATH:C:\Program Files (x86)\Windows Kits\10\lib\10.0.22621.0\ucrt\x64" "/LIBPATH:C:\Program Files (x86)\Windows Kits\10\\lib\10.0.22621.0\\um\x64" c10.lib torch.lib torch_cpu.lib torch_python.lib cudart.lib c10_cuda.lib torch_cuda.lib /EXPORT:PyInit_leapct build\temp.win-amd64-cpython-312\Release\src\analytic_ray_tracing.obj build\temp.win-amd64-cpython-312\Release\src\bilateral_filter.obj build\temp.win-amd64-cpython-312\Release\src\cpu_utils.obj build\temp.win-amd64-cpython-312\Release\src\cuda_utils.obj build\temp.win-amd64-cpython-312\Release\src\file_io.obj build\temp.win-amd64-cpython-312\Release\src\filtered_backprojection.obj build\temp.win-amd64-cpython-312\Release\src\find_center_cpu.obj build\temp.win-amd64-cpython-312\Release\src\list_of_tomographic_models.obj build\temp.win-amd64-cpython-312\Release\src\matching_pursuit.obj build\temp.win-amd64-cpython-312\Release\src\noise_filters.obj build\temp.win-amd64-cpython-312\Release\src\parameters.obj build\temp.win-amd64-cpython-312\Release\src\phantom.obj build\temp.win-amd64-cpython-312\Release\src\projectors.obj build\temp.win-amd64-cpython-312\Release\src\projectors_Joseph.obj build\temp.win-amd64-cpython-312\Release\src\projectors_Joseph_cpu.obj build\temp.win-amd64-cpython-312\Release\src\projectors_SF.obj build\temp.win-amd64-cpython-312\Release\src\projectors_SF_cpu.obj build\temp.win-amd64-cpython-312\Release\src\projectors_Siddon.obj build\temp.win-amd64-cpython-312\Release\src\projectors_Siddon_cpu.obj build\temp.win-amd64-cpython-312\Release\src\projectors_attenuated.obj build\temp.win-amd64-cpython-312\Release\src\projectors_extendedSF.obj build\temp.win-amd64-cpython-312\Release\src\projectors_symmetric.obj build\temp.win-amd64-cpython-312\Release\src\projectors_symmetric_cpu.obj build\temp.win-amd64-cpython-312\Release\src\ramp_filter.obj build\temp.win-amd64-cpython-312\Release\src\ramp_filter_cpu.obj build\temp.win-amd64-cpython-312\Release\src\ray_weighting.obj build\temp.win-amd64-cpython-312\Release\src\ray_weighting_cpu.obj build\temp.win-amd64-cpython-312\Release\src\rebin.obj build\temp.win-amd64-cpython-312\Release\src\scatter_models.obj build\temp.win-amd64-cpython-312\Release\src\sensitivity.obj build\temp.win-amd64-cpython-312\Release\src\sensitivity_cpu.obj build\temp.win-amd64-cpython-312\Release\src\tomographic_models.obj build\temp.win-amd64-cpython-312\Release\src\tomographic_models_c_interface.obj build\temp.win-amd64-cpython-312\Release\src\total_variation.obj /OUT:build\lib.win-amd64-cpython-312\leapct.cp312-win_amd64.pyd /IMPLIB:build\temp.win-amd64-cpython-312\Release\src\leapct.cp312-win_amd64.lib
  LINK : error LNK2001: unresolved external symbol PyInit_leapct
  build\temp.win-amd64-cpython-312\Release\src\leapct.cp312-win_amd64.lib : fatal error LNK1120: 1 unresolved externals
  error: command 'C:\\Program Files\\Microsoft Visual Studio\\2022\\Community\\VC\\Tools\\MSVC\\14.39.33519\\bin\\HostX86\\x64\\link.exe' failed with exit code 1120
  error: subprocess-exited-with-error

I also tried building the project on WSL on the same machine and it worked, so I can use the tool, but with sacrificing a little bit of performance to the virtualization.

About the the differentiability of BackProjection

Hi, thank you for the amazing work! I really enjoy to use it.
I have a question about the differentiability of Backprojection.

When I tried to calculate the gradient after BackProjection, RuntimeError: One of the differentiated Tensors appears to not have been used in the graph. has occured.

Below is a simple script to reproduce the issue. If “forward_project ” is set to True, there is no problem, but if it is set to False (=using BackProjcection), an error occurs in my environment.
(I am aware that this code, which calculates the loss of the projected image itself, is pointless, but I put it because I thought it would make it easy to understand the problem.)

At first I thought that the computational graph was not generated correctly through backprojection function. Or am I using the library incorrectly?

I would really appreciate it if you could let me know.
Thanks in advance!

import numpy as np
import torch
from leaptorch import Projector

device = torch.device("cuda:0")

# define projector
forward_project = True #Or False
proj = Projector(
    forward_project=forward_project, use_static=True, 
    use_gpu=True, gpu_device=device, batch_size=1)

# set geometory
numCols = 256
numAngles = 2*int(360*numCols/1024)
pixelSize = 0.5*512/numCols
numRows = 1
proj.leapct.set_parallelbeam(
    numAngles, numRows, numCols, pixelSize, pixelSize, 
    0.5*(numRows-1), 0.5*(numCols-1),
    proj.leapct.setAngleArray(numAngles, 180.0))
proj.leapct.set_default_volume()
proj.allocate_batch_data()

# generate input(dummy)
rng = np.random.default_rng()
if forward_project:
  x = rng.random(size=(1,256,256)) # size of the image
else:
  x = rng.random(size=(180,1,256)) # size of the projection 
x = torch.from_numpy(x).to(device).unsqueeze(0)
x.requires_grad = True

# get projection
y = proj(x)

# calculate gradient
grad = torch.autograd.grad(y.mean(), x, retain_graph=True, create_graph=True)[0]
print(grad, grad.sum())

meet some problem when install LEAP

Hello developer, I want to install the LEAP function library on my X86-linux server. After downloading and changing the current folder with cd, I executed pip install . and encountered an error: ERROR: Could not build wheels for leapct, which is required to install pyproject.toml-based projects. Can you provide a solution?

Failed installation on window

CMake Deprecation Warning at CMakeLists.txt:1 (cmake_minimum_required):
Compatibility with CMake < 3.5 will be removed from a future version of
CMake.

Update the VERSION argument value or use a ... suffix to tell
CMake that the project does not need compatibility with older versions.

CMake Error at CMakeLists.txt:3 (project):
Generator

Visual Studio 16 2019

could not find any instance of Visual Studio.

-- Configuring incomplete, errors occurred!

C:\Users\HL277986\Documents\GitHub\LEAP\win_build>cmake --build . --config Release
CMake Error:
Generator

Visual Studio 16 2019

could not find any instance of Visual Studio.

Failed to compile!

LEAP installation Issue

Hi!

I am working for building a docker file and a docker image for LEAP.

But things got wrong on installation on LEAP on the current version.

Despite I'm using Leap very well on my local computer, and I'm testing for several environments on several server, the installation of LEAP occurs the same error on every attempts.

When executing pip install . , the error below outcame.

root@workspace-2hw8cpn8igqf-0:~/LEAP# pip install .
Processing /root/LEAP
Preparing metadata (setup.py) ... done
ERROR: No .egg-info directory found in /tmp/pip-pip-egg-info-bok50650

Earlier versions of Leap (like versions from Feb 2024 or late 2023) didn't shown this kind of error, so something might gone wrong when fixing other problem from then.

I hope this problem will go well.

Questions about rebin function used in the curved detector

Thanks for your excellent work, it is the first work I've seen support curved detector for FDK reconstruction. So I tried it but the result seems not good. I'm wondering if I missed some part. The following is the code I'm used to reconstructed,

If I didn't add the following lines, it worked but with some artifacts.
'leapct.set_curvedDetector()
fanAngles = np.array(range(numCols),dtype=np.float32)
fanAngles *= 180.0/np.pi
leapct.rebin_curved(g, fanAngles, order=6)
leapct.FBP(g,f)`
After added the rebin code, the final image will be blanked like the following:

Is there any way to solve this problem?
Thanks very much :)

long height sample case with leap's function

@kylechampley. At the case of long length sample as like cylinder battery, there is some tile issue of pivot which is the point of rotation axis at the sample. In other words, the true CoR of the center of detector is different from the bottom one. Such case, I can not overcome this issue with your new detector rotation function. Could you give me any alternative solution at this case?
As you can see below, the center slice is good but the bottom slice has wrong points at the boundary of circle.

1. center slice:
   

2. bottom slice:
   

Cylinder type Sample Auto-calibration Issue

I have checked your new v1.13 auto-cal function with my 4680 battery projection. And I can see that there is some issue at that case.

1. Center slice:

2. Bottom slice:

3. console log:

I will send my full test script code and my sample image location by your e-mail. I guess that the issue is related to modular-beam converting action because the if np.abs(res.x[2]) > 0.1: condition matched and go into below sub-routine, then this issue happens.

Best regards.

Will the LEAP add more iterative algorithms like TIGRE?

Thanks for your work. I just wondering if LEAP will support more iterative algorithms like OSART-TV and the following in the future.

Iterative algorithms

Gradient-based algorithms (SART, OS-SART, SIRT, ASD-POCS, OS-ASD-POCS, B-ASD-POCS-β, PCSD, AwPCSD, Aw-ASD-POCS) with multiple tuning parameters (Nesterov acceleration, initialization, parameter reduction, ...)

Krylov subspace algorithms (CGLS, LSQR, hybrid LSQR, LSMR, IRN-TV-CGLS, hybrid-fLSQR-TV, AB/BA-GMRES)

Statistical reconstruction (MLEM)

Variational methods (FISTA, OSSART-TV)

Detector rotation functions for three cases(skew, tilt, rolling)

@kylechampley I checked your detector rotation function . And I am wondering why this function combined with two more factors, for example, tilt+skew, or skew+roll.

Detector is fixed into a equipment mechanically at most case. But there are sometime some mismatching angles of three factors. So from the real projection image, we need to extract these error angles for exact reconstruction. For this case, we may needs three separated functions for geometry correction.

I think that at least there needs to exist of three separated functions for tilt-rotate(), skew-rotate() and roll-rotate() of detector.
How about your opinion ?

.

Bug related to subsets partition

Hello,
In class subsetParameters and function breakIntoSubsets, subsets are defined like g[m:-1:numSubsets,:,:].
However, for the last subset, the last element will be dropped in python, which makes an inconsistent subset.
For example,

x = range(16)
print(list(x[3:-1:4]))

gives [3, 7, 11].

Geometry of the cone-beam projection

Hi I tried the cone-beam projection in LEAP, the detector had 200 rows and the centerRow is 100.5 (see below):

proj = Projector(use_gpu=True, gpu_device=device)
proj.leapct.set_volume(256, 256, 100, 1, 1, 0, 0, 0)
proj.leapct.set_conebeam(360, 200, 256, 1, 1.5, 100.5, 128.5, proj.leapct.setAngleArray(360, 360.0), 1000, 1500, tau=0.0, helicalPitch=0.0)
...
sino = proj.leapct.project(sino,img)
plt.imshow(sino[100,:,:].cpu().numpy(), cmap='gray')

However, the displayed projection in the 101 view looks like:

Seems that the parameter of centerRow doesn't work properly?