slm: OpenCL code base  0.1
integrationfns.cl
Go to the documentation of this file.
1 
14 
15 #ifdef KERNEL_INTEGRATE_TRAJECTORY
16 static inline bool runge_kutta_step_record(float *dt, float *dl, float *l_trajectory,
50  float2 *dxy1_vec, float2 *dxy2_vec,
51  float2 *vec, float2 *prev_vec,
52  float2 *next_vec,
53  uint *n_steps, uint *idx, uint *prev_idx,
54  __global char2 *trajectory_vec)
55 {
56  const float step_error = fast_length((*dxy2_vec-*dxy1_vec)/GRID_SCALE);
57 
58  *dl = fast_length(*dxy2_vec);
59  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
60  *next_vec = *vec+*dxy1_vec;
61  *dl = fast_length(*dxy1_vec);
62  }
63  *vec = *next_vec;
64  *idx = get_array_idx(*next_vec);
65  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
66  update_record_trajectory(*dl,l_trajectory,*vec,*prev_vec,n_steps,trajectory_vec);
67  return true;
68  }
69  *dt = select( fmin(DT_MAX,(ADJUSTED_MAX_ERROR*(*dt))/(step_error)), DT_MAX,
70  isequal(step_error,0.0f) );
71  update_record_trajectory(*dl,l_trajectory,*vec,*prev_vec,n_steps,trajectory_vec);
72  *prev_vec = *vec;
73 // *prev_idx = select(*prev_idx,*idx,isnotequal(*idx,*prev_idx));
74  *prev_idx = *idx;
75  return false;
76 }
77 #endif
78 
79 #ifdef KERNEL_INTEGRATE_TRAJECTORY
80 static inline void euler_step_record(float *dt, float *dl,
104  float *l_trajectory, const float2 uv_vec,
105  float2 *vec, const float2 prev_vec,
106  uint *n_steps, __global char2 *trajectory_vec)
107 {
108  const float2 sgnd_uv_vec = uv_vec*(float2)(DOWNUP_SIGN,DOWNUP_SIGN);
109  const float dt_x = dt_to_nearest_edge((*vec)[0], sgnd_uv_vec[0]);
110  const float dt_y = dt_to_nearest_edge((*vec)[1], sgnd_uv_vec[1]);
111 
112  *dt = minmag(dt_x,dt_y);
113  *vec += sgnd_uv_vec*(*dt);
114  *vec = approximate(*vec);
115  *vec = (float2)( fmin(fmax((*vec)[0],-0.5f),X_MAX),
116  fmin(fmax((*vec)[1],-0.5f),Y_MAX) );
117  *dl = fast_length(*vec-prev_vec);
118  update_record_trajectory(*dl,l_trajectory,*vec,prev_vec,n_steps,trajectory_vec);
119  }
120 #endif
121 
122 #ifdef KERNEL_INTEGRATE_TRAJECTORY
123 static inline bool runge_kutta_step_write_sl_data(
163  float *dt, float *dl, float *l_trajectory,
164  float2 *dxy1_vec, float2 *dxy2_vec,
165  float2 *vec, float2 *prev_vec, const float2 next_vec,
166  uint *n_steps, uint *idx, uint *prev_idx,
167  __global const bool *mask_array,
168  __global uint *slt_array, __global uint *slc_array)
169 {
170  const float step_error = fast_length((*dxy2_vec-*dxy1_vec)/GRID_SCALE);
171 
172  *dl = fast_length(*dxy2_vec);
173  *vec = next_vec;
174  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
175  update_trajectory_write_sl_data(*dl,l_trajectory,*vec,*prev_vec,n_steps,
176  idx, prev_idx, mask_array, slt_array, slc_array);
177  return true;
178  }
179  *dt = select( fmin(DT_MAX,(ADJUSTED_MAX_ERROR*(*dt))/(step_error)), DT_MAX,
180  isequal(step_error,0.0f) );
181  update_trajectory_write_sl_data(*dl,l_trajectory,*vec,*prev_vec,n_steps,
182  idx, prev_idx, mask_array, slt_array, slc_array);
183  *prev_vec = *vec;
184  return false;
185 }
186 #endif
187 
188 #ifdef KERNEL_INTEGRATE_TRAJECTORY
189 static inline void euler_step_write_sl_data(float *dt, float *dl,
220  float *l_trajectory,
221  const float2 uv_vec, float2 *vec,
222  const float2 prev_vec,
223  uint *n_steps,
224  uint *idx, uint *prev_idx,
225  __global const bool *mask_array,
226  __global uint *slt_array,
227  __global uint *slc_array)
228 {
229  __private float2 sgnd_uv_vec;
230  __private float dt_x, dt_y;
231 
232  sgnd_uv_vec = uv_vec*(float2)(DOWNUP_SIGN,DOWNUP_SIGN);
233  dt_x = dt_to_nearest_edge((*vec)[0], sgnd_uv_vec[0]);
234  dt_y = dt_to_nearest_edge((*vec)[1], sgnd_uv_vec[1]);
235  *dt = minmag(dt_x,dt_y);
236  *vec += sgnd_uv_vec*(*dt);
237  *vec = approximate(*vec);
238  *vec = (float2)( fmin(fmax((*vec)[0],-0.5f),X_MAX),
239  fmin(fmax((*vec)[1],-0.5f),Y_MAX) );
240  *dl = fast_length(*vec-prev_vec);
241  update_trajectory_write_sl_data(*dl,l_trajectory,*vec,prev_vec,n_steps,
242  idx, prev_idx, mask_array, slt_array, slc_array);
243  }
244 #endif
245 
246 #ifdef KERNEL_CONNECT_CHANNELS
247 static inline bool connect_runge_kutta_step_record(float *dt, float *dl,
281  float *l_trajectory,
282  float2 *dxy1_vec, float2 *dxy2_vec,
283  float2 *vec, float2 *prev_vec,
284  float2 *next_vec,
285  uint *n_steps,
286  uint *idx, uint *prev_idx,
287  __private char2 *trajectory_vec)
288 {
289  const float step_error = fast_length((*dxy2_vec-*dxy1_vec)/GRID_SCALE);
290 
291  *dl = fast_length(*dxy2_vec);
292  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
293  *next_vec = *vec+*dxy1_vec;
294  *dl = fast_length(*dxy1_vec);
295  }
296  *vec = *next_vec;
297  *idx = get_array_idx(*next_vec);
298  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
299  update_record_private_trajectory(*dl,l_trajectory,*vec,*prev_vec,n_steps,
300  trajectory_vec);
301  return true;
302  }
303  *dt = select( fmin(DT_MAX,(ADJUSTED_MAX_ERROR*(*dt))/(step_error)), DT_MAX,
304  isequal(step_error,0.0f) );
305  update_record_private_trajectory(*dl,l_trajectory,*vec,*prev_vec,n_steps,
306  trajectory_vec);
307  *prev_vec = *vec;
308  return false;
309 }
310 #endif
311 
312 #ifdef KERNEL_MAP_CHANNEL_HEADS
313 static inline void channelheads_runge_kutta_step(float *dt, float *dl,
336  float2 *dxy1_vec, float2 *dxy2_vec,
337  float2 *vec, float2 *next_vec,
338  uint *n_steps, uint *idx)
339 {
340  const float step_error = fast_length((*dxy2_vec-*dxy1_vec)/GRID_SCALE);
341 
342  *dl = fast_length(*dxy2_vec);
343  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
344  *next_vec = *vec+*dxy1_vec;
345  *dl = fast_length(*dxy1_vec);
346  }
347  *vec = *next_vec;
348  *idx = get_array_idx(*next_vec);
349  *n_steps += 1u;
350  *dt = select( fmin(DT_MAX,(ADJUSTED_MAX_ERROR*(*dt))/(step_error)), DT_MAX,
351  isequal(step_error,0.0f) );
352  return;
353 }
354 #endif
355 
356 #if defined(KERNEL_COUNT_DOWNCHANNELS) || defined(KERNEL_LINK_HILLSLOPES)
357 static inline bool countlink_runge_kutta_step(float *dt, float *dl,
383  float2 *dxy1_vec, float2 *dxy2_vec,
384  float2 *vec, float2 *next_vec, uint *idx,
385  __global uint *mapping_array)
386 {
387  const float step_error = fast_length((*dxy2_vec-*dxy1_vec)/GRID_SCALE);
388 
389  *dl = fast_length(*dxy2_vec);
390  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
391  *next_vec = *vec+*dxy1_vec;
392  *dl = fast_length(*dxy1_vec);
393  }
394  *vec = *next_vec;
395  *idx = get_array_idx(*next_vec);
396  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
397  printf("Count-link @ %g,%g: stuck\n",(*vec)[0],(*vec)[1]);
398  atomic_or(&mapping_array[*idx],IS_STUCK);
399  return true;
400  }
401  *dt = select( fmin(DT_MAX,(ADJUSTED_MAX_ERROR*(*dt))/(step_error)), DT_MAX,
402  isequal(step_error,0.0f) );
403  return false;
404 }
405 #endif
406 
407 #if defined(KERNEL_SEGMENT_HILLSLOPES) || defined(KERNEL_SUBSEGMENT_FLANKS)
408 static inline bool segment_runge_kutta_step(float *dt, float *dl,
433  float2 *dxy1_vec, float2 *dxy2_vec,
434  float2 *vec, float2 *next_vec,
435  uint *n_steps, uint *idx)
436 {
437  const float step_error = fast_length((*dxy2_vec-*dxy1_vec)/GRID_SCALE);
438 
439  *dl = fast_length(*dxy2_vec);
440  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
441  *next_vec = *vec+*dxy1_vec;
442  *dl = fast_length(*dxy1_vec);
443  }
444  *vec = *next_vec;
445  *idx = get_array_idx(*next_vec);
446  *n_steps += 1u;
447  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
448  printf("Segment @ %g,%g: stuck\n",(*vec)[0],(*vec)[1]);
449  return true;
450  }
451  *dt = select( fmin(DT_MAX,(ADJUSTED_MAX_ERROR*(*dt))/(step_error)), DT_MAX,
452  isequal(step_error,0.0f) );
453  return false;
454 }
455 #endif
456 
457 #ifdef KERNEL_HILLSLOPE_LENGTHS
458 static inline bool lengths_runge_kutta_step(float *dt, float *dl, float *l_trajectory,
488  float2 *dxy1_vec, float2 *dxy2_vec,
489  float2 *vec, float2 *prev_vec,
490  float2 *next_vec,
491  uint *n_steps, uint *idx)
492 {
493  const float step_error = fast_length((*dxy2_vec-*dxy1_vec)/GRID_SCALE);
494 
495  *dl = fast_length(*dxy2_vec);
496  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
497  *next_vec = *vec+*dxy1_vec;
498  *dl = fast_length(*dxy1_vec);
499  }
500  *vec = *next_vec;
501  *idx = get_array_idx(*next_vec);
502  if (*dl<(INTEGRATION_HALT_THRESHOLD)) {
503  update_trajectory(*dl,l_trajectory,n_steps);
504  return true;
505  }
506  *dt = select( fmin(DT_MAX,(ADJUSTED_MAX_ERROR*(*dt))/(step_error)), DT_MAX,
507  isequal(step_error,0.0f) );
508  update_trajectory(*dl,l_trajectory,n_steps);
509  *prev_vec = *vec;
510  return false;
511 }
512 #endif
#define DT_MAX
Definition: info.h:48
static void euler_step_record(float *dt, float *dl, float *l_trajectory, const float2 uv_vec, float2 *vec, const float2 prev_vec, uint *n_steps, __global char2 *trajectory_vec)
Compute a single Euler integration step of of a streamline.
static void update_trajectory(float dl, float *l_trajectory, uint *n_steps)
Update variables tracking trajectory length and integration step counter.
#define IS_STUCK
Definition: info.h:109
static void update_trajectory_write_sl_data(float dl, float *l_trajectory, float2 vec, float2 prev_vec, uint *n_steps, uint *idx, uint *prev_idx, __global const bool *mask_array, __global uint *slt_array, __global uint *slc_array)
Update variables tracking trajectory length and integration step counter.
static void update_record_private_trajectory(float dl, float *l_trajectory, float2 vec, float2 prev_vec, uint *n_steps, __private char2 *trajectory_vec)
Update variables tracking trajectory length and integration step counter.
static bool lengths_runge_kutta_step(float *dt, float *dl, float *l_trajectory, float2 *dxy1_vec, float2 *dxy2_vec, float2 *vec, float2 *prev_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...
#define INTEGRATION_HALT_THRESHOLD
Definition: info.h:45
static bool runge_kutta_step_record(float *dt, float *dl, float *l_trajectory, float2 *dxy1_vec, float2 *dxy2_vec, float2 *vec, float2 *prev_vec, float2 *next_vec, uint *n_steps, uint *idx, uint *prev_idx, __global char2 *trajectory_vec)
Compute a single step of 2nd-order Runge-Kutta numerical integration of a streamline given precompute...
static float dt_to_nearest_edge(float x, float u)
In Euler streamline integration (which is the last step), this function provides the delta time requi...
Definition: essentials.cl:181
static float2 approximate(float2 raw_position)
Approximate a float vector at the resolution provided by a scaled byte vector.
Definition: essentials.cl:164
static bool runge_kutta_step_write_sl_data(float *dt, float *dl, float *l_trajectory, float2 *dxy1_vec, float2 *dxy2_vec, float2 *vec, float2 *prev_vec, const float2 next_vec, uint *n_steps, uint *idx, uint *prev_idx, __global const bool *mask_array, __global uint *slt_array, __global uint *slc_array)
Compute a single step of 2nd-order Runge-Kutta numerical integration of a streamline given precompute...
#define Y_MAX
Definition: info.h:72
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 X_MAX
Definition: info.h:71
static bool segment_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...
static void euler_step_write_sl_data(float *dt, float *dl, float *l_trajectory, const float2 uv_vec, float2 *vec, const float2 prev_vec, uint *n_steps, uint *idx, uint *prev_idx, __global const bool *mask_array, __global uint *slt_array, __global uint *slc_array)
Compute a single Euler integration step of of a streamline.
#define ADJUSTED_MAX_ERROR
Definition: info.h:43
static void update_record_trajectory(float dl, float *l_trajectory, float2 vec, float2 prev_vec, uint *n_steps, __global char2 *trajectory_vec)
Update variables tracking trajectory length and integration step counter.
#define GRID_SCALE
Definition: info.h:46
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 DOWNUP_SIGN
Definition: info.h:40
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...
static bool connect_runge_kutta_step_record(float *dt, float *dl, float *l_trajectory, float2 *dxy1_vec, float2 *dxy2_vec, float2 *vec, float2 *prev_vec, float2 *next_vec, uint *n_steps, uint *idx, uint *prev_idx, __private char2 *trajectory_vec)
Compute a single step of 2nd-order Runge-Kutta numerical integration of a streamline given precompute...