/// 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); }