SFR算法详解（三）——浅析关于SFR的一些物理意义及原理

本文内容大量援引自如https://wenku.baidu.com/view/c56f2358804d2b160b4ec04a.html等文章和博客的内容，并加入适当个人补充观点，相关内容仅作技术交流讨论，禁止商用，未经授权产生的相关商业纠纷，本人概不负责，若本文观点存在错误，欢迎联系并指正。

1.什么是点扩散函数？

        点扩散函数（point spread function (PSF) 以下均使用PSF缩写）描述了一个成像系统对一个点光源（物体）的响应。PSF的一般术语就是系统响应，PSF是一个聚焦光学系统的冲击响应。

上图水平面坐标x,y为像面空间坐标，z轴为像面该点处亮度。

以上摘自博客https://blog.csdn.net/qq254612999/article/details/50509793?locationNum=2&fps=1，

https://blog.csdn.net/glorydream2015/article/details/44966369

由上图我们不难发现点扩散函数的切面图具有如下特性：

       在相干光照明时，一个点物在像面上造成的强度分布即为点扩散函数h(xi,yi)。由于像是中心对称的，因此，通过一条过中心的狭缝观察像斑，我们便可以得到如上图b所示的强度分布曲线h(xi),作为沿xi的点扩散函数。

 

2.什么是线扩散函数？

       线扩散函数顾名思义，是由线光源产生，当然你也可以认为是由一排点光源叠加产生。

参见文章http://www.doc88.com/p-6931137520466.html【这篇文章写的很好】，有如下图，小棋认为比较直观的反应了以上内容：

 

       如果用亮狭缝或一个亮线通过光学系统成像，取直线像的长度方向为yi，则沿xi方向的光强分布L(xi)就叫做线扩散函数，线扩散函数是由点扩散函数叠加而成。这里小棋不想深入去研究数学模型，读者可以自行阅读小棋分享的几个链接的资料去了解相关内容，以上介绍的成像模型可参考如下图，由此我们发现，其实点扩散函数和线扩散函数从函数外观上是几乎一样的，但注意它们的物理含义是不同的，点扩散函数实际是描述二维的，而线扩散函数侧重于描述一维。

 

3.什么是边沿扩散函数？

        有了对线扩散函数理解的基础，请试想一下，假设用一与狭缝方向平行的刀片放置在像平面上。开始时刀片完全挡住狭缝像，逐渐沿xi方向移动刀片，刀片后的探测器接收到的全部光信号与刀片位置关系可用下图的阴影变化表示。光通量随xi的变化称边缘扩散函数E(xi)，与线扩散函数关系为：

 

        对照上图，所谓探测器接收到的总光通量就是图中阴影部分的面积，它是对线扩散函数左半区间内的无穷积分【xi的右侧被刀片挡住，光无法通过】。也即是该处的边缘扩散函数值。因此对线扩散函数积分以后我们有类似如下图的边缘扩散函数：

 

        接下来，前述文章对线扩散函数做了进一步展开，试想一下，假如我们有两个一强一弱的线光源，这时候发生在像平面上的响应将是这两个线光源的叠加，如下图：

但如果我们的实景是如下的样子呢？假设我们的面光源是中心区域亮，两侧区域逐渐变暗。

这时候我们可以把面中的每一条竖线当成一个假想的线光源，于是我们可以抽象出如下图所示，无限个线光源产生的线扩散函数叠加的像面情形。

 

        接下来，我们回归正题，探讨一下以上内容对于指导斜棱法的意义，由于小棋在上一篇博文中已经提过了ISO12233中是如何通过4倍超采样提升系统亚像素级的采样能力，这里不再赘述，在上一篇中，通过4倍超采样，实际我们对每一行上的像素到斜边矩心距进行了一个加权平均并归档，而归档的过程实际也是一个相位调整的过程，它将原来每一行间互相偏移的相位，又拉回了同一相位上。因此，我们实际上需要探讨的应该是垂直刀口模型。

同样的抽象方法，我们有如下图的物面和像面响应曲线，我们依旧把物面拆解成无限个线光源的线扩散函数来分析。

       上图实际上反应到光学系统中，就是对线扩散函数积分的一个过程，因此反应到4倍超采样的结果上，就是ESF曲线，为什么这么说，下面，小棋进一步图解一下上述过程。

       如上图，为简化模型，假定线扩散函数间的间距就是一个像素宽度，于是像素点上面包含的所有线扩散函数的面积就是将投影到这个像素点内的光通量的总和，但是它可能来自多个线光源的叠加。

         如上图①区间对应的a像素，仅获取到Ⅰ号线光源的部分能量，而②号区间对应的b号像素，则既获取到了Ⅰ号线光源的部分能量，又获取到了Ⅱ号线光源的部分能量，而我们不难发现b像素获取到的Ⅱ号线光源的能量应该等于a像素获取的对应Ⅰ号线光源的能量，原因是我们假定物面是均匀的，即各线光源具有相同的能量密度。

        于是我们可以把b像素的两部分能量拼接，也即等价于在②区间左边的所有来自Ⅰ号线光源的能量，显然，后面的情况可以依此类推。

         于是投射到各像素上的能量等价于Ⅰ号线扩散函数覆盖区域中，该像素边界处左边所有Ⅰ号线光源能量的积分和，即面积。【两个问题：一、上图最后面有的区间好像缺失了部分能量，那是因为我只画了7条曲线，右边还有无限个，实际情况中也可以再进一步细分出无限个，这里简化了模型。二、像素点到了g以后，就越过了Ⅰ号曲线覆盖范围，显然在分析g像素点时，Ⅰ号曲线的能量已经截止了，这时分析时，应该变更为等价于对应Ⅱ号线扩散函数左边能量的积分和。于是，不难发现，从f开始，后续所有像素均已达到了最大能量密度，即包含整一个线扩散函数的能量密度，f对应Ⅰ号，g对应Ⅱ号，依此类推。】

 

        通过以上分析，请读者再倒回去本文介绍边沿扩散函数时，列举的刀片示例，不难发现，超采样后获得的序列正是边沿扩散函数ESF。

 

       此外，上篇博客中提到，ESF卷积一个有限差分滤波器[-1,1]【即微分过程】后，自然得到的也就是线扩散函数LSF。

 

