mmuckley / torchkbnufft

A high-level, easy-to-deploy non-uniform Fast Fourier Transform in PyTorch.
https://torchkbnufft.readthedocs.io/
MIT License
201 stars 44 forks source link

Toeplitz Kernel has boundary artifacts #85

Closed roflmaostc closed 1 year ago

roflmaostc commented 1 year ago

Hi,

I'm not quite familiar with the mathematics of the kernel but it looks like that the kernel has some boundary issues. I would have expected that it decreases with intensities but it doesn't.

This is for the use-case of radial polar patterns.

What you see is the real part of the kernel, after fftshift.

Did you observe that before?

Best,

Felix

Screenshot_2023-04-18_18-32-16

mmuckley commented 1 year ago

Hello @roflmaostc, I think there are some zeros that are manually inserted in toep.py. It's been quite some time since I worked with this to remember the justification, but it is used in other reference implementations. For details I think you'll need to look at the Feichtinger reference.

If you're talking about the ring boundary, that should be due to using a radial trajectory. Those are the k-space locations that are not sampled.

roflmaostc commented 1 year ago

Thanks, I'll check that out!

The kernel is in real space. I'm more worried about the strong ringings in the corners. As you said, it might come from the hard cut off in the Fourier space. But in my application it causes some background blur which smears into the whole image, not so nice.

I think I might just zero them out in real space, because those are exactly the regions which are also not measured at all.

mmuckley commented 1 year ago

Hmm, yes, it can be a little bit non-intuitive since the matrix A is being embedded into a larger 2-factor size matrix. Considering cropping, I wonder if convolving with the kernel you show will touch the ringing locations. Are you using iterative reconstruction?

If you notice anything inconsistent about toep.py with Toeplitz embedding feel free to open a PR.

roflmaostc commented 1 year ago

So my input was first quadratic and I used iteratively. Then the ringing locations will smear data into a circle of a diameter with array size, so will mess up things.

But if you set all data outside of this circle to zero in every iteration, then it is fine.

This actually kind of makes sense since a CT system would not be able to measure the corners of the array anyways (because of finite detector size).

mmuckley commented 1 year ago

I see. Did you use regularization? Technically you should get ringing because those corners are in the nullspace of the operator and you're trying to invert 0s without regularization.

Setting data outside the circle to zero is an equivalent form of regularization via constraining, I guess applying projected gradient descent.

roflmaostc commented 1 year ago

Yes, that makes sense for me.

No, I didn't regularize. My approach is similar t IFTA with manually setting constraints in each iteration.

roflmaostc commented 1 year ago

Btw, because of the strong ringing, convolving with the Toeplitz kernel is also not energy preserving.

A lot of energy is lost in the ringing and hence cut-off.

mmuckley commented 1 year ago

By "energy-preserving" do you mean norm(x) = norm(toep(x))? That is also expected. The NUFFT is not an energy-preserving operation. The Toeplitz approximates A'A where A is a NUFFT. The NUFFT can be made approximately energy-preserving if you apply density compensation (analogous to ramp filters in CT).

roflmaostc commented 1 year ago

Yeah, or simply sum(conv(x, toeplitz_matrix)).

Even if you normalize so that sum(toeplitz_matrix)=1, the operation is not energy preserving because a lot of the energy is distributed in the ringing which is then cut-off