slm: OpenCL code base  0.1
channelheads.cl
Go to the documentation of this file.
1 
14 
15 #ifdef KERNEL_MAP_CHANNEL_HEADS
16 __kernel void map_channel_heads(
34  __global const float2 *seed_point_array,
35  __global const bool *mask_array,
36  __global const float2 *uv_array,
37  __global uint *mapping_array
38  )
39 {
40  // For every non-masked pixel...
41 
42  const uint global_id = get_global_id(0u)+get_global_id(1u)*get_global_size(0u);
43  __private uint idx, prev_idx, n_steps = 0u;
44  __private float dl = 0.0f, dt = DT_MAX;
45  __private float2 uv1_vec, uv2_vec, dxy1_vec, dxy2_vec,
46  vec = seed_point_array[global_id], next_vec;
47 
48  // Remember here
49  idx = get_array_idx(vec);
50  prev_idx = idx;
51  // Integrate downstream one pixel
52  while (prev_idx==idx && !mask_array[idx] && n_steps<(MAX_N_STEPS-1)) {
53  compute_step_vec(dt, uv_array, &dxy1_vec, &dxy2_vec, &uv1_vec, &uv2_vec,
54  vec, &next_vec, &idx);
55  channelheads_runge_kutta_step(&dt, &dl, &dxy1_vec, &dxy2_vec,
56  &vec, &next_vec, &n_steps, &idx);
57  }
58  // If need be, integrate further downstream until a IS_THINCHANNEL pixel is reached
59  n_steps = 0u;
60  while (!mask_array[idx] && !(mapping_array[idx] & IS_THINCHANNEL)
61  && n_steps<(MAX_N_STEPS-1)) {
62  compute_step_vec(dt, uv_array, &dxy1_vec, &dxy2_vec, &uv1_vec, &uv2_vec,
63  vec, &next_vec, &idx);
64  channelheads_runge_kutta_step(&dt, &dl, &dxy1_vec, &dxy2_vec,
65  &vec, &next_vec, &n_steps, &idx);
66  }
67  if (n_steps>=(MAX_N_STEPS-1)) {
68 // printf("stuck...");
69  return;
70  }
71  // Unset the channel head flag unless we're really at a channel head maybe
72  if (!mask_array[idx]) {
73  idx = get_array_idx(vec);
74  // If here is not channel...
75  if ((~mapping_array[idx]) & IS_THINCHANNEL){
76  // ...flag here as not channel head
77  atomic_and(&mapping_array[idx],~IS_CHANNELHEAD);
78  } else {
79  // Here is a channel
80  // If previous pixel was channel...
81  if (mapping_array[prev_idx] & IS_THINCHANNEL) {
82  // ...flag here as not channel head
83  atomic_and(&mapping_array[idx],~IS_CHANNELHEAD);
84  }
85  }
86  }
87  return;
88 }
89 #endif
90 
91 #ifdef KERNEL_PRUNE_CHANNEL_HEADS
92 
93 // Check if this nbr is a thin channel pixel and not masked
94 // If so, add one to the 'flag'.
95 // Add 16 if it's masked, thus recording if *any* nbr is masked.
96 #define CHECK_IS_THINCHANNEL(idx) ((mapping_array[idx] & IS_THINCHANNEL)>0)
97 #define CHECK_IS_MASKED(idx) (mask_array[idx])
98 #define CHECK_THINCHANNEL(nbr_vec) { \
99  idx = get_array_idx(nbr_vec); \
100  flag += (CHECK_IS_THINCHANNEL(idx) | CHECK_IS_MASKED(idx)*16); \
101  }
102 // Check all eight pixel-nbr directions
103 #define CHECK_E(vec) CHECK_THINCHANNEL((float2)( vec[0]+1.0f, vec[1] ))
104 #define CHECK_NE(vec) CHECK_THINCHANNEL((float2)( vec[0]+1.0f, vec[1]+1.0f ))
105 #define CHECK_N(vec) CHECK_THINCHANNEL((float2)( vec[0] , vec[1]+1.0f ))
106 #define CHECK_NW(vec) CHECK_THINCHANNEL((float2)( vec[0]-1.0f, vec[1]+1.0f ))
107 #define CHECK_W(vec) CHECK_THINCHANNEL((float2)( vec[0]-1.0f, vec[1] ))
108 #define CHECK_SW(vec) CHECK_THINCHANNEL((float2)( vec[0]-1.0f, vec[1]-1.0f ))
109 #define CHECK_S(vec) CHECK_THINCHANNEL((float2)( vec[0] , vec[1]-1.0f ))
110 #define CHECK_SE(vec) CHECK_THINCHANNEL((float2)( vec[0]+1.0f, vec[1]-1.0f ))
111 
130 __kernel void prune_channel_heads(
131  __global const float2 *seed_point_array,
132  __global const bool *mask_array,
133  __global const float2 *uv_array,
134  __global uint *mapping_array
135  )
136 {
137  // For every provisional channel head pixel...
138 
139  const uint global_id = get_global_id(0u)+get_global_id(1u)*get_global_size(0u);
140  // Use an unsigned integer for the 'flag' var, because it needs it to be big
141  __private uint idx, flag = 0;
142  __private float2 vec = seed_point_array[global_id];
143 
144  CHECK_N(vec);
145  CHECK_S(vec);
146  CHECK_E(vec);
147  CHECK_W(vec);
148  CHECK_NE(vec);
149  CHECK_SE(vec);
150  CHECK_NW(vec);
151  CHECK_SW(vec);
152  // If flag==1, one and only one nbr is a thin channel pixel
153  // Otherwise, remove this provisional channel head.
154  if (flag!=1) {
155  idx = get_array_idx(vec);
156  atomic_and(&mapping_array[idx],~IS_CHANNELHEAD);
157  // If there are no thin channel neighbors AT ALL,
158  // we must be at an isolated pixel.
159  // Thus redesignate this pixel as 'not channelized at all'.
160  if (flag==0 || flag>=16) {
161  atomic_and(&mapping_array[idx],IS_THINCHANNEL);
162  atomic_and(&mapping_array[idx],IS_CHANNEL);
163  }
164  }
165  return;
166 }
167 #endif
#define DT_MAX
Definition: info.h:48
#define IS_CHANNELHEAD
Definition: info.h:99
__kernel void prune_channel_heads(__global const float2 *seed_point_array, __global const bool *mask_array, __global const float2 *uv_array, __global uint *mapping_array)
Keep only those provisional channel heads that lie on the &#39;thin channel&#39; skeletonized network and hav...
__kernel void map_channel_heads(__global const float2 *seed_point_array, __global const bool *mask_array, __global const float2 *uv_array, __global uint *mapping_array)
Map provisional channel heads, even including those not on an IS_THINCHANNEL pixel and thus extraneou...
Definition: channelheads.cl:33
#define CHECK_E(vec)
#define CHECK_SW(vec)
#define CHECK_N(vec)
#define CHECK_S(vec)
static void compute_step_vec(const float dt, const __global float2 *uv_array, float2 *dxy1_vec, float2 *dxy2_vec, float2 *uv1_vec, float2 *uv2_vec, const float2 vec, float2 *next_vec, uint *idx)
Compute a 2nd-order Runge-Kutta integration step along a streamline.
Definition: computestep.cl:41
static uint get_array_idx(float2 vec)
Compute the array index of the padded grid pixel pointed to by a float2 grid position vector (choice ...
Definition: essentials.cl:62
#define CHECK_NE(vec)
#define CHECK_NW(vec)
#define IS_THINCHANNEL
Definition: info.h:97
#define CHECK_W(vec)
#define MAX_N_STEPS
Definition: info.h:49
#define CHECK_SE(vec)
#define IS_CHANNEL
Definition: info.h:96
static void channelheads_runge_kutta_step(float *dt, float *dl, float2 *dxy1_vec, float2 *dxy2_vec, float2 *vec, float2 *next_vec, uint *n_steps, uint *idx)
Compute a single step of 2nd-order Runge-Kutta numerical integration of a streamline given precompute...