4.空间频率响应函数SFR的物理意义

        我们进一步来分析一下SFR，以及它的物理意义。这里我先给出结论，结论就是对LSF曲线做一维傅里叶变换，就能得到空间频率响应函数e-SFR。

        我们知道傅里叶变换可以用来分析时域或空间域的幅频响应特性，显然这里我们分析的是空间域不是时域，我们知道在一幅图像中，有的场景光亮变化很频密，可能在很小的空间内，就发生了频繁的明暗变化，有的场景又变化很缓慢，光处于一个比较柔和的变化中，而我们在此处使用傅里叶变换，就是来分析它针对以上不同空间频率下的一个响应情况，假设它在所有频率上响应都没有衰减，那么它将真实反应实际场景在各个频率下的所有细节。但由于真实的光学系统幅频响应曲线【Fourier Translate{LSF}】逐渐衰减，因此,随着空间频率的升高，最后高频细节将由于混叠而丧失。而通过LSF曲线经傅里叶变换后得到的幅频响应曲线SFR，则正是反应了光学系统对于各种空间频率下的响应能力。

        如上图，我们控制物面输入的光能量正弦变化【正弦光栅】，且频率越来越高，幅度不变，则得到以上响应输出，以上这种光能量密度正弦变化的标靶，很具有代表性。

        我们保持输入为一个同频w的正弦波【对每个单波，波峰一半处为原点】，正弦波sinwt的频谱为两个单位冲激响应【其证明利用了欧拉公式，详见链接文章https://wenku.baidu.com/view/0b7e6f0b783e0912a2162ac3.html?pn=51】，结合之前的分析，该正弦输入会与LSF函数卷积并得到输出图像，而空间域的卷积又等于频域的相乘，正弦波的频谱，仅f0和-f0处有值，其余地方均为0。于是在频域输出上，我们得到的输出频谱将仍是两个冲激响应，只是它的幅度受到了SFR曲线的幅度调制，其逆变换以后，则是一个幅度受LSF影响，产生了一定衰减的同频正弦输出，它的幅度衰减情况与SFR的幅频曲线衰减情况正相关【注意我们这里只讨论了幅频响应，相频变化及相移情况此处未做研究】。因此，我们通过不断调整正弦标靶的空间频率，就可以把SFR曲线上的各个点都测算出来。

 

      但由于上述正弦标靶难于实现，工业上我们不得不妥协，而采用了如下方波输入型【矩形光栅】标靶来近似测算。

 

        最后，结合上一篇博文中计算得到的几条结果曲线，解释一下其横轴的表示方法和单位问题。对于ESF和LSF，横轴坐标可以用水平方向跨越实际像素点的个数表示，又或者是像素点个数的4倍，抑或是用pixel size转换为实际的物理距离。

        这里主要讨论一下SFR曲线横轴的几种表现方法，一种是用奈奎斯特频率的倍数，采用这种方式右侧截止位置值为1倍Ny，即1，另一种是采用LP/mm【黑白线对/毫米】的形式【假定pixel size为1..12um,且pixel是正方形的】，有：

                奈奎斯特频率Ny = 1/2 * 采样频率 = 1/2 * [1 / (1.12/1000 mm)] = 446.429 (LP/mm)

        于是，我们可以根据奈奎斯特频率除以SFR曲线样点宽度，得到各细分频率上的SFR值，如SFR宽度为120样点，则SFR曲线横轴坐标步距为：

                步距=采样频率/120= 1000/1.12/120 (LP/mm)= 7.44(LP/mm)

        显然此处的曲线最右侧取样点的频率为 119 * 7.44(LP/mm) = 885.36(LP/mm)，但实际上我们只关心截止频率Ny之前的曲线，即小于446.429(LP/mm)之前的数据。

以上，为什么小棋可以这样计算，其实这是离散傅里叶变换的物理意义，建议以离散傅里叶变换DFT或FFT的物理意义为关键词搜索相关内容，或参看以下相关链接，由于控制篇幅，及网上相关时域分析上的介绍资源较多，这里小棋就略去这部分内容了。

http://blog.sina.com.cn/s/blog_90b4c7ff010157gv.html

https://blog.csdn.net/iloveyoumj/article/details/53308142