Commit: 31d7dcd16a2e31e040cafdba26966d93385449bf
Parent: d72274eb5a685c634e5129fbdfdea6f618c22f08
Author: Randy Palamar
Date: Fri, 28 Aug 2026 06:23:06 -0700
core/das: allow das to handle F16 inputs
This is useful for VLS/TPW since we can just let the demodulation
data kind flow through without issue. There may be precision loss
by storing the filter result in F16 but in practice I haven't seen
an issue (and we have already been doing it for FORCES/HERCULES).
I could see it mattering for a PSF simulation study but in that
case you can just pass F32 data to the beamformer and everything
will happen in F32.
Diffstat:
| M | beamformer_core.c | | | 98 | +++++++++++++++++++++++++++++++++++++++++++++---------------------------------- |
| M | shaders/das.glsl | | | 73 | +++++++++++++++++++++++++++++++++++++++++++++---------------------------- |
2 files changed, 101 insertions(+), 70 deletions(-)
diff --git a/beamformer_core.c b/beamformer_core.c
@@ -1,5 +1,9 @@
/* See LICENSE for license details. */
/* TODO(rnp):
+ * [ ]: lib: perform conversion from S16->F16 while copying
+ * [ ]: bug: out of bounds reads and extra computation when ChunkChannelCount is not a multiple
+ * of (receive) ChannelCount.
+ * - probably won't show up with any of our arrays anytime soon but should be fixed
* [ ]: backtrace dumping on SIGSEGV
* [ ]: cooperative shared memory loading in decode shader
* [ ]: refactor: save filter parameters with rest of parameters, whole slot thing is dumb
@@ -39,18 +43,6 @@
* 1. Maximize value of wire target
* 2. Maximize value of single cyst contrast
*
- * [ ]: Need somewhere to store display image temporaries
- * - Coherency Weighting needs somewhere to put its incoherent sum
- * - Image Averaging needs to put its output somewhere
- * - Recursive imaging needs to be able to sum/subtract to update on each recursion step
- * - Power doppler needs somewhere to store current power map
- * - We don't need a backlog of these but it may be useful to have a front and back buffer
- * so that the UI always sees consistent state
- * - Really what we want is a separate temp (GPU) arena that gets recreated on
- * pipeline creation (in plan_compute_pipeline()).
- * - NOTE: image offsets/etc needed below don't have a predefined sizes so they would
- * also be a candidates for storage in GPU temp arena.
- *
* [ ]: Recursive Imaging Handling
* [ ]: add array of image offsets to Compute Array Parameters
* [ ]: add array of weighting coeffiecents to Compute Array Parameters
@@ -75,33 +67,6 @@
#include "beamformer_internal.h"
-typedef struct BeamformerComputeGraphNode BeamformerComputeGraphNode;
-struct BeamformerComputeGraphNode {
- // NOTE(rnp): will be BeamformerShaderKind_Count for root node
- BeamformerShaderKind kind;
-
- // NOTE(rnp): when any of input or output stride is assigned it is assumed that
- // the shader requires a fixed layout for input, output, or both. When two adjacent
- // nodes require incompatible layouts the second pass over the graph will insert
- // Reshape shaders in between.
- BeamformerDataKind input_data_kind;
- iv3 input_stride;
-
- BeamformerDataKind output_data_kind;
- iv3 output_stride;
-
- i32 user_pipeline_index;
-
- BeamformerComputeGraphNode *prev;
- BeamformerComputeGraphNode *next;
-};
-
-typedef struct {
- BeamformerComputeGraphNode *first;
- BeamformerComputeGraphNode *last;
- u64 count;
-} BeamformerComputeGraph;
-
#define GPU_RESOURCE_HASH_TABLE_COUNT 256
typedef struct U64ReferenceNode U64ReferenceNode;
struct U64ReferenceNode {u64 *v; U64ReferenceNode *next;};
@@ -502,6 +467,35 @@ dispatch_for_output(uv3 layout, iv3 points)
return result;
}
+typedef struct BeamformerComputeGraphNode BeamformerComputeGraphNode;
+struct BeamformerComputeGraphNode {
+ // NOTE(rnp): will be BeamformerShaderKind_Count for root node
+ BeamformerShaderKind kind;
+
+ // NOTE(rnp): when any of input or output stride is assigned it is assumed that
+ // the shader requires a fixed layout for input, output, or both. When two adjacent
+ // nodes require incompatible layouts the second pass over the graph will insert
+ // Reshape shaders in between.
+ BeamformerDataKind input_data_kind;
+ iv3 input_stride;
+
+ BeamformerDataKind output_data_kind;
+ iv3 output_stride;
+
+ b32 requires_fp_input;
+
+ i32 user_pipeline_index;
+
+ BeamformerComputeGraphNode *prev;
+ BeamformerComputeGraphNode *next;
+};
+
+typedef struct {
+ BeamformerComputeGraphNode *first;
+ BeamformerComputeGraphNode *last;
+ u64 count;
+} BeamformerComputeGraph;
+
function b32
compute_plan_push_shader(BeamformerComputePlan *p, BeamformerComputeGraphNode *node, BeamformerShaderParameters *sp)
{
@@ -609,6 +603,15 @@ plan_compute_pipeline(BeamformerComputePlan *cp, BeamformerParameterBlock *pb, A
[BeamformerDataKind_Float32Complex] = BeamformerDataKind_Float32,
};
+ read_only local_persist BeamformerDataKind data_kind_to_fp_kind[] = {
+ [BeamformerDataKind_Int16] = BeamformerDataKind_Float16,
+ [BeamformerDataKind_Float16] = BeamformerDataKind_Float16,
+ [BeamformerDataKind_Float32] = BeamformerDataKind_Float32,
+ [BeamformerDataKind_Int16Complex] = BeamformerDataKind_Float16Complex,
+ [BeamformerDataKind_Float16Complex] = BeamformerDataKind_Float16Complex,
+ [BeamformerDataKind_Float32Complex] = BeamformerDataKind_Float32Complex,
+ };
+
//////////////////////////////////////
// NOTE(rnp): First Pass: build initial graph and insert hard layout constraints
BeamformerComputeGraph graph = {0};
@@ -652,9 +655,9 @@ plan_compute_pipeline(BeamformerComputePlan *cp, BeamformerParameterBlock *pb, A
(acquisition_count % 16 == 0) &&
(chunk_channel_count % 16 == 0);
+ node->requires_fp_input = 1;
+
// NOTE(rnp): fixed input layout required for reasonable performance
- if (low_precision && beamformer_data_kind_complex[input_data_kind])
- node->input_data_kind = BeamformerDataKind_Float16Complex;
node->input_stride.x = chunk_channel_count * acquisition_count;
node->input_stride.y = acquisition_count;
node->input_stride.z = 1;
@@ -667,7 +670,12 @@ plan_compute_pipeline(BeamformerComputePlan *cp, BeamformerParameterBlock *pb, A
}break;
case BeamformerShaderKind_DAS:{
- node->input_data_kind = das_data_kind;
+ // NOTE(rnp): if we aren't decoding we can let the previous stage's
+ // data kind pass through as long as it is floating point
+ node->requires_fp_input = 1;
+ if (pb->parameters.decode_mode != BeamformerDecodeMode_None)
+ node->input_data_kind = das_data_kind;
+
node->input_stride.x = 1; // Sample Stride
node->input_stride.y = input_sample_count * acquisition_count; // Channel Stride
node->input_stride.z = input_sample_count; // Receive Event Stride
@@ -721,6 +729,12 @@ plan_compute_pipeline(BeamformerComputePlan *cp, BeamformerParameterBlock *pb, A
if (prev_output_dont_care && input_dont_care)
node->input_data_kind = node->prev->output_data_kind = node->prev->input_data_kind;
+ if (node->requires_fp_input) {
+ node->input_data_kind = data_kind_to_fp_kind[node->input_data_kind];
+ if (prev_output_dont_care)
+ node->prev->output_data_kind = node->input_data_kind;
+ }
+
needs_reshape |= node->input_data_kind != node->prev->output_data_kind;
}
diff --git a/shaders/das.glsl b/shaders/das.glsl
@@ -1,12 +1,12 @@
/* See LICENSE for license details. */
-#if InputDataKind == DataKind_Float32
+#if InputDataKind == DataKind_Float32 || InputDataKind == DataKind_Float16
#if CoherencyWeighting
#define RESULT_TYPE vec2
#define RESULT_COHERENT_CAST(a) (a).x
#define RESULT_INCOHERENT_CAST(a) (a).y
#endif
#define SAMPLE_TYPE f32
-#elif InputDataKind == DataKind_Float32Complex
+#elif InputDataKind == DataKind_Float32Complex || InputDataKind == DataKind_Float16Complex
#if CoherencyWeighting
#define RESULT_TYPE vec3
#define RESULT_COHERENT_CAST(a) (a).xy
@@ -38,19 +38,22 @@ layout(std430, buffer_reference) buffer IncoherentOutput {
f32 x[];
};
-layout(std430, buffer_reference) buffer F16 { f16 x[]; };
-layout(std430, buffer_reference) buffer F32 { f32 x[]; };
-layout(std430, buffer_reference) buffer S16 { s16 x[]; };
-layout(std430, buffer_reference) buffer U8 { u8 x[]; };
-layout(std430, buffer_reference) buffer V2 { vec2 x[]; };
-layout(std430, buffer_reference) buffer V4 { vec4 x[]; };
+layout(std430, buffer_reference) buffer F16 { f16 x[]; };
+layout(std430, buffer_reference) buffer F32 { f32 x[]; };
+layout(std430, buffer_reference) buffer S16 { s16 x[]; };
+layout(std430, buffer_reference) buffer U8 { u8 x[]; };
+layout(std430, buffer_reference) buffer U32V4 { u32vec4 x[]; };
+layout(std430, buffer_reference) buffer F32V2 { f32vec2 x[]; };
+layout(std430, buffer_reference) buffer F32V4 { f32vec4 x[]; };
+layout(std430, buffer_reference) buffer F16V2 { f16vec2 x[]; };
+layout(std430, buffer_reference) buffer F16V4 { f16vec4 x[]; };
#define RX_ORIENTATION(tx_rx) bitfieldExtract((tx_rx), 0, 4)
#define TX_ORIENTATION(tx_rx) bitfieldExtract((tx_rx), 4, 4)
#define C_SPLINE 0.5
-#if InputDataKind == DataKind_Float32Complex
+#if InputDataKind == DataKind_Float32Complex || InputDataKind == DataKind_Float16Complex
vec2 rotate_iq(const vec2 iq, const float time)
{
float arg = radians(360) * DemodulationFrequency * time;
@@ -74,29 +77,37 @@ SAMPLE_TYPE cubic(const u64 rf_pointer, const f32 t)
);
#if InputDataKind == DataKind_Float32
- vec4 samples = V4(rf_pointer).x[0];
+ f32vec4 samples = F32V4(rf_pointer).x[0];
+ #elif InputDataKind == DataKind_Float16
+ f16vec4 samples = F16V4(rf_pointer).x[0];
+ #elif InputDataKind == DataKind_Float16Complex
+ f32vec2 samples[4];
+ uvec4 load = U32V4(rf_pointer).x[0];
+ samples[0] = unpackHalf2x16(load[0]);
+ samples[1] = unpackHalf2x16(load[1]);
+ samples[2] = unpackHalf2x16(load[2]);
+ samples[3] = unpackHalf2x16(load[3]);
#else
- SAMPLE_TYPE samples[4];
- vec4 load1 = V4(rf_pointer).x[0];
- vec4 load2 = V4(rf_pointer).x[1];
- samples[0] = load1.xy;
- samples[1] = load1.zw;
- samples[2] = load2.xy;
- samples[3] = load2.zw;
+ f32vec2 samples[4];
+ vec4 load1 = F32V4(rf_pointer).x[0];
+ vec4 load2 = F32V4(rf_pointer).x[1];
+ samples[0] = load1.xy;
+ samples[1] = load1.zw;
+ samples[2] = load2.xy;
+ samples[3] = load2.zw;
#endif
- vec4 S = vec4(t * t * t, t * t, t, 1);
+ vec4 Sh = vec4(t * t * t, t * t, t, 1) * h;
SAMPLE_TYPE P1 = samples[1];
SAMPLE_TYPE P2 = samples[2];
SAMPLE_TYPE T1 = C_SPLINE * (P2 - samples[0]);
SAMPLE_TYPE T2 = C_SPLINE * (samples[3] - P1);
- #if InputDataKind == DataKind_Float32
- vec4 C = vec4(P1.x, P2.x, T1.x, T2.x);
- SAMPLE_TYPE result = dot(S, h * C);
- #elif InputDataKind == DataKind_Float32Complex
- mat2x4 C = mat2x4(vec4(P1.x, P2.x, T1.x, T2.x), vec4(P1.y, P2.y, T1.y, T2.y));
- SAMPLE_TYPE result = S * h * C;
+ #if InputDataKind == DataKind_Float32 || InputDataKind == DataKind_Float16
+ SAMPLE_TYPE result = dot(Sh, vec4(P1, P2, T1, T2));
+ #else
+ mat2x4 C = mat2x4(vec4(P1.x, P2.x, T1.x, T2.x), vec4(P1.y, P2.y, T1.y, T2.y));
+ SAMPLE_TYPE result = Sh * C;
#endif
return result;
}
@@ -113,11 +124,17 @@ SAMPLE_TYPE sample_rf(const u64 rf_pointer, const f32 index)
case InterpolationMode_Linear:{
if (index >= 0.f && index < f32(SampleCount - 1)) {
#if InputDataKind == DataKind_Float32
- vec2 rf = V2(rf_pointer + InputDataKindByteSize * u32(index)).x[0];
+ f32vec2 rf = F32V2(rf_pointer + InputDataKindByteSize * u32(index)).x[0];
+ #elif InputDataKind == DataKind_Float16
+ f16vec2 rf = F16V2(rf_pointer + InputDataKindByteSize * u32(index)).x[0];
+ #elif InputDataKind == DataKind_Float16Complex
+ f16vec4 load = F16V4(rf_pointer + InputDataKindByteSize * u32(index)).x[0];
+ f16vec2 rf[2] = {load.xy, load.zw};
#else
- vec4 load = V4(rf_pointer + InputDataKindByteSize * u32(index)).x[0];
- vec2 rf[2] = {load.xy, load.zw};
+ f32vec4 load = F32V4(rf_pointer + InputDataKindByteSize * u32(index)).x[0];
+ f32vec2 rf[2] = {load.xy, load.zw};
#endif
+
f32 t = fract(index);
result = (1 - t) * rf[0] + t * rf[1];
result = rotate_iq(result, index / SamplingFrequency);
@@ -186,7 +203,7 @@ u8 tx_rx_orientation_for_acquisition(const s32 acquisition)
f32vec2 focal_vector_for_acquisition(const s32 acquisition)
{
- f32vec2 result = SingleFocus ? f32vec2(TransmitAngle, FocusDepth) : V2(FocalVectors).x[acquisition];
+ f32vec2 result = SingleFocus ? f32vec2(TransmitAngle, FocusDepth) : F32V2(FocalVectors).x[acquisition];
return result;
}