Fraxinus  16.5.0-fx-rc1
An IGT application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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) 2008-2014, SINTEF Department of Medical Technology
5 All rights reserved.
6 
7 Redistribution and use in source and binary forms, with or without
8 modification, are permitted provided that the following conditions are met:
9 
10 1. Redistributions of source code must retain the above copyright notice,
11  this list of conditions and the following disclaimer.
12 
13 2. Redistributions in binary form must reproduce the above copyright notice,
14  this list of conditions and the following disclaimer in the documentation
15  and/or other materials provided with the distribution.
16 
17 3. Neither the name of the copyright holder nor the names of its contributors
18  may be used to endorse or promote products derived from this software
19  without specific prior written permission.
20 
21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
24 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
25 FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
27 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
29 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 =========================================================================*/
32 
33 #ifndef KERNELS_CLH_
34 #define KERNELS_CLH_
35 
40 /*******************/
41 /* Begin compile time definitions */
42 /*******************/
43 //MAX_PLANES - number of planes to include in closest planes
44 //N_PLANES - number of planes/images (the z dimension)
45 //METHOD - the reconstruction method
46 //PLANE_METHOD - the method used to find the closes planes
47 //MAX_MULTISTART_STARTS - how many guesses the algorithm should produce as starting points for finding the closes planes
48 //NEWNESS_FACTOR - Newness weight, should newer (pixels in newer planes) be preferred
49 //BRIGHTNESS_FACTOR - Brightness weight, when selecting pixrls, should similarity in intensity count
50 
51 /*******************/
52 /* End compile time definitions */
53 /*******************/
54 
55 /*******************/
56 /* Begin constants */
57 /*******************/
58 
59 #define CUBE_SIZE 4
60 
61 #define N_BLOCKS 10
62 
63 // Reconstruction methods
64 #define METHOD_VNN 0
65 #define METHOD_VNN2 1
66 #define METHOD_DW 2
67 #define METHOD_ANISOTROPIC 3
68 
69 // Plane searching methods
70 #define PLANE_HEURISTIC 0
71 #define PLANE_CLOSEST 1
72 
73 // Anisotropic method specific constants
74 #define ANISOTROPIC_SIZE 3
75 #define ANISOTROPIC_GAUSSIAN_SIZE 3
76 
77 #define ANISOTROPIC_METHOD_INCLUDE_ALL 0
78 #define ANISOTROPIC_METHOD_BILINEAR_ON_PLANE 1
79 #define ANISOTROPIC_METHOD_GAUSSIAN_ON_PLANE 2
80 
81 #define ANISOTROPIC_WEIGHT_METHOD_DISTANCE 0
82 #define ANISOTROPIC_WEIGHT_METHOD_BRIGHTNESS 1
83 #define ANISOTROPIC_WEIGHT_METHOD_LATENESS 2
84 #define ANISOTROPIC_WEIGHT_METHOD_BOTH 3
85 
86 #define ANISOTROPIC_METHOD ANISOTROPIC_METHOD_BILINEAR_ON_PLANE
87 
88 // Local search related parameters
89 
90 //#define LOCAL_SEARCH_DISTANCE (N_PLANES/(MAX_MULTISTART_STARTS))
91 //#define LOCAL_SEARCH_DISTANCE 10
92 /*****************/
93 /* End constants */
94 /*****************/
95 
96 /*****************/
97 /* Begin structs */
98 /*****************/
99 
100 typedef struct _close_plane
101 {
102  float dist;
103  short plane_id;
104  unsigned char intensity;
105  unsigned char padding; // Align with 4
106 } close_plane_t;
107 
108 typedef struct _output_volume_type
109 {
110  int3 size;
111  float3 spacing;
112  __global unsigned char* volume;
114 
115 /***************/
116 /* End structs */
117 /***************/
118 
119 /****************/
120 /* Begin macros */
121 /****************/
122 
123 //#define DEBUG
124 #define CHECK_PLANE_INDICES
125 
126 #ifdef DEBUG
127  #define DEBUG_PRINTF(...) if((get_global_id(0) % 5000) == 0) printf(##__VA_ARGS__)
128  //#define DEBUG_PRINTF(...) printf(##__VA_ARGS__)
129  //#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)
130  #define BOUNDS_CHECK(x, min, max)
131 #else
132 // #define DEBUG_PRINTF(...)
133  #define BOUNDS_CHECK(x, min, max)
134 #endif
135 
136 //#define PLANE_DIST(voxel, matrix) (dot(matrix.s26AE,voxel) - dot(matrix.s26AE, matrix.s37BF))
137 
138 #define EUCLID_DIST(a, b, c) sqrt((a)*(a) + (b)*(b) + (c)*(c))
139 
140 #define PROJECTONTOPLANE(voxel, matrix, dist) (voxel - dist*(matrix.s26AE))
141 
142 #define PROJECTONTOPLANEEQ(voxel, eq, dist) (voxel - dist*(eq))
143 
144 #define ISINSIDE(x, size) ((x) >= 0 && (x) < (size))
145 #define ISNOTMASKED(x, y, mask, xsize) ((mask)[(x) + (y)*(xsize)] > 0)
146 //#define ISNOTMASKED(x, y, mask, xsize) true
147 
148 #define VOXEL(v,x,y,z) v[x + y*volume_xsize + z*volume_ysize*volume_xsize]
149 
150 #define WEIGHT_INV(x) (1.0f/fabs(x))
151 #define WEIGHT_INV2(x) (1.0f/fabs(x*x))
152 #define WEIGHT_INV4(x) (1.0f/fabs(x*x*x*x))
153 #define WEIGHT_SUB(x) (1.0f - fabs(x))
154 
155 #define WEIGHT_TERNARY(val, mean, factor) \
156  ((val) >= (mean) ? (factor) : 0.0f)
157 
158 #define ANISOTROPIC_GAUSS_WEIGHT(px, var, mean, mean_id, sigma) WEIGHT_GAUSS(px.dist, sigma)
159 
160 #ifndef ANISOTROPIC_WEIGHT_METHOD
161  #define ANISOTROPIC_WEIGHT_METHOD ANISOTROPIC_WEIGHT_METHOD_BOTH
162 #endif
163 
164 #ifndef BRIGHTNESS_FACTOR
165  #define BRIGHTNESS_FACTOR 5.0f
166 #endif
167 
168 #ifndef NEWNESS_FACTOR
169  #define NEWNESS_FACTOR 5.0f
170 #endif
171 
172 #define ANISOTROPIC_WEIGHT_BRIGHTNESS(px, var, mean, mean_id, sigma) \
173  ((WEIGHT_GAUSS(px.dist, sigma)) + (WEIGHT_TERNARY(px.intensity, mean, BRIGHTNESS_FACTOR)))
174 
175 #define ANISOTROPIC_WEIGHT_LATENESS(px, var, mean, mean_id, sigma) \
176  ((WEIGHT_GAUSS(px.dist, sigma)) + (WEIGHT_TERNARY(px.plane_id, mean_id, NEWNESS_FACTOR)))
177 
178 #define ANISOTROPIC_WEIGHT_BOTH(px, var, mean, mean_id, sigma) \
179  ((WEIGHT_GAUSS(px.dist, sigma)) \
180  + (WEIGHT_TERNARY(px.plane_id, mean_id, NEWNESS_FACTOR)) \
181  + (WEIGHT_TERNARY(px.intensity, mean, BRIGHTNESS_FACTOR)))
182 
183 #if ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_DISTANCE
184  #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_GAUSS_WEIGHT(px, var, mean, mean_id, sigma)
185 #elif ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_BRIGHTNESS
186  #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_WEIGHT_BRIGHTNESS(px, var, mean, mean_id, sigma)
187 #elif ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_LATENESS
188  #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_WEIGHT_LATENESS(px, var, mean, mean_id, sigma)
189 #elif ANISOTROPIC_WEIGHT_METHOD == ANISOTROPIC_WEIGHT_METHOD_BOTH
190  #define ANISOTROPIC_WEIGHT(px, var, mean, mean_id, sigma) ANISOTROPIC_WEIGHT_BOTH(px, var, mean, mean_id, sigma)
191 #endif
192 
193 // Gaussian weight function
194 #define WEIGHT_GAUSS_SIGMA (0.05f)
195 
196 #define WEIGHT_GAUSS_SQRT_2PI 2.506628275f
197 
198 #define WEIGHT_GAUSS_NONEXP_PART(sigma) (1.0f/((sigma)*WEIGHT_GAUSS_SQRT_2PI))
199 #define WEIGHT_GAUSS_EXP_PART(dist, sigma) exp(-((dist)*(dist))/(2*(sigma)*(sigma)))
200 
201 #define WEIGHT_GAUSS(x, sigma) (WEIGHT_GAUSS_NONEXP_PART(sigma)*WEIGHT_GAUSS_EXP_PART(x, sigma))
202 
203 #define DW_WEIGHT(x) WEIGHT_INV(x)
204 #define VNN2_WEIGHT(x) WEIGHT_INV(x)
205 
206 #define CLOSE_PLANE_IDX(p, i) p[get_local_id(0)*(MAX_PLANES+1)+(i)]
207 
208 #define CREATE_OUTPUT_VOLUME_TYPE(name, in_size, in_spacing, in_volume) \
209  output_volume_type name; \
210  name.size = in_size; \
211  name.spacing = in_spacing; \
212  name.volume = in_volume;
213 
214 /**************/
215 /* End macros */
216 /**************/
217 
218 /********************/
219 /* Begin prototypes */
220 /********************/
221 
222 // Declare all the functions, as Apple seems to need that
223 //---------------------DEBUGGING-FUNCTIONALITY---------------------
224 #ifdef DEBUG
225  void printMatrix(float16 matrix);
226 #endif /* DEBUG */
227 
228 //---------------------DEBUGGING-FUNCTIONALITY---------------------
229 
230 int isValidPixel(int2 point, const __global unsigned char* mask, int2 in_size);
231 
232 int findHighestIdx(__local close_plane_t *planes, int n);
233 
234 int2 findClosestPlanes_heuristic(__local close_plane_t *close_planes,
235  __local float4* const plane_eqs,
236  __global float16* const plane_matrices,
237  const float4 voxel,
238  const float radius,
239  int guess,
240  bool doTermDistance,
241  __global const unsigned char* mask,
242  int2 in_size,
243  float2 in_spacing);
244 
245 int2 findClosestPlanes_multistart(__local close_plane_t *close_planes,
246  __local float4* const plane_eqs,
247  __global float16* const plane_matrices,
248  const float4 voxel,
249  const float radius,
250  int *multistart_guesses,
251  int n_multistart_guesses,
252  bool doTermDistance,
253  __global const unsigned char* mask,
254  int2 in_size,
255  float2 in_spacing);
256 
257 #if PLANE_METHOD == PLANE_EXACT
258  #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)
259 #elif PLANE_METHOD == PLANE_CLOSEST
260  #ifdef MAX_MULTISTART_STARTS
261  #undef MAX_MULTISTART_STARTS
262  #define MAX_MULTISTART_STARTS 1
263  #endif
264 
265  #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)
266 #endif
267 
268 __global const unsigned char* getImageData(int plane_id,
269  __global const unsigned char* bscans_blocks[],
270  int2 in_size);
271 
272 float4 transform(float16 matrix, float4 voxel);
273 
274 float4 transform_inv(float16 matrix, float4 voxel);
275 
276 float2 transform_inv_xy(float16 matrix, float4 voxel);
277 
278 int2 round_int(float2 value);
279 
280 int2 toImgCoord_int(float4 voxel, float16 plane_matrix, float2 in_spacing);
281 
282 float2 toImgCoord_float(float4 voxel, float16 plane_matrix, float2 in_spacing);
283 
284 float bilinearInterpolation(float x,
285  float y,
286  const __global unsigned char* image,
287  int in_xsize);
288 
289 #if METHOD == METHOD_VNN
290  #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i) \
291  performInterpolation_vnn(a, b, c, d, e, f, g, h, i)
292 #elif METHOD == METHOD_VNN2
293  #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i) \
294  performInterpolation_vnn2(a, b, c, d, e, f, g, h, i)
295 #elif METHOD == METHOD_DW
296  #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i) \
297  performInterpolation_dw(a, b, c, d, e, f, g, h, i)
298 #elif METHOD == METHOD_ANISOTROPIC
299  #define PERFORM_INTERPOLATION(a, b, c, d, e, f, g, h, i) \
300  performInterpolation_anisotropic(a, b, c, d, e, f, g, h, i)
301 #endif
302 
303 unsigned char performInterpolation_vnn(__local close_plane_t *close_planes,
304  int n_close_planes,
305  __global const float16 *plane_matrices,
306  __local const float4 *plane_eqs,
307  __global const unsigned char* bscans_blocks[],
308  int2 in_size,
309  float2 in_spacing,
310  __global const unsigned char* mask,
311  float4 voxel);
312 
313 unsigned char performInterpolation_vnn2(__local close_plane_t *close_planes,
314  int n_close_planes,
315  __global const float16 *plane_matrices,
316  __local const float4 *plane_eqs,
317  __global const unsigned char* bscans_blocks[],
318  int2 in_size,
319  float2 in_spacing,
320  __global const unsigned char* mask,
321  float4 voxel);
322 
323 unsigned char performInterpolation_dw(__local close_plane_t *close_planes,
324  int n_close_planes,
325  __global const float16 *plane_matrices,
326  __local const float4 *plane_eqs,
327  __global const unsigned char* bscans_blocks[],
328  int2 in_size,
329  float2 in_spacing,
330  __global const unsigned char* mask,
331  float4 voxel);
332 
333 unsigned char performInterpolation_anisotropic(__local close_plane_t *close_planes,
334  int n_close_planes,
335  __global const float16 *plane_matrices,
336  __local const float4 *plane_eqs,
337  __global const unsigned char* bscans_blocks[],
338  int2 in_size,
339  float2 in_spacing,
340  __global const unsigned char* mask,
341  float4 voxel);
342 
343 unsigned char anisotropicFilter(__local const close_plane_t *pixels,
344  int n_planes);
345 
346 void prepare_plane_eqs(__global float16 *plane_matrices,
347  __local float4 *plane_eqs);
348 
349 int findLocalMinimas(int *guesses,
350  __local const float4 *plane_eqs,
351  float radius,
352  float4 voxel,
353  float3 out_spacing,
354  __global const float16 *plane_matrices,
355  __global const unsigned char *mask);
356 
357 __kernel void voxel_methods(int volume_xsize,
358  int volume_ysize,
359  int volume_zsize,
360  float volume_xspacing,
361  float volume_yspacing,
362  float volume_zspacing,
363  int in_xsize,
364  int in_ysize,
365  float in_xspacing,
366  float in_yspacing,
367  // TODO: Wouldn't it be kind of nice if the bscans was an image sampler object?
368  __global unsigned char* in_bscans_b0,
369  __global unsigned char* in_bscans_b1,
370  __global unsigned char* in_bscans_b2,
371  __global unsigned char* in_bscans_b3,
372  __global unsigned char* in_bscans_b4,
373  __global unsigned char* in_bscans_b5,
374  __global unsigned char* in_bscans_b6,
375  __global unsigned char* in_bscans_b7,
376  __global unsigned char* in_bscans_b8,
377  __global unsigned char* in_bscans_b9,
378  __global unsigned char* out_volume,
379  __global float16 *plane_matrices,
380  __global unsigned char* mask,
381  __local float4 *plane_eqs,
382  __local close_plane_t *planes,
383  float radius);
384 
385 /******************/
386 /* End prototypes */
387 /******************/
388 
389 #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:104
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)
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:103
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:105
float2 transform_inv_xy(float16 matrix, float4 voxel)
__global unsigned char * volume
Definition: kernels.cl.h:112
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)