一种图卡图像的SFR自动提取与分析方法与流程

本发明涉及客观图像质量评价领域，更具体地，涉及一种图卡图像的sfr自动提取与分析方法。

背景技术：

sfr(spatialfrequencyresponse)算法衡量影像的锐利度(sharpness)。sfr数值越大，代表图像越锐利，此时图片也越清晰。本领域现阶段进行sfr提取的流程包括：

1)手动或者自动标定待评估区域(存在标记点的情况下)；

2)手动微调待评估区域；

3)基于刃边法或者改进的刃边法提取sfr。

3.1)提取评估区域的esf；

3.2)esf转换为lsf；

3.3).lsf进行傅里叶变换得到mtf；

3.4)对mtf做去噪、平滑，得到最终的sfr；

3.5)当评估区域多于1时，利用平均的sfr做图像整体sfr。

4)利用求得的sfr对图像锐度做评估。

以上方式在人机交互、评估图像数量较少的情况下，可以简单并且形象地分析sfr，但缺点是：

1)自动标定待评估区域方法不稳定，尤其是图像较多且评估图像中没有标记点时；

2)手动、自动标定评估区域过程中缺少自动微调功能，这会大大降低评估效率；

3)提取sfr噪声缺少实际图像中的噪声功率谱，导致得到的sfr中包含噪声，使评估结果并不准确。

本发明的参考文献为

文献1：陈起行,何斌,王俊琦.刃边法边缘扩散函数最优提取方法[j].激光与光电子学进展,2015,52(11):102-108.

文献2：http://www.imatest.com/products/imatest-master/。

技术实现要素：

本发明的目的在于克服现有技术存在的上述缺陷，提供一种图卡图像的sfr自动提取与分析方法，通过对图卡色彩及灰度的计算，自动判断图卡类别；对图卡的标识区域进行预估和优化，将人工主管微调过程转换为自动量化的优化过程，可以保证后续提取sfr的精度，将每个区域的sfr与噪声分布均做细致分离，在完成去噪后再进行平滑得到修正后的sfr，从而大大提升了提取所得sfr的可信度。

为实现上述目的，本发明的技术方案如下：

一种图卡图像的sfr自动提取与分析方法，其特征在于，包括以下步骤：

s1：使用待测试的成像系统对图卡进行拍摄，并保存拍摄图像数据；

s2：判定拍摄图像所拍摄的图卡的类别；

s3：基于步骤s2中判定的图卡类型，对标识区域做预先估计；

s4：进行待评估区域的优化，以步骤s3中确定的矩形评估区域初始值为基础，以下列公式为目标函数

其中r(x,y)代表的是以横纵坐标x,y为中心的外接矩形区域坐标，gdark(x,y)是深色图像在x,y处的灰度，glight(x,y)是浅色图像在x,y处的灰度；x0,y0为区域的起始中心坐标，△x,△y是区域的调整坐标；

s5：使用刃边法计算成像系数所拍摄图像直接的自动提取的待评估区域的调制传递函数mtf；

s6：计算步骤s4中自动提取的待评估区域的噪声分布；

s7：对步骤s5中计算得到的调制传递函数mtf进行降噪；

s8：利用数值方法对mtfcorr(r)做平滑，得到最终修正后的sfr。

进一步地，所述步骤s2具体包括以下步骤：

s21：基于lab空间的图像色彩能量评估，对图像做lab变换，然后计算ab空间的平均能量，设定阈值，阈值的参考范围是0.1～0.9，当平均能量大于该阈值时，拍摄图像所拍摄的图卡是iso12233：2000或者iso12233：2014，进入步骤s3；当平均能量小于该阈值时，进入步骤s22；

s22：计算图像的灰度直方图，通过大津法求取自适应二值化阈值，统计直方图小于该阈值的频率，若频率小于设定的频率阈值，则认为拍摄图像所拍摄的图卡为sfrreg图卡，否则为checkerboard图卡，进入步骤s3。

进一步地，所述步骤s3具体包括以下步骤：

