【基于python的图像质量测试】ISO12233 eSFR MTF (二)

接前文
#include 【基于python的图像质量测试】ISO12233 eSFR MTF (一)

3.3，综合每一行图像计算拟合边缘位置的斜率和截距(续)

得到拟合直线后第二次滤波
第二次滤波时将以每行的黑白边缘位置作为汉明窗的中点对求导后的每一行重新进行滤波
这个分界点即为上一轮的拟合直线所确认的点。
注意：这里并不是在上一次滤波后的数据上进行，而是重新对求导后的每一行进行滤波

以第一行为例：图示蓝点视为第一行中心点

下图为对第一行图像的导数以及汉明窗的应用

第一行滤波后的图像如下图所示：

之后就是重复之前的操作：计算矩心+拟合直线
代码示例：
声明x1,y1用于存放第二次滤波后得到的点

x1 = []
y1 = []
for k in range(h):
    centroid = int((k-intercept)/slope)#矩心
    diff_row = np.convolve(img_lum[k].astype(int),kernel,mode="valid")#卷积自定义的微分
    if centroid <= len(diff_row)/2:#如果矩心位于左侧，则在汉明窗右侧进行补全
        hamming_window_2 = np.hamming(2*centroid)#汉明窗大小
        completion_list = np.array([hamming_window_2[0] for _ in range(len(diff_row)-len(hamming_window_2))])#使用汉明窗单侧的最小值进行补全
        hamming_window_2_fixed = np.append(hamming_window_2,completion_list)#把自定义的滤波补全到整行大小
        diff_row = diff_row*hamming_window_2_fixed#对导数应用滤波
    else:#如果矩心位于右侧，则在汉明窗左侧进行补全
        hamming_window_2 = np.hamming((w-centroid)*2)
        completion_list = np.array([hamming_window_2[0] for _ in range(len(diff_row)-len(hamming_window_2))])#生成一个列表用于补全汉明窗
        hamming_window_2_fixed = np.append(completion_list,hamming_window_2)
        deff_row = diff_row*hamming_window_2_fixed


    #以下代码和第一次滤波一样
    dt1 = np.sum(diff_row)
    dt = result = np.sum(diff_row * (np.arange(len(diff_row))+1))#每个元素*（元素索引+1），因为索引从0开始所以需要+1
    shift = dt/dt1
    x1.append(int(shift))
    y1.append(k)
coefficients = np.polyfit(x1, y1, 1)#线的拟合
slope = coefficients[0]#斜率
intercept = coefficients[1]#截距

最后得到的直线：

看起来好像没啥区别

4、边缘扩展函数ESF

4.1、4倍超采样

关于超采样，贴一下关于ISO:12233的原文：

贴公式：

目前似乎对 $\alpha （p,r,j）$的处理存在一些分歧，我在查阅一些资料后尝试采取以下方法：
通过excel举一个例子
以下近似看作原始的图像数据

扩增4倍后

尝试将没有函数值的横坐标位置，向前寻找非零的函数值进行替换

然后将以拟合直线为轴进行归档

通过拟合直线计算出与每一行像素的交点，这个焦点坐标x4并向下取整（因为是扩增4倍的）。

for j in range(h):
	centroid_4 = int(((j-intercept)/slope)*4)

4.2、计算ESF

以拟合直线对应的坐标为轴，向左向右分别以一个像素为单位计算每一列的平均值，如上图
以此将二维图像计算为一维的ESF

#创建两个足够大的容器用于对图像进行归档
input_arr_n = [[] for _ in range(400)]
input_arr_p = [[] for _ in range(400)] 
#4倍超采样扩增数组
img_lum_4 = np.repeat(img_lum,4,axis=1)
#获得4倍超采样后的宽和高
h,w = img_lum_4.shape

for j in range(h):
    #4倍超采样的边缘位置，向下取整
    centroid_4 = int(((j-intercept)/slope)*4)
    #以边缘为0点开始归档,分为左侧和右侧，0的左侧存入input_arr_n,右侧存入input_arr_p
    for left_cnt in range(0,int(centroid_4),1):
        input_arr_n[left_cnt].append(img_lum_4[j][centroid_4-left_cnt])
        
    for right_cnt in range(centroid_4 + 1,w,1):
        input_arr_p[right_cnt-centroid_4+1].append(img_lum_4[j][right_cnt])
        
