Improved Resolution from Subpixel Shifted Pictures 论文理解

references【前置知识】

Theoretical characterizations of images that can be totally refined using the transducer oscillation algorithm

这篇文章据说就是用解方程做的，但是根本找不到文章本体x

Improving image resolution using subpixel motion

用local search+退火做的

Generalized Sampling Expansion（比较有用）

数学上证明了：对于\(f(t)\)，如果傅里叶变换满足：

\[\begin{gathered} F(\omega)=0, \quad \text { for } \quad|\omega| \geqslant \sigma \\ \int_{-\infty}^{\infty}|f(t)|^{2} d t=\frac{1}{2 \pi} \int_{-\sigma}^{\sigma}|F(\omega)|^{2} d \omega<\infty \end{gathered} \]

当用m个线性系统\(H_{1}(\omega), \cdots, H_{m}(\omega)\)去处理\(f(t)\)得到：

\[g_{k}(t)=\frac{1}{2 \pi} \int_{-\sigma}^{\sigma} F(\omega) H_k(\omega) e^{j \omega t} d \omega \]

然后对\(g_k(t)\)进行采样得到\(g_k(nT)\)，其中\(T=\frac{2 \pi}{c} \quad c=\frac{2 \sigma}{m}\)，则\(f(t)\)可以用\(g_k(nT)\)来表示。过程为：

先列一个方程组：

\[\begin{aligned} &H_{1}(\omega) Y_{1}(\omega, t)+\cdots+H_{m}(\omega) Y_{m}(\omega, t)=1 \\ &\ldots \\ &H_{1}(\omega+r c) Y_{1}(\omega, t)+\cdots+H_{m}(\omega+r c) Y_{m}(\omega, t)=e^{j r c t} \\ &\ldots \\ &H_{1}(\omega+m c-c) Y_{1}(\omega, t)+\cdots+H_{m}(\omega+m c-c) Y_{m}(\omega, t)=e^{j(m-1) c t} \end{aligned} \]

where \(t\) is any number and \(\omega\) is between \(-\sigma\) and \(-\sigma+c\).

这样可以定义m个函数：\(Y_{1}(\omega, t), \cdots, Y_{m}(\omega, t)\)

然后可以用它们来表示\(f(t)\)：

\[\begin{array}{r} f(t)=\sum_{n=-\infty}^{\infty}\left[g_{1}(n T) y_{1}(t-n T)+\right. \left.\cdots+g_{m}(n T) y_{m}(t-n T)\right] ,\\ where \ \ \ \ \ y_{k}(t)=\frac{1}{c} \int_{-\sigma}^{-\sigma+c} Y_{k}(\omega, t) e^{j \omega t} d \omega, \quad k=1, \cdots, m \end{array} \]

论文本体

重建过程

(1) Determining the model parameters, i.e., the PSF and the relative subpixel shifts in the ensemble \(\left(g_{k}\right)\).
(2) Registering and merging the ensemble over a finer grid to obtain \(\hat{g}\).
(3) Deconvolving the combined picture \((\hat{g})\) with the restoration filter \(\left(y_{I}\right)\) to obtain an improved picture \((\hat{f})\).
如何得到 \(\hat{g}\)？通过下面的分析：

the sampling

\(h_{I}\)和\(h_{II}\)是filter的分解，分别代表了shift和blur，由于是linear space invariant operator，a和b是等价的。对channel的response用傅里叶变换：

\[\begin{aligned} H_{k}(\omega) &=\frac{1}{j \omega}\left(e^{j \omega A / 2}-e^{-j \omega A / 2}\right) e^{-j \omega x_{k}} \\ &=\underbrace{A \frac{\sin (\omega A / 2)}{(\omega A / 2)}}_{H_{l}(\omega)} \cdot \underbrace{e^{-j \omega A / k}}_{H_{I I}(\omega, k)} \end{aligned} \]

A是detector长度。

reference中采样定理的channel：

k个线性channel的output（\(g_k\)）sampled at a rate of \(\omega=2 \sigma / K\)（\(1/K\)of the Nyquist rate）。如何重建：passing the sampled signals through \(K\) linear filters whose impulse responses are \(y_{k}\). The output \(\hat{f}\) is equal to \(f\), and

\[\hat{f}(x)=\sum_{k=1}^{K} \hat{f}_{k}(x)=\sum_{k=1}^{K} \sum_{m=-\infty}^{\infty} g_{k}(m D) \cdot y_{k}(x-m D) \]

\[\begin{aligned} &H(\omega)={\left[\begin{array}{ccc} H_{1}(\omega) & \cdots & H_{K}(\omega) \\ H_{1}(\omega+\Omega) & \cdots & H_{K}(\omega+\Omega) \\ \vdots & & \vdots \\ H_{1}(\omega+(K-1) \Omega) & \cdots & H_{K}(\omega+(K-1) \Omega) \end{array}\right] } \end{aligned} \]

