Toeplitz Kernel has boundary artifacts about torchkbnufft

[roflmaostc]

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]

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]

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]

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]

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]

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]

Yes, that makes sense for me.

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

[mmuckley]

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]

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