In the previous article, we have analyzed the algorithm principle and steps of SFR, let’s directly analyze the source code, the main function in the source code is mainly divided into the following several:

1. Locate_centroids is used to locate the center of moment of each row of pixels

unsigned short locate_centroids(double *farea, double *temp, double *shifts,unsigned short size_x, unsigned short size_y,double *offset)
Copy the code
 unsigned long i, j;
  double dt, dt1, dt2;

  /* Compute the first difference on each line. Interpolate to find the centroid of the first derivatives. */
  // Compute difference compute moment center position
  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;
  }
// Prevent the center of moment from getting too close to the image boundary
  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;
Copy the code

2. Fit function: According to the results of locate_centroid function above, the least square method is used to fit and a fitted centroid line is obtained.

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 of the stack
    sy += y[i];/ / y of the stack
  }
  // Find the average value
  ss=(double)ndata;
  sxoss=sx/ss;   
  syoss=sy/ss;
  for ( i=0; i < ndata; i++ ) {
    t = x[i] - sxoss;  //
    st2 += t*t; / / variance
    *b += t * y[i];  
  }
  *b /= st2;         /* slope *// / the slope
  *a =(sy-sx*(*b))/ss; /* intercept *// / 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/(ndata2 -));
  siga *= sigdat;
  sigb *= sigdat;
  *R2 = 1.0 - chi2/sst;// The degree of fitting
  *avar = siga;
  *bvar = sigb;
  return 0;
}
Copy the code

3. Functions of Bin_TO_REGULar_Xgrid: ESF is obtained by quadruple oversampling

unsigned short bin_to_regular_xgrid(unsigned short alpha,Alpha -> refers to the multiple of the oversampling
				    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;  // Quadruple

  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]);// round down
    if (bin_number >= 0) {
      if (bin_number <= (bin_len - 1) ) {
        AveEdge[bin_number] = AveEdge[bin_number] + Signal[i];// Add the signals in each row that are as far away from the edge X-axis
        counts[bin_number] = (counts[bin_number])+1;// Record how many signals are added to the corresponding position
      }
    }
  }

  nzeros = 0;
  for (i=0; i<bin_len; i++) {
    j = 0;
    k = 1;
    // It's a bit complicated
    if (counts[i] == 0) {
      nzeros++; // Record how many positions are empty, i.e. there is no signal
      Since the signal here is 0, find the first non-zero value and assign it to the current zero
      if (i == 0) {
        while(! j) {// if j==0, the signal is 0
          if(counts[i+k] ! =0) {// The element in the first line
	    AveEdge[i] = AveEdge[i+k]/((double) counts[i+k]);
            j = 1;// act as flag... Why not use a Boolean type
	  }
          else k++;// Stop until the first non-zero number is found}}else {
        while(! j && ((i-k) >=0)) {//j==0&& I -k>=0 j==0 counts[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 ) {// count [I] counts = 0 so AveEdge[I] counts = 0
	  k = 1;
	  while(! j) {if(counts[i+k] ! =0) {
	      AveEdge[i] = AveEdge[i+k]/((double) counts[i+k]);
	      j = 1;
	    } 
	    elsek++; }}}}else 
      AveEdge[i] = (AveEdge[i])/ ((double) counts[i]);// If this is not zero, just find the average value
  }

  if (nzeros > 0) {// Prompt message
    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;
}
Copy the code

Okay, so I’m going to write that today, and then I’ll fill in the next few functions when I have time.

 

Welcome to pay attention to my personal public number, there are more good health oh ~