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

标准步骤

操作步骤

1，选择ROI
2，图像线性化转化为有效曝光（gamma反变换，将相机光电转换功能（OECF）转换为焦平面曝光或场景亮度的数字信号转换）
3，边缘位置估计

• 3.1，计算每一行的一维导数(每行的ESF)
• 3.2，计算每一行的边缘位置
• 3.3，综合每一行图像计算拟合边缘位置的斜率和截距

4，边缘扩展函数ESF
5，线扩散函数LSF
6，傅里叶变换DFT

步骤详解

一、选择ROI

• 水平方向斜边
• 左暗右亮

引入我所用到的包

import cv2
from matplotlib.lines import lineStyles
import numpy as np
import matplotlib.pyplot as plt

导入图像

img = cv2.imread(r"C:\Users\Thinkpad\Desktop\123.jpg")
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)#使用灰度图像

二、线性化处理

gamma变换公式：
$$O=c\ast I^{\gamma }$$

• O:输出灰度级别
• I:输入灰度级别
• c:参数
• γ：调整参数，较小时降低对比度，较高时增加对比度

根据定义，线性化即为将图像还原
我这边简化了操作，把系数c设为1，那么将上面的输出O代入下面的公式
推导线性化公式即为：
$$Luminance=O^{1/\gamma }=(I^{\gamma }{)}^{1/\gamma }=I$$
不确定这线性化公式对不对，没有找到相关资料

def luminance_trance(image,gamma = 0.5):#这个gamma是可调整的，通常设为0.5左右
    return np.power(image,1/gamma)
img_lum = luminance_trance(gray)
cv2.imshow("img_correct ",img_lum )
cv2.waitKey(0)

得到的图像结果：

三、边缘位置估计

边缘位置公式贴图：

3.1 计算每一行的一维导数

对于得到的阵列中的每一行像素，线性化图像数据的导数是使用有限冲激响应（FIR）滤波器。结果是一个数组，它的大小与输入ROI大小相同。该导数矩阵的一维质心是逐行计算的，以确定每一条线上的边的位置。
将上述公式进行步骤拆分，先进行求导:
$$\Phi (p+1,r)-\Phi (p,r)$$
这里对于这个公式的使用我找到了不同的版本，基本是基于不同的滤波器,
当滤波器为[-1,1],将滤波器应用到每一行像素

这里介绍python函数：
np.diff()☞函数用法
这是个用于计算差分的函数，对于一维数组的差分则近似于求导

h,w = gray.shape#获取图像的高和宽
for i in range(h):#逐行求导
    single_row= img_lum[i].astype(int)#因为差分后会出现负数，所以将uint8转化为int
    diff_row = np.diff(single_row)
    hamming_window_1 = np.hamming(len(diff_row))#汉明窗
    diff_row = diff_row*hamming_window_1#将每一行的数据进行汉明滤波

当系数Φ=1/2,即滤波器为[-1/2,1/2]时np.diff()将不再实用，这里我没有找到合适的函数，尝试使用如下函数:
np.convolve()☞函数用法
这是一个处理卷积的函数
根据离散线性卷积的步骤，该函数会对卷积核做180度反转，然后再进行滑动滤波。因此，当滤波器为[-1/2,1/2]时，卷积核即为[1/2,-1/2]

kernel = np.array([0.5,-0.5])
for i in range(h):
    single_row= img_lum[i].astype(int)#因为差分后会出现负数，所以将uint8转化为int
    #diff_row = np.diff(single)
    diff_row = np.convolve(single_row,kernel,mode="valid")#自定义的求导
    hamming_window_1 = np.hamming(len(diff_row))#汉明窗
    diff_row = diff_row*hamming_window_1#将每一行的数据进行汉明滤波

在上述的方法中，我不仅进行了求导，还应用了汉明窗滤波

下面展示的时两种求导方式以及运用汉明窗的区别

[-1,1]滤波求导：

[-1/2,1/2]滤波求导

这里提一嘴汉明窗：☞内容参考
在ISO12233:2017标准中使用的是汉明窗，而ISO12233:2023标准中使用的是海宁窗
区别参见下图

3.2，计算每一行的边缘位置

每一行的边缘位置即为一阶导数的矩心
贴图公式简化为：
$$c(r)=\frac{{\sum }_{p=1}^{p}p\ast diff(p,r)}{{\sum }_{p=1}^{p}diff(p,r)}$$
其中diff即为上面求的导数
设置两个列表x = [],y= []分别存入边缘点横坐标和纵坐标，这里我使用[-1/2,1/2]滤波+汉明窗。
接之前的代码：

kernel = np.array([-0.5,0.5])
x = []
y = []
for i in range(h):
    single_row= img_lum[i].astype(int)#因为差分后会出现负数，所以将uint8转化为int
    #diff_row = np.diff(single)
    diff_row = np.convolve(single_row,kernel,mode="valid")#自定义的求导
    hamming_window_1 = np.hamming(len(diff_row))#汉明窗
    diff_row = diff_row*hamming_window_1#将每一行的数据进行汉明滤波
    dt1 = np.sum(diff_row)
    dt = result = np.sum(diff_row * (np.arange(len(diff_row))+1))#每个元素*（元素索引+1），因为索引从0开始所以需要+1
    shift = dt/dt1
    x.append(int(shift))
    y.append(i)

x,y中存放的是所有边缘点的合集

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

在获得x和y组成的坐标集合后，使用最小二乘法拟合直线
$$斜率b=\frac{{\sum }_{i=1}^n(x[i]-average(x))\ast (y[i]-average(y))}{{\sum }_{i=1}^n(x[i]-average(x){)}^2}$$
$$截距a=average(y)-b\ast average(x)$$
python有现成的公式np.polyfit()☞函数参考
获得的直线代码及直线效果图，x,y为前面获取的点坐标合集

coefficients = np.polyfit(x, y, 1)#线的拟合
slope = coefficients[0]#斜率
intercept = coefficients[1]#截距

cv2.line(img,(0,int(intercept)),(int(-intercept/slope),0),(0,0,255),2)#这个方法在斜率为正时画出的线段将不可见
cv2.imshow("2",img)
cv2.waitKey(0)

红色为拟合直线：

未完待续。。。