slm: OpenCL code base  0.1
countlink.cl
Go to the documentation of this file.
1 
10 #ifdef KERNEL_COUNT_DOWNCHANNELS
11 __kernel void count_downchannels(
35  __global const float2 *seed_point_array,
36  __global const bool *mask_array,
37  __global const float2 *uv_array,
38  __global uint *mapping_array,
39  __global uint *count_array,
40  __global uint *link_array
41  )
42 {
43  // For every channel head pixel...
44 
45  const uint global_id = get_global_id(0u)+get_global_id(1u)*get_global_size(0u);
46  __private uint idx, prev_idx, n_steps=0u, counter=1u;
47  __private float dl=0.0f, dt=DT_MAX;
48  __private float2 uv1_vec, uv2_vec, dxy1_vec, dxy2_vec,
49  vec = seed_point_array[global_id], next_vec;
50 
51  // Remember here
52  idx = get_array_idx(vec);
53  prev_idx = idx;
54  // Counter is already set=1 above
55  atomic_xchg(&count_array[idx],counter++);
56  atomic_or(&mapping_array[idx],IS_THINCHANNEL);
57  // Integrate downstream until the masked boundary is reached or n_steps too big
58  // HACK: factor 100x
59  while (!mask_array[idx] && n_steps<100*(MAX_N_STEPS-1)) {
60  compute_step_vec(dt, uv_array, &dxy1_vec, &dxy2_vec, &uv1_vec, &uv2_vec,
61  vec, &next_vec, &idx);
62  if (countlink_runge_kutta_step(&dt, &dl, &dxy1_vec, &dxy2_vec,
63  &vec, &next_vec, &idx, mapping_array)) {
64  break;
65  }
66  n_steps++;
67  // If at a new pixel
68  if (prev_idx!=idx) {
69  // Redesignate as a thin channel pixel
70  atomic_or(&mapping_array[idx],IS_THINCHANNEL);
71  // Link to here from the last pixel,
72  // i.e., point the previous pixel to this its downstream neighbor
73  atomic_xchg(&link_array[prev_idx],idx);
74  // If we've landed on a pixel whose channel length count
75  // exceeds our counter, we must have stepped off a minor onto a major
76  // channel, and thus need to stop
77  if (counter++<count_array[idx]) {
78  break;
79  }
80  prev_idx = idx;
81  }
82  }
83  if (!mask_array[prev_idx]) atomic_xchg(&link_array[prev_idx],idx);
84  return;
85 }
86 #endif
87 
88 #ifdef KERNEL_FLAG_DOWNCHANNELS
89 __kernel void flag_downchannels(
106  __global const float2 *seed_point_array,
107  __global const bool *mask_array,
108  __global const float2 *uv_array,
109  __global uint *mapping_array,
110  __global uint *count_array,
111  __global uint *link_array
112  )
113 {
114  // For every channel head pixel...
115 
116  const uint global_id = get_global_id(0u)+get_global_id(1u)*get_global_size(0u);
117  __private uint idx, prev_idx, counter=1u;
118  __private float2 vec = seed_point_array[global_id];
119 
120  // Remember here
121  idx = get_array_idx(vec);
122  prev_idx = idx+1;
123  // Counter=1 at channel head (set by count_downchannels)
124  // atomic_xchg(&count_array[idx],counter++);
125  atomic_or(&mapping_array[idx],IS_THINCHANNEL);
126  // Integrate downstream until the masked boundary is reached
127  while (!mask_array[idx] && prev_idx!=idx) {
128  prev_idx = idx;
129  idx = link_array[idx];
130  // Assume this idx is on the grid?
131  if (!mask_array[idx]) {
132  atomic_or(&mapping_array[idx],IS_THINCHANNEL);
133  atomic_max(&count_array[idx],counter++);
134  } else {
135  break;
136  }
137  }
138  // We have just stepped onto a masked pixel, so let's tag the previous pixel
139  // as a channel tail
140  if (!mask_array[prev_idx]) atomic_or(&mapping_array[prev_idx],IS_CHANNELTAIL);
141  return;
142 }
143 #endif
144 
145 #ifdef KERNEL_LINK_HILLSLOPES
146 __kernel void link_hillslopes(
167  __global const float2 *seed_point_array,
168  __global const bool *mask_array,
169  __global const float2 *uv_array,
170  __global uint *mapping_array,
171  __global uint *count_array,
172  __global uint *link_array
173  )
174 {
175  // For every hillslope aka non-thin-channel pixel...
176 
177  const uint global_id = get_global_id(0u)+get_global_id(1u)*get_global_size(0u);
178  __private uint idx, prev_idx, n_steps=0u;
179  __private float dl=0.0f, dt=DT_MAX;
180  __private float2 uv1_vec, uv2_vec, dxy1_vec, dxy2_vec,
181  vec = seed_point_array[global_id], next_vec;
182 
183  // Remember here
184  idx = get_array_idx(vec);
185  prev_idx = idx;
186  // Integrate downstream one pixel
187  // HACK: factor 100x
188  while (prev_idx==idx && !mask_array[idx] && n_steps<100*(MAX_N_STEPS-1)) {
189  prev_idx = idx;
190  compute_step_vec(dt, uv_array, &dxy1_vec, &dxy2_vec, &uv1_vec, &uv2_vec,
191  vec, &next_vec, &idx);
192  if (countlink_runge_kutta_step(&dt, &dl, &dxy1_vec, &dxy2_vec,
193  &vec, &next_vec, &idx, mapping_array)) {
194  printf("stuck\n");
195  break;
196  }
197  n_steps++;
198  }
199  // If we're on a new pixel, link previous pixel to here
200  if (prev_idx!=idx && !mask_array[idx]) {
201  atomic_xchg(&link_array[prev_idx],idx);
202  }
203  return;
204 }
205 #endif
#define DT_MAX
Definition: info.h:48
__kernel void count_downchannels(__global const float2 *seed_point_array, __global const bool *mask_array, __global const float2 *uv_array, __global uint *mapping_array, __global uint *count_array, __global uint *link_array)
REVISE? Integrate downstream from all channel heads until either a masked boundary pixel is reached o...
Definition: countlink.cl:34
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 IS_THINCHANNEL
Definition: info.h:97
#define MAX_N_STEPS
Definition: info.h:49
__kernel void link_hillslopes(__global const float2 *seed_point_array, __global const bool *mask_array, __global const float2 *uv_array, __global uint *mapping_array, __global uint *count_array, __global uint *link_array)
TBD.
Definition: countlink.cl:166
static bool countlink_runge_kutta_step(float *dt, float *dl, float2 *dxy1_vec, float2 *dxy2_vec, float2 *vec, float2 *next_vec, uint *idx, __global uint *mapping_array)
Compute a single step of 2nd-order Runge-Kutta numerical integration of a streamline given precompute...
#define IS_CHANNELTAIL
Definition: info.h:100