Helical reconstruction

[MshGhazale]

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.

[kylechampley]

Hello.

Yes, that demo script only runs the helical FBP, but all iterative reconstruction algorithm in LEAP can be applied to all geometries. If you look further down the script, you will see a bunch of iterative reconstruction algorithms that are commented out. You can remove the comment characters and try any of these out.

The helical FBP algorithm is scattered through multiple cu and cpp files. For example, there are files ramp_filter.cu and ramp_filter_cpu.cpp that perform GPU-based and CPU-based ramp filter for each geometry. The ray weighting is in the ray_weighting.cu and ray_weighting_cpu.cpp files. The backprojection is in projectors_SF.cu, projectors_SF_cpu.cpp, or projectors_extendedSF.cu.

The derivative step in helical FBP is implemented in the function parallelRay_derivative in ramp_filter.cu And the Hilbert filter step is implemented in the function: Hilbert1D is ramp_filter.cu

It is fairly widely accepted that the best helical FBP reconstruction algorithm is the WFBP algorithm. But this algorithm applies to cone-parallel coords which is a rebinning of the native cone-beam coordinates where the transaxial rays are rebinned to be parallel.

I did not want to perform this rebinning, so I used the "weighted backprojection trick" from the WFBP algorithm and used it in a helical FBP algorithm in native coordinates given here.

I very briefly mention this on pages 14-15 in the LEAP manual.

[oooo1114]

Hello, thanks for the excellent library. I have tried FBP with an sinogram, the geometry looks fine since all details are sharp. But I have encountered with some helical related artifacts related with in-homogeneous contrast. As shown in the following image, the right part seems darker while the left part seems brighter. And this effect is in a helical manner as I scroll up and down slices. Is the artifact related to the algorithm in native coordinates? Thanks a lot.

[kylechampley]

Medical CT vendors usually perform helical FBP in cone-parallel coordinates because it provides reconstructions where the noise and resolution are more uniform across the field of view. There is nothing inherent about helical FBP in native coordinates that would give stronger helical artifacts.

I do see the intensity difference you speak about. My best guess is that this is caused by 1) a bug in LEAP 2) something is wrong with the way you specified the geometry 3) something is wrong with the data

If you send me your data and a script to reproduce the result, I'll look into it.

[oooo1114]

Thanks for your reply. The data is too big to upload. I reproduce the bug with my geometry on the FORBILD phantom. Here are the codes:

import time
import numpy as np
from leapctype import *
leapct = tomographicModels()

numCols = 860
numTurns = 6
numAngles = 8640
numRows = 64
pixelHeight = 1.0947
pixelWidth = 1.0178
centerRow = 0.5*(numRows-1)
centerCol = 0.5*(numCols-1)+4.25048
phis = leapct.setAngleArray(numAngles, 360.0*numTurns)
sod = 570
sdd = 990
pitch = 0.843
leapct.set_conebeam(numAngles, numRows, numCols, pixelHeight, pixelWidth, centerRow, centerCol, phis, sod, sdd)
leapct.set_normalizedHelicalPitch(pitch)
leapct.set_curvedDetector()

numX = 768
numY = 768
numZ = 137
voxelWidth = 0.77142188
voxelHeight = 1.26
offsetX = 0
offsetY = 0
offsetZ = 0
leapct.set_volume(numX, numY, numZ, voxelWidth, voxelHeight, offsetX, offsetY, offsetZ)
leapct.set_diameterFOV(640)
leapct.print_parameters()

g = leapct.allocateProjections()
f = leapct.allocateVolume()
leapct.set_FORBILD(f,True)
startTime = time.time()
leapct.project(g,f)
print('Forward Projection Elapsed Time: ' + str(time.time()-startTime))
startTime = time.time()
leapct.FBP(g,f)
print('Reconstruction Elapsed Time: ' + str(time.time()-startTime))

Here are the reconstruction volume in three views to better visualized the artifacts with a narrow window near the air.

[kylechampley]

I really appreciate the self-contained script that shows the artifact! Yes, that is definitely a bug. I’ll resolve it and let you know when it’s fixed.

In the meantime, you could probably mitigate that artifact with just a few iterations of an iterative reconstruction method because that is a low frequency error.

[oooo1114]

Thanks, I will try.

[kylechampley]

I was able to find and fix that bug in the helical FBP reconstruction. I'll merge the changes into the main branch and make a new release sometime this weekend, but if you want early access to the bug fix, you can checkout the code from the champley_dev branch.

[MshGhazale]

I encountered an issue while using the FBP function in the d03_helical.py script, which runs on the CPU. The source file filtered_backprojection.cpp has an execute function containing the following code snippet: bool filteredBackprojection::execute(float g, float f, parameters* params, bool data_on_cpu) { ... if (params->geometry == parameters::CONE && params->helicalPitch != 0.0 && params->whichGPU < 0) { printf("Error: CPU-based FBP not yet implemented for helical cone-beam geometry\n"); return false; } ... } As indicated by the error message: Error: CPU-based FBP not yet implemented for helical cone-beam geometry Since my geometry is CONE BEAM and involves helical pitch, this message suggests that the CPU-based FBP is not yet supported for this configuration. Could you please provide guidance on how to proceed with CPU-based FBP for helical cone-beam geometry?

[kylechampley]

At this time LEAP does not support CPU-based helical FBP. I forgot that this never got finished. I can put this on a list of priorities for improvements to LEAP, but honestly it will take me some time to get to this. Helical FBP is a non-trivial algorithm and it will take some time to write up a CPU-based algorithm that runs efficiently.

[kylechampley]

The bug fixes to the helical FBP algorithm have been merged into the main branch and release v1.14 generated.

Let me know how it goes.

[MshGhazale]

Does LEAP currently offer support for a GPU-based implementation of the helical FBP algorithm in its C++ functions?

[kylechampley]

Yes, there are GPU-based helical FBP algorithms in LEAP. Maybe I don't understand your question because I thought I already mentioned this.

[oooo1114]

The bug fixes to the helical FBP algorithm have been merged into the main branch and release v1.14 generated.

Let me know how it goes.

Thanks, I have tested and the bug is fixed.

[kylechampley]

Great! I'm going to close this issue, but feel free to open a new one if something comes up.

[MshGhazale]

I am a beginner in CT image reconstruction and encountered a CUDA error while using the FBP_gpu(g, f) function for helical reconstruction. Environment: OS: Windows 10 CUDA Version:11.8 PyTorch Version: 2.0.1 import sys import os import time import numpy as np from leapctype import * import torch import nibabel as nib

leapct = tomographicModels() os.environ["CUDA_LAUNCH_BLOCKING"] = "1"

numCols = 512 numTurns = 10 numAngles = 2 2 int(360 numCols / 1024) numTurns pixelSize = 0.65 * 512 / numCols

numRows = numCols // 4 leapct.set_GPU(0) leapct.set_conebeam(numAngles, numRows, numCols, pixelSize, pixelSize, 0.5 (numRows - 1), 0.5 (numCols - 1), leapct.setAngleArray(numAngles, 360.0 * numTurns), 1100, 1400) leapct.set_curvedDetector() leapct.set_normalizedHelicalPitch(0.5) leapct.set_default_volume()

f = leapct.allocateVolume() g = leapct.allocateProjections_gpu() leapct.set_FORBILD(f, True) f = torch.from_numpy(f)

startTime = time.time() leapct.project(g, f)

try: leapct.FBP_gpu(g, f) except RuntimeError as e: print(f"Error during FBP_gpu: {e}") print('Reconstruction Elapsed Time: ' + str(time.time() - startTime))

Error Message: Torch is on cuda: True leapct.print_parameters() cudaMemcpy3D Error: invalid argument kernel failed! error name: cudaErrorIllegalAddress error msg: an illegal memory access was encountered error during FBP_gpu: CUDA error: an illegal memory access was encountered Compile with TORCH_USE_CUDA_DSA to enable device-side assertions.

Traceback (most recent call last): File "E:\helical\leap\LEAP-1.12\test_helical1.py", line 85, in leapct.FBP_gpu(g, f) File "E:\helical\leap\LEAP-1.12\src\leapctype.py", line 1965, in FBP_gpu q = self.copyData(g) File "E:\helical\leap\LEAP-1.12\src\leapctype.py", line 1407, in copyData x_copy = x.clone() RuntimeError: CUDA error: an illegal memory access was encountered Compile with TORCH_USE_CUDA_DSA to enable device-side assertions. Could you guide me on how to resolve this illegal memory access error during the FBP_gpu function call?

[kylechampley]

You got an error because LEAP requires that the projection and volume data must either both be on the CPU or both be on the GPU. In your case f was on the CPU, so all you have to do is copy it to the GPU before calling any LEAP functions.

[kylechampley]

@MshGhazale has this issue been resolved for you? If so, let me know so I can close this issue.

[MshGhazale]

Thank you. Yes, I moved volume data to the GPU and the error was resolved.