/// Reconstruct a 1D function from a 2D symmetric function. This function takes a 2D image f(x,y) as input and
/// builds a 1D function f(r) where r = sqrt(x^2 + y^2) to approximate this 2D function.
/// This is useful for several applications, such as:
/// 1) Calculating a 1D function from a noisy 2D image, when you know the 2D image is supposed to be symmetric
/// 2) Calculating the average value for every r = sqrt(x^2 + y^2)
/// Given a set of function samples equally spaced by dx, calculate the two samples closest to x and the proximity ratio alpha.
/// This can be used to linearly interpolate between an array of equally spaced values. Given the query value x, the
/// interpolated value can be calculated as r = values[sample] * alpha + values[sample + 1] * (1 - alpha)
/// @param sample is the lowest bin closest to the query point x
/// @param alpha is the ratio of x between [sample, sample + 1]
/// @param dx is the spacing between values
/// @param x is the query point
template
void lerp_alpha(T& sample, T& alpha, T dx, T x){
sample = std::floor(x/dx);
alpha = 1 - (x - (b * dx)) / dx;
}
/// This function assumes that the input image is square, that the # of samples are odd, and that r=0 is at the center
/// @param fr is an array of X elements that will store the reconstructed function
/// @param dr is the spacing (in pixels) between samples in fr
template
void cpu_func1_from_symmetric2(T* fr, T& dr, T* fxy, size_t X){
if(X%2 == 0){ //the 2D function must be odd (a sample must be available for r=0)
std::err<<"Error, X = "< C) xii = 2 * C - xi - 1; //calculate the folded index of x
if(yi > C) yii = 2 * C - yi - 1; //calculate the folded index of y
if(xii < yii) std::swap(xii, yii); //fold the function again along the 45-degree line
folded[yii * C + xii] += v; //add the value to the folded function
count[yii * C + xii] += 1; //add a counter to the counter table
}
}
//divide out the counter to correct the folded function
for(size_t i = 0; i < N){
folded[i] /= (T)count[i]; //divide out the counter
}
T max_r = sqrt(X * X + Y * Y); //calculate the maximum r value, which will be along the image diagonal
T dr = max_r / (X - 1); //spacing between samples in the output function f(r)
T* fA = (T*) malloc( sizeof(T) * X); //allocate space for a counter function storing alpha weights
memset(fA, 0, sizeof(T) * X); //zero out the alpha array
memset(fr, 0, sizeof(T) * X); //zero out the output function
T r; //register to store the value of r at each point
size_t sample;
T alpha;
for(xi = 0; xi < C; xi++){
for(yi = 0; yi < xi; yi++){
r = sqrt(xi*xi + yi*yi); //calculate the value of r for the current (x, y)
lerp_alpha(sample, alpha, dr, r); //calculate the lowest nearby sample index and the associated alpha weight
fr[sample] += folded[yi * C + xi] * alpha; //sum the weighted value from the folded function
fA[sample] += alpha; //sum the weight
if(sample < X - 1){ //if we aren't dealing with the last bin
fr[sample + 1] += folded[yi * C + xi] * (1.0 - alpha); //calculate the weighted value for the second point
fA[sample + 1] += 1 - alpha; //add the second alpha value
}
}
}
//divide out the alpha values
for(size_t i = 0; i < X; i++)
fr[i] /= fA[i];
//free allocated memory
free(folded);
free(count);
free(fA);
}