In the last image analysis force algorithm – SFR (Spatial Frequency Response) source analysis (A) introduced several important functions of SFR, next we look at the main process and other functions.
SfrProc is the main process function that calculates SFR values
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 ) {
// Position the torque center
err = locate_centroids(farea, temp, shifts, size_x, size_y, &offset1);
if(err ! =0) { return 2; }
/* Calculate the best fit line to the centroids */
/* Calculate the fitting curve of the moment center */
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 */
// Check the edge slope to ensure the quality of the following oversampling
err = check_slope(*slope, &size_y, numcycles, cyclelimit, 1);
/* Start image at new location, so that same row is center */
// Adjust the ROW of 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. */
// Update shifts based on the fitting results
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 */
// Do not normalize the data. The data will be normalized during DFT processing
/* To normalize here, set dt = min and dt1 = max of farea data */
// If you want to initialize here, set dt to minimum and dT1 to maximum
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];// Calculate the distance from each point to the knife edge
Signal[pcnt] = ((farea[((j*(long)size_x)+i)]) - dt)/(dt1-dt); // Normalize the brightness of each point
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 */
// Oversampling the current line image (ESF) of size_x*ALPHA(4) and saving it in 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. */
// Convert ESF to difference graph LSF, and calculate the moment center of 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 */
// Find the maximum value of LSF
locate_max_PSF( bin_len, AveEdge, &pcnt);// Here we get the center of the maximum
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)// Ignore the maximum value after difference
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. */
// Move the maximum value of the LSF to the center and add the Hamming window
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); // The frequency corresponding to each point
(*sfr)[i] = (double) (sfrc/AveTmp[0]); // Normalize SFR
}
/* Free */
free(AveEdge);
free(AveTmp);
*nrows = size_y;
*pcnt2 = pcnt;
return(0);
}
Copy the code5. Apply_hamming_window function: Apply hamming window to LSF to reduce noise
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 */
// Move the Avgedge to the center, and fill the space caused by the movement with 0
edge_offset = (*pcnt2) - (oldlen/2);// Pcnt2 must be less than oldlen/2.
if(edge_offset ! =0) { //cer=6
if (edge_offset < 0 ) { // Edge_offset should be less than 0. left
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]); / / write write write
for (i=0; i < -edge_offset; i++) // left center=3 right
AveEdge[i] = 0.00; /* last operation */ //offset = center-cer=-3
} else { // cer=6
/ / left
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]); / / write write write
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 */
// Fill the values for begin and end with 0, but feel useless
begin = (oldlen/2)-(newxwidth*alpha/2);// Begin will only equal 0, because 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;
// Add hamming window to the data between begin and end
/ / hamming window W (n, alpha) = (1 - alpha - alpha cos (2 * PI * n/(n - 1)), (0, n, n - 1 or less or less)
// Generally, the value of α is 0.46
// The window has been shifted (so the sign has changed) and the result is the same
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; }// Move lsfBEGIN left to index0
// This should not work in code,
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];
}
Copy the code6. Function of FTWOS: calculate DFT and convert airspace to frequency domain.
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 first half
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;
}
Copy the codeWell, the whole SFR important code notes to here also ended, drag quite a long time, mainly the previous period of time is also a little busy, we have what is not clear or error correction can leave a message. Let’s talk about it