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