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

图像传感器与信号处理——SFR算法/ISO 12233解读
图像传感器与信号处理——SFR算法/ISO 12233解读
1. 前言
2. 基于视觉的分辨率测量方法
3. 基于边界的SFR算法
4. 基于正弦波的SFR算法
图像传感器与信号处理——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校正对非线性的图像编码值(若为8 88位图像，即灰度值0 − 255 0-2550−255)进行逆向线性化处理，将其还原为一个类CCD原始光电采样信号的序列。关于Gamma校正的原理可以参考色彩中的Gamma值是什么，在ISO 12233中这一步表示为：
φ ( p , r ) = a OECF ⁡ [ D N r e d ] + b OECF ⁡ [ D N g r e e n ] + c OECF ⁡ [ D N b l u e ] \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[DN
red

]+bOECF[DN
green

]+cOECF[DN
blue

]
其中D N r e d DN_{\mathrm{red}}DN
red

，D N g r e e n DN_{\mathrm{green}}DN
green

，D N b l u e DN_{\mathrm{blue}}DN
blue

分别为各个通道的数字输出，φ ( p , r ) \varphi(p, r)φ(p,r)为CCD上第r rr行，第p pp列原始光电采样信号，a aa，b bb，c cc为各个通道的加权值

（3）计算每一行的矩心
在ISO 12233中通过如下公式进行计算：
C ( r ) = ∑ p = 1 P − 1 p [ ϕ ( p + 1 , r ) − ϕ ( p , r ) ] ∑ p = 1 P − 1 [ ϕ ( p + 1 , r ) − ϕ ( p , r ) ] − 0.5 C(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=1
P−1

[ϕ(p+1,r)−ϕ(p,r)]
∑
p=1
P−1

p[ϕ(p+1,r)−ϕ(p,r)]

−0.5
其中，P PP为ROI区域的总列数，C ( r ) C(r)C(r)为第r rr行的矩心位置，这里我们以SFR算法详解（二）——斜棱法中的例子对公式进行说明，我们假定ROI区域光电采集信号如下表所示（为了方便解释，这里仍然用图像编码值进行举例说明）：
0 1 2 15 83 236 255 254 0 0 3 13 97 244 253 254 2 1 1 18 110 241 252 254 0 0 3 18 124 249 253 252 1 1 2 21 140 251 252 255 2 1 7 26 160 253 255 254 3 1 1 34 171 252 252 254 1 2 7 45 189 249 254 255
00201231101011122313271715131818212634458397110124140160171189236244241249251253252249255253252253252255252254254254254252255254254255

然后对每一行求加权平均获得矩心，对于第一行来说即
C ( 1 ) = 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)=\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−1
1∗1+2∗1+3∗13+4∗68+5∗153+6∗19+7∗−1

=4.669291
对每一行都进行同样的操作及获得每一行矩心的位置：
4.669291 4.598425 4.575397 4.43254 4.397638 4.261905 4.239044 4.086614
4.6692914.5984254.5753974.432544.3976384.2619054.2390444.086614
4.669291
4.598425
4.575397
4.43254
4.397638
4.261905
4.239044
4.086614
4.669291
4.598425
4.575397
4.43254
4.397638
4.261905
4.239044
4.086614

（4）对矩心进行直线拟合
在ISO 12233中公式如下：
m = [ Δ r Δ C ( r ) ] ‾ m=\overline{\left[\frac{\Delta r}{\Delta C(r)}\right]}
m=
[
ΔC(r)
Δr

]

这个公式的意思就是求相邻两行之间的矩心的斜率，然后对所有斜率求均值m mm，最后获得函数：
S ( r ) = R / 2 − r m S(r)=\frac{R / 2-r}{m}
S(r)=
m
R/2−r

其中R RR是总ROI区域的总行数，S ( r ) S(r)S(r)并不是拟合的直线，但是在后面和其他公式结合起来之后起到了相当于矩心直线的作用。

