少有人走的路

在前面的文章中，我们已经分析了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，转载请注明出处！讨论可扫码加群：