编程语言
首页 > 编程语言> > 将ista算法转化为图像处理代码(压缩感知磁共振)

将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的示意图:

image-20220518171420329

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

这里参考了论文 “基于迭代网络的快速磁共振成像 南昌大学 刘沂玲 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} $ 是 ifftO = 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.

标签:ista,Phi,ifft,迭代,align,mask,图像处理,磁共振,mathrm
来源: https://www.cnblogs.com/paprikatree/p/16288242.html