input_arr_n = [sublist for sublist in input_arr_n if sublist]#清空空列表
input_arr_p = [sublist for sublist in input_arr_p if sublist]#清空空列表
averages_n = [sum(sublist)/len(sublist) for sublist in input_arr_n]#计算归档后的平均值
averages_p = [sum(sublist)/len(sublist) for sublist in input_arr_p]
#负向列表里的数据顺序进行反转
averages_n.reverse()
ESF = averages_n + averages_p#边缘扩散函数
plt.plot([i for i in range(len(ESF))],ESF)
plt.title("ESF")
plt.show()

*众所周知，python并不适合直接处理大数据量的迭代。通常会借用numpy等python包协助处理。在ROI较大时，上述代码的运行将会变得非常缓慢。因此这段代码还有很大的优化空间。
同时我还有个不成熟的想法，就是把上面这段代码使用C++编写并编译成dll文件，在使用python的ctypes库去调用这个dll。 *

ESF效果图：

5、线扩展函数LSF

$$LSF(J)=\frac{ESF(J+1)-ESF(J-1)}{2},forj=2,...,N-1$$
使用[-1/2,0,1/2]滤波器进行求导
这里再次尝试使用np.convolve()来达成求导效果

#反转滤波器以达到求导效果
ESF_kernel = np.array([0.5,0,-0.5])
#LSF = np.diff(ESF)#线扩散函数
LSF = np.convolve(ESF,ESF_kernel,mode="valid")

之后便是LSF对横轴居中。资料上没说是如何居中的，我这边做法是将LSF对最高点非对称的部分截断后居中
加汉明窗滤波，这个汉明依旧是要以居中后的LSF曲线的最高点为中心

#应用汉明窗
#按前面的矩心计算公式求出最高点
dt1 = np.sum(LSF)
dt = result = np.sum(LSF * (np.arange(len(LSF))+1))
#汉明窗中点坐标
shift = dt/dt1
#将LSF居中
LSF = LSF[len(LSF)-2*(len(LSF)-int(shift)):]#截断，取以最高点对称的曲线
hamming_window = np.hamming(len(LSF))
LSF = LSF*hamming_window
plt.plot([i for i in range(len(LSF))],LSF)
plt.title("LSF")
plt.show()

LSF曲线以及对应的汉明窗

6、傅里叶变换DFT

贴公式：

将公式分解处理
其中$D(k)$是矫正函数

LSF'w是从图表图像的选定区域形成的加窗、平均、居中、超采样线扩展函数,可以理解为这里的LSF'w与前面求得的LSF等价
那么，公式上半部分即为$f1=|{\sum }_{J=1}^{N}LSF(J)e^{-i2\pi kj/N}|$
对LSF曲线取绝对值并做傅里叶变换

SFR = np.fft.fft(LSF)
abs_SFR = np.abs(SFR)

公式下半部分$f2=|{\sum }_{j=1}^{N}LSF(j)|$的作用是归一化
那么$f1/f2$即为

freq = np.fft.fftfreq(len(SFR),0.01)
normalize_SFR = abs_SFR/abs(sum(LSF))
plt.plot(freq[0:int(len(freq)/2)],normalize_SFR[0:int(len(normalize_SFR)/2)])#只绘制频率为正的部分
plt.title("e-SFR(not correct)")
plt.xlabel("X")
plt.ylabel("Y")
plt.show()

得到一个未矫正的图像：

这是一个跟频率相关的FFT归一化曲线。
尝试处理矫正函数$D(k)$：

Dk = [min([1/math.sin(2*x*math.pi/len(LSF)),10]) for x in range(1,int(len(LSF)/2)+1,1)]

矫正后的图像可视为
e-sfr = e-sfr(not correct)*Dk
获得矫正后的图像：

由于资料有限，这篇文章中可能存在一些错误，大致罗列一下我目前不是很确定的点：
1，4倍超采样：查阅一些文章后发现方法并不是很统一，我这边选择了一个比较好实现的方法
2，关于求LSF的汉明窗的处理，我这边选择截断可能导致最后做DFT的时候数据出现偏差
3，上一篇的文章的反gamma变换，那个公式是我自己推导的，我并没有找到参考公式
4，最后的e-SFR曲线。我在横轴上标明了???是因为真正的横轴可能和频率并不相关，也没有找到绘制这个图像的具体方法

希望各位看官不吝赐教

参考文章：图像分辨率测试ISO12233 - 2017中文翻译