图像传感器与信号处理——SFR算法/ISO 12233解读

SFR的全称是Spatial Frequency Response，空间频率响应。介绍SFR算法的参考文献主要是ISO 12233这篇英文文档。和读一般的论文感觉不同，文档读起来更加晦涩难懂，可能是因为太专业，反正第一遍我是没有看懂的，后来参考了一些前辈的资料才搞明白SFR算法到底在干嘛。这篇博客我打算结合ISO 12233的部分内容以及一些参考资料，加入一些我自己的理解，对SFR算法进行一个总结。

1. 前言

以我的理解，ISO 12233这篇文章主要是规定了摄像头分辨率的测量标准，一共有三种方法：基于视觉的分辨率测量方法，基于边界的SFR算法以及基于正弦波的SFR算法。本博客以基于边界的SFR算法为主，考虑到全面性，也简单介绍下基于视觉的分辨率测量方法以及基于正弦波的SFR算法。

（1）那么首先我们得知道什么是分辨率？
在生活中，我们经常将分辨率定义为相机传感器像素的数量，而在ISO 12233中，将分辨率定义为相机捕捉空间细节的能力，可以分为单值标准和多值标准，单值标准即通过基于视觉的分辨率测量方法测得，而多值标准即通过两种SFR算法获得。

（2）单值标准好理解，那么通过SFR算法获得的多值标准是如何描述分辨率的呢？
这个通过SFR算法获得的多值标准其实指的就是MTF，全称是Modulation Transfer Function，调制传递函数。该函数是空间频率的传输函数，典型的调制传递函数的示意图如下图所示：

如图所示，典型的调制传递函数的示意图都会表现为一条向右衰减的曲线，这是因为空间频率越高（物体图案中的细节越多〉，镜头对物体的再现能力越弱。可以认为高频率下的高调制传递函数值指的是分辨率，而中等频率下的调制传递函数值指的是对比度，下图展示了该结论：

（3）最后，既然是测量相机的分辨率，那么具体是如何测量的呢？
在ISO 12233中给出了相机分辨率测量的条件，然后规定了一些列对于光照、白平衡、聚焦的设置要求，这里就不进行赘述，反正文档里都有。当搭建的环境满足以上要求后，将Test Chart放在规定位置，用待测的相机拍摄一张图片，然后用基于视觉的分辨率测量方法或者SFR算法对图片进行分析，最终获得该相机的分辨率。当然，对于不同的算法所用的Test Chart也不相同，下面就开始对各个算法进行介绍。

2. 基于视觉的分辨率测量方法

如下图所示是基于视觉的分辨率测量方法所用的Test Chart：

基于视觉的分辨率测量方法是利用图像空间频率过高而产生的混叠效应来进行判定的，如下图所示待测相机在前文所述的测试环境下拍摄的Test Chart图像的一部分，可以看到在图像中箭头所指的位置，最左侧两条线段出现了混叠效应，根据混叠效应产生位置，通过查表就可以获得相机分辨率：

基于视觉的分辨率测量方法我觉得需要注意如下两点：
（1）这种测试方法获得的结果单位是LW/PH，即line widths per picture height；
（2）文中提到，相机拍摄的Test Chart可能并不是占据整张像，因此获得的结果也需要按照Test Chart的大小比例进行放缩，例如，图像中的Test Chart只占据整张图像的1/2，按照Test Chart获得的结果是2000LW/PH，那么相机的实际分辨率是4000LW/P。.

3. 基于边界的SFR算法

如下图所示是基于边界的SFR算法所用的Test Chart（ISO 12233 2017版）：

基于边界的SFR算法的介绍可以参考博客SFR算法详解（二）——斜棱法，这里结合前辈的博客和ISO 12233说一下我的理解：
算法的流程图如下图所示：

具体如下：
（1）选取ROI区域
ROI区域的选取如下图红色边框区域所示，选择的ROI区域由黑白两部分组成，中间为分界线。分界线与的ROI区域的上下边框交接点距离左右边框不得小于两个像素（其实就是所选取的ROI区域中的分界线不能太过倾斜），此外，ROI区域亮度调制比不得小于20%（其实就是黑白对比明显）。

