Source browser

examples/shaders/evt3_flow.cpp

Read the complete file below or download its original bytes. Line numbers provide stable links to a declaration or processing step.

1// Estimate normal optical flow by fitting a local plane to event timestamps.
2#ifdef __cplusplus // Include the runtime only when compiling the CPU adapter.
3#include "shdr.h" // Share the portable resource and math declarations.
4#endif // Finish the C++ declaration include guard.
5layout(local_size_x = 64) in; // Assign one decoded event to each flow-estimation invocation.
6layout(std430, binding = 0) readonly buffer Events { // Access decoded records in their original event order.
7    uvec4 events[]; // Use decoded pixel or trigger records with paired timestamp words.
8}; // Finish this shader storage resource declaration.
9layout(std430, binding = 1) readonly buffer Offsets { // Read each block's exclusive output interval and the final total.
10    uint offsets[]; // Read exclusive event offsets and the final total event count.
11}; // Finish this shader storage resource declaration.
12layout(std430, binding = 2) writeonly buffer Flow { // Write one normal-flow result per decoded event.
13    vec4 flow[]; // Store normal velocity in pixels per second, neighbor count and validity.
14}; // Finish this shader storage resource declaration.
15uniform uint blocks; // Receive the number of summarized input blocks.
16uniform float radius; // Receive the spatial neighborhood radius in pixels.
17uniform uint history; // Receive the maximum preceding records to inspect, capped at 128.
18uniform uint time_window; // Receive the maximum neighbor age in microseconds.
19
20void main() { // Fit a local timestamp plane and publish its observable normal velocity.
21    uint index = gl_GlobalInvocationID.x; // Select the decoded event whose flow is being estimated.
22    if (index >= offsets[blocks]) // Ignore invocations beyond the decoder's final compact event count.
23        return; // Stop dispatch padding at the compact event-count boundary.
24    uvec4 center = events[index]; // Load the current event's pixel metadata and complete timestamp.
25    flow[index] = vec4(0.0); // Default to an invalid zero result until a usable plane is fitted.
26    // Fit spatial neighborhoods around decoded CD pixel events.
27    if ((center.y & 2147483648u) != 0u) // Recognize a trigger record before selecting a pixel neighborhood.
28        return; // Retain the initialized zero flow for trigger records.
29    float sx = 0.0, sy = 0.0, st = 0.0, sxx = 0.0, sxy = 0.0, syy = 0.0, sxt = 0.0, syt = 0.0; // Initialize first and second moments for a centered least-squares plane.
30    uint neighbors = 0u, lookback = min(history, 128u); // Count accepted neighbors and bound the per-event search cost.
31    ivec2 reference_offset = ivec2(0); // Remember the first nonzero integer pixel direction from the center.
32    bool spatial_rank_two = false; // Require two independent directions before fitting a two-dimensional plane.
33    for (uint step = 1u; step <= lookback && step <= index; ++step) { // Inspect only available preceding events within the bounded history.
34        uvec4 previous = events[index - step]; // Load one candidate neighbor from the decoded stream.
35        if ((previous.y & 2147483648u) != 0u || (previous.y & 6144u) != (center.y & 6144u)) // Reject triggers and events with a different polarity or camera.
36            continue; // Exclude this candidate while preserving the accumulated accepted neighbors.
37        // Subtract the timestamp pair exactly, including a low-word borrow.
38        // Only differences fitting uint32 and the requested window become float.
39        if (previous.w > center.w) // Reject timestamps newer than the center in their high word.
40            continue; // Exclude this candidate while preserving the accumulated accepted neighbors.
41        uint high = center.w - previous.w, borrow = center.z < previous.z ? 1u : 0u; // Subtract high words and determine the low-word borrow exactly.
42        if (high < borrow) // Reject a future timestamp whose low-word borrow would underflow.
43            continue; // Exclude this candidate while preserving the accumulated accepted neighbors.
44        high -= borrow; // Apply the borrow to the remaining high-word difference.
45        uint elapsed = center.z - previous.z; // Compute the low-word age with the defined unsigned wrapping subtraction.
46        if (high != 0u || elapsed > time_window) // Reject ages beyond uint32 or the requested temporal neighborhood.
47            continue; // Exclude this candidate while preserving the accumulated accepted neighbors.
48        ivec2 offset = ivec2(int(previous.x) - int(center.x), int(previous.y & 2047u) - int(center.y & 2047u)); // Subtract the decoded eleven-bit pixel coordinates exactly as signed integers.
49        float dx = float(offset.x), dy = float(offset.y); // Convert the small integer offsets for the floating-point moment calculation.
50        if (dx * dx + dy * dy > radius * radius) // Reject candidates outside the circular pixel neighborhood.
51            continue; // Exclude this candidate while preserving the accumulated accepted neighbors.
52        if (!spatial_rank_two) { // Stop testing directions once the neighborhood is known to span a plane.
53            if (reference_offset.x == 0 && reference_offset.y == 0) { // Find the first nonzero pixel direction from the center.
54                reference_offset = offset; // Keep this direction for comparison with later accepted neighbors.
55            } else if (offset.x * reference_offset.y != offset.y * reference_offset.x) { // Test independence with exact integer products bounded by eleven-bit coordinates.
56                spatial_rank_two = true; // Record exact two-dimensional spatial rank from the integer cross product.
57            } // Finish the exact spatial direction test.
58        } // Skip later direction tests once two independent offsets are found.
59        float dt = -float(elapsed); // Express an accepted earlier timestamp relative to the center in microseconds.
60        sx += dx; // Accumulate the horizontal first moment.
61        sy += dy; // Accumulate the vertical first moment.
62        st += dt; // Accumulate the relative timestamp first moment.
63        sxx += dx * dx; // Accumulate the horizontal squared moment.
64        sxy += dx * dy; // Accumulate the horizontal and vertical cross moment.
65        syy += dy * dy; // Accumulate the vertical squared moment.
66        sxt += dx * dt; // Accumulate the horizontal and timestamp cross moment.
67        syt += dy * dt; // Accumulate the vertical and timestamp cross moment.
68        ++neighbors; // Count this accepted spatial and temporal neighbor.
69    } // Finish this decoder or fitting branch.
70    if (neighbors < 6u || !spatial_rank_two) // Require six neighbors and exact two-dimensional spatial rank.
71        return; // Return the initialized zero result for insufficient or collinear neighbors.
72    // Include the center sample and eliminate the intercept with count-scaled moments.
73    // Subtract count-scaled integer moments before the final slope divisions.
74    // The common sample-count factors cancel in the slope solve.
75    float n = float(neighbors + 1u); // Include the center event at zero offset and zero relative time.
76    float xx = n * sxx - sx * sx, xy = n * sxy - sx * sy, yy = n * syy - sy * sy; // Form spatial covariances with the common sample-count scale.
77    float xt = n * sxt - sx * st, yt = n * syt - sy * st; // Form matching count-scaled space-time covariances.
78    float determinant = xx * yy - xy * xy; // Measure invertibility with both covariance factors scaled by the sample count.
79    if (determinant <= 0.0001 * n * n) // Preserve the original determinant threshold under the squared count scaling.
80        return; // Return the initialized zero result for degenerate spatial covariance.
81    float a = (xt * yy - yt * xy) / determinant, b = (yt * xx - xt * xy) / determinant; // Solve the horizontal and vertical timestamp slopes in microseconds per pixel.
82    float gradient = a * a + b * b; // Compute the squared timestamp-gradient magnitude.
83    if (gradient <= 0.000001) // Reject a plane with too little temporal slope to define finite velocity.
84        return; // Return the initialized zero result for a timestamp slope below the velocity threshold.
85    // Timestamps are microseconds, publish pixels/second, neighbor count, validity.
86    flow[index] = vec4(vec2(a, b) * (1000000.0 / gradient), float(neighbors), 1.0); // Convert the observable normal velocity to pixels per second and mark this fit valid.
87} // Finish this decoder or fitting branch.