\[y_{k}(k)=\frac{1}{\sigma} \sum_{\nu=1}^{K} \int_{I_{1}} y_{k \nu} e^{j(\omega+(\nu-1) \sigma) x} d \omega \]

where \(y_{k \nu}\) is the \(k \nu^{\text {th }}\) element of the inverse matrix of \(H(\omega)\)
【变量对应采样定理】

我的理解：这里解出的就是\(\hat{g}\)。
论文后面就直接说，shift的closed form表达式：

\[\begin{aligned} &\hat{g}(x)=\sum_{k=1}^{x} \sum_{n=-x}^{\infty} g\left(m D+x_{k}\right) \cdot \frac{\Pi_{n=1}^{K} \sin \left(\Omega\left[\left(x-m D-x_{n}\right) / 2\right]\right)}{\left(\Omega\left[\left(x-m D-x_{k}\right) / 2\right]\right) \cdot \prod_{n=1 \atop n \neq k}^{K} \sin \Omega\left[\left(x_{k}-x_{n}\right) / 2\right]} \end{aligned} \]

也可以写作（是通过另一个定理还是直接写，没看懂）：

\[\begin{aligned} \hat{g}(x)=& \sum_{k=1}^{\infty} \sum_{m=-\infty}^{\infty} g_{k}(m D) \varphi_{k m}\left(x, m, x_{1}, \ldots ., x_{k}, \ldots, x_{K}\right) \end{aligned} \]

为了计算方便，regroup：

\[\hat{g}(x)=\sum_{k=1}^{K} \sum_{m=-\infty}^{\infty} g\left(m D+x_{k}\right) \frac{\sin \Omega\left[\left(x-m D-x_{k}\right) / 2\right]}{\Omega\left[\left(x-m D-x_{k}\right) / 2\right]} \cdot \frac{\prod_{n=1 \atop n \neq k}^{K} \sin \left(\Omega\left[\left(x-m D-x_{n}\right) / 2\right]\right)}{\left.\prod_{n=1 \atop n \neq k}^{K} \sin \Omega\left[\left(x_{k}-x_{n}\right) / 2\right]\right)} . \]

则：\(\varphi_{mk}(x)=\frac{\sin \Omega\left[\left(x-m D-x_{k}\right) / 2\right]}{\Omega\left[\left(x-m D-x_{k}\right) / 2\right]} \cdot \frac{\prod_{n=1 \atop n \neq k}^{K} \sin \left(\Omega\left[\left(x-m D-x_{n}\right) / 2\right]\right)}{\left.\prod_{n=1 \atop n \neq k}^{K} \sin \Omega\left[\left(x_{k}-x_{n}\right) / 2\right]\right)}\)

\(sinc= \frac{\sin \Omega\left[\left(x-m D-x_{k}\right) / 2\right]}{\Omega\left[\left(x-m D-x_{k}\right) / 2\right]}\), \(\varphi_{m k}^{*}(x)=\frac{\prod_{n=1 \atop n \neq k}^{K} \sin \left(\Omega\left[\left(x-m D-x_{n}\right) / 2\right]\right)}{\left.\prod_{n=1 \atop n \neq k}^{K} \sin \Omega\left[\left(x_{k}-x_{n}\right) / 2\right]\right)}\)

我的理解：只是为了方便计算变形了一下。

这篇文章大概的思路了解了，但是没有解决最关心的问题：为什么这个问题要到频域上做，和退火那个做法的根本出发点区别在哪里。不是从数学上知道“噢有这么一个定理可以在频域上做”，而是更有逻辑的直观的“为什么频域上更好”。以及有没有别的正交空间上也有这样的好的性质呢？

说实话我傅里叶变换这节学得很差...其实最近看这些都挺难理解挺痛苦的...

顺便我二刷也没看懂前面说的一些看起来很关键的想法和后面的解法有什么关系。大概是这样的内容：
bandwidth is limited by sampling rate（单位面积的detector数量），限制：intersample distance cannot be smaller than the detector size. -> shifted picture

temporal和special bandwidth的trade-off：对于resolving power超过传统限制的光学系统，fundamental invariant是系统传递的光学图像的N个自由度，并不仅仅与spatial bandwidth有关（是product of the spatial x, y and temporal dimensions）。因此可以通过减小temporal bandwidth来增大spatial bandwidth。

这篇文章考虑的sample theorem：Papoulis considers a \(\sigma\)-Band-Limited (BL) signal that passes through \(K\) linear channels and then sampled at \(1/K\) the Nyquist rate. 这里的linear channel对应的是这篇文章中得到低分辨率图像的过程和重建过程。

把低分辨率图像merge到一张大的图的过程可以认为是周期非均匀采样，有插值方法。