少有人走的路在上一篇图像解析力算法—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的重要代码的注释到这里也就结束了，拖了挺久的，主要是前段时间也是有点忙，大家有啥不清楚或者错误指正可以留言一下。大家互相探讨一波吧~
————————————————
版权声明：本文为CSDN博主「Dylan_young」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/weixin_38419133/article/details/100726062
本文出自勇哥的网站《少有人走的路》wwww.skcircle.com，转载请注明出处！讨论可扫码加群：
