SFR解析算法 - SFR_Calculation （C语言）

目录

一. 简介：

二. 软件操作流程：

三. 编译

四. 算法流程

五. 关键算法

六. 鸣谢

---

一. 简介：

SFR是"Spatial Frequency Response"（空间频率响应）的简称，用来表示图形清晰度的算法。本篇文章介绍的的代码实现来自https://github.com/RayXie29/SFR_Calculation，虽然严谨性无法和mitre_sfr算法媲美，但是算法遵循ISO12233，依然可以作为初学者学习的资料。

二. 软件操作流程：

如果读者没有自己的斜边图形，该算法的文件夹imgs包含了一个图形original_img.bmp可以使用，如下图。

下面是选择ROI区域，如下图示，

切割后的图形如下，

最后会生成一个csv文件，保存不同空间频率下的MTF值，如下，

0,1
0.0277778,0.997496
0.0555556,0.990619
0.0833333,0.979686
0.111111,0.965329
0.138889,0.947379
0.166667,0.928097
0.194444,0.909256
0.222222,0.889783
0.25,0.867951
0.277778,0.843465
0.305556,0.818012
0.333333,0.792307
0.361111,0.764379
0.388889,0.73298
0.416667,0.698145
0.444444,0.661594
0.472222,0.625916
0.5,0.590796
0.527778,0.556283
0.555556,0.522798
0.583333,0.490742
0.611111,0.457425
0.638889,0.422138
0.666667,0.38516
0.694444,0.351526
0.722222,0.322769
0.75,0.294444
0.777778,0.262373
0.805556,0.228476
0.833333,0.200449
0.861111,0.18049
0.888889,0.159404
0.916667,0.134192
0.944444,0.106814
0.972222,0.0801231
1,0.0798865

三. 编译

请参考下面的blog搭建编译环境，(134条消息) VS2019+OpenCV2环境搭建_jgw2008的专栏-CSDN博客;
代码下载地址为 https://github.com/RayXie29/SFR_Calculation, 注意在sfr.cpp中将头文件"ISOsfr.h"改为“sfr.h”。

四. 算法流程

下面是引用一篇mitre_sfr算法流程的图片，该算法和其近似度十分高。
引用地址1：https://blog.csdn.net/weixin_38419133/article/details/100019551, 
引用地址2：https://blog.csdn.net/jaych/article/details/50700576

1）在相机获取到的RGB图像数据中选择对应垂直斜边ROI，并将其转化为Gray ROI;
2）对得到的Gray ROI数据进行线性化，得到linear ROI数据；
3）计算每一行的centroid，即图像重心，这样能够得到每一行数据的边缘位置。
4）获取到的centroid进行线性回归，得到对应的边缘数据，包括斜率、截距等。
5）获取得到重新定位ROI，根据图像数据及其到边缘的距离，获得ESF。
6）对ESF进行4倍超采样，根据每个数据与它同一行的边缘位置的距离，分别放入4个容器中，获得4xESF。
7）对4xESF进行差分运算，获得LSF。
8）对LSF应用汉明窗，减少图像边缘及噪声影响。
9）进行DFT运算，获得OTF。
10）从OTF的实部获得SFR。

五. 关键算法

此份代码的优点是简洁，可以和ISO12233算法流程匹配度高，但严谨性缺失，读者可以在使用过程中参考经典的mitre_sfr算法进行完善。

int SFRCalculation(cv::Mat &ROI, double gamma)
{
	if (ROI.empty())
	{
		std::cerr << "Open the ROI image error" << std::endl;
		return 0;
	}

	int height = ROI.rows, width = ROI.cols;
	//Do the gamma decoding to eliminate the gamma encoded by camera device 
	de_Gamma(ROI, gamma);
	int i, j;

	double slope = 0, intercept = 0;
	
	//Center centroid offset
	double CCoffset = 0;
	std::vector<double> y_shifts(height);
	
	Calculate the shifts between Centroids and Centroid of image center	
	std::vector<double> Cen_Shifts = CentroidFind(ROI, y_shifts, &CCoffset);
	if (Cen_Shifts.empty()) { return 0; }

	//simple linear regression for slanted edge fitting
	SLR(Cen_Shifts, y_shifts, &intercept, &slope);
	std::cout << "Slope is " << slope << std::endl;

	//Truncate the number of rows of data to largest slope cycle which will have an integer number of full phase rotations
	ReduceRows(slope, &height);

	//update the CCoffset to the offset between original mid point of image and reference mid point we calculated
	CCoffset = CCoffset + 0.5 + intercept - width / 2;

	//Mapping the pixel value of original image into a sampling data which the length is 4 times of original image width
	//This step is for concentrating the amount of change of original pixel values
	int SamplingLen = width * 4;
	std::vector<double> OverSamplingData = OverSampling(ROI, slope, CCoffset, height, width, &SamplingLen);

	//Using hamming window to filter the ripple signal of two side of data
	OverSamplingData = HammingWindows(OverSamplingData, SamplingLen);

	//decrete four transform
	DFT(OverSamplingData, SamplingLen);
	width = int(SamplingLen / 4);
	double maxData = 0;
	for (i = 0; i < SamplingLen; ++i)
	{
		if (OverSamplingData[i] != 0)
		{
			maxData = OverSamplingData[i];
			break;
		}
	}

	for (int i = 0; i < SamplingLen; ++i) { OverSamplingData[i] /= maxData; }

	std::clog << "create ./ref/mtf.csv file" << std::endl;
	std::fstream mtf_file("./ref/mtf.csv", std::ios::out);
	
	for (i = 0; i <= width; ++i)
	{
		double frequency = (double)i / width;
		mtf_file << frequency << "," << OverSamplingData[i] << "\n";
	}

5.1 线性化ROI数据

在Sensor获得图像之后，呈现出来的图像由于要符合人眼的感觉，会对图像像素进行伽马变换，使得其变成非线性的像素数据，从而使图像的显示更加符合人眼的感受。所以，当我们要进行sensor的成像解析力分析时，要先将图像处理成没有经过伽马变换前的。

一般sensor会对raw图像进行一个2.2的gamma变换，若我们想恢复原始图像时，我们只需要进行一个1/2.2的gamma变换即可。

代码如下，首先检查图形是否单通道，然后进行变换
变换公式：目标值 = 255*((原始值/255) ** (1/gamma)).

注意，代码默认赋予的gamma值为1，这样图形数据没有进行任何变换；比较合理的gamma值为2.2

void de_Gamma(cv::Mat &Src, double gamma)
{
	std::clog << "de_Gamma ++" << std::endl;

	if (Src.channels() != 1) { return; }

	for (int i = 0; i < Src.rows; ++i)
	{
		uchar *SrcP = Src.ptr(i);
		for (int j = 0; j < Src.cols; ++j)
		{
			SrcP[j] = 255 * (pow((double)SrcP[j] / 255, 1 / gamma));
		}
	}
}

六. 鸣谢

(133条消息) 图像解析力算法—SFR（Spatial Frequency Response）原理分析(一)_Dylan的博客-CSDN博客

(134条消息) 【图像处理】SFR算法详解1_Image-CSDN博客_sfr算法

(134条消息) 【图像处理】SFR算法详解2_Image-CSDN博客_sfr算法详解

(134条消息) 【图像处理】SFR算法详解3_Image-CSDN博客_imatest sfr