func1_from_symmetric2.h
4.38 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
/// 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);
}