将ista算法转化为图像处理代码（压缩感知磁共振）

将ISTA针对图像问题进行展开，并用Pytorch编程的过程进行说明。

针对图像压缩感知（逆问题）问题（1）：

\[\underset{x}\min{\frac{1}{2}|| \Phi x - y ||_2^2} + \lambda|| \Psi x || \]

这里给出迭代解ISTA的表示（2）：

\[\left\{{\begin{align} & r^{(k)}=x^{(k-1)}-\rho \Phi^{\mathrm{T}}(\Phi x^{(k-1)}-y) \\ & x^{(k)}=\underset{x}{\arg\min}{\frac{1}{2}|| x - r^{(k)} ||_2^2} + \lambda|| \Psi x || \end{align}}\right. \]

然后给出软阈值优化迭代算法（3）：

\[\left\{{\begin{align} & r^{(k)}=x^{(k-1)}-\rho \Phi^{\mathrm{T}}(\Phi x^{(k-1)}-y) \\ & x^{(k)}=G(\mathrm{soft}(F(r^{(k)},\theta))) \end{align}}\right. \]

将（3）展开成迭代网络，图下给出ISTA-Net的示意图：

可以看出结构还是很清晰的。

这里参考了论文 “基于迭代网络的快速磁共振成像 南昌大学 刘沂玲 2019” 和 ISTA-NET 原文及其代码。

将（3）变形成（4）：

\[\left\{{\begin{align} & r^{(k)}=x^{(k-1)}-\rho \Phi^{\mathrm{T}}\Phi x^{(k-1)}+\rho \Phi^{\mathrm{T}} y \\ & x^{(k)}=G(\mathrm{soft}(F(r^{(k)},\theta))) \end{align}}\right. \]

这样就得到了我们编程中用的表达式。下面将重点讲解里面的每一个符号。\(\rho\) 表示一个学习率，通常迭代中取常数利普西斯 \(L\) 的导数，但这里的展开网络可以随意一定。\(\Phi^{\mathrm{T}}\) 表示观测矩阵的转置，\(\Phi\) 表示观测矩阵。假设我们的问题是CS-MRI （compressed sensing MRI），压缩感知磁共振成像问题的话，那么观察矩阵 \(\Phi\) 就意味着在傅里叶空间进行下采样。\(\Phi=\mathrm{MF}\)，\(\mathrm{M}\) 表示下采样矩阵，\(\mathrm{F}\) 表示傅里叶变换矩阵。

但是，如果这里很老实的用 \(\Phi=\mathrm{MF}\) 这种表达形式来做压缩感知问题，计算量会特别大，因为此时 \(y\) 会被拉成一个列向量。所以还是使用fft 函数最合适。下面一一来分析怎么使用fft 来对应 \(\Phi=\mathrm{MF}\)。

对于 \(\Phi x^{(k-1)}\)，我们知道 \(x^{(k-1)}\) 就是一个迭代中的图像，那么在一个观察矩阵 \(\Phi\) 下，\(\Phi x^{(k-1)}\) 就是迭代图像的\(k\) 空间数据乘以观测\(\mathrm{M}\)，\(\Phi^{\mathrm{T}}\Phi x^{(k-1)}\) 就是用重建出的k空间数据得到的图像。所以， $\Phi x^{(k-1)} $ = fft(x(k-1)*mask), \(\Phi^{\mathrm{T}}\Phi x^{(k-1)}\) = ifft(fft(x(k-1))*mask). 对于 \(\Phi^{\mathrm{T}} y\) 其实也一样，它的值其实就是我们常说的zero-filling image。不同的是，\(y\) 是已观察的数据，是一个常数了，那么 \(\Phi^{\mathrm{T}} y\) 也是一个常数了，所以编程里面，常常就在初试阶段就已经求好，不需要重复求值。在ISTANET 的官方代码里面，用的是 \(\Phi^{\mathrm{T}} b\) 来表示的。假设有全采样k空间数据：f_k。那么：\(\Phi^{\mathrm{T}} b\) = ifft(mask * f_k)。

最后，我们来对一下ISTANET 中的官方代码和公式里面的符号：

\[\left\{{\begin{align} & \mathrm{PhiTPhi}= \Phi^{\mathrm{T}}\Phi \\ & \mathrm{FFT-Mask-ForBack}=zerofilingimg=\Phi^{\mathrm{T}} b \\ & \mathrm{\Phi}=fft(mask*) \\ & \mathrm{\Phi^{\mathrm{T}}}=ifft(mask) \\ & \mathrm{\Phi^{-1}}=ifft() \end{align}}\right. \]

对于\(\mathrm{\Phi^{-1}}=ifft()\) 的原因，1. FH 是逆傅里叶变换的定义 这里给出了一些定义。2. denotes the Hermitian of F, 3. 证明F-1 = FH。为什么要用 \(\mathrm{\Phi^{-1}}\)呢？因为在 ADMM-LIKE 的网络中，如 “基于迭代网络的快速磁共振成像 南昌大学 刘沂玲 2019” 【IFR-Net】 会由这个 \(\mathrm{\Phi^{-1}}\) 的中间结果。在 IFR-Net 的代码中，https://github.com/yqx7150/IFR-Net-Code，In IFR-Net-Code/layersfunction/xorg.m 中，line15 也说明了$F^{-1} $ 是 ifft。O = real(ifft2(Y .* Q1));。

对于\(\mathrm{\Phi^{-1}}=ifft()\) 和 \(\mathrm{\Phi^{\mathrm{T}}}=ifft(mask)\) 的区别，个人认为有以下不同：\(\mathrm{\Phi^{-1}}\) 是直接的傅里叶逆变换，而 \(\mathrm{\Phi^{\mathrm{T}}}\) 是观测的傅里叶逆变换。

所以对于公式（4），可以用代码写出:

\[\left\{{\begin{align} & r^{(k)}=r^{(k-1)}-\rho\cdot ifft(fft(r^{(k)})\cdot mask)+\rho\cdot ifft(fullk\cdot mask) \\ & x^{(k)}=G(\mathrm{soft}(F(r^{(k)},\theta))) \end{align}}\right. \]

G 和F 分别是两个一样的两层卷积神经网络，并且是可逆的。需要用loss来约束G 和F，来使得他们互为反函数。

End.