[Implementation question] XPDnet

[zwep]

Sorry for not sticking to the Issue convention. I am a little bit hesitant to ask for help like this, because I think Im making a silly mistake somewhere... But Im curious if you can remark on the implementation of the CrossDomainNetwork class.

There, in the kspace_correction method of the CrossDomainNetwork class, I see the following update

            kspace_buffer = torch.cat([kspace_buffer, forward_buffer, masked_kspace], self._complex_dim)

            ...

            kspace_buffer = kspace_buffer[..., :2] - kspace_buffer[..., 2:4]

I have read the paper on the DPHG algorithm (https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=8271999) and the XPDNet (https://arxiv.org/pdf/2010.07290.pdf), but I couldnt see where this difference is coming from.

If I understand it correctly, we are taking the difference between the kspace_buffer (== initial kspace at the first iteration) and the forward_buffer (==kspace transformed from image space input at the first iteration). So this would end up being a very noisy image, since both the kspace_buffer and forward_buffer are very similar, if not identical. And even if this is not the case for the first iteration, my understanding is that ultimately the kspace_buffer and forward_buffer should start to convergence to each other.

How does this then become an appropriate input for the image_correction method? I believe that I am making an error somewhere with the interpretation of the variables, but I cannot see my error at this point. So maybe one of you can help me with this. If you think responding this will take too much time please say so, then I'll think some more about it

[zwep]

I closed this because the numerical experiments showed at least proper image creation. So in practice it works fine I guess? But I am not able (limited time as well) to explain it why it works the way it does.

[zwep]

Okay I went over it one more time. These are my findings

1. I am not certain that taking the difference between the k-space is an appropriate replacement of the learnable operator in the original algorithm. What do you think about taking the average between the two spaces instead?
2. Taking the difference between the k-space arrays does result in images with only high frequency components. Meaning that the difference is relatively small. However, in further iterations this effect diminishes and the whole process corrects for this.
3. Thus, taking the difference between the k-space arrays does not result in an array with zeros (or even machine precision). Despite the fact that the following operations are applied (I have simplified it to a mathematical notation)

$$X^{(1)} = X^{(0)} - \hat{X}^{(1)}$$ $$\text{This is meant to be identical to: kspace\_buffer = kspace\_buffer[..., :2] - kspace\_buffer[..., 2:4]}$$ Where

$$\hat{X}^{(1)} = F(Y^{(0)}, S, P)$$ Where

$$Y^{(0)} = F^{-1}(X^{(0)}, S, P)$$

Note here that $X^{(i)}$ represent the k-space array at iteration $i$, analogous for the image array $Y^{(i)}$. The operator $F(..., S, P)$ represents the Fourier transform using sensitivity map $S$ and sampling mask $P$.

In my idea, the above equations should simplify to: $$X^{(1)} = X^{(0)} - \hat{X}^{(1)} \rightarrow X^{(0)} - F(Y^{(0)}, S, P) \rightarrow X^{(0)} - F(F^{-1}(X^{(0)}, S, P), S, P) \rightarrow X^{(0)} - X^{(0)} = 0$$

But in the implementation I experience that applying the _backward and _forward operator in succession does not result in the identity operator.

Please correct me where I am wrong if you have the time.

[georgeyiasemis]

Hi @zwep. Thanks for opening this. ATM I am not sure I fully understand your question, but in the next weeks I will try to go through the code to see if there is indeed some bug, or something.

A few notes:

• _backward_operator and _forward operator should not result in the identity. They are adjoints, not inverses. In accelerated MRI, an inverse to the forward operator does not exist (ill-posed problem due to sub-sampling).
• Code for XPDNet was adapted from https://github.com/zaccharieramzi/fastmri-reproducible-benchmark/blob/master/fastmri_recon/models/subclassed_models/cross_domain.py.

Hope the above helps

[zwep]

Ah okay thanks, I get it. The goal is really to use the residual of the k-space. I thought that the k-space refinement stage should eventually produce an estimate of the fullly sampled k-space.... but I think I understand now why approximating the residual can also be useful