（2）进行数据归一化
数据归一化处理通过逆光电转换函数(OECF)/逆Gamma校正对非线性的图像编码值(若为888位图像，即灰度值0−2550-2550−255)进行逆向线性化处理，将其还原为一个类CCD原始光电采样信号的序列。关于Gamma校正的原理可以参考色彩中的Gamma值是什么，在ISO 12233中这一步表示为：φ(p,r)=aOECF⁡[DNred]+bOECF⁡[DNgreen]+cOECF⁡[DNblue]\varphi(p, r)=a \operatorname{OECF}\left[DN_{\mathrm{red}}\right]+b \operatorname{OECF}\left[DN_{\mathrm{green}}\right]+c \operatorname{OECF}\left[DN_{\mathrm{blue}}\right] φ(p,r)=aOECF[DNred]+bOECF[DNgreen]+cOECF[DNblue]其中DNredDN_{\mathrm{red}}DNred，DNgreenDN_{\mathrm{green}}DNgreen，DNblueDN_{\mathrm{blue}}DNblue分别为各个通道的数字输出，φ(p,r)\varphi(p, r)φ(p,r)为CCD上第rrr行，第ppp列原始光电采样信号，aaa，bbb，ccc为各个通道的加权值

（3）计算每一行的矩心
在ISO 12233中通过如下公式进行计算：C(r)=∑p=1P−1p[ϕ(p+1,r)−ϕ(p,r)]∑p=1P−1[ϕ(p+1,r)−ϕ(p,r)]−0.5C(r)=\frac{\sum_{p=1}^{P-1} p[\phi(p+1, r)-\phi(p, r)]}{\sum_{p=1}^{P-1}[\phi(p+1, r)-\phi(p, r)]}-0.5 C(r)=∑p=1P−1[ϕ(p+1,r)−ϕ(p,r)]∑p=1P−1p[ϕ(p+1,r)−ϕ(p,r)]−0.5其中，PPP为ROI区域的总列数，C(r)C(r)C(r)为第rrr行的矩心位置，这里我们以SFR算法详解（二）——斜棱法中的例子对公式进行说明，我们假定ROI区域光电采集信号如下表所示（为了方便解释，这里仍然用图像编码值进行举例说明）：01215832362552540031397244253254211181102412522540031812424925325211221140251252255217261602532552543113417125225225412745189249254255\begin{array}{|c|c|c|c|c|c|c|c|} \hline 0 & 1 & 2 & 15 & 83 & 236 & 255 & 254 \\ \hline 0 & 0 & 3 & 13 & 97 & 244 & 253 & 254 \\ \hline 2 & 1 & 1 & 18 & 110 & 241 & 252 & 254 \\ \hline 0 & 0 & 3 & 18 & 124 & 249 & 253 & 252 \\ \hline 1 & 1 & 2 & 21 & 140 & 251 & 252 & 255 \\ \hline 2 & 1 & 7 & 26 & 160 & 253 & 255 & 254 \\ \hline 3 & 1 & 1 & 34 & 171 & 252 & 252 & 254 \\ \hline 1 & 2 & 7 & 45 & 189 & 249 & 254 & 255 \\ \hline \end{array} 00201231101011122313271715131818212634458397110124140160171189236244241249251253252249255253252253252255252254254254254252255254254255按照公式，首先用每一行相邻两列相减，得到下表：11136815319−103108414791−10179213111203151061254−1011911911113−1619134932−1−2033137810215381446051\begin{array}{|c|c|c|c|c|c|c|} \hline 1 & 1 & 13 & 68 & 153 & 19 & -1 \\ \hline 0 & 3 & 10 & 84 & 147 & 9 & 1 \\ \hline-1 & 0 & 17 & 92 & 131 & 11 & 2 \\ \hline 0 & 3 & 15 & 106 & 125 & 4 & -1 \\ \hline 0 & 1 & 19 & 119 & 111 & 1 & 3 \\ \hline-1 & 6 & 19 & 134 & 93 & 2 & -1 \\ \hline-2 & 0 & 33 & 137 & 81 & 0 & 2 \\ \hline 1 & 5 & 38 & 144 & 60 & 5 & 1 \\ \hline \end{array} 10−100−1−211303160513101715191933386884921061191341371441531471311251119381601991141205−112−13−121然后对每一行求加权平均获得矩心，对于第一行来说即C(1)=1∗1+2∗1+3∗13+4∗68+5∗153+6∗19+7∗−11+1+13+68+153+19−1=4.669291C(1)=\frac{1*1+2*1+3*13+4*68+5*153+6*19+7*-1}{1+1+13+68+153+19-1}=4.669291 C(1)=1+1+13+68+153+19−11∗1+2∗1+3∗13+4∗68+5∗153+6∗19+7∗−1=4.669291对每一行都进行同样的操作及获得每一行矩心的位置：4.6692914.5984254.5753974.432544.3976384.2619054.2390444.086614\begin{array}{|r|} \hline 4.669291 \\ \hline 4.598425 \\ \hline 4.575397 \\ \hline 4.43254 \\ \hline 4.397638 \\ \hline 4.261905 \\ \hline 4.239044 \\ \hline 4.086614 \\ \hline \end{array} 4.6692914.5984254.5753974.432544.3976384.2619054.2390444.086614