s31：如果拍摄图像所拍摄的图卡是iso12233：2000或者iso12233：2014，拍摄图像的分割标识与过程包括以下步骤：

s311：计算拍摄图像每个像素位置的ab能量，得到拍摄图像的能量图；

s312：对能量图做最大值归一化，得到归一化图像，然后做归一化图像的直方图统计，以直方图累积统计频率为70％～99％范围中一具体值的能量作为阈值，对能量图做二值化，然后进行形态学处理，对二值图去噪，得到彩色区域块的蒙板图，并计算蒙板图的最小外接矩形，若矩形中心距离蒙板图的图像中心小于图像比例长度，比例长度参考范围0.125～0.25，则认为矩形中心接近蒙板图的图像中心，进入步骤s313，否则进入步骤s314；

s313：若矩形中心接近蒙板图的图像中心，则认为拍摄图像所拍摄的图卡为iso12233：2014，基于先验知识，利用该图卡中用于评估sfr的区域与中心、矩形的几何关系，估计出所有评估sfr区域的外接矩形，进入步骤s4；

s314：计算矩形中心的横坐标与蒙板图图像中心的横坐标距离，若距离与矩形长度比值小于设定比例阈值，该阈值参考范围是2.5～10，且矩形中心与蒙板图图像边缘棋盘距离与蒙版图图像长度比例小于比例阈值，该阈值参考范围是0.125～0.25，则认为拍摄图像所拍摄的图卡为iso12233：2000，基于先验知识，估计出所有评估sfr的外接矩形，进入步骤s4；若中心坐标均不满足s313与s314的条件，认为拍摄图像所拍摄的图卡非标准图卡，拒绝进行sfr自动分析；

s32：如果拍摄图像所拍摄的图卡是checkerboard图卡，拍摄图像的分割标识与过程包括以下步骤：

s321：进行拍摄图像边界提取与细化；

s322：基于霍夫变换做直线估计，统计每条直线的斜率与截距值；

s323：基于斜率做行方向与列方向的聚类，根据不同方向的直线描述，得到每条直线的交点，每相邻的4个交点，确定一个可以用来评估sfr的最小外接矩形区域，进入步骤s4；计算每个相邻角点间距离的平均值与标准差,若标准差与平均值比值大于比例阈值，该阈值参考范围为0.05～0.1，则认为该图像是非标准图卡或几何畸变过大的图卡，均拒绝做sfr自动分析；

s33：如果拍摄图像所拍摄的图卡是sfrreg图卡，拍摄图像的分割标识与过程包括以下步骤：

s331：基于mser方法统计sfr评估候选区域；

s332：对区域做基于空间几何位置阈值的合并，若合并后数量多于9个或者小于5个均认为图像是非标准图卡，拒绝做sfr自动分析，否则计算每个区域的外接矩形，用来做区域内sfr评估，进入步骤s4。

进一步地，所述步骤s4的具体优化方法包括以下步骤：

s41：sfrreg图卡以矩形的4个顶点坐标做起始中心坐标，checkerboard图卡、iso12233：2000图卡和iso12233：2014图卡以候选区域的外接矩形四条边中心做起始中心坐标，计算得到目标初始值δ0；

s42：对于矩形不同的边，在进行第k次优化时，k为优化循环次数，进行以下操作

对于矩形上边，先做纵坐标的优化，△y＝△y+1，然后计算δk，若δk>δk-1，△y＝△y-1；然后进行横坐标优化，△x＝△x+1，然后计算δk，若δk>δk-1，下次调整时按照△x＝△x-2进行调整，令δk＝δk-1，每次δk>δk-1均做一次优化方向的改变；当不论增加或者减少△x，均满足δk>δk-1，或者δk＝δk-1，则保留最小△x，迭代结束；

对于矩形下边，先做纵坐标的优化，△y＝△y-1，然后计算δk，若δk>δk-1，△y＝△y+1；然后进行横坐标优化，△x＝△x+1，然后计算δk，若δk>δk-1，下次调整时按照△x＝△x-2进行调整，令δk＝δk-1，每次δk>δk-1均做一次优化方向的改变；当不论增加或者减少△x，均满足δk>δk-1，或者δk＝δk-1，则保留最小△x，迭代结束；

