func1_from_symmetric2.h 4.38 KB
``````/// 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<typename T>
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<typename T>
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 = "<<X<<" must be odd."<<std::endl;
exit(1);
}
size_t C = X/2+1;												//calculate the center pixel coordinate
size_t N = C * C;												//number of values in the folded function

// The first step is to fold the function 8 times to take advantage of symmetry in the grid
T* folded = (T*) malloc(sizeof(T) * N );					//allocate space for the folded function
memset(folded, 0, sizeof(T) * N);
char* count = (char*) malloc( N );								//allocate space for a counter for the folded function
memset(count, 0, sizeof(T) * N);
size_t xi, yi;													//indices into the image f(xi, yi)
size_t xii, yii;												//indices into the folded image
T v;															//register to store the value at point (xi, yi)
for(xi = 0; xi < X; xi++){
for(yi = 0; yi < X; yi++){
v = fxy[yi * X + xi];									//retrieve f(x, y)

xii = xi;
yii = yi;												//initialize the indices into the folded image

//fold the function along the x and y axes
if(xi > 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<T>(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);
}``````