（4）对矩心进行直线拟合
在ISO 12233中公式如下：m=[ΔrΔC(r)]‾m=\overline{\left[\frac{\Delta r}{\Delta C(r)}\right]} m=[ΔC(r)Δr]这个公式的意思就是求相邻两行之间的矩心的斜率，然后对所有斜率求均值mmm，最后获得函数：S(r)=R/2−rmS(r)=\frac{R / 2-r}{m} S(r)=mR/2−r其中RRR是总ROI区域的总行数，S(r)S(r)S(r)并不是拟合的直线，但是在后面和其他公式结合起来之后起到了相当于矩心直线的作用。

（5）获得边缘扩散函数ESF
这一步相对来说比较复杂，在ISO 12233中公式如下：ESF′(j)=∑r=1R∑p=1Pϕ(p,r)⋅α(p,r,j)∑r=1R∑p=1Pα(p,r,j)E S F^{\prime}(j)=\frac{\sum_{r=1}^{R} \sum_{p=1}^{P} \phi(p, r) \cdot \alpha(p, r, j)}{\sum_{r=1}^{R} \sum_{p=1}^{P} \alpha(p, r, j)} ESF′(j)=∑r=1R∑p=1Pα(p,r,j)∑r=1R∑p=1Pϕ(p,r)⋅α(p,r,j)其中函数α(p,r,j)\alpha(p, r, j)α(p,r,j)作用是将上述ROI光电采集信号归并到边缘扩散函数ESF中，其定义如下α(p,r,j)={1,−0,125≤[p−S(r)−j]<0,1250,otherwise \alpha(p, r, j)=\left\{\begin{array}{ll} 1, & -0,125 \leq[p-S(r)-j]<0,125 \\ 0, & \text { otherwise } \end{array}\right. α(p,r,j)={1,0,−0,125≤[p−S(r)−j]<0,125 otherwise 其中jjj为整数，其实也就是边缘扩散函数的下标，这个公式看起来其实不是很好理解，下面通过继续演示前文的例子来说明这个过程。上面已经可以得到每行矩心的位置，那么记下来计算各列到矩心的距离，例如第一行第一列有：1−4.669291=−3.6692911-4.669291=-3.669291 1−4.669291=−3.669291则所有距离如下表所示：−3.66929−2.66929−1.66929−0.669290.3307091.3307092.3307093.330709−3.59843−2.59843−1.59843−0.598430.4015751.4015752.4015753.401575−3.5754−2.5754−1.5754−0.57540.4246031.4246032.4246033.424603−3.43254−2.43254−1.43254−0.432540.567461.567462.567463.56746−3.39764−2.39764−1.39764−0.397640.6023621.6023622.6023623.602362−3.2619−2.2619−1.2619−0.26190.7380951.7380952.7380953.738095−3.23904−2.23904−1.23904−0.239040.7609561.7609562.7609563.760956−3.08661−2.08661−1.08661−0.086610.9133861.9133862.9133863.913386\begin{array}{|r|r|r|r|r|r|r|} \hline-3.66929 & -2.66929 & -1.66929 & -0.66929 & 0.330709 & 1.330709 & 2.330709 & 3.330709 \\ \hline-3.59843 & -2.59843 & -1.59843 & -0.59843 & 0.401575 & 1.401575 & 2.401575 & 3.401575 \\ \hline-3.5754 & -2.5754 & -1.5754 & -0.5754 & 0.424603 & 1.424603 & 2.424603 & 3.424603 \\ \hline-3.43254 & -2.43254 & -1.43254 & -0.43254 & 0.56746 & 1.56746 & 2.56746 & 3.56746 \\ \hline-3.39764 & -2.39764 & -1.39764 & -0.39764 & 0.602362 & 1.602362 & 2.602362 & 3.602362 \\ \hline-3.2619 & -2.2619 & -1.2619 & -0.2619 & 0.738095 & 1.738095 & 2.738095 & 3.738095 \\ \hline-3.23904 & -2.23904 & -1.23904 & -0.23904 & 0.760956 & 1.760956 & 2.760956 & 3.760956 \\ \hline-3.08661 & -2.08661 & -1.08661 & -0.08661 & 0.913386 & 1.913386 & 2.913386 & 3.913386 \\ \hline \end{array} −3.66929−3.59843−3.5754−3.43254−3.39764−3.2619−3.23904−3.08661−2.66929−2.59843−2.5754−2.43254−2.39764−2.2619−2.23904−2.08661−1.66929−1.59843−1.5754−1.43254−1.39764−1.2619−1.23904−1.08661−0.66929−0.59843−0.5754−0.43254−0.39764−0.2619−0.23904−0.086610.3307090.4015750.4246030.567460.6023620.7380950.7609560.9133861.3307091.4015751.4246031.567461.6023621.7380951.7609561.9133862.3307092.4015752.4246032.567462.6023622.7380952.7609562.9133863.3307093.4015753.4246033.567463.6023623.7380953.7609563.913386然后以[−0.125，0.125][-0.125，0.125][−0.125，0.125]为第000区间，以0.250.250.25为区间长度朝正负方向进行延伸，对以上距离进行区间划分，例如[0.125，0.375][0.125，0.375][0.125，0.375]为第111区间，[−0.375，−0.125][-0.375，-0.125][−0.375，−0.125]为第−1-1−1区间，由此可获得下表：−15−11−7−315913−14−10−6−2261014−14−10−6−2261014−14−10−6−2261014−14−10−6−2261014−13−9−5−1371115−13−9−5−1371115−12−8−40481216\begin{array}{|c|c|c|c|c|c|c|c|} \hline-15 & -11 & -7 & -3 & 1 & 5 & 9 & 13 \\ \hline-14 & -10 & -6 & -2 & 2 & 6 & 10 & 14 \\ \hline-14 & -10 & -6 & -2 & 2 & 6 & 10 & 14 \\ \hline-14 & -10 & -6 & -2 & 2 & 6 & 10 & 14 \\ \hline-14 & -10 & -6 & -2 & 2 & 6 & 10 & 14 \\ \hline-13 & -9 & -5 & -1 & 3 & 7 & 11 & 15 \\ \hline-13 & -9 & -5 & -1 & 3 & 7 & 11 & 15 \\ \hline-12 & -8 & -4 & 0 & 4 & 8 & 12 & 16 \\ \hline \end{array} −15−14−14−14−14−13−13−12−11−10−10−10−10−9−9−8−7−6−6−6−6−5−5−4−3−2−2−2−2−1−1012222334566667789101010101111121314141414151516上表中的区间值其实就是公式中定义的参数jjj，最后一步就是按照上表中的区间值对位于同一区间的光电信号进行平均，以第−15-15−15区间为例，即0/1=10/1=10/1=1，以第−14-14−14区间为例，即(0+2+0+1)/4=0.75(0+2+0+1)/4=0.75(0+2+0+1)/4=0.75，那么将获得边缘扩散函数ESF如下表所示：x−15−14−13−12−11−10−9−8y00.750210.512x−7−6−5−4−3−2−10y22.25471517.53045x12345678y83117.75165.5189236246.25252.5249x910111213141516y255252.5253.5254254253.75254255\begin{array}{|c|c|c|c|c|c|c|c|c|} \hline \mathrm{x} & -15 & -14 & -13 & -12 & -11 & -10 & -9 & -8 \\ \hline \mathrm{y} & 0 & 0.75 & 0 & 2 & 1 & 0.5 & 1 & 2 \\ \hline \mathrm{x} & -7 & -6 & -5 & -4 & -3 & -2 & -1 & 0 \\ \hline \mathrm{y} & 2 & 2.25 & 4 & 7 & 15 & 17.5 & 30 & 45 \\ \hline \mathrm{x} & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ \hline \mathrm{y} & 83 & 117.75 & 165.5 & 189 & 236 & 246.25 & 252.5 & 249 \\ \hline \mathrm{x} & 9 & 10 & 11 & 12 & 13 & 14 & 15 & 16 \\ \hline \mathrm{y} & 255 & 252.5 & 253.5 & 254 & 254 & 253.75 & 254 & 255 \\ \hline \end{array} xyxyxyxy−150−721839255−140.75−62.252117.7510252.5−130−543165.511253.5−122−47418912254−111−315523613254−100.5−217.56246.2514253.75−91−1307252.515254−82045824916255进而可以获得边缘扩散函数ESF的曲线图（关于这个下标jjj其实我有点纠结，按照公式推导横坐标确实会有负值，但是下文中LSF又是从222开始）

（7）四倍超采样
在获得边缘扩散函数的过程中采用的离散距离单位缩小为0.250.250.25个像素，我们将边缘扩散函数的横坐标乘以444倍，然后将没有函数值的横坐标位置，向前寻找非零的函数值进行替换。这样将离散距离单位归为111，但它代表的还0.250.250.25个像素，而这也就是我们说的四倍倍超采样。

（8）生成线性扩散序列LSF
在ISO 12233中公式如下：LSF⁡W′(j)=W(j)ESF⁡′(j+1)−ESF⁡′(j−1)2,for j=2,…,N−1\operatorname{LSF}_{W}^{\prime}(j)=W(j) \frac{\operatorname{ESF}^{\prime}(j+1)-\operatorname{ESF}^{\prime}(j-1)}{2}, \text { for } j=2, \dots, \mathrm{N}-1 LSFW′(j)=W(j)2ESF′(j+1)−ESF′(j−1), for j=2,…,N−1其中W(j)=0,54+0,46cos⁡[2π(j−2X)/4X]W(j)=0,54+0,46 \cos [2 \pi(j-2 X) / 4 X] W(j)=0,54+0,46cos[2π(j−2X)/4X]其实线性扩散序列LSF就是边缘扩散函数ESF均值滤波后再进行微分，在此基础上通过汉明窗滤波去除噪声，线性扩散序列LSF的图像大致如下图所示：
上面是根据前面的例子计算得到的线性扩散函数LSF，而理想的线性扩散函数LSF应该如下图所示
（9）归一化离散傅里叶变换
在ISO 12233中公式如下：e−SFR(k)=D(k)∣∑j=1NLSFW′(j)e−i2πkj/N∑j=1NLSFW′(j)∣,for k=0,1,2,…,N/2,or (N+1)/2if Nis odd e-S F R(k)=D(k)\left|\frac{\sum_{j=1}^{N} L S F_{W}^{\prime}(j) e^{-i 2 \pi k j / N}}{\sum_{j=1}^{N} L S F_{W}^{\prime}(j)}\right|, \text { for } k=0,1,2, \ldots, N / 2, \text { or }(N+1) / 2 \text { if } N \text { is odd } e−SFR(k)=D(k)∣∣∣∣∣∑j=1NLSFW′(j)∑j=1NLSFW′(j)e−i2πkj/N∣∣∣∣∣, for k=0,1,2,…,N/2, or (N+1)/2 if N is odd 即对上述线性扩散序列应用归一化离散傅里叶变换即可过得最后的结果，大致如下图所示：

4. 基于正弦波的SFR算法

如下图所示是基于正弦波的SFR算法所用的Test Chart：

但是关于这一部分算法，网上能够参考的资料较少，我反复看了ISO 12233文档好几遍，始终没能理解算法是怎么实现的，因为之后不一定会用到这一部分知识，所以不打算再话更多时间来研究这个算法，如果有同学刚好看到了这篇博客，和我有相同的问题，希望可以相互交流下~当然，如果有大佬能够直接指点下就最好了！

这篇博客先写到这里了，欢迎指正交流。