对于矩形左边，先做横坐标的优化，△x＝△x+1，然后计算δk，若δk>δk-1，△x＝△x-1；然后进行纵坐标优化，△y＝△y+1，然后计算δk，若δk>δk-1，下次调整时按照△y＝△y-2进行调整，令δk＝δk-1，每次δk>δk-1均做一次优化方向的改变；当不论增加或者减少△y，均满足δk>δk-1，或者δk＝δk-1，则保留最小△y，迭代结束；

对于矩形右边，先做横坐标的优化，△x＝△x-1，然后计算δk，若δk>δk-1，△x＝△x+1；然后进行纵坐标优化，△y＝△y+1，然后计算δk，若δk>δk-1，下次调整时按照△y＝△y-2进行调整，令δk＝δk-1，每次δk>δk-1均做一次优化方向的改变；当不论增加或者减少△y，均满足δk>δk-1，或者δk＝δk-1，则保留最小△y，迭代结束；

s43：完成优化后，得到拍摄图像的自动提取的待评估区域。

进一步地，所述步骤s6具体包括以下步骤：

s61：使用大津法对自动提取的待评估区域进行二值化；

s62：规定深色、浅色区域分别为gdark(x,y)、glight(x,y)，其中，x∈[1,m]，y∈[1,n]，做傅里叶变换，得到二维的幅度谱gd(u,v)与gl(u,v)，其中u∈[1,m]，v∈[1,n]，其中m,n是傅里叶变换后的频谱行列值；

s63：根据下式算基于欧式距离的半径值r，

s64：按照半径值r与(u,v)的对应关系，将二维的幅度谱gd(u,v)和gl(u,v)转换为一维幅度谱gl(r)和gl(r)，

s65：通过求平均的方式将所得的噪声分布取均值，g(r)为对应区域相对精准的噪声能量分布：

进一步地，所述步骤s7具体包括以下步骤：

s71：每个所拍摄图像直接的自动提取的待评估区域的mtf(r)与对应区域的噪声分布g(r)相减去除0频分量的差值，得到去噪后的mtfcorr(r)：

s72：对于mtfcorr(r)小于0的全部置为0；

s73：对mtfcorr(r)做最大值归一化：

进一步地，所述步骤s8中所述的数值方法是三次样条插值或者样条拟合。

从上述技术方案可以看出，本发明通过自动判断图卡类别、自动最优化sfr提取区域和对噪声分布的细致分离，确保了sfr的精度和可信度。因此，本发明具有极大地提升了现阶段人机交互提取sfr的效率的显著特点。

附图说明

图1是本发明的方法流程图。

具体实施方式

下面结合附图，对本发明的具体实施方式作进一步的详细说明。

需要说明的是，在下述的具体实施方式中，在详述本发明的实施方式时，为了清楚地表示本发明的结构以便于说明，特对附图中的结构不依照一般比例绘图，并进行了局部放大、变形及简化处理，因此，应避免以此作为对本发明的限定来加以理解。

在以下本发明的具体实施方式中，请参阅图1，图1是本发明的流程图。

一种图卡图像的sfr自动提取与分析方法，其特征在于，包括以下步骤

使用待测试的成像系统对图卡进行拍摄，并保存拍摄图像数据；

s2：判定拍摄图像所拍摄的图卡的类别，具体地，包括以下步骤

s21：基于lab空间的图像色彩能量评估，对图像做lab变换，然后计算ab空间的平均能量，设定阈值，阈值参考值范围0.1～0.9，当平均能量大于该阈值时，拍摄图像所拍摄的图卡是iso12233：2000或者iso12233：2014，进入步骤s3；当平均能量小于该阈值时，进入步骤s22；

s22：计算图像的灰度直方图，通过大津法求取自适应二值化阈值，统计直方图小于该阈值的频率，若频率小于设定的频率阈值，则认为拍摄图像所拍摄的图卡为sfrreg图卡，否则为checkerboard图卡，进入步骤s3；

