在上一篇
图像解析力算法—SFR（Spatial Frequency Response）源码分析(一）
中介绍了SFR的几个重要函数，接下来我们看一下主流程和其他函数。
4、sfrProc作用：计算SFR数值的主流程函数
-
short sfrProc (double **freq, double **sfr,
-
int *len,
-
double *farea,
-
unsigned
short size_x,
int *nrows,
-
double *slope,
int *numcycles,
int *pcnt2,
double *off,
double *R2,
-
int version,
int iterate,
int user_angle)
-
{
-
unsigned
short i, j, col, err =
0;
-
long pcnt;
-
double dt, dt1, sfrc, tmp, tmp2;
-
double *temp=
NULL, *shifts=
NULL, *edgex=
NULL, *Signal=
NULL;
-
double *AveEdge=
NULL, *AveTmp=
NULL;
-
long *counts=
NULL;
-
int nzero;
-
unsigned
short size_y;
-
unsigned
int bin_len;
-
double avar, bvar, offset1, offset2, offset;
-
double centroid;
-
int start_row, center_row;
-
double *farea_old;
-
double cyclelimit;
-
FILE *fp =
NULL;
-
-
size_y = *nrows;
-
-
/* Verify input selection dimensions are EVEN */
-
if (size_x%
2 !=
0) {
-
fprintf(
stderr,
"ROI width is not even. Does this really matter???\n");
-
return
1;
-
}
-
-
/* At least this many cycles required. */
-
/* For special iterative versions of the code, it can go lower */
-
if (iterate) cyclelimit =
1.0;
-
else cyclelimit =
5.0;
-
-
/* Allocate memory */
-
shifts = (
double *)
malloc(size_y*
sizeof(
double));
-
temp = (
double *)
malloc(size_y*
sizeof(
double));
-
edgex = (
double *)
malloc(size_y*size_x*
sizeof(
double));
-
Signal = (
double *)
malloc(size_y*size_x*
sizeof(
double));
-
-
if( !user_angle ) {
-
//定位矩心位置
-
err = locate_centroids(farea, temp, shifts, size_x, size_y, &offset1);
-
if (err !=
0) {
return
2; }
-
-
/* Calculate the best fit line to the centroids */
-
/*计算矩心的拟合曲线*/
-
err = fit(size_y, temp, shifts, slope, &offset2, R2, &avar, &bvar);
-
if (err !=
0) {
return
3; }
-
-
if (version)
-
MTFPRINT4(
"\nLinear Fit: R2 = %.3f SE_intercept = %.2f SE_angle = %.3f\n",
-
*R2, avar,
atan(bvar)*(
double)(
180.0/M_PI))
-
}
-
-
/* Check slope is OK, and set size_y to be full multiple of cycles */
-
//检查刀口斜率，以确保后面超采样的质量
-
err = check_slope(*slope, &size_y, numcycles, cyclelimit,
1);
-
-
/* Start image at new location, so that same row is center */
-
//调整ROI的行
-
center_row = *nrows/
2;
-
start_row = center_row - size_y/
2;
-
farea_old = farea;
-
farea = farea + start_row*size_x;
-
/* On center row how much shift to get edge centered in row. */
-
/* offset = 0.; Original code effectively used this (no centering)*/
-
if(user_angle)
-
offset = *off - size_x/
2;
-
else
-
offset = offset1 +
0.5 + offset2 - size_x/
2;
-
-
*off = offset;
-
if(version & ROUND || version & DER3)
-
offset +=
0.125;
-
-
if (err !=
0) {
/* Slopes are bad. But send back enough
-
data, so a diagnostic image has a chance. */
-
*pcnt2 =
2*size_x;
/* Ignore derivative peak */
-
return
4;
-
}
-
-
/* reference the temp and shifts values to the new y centre */
-
/* Instead of using the values in shifts, synthesize new ones based on
-
the best fit line. */
-
// 基于拟合结果更新shifts
-
col = size_y/
2;
-
for (i=
0; i < size_y; i++) {
-
shifts[i] = (*slope) * (
double)(i-col) + offset;
-
}
-
-
/* Don't normalize the data. It gets normalized during dft process */
-
//不要对数据进行归一化，数据在DFT处理过程中会被归一化
-
/* To normalize here, set dt = min and dt1 = max of farea data */
-
//如果要在这里初始化，将dt设置为最小值，dt1设置为最大值
-
dt =
0.0;
-
dt1 =
1.0;
-
-
if (version & ESFFILE)
-
fp = fopen(
"esf.txt",
"w");
-
-
/* Calculate a long paired list of x values and signal values */
-
pcnt =
0;
-
for (j =
0; j < size_y; j++) {
-
for (i =
0; i < size_x; i++) {
-
edgex[pcnt] = (
double)i - shifts[j];
//计算每个点到刀口的距离
-
Signal[pcnt] = ((farea[((j*(
long)size_x)+i)]) - dt)/(dt1-dt);
//归一化每个点的亮度
-
if ((version & ESFFILE) && edgex[pcnt] < size_x/
2 +
3 &&
-
edgex[pcnt] > size_x/
2 -
3)
-
fprintf(fp,
"%f %f\n", edgex[pcnt], Signal[pcnt]);
-
pcnt++;
-
}
-
}
-
-
if (version & ESFFILE)
-
fclose(fp);
-
-
/* Allocate more memory */
-
bin_len = (
unsigned
int)(ALPHA*size_x);
-
AveEdge = (
double *)
malloc(bin_len*
sizeof(
double));
-
AveTmp = (
double *)
malloc(bin_len*
sizeof(
double));
-
counts = (
long *)
malloc(bin_len*
sizeof(
long));
-
-
/* Project ESF data into supersampled bins */
-
//进行超采样，生成长度为size_x*ALPHA(4)的当行图像（ESF），保存在AveEdge中
-
nzero = bin_to_regular_xgrid((
unsigned
short)ALPHA, edgex, Signal,
-
AveEdge, counts,
-
size_x, size_y);
-
free(counts);
-
free(Signal);
-
free(edgex);
-
-
/* Compute LSF from ESF. Not yet centered or windowed. */
-
// 将ESF转换为差分图LSF， 并计算LSF的矩心
-
if ( (version&DER3) )
-
calculate_derivative( bin_len, AveTmp, AveEdge, ¢roid,
1);
-
else
-
calculate_derivative( bin_len, AveTmp, AveEdge, ¢roid,
0);
-
-
if (iterate ==
0) {
-
/* Copy ESF to output area */
-
for ( i=
0; i<bin_len; i++ )
-
farea_old[i] = AveTmp[i];
-
}
-
-
/* Find the peak/center of LSF */
-
//寻找LSF的最大值
-
locate_max_PSF( bin_len, AveEdge, &pcnt);
//这里获得了最大值的中心
-
-
free(shifts);
-
free(temp);
-
-
if(version)
-
MTFPRINT3(
"Off center distance (1/4 pixel units): Peak %ld Centroid %.2f\n",
-
pcnt-bin_len/
2, centroid-bin_len/
2)
-
-
if((version & PEAK) ==
0)
//忽略差分后的最大值
-
pcnt = bin_len/
2;
/* Ignore derivative peak */
-
else
-
MTFPRINT(
"Shifting peak to center\n")
-
-
/* Here the array length is shortened to ww_in_pixels*ALPHA,
-
and the LSF peak is centered and Hamming windowed. */
-
//将LSF的最大值移到中心位置，并加上汉明窗
-
apply_hamming_window((
unsigned
short)ALPHA, bin_len,
-
(
unsigned
short)size_x,
-
AveEdge, &pcnt);
-
-
/* From now on this is the length used. */
-
*len = bin_len/
2;
-
-
if (iterate ==
0) {
-
/* Copy LSF_w to output area */
-
for ( i=
0; i<bin_len; i++ )
-
farea_old[size_x*(
int)ALPHA+i] = AveEdge[i];
-
}
-
-
tmp =
1.0;
-
tmp2 =
1.0/((
double)bin_len) ;
-
-
/* Now perform the DFT on AveEdge */
-
/* ftwos ( nx, dx, lsf(x), nf, df, sfr(f) ) */
-
//ftwos( number, dx, *lsf, ns, ds, *sfr)
-
(
void) ftwos(bin_len, tmp, AveEdge, (
long)(*len),
-
tmp2, AveTmp);
-
-
if(*freq==
NULL)
-
(*freq) = (
double *)
malloc((*len)*
sizeof(
double));
-
if(*sfr==
NULL)
-
(*sfr) = (
double *)
malloc((*len)*
sizeof(
double));
-
-
for (i=
0; i<(*len); i++) {
-
sfrc = AveTmp[i];
-
(*freq)[i]= ((
double)i/(
double)size_x);
//每个点对应的频率
-
(*sfr)[i] = (
double) (sfrc/AveTmp[
0]);
//归一化sfr
-
}
-
-
/* Free */
-
free(AveEdge);
-
free(AveTmp);
-
-
*nrows = size_y;
-
*pcnt2 = pcnt;
-
-
return(
0);
-
}
5、apply_hamming_window作用：对lsf应用汉明窗，减少噪声
-
void apply_hamming_window( unsigned short alpha,
-
unsigned
int oldlen,
//size_x*4
-
unsigned
short newxwidth,
//size_x
-
double *AveEdge,
long *pcnt2)
-
{
-
long i,j,k, begin, end, edge_offset;
-
double sfrc;
-
-
/* Shift the AvgEdge[i] vector to centre the lsf in the transform window */
-
// 将Avgedge移到中心位置， 两边由于移动造成的空位由0补齐
-
edge_offset = (*pcnt2) - (oldlen/
2);
//不能理解为什么一定要反着计算，pcnt2肯定是小于oldlen/2吧。。
-
if (edge_offset !=
0) {
//cer=6
-
if (edge_offset <
0 ) {
//这里根据分析的话，应该edge_offset只会小于0啊。。 ↓
-
for (i=oldlen
-1; i > -edge_offset
-1; i--)
// [l l l l l l l max l l l l l]
-
AveEdge[i] = (AveEdge[i+edge_offset]);
// ↑ ↑ ↑
-
for (i=
0; i < -edge_offset; i++)
// left center=3 right
-
AveEdge[i] =
0.00;
/* last operation */
//offset = center-cer=-3
-
}
else {
// cer=6
-
// ↓
-
for (i=
0; i < oldlen-edge_offset; i++)
// [0 0 0 l l l l l l l max l l]
-
AveEdge[i] = (AveEdge[i+edge_offset]);
// ↑ ↑ ↑
-
for (i=oldlen-edge_offset; i < oldlen; i++)
// center=3
-
AveEdge[i] =
0.00;
-
}
-
}
-
/* Multiply the LSF data by a Hamming window of width NEWXWIDTH*alpha */
-
//将begin和end两侧的值用0填充,但是感觉没啥用
-
begin = (oldlen/
2)-(newxwidth*alpha/
2);
//上下看来，begin只会等于0，因为oldlen=bin_len, bin_len=size_x*alpha == newxwidth*alpha
-
if (begin <
0) begin =
0;
-
end = (oldlen/
2)+(newxwidth*alpha/
2);
-
if (end > oldlen ) end = oldlen;
-
for (i=
0; i< begin; i++)
-
AveEdge[i] =
0.0;
-
for (i=end; i< oldlen; i++)
-
AveEdge[i] =
0.0;
-
-
// 给begin和end之间的数据加上汉明窗
-
// 汉明窗 W(n,α ) = (1 -α ) - α cos(2*PI*n/(N-1)) ，（0≤n≤N-1）
-
// 一般情况下，α取0.46
-
// 下面计算方法等于窗发生了平移（故符号发生了变化）， 结果是一样的
-
for (i=begin,j = -newxwidth*alpha/
2; i < end; i++,j++) {
-
sfrc =
0.54 +
0.46*
cos( (M_PI*(
double)j)/(newxwidth*alpha/
2) );
-
AveEdge[i] = (AveEdge[i])*sfrc;
-
}
-
//将lsfbegin的位置左移到index0的位置
-
//但在代码中应该是不会起作用的，
-
if (begin !=
0)
/* Shift LSF to begin at index 0 (rather than begin) */
-
for (k=
0, i=begin; k<newxwidth*alpha; i++,k++)
-
AveEdge[k] = AveEdge[i];
-
-
}
6、ftwos作用：计算DFT，将空域转换到频域范围上。
-
unsigned short ftwos(long number, double dx, double *lsf,
-
long ns,
double ds,
double *sfr)
-
{
-
double a, b, twopi, g;
-
long i,j;
-
-
// n-1 k
-
// DFT ==> X[k] = Σ x[n]e^(-j2π - n)
-
// n=0 N
-
-
twopi =
2.0 * M_PI;
-
for (j =
0; j < ns; j++){
//ns=1/2*bin_len 前一半
-
g = twopi * dx * ds * (
double)j;
-
for (i =
0, a =
0, b =
0; i < number; i++) {
-
a += lsf[i] *
cos(g * (
double)(i));
-
b += lsf[i] *
sin(g * (
double)(i));
-
}
-
sfr[j] =
sqrt(a * a + b * b);
-
}
-
return
0;
-
}
好了，整个SFR的重要代码的注释到这里也就结束了，拖了挺久的，主要是前段时间也是有点忙，大家有啥不清楚或者错误指正可以留言一下。大家互相探讨一波吧~