从零开始学习SFR--4.0
前言：因为课题涉及镜头质量检测，而现在镜头检测最普遍的方法便是MTF曲线作为检测镜头质量的标准。网上相关的学习资料并不多，也有一些大佬做了相关算法的研究，不过零零散散，难以成系统。为了学习并实现相关算法，参考各大佬的文章，对整个学习思路进行整理，特开此贴作为学习笔记。
从这篇继上一篇对sfr_iso.c文件注解之后，继续对find_area.c和mitre_sfr.c文件中的代码进行注解，当中参考了很多大佬的文章。
代码注解参考：
链接: Mitre_sfr代码注解（一）main函数
链接: Mitre_sfr代码注解（二） sfrproc函数
链接: Mitre_sfr代码注解（三） 将ROI超采样为ESF（locate_centroids函数/fit函数和bin_to_regular_xgrid函数）
链接: Mitre_sfr代码注解（四） LSF / 汉明窗 / SFR（DFT）计算
find_area.c
这里从拍摄的图像（一般为RGB）中截取ROI区域，并将其转换为Gray。
有个问题就是：用于SFR算法的图像大小多大比较合适？ROI的区域大小信息并没有严格定义，但会有影响，具体各个方面的要求可参考如下：
对于定义ROI，有以下几个要求：（针对垂直的斜边）
1） 对于垂直的斜边，ROI的宽度要小于高度。
2）在ROI中，只允许出现一个斜边，并且大致居中。
3）在底部或者顶部，比例小的那一部分不能低于5个像素
4）在中心位置，边界到中间斜边的距离需要在20至60个像素之间
5）高度需要在80到300之间。
6）ROI中避免污点，保证图像边界基本连续。
/*****************************************************************************/
/* External functions */
unsigned short fit(unsigned long, double *, double *, double *, double *,
double *, double *, double *);
unsigned short locate_centroids(double *, double *, double *,
unsigned short, unsigned short, double *);
unsigned short check_slope(double, unsigned short *, int *, double, int);
/* Internal function */
void find_10percent(double *, int, double, int *, int *);
/*****************************************************************************/
/* Data passed to this function is assumed to be radiometrically corrected, */
/* and oriented vertically, with black on left, white on right. The black to */
/* white orientation doesn't appear to be terribly important in this code. */
//传递给该函数的数据被认为是经过辐射校正的，并垂直定向，黑色在左边，白色在右边。
//在这段代码中，黑白方向似乎并不十分重要。
/* Parameters:
Input: farea = radiometric image data in original ROI
size_x = number of columns in original ROI
size_y = number of rows in original ROI
verbose Whether ROI search messages are desired
0 = No printout
1 = Print info as better ROIs are found
Output: ncols = number of columns in recommended ROI
nrows = number of rows in recommended ROI
center_x = column center of recommended ROI
center_y = row center of recommended ROI
*/
/*
参数：
输入：farea=原始ROI中的辐射图像数据
size_x=原始ROI中的列数
size_y=原始ROI中的行数
是否需要ROI搜索消息
0=无打印输出
1=找到更好的ROI时打印信息
输出：ncols=建议ROI中的列数
nrows=建议ROI中的行数
center_x=建议ROI的列中心
center_y=建议ROI的行中心
*/
/*****************************************************************************/
short find_area (double *farea, unsigned short size_x, unsigned short size_y,
int *ncols, int *nrows, int *center_x, int *center_y,
unsigned char verbose)
{
int numcycles;
unsigned short len, i, dummy, err = 0;
double *temp=NULL, *shifts=NULL;
double R2, avar, bvar, offset1, offset2;
double slope, min, expect_edge;
int minlen, mincenter, center, signflag, cnt;
int center_row, check_row, half_width;
int min_x, max_x, first1, first2, last1, last2;
double *tptr, *sptr;
*nrows = size_y;
*ncols = size_x;
/* Allocate memory */
shifts = (double *)malloc(size_y*sizeof(double));
temp = (double *)malloc(size_y*sizeof(double));
center_row = *nrows/2;
/* These values should change during correct operation */
mincenter = -1;
minlen = -1;
err = locate_centroids(farea, temp, shifts, size_x, size_y, &offset1);
min = 0.0;
slope = -1./10.;
signflag = 1;
for (i = 0; i<center_row - 1.0/fabs(slope)/2;) {
center = center_row + signflag*i;
len = (center< size_y-center)? center: size_y-center;
sptr = shifts + center - len;
tptr = temp + center - len;
len = 2*len+1;
numcycles = 5;
cnt = 0;
for (; numcycles>4; len--, sptr++, tptr++) {
/* Calculate the best fit line to the centroids */\
//计算与质心的最佳拟合线
if (len <= size_y) {
err = fit(len, tptr, sptr, &slope, &offset2, &R2, &avar, &bvar);
/* Check slope is OK, and set size_y to be full multiple of cycles */
//检查斜率是否正常，并将size_y设置为周期的全倍数
dummy = len;
numcycles = 0;
err = check_slope(slope, &dummy, &numcycles,5.0, 0);
if (numcycles > 4 && R2 > min) {
if(verbose)
MTFPRINT7("Center row: %d Rows: %d Angle %.3f R2 = %f #Cycles=%d SE = %.3f\n",
center, len, atan(fabs(slope))*180.0/M_PI,
R2, numcycles, atan(bvar)*180/M_PI)
min = R2;
minlen = len;
mincenter = center;
}
if (numcycles > 4) cnt++;
}
len--;
/* Calculate the best fit line to the centroids */
//计算与质心的最佳拟合线
err = fit(len, tptr, sptr, &slope, &offset2, &R2, &avar, &bvar);
/* Check slope is OK, and set size_y to be full multiple of cycles */
//检查斜率是否正常，并将size_y设置为周期的全倍数
dummy = len;
numcycles = 0;
err = check_slope(slope, &dummy, &numcycles, 5.0, 0);
if (numcycles > 4 && R2 > min) {
if(verbose)
MTFPRINT7("Center row: %d Rows: %d Angle %.3f R2 = %f #Cycles=%d SE = %.3f\n",
center, len, atan(fabs(slope))*180.0/M_PI,
R2, numcycles, atan(bvar)*180/M_PI)
min = R2;
minlen = len;
mincenter = center;
}
if (numcycles > 4) cnt++;
}
signflag *= -1;
if (signflag == -1) i++;
if (signflag == -1 && cnt == 0) break;
}
sptr = shifts + mincenter - minlen/2;
tptr = temp + mincenter - minlen/2;
err = fit(minlen, tptr, sptr, &slope, &offset2, &R2, &avar, &bvar);
/* Check slope is OK, and set size_y to be full multiple of cycles */
dummy = minlen;
numcycles = 0;
err = check_slope(slope, &dummy, &numcycles,5.0, 1);
if(verbose)
MTFPRINT7("Best Center: %d Hgt: %d Angle %.3f R2 = %f #Cycles=%d SE = %.3f\n",
mincenter, minlen, atan(fabs(slope))*180.0/M_PI,
R2, numcycles, atan(bvar)*180/M_PI)
*center_x = *ncols/2;
*center_y = mincenter;
*nrows = minlen;
/* Figure out where edge is on first/last line, and then find
where derivative goes consistently (5pixels) below 10% that value */
//计算出边在第一行/最后一行的位置，然后找出导数一致（5像素）低于该值10%的位置
check_row = mincenter - dummy/2;
expect_edge = (double)(check_row-center_row)*slope + offset2 + offset1;
find_10percent(farea+check_row*size_x,size_x,expect_edge,&first1, &last1);
check_row = mincenter + dummy/2;
expect_edge = (double)(check_row-center_row)*slope + offset2 + offset1;
find_10percent(farea+check_row*size_x,size_x,expect_edge,&first2, &last2);
max_x = (last1 < last2)? last2: last1;
min_x = (first1 < first2)? first1: first2;
if(verbose)
MTFPRINT4("Needs from: %d to %d = %d cols \n",min_x, max_x, max_x-min_x)
if (max_x-min_x+31 < *ncols) { /* Try adding 15 pixels beyond ends */
*ncols = max_x-min_x+30;
*center_x = (max_x+min_x)/2;
half_width = (*ncols+1)/2;
if (*center_x < half_width) *center_x = half_width;
if (size_x - *center_x < half_width) *center_x = size_x - half_width;
if ((*ncols)%2 != 0) (*ncols)++;
}
else if (max_x-min_x > 32) {
*ncols = max_x-min_x;
*center_x = (max_x+min_x)/2;
if ((*ncols)%2 != 0) (*ncols)++;
}
else { /* Can't add 15 at each end, but not as big as 32 yet */
/* Need to add enough to get up to 32 */
*center_x = (max_x + min_x)/2;
if (*center_x < 16) *center_x = 16;
else if (*center_x > size_x-16) *center_x = size_x - 16;
if (size_x>= 32) *ncols = 32;
}
free(temp);
free(shifts);
return(0);
}
mitre_sfr.c
在获取图像以后还要对图像进行线性化处理。
在Sensor获得图像之后，呈现出来的图像由于要符合人眼的感觉，会对图像像素进行伽马变换，使得其变成非线性的像素数据，从而使图像的显示更加符合人眼的感受。所以，当我们要进行sensor的成像解析力分析时，要先将图像处理成没有经过伽马变换前的。
主要是对获取到的斜边数据进行伽马变换，获取到线性的图像数据。
为什么要获取线性的图像数据？
在SFR的标准算法中，建议的算法是用OECF将图像数据转化为线性数据。因此需要讨论下我们获取到的数据是否线性数据。线性跟非线性应该怎样去理解。
通常，sensor获取到的数据是线性的，sensor上每个pixel都是根据物体所呈现出来的颜色强度转化为对应的数值。但是，由于人眼对颜色的感知不是线性的，因此，通常在相机内部处理中，都有对转化后的图像进行伽马变换，变换后的图像就不是线性的，使得图像数据能够适应人眼。而通常取值为2.2的伽马变换。
为了获取到线性的图像数据，需要抵消这个伽马变换，在算法上，我们可以采取1/2.2的伽马变换来抵消。
int main(int argc, char **argv)
{
char problem_string[82];
unsigned char rotation;
int i;
double *farea;
double slope, scale, b;
int size_x, size_y;
int len, err, bin_len;
int rgt_side, left_side;
int center;
int numcycles=0;
double *Freq=NULL;
double *disp=NULL;
double *ref_lut;
double off, R2;
int center_x, center_y, new_x, new_y;
int lowest_val, highest_val, grey_level, first, last;
int piv_err;
TIFF* tif = NULL;
#if !defined(MSDOS)
char *fnamep;
fnamep = dirname(argv[0]);
if (argc==1 && fnamep != NULL && fnamep[0] == '/') {
printf ("Changed to directory of executable: %s\n\n", fnamep);
chdir((const char *)fnamep);
}
#else
atexit(wait_to_exit);
#endif
/* Once-only initializations */
g_scan_image_file_id = 0;
g_problem_count = 0; /* # of problems in our problem report */
g_IQS_problem_count = 0;
get_switches(argc, argv, image_filename, data_filename);
/* always append output to this file */
if ((g_mtfout = fopen(MTFOUTNAME,"a")) == NULL) {
fprintf(stderr,"Can't open %s\n",MTFOUTNAME);
exit(1);
}
/* command line args */
get_args(image_filename,data_filename,&tif); // 获取源文件名并读取文件头， 获取ROI位置， 如果需要读取计算角度的两个点位
/* print the user entered data (corner coords, etc.) in the output */
/* 时间， 版本号， 以及get_args函数里面获取到的信息*/
print_header(image_filename,data_filename);
g_target_res = 500;
if (abs((int)g_ppi - g_target_res) > 10) {
sprintf(problem_string,
"Resolution outside the %d-%d ppi PIV spec range\n",
g_target_res-10,g_target_res+10);
put_problem(problem_string, IQS);
put_problem("\n", IQS);
}
// 从ROI 4个角落各读取4个像素点， 用于判断刀口方向， 并确认ROI区域的合法性
input_area(tif, &rotation,
g_test_pattern_yul, g_test_pattern_ylr,
g_test_pattern_xul, g_test_pattern_xlr);
/* read the image in - we store it internally in a standard way */
// 读取ROI图像， 如果并将水平刀口转换为垂直刀口
read_in_image(tif,rotation,image_filename,&size_x,&size_y);
/* get reflectance */
// 创建反射查找表， 将图像灰阶映射在0~1范围内
ref_lut = read_in_ref_lut(data_filename, &lowest_val, &highest_val);
// 将原图ROI映射到反射ROI
len = size_x*size_y;
farea = (double *)malloc(len*sizeof(double));
for(i=0; i<len; i++) {
grey_level = g_image_array[i];
if (grey_level < lowest_val || grey_level > highest_val) {
MTFPRINT("ERROR: OECF range does not cover the entire image\n")
MTFPRINT2(" Greylevel %d found in image, but not OECF\n", grey_level)
exit(-1);
}
farea[i] = (double)ref_lut[grey_level];
}
if(g_autorefine) {
find_area(farea, (unsigned short)size_x, (unsigned short)size_y, &new_x, &new_y, ¢er_x, ¢er_y, g_extended);
if (new_x != size_x || new_y != size_y ||
center_x != size_x/2 || center_y != size_y/2) {
center_x += (g_test_pattern_xlr + g_test_pattern_xul)/2 - size_x/2;
center_y += (g_test_pattern_ylr + g_test_pattern_yul)/2 - size_y/2;
g_test_pattern_xul = center_x - new_x/2;
g_test_pattern_yul = center_y - new_y/2;
g_test_pattern_xlr = g_test_pattern_xul + new_x;
g_test_pattern_ylr = g_test_pattern_yul + new_y;
read_in_image(tif,rotation,image_filename,&size_x,&size_y);
len = size_x*size_y;
for(i=0; i<len; i++)
farea[i] = (double)ref_lut[g_image_array[i]];
}
MTFPRINT("Refined coordinates of slanted edge area:\n")
MTFPRINT3("Upper Left Col =\t\t%d\tUpper Left Row =\t%d\n",
g_test_pattern_xul,g_test_pattern_yul)
MTFPRINT3("Center Col =\t\t\t%d\tCenter Row =\t\t%d\n",
(g_test_pattern_xlr+g_test_pattern_xul)/2,
(g_test_pattern_ylr+g_test_pattern_yul)/2)
MTFPRINT3("Width =\t\t\t\t%d\tHeight =\t\t%d\n",
(g_test_pattern_xlr-g_test_pattern_xul),
(g_test_pattern_ylr-g_test_pattern_yul))
MTFPRINT("\n\n")
}
if (g_userangle) {
double fslope;
slope = (g_pt1x-g_pt2x)/(double)(g_pt1y-g_pt2y);
b = g_pt1x - slope*g_pt1y;
if (rotation == RIGHT)
off = (g_test_pattern_yul + size_y/2)*slope + b - g_test_pattern_xul;
else {
off = (g_test_pattern_xlr - 1 - size_y/2 - b ) / slope - g_test_pattern_yul;
slope = -1.0/slope;
}
fslope = fabs(slope)*size_y/2;
if( off-fslope < 2 || off+fslope >= size_x-2 ) {
MTFPRINT("** ERROR: User specified edge does not work with ROI **\n");
if (off-fslope < 2)
MTFPRINT3("Edge Angle = %.3f degrees located %d pixels before ROI\n",
(atan(slope)*(double)(180.0/M_PI)), (int)rint(off))
else
MTFPRINT3("Edge Angle = %.3f degrees located %d pixels after ROI\n",
(atan(slope)*(double)(180.0/M_PI)), (int)rint(off-size_x))
if (g_debug) {
MTFPRINT("See location of edge points in diagnostic image\n")
write_debug_image(g_debug_array,g_debug_fullwidth, g_debug_fullheight);
}
exit(-1);
}
}
// 计算ROI中超出取值范围[0, 255]外的像素比例
// 前面像素映射的时候已经逐点确认过是否所有像素在oecf范围内， 感觉这个函数没什么实际作用
clipping (0, 255, 0.02, g_image_array, len);
/* calculate the sfr on this area */
err = sfrProc(&Freq, &disp, &bin_len, farea, (unsigned short)size_x, &size_y, &slope, &numcycles,¢er,&off, &R2, g_version, 0, g_userangle);
/* Add messages to problem report */
slope_bounds( slope, size_y, numcycles, 5.0, problem_string);
if(g_debug) {
draw_lines(g_debug_array, slope, rotation, center, size_y, off, size_x);
write_debug_image(g_debug_array,g_debug_fullwidth, g_debug_fullheight);
}
if (!g_userangle)
MTFPRINT2("R2 of linear edge fit = %.3f\n", R2)
else
MTFPRINT("User input edge location: No R2 available\n")
MTFPRINT4("Edge Angle = %.3f degrees \nCycle Length = %.3f \t#Cycles = %d\n",
(atan(slope)*(double)(180.0/M_PI)), fabs(1.0/slope), numcycles )
left_side = size_x/2+off; rgt_side = size_x/2-off;
if (left_side > bin_len/4) left_side = bin_len/4;
if (rgt_side > bin_len/4) rgt_side = bin_len/4;
MTFPRINT2("SFR computed using ~%d pixels on left or top side of the edge,\n", left_side)
MTFPRINT2(" and ~%d pixels on right or bottom side of the edge,\n", rgt_side)
MTFPRINT(" all within the ROI\n")
if (rotation == RIGHT)
MTFPRINT3(" over %d rows, centered at row %d\n", size_y, (g_test_pattern_yul + g_test_pattern_ylr)/2 )
else
MTFPRINT3(" over %d cols, centered at col %d\n", size_y, (g_test_pattern_xul + g_test_pattern_xlr)/2 )
if (center/4 < 16) {
MTFPRINT("\nERROR: Too few pixels across the edge for valid SFR computation\n\n")
exit(-1);
}
if ( left_side < 20 || rgt_side < 20)
MTFPRINT("Warning: Low width across the edge. SFR values may be suspect.\n")
/* log any problems we encountered into output */
if (g_reversepolarity)
MTFPRINT("\nNOTE: Original image polarity was reversed by user request.\n")
if (err) {
MTFPRINT ("** ERROR in computation. SFR values unknown **\n\n")
exit(-1);
}
scale = g_ppi/MM_PER_INCH;
MTFPRINT("\n\ncy/mm \t SFR ")
if(rotation == RIGHT)
MTFPRINT("\t Vert ")
else
MTFPRINT("\t Horz ")
if(farea[4*size_x-1] - farea[0] > 0)
MTFPRINT("black-to-white** edge\n(obj.plane)\n")
else
MTFPRINT("white-to-black** edge\n(obj.plane)\n")
piv_err=0;
for( i=0; i<bin_len/2; i++) {
double f, sfr;
double freq, fd_scale;
// 对加海明窗导致的功率损失进行补偿？
freq = M_PI*Freq[i];
/* [-1 0 0 0 1] freq /= 1.0; */
if (g_version & 4) /* [-1 0 1] */
freq /= 2.0;
else freq /= 4.0; /* [-1 1] */
if (freq == 0.0) fd_scale = 1.0;
else fd_scale = freq / sin(freq);
f = Freq[i]*scale;
sfr = disp[i]*fd_scale;
MTFPRINT3("%9.6f \t%f ",f, sfr)
if(Freq[i] > 0.5)
MTFPRINT("\t#\n") /* Mark frequencies above Nyquist */
else {
double lower = 1.02829 - 1.67473e-1*f + 1.06255e-2*f*f - 2.80874e-4*f*f*f;
double upper = 1.12;
if ((f<=10 && f >= 0.998) && (sfr < lower || sfr > upper)) {
if (g_nocompare)
MTFPRINT("\n")
else {
piv_err = 1;
MTFPRINT("\t*\n") /* Mark out of PIV spec values */
}
sprintf(problem_string,
"Computed SFR at %.3f cy/mm is outside PIV spec range of %.3f - %.2f\n",
f,lower,upper);
put_problem(problem_string, IQS);
}
else MTFPRINT("\n")
}
}
MTFPRINT("\n# Frequencies above Nyquist\n")
if (piv_err)
MTFPRINT("\n* SFR is outside PIV spec. See problem report for details.\n")
first = reverse_lut(ref_lut, farea[0], lowest_val, highest_val);
last = reverse_lut(ref_lut, farea[4*size_x-1], lowest_val, highest_val);
if (rotation == RIGHT) {
MTFPRINT2("** Edge type based on leftside gray level (%d)\n", first)
MTFPRINT2(" and rightside gray level (%d)\n", last)
}
else {
MTFPRINT2("** Edge type based on top gray level (%d)\n", first)
MTFPRINT2(" and bottom gray level (%d)\n", last)
}
MTFPRINT(" If type isn't correct, use 'f' option to adjust polarity\n\n")
print_problems();
if (g_extended){
int offset, begin, end;
int j, cnt;
double threshold,max_edge,current;
/* Derivatives below this value is considered insignificant for printing */
max_edge = fabs(farea[4*size_x+bin_len]);
threshold = max_edge / 25.0;
/* Find 5 contiguous insignificant values first occur before center */
offset = center - bin_len;
j = bin_len; cnt = 0;
while (cnt < 5 && j>=0) {
current = fabs(farea[4*size_x+j]);
if ( current < threshold ) cnt ++;
else cnt = 0;
if (current > max_edge) max_edge = current;
j--;
}
if (j<0) j=0;
begin = j + offset;
/* Find 5 contiguous insignificant values first occur after center */
j = bin_len; cnt = 0;
while (cnt < 5 && j<bin_len*2) {
current = fabs(farea[4*size_x+j]);
if ( current < threshold ) cnt ++;
else cnt = 0;
if (current > max_edge) max_edge = current;
j++;
}
if (j==2*bin_len) j=2*bin_len-1;
end = j + offset;
if ( farea[4*size_x+bin_len] < 0 ) max_edge *= -1;
MTFPRINT("\n\nPixel+ \t\t LSF_w \t\t AvgESF")
if (g_version&4)
MTFPRINT(" \t \n")
else
MTFPRINT("(with -0.125 offset)\n")
for( i=begin, j=begin-offset; i<0; j++, i++)
MTFPRINT3("% .2f \t\t% f\n", (i-center)/4.0, farea[4*size_x + j]/max_edge)
for( ; i<end && i<size_x*4; j++, i++) {
MTFPRINT4("% .2f \t\t% f \t% f\n", (i-center)/4.0, farea[4*size_x + j]/max_edge, farea[i])
}
for(; i<end; j++, i++)
MTFPRINT3("% .2f \t\t% f\n", (i-center)/4.0, farea[4*size_x + j]/max_edge)
MTFPRINT("\n")
MTFPRINT("+ Pixel values are centered on the computed line fit\n")
MTFPRINT(" LSF_w values are scaled to a maximum of 1\n")
if (!(g_version&4))
MTFPRINT(" Actual pixel values for the AvgESF are shifted up +0.125 from values here\n")
MTFPRINT("\n")
}
return(EXIT_SUCCESS);
} /* end of main */
1133
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
