CN104933713A - 一种利用边缘分析的图像mtf估计方法 - Google Patents

Abstract

本发明公开了一种利用边缘分析的图像MTF估计方法，包括如下步骤：（1）盲估计点扩散函数PSF；（2）从PSF计算获取MTF。本发明方法针对单幅模糊退化图，首先利用双边滤波去除噪声与波纹等效应，进而采用WLS滤波提取边缘从而充分利用边缘信息，再次使用数学方法构建盲估计PSF的算法，进而换算到频域，实现单幅图像MTF的估计。在本发明方法中，只要一幅输入模糊图像，即可获取较好的MTF估计结果。本发明方法可应用于遥感成像质量、普通相机成像质量的评价，以实现对整个成像系统质量的评估。

Description

一种利用边缘分析的图像MTF估计方法

技术领域

本发明涉及图像处理与图像评价技术，尤其涉及一种利用边缘分析的图像MTF估计方法。

背景技术

MTF，全称Modulation Transfer Function，常称为调制传递函数，它本身是属于物理光学理论范畴的一个概念，同时也是衡量光学系统的质量的一个重要指标。MTF以空间频率的函数形式存在，较之仅凭借某一个数字量(分辨率、清晰度等)来对成像系统进行质量评价更具有可分析性、权威性。经过多年实践，用图像MTF来评价光学系统的成像质量已经获得普遍承认并得到广泛应用。

目前，国内外常用的MTF计算方法主要包括：冲击输入法、正弦输入法、相片判读法、脉冲法以及刃边法。其中，后两种方法适合后处理计算，比如在轨卫星图像MTF的计算。实际上，在很多应用情况下，系统的MTF并不只是由单纯的光学系统决定，后端的电路系统等都对MTF有影响，而图像是最终输出，从图像计算MTF的方法能够衡量整个成像系统的成像质量。平常使用的脉冲法以及刃边法，需要从图像中寻找合适的位置——即脉冲存在的位置与刃边存在的位置，有时候图像中可能不存在这样的位置与信息，而且有时候需要人工干预找出这些位置，这都将导致无法自动产生准确的MTF。

因此，对于从图像获取MTF，需要开发一种自动化的方法，即输入图像就可以输出MTF。

发明内容

本发明提出一种利用边缘分析的图像MTF估计方法，从单幅图角度，充分利用边缘信息，构建盲估计点扩散函数PSF(Point Spread Function)，进而换算到频域，实现单幅图像MTF的估计。

本发明的主要思路是：

(1)盲估计点扩散函数

对于图像，其常用的退化模型表示为：

$I=f\⊗h+N$

上式中观测退化图像I，f是原始场景即清晰景物，h即为点扩散函数PSF，N为噪声，是卷积符号；PSF即为Point Spread Function；

在输入观测退化图像I的情况下，进行点扩散函数PSF的估计即h的估计，忽略N的影响；

(1-a)：首先，利用双边滤波器进行图像滤波，观测退化图像I经双边滤波器后，得到去除噪声和波纹的图像g；

在这里，双边滤波器是一个从局部操作的过程，其主要方程如下：