s3：基于步骤s2中判定的图卡类型，对标识区域做预先估计，具体地，包括以下步骤：

s31：如果拍摄图像所拍摄的图卡是iso12233：2000或者iso12233：2014，拍摄图像的分割标识与过程包括以下步骤：

s311：计算拍摄图像每个像素位置的ab能量，得到拍摄图像的能量图；

s312：对能量图做最大值归一化，得到归一化图像，然后做归一化图像的直方图统计，以直方图累积统计频率为95％的能量作为阈值，对能量图做二值化，然后进行形态学处理，对二值图去噪，得到彩色区域块的蒙板图，并计算蒙板图的最小外接矩形，若矩形中心距离蒙板图的图像中心小于图像比例长度，比例长度参考范围0.125～0.25，进入步骤s313，否则进入步骤s314；

s313：若矩形中心接近蒙板图的图像中心，则认为拍摄图像所拍摄的图卡为iso12233：2014，基于先验知识，可以利用该图卡中用于评估sfr的区域与中心、矩形的几何关系，即每个区块与提取矩形中心的距离比值、每个区块与提取矩形大小的比例比值是几乎不变的，则可以利用这些固定值估计出所有评估sfr区域的外接矩形，进入步骤s4，这里以实例举例说明计算方法，假定矩形中心为(x,y)，尺寸长度为h×w，其中某个sfr区域距离中心的比例为(△x,△y)，区域大小与矩形尺寸大小比例为(δx,δy)，则该区域的实际中心计算方式是，

(xsfr,ysfr)＝(△xw+x,△yh+y)

sfr区域实际大小的计算方式是，

(hsfr,wsfr)＝(δyh,δxw)

s314：计算矩形中心的横坐标与蒙板图图像中心的横坐标距离，若距离与矩形长度比值小于设定比例阈值，该阈值参考范围是2.5～10，且矩形中心与蒙板图图像边缘棋盘距离与蒙版图图像长度比例小于比例阈值，该阈值参考范围是0.125～0.25，则认为拍摄图像所拍摄的图卡为iso12233：2000，基于先验知识，估计出所有评估sfr的外接矩形，进入步骤s4；若中心坐标均不满足s313与s314的条件，认为拍摄图像所拍摄的图卡非标准图卡，拒绝进行sfr自动分析；

s32：如果拍摄图像所拍摄的图卡是checkerboard图卡，拍摄图像的分割标识与过程包括以下步骤：

s321：进行拍摄图像边界提取与细化；

s322：基于霍夫变换做直线估计，统计每条直线的斜率与截距值；

s323：基于斜率做行方向与列方向的聚类，根据不同方向的直线描述，得到每条直线的交点，每相邻的4个交点，确定一个可以用来评估sfr的最小外接矩形区域，进入步骤s4；计算每个相邻角点间距离的平均值与标准差,若标准差与平均值比值大于比例阈值，该阈值参考范围为0.05～0.1，则认为该图像是非标准图卡，或者几何畸变过大的图卡，均拒绝做sfr自动分析；

s33：如果拍摄图像所拍摄的图卡是sfrreg图卡，拍摄图像的分割标识与过程包括以下步骤：

s331：基于mser方法统计sfr评估候选区域，

s332：对区域做基于空间几何位置阈值的合并，若合并后数量多于9个或者小于5个均认为图像是非标准图卡，拒绝做sfr自动分析，否则计算每个区域的外接矩形，用来做区域内sfr评估，进入步骤s4；

s4：进行待评估区域的优化，以步骤s3中确定的矩形评估区域初始值为基础，以下列公式为目标函数

其中r(x,y)代表的是以横纵坐标x,y为中心的外接矩形区域坐标，gdark(x,y)是深色图像在x,y处的灰度，glight(x,y)是浅色图像在x,y处的灰度；x0,y0为区域的起始中心坐标，△x,△y是区域的调整坐标；

具体地优化方法包括以下步骤：

