|
12 | 12 | #define SIZE ((SHAPE_Z * SLICE)) |
13 | 13 |
|
14 | 14 |
|
15 | | -kernel |
16 | | -void rotate_grid3d( |
17 | | - global float *grid, float16 rotmat, global float *out, int nearest |
18 | | - ) |
19 | | -{ |
20 | | - // Rotate grid around the origin. Only grid points within LLENGTH of the |
21 | | - // origin are rotated. |
22 | | - |
23 | | - int zid = get_global_id(0); |
24 | | - int yid = get_global_id(1); |
25 | | - int xid = get_global_id(2); |
26 | | - int zstride = get_global_size(0); |
27 | | - int ystride = get_global_size(1); |
28 | | - int xstride = get_global_size(2); |
29 | | - |
30 | | - int z, y, x, x0, y0, z0, x1, y1, z1, offset0, offset1, grid_ind; |
31 | | - float dx, dy, dz, dx1, dy1, dz1, c00, c10, c01, c11, c0, c1, c; |
32 | | - float4 dist2, coor_z, coor_zy, coor_zyx; |
33 | | - int4 out_ind; |
34 | | - |
35 | | - for (z = zid - LLENGTH; z <= LLENGTH; z += zstride) { |
36 | | - dist2.s2 = SQUARE(z); |
37 | | - |
38 | | - coor_z.s0 = rotmat.s6 * z; |
39 | | - coor_z.s1 = rotmat.s7 * z; |
40 | | - coor_z.s2 = rotmat.s8 * z; |
41 | | - |
42 | | - out_ind.s0 = z * SLICE; |
43 | | - if (z < 0) |
44 | | - out_ind.s0 += SIZE; |
45 | | - |
46 | | - for (y = yid - LLENGTH; y <= LLENGTH; y += ystride) { |
47 | | - dist2.s1 = SQUARE(y) + dist2.s2; |
48 | | - if (dist2.s1 > LLENGTH2) |
49 | | - continue; |
50 | | - |
51 | | - coor_zy.s0 = rotmat.s3 * y + coor_z.s0; |
52 | | - coor_zy.s1 = rotmat.s4 * y + coor_z.s1; |
53 | | - coor_zy.s2 = rotmat.s5 * y + coor_z.s2; |
54 | | - |
55 | | - out_ind.s1 = out_ind.s0 + y * SHAPE_X; |
56 | | - if (y < 0) |
57 | | - out_ind.s1 += SLICE; |
58 | | - |
59 | | - for (x = xid - LLENGTH; x <= LLENGTH; x += xstride) { |
60 | | - dist2.s0 = SQUARE(x) + dist2.s1; |
61 | | - if (dist2.s0 > LLENGTH2) |
62 | | - continue; |
63 | | - coor_zyx.s0 = rotmat.s0 * x + coor_zy.s0; |
64 | | - coor_zyx.s1 = rotmat.s1 * x + coor_zy.s1; |
65 | | - coor_zyx.s2 = rotmat.s2 * x + coor_zy.s2; |
66 | | - |
67 | | - out_ind.s2 = out_ind.s1 + x; |
68 | | - if (x < 0) |
69 | | - out_ind.s2 += SHAPE_X; |
70 | | - |
71 | | - if (nearest > 0) { |
72 | | - |
73 | | - x0 = (int) round(coor_zyx.s0); |
74 | | - y0 = (int) round(coor_zyx.s1); |
75 | | - z0 = (int) round(coor_zyx.s2); |
76 | | - |
77 | | - grid_ind = z0 * SLICE + y0 * SHAPE_X + x0; |
78 | | - if (x0 < 0) |
79 | | - grid_ind += SHAPE_X; |
80 | | - if (y0 < 0) |
81 | | - grid_ind += SLICE; |
82 | | - if (z0 < 0) |
83 | | - grid_ind += SIZE; |
84 | | - |
85 | | - out[out_ind.s2] = grid[grid_ind]; |
86 | | - |
87 | | - } else { |
88 | | - x0 = (int) floor(coor_zyx.s0); |
89 | | - y0 = (int) floor(coor_zyx.s1); |
90 | | - z0 = (int) floor(coor_zyx.s2); |
91 | | - x1 = x0 + 1; |
92 | | - y1 = y0 + 1; |
93 | | - z1 = z0 + 1; |
94 | | - |
95 | | - // Grid index |
96 | | - grid_ind = z0 * SLICE + y0 * SHAPE_X + x0; |
97 | | - if (x0 < 0) |
98 | | - grid_ind += SHAPE_X; |
99 | | - if (y0 < 0) |
100 | | - grid_ind += SLICE; |
101 | | - if (z0 < 0) |
102 | | - grid_ind += SIZE; |
103 | | - |
104 | | - dx = coor_zyx.s0 - x0; |
105 | | - dy = coor_zyx.s1 - y0; |
106 | | - dz = coor_zyx.s2 - z0; |
107 | | - dx1 = 1 - dx; |
108 | | - dy1 = 1 - dx; |
109 | | - dz1 = 1 - dx; |
110 | | - |
111 | | - offset1 = 1; |
112 | | - if (x1 == 0) |
113 | | - offset1 -= SHAPE_X; |
114 | | - c00 = grid[grid_ind] * dx1 + |
115 | | - grid[grid_ind + offset1] * dx; |
116 | | - |
117 | | - offset0 = SHAPE_X; |
118 | | - if (y1 == 0) |
119 | | - offset0 -= SLICE; |
120 | | - offset1 = offset0 + 1; |
121 | | - if (x1 == 0) |
122 | | - offset1 -= SHAPE_X; |
123 | | - c10 = grid[grid_ind + offset0] * dx1 + |
124 | | - grid[grid_ind + offset1] * dx; |
125 | | - |
126 | | - offset0 = SLICE; |
127 | | - if (z1 == 0) |
128 | | - offset0 -= SIZE; |
129 | | - offset1 = offset0 + 1; |
130 | | - if (x1 == 0) |
131 | | - offset1 -= SHAPE_X; |
132 | | - c01 = grid[grid_ind + offset0] * dx1 + |
133 | | - grid[grid_ind + offset1] * dx; |
134 | | - |
135 | | - offset0 = SLICE + SHAPE_X; |
136 | | - if (z1 == 0) |
137 | | - offset0 -= SIZE; |
138 | | - if (y1 == 0) |
139 | | - offset0 -= SLICE; |
140 | | - offset1 = offset0 + 1; |
141 | | - if (x1 == 0) |
142 | | - offset1 -= SHAPE_X; |
143 | | - c01 = grid[grid_ind + offset0] * dx1 + |
144 | | - grid[grid_ind + offset1] * dx; |
145 | | - |
146 | | - c0 = c00 * dy1 + c10 * dy; |
147 | | - c1 = c01 * dy1 + c11 * dy; |
148 | | - // TODO fix why nan = 1 * 1 + 0 * 0, instead of 1 |
149 | | - printf("c1 = c01 * dy1 + c11 * dy: %f = %f * %f + %f * %f\n", c1, c01, dy1, c11, dy); |
150 | | - |
151 | | - c = c0 * dz1 + c1 * dz; |
152 | | - |
153 | | - //printf("c = c0 * dz1 + c1 * dz: %f = %f * %f + %f * %f\n", c, c0, dz1, c1, dz); |
154 | | - out[out_ind.s2] = c; |
155 | | - } |
156 | | - } |
157 | | - } |
158 | | - } |
159 | | -} |
160 | | - |
161 | | - |
162 | 15 | kernel |
163 | 16 | void rotate_image3d( |
164 | 17 | read_only image3d_t image, sampler_t sampler, float16 rotmat, |
|
0 commit comments