在前面的文章中，我们已经分析了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;
-
}
好了，今天先写到这，接下来有空会把接下来几个函数补上。
1636
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