（5）获得边缘扩散函数ESF
这一步相对来说比较复杂，在ISO 12233中公式如下：
E S F ′ ( j ) = ∑ r = 1 R ∑ p = 1 P ϕ ( p , r ) ⋅ α ( p , r , j ) ∑ r = 1 R ∑ p = 1 P α ( 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=1
R

∑
p=1
P

α(p,r,j)
∑
r=1
R

∑
p=1
P

ϕ(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 , 125 0 , otherwise \alpha(p, r, j)=\left\{
1,0,−0,125≤[p−S(r)−j]<0,125 otherwise
1
,
−
0
,
125
≤
[
�
−
�
(
�
)
−
�
]
<
0
,
125
0
,
otherwise
\right.
α(p,r,j)={
1,
0,


−0,125≤[p−S(r)−j]<0,125
otherwise

其中j jj为整数，其实也就是边缘扩散函数的下标，这个公式看起来其实不是很好理解，下面通过继续演示前文的例子来说明这个过程。上面已经可以得到每行矩心的位置，那么记下来计算各列到矩心的距离，例如第一行第一列有：
1 − 4.669291 = − 3.669291 1-4.669291=-3.669291
1−4.669291=−3.669291
则所有距离如下表所示：
− 3.66929 − 2.66929 − 1.66929 − 0.66929 0.330709 1.330709 2.330709 3.330709 − 3.59843 − 2.59843 − 1.59843 − 0.59843 0.401575 1.401575 2.401575 3.401575 − 3.5754 − 2.5754 − 1.5754 − 0.5754 0.424603 1.424603 2.424603 3.424603 − 3.43254 − 2.43254 − 1.43254 − 0.43254 0.56746 1.56746 2.56746 3.56746 − 3.39764 − 2.39764 − 1.39764 − 0.39764 0.602362 1.602362 2.602362 3.602362 − 3.2619 − 2.2619 − 1.2619 − 0.2619 0.738095 1.738095 2.738095 3.738095 − 3.23904 − 2.23904 − 1.23904 − 0.23904 0.760956 1.760956 2.760956 3.760956 − 3.08661 − 2.08661 − 1.08661 − 0.08661 0.913386 1.913386 2.913386 3.913386
−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
−
3.66929
−
2.66929
−
1.66929
−
0.66929
0.330709
1.330709
2.330709
3.330709
−
3.59843
−
2.59843
−
1.59843
−
0.59843
0.401575
1.401575
2.401575
3.401575
−
3.5754
−
2.5754
−
1.5754
−
0.5754
0.424603
1.424603
2.424603
3.424603
−
3.43254
−
2.43254
−
1.43254
−
0.43254
0.56746
1.56746
2.56746
3.56746
−
3.39764
−
2.39764
−
1.39764
−
0.39764
0.602362
1.602362
2.602362
3.602362
−
3.2619
−
2.2619
−
1.2619
−
0.2619
0.738095
1.738095
2.738095
3.738095
−
3.23904
−
2.23904
−
1.23904
−
0.23904
0.760956
1.760956
2.760956
3.760956
−
3.08661
−
2.08661
−
1.08661
−
0.08661
0.913386
1.913386
2.913386
3.913386
−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.08661


0.330709
0.401575
0.424603
0.56746
0.602362
0.738095
0.760956
0.913386


1.330709
1.401575
1.424603
1.56746
1.602362
1.738095
1.760956
1.913386


2.330709
2.401575
2.424603
2.56746
2.602362
2.738095
2.760956
2.913386


3.330709
3.401575
3.424603
3.56746
3.602362
3.738095
3.760956
3.913386

然后以[ − 0.125 ， 0.125 ] [-0.125，0.125][−0.125，0.125]为第0 00区间，以0.25 0.250.25为区间长度朝正负方向进行延伸，对以上距离进行区间划分，例如[ 0.125 ， 0.375 ] [0.125，0.375][0.125，0.375]为第1 11区间，[ − 0.375 ， − 0.125 ] [-0.375，-0.125][−0.375，−0.125]为第− 1 -1−1区间，由此可获得下表：
− 15 − 11 − 7 − 3 1 5 9 13 − 14 − 10 − 6 − 2 2 6 10 14 − 14 − 10 − 6 − 2 2 6 10 14 − 14 − 10 − 6 − 2 2 6 10 14 − 14 − 10 − 6 − 2 2 6 10 14 − 13 − 9 − 5 − 1 3 7 11 15 − 13 − 9 − 5 − 1 3 7 11 15 − 12 − 8 − 4 0 4 8 12 16
−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

进而可以获得边缘扩散函数ESF的曲线图（关于这个下标j jj其实我有点纠结，按照公式推导横坐标确实会有负值，但是下文中LSF又是从2 22开始）

（7）四倍超采样
在获得边缘扩散函数的过程中采用的离散距离单位缩小为0.25 0.250.25个像素，我们将边缘扩散函数的横坐标乘以4 44倍，然后将没有函数值的横坐标位置，向前寻找非零的函数值进行替换。这样将离散距离单位归为1 11，但它代表的还0.25 0.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
LSF
W
′

(j)=W(j)
2
ESF
′
(j+1)−ESF
′
(j−1)

, for j=2,…,N−1
其中
W ( j ) = 0 , 54 + 0 , 46 cos ⁡ [ 2 π ( j − 2 X ) / 4 X ] 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 − S F R ( k ) = D ( k ) ∣ ∑ j = 1 N L S F W ′ ( j ) e − i 2 π k j / N ∑ j = 1 N L S F W ′ ( j ) ∣ , for k = 0 , 1 , 2 , … , N / 2 , or ( N + 1 ) / 2 if N is 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=1
N

LSF
W
′

(j)
∑
j=1
N

LSF
W
′

(j)e
−i2πkj/N

∣
∣
∣
∣
∣

, for k=0,1,2,…,N/2, or (N+1)/2 if N is odd
即对上述线性扩散序列应用归一化离散傅里叶变换即可过得最后的结果，大致如下图所示：