s41：sfrreg图卡以矩形的4个顶点坐标做起始中心坐标，checkerboard图卡、iso12233：2000图卡和iso12233：2014图卡以候选区域的外接矩形四条边中心做起始中心坐标，计算得到目标初始值δ0

s42：对于矩形不同的边，在进行第k次优化时，k为优化循环次数(为了使优化过程不出现死循环的情况，可以对最大优化次数做规定，比如100，当优化次数超过该数则也认为优化结束)，进行以下操作

对于矩形上边，先做纵坐标的优化，△y＝△y+1，然后计算δk，若δk>δk-1，△y＝△y-1；然后进行横坐标优化，△x＝△x+1，然后计算δk，若δk>δk-1，下次调整时按照△x＝△x-2进行调整，令δk＝δk-1，每次δk>δk-1均做一次优化方向的改变；当不论增加或者减少△x，均满足δk>δk-1，或者δk＝δk-1，则保留最小△x，迭代结束；

对于矩形下边，先做纵坐标的优化，△y＝△y-1，然后计算δk，若δk>δk-1，△y＝△y+1；然后进行横坐标优化，△x＝△x+1，然后计算δk，若δk>δk-1，下次调整时按照△x＝△x-2进行调整，令δk＝δk-1，每次δk>δk-1均做一次优化方向的改变；当不论增加或者减少△x，均满足δk>δk-1，或者δk＝δk-1，则保留最小△x，迭代结束；

对于矩形左边，先做横坐标的优化，△x＝△x+1，然后计算δk，若δk>δk-1，△x＝△x-1；然后进行纵坐标优化，△y＝△y+1，然后计算δk，若δk>δk-1，下次调整时按照△y＝△y-2进行调整，令δk＝δk-1，每次δk>δk-1均做一次优化方向的改变；当不论增加或者减少△y，均满足δk>δk-1，或者δk＝δk-1，则保留最小△y，迭代结束；

对于矩形右边，先做横坐标的优化，△x＝△x-1，然后计算δk，若δk>δk-1，△x＝△x+1；然后进行纵坐标优化，△y＝△y+1，然后计算δk，若δk>δk-1，下次调整时按照△y＝△y-2进行调整，令δk＝δk-1，每次δk>δk-1均做一次优化方向的改变；当不论增加或者减少△y，均满足δk>δk-1，或者δk＝δk-1，则保留最小△y，迭代结束；

s43：完成优化后，得到拍摄图像的自动提取的待评估区域；

s5：使用刃边法计算成像系数所拍摄图像直接的自动提取的待评估区域的调制传递函数mtf；

s6：计算步骤s4中自动提取的待评估区域的噪声分布，具体的，包括以下步骤：

s61：使用大津法对自动提取的待评估区域进行二值化；

s62：规定深色、浅色区域分别为gdark(x,y)、glight(x,y)，其中，x∈[1,m]，y∈[1,n]，做傅里叶变换，得到二维的幅度谱gd(u,v)与gl(u,v)，其中u∈[1,m]，v∈[1,n]，其中m,n是傅里叶变换后的频谱行列值；

s63：根据下式算基于欧式距离的半径值r，

s64：按照半径值r与(u,v)的对应关系，将二维的幅度谱gd(u,v)和gl(u,v)转换为一维幅度谱gl(r)和gl(r)，

s65：通过求平均的方式将所得的噪声分布取均值，g(r)为对应区域相对精准的噪声能量分布：

s7：对步骤s5中计算得到的调制传递函数mtf进行降噪，具体地包括以下步骤：

s71：每个所拍摄图像直接的自动提取的待评估区域的mtf(r)与对应区域的噪声分布g(r)相减去除0频分量的差值，得到去噪后的mtfcorr(r)：

s72：对于mtfcorr(r)小于0的全部置为0；

s73：对mtfcorr(r)做最大值归一化：

s8：利用数值方法对mtfcorr(r)进行做平滑，得到最终修正后的sfr。

数值方法包括三次样条插值或者样条拟合。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。