Use the projection residual to get the value that the current reconstructed image should update

[Kwong-Yam-Lie]

Hello,

I'm implementing an algorithm .

I want to "Use the projection residual to get the value that the current reconstructed image should update" , but I falied. For example, with SART-based algs. I find that all algorithms return non-negative results .

So when I want to get the value that the current reconstructed image should update, I got the wrong result.

# Load thorax phatom data
head = sample_loader.load_head_phantom(geo.nVoxel)
# generate projections
projections = tigre.Ax(head, geo, angles)
# SART
imgSART = algs.sart(-projections, geo, angles, niter = 10) 
## Just a test, the 'projections' is  negative but 'imgSART' is filled with zeros.  

Excuse me how can I realize my requirements by calling this software?

Does the software have this function?
I need your help, thanks!

[AnderBiguri]

SART does that for you inside, the algorithms are there so you don't need to write what you want.

In fact, SART is an algorithm that does exactly that, for you, for niter updates. If you want to reproduce it yourself, I suggest you read the source code of SART.

An example of computing the projection residual is

x=#some initial guess
projection_residual = tigre.Ax(x,geo,angles)-projections

Updating the image with that requires a bit more. The code is inside the algorithms.

[Kwong-Yam-Lie]

Thank you for your patience. 🔥

I know that backprojection needs to use texture memory. I have also implemented some 2D algorithms using CUDA and Cplusplus.

Usually in order to speed up the iteration, we will zero the non-negative voxels in the result of each iteration. This is very important ！But now I just want to get the update value ($img_{residual}^{k}$) of the voxel in 3D scene, which may be a negative value.🏁
$$img^{k+1} = img^{k} + \lambda*img^{k}_{residual}$$
$$img^{k+1}[img^{k+1}<0] = 0$$

I guess there may not be a call interface for this small module in the python code.
So I think I may need to directly call the Cplusplus source code or write it myself. 🤔

I'm working on X-Ray dual spectral CT. This software provides me with great convenience👍. Thanks again！

[AnderBiguri]

I believe imgSART = algs.sart(projections, geo, angles, niter = 10, nonneg= False) should remove the non-negativity constraint, if that is what you are looking for.

But again, getting img_residual is just a line that uses the backprojection and forward projection.
You can find the exact line here:

TIGRE/Python/tigre/algorithms/iterative_recon_alg.py

Lines 385 to 397 in 8df6326

	self.res += (
	self.lmbda
	* 1.0
	/ self.V[iteration]
	* Atb(
	self.W[ang_index]
	* (self.proj[ang_index] - Ax(self.res, geo, angle, "Siddon", gpuids=self.gpuids)),
	geo,
	angle,
	"FDK",
	gpuids=self.gpuids,
	)
	)

or you can have a simple non-weighted one by doing:

img_res=tigre.Atb(projections-tigre.Ax(x,geo,angles),geo,angles)

The purpose of TIGRE is to make you not use your cpp code, because we provide it. Atb is backprojection, Ax is projection.

[Kwong-Yam-Lie]

With your suggestion, my problem has been solved！ 😄