-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathkernel_3d.cl
121 lines (113 loc) · 4.32 KB
/
kernel_3d.cl
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
111
112
113
114
115
116
117
118
119
120
__kernel void SFEGO_3d_gradient(__global float *data, __global float *diff,
__global float *direct_x, __global float *direct_y, __global float *direct_z,
__global int *list_x, __global int *list_y, __global int *list_z,
__global float *list_unit_x, __global float *list_unit_y, __global float *list_unit_z,
__global int *surface_indexs,
__global int *hemisphere_dp_pos_add, __global int *hemisphere_dp_pos_sub,
__global int *hemisphere_dp_neg_add, __global int *hemisphere_dp_neg_sub,
__global int *dp_pos_add_start_idxs, __global int *dp_pos_sub_start_idxs,
__global int *dp_neg_add_start_idxs, __global int *dp_neg_sub_start_idxs,
__global int *dp_pos_add_start_lens, __global int *dp_pos_sub_start_lens,
__global int *dp_neg_add_start_lens, __global int *dp_neg_sub_start_lens,
int list_len, int dp_len, int dim_x, int dim_y, int dim_z
) {
//kernel index
int x=get_global_id(0); //x
int y=get_global_id(1); //y
int z=get_global_id(2); //z
//np_data[z*dim_y*dim_x+y*dim_x+x]
if(x<dim_x && y<dim_y && z<dim_z){
//init variable
float pos_sum=0.0, neg_sum=0.0;
int pos_count=0, neg_count=0;
float avg_pos, avg_neg, avg_diff;
float max_diff=-99999.9;
int max_index;
float current=data[z*dim_y*dim_x+y*dim_x+x];
for(int i=0;i<dp_len;++i){
for(int j=0;j<dp_pos_add_start_lens[i];++j){
int idx=j+dp_pos_add_start_idxs[i];
int index=hemisphere_dp_pos_add[idx];
int tx=x+list_x[index];
int ty=y+list_y[index];
int tz=z+list_z[index];
if( (tx>=0 && tx<dim_x) && (ty>=0 && ty<dim_y) && (tz>=0 && tz<dim_z)){
pos_sum+=data[tz*dim_y*dim_x+ty*dim_x+tx];
pos_count++;
}
}
for(int j=0;j<dp_pos_sub_start_lens[i];++j){
int idx=j+dp_pos_sub_start_idxs[i];
int index=hemisphere_dp_pos_sub[idx];
int tx=x+list_x[index];
int ty=y+list_y[index];
int tz=z+list_z[index];
if( (tx>=0 && tx<dim_x) && (ty>=0 && ty<dim_y) && (tz>=0 && tz<dim_z)){
pos_sum-=data[tz*dim_y*dim_x+ty*dim_x+tx];
pos_count--;
}
}
for(int j=0;j<dp_neg_add_start_lens[i];++j){
int idx=j+dp_neg_add_start_idxs[i];
int index=hemisphere_dp_neg_add[idx];
int tx=x+list_x[index];
int ty=y+list_y[index];
int tz=z+list_z[index];
if( (tx>=0 && tx<dim_x) && (ty>=0 && ty<dim_y) && (tz>=0 && tz<dim_z)){
neg_sum+=data[tz*dim_y*dim_x+ty*dim_x+tx];
neg_count++;
}
}
for(int j=0;j<dp_neg_sub_start_lens[i];++j){
int idx=j+dp_neg_sub_start_idxs[i];
int index=hemisphere_dp_neg_sub[idx];
int tx=x+list_x[index];
int ty=y+list_y[index];
int tz=z+list_z[index];
if( (tx>=0 && tx<dim_x) && (ty>=0 && ty<dim_y) && (tz>=0 && tz<dim_z)){
neg_sum-=data[tz*dim_y*dim_x+ty*dim_x+tx];
neg_count--;
}
}
if(pos_count>0) avg_pos=pos_sum/pos_count;
else avg_pos=current;
if(neg_count>0) avg_neg=neg_sum/neg_count;
else avg_neg=current;
avg_diff=avg_pos-avg_neg;
if(avg_diff>max_diff){
max_diff=avg_diff;
max_index=surface_indexs[i];
}
}
//printf("%d %d %d => %d %f\n", x, y, z, max_index, max_diff);
direct_x[z*dim_y*dim_x+y*dim_x+x]=list_unit_x[max_index];
direct_y[z*dim_y*dim_x+y*dim_x+x]=list_unit_y[max_index];
direct_z[z*dim_y*dim_x+y*dim_x+x]=list_unit_z[max_index];
diff[z*dim_y*dim_x+y*dim_x+x]=max_diff;
}
}
__kernel void SFEGO_3d_integral(__global float *result, __global float *diff,
__global float *direct_x, __global float *direct_y, __global float *direct_z,
__global int *list_x, __global int *list_y, __global float *list_z,
__global float *list_unit_x, __global float *list_unit_y, __global float *list_unit_z,
int list_len, int dim_x, int dim_y, int dim_z) {
//kernel index
int x=get_global_id(0); //x
int y=get_global_id(1); //y
int z=get_global_id(2); //z
if(x<dim_x && y<dim_y && z<dim_z){
int tx,ty,tz;
result[z*dim_y*dim_x+y*dim_x+x]=0.0;
for(int i=0;i<list_len;++i){
tx=x+list_x[i];
ty=y+list_y[i];
tz=z+list_z[i];
if( (tx>=0 && tx<dim_x) && (ty>=0 && ty<dim_y) && (tz>=0 && tz<dim_z) ){
float dot_product= (list_unit_x[i]*direct_x[tz*dim_y*dim_x+ty*dim_x+tx])+
(list_unit_y[i]*direct_y[tz*dim_y*dim_x+ty*dim_x+tx])+
(list_unit_z[i]*direct_z[tz*dim_y*dim_x+ty*dim_x+tx]);
result[z*dim_y*dim_x+y*dim_x+x]+=dot_product*diff[tz*dim_y*dim_x+ty*dim_x+tx];
}
}
}
}