search

LLNL / LEAP

comprehensive library of 3D transmission Computed Tomography (CT) algorithms with Python API and fully integrated with PyTorch
https://leapct.readthedocs.io
MIT License
74 stars 8 forks source link

Pytorch optimization through FBP #25

Closed shaoyanpan closed 3 months ago

shaoyanpan commented 3 months ago

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.

kylechampley commented 3 months ago

Hello.

Before I get into your main questions, I'd like to discuss some other things.

First of all, it seems that in your set_conebeam function you are setting your center row to 256, but you only have 4 detector rows. You can do this, but this would mean that your detector is shifted very far up and I assume this is not what you wanted. A perfectly centered detector would have this value at (4-1)/2.

Next, we usually use the variable "f" for the reconstruction volume and "g" for the projection data. It seems like you have these switched in your code above which is a little confusing to me. You can use whatever variable names you want, but if you switch these names it might be hard to read our code and it may contribute to some of the problems you are experiencing.

Lastly, I am not exactly sure what you are trying to do, but I'll make some comments, but feel free to write back if I did not answer your question. OK, we do provide an implementation of FBP, but it is not setup to be a differentiable operator in our package. If you check our Projector class, only our forward and backprojections are differentiable. If you are interested in the derivative of FBP, you will have to implement this yourself, but it is not trivial to do. You could implement the adjoint of FBP as the following.

proj.leapct.project(...) proj.leapct.filterProjections(...)

shaoyanpan commented 3 months ago

Hi:

I think your answer is perfect actually. I will try the adjoint of FBP (I see them in the instruction) so I can calculate the derivative of FBP. Thank you so much!

Also I notice your suggestions and they are very valuable. Thank you so much!