少有人走的路在前面的文章中，我们已经分析了SFR的算法原理与步骤，下面我们直接来分析源码，源码中主要的函数主要分为一下几个：
1、locate_centroids作用：定位每一行像素的矩心位置
unsigned short locate_centroids(double *farea, double *temp, double *shifts,unsigned short size_x, unsigned short size_y,double *offset)
unsigned long i, j;
double dt, dt1, dt2;
/* Compute the first difference on each line. Interpolate to find the
centroid of the first derivatives. */
//计算差分 计算矩心位置
for (j = 0; j < size_y; j++) {
dt = 0.0;
dt1 = 0.0;
for (i = 0; i < size_x-1; i++) {
dt2 = farea[(j*(long)size_x)+(i+1)] - farea[(j*(long)size_x)+i];
dt += dt2 * (double)i;
dt1 += dt2;
}
shifts[j] = dt / dt1;
printf("=========\n");
printf("%f\n", shifts[j]);
}
/* check again to be sure we aren't too close to an edge on the corners.
If the black to white transition is closer than 2 pixels from either
side of the data box, return an error of 5; the calling program will
display an error message (the same one as if there were not a difference
between the left and right sides of the box ) */
if (shifts[size_y-1] < 2 || size_x - shifts[size_y-1] < 2) {
fprintf(stderr,"** WARNING: Edge comes too close to the ROI corners.\n");
// return 5;
}
//防止矩心过于靠近图像的边界
if (shifts[0] < 2 || size_x - shifts[0] < 2){
fprintf(stderr,"** WARNING: Edge comes too close to the ROI corners.\n");
// return 5;
}
/* Reference rows to the vertical centre of the data box */
j = size_y/2;
dt = shifts[j];
for (i = 0; i < size_y; i++) {
temp[i] = (double)i - (double)j;
shifts[i] -= dt;
}
*offset = dt;2、fit函数：根据上面locate_centroid函数的结果，利用最小二乘法进行拟合，得到一条拟合后的质心直线。
unsigned short fit(unsigned long ndata, double *x, double *y, double *b,
double *a, double *R2, double *avar, double *bvar)
{
unsigned long i;
double t,sxoss,syoss,sx=0.0,sy=0.0,st2=0.0;
double ss,sst,sigdat,chi2,siga,sigb;
*b=0.0;
for ( i=0; i < ndata; i++ ) {
sx += x[i];//x的叠加
sy += y[i];//y的叠加
}
//求平均值
ss=(double)ndata;
sxoss=sx/ss;
syoss=sy/ss;
for ( i=0; i < ndata; i++ ) {
t = x[i] - sxoss; //
st2 += t*t; //方差
*b += t * y[i];
}
*b /= st2; /* slope *///斜率
*a =(sy-sx*(*b))/ss; /* intercept *///截距
siga=sqrt((1.0+sx*sx/(ss*st2))/ss);
sigb=sqrt(1.0/st2);
chi2=0.0;
sst=0.0;
for (i=0; i < ndata; i++) {
chi2 += SQR( y[i] - (*a) - (*b) * x[i]);
sst += SQR( y[i] - syoss);
}
sigdat=sqrt(chi2/(ndata-2));
siga *= sigdat;
sigb *= sigdat;
*R2 = 1.0 - chi2/sst;//拟合程度
*avar = siga;
*bvar = sigb;
return 0;
}3、bin_to_regular_xgrid的作用：进行四倍的超采样得到ESF
unsigned short bin_to_regular_xgrid(unsigned short alpha,//alpha->指的是超取样的倍数
double *edgex, double *Signal,
double *AveEdge, long *counts,
unsigned short size_x,
unsigned short size_y)
{
long i, j, k,bin_number;
long bin_len;
int nzeros;
bin_len = size_x * alpha; //扩大四倍
for (i=0; i<bin_len; i++) {
AveEdge[i] = 0;
counts[i] = 0;
}
for (i=0; i<(size_x*(long)size_y); i++) {
bin_number = (long)floor((double)alpha*edgex[i]);//向下取整
if (bin_number >= 0) {
if (bin_number <= (bin_len - 1) ) {
AveEdge[bin_number] = AveEdge[bin_number] + Signal[i];//把每一行的距离边缘x轴一样远的信号相加
counts[bin_number] = (counts[bin_number])+1;//记录下对应位置有多少个信号相加
}
}
}
nzeros = 0;
for (i=0; i<bin_len; i++) {
j = 0;
k = 1;
//感觉写的有点复杂
if (counts[i] == 0) {
nzeros++; //记录有多少个位置是空的，即没有信号
//K的作用：因为这里的信号为0，找到后面第一个不为零的值，赋给当前这个零
if (i == 0) {
while (!j) { //当j==0时，表示此处的信号为0
if (counts[i+k] != 0) {//第一行的元素
AveEdge[i] = AveEdge[i+k]/((double) counts[i+k]);
j = 1;//充当flag。。。为啥不用布尔类型
}
else k++;//直到找到第一个不为零的数字才会停止
}
} else {
while (!j && ((i-k) >= 0) ) { //j==0&&i-k>=0 j==0说明 counts[i]是0 i-k>0说明 k在i前面，找前面不为零的数值赋给AveEdge[i]
if ( counts[i-k] != 0) {
AveEdge[i] = AveEdge[i-k]; /* Don't divide by counts since it already happened in previous iteration */
j = 1;
}
else k++;
}
if ( (i-k) < 0 ) {//k>i,其实联合上面那段，就是：此处的counts[i]累加次数为零，所以AveEdge[i]也就是0，所以要找到附近一个不为零的值赋给AveEdge[i]
k = 1;
while (!j) {
if (counts[i+k] != 0) {
AveEdge[i] = AveEdge[i+k]/((double) counts[i+k]);
j = 1;
}
else k++;
}
}
}
}
else
AveEdge[i] = (AveEdge[i])/ ((double) counts[i]);//如果此处不为零，直接就求个平均值
}
if (nzeros > 0) {//提示信息
fprintf(stderr, "\nWARNING: %d Zero counts found during projection binning.\n", nzeros);
fprintf(stderr, "The edge angle may be large, or you may need more lines of data.\n\n");
}
return nzeros;
}好了，今天先写到这，接下来有空会把接下来几个函数补上。
本文出自勇哥的网站《少有人走的路》wwww.skcircle.com，转载请注明出处！讨论可扫码加群：
