Fraxinus  2023.01.05-dev+develop.0da12
An IGT application
kernels.cl.h
Go to the documentation of this file.
1 /*=========================================================================
2 This file is part of CustusX, an Image Guided Therapy Application.
3 
4 Copyright (c) SINTEF Department of Medical Technology.
5 All rights reserved.
6 
7 CustusX is released under a BSD 3-Clause license.
8 
9 See Lisence.txt (https://github.com/SINTEFMedtek/CustusX/blob/master/License.txt) for details.
10 =========================================================================*/
11 
12 #ifndef KERNELS_CLH_
13 #define KERNELS_CLH_
14 
19 /*******************/
20 /* Begin compile time definitions */
21 /*******************/
22 //MAX_PLANES - number of planes to include in closest planes
23 //N_PLANES - number of planes/images (the z dimension)
24 //METHOD - the reconstruction method
25 //PLANE_METHOD - the method used to find the closes planes
26 //MAX_MULTISTART_STARTS - how many guesses the algorithm should produce as starting points for finding the closes planes
27 //NEWNESS_FACTOR - Newness weight, should newer (pixels in newer planes) be preferred
28 //BRIGHTNESS_FACTOR - Brightness weight, when selecting pixrls, should similarity in intensity count
29 
30 /*******************/
31 /* End compile time definitions */
32 /*******************/
33 
34 /*******************/
35 /* Begin constants */
36 /*******************/
37 
38 #define CUBE_SIZE 4
39 
40 #define N_BLOCKS 10
41 
42 // Reconstruction methods
43 #define METHOD_VNN 0
44 #define METHOD_VNN2 1
45 #define METHOD_DW 2
46 #define METHOD_ANISOTROPIC 3
47 
48 // Plane searching methods
49 #define PLANE_HEURISTIC 0
50 #define PLANE_CLOSEST 1
51 
52 // Anisotropic method specific constants
53 #define ANISOTROPIC_SIZE 3
54 #define ANISOTROPIC_GAUSSIAN_SIZE 3
55 
56 #define ANISOTROPIC_METHOD_INCLUDE_ALL 0
57 #define ANISOTROPIC_METHOD_BILINEAR_ON_PLANE 1
58 #define ANISOTROPIC_METHOD_GAUSSIAN_ON_PLANE 2
59 
60 #define ANISOTROPIC_WEIGHT_METHOD_DISTANCE 0
61 #define ANISOTROPIC_WEIGHT_METHOD_BRIGHTNESS 1
62 #define ANISOTROPIC_WEIGHT_METHOD_LATENESS 2
63 #define ANISOTROPIC_WEIGHT_METHOD_BOTH 3
64 
65 #define ANISOTROPIC_METHOD ANISOTROPIC_METHOD_BILINEAR_ON_PLANE
66 
67 // Local search related parameters
68 
69 //#define LOCAL_SEARCH_DISTANCE (N_PLANES/(MAX_MULTISTART_STARTS))
70 //#define LOCAL_SEARCH_DISTANCE 10
71 /*****************/
72 /* End constants */
73 /*****************/
74 
75 /*****************/
76 /* Begin structs */
77 /*****************/
78 
79 typedef struct _close_plane
80 {
81  float dist;
82  short plane_id;
83  unsigned char intensity;
84  unsigned char padding; // Align with 4
86 
87 typedef struct _output_volume_type
88 {
89  int3 size;
90  float3 spacing;
91  __global unsigned char* volume;
93 
94 /***************/
95 /* End structs */
96 /***************/
97 
98 /****************/
99 /* Begin macros */
100 /****************/
101 
102 //#define DEBUG
103 #define CHECK_PLANE_INDICES
104 
105 #ifdef DEBUG
106  #define DEBUG_PRINTF(...) if((get_global_id(0) % 5000) == 0) printf(##__VA_ARGS__)
107  //#define DEBUG_PRINTF(...) printf(##__VA_ARGS__)
108  //#define BOUNDS_CHECK(x, min, ma x) if(x < min || x >= max) printf("Line %d: %s out of range: %d min: %d max: %d\n", __LINE__, #x, x, min, max)
109  #define BOUNDS_CHECK(x, min, max)
110 #else
111 // #define DEBUG_PRINTF(...)
112  #define BOUNDS_CHECK(x, min, max)
113 #endif
114 
115 //#define PLANE_DIST(voxel, matrix) (dot(matrix.s26AE,voxel) - dot(matrix.s26AE, matrix.s37BF))
116 
117 #define EUCLID_DIST(a, b, c) sqrt((a)*(a) + (b)*(b) + (c)*(c))
118 
119 #define PROJECTONTOPLANE(voxel, matrix, dist) (voxel - dist*(matrix.s26AE))
120 
121 #define PROJECTONTOPLANEEQ(voxel, eq, dist) (voxel - dist*(eq))
122 
123 #define ISINSIDE(x, size) ((x) >= 0 && (x) < (size))
124 #define ISNOTMASKED(x, y, mask, xsize) ((mask)[(x) + (y)*(xsize)] > 0)
125 //#define ISNOTMASKED(x, y, mask, xsize) true
126 
127 #define VOXEL(v,x,y,z) v[x + y*volume_xsize + z*volume_ysize*volume_xsize]
128 
129 #define WEIGHT_INV(x) (1.0f/fabs(x))
130 #define WEIGHT_INV2(x) (1.0f/fabs(x*x))
131 #define WEIGHT_INV4(x) (1.0f/fabs(x*x*x*x))
132 #define WEIGHT_SUB(x) (1.0f - fabs(x))
133 
134 #define WEIGHT_TERNARY(val, mean, factor) \
135  ((val) >= (mean) ? (factor) : 0.0f)
136 
137 #define ANISOTROPIC_GAUSS_WEIGHT(px, var, mean, mean_id, sigma) WEIGHT_GAUSS(px.dist, sigma)
138 
139 #ifndef ANISOTROPIC_WEIGHT_METHOD
140  #define ANISOTROPIC_WEIGHT_METHOD ANISOTROPIC_WEIGHT_METHOD_BOTH
141 #endif
142 
143 #ifndef BRIGHTNESS_FACTOR
144  #define BRIGHTNESS_FACTOR 5.0f
145 #endif
146 
147 #ifndef NEWNESS_FACTOR
148  #define NEWNESS_FACTOR 5.0f
149 #endif
150 
151 #define ANISOTROPIC_WEIGHT_BRIGHTNESS(px, var, mean, mean_id, sigma) \
152  ((WEIGHT_GAUSS(px.dist, sigma)) + (WEIGHT_TERNARY(px.intensity, mean, BRIGHTNESS_FACTOR)))
153 
154 #define ANISOTROPIC_WEIGHT_LATENESS(px, var, mean, mean_id, sigma) \
155  ((WEIGHT_GAUSS(px.dist, sigma)) + (WEIGHT_TERNARY(px.plane_id, mean_id, NEWNESS_FACTOR)))
156 
157 #define ANISOTROPIC_WEIGHT_BOTH(px, var, mean, mean_id, sigma) \
158  ((WEIGHT_GAUSS(px.dist, sigma)) \
159  + (WEIGHT_TERNARY(px.plane_id, mean_id, NEWNESS_FACTOR)) \
160  + (WEIGHT_TERNARY(px.intensity, mean, BRIGHTNESS_FACTOR)))
161 
162 #if ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_DISTANCE
163  #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_GAUSS_WEIGHT(px, var, mean, mean_id, sigma)
164 #elif ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_BRIGHTNESS
165  #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_WEIGHT_BRIGHTNESS(px, var, mean, mean_id, sigma)
166 #elif ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_LATENESS
167  #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_WEIGHT_LATENESS(px, var, mean, mean_id, sigma)
168 #elif ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_BOTH
169  #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_WEIGHT_BOTH(px, var, mean, mean_id, sigma)
170 #endif
171 
172 // Gaussian weight function
173 #define WEIGHT_GAUSS_SIGMA (0.05f)
174 
175 #define WEIGHT_GAUSS_SQRT_2PI 2.506628275f
176 
177 #define WEIGHT_GAUSS_NONEXP_PART(sigma) (1.0f/((sigma)*WEIGHT_GAUSS_SQRT_2PI))
178 #define WEIGHT_GAUSS_EXP_PART(dist, sigma) exp(-((dist)*(dist))/(2*(sigma)*(sigma)))
179 
180 #define WEIGHT_GAUSS(x, sigma) (WEIGHT_GAUSS_NONEXP_PART(sigma)*WEIGHT_GAUSS_EXP_PART(x, sigma))
181 
182 #define DW_WEIGHT(x) WEIGHT_INV(x)
183 #define VNN2_WEIGHT(x) WEIGHT_INV(x)
184 
185 #define CLOSE_PLANE_IDX(p, i) p[get_local_id(0)*(MAX_PLANES+1)+(i)]
186 
187 #define CREATE_OUTPUT_VOLUME_TYPE(name, in_size, in_spacing, in_volume) \
188  output_volume_type name; \
189  name.size = in_size; \
190  name.spacing = in_spacing; \
191  name.volume = in_volume;
192 
193 /**************/
194 /* End macros */
195 /**************/
196 
197 /********************/
198 /* Begin prototypes */
199 /********************/
200 
201 // Declare all the functions, as Apple seems to need that
202 //---------------------DEBUGGING-FUNCTIONALITY---------------------
203 #ifdef DEBUG
204  void printMatrix(float16 matrix);
205 #endif /* DEBUG */
206 
207 //---------------------DEBUGGING-FUNCTIONALITY---------------------
208 
209 int isValidPixel(int2 point, const __global unsigned char* mask, int2 in_size);
210 
211 int findHighestIdx(__local close_plane_t *planes, int n);
212 
213 int2 findClosestPlanes_heuristic(__local close_plane_t *close_planes,
214  __local float4* const plane_eqs,
215  __global float16* const plane_matrices,
216  const float4 voxel,
217  const float radius,
218  int guess,
219  bool doTermDistance,
220  __global const unsigned char* mask,
221  int2 in_size,
222  float2 in_spacing);
223 
224 int2 findClosestPlanes_multistart(__local close_plane_t *close_planes,
225  __local float4* const plane_eqs,
226  __global float16* const plane_matrices,
227  const float4 voxel,
228  const float radius,
229  int *multistart_guesses,
230  int n_multistart_guesses,
231  bool doTermDistance,
232  __global const unsigned char* mask,
233  int2 in_size,
234  float2 in_spacing);
235 
236 #if PLANE_METHOD == PLANE_EXACT
237  #define FIND_CLOSE_PLANES(a, b, c, d, e, f, g, h, i, j) findClosestPlanes_multistart(a, b, c, d, e, f, g, 1, h, i, j)
238 #elif PLANE_METHOD == PLANE_CLOSEST
239  #ifdef MAX_MULTISTART_STARTS
240  #undef MAX_MULTISTART_STARTS
241  #define MAX_MULTISTART_STARTS 1
242  #endif
243 
244  #define FIND_CLOSE_PLANES(a, b, c, d, e, f, g, h, i, j) findClosestPlanes_multistart(a, b, c, d, e, f, g, 0, h, i, j)
245 #endif
246 
247 __global const unsigned char* getImageData(int plane_id,
248  __global const unsigned char* bscans_blocks[],
249  int2 in_size);
250 
251 float4 transform(float16 matrix, float4 voxel);
252 
253 float4 transform_inv(float16 matrix, float4 voxel);
254 
255 float2 transform_inv_xy(float16 matrix, float4 voxel);
256 
257 int2 round_int(float2 value);
258 
259 int2 toImgCoord_int(float4 voxel, float16 plane_matrix, float2 in_spacing);
260 
261 float2 toImgCoord_float(float4 voxel, float16 plane_matrix, float2 in_spacing);
262 
263 float bilinearInterpolation(float x,
264  float y,
265  const __global unsigned char* image,
266  int in_xsize);
267 
268 #if METHOD == METHOD_VNN
269  #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i) \
270  performInterpolation_vnn(a, b, c, d, e, f, g, h, i)
271 #elif METHOD == METHOD_VNN2
272  #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i) \
273  performInterpolation_vnn2(a, b, c, d, e, f, g, h, i)
274 #elif METHOD == METHOD_DW
275  #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i) \
276  performInterpolation_dw(a, b, c, d, e, f, g, h, i)
277 #elif METHOD == METHOD_ANISOTROPIC
278  #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i) \
279  performInterpolation_anisotropic(a, b, c, d, e, f, g, h, i)
280 #endif
281 
282 unsigned char performInterpolation_vnn(__local close_plane_t *close_planes,
283  int n_close_planes,
284  __global const float16 *plane_matrices,
285  __local const float4 *plane_eqs,
286  __global const unsigned char* bscans_blocks[],
287  int2 in_size,
288  float2 in_spacing,
289  __global const unsigned char* mask,
290  float4 voxel);
291 
292 unsigned char performInterpolation_vnn2(__local close_plane_t *close_planes,
293  int n_close_planes,
294  __global const float16 *plane_matrices,
295  __local const float4 *plane_eqs,
296  __global const unsigned char* bscans_blocks[],
297  int2 in_size,
298  float2 in_spacing,
299  __global const unsigned char* mask,
300  float4 voxel);
301 
302 unsigned char performInterpolation_dw(__local close_plane_t *close_planes,
303  int n_close_planes,
304  __global const float16 *plane_matrices,
305  __local const float4 *plane_eqs,
306  __global const unsigned char* bscans_blocks[],
307  int2 in_size,
308  float2 in_spacing,
309  __global const unsigned char* mask,
310  float4 voxel);
311 
312 unsigned char performInterpolation_anisotropic(__local close_plane_t *close_planes,
313  int n_close_planes,
314  __global const float16 *plane_matrices,
315  __local const float4 *plane_eqs,
316  __global const unsigned char* bscans_blocks[],
317  int2 in_size,
318  float2 in_spacing,
319  __global const unsigned char* mask,
320  float4 voxel);
321 
322 unsigned char anisotropicFilter(__local const close_plane_t *pixels,
323  int n_planes);
324 
325 void prepare_plane_eqs(__global float16 *plane_matrices,
326  __local float4 *plane_eqs);
327 
328 int findLocalMinimas(int *guesses,
329  __local const float4 *plane_eqs,
330  float radius,
331  float4 voxel,
332  float3 out_spacing,
333  __global const float16 *plane_matrices,
334  __global const unsigned char *mask);
335 
336 __kernel void voxel_methods(int volume_xsize,
337  int volume_ysize,
338  int volume_zsize,
339  float volume_xspacing,
340  float volume_yspacing,
341  float volume_zspacing,
342  int in_xsize,
343  int in_ysize,
344  float in_xspacing,
345  float in_yspacing,
346  // TODO: Wouldn't it be kind of nice if the bscans was an image sampler object?
347  __global unsigned char* in_bscans_b0,
348  __global unsigned char* in_bscans_b1,
349  __global unsigned char* in_bscans_b2,
350  __global unsigned char* in_bscans_b3,
351  __global unsigned char* in_bscans_b4,
352  __global unsigned char* in_bscans_b5,
353  __global unsigned char* in_bscans_b6,
354  __global unsigned char* in_bscans_b7,
355  __global unsigned char* in_bscans_b8,
356  __global unsigned char* in_bscans_b9,
357  __global unsigned char* out_volume,
358  __global float16 *plane_matrices,
359  __global unsigned char* mask,
360  __local float4 *plane_eqs,
361  __local close_plane_t *planes,
362  float radius);
363 
364 /******************/
365 /* End prototypes */
366 /******************/
367 
368 #endif /* KERNELS_CLH_ */
unsigned char performInterpolation_anisotropic(__local close_plane_t *close_planes, int n_close_planes, __global const float16 *plane_matrices, __local const float4 *plane_eqs, __global const unsigned char *bscans_blocks[], int2 in_size, float2 in_spacing, __global const unsigned char *mask, float4 voxel)
unsigned char intensity
Definition: kernels.cl.h:83
int isValidPixel(int2 point, const __global unsigned char *mask, int2 in_size)
int2 round_int(float2 value)
__global const unsigned char * getImageData(int plane_id, __global const unsigned char *bscans_blocks[], int2 in_size)
float2 toImgCoord_float(float4 voxel, float16 plane_matrix, float2 in_spacing)
int findHighestIdx(__local close_plane_t *planes, int n)
float dist
Definition: kernels.cl.h:81
void prepare_plane_eqs(__global float16 *plane_matrices, __local float4 *plane_eqs)
int2 findClosestPlanes_heuristic(__local close_plane_t *close_planes, __local float4 *const plane_eqs, __global float16 *const plane_matrices, const float4 voxel, const float radius, int guess, bool doTermDistance, __global const unsigned char *mask, int2 in_size, float2 in_spacing)
__kernel void voxel_methods(int volume_xsize, int volume_ysize, int volume_zsize, float volume_xspacing, float volume_yspacing, float volume_zspacing, int in_xsize, int in_ysize, float in_xspacing, float in_yspacing, __global unsigned char *in_bscans_b0, __global unsigned char *in_bscans_b1, __global unsigned char *in_bscans_b2, __global unsigned char *in_bscans_b3, __global unsigned char *in_bscans_b4, __global unsigned char *in_bscans_b5, __global unsigned char *in_bscans_b6, __global unsigned char *in_bscans_b7, __global unsigned char *in_bscans_b8, __global unsigned char *in_bscans_b9, __global unsigned char *out_volume, __global float16 *plane_matrices, __global unsigned char *mask, __local float4 *plane_eqs, __local close_plane_t *planes, float radius)
short plane_id
Definition: kernels.cl.h:82
int2 toImgCoord_int(float4 voxel, float16 plane_matrix, float2 in_spacing)
unsigned char anisotropicFilter(__local const close_plane_t *pixels, int n_planes)
unsigned char performInterpolation_dw(__local close_plane_t *close_planes, int n_close_planes, __global const float16 *plane_matrices, __local const float4 *plane_eqs, __global const unsigned char *bscans_blocks[], int2 in_size, float2 in_spacing, __global const unsigned char *mask, float4 voxel)
int findLocalMinimas(int *guesses, __local const float4 *plane_eqs, float radius, float4 voxel, float3 out_spacing, __global const float16 *plane_matrices, __global const unsigned char *mask)
unsigned char performInterpolation_vnn(__local close_plane_t *close_planes, int n_close_planes, __global const float16 *plane_matrices, __local const float4 *plane_eqs, __global const unsigned char *bscans_blocks[], int2 in_size, float2 in_spacing, __global const unsigned char *mask, float4 voxel)
int2 findClosestPlanes_multistart(__local close_plane_t *close_planes, __local float4 *const plane_eqs, __global float16 *const plane_matrices, const float4 voxel, const float radius, int *multistart_guesses, int n_multistart_guesses, bool doTermDistance, __global const unsigned char *mask, int2 in_size, float2 in_spacing)
struct _close_plane close_plane_t
float4 transform_inv(float16 matrix, float4 voxel)
float bilinearInterpolation(float x, float y, const __global unsigned char *image, int in_xsize)
struct _output_volume_type output_volume_type
unsigned char padding
Definition: kernels.cl.h:84
float2 transform_inv_xy(float16 matrix, float4 voxel)
__global unsigned char * volume
Definition: kernels.cl.h:91
unsigned char performInterpolation_vnn2(__local close_plane_t *close_planes, int n_close_planes, __global const float16 *plane_matrices, __local const float4 *plane_eqs, __global const unsigned char *bscans_blocks[], int2 in_size, float2 in_spacing, __global const unsigned char *mask, float4 voxel)
float4 transform(float16 matrix, float4 voxel)