-
Notifications
You must be signed in to change notification settings - Fork 0
/
etfWeight.cu
executable file
·110 lines (84 loc) · 2.92 KB
/
etfWeight.cu
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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
#include <math.h>
const int LEN = 4;
// get a line centered at *ctrp, along the direction of *vec
__device__ void line2D(double out[][2*LEN+1], double *ctrp, double *vec) {
double x = ctrp[1];
double y = ctrp[0];
double gx = vec[0];
double gy = vec[1];
for (int i = -LEN; i <= LEN; i++){
out[0][i+LEN] = x+i*gx;
out[1][i+LEN] = y+i*gy;
}
}
__global__ void etfWeight(
double *xout, double *yout, double *outmag,
double *tx, double *ty,
double *im, double *gmag,
int height, int width) {
// calculate pixels' location
int x = (blockIdx.x * blockDim.x) + threadIdx.x;
int y = (blockIdx.y * blockDim.y) + threadIdx.y;
// Only execute valid pixels
if (x>=width || y>=height) {
return;
}
/*
// get parallel line centered at current point
double ctrp[2] = {y,x};
// direction vector
double vec[2] = {tx[y*width+x], ty[y*width+x]};
double p_line[2][2*LEN+1] = {0.0};
line2D(p_line, ctrp, vec);
double sum_inten = 0.0;
for(int i = 0; i < 2*LEN+1; i++) {
if (p_line[0][i]>=0 && p_line[0][i]<width &&
p_line[1][i]>=0 && p_line[1][i]<height){
int posx = (int)(p_line[0][i]+0.5);
int posy = (int)(p_line[1][i]+0.5);
int ind = posy*width+posx;
if (ind<width*height){
sum_inten += im[ind];
}
}
}
*/
// New try: trace forward and backward for points on the line defined
// by the ETF calculation.
double p_line[2][2*LEN+1] = {0.0};
double cur_ctrp[2] = {y, x};
double cur_pos = cur_ctrp[0]*width+cur_ctrp[1];
double cur_vec[2] = {tx[cur_pos], ty[cur_pos]}
for(int i = 1; i <= LEN; i++) {
int posx = (int)(cur_ctrp[1]+cur_vec[0]+0.5);
int posy = (int)(cur_ctrp[0]+cur_vec[1]+0.5);
p_line[0][LEN+i] = posx;
p_line[1][LEN+i] = posy;
// go to next point
cur_ctrp = {posy, posx};
cur_pos = cur_ctrp[0]*width+cur_ctrp[1];
cur_vec = {tx[cur_pos], ty[cur_pos]};
}
// reset starting point
cur_ctrp = {y, x};
cur_pos = cur_ctrp[0]*width+cur_ctrp[1];
cur_vec = {tx[cur_pos], ty[cur_pos]};
for(int i = -1; i >= -LEN; i--) {
int posx = (int)(cur_ctrp[1]-cur_vec[0]+0.5);
int posy = (int)(cur_ctrp[0]-cur_vec[1]+0.5);
p_line[0][LEN+i] = posx;
p_line[1][LEN+i] = posy;
// go to next point
cur_ctrp = {posy, posx};
cur_pos = cur_ctrp[0]*width+cur_ctrp[1];
cur_vec = {tx[cur_pos], ty[cur_pos]};
)
// for each point on the new line, calculate the orthogonal line
for(int i = 0; i < 2*LEN+1; i++) {
cur_ctrp = {p_line[1][i], ;
}
__syncthreads();
outmag[y*width+x] = sum_inten/(2*LEN+1);
xout[y*width+x] = tx[y*width+x];
yout[y*width+x] = ty[y*width+x];
}