67 #define METHOD_ANISOTROPIC 3
70 #define PLANE_HEURISTIC 0
71 #define PLANE_CLOSEST 1
74 #define ANISOTROPIC_SIZE 3
75 #define ANISOTROPIC_GAUSSIAN_SIZE 3
77 #define ANISOTROPIC_METHOD_INCLUDE_ALL 0
78 #define ANISOTROPIC_METHOD_BILINEAR_ON_PLANE 1
79 #define ANISOTROPIC_METHOD_GAUSSIAN_ON_PLANE 2
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
86 #define ANISOTROPIC_METHOD ANISOTROPIC_METHOD_BILINEAR_ON_PLANE
124 #define CHECK_PLANE_INDICES
127 #define DEBUG_PRINTF(...) if((get_global_id(0) % 5000) == 0) printf(##__VA_ARGS__)
130 #define BOUNDS_CHECK(x, min, max)
133 #define BOUNDS_CHECK(x, min, max)
138 #define EUCLID_DIST(a, b, c) sqrt((a)*(a) + (b)*(b) + (c)*(c))
140 #define PROJECTONTOPLANE(voxel, matrix, dist) (voxel - dist*(matrix.s26AE))
142 #define PROJECTONTOPLANEEQ(voxel, eq, dist) (voxel - dist*(eq))
144 #define ISINSIDE(x, size) ((x) >= 0 && (x) < (size))
145 #define ISNOTMASKED(x, y, mask, xsize) ((mask)[(x) + (y)*(xsize)] > 0)
148 #define VOXEL(v,x,y,z) v[x + y*volume_xsize + z*volume_ysize*volume_xsize]
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))
155 #define WEIGHT_TERNARY(val, mean, factor) \
156 ((val) >= (mean) ? (factor) : 0.0f)
158 #define ANISOTROPIC_GAUSS_WEIGHT(px, var, mean, mean_id, sigma) WEIGHT_GAUSS(px.dist, sigma)
160 #ifndef ANISOTROPIC_WEIGHT_METHOD
161 #define ANISOTROPIC_WEIGHT_METHOD ANISOTROPIC_WEIGHT_METHOD_BOTH
164 #ifndef BRIGHTNESS_FACTOR
165 #define BRIGHTNESS_FACTOR 5.0f
168 #ifndef NEWNESS_FACTOR
169 #define NEWNESS_FACTOR 5.0f
172 #define ANISOTROPIC_WEIGHT_BRIGHTNESS(px, var, mean, mean_id, sigma) \
173 ((WEIGHT_GAUSS(px.dist, sigma)) + (WEIGHT_TERNARY(px.intensity, mean, BRIGHTNESS_FACTOR)))
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)))
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)))
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)
194 #define WEIGHT_GAUSS_SIGMA (0.05f)
196 #define WEIGHT_GAUSS_SQRT_2PI 2.506628275f
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)))
201 #define WEIGHT_GAUSS(x, sigma) (WEIGHT_GAUSS_NONEXP_PART(sigma)*WEIGHT_GAUSS_EXP_PART(x, sigma))
203 #define DW_WEIGHT(x) WEIGHT_INV(x)
204 #define VNN2_WEIGHT(x) WEIGHT_INV(x)
206 #define CLOSE_PLANE_IDX(p, i) p[get_local_id(0)*(MAX_PLANES+1)+(i)]
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;
225 void printMatrix(float16 matrix);
230 int isValidPixel(int2 point,
const __global
unsigned char* mask, int2 in_size);
235 __local float4*
const plane_eqs,
236 __global float16*
const plane_matrices,
241 __global
const unsigned char* mask,
246 __local float4*
const plane_eqs,
247 __global float16*
const plane_matrices,
250 int *multistart_guesses,
251 int n_multistart_guesses,
253 __global
const unsigned char* mask,
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
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)
268 __global
const unsigned char*
getImageData(
int plane_id,
269 __global
const unsigned char* bscans_blocks[],
272 float4
transform(float16 matrix, float4 voxel);
280 int2
toImgCoord_int(float4 voxel, float16 plane_matrix, float2 in_spacing);
282 float2
toImgCoord_float(float4 voxel, float16 plane_matrix, float2 in_spacing);
286 const __global
unsigned char* image,
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)
305 __global
const float16 *plane_matrices,
306 __local
const float4 *plane_eqs,
307 __global
const unsigned char* bscans_blocks[],
310 __global
const unsigned char* mask,
315 __global
const float16 *plane_matrices,
316 __local
const float4 *plane_eqs,
317 __global
const unsigned char* bscans_blocks[],
320 __global
const unsigned char* mask,
325 __global
const float16 *plane_matrices,
326 __local
const float4 *plane_eqs,
327 __global
const unsigned char* bscans_blocks[],
330 __global
const unsigned char* mask,
335 __global
const float16 *plane_matrices,
336 __local
const float4 *plane_eqs,
337 __global
const unsigned char* bscans_blocks[],
340 __global
const unsigned char* mask,
347 __local float4 *plane_eqs);
350 __local
const float4 *plane_eqs,
354 __global
const float16 *plane_matrices,
355 __global
const unsigned char *mask);
360 float volume_xspacing,
361 float volume_yspacing,
362 float volume_zspacing,
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,
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)
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)
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
float2 transform_inv_xy(float16 matrix, float4 voxel)
__global unsigned char * volume
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)