$g(i,j)=\frac{1}{M(i,j)}\underset{p\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left\{G_s\right[D_s\left(p\right),{\mathrm{\σ}}_s\left]G_v\right[D_v\left(p\right),{\mathrm{\σ}}_v\left]I_p\right\}$

上式中，Ω是局部区域，也可称作支持域；(i,j)是像素坐标，且它是局部区域Ω的中心，p代表局部区域Ω中的任意像素；在这个公式中，I_p表示p处的像素强度值，并且G_s[D_s(p),σ_s]与G_v[D_v(p),σ_v]是两个高斯函数；归一化的因子M(i,j)定义为：

$M_{\mathrm{ij}}=\underset{p\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left\{G_s\right[D_s\left(p\right),{\mathrm{\σ}}_s\left]G_v\right[D_v\left(p\right),{\mathrm{\σ}}_v\left]\right\}$

高斯函数G_s[D_s(p),σ_s]是为空间约束设计的，D_s(p)是尺寸，σ_s是标准差或偏差；D_s(p)代表(i,j)与p之间的空间距离，用像素表示；而G_v[D_v(p),σ_v]为像素强度约束，D_v(p)是尺寸,而σ_v是标准差或偏差；D_v(p)是像素(i,j)与p之间的强度距离：

D_v(p)＝|I(i,j)-I_p|

(1-b)：对步骤(1-a)中的结果，进而利用加权均方最小二乘滤波WLS方法滤波，得到具有强烈边缘的图像g₁；WLS即为Weighted Least Square；

WLS滤波器对g进行滤波，滤波后图像g₁：

g₁＝(κ+λH)^-1g

这里D_x与D_y代表差分算子，κ表示单位矩阵；平滑加权因子w_x和w_y定义为：

$w_x\left(a\right)=\frac{1}{({\left|D_x{\left(g\right)}_L\right|}^a+\mathrm{\ϵ})}$

$w_y\left(a\right)=\frac{1}{({\left|D_y{\left(g\right)}_L\right|}^a+\mathrm{\ϵ})}$

其中(g)_L定义为g的log操作，参数a可以控制梯度的变化能力，而ε是一个正常数，为了避免计算中的病态；

(1-c)：从步骤(1-b)中，获取具有强烈边缘的图像g₁，结合观测退化图像I及g₁，用PSF估计算法得到相对精确的PSF；

PSF求取式子如下：

$h={\arg \min }_h{\left|\right|h\⊗g_1-I\left|\right|}^2+\mathrm{\β}{\left|\right|h\left|\right|}^2$

在上式中，β为Tikhonov正则化能量约束的系数；

于是，可以将PSF求取式子写成矩阵的形式：

h＝argm_hin(||Ah-b||²+β||h||²)

其中A是上式中g₁的Toeplitz矩阵形式，b为I的列堆砌的向量形式，求解上式可以得到：

(A^TA+βκ)h＝A^Tb

用共轭梯度法求上式中的h；初始的PSF可用尺寸为S×S,标准差为0.0001的高斯核；

于是，获得了PSF的估计，即求解获得了h；

(2)从PSF计算获取MTF；

从步骤(1)中获取了PSF h，对h进行傅里叶变换：

H＝FFT(h)

h对应到频域为H，而FFT为傅里叶变换；

假设坐标(x,y)为H的中心，截取主对角线；从左上角到(x,y)刚好是主对角线的一半，其分布即为MTF的分布；

画出MTF曲线图，纵坐标为值即左上角到(x,y)的一系列数值，而横坐标为数值个数，可归一化。

本发明方法针对单幅模糊退化图，首先利用双边滤波去除噪声与波纹等效应，进而采用WLS滤波提取边缘从而充分利用边缘信息，再次使用数学方法构建盲估计PSF的算法，进而换算到频域，实现单幅图像MTF的估计。在本发明方法中，只要一幅输入模糊图像，即可获取较好的MTF估计结果。本发明方法可应用于遥感成像质量、普通相机成像质量的评价，以实现对整个成像系统质量的评估。

附图说明

图1为本发明方法的操作流程框图；

图2为具体实施例图的原始观测模糊图像；

图3为具体实施例图的MTF曲线，即图2经本方面方法处理后获取的结果。

具体实施方式

面结合附图，通过具体实施例，对本发明的技术方案进行清楚、完整的描述。

利用本发明方法处理图像，如图1所示，输入观测模糊图像，即可得到MTF估计的结果。以图2(观测图像)为例(以下称为I),其主要操作步骤如下：

(1)盲估计点扩散函数

对于图像，其常用的退化模型表示为：

$I=f\⊗h+N$

上式中观测退化图像I，f是原始场景即清晰景物，h即为点扩散函数PSF，N为噪声，是卷积符号；PSF即为Point Spread Function；

在输入观测退化图像I的情况下，进行点扩散函数PSF的估计即h的估计，忽略N的影响；

(1-a)：首先，利用双边滤波器进行图像滤波，观测退化图像I经双边滤波器后，得到去除噪声和波纹的图像g；

在这里，双边滤波器是一个从局部操作的过程，其主要方程如下：

$g(i,j)=\frac{1}{M(i,j)}\underset{p\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left\{G_s\right[D_s\left(p\right),{\mathrm{\σ}}_s\left]G_v\right[D_v\left(p\right),{\mathrm{\σ}}_v\left]I_p\right\}$

上式中，Ω是局部区域，也可称作支持域；(i,j)是像素坐标，且它是局部区域Ω的中心，p代表局部区域Ω中的任意像素；在这个公式中，I_p表示p处的像素强度值，并且G_s[D_s(p),σ_s]与G_v[D_v(p),σ_v]是两个高斯函数；归一化的因子M(i,j)定义为：

$M_{\mathrm{ij}}=\underset{p\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left\{G_s\right[D_s\left(p\right),{\mathrm{\σ}}_s\left]G_v\right[D_v\left(p\right),{\mathrm{\σ}}_v\left]\right\}$

高斯函数G_s[D_s(p),σ_s]是为空间约束设计的，D_s(p)是尺寸，σ_s是标准差或偏差；D_s(p)代表(i,j)与p之间的空间距离，用像素表示；而G_v[D_v(p),σ_v]为像素强度约束，D_v(p)是尺寸,而σ_v是标准差或偏差；D_v(p)是像素(i,j)与p之间的强度距离：

D_v(p)＝|I(i,j)-I_p|

(1-b)：对步骤(1-a)中的结果，进而利用加权均方最小二乘滤波WLS方法滤波，得到具有强烈边缘的图像g₁；WLS即为Weighted Least Square；

WLS滤波器对g进行滤波，滤波后图像g₁：

g₁＝(κ+λH)^-1g

这里D_x与D_y代表差分算子，κ表示单位矩阵；平滑加权因子w_x和w_y定义为：

$w_x\left(a\right)=\frac{1}{({\left|D_x{\left(g\right)}_L\right|}^a+\mathrm{\ϵ})}$

$w_y\left(a\right)=\frac{1}{({\left|D_y{\left(g\right)}_L\right|}^a+\mathrm{\ϵ})}$

其中(g)_L定义为g的log操作，参数a可以控制梯度的变化能力，而ε是一个正常数，为了避免计算中的病态；

(1-c)：从步骤(1-b)中，获取具有强烈边缘的图像g₁，结合观测退化图像I及g₁，用PSF估计算法得到相对精确的PSF；

PSF求取式子如下：

$h={\arg \min }_h{\left|\right|h\⊗g_1-I\left|\right|}^2+\mathrm{\β}{\left|\right|h\left|\right|}^2$

在上式中，β为Tikhonov正则化能量约束的系数；

于是，可以将PSF求取式子写成矩阵的形式：

h＝argm_hin(||Ah-b||²+β||h||²)

其中A是上式中g₁的Toeplitz矩阵形式，b为I的列堆砌的向量形式，求解上式可以得到：

(A^TA+βκ)h＝A^Tb

用共轭梯度法求上式中的h；初始的PSF可用尺寸为S×S,标准差为0.0001的高斯核；

于是，获得了PSF的估计，即求解获得了h；

(2)、从PSF计算获取MTF。

从步骤(1)中获取了PSF h，本发明对h进行傅里叶变换：

H＝FFT(h)

H为h在频域，而FFT为傅里叶变换。

假设(x,y)为H的中心，本发明截取主对角线，从左上角到(x,y)刚好是主对角线的一半，其分布即为MTF的分布。

可以画出MTF曲线图，纵坐标为值即左上角到(x,y)的一系列数值，而横坐标为数值个数，归一化。MTF曲线如图3所示。在实例中，采用的参数如下：

本发明虽然已以较佳实施例公开如上，但其并不是用来限定本发明，任何本领域技术人员在不脱离本发明的精神和范围内，都可以利用上述揭示的方法和技术内容对本发明技术方案做出可能的变动和修改，因此，凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化及修饰，均属于本发明技术方案的保护范围。

Claims (1)

1.一种利用边缘分析的图像MTF估计方法，其特征在于，包括如下步骤：

(1)盲估计点扩散函数

对于图像，其常用的退化模型表示为：

$I=f\⊗h+N$

上式中观测退化图像I，f是原始场景即清晰景物，h即为点扩散函数PSF，N为噪声，是卷积符号；PSF即为Point Spread Function；

在输入观测退化图像I的情况下，进行点扩散函数PSF的估计即h的估计，忽略N的影响；

(1-a)：首先，利用双边滤波器进行图像滤波，观测退化图像I经双边滤波器后，得到去除噪声和波纹的图像g；

在这里，双边滤波器是一个从局部操作的过程，其主要方程如下：

$g(i,j)=\frac{1}{M(i,j)}\underset{p\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left\{G_s\right[D_s\left(p\right),{\mathrm{\σ}}_s\left]G_v\right[D_v\left(p\right),{\mathrm{\σ}}_v\left]I_p\right\}$

上式中，Ω是局部区域，也可称作支持域；(i,j)是像素坐标，且它是局部区域Ω的中心，p代表局部区域Ω中的任意像素；在这个公式中，I_p表示p处的像素强度值，并且G_s[D_s(p),σ_s]与G_v[D_v(p),σ_v]是两个高斯函数；归一化的因子M(i,j)定义为：

$M_{\mathrm{ij}}=\underset{p\∈\mathrm{\Ω}}{\mathrm{\Σ}}\left\{G_s\right[D_s\left(p\right),{\mathrm{\σ}}_s\left]G_v\right[D_v\left(p\right),{\mathrm{\σ}}_v\left]\right\}$

高斯函数G_s[D_s(p),σ_s]是为空间约束设计的，D_s(p)是尺寸，σ_s是标准差或偏差；D_s(p)代表(i,j)与p之间的空间距离，用像素表示；而G_v[D_v(p),σ_v]为像素强度约束，D_v(p)是尺寸,而σ_v是标准差或偏差；D_v(p)是像素(i,j)与p之间的强度距离：

D_v(p)＝|I(i,j)-I_p|

(1-b)：对步骤(1-a)中的结果，进而利用加权均方最小二乘滤波WLS方法滤波，得到具有强烈边缘的图像g₁；WLS即为Weighted Least Square；

WLS滤波器对g进行滤波，滤波后图像g₁：

g₁＝(κ+λH)^-1g

这里D_x与D_y代表差分算子，κ表示单位矩阵；平滑加权因子w_x和w_y定义为：

$w_x\left(a\right)=\frac{1}{({\left|D_x{\left(g\right)}_L\right|}^a+\mathrm{\ϵ})}$

$w_y\left(a\right)==\frac{1}{({\left|D_y{\left(g\right)}_L\right|}^a+\mathrm{\ϵ})}$

其中(g)_L定义为g的log操作，参数a可以控制梯度的变化能力，而ε是一个正常数，为了避免计算中的病态；

(1-c)：从步骤(1-b)中，获取具有强烈边缘的图像g₁，结合观测退化图像I及g₁，用PSF估计算法得到相对精确的PSF；

PSF求取式子如下：

$h=\arg {\min }_h{\left|\right|h\⊗g_1-I\left|\right|}^2+\mathrm{\β}{\left|\right|h\left|\right|}^2$

在上式中，β为Tikhonov正则化能量约束的系数；

于是，可以将PSF求取式子写成矩阵的形式：

$h=\arg \underset{h}{\min }({\left|\right|\mathrm{Ah}-b\left|\right|}^2+\mathrm{\β}{\left|\right|h\left|\right|}^2)$

其中A是上式中g₁的Toeplitz矩阵形式，b为I的列堆砌的向量形式，求解上式可以得到：

(A^TA+βκ)h＝A^Tb

用共轭梯度法求上式中的h；初始的PSF可用尺寸为S×S,标准差为0.0001的高斯核；

于是，获得了PSF的估计，即求解获得了h；

(2)从PSF计算获取MTF；

从步骤(1)中获取了PSF h，对h进行傅里叶变换：

H＝FFT(h)

h对应到频域为H，而FFT为傅里叶变换；

假设坐标(x,y)为H的中心，截取主对角线；从左上角到(x,y)刚好是主对角线的一半，其分布即为MTF的分布；

画出MTF曲线图，纵坐标为值即左上角到(x,y)的一系列数值，而横坐标为数值个数，可归一化。