在上一篇圖像解析力算法—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的重要代碼的注釋到這裡也就結束了，拖了挺久的，主要是前段時間也是有點忙，大家有啥不清楚或者錯誤指正可以留言一下。大家互相探讨一波吧~
![opencv028-輪廓發現[圖]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)