ogl_beamforming

Ultrasound Beamforming Implemented with OpenGL
git clone anongit@rnpnr.xyz:ogl_beamforming.git
Log | Files | Refs | Feed | Submodules | README | LICENSE

Commit: 1643d83dfb2a08cea4f98c74a0064d8cbd5d08ea
Parent: 64fcfa5c6ac52b400146906fd76c3cf0488251e5
Author: Randy Palamar
Date:   Thu,  5 Nov 2026 17:54:04 -0700

shaders: filter/das: add SamplesPerWave

This should be more accurate than DemodulationFrequency /
SamplingFrequency when we know the value is either 2 or 1 (in IQ
samples).

No difference to compiled DAS on my test dataset but does give
better code for compiled filter shader.

Diffstat:
Mbeamformer.meta | 38+++++++++++++++++++++-----------------
Mbeamformer_core.c | 11++++++-----
Mgenerated/beamformer.c | 20+++++++++++---------
Mshaders/das.glsl | 8++++----
Mshaders/filter.glsl | 4++--
5 files changed, 44 insertions(+), 37 deletions(-)

diff --git a/beamformer.meta b/beamformer.meta @@ -20,11 +20,12 @@ Columns } -@Enumeration SamplingMode +@Table([name samples]) SamplingModeTable { - 2X - 4X + [2X 2] + [4X 4] } +@Expand(SamplingModeTable) @Enumeration(`$(name)`) SamplingMode @Table([name size elements complex glsl]) DataKindTable { @@ -294,6 +295,10 @@ @Expand(DataKindTable) ` str8_comp("$(name)"),` `};` `` + `read_only f32 beamformer_sampling_mode_samples_per_wave[] = {` + @Expand(SamplingModeTable) ` $(samples),` + `};` + `` `read_only u32 beamformer_acquisition_kind_flags[] = {` @Expand(AcquisitionKindTable) ` ($(fixed_transmits)*BeamformerAcquisitionKindFlags_FixedTransmits)|($(sparse)*BeamformerAcquisitionKindFlags_Sparse),` `};` @@ -373,19 +378,18 @@ @Bake { - [FilterCoefficients U32] - [FilterLength U32] - [SamplingFrequency F32] - [DemodulationFrequency F32] - [DecimationRate U32] - [SampleCount U32] - [BatchSampleCount U32] - [InputChannelStride U32] - [InputSampleStride U32] - [InputTransmitStride U32] - [OutputChannelStride U32] - [OutputSampleStride U32] - [OutputTransmitStride U32] + [FilterCoefficients U32] + [FilterLength U32] + [SamplesPerWave F32] + [DecimationRate U32] + [SampleCount U32] + [BatchSampleCount U32] + [InputChannelStride U32] + [InputSampleStride U32] + [InputTransmitStride U32] + [OutputChannelStride U32] + [OutputSampleStride U32] + [OutputTransmitStride U32] } @PushConstants @@ -433,7 +437,7 @@ [ReadiGroupCount U16] [SamplingFrequency F32] - [DemodulationFrequency F32] + [SamplesPerWave F32] [SpeedOfSound F32] [TimeOffset F32] [FNumber F32] diff --git a/beamformer_core.c b/beamformer_core.c @@ -560,6 +560,10 @@ plan_compute_pipeline(BeamformerComputePlan *cp, BeamformerParameterBlock *pb, A u32 input_sample_count = pb->parameters.sample_count; u32 acquisition_count = pb->parameters.acquisition_count; u32 decimation_rate = Max(pb->parameters.decimation_rate, 1); + // NOTE(rnp): only relevant when demodulating + f32 samples_per_wave = (pb->parameters.sampling_frequency / 2.0f) / pb->parameters.demodulation_frequency; + if ((u32)pb->parameters.sampling_mode < (u32)BeamformerSamplingMode_Count) + samples_per_wave = beamformer_sampling_mode_samples_per_wave[pb->parameters.sampling_mode] / 2.0f; cp->raw_channel_byte_stride = pb->parameters.sample_count * pb->parameters.acquisition_count * beamformer_data_kind_byte_size[pb->pipeline.data_kind]; @@ -902,10 +906,7 @@ plan_compute_pipeline(BeamformerComputePlan *cp, BeamformerParameterBlock *pb, A * Q[n] = RF[n + 1] * IQ[n] = I[n] - j*Q[n] */ - if (demod) { - fb->DemodulationFrequency = pb->parameters.demodulation_frequency; - fb->SamplingFrequency = pb->parameters.sampling_frequency / 2; - } + fb->SamplesPerWave = samples_per_wave; sd->layout = (uv3){{subgroup_size, 1, 1}}; sd->dispatch.x = DivRoundUp(node->rf_dimensions.E[BeamformerRFDimension_Samples], sd->layout.x); @@ -918,7 +919,7 @@ plan_compute_pipeline(BeamformerComputePlan *cp, BeamformerParameterBlock *pb, A BeamformerDASBakeParameters *db = &sd->bake.DAS; db->SamplingFrequency = sampling_frequency; - db->DemodulationFrequency = pb->parameters.demodulation_frequency; + db->SamplesPerWave = samples_per_wave / decimation_rate; db->SpeedOfSound = pb->parameters.speed_of_sound; db->TimeOffset = time_offset; db->FNumber = pb->parameters.f_number; diff --git a/generated/beamformer.c b/generated/beamformer.c @@ -195,8 +195,7 @@ typedef struct { typedef struct { u32 FilterCoefficients; u32 FilterLength; - f32 SamplingFrequency; - f32 DemodulationFrequency; + f32 SamplesPerWave; u32 DecimationRate; u32 SampleCount; u32 BatchSampleCount; @@ -229,7 +228,7 @@ typedef struct { u16 InterpolationMode; u16 ReadiGroupCount; f32 SamplingFrequency; - f32 DemodulationFrequency; + f32 SamplesPerWave; f32 SpeedOfSound; f32 TimeOffset; f32 FNumber; @@ -469,7 +468,7 @@ read_only MetaStructMember *meta_struct_members_by_id[] = { {18, 0, 1, 0}, {18, 4, 1, 0}, {8, 8, 1, 0}, - {8, 12, 1, 0}, + {18, 12, 1, 0}, {18, 16, 1, 0}, {18, 20, 1, 0}, {18, 24, 1, 0}, @@ -478,7 +477,6 @@ read_only MetaStructMember *meta_struct_members_by_id[] = { {18, 36, 1, 0}, {18, 40, 1, 0}, {18, 44, 1, 0}, - {18, 48, 1, 0}, }, (MetaStructMember []){ {18, 0, 1, 0}, @@ -544,8 +542,7 @@ read_only str8 *meta_struct_member_names_by_id[] = { (str8 []){ str8_comp("FilterCoefficients"), str8_comp("FilterLength"), - str8_comp("SamplingFrequency"), - str8_comp("DemodulationFrequency"), + str8_comp("SamplesPerWave"), str8_comp("DecimationRate"), str8_comp("SampleCount"), str8_comp("BatchSampleCount"), @@ -577,7 +574,7 @@ read_only str8 *meta_struct_member_names_by_id[] = { str8_comp("InterpolationMode"), str8_comp("ReadiGroupCount"), str8_comp("SamplingFrequency"), - str8_comp("DemodulationFrequency"), + str8_comp("SamplesPerWave"), str8_comp("SpeedOfSound"), str8_comp("TimeOffset"), str8_comp("FNumber"), @@ -605,7 +602,7 @@ read_only str8 *meta_struct_member_names_by_id[] = { read_only MetaStructInfo meta_struct_info_by_id[] = { {str8_comp("DecodeBakeParameters"), 11, 44, 0}, - {str8_comp("FilterBakeParameters"), 13, 52, 0}, + {str8_comp("FilterBakeParameters"), 12, 48, 0}, {str8_comp("DASBakeParameters"), 27, 92, 0}, {str8_comp("CoherencyWeightingBakeParameters"), 3, 12, 0}, {str8_comp("ReshapeBakeParameters"), 9, 36, 0}, @@ -940,6 +937,11 @@ read_only str8 beamformer_data_kind_str8[] = { str8_comp("Float16Complex"), }; +read_only f32 beamformer_sampling_mode_samples_per_wave[] = { + 2, + 4, +}; + read_only u32 beamformer_acquisition_kind_flags[] = { (1*BeamformerAcquisitionKindFlags_FixedTransmits)|(0*BeamformerAcquisitionKindFlags_Sparse), (0*BeamformerAcquisitionKindFlags_FixedTransmits)|(1*BeamformerAcquisitionKindFlags_Sparse), diff --git a/shaders/das.glsl b/shaders/das.glsl @@ -57,7 +57,7 @@ const f32 C_SPLINE = 0.5f; #if InputDataKind == DataKind_Float32Complex || InputDataKind == DataKind_Float16Complex vec2 rotate_iq(const vec2 iq, const float time) { - float arg = radians(360) * DemodulationFrequency * time; + float arg = radians(360) * time / SamplesPerWave; mat2 phasor = mat2( cos(arg), sin(arg), -sin(arg), cos(arg)); vec2 result = phasor * iq; @@ -153,7 +153,7 @@ SAMPLE_TYPE sample_rf(const u64 rf_pointer, const f32 index) switch (u32(InterpolationMode)) { case InterpolationMode_Nearest:{ if (index >= 0.f && index < (f32(SampleCount) - 0.5f)) - result = rotate_iq(Input(rf_pointer + InputDataKindByteSize * u32(round(index))).x[0], index / SamplingFrequency); + result = rotate_iq(Input(rf_pointer + InputDataKindByteSize * u32(round(index))).x[0], index); }break; case InterpolationMode_Linear:{ if (index >= 0.f && index < f32(SampleCount - 1)) { @@ -171,12 +171,12 @@ SAMPLE_TYPE sample_rf(const u64 rf_pointer, const f32 index) f32 t = fract(index); result = (1 - t) * rf[0] + t * rf[1]; - result = rotate_iq(result, index / SamplingFrequency); + result = rotate_iq(result, index); } }break; case InterpolationMode_Cubic:{ if (index >= 1.f && index < f32(SampleCount - 2)) - result = rotate_iq(cubic(rf_pointer + InputDataKindByteSize * u32(index), fract(index)), index / SamplingFrequency); + result = rotate_iq(cubic(rf_pointer + InputDataKindByteSize * u32(index), fract(index)), index); }break; } return result; diff --git a/shaders/filter.glsl b/shaders/filter.glsl @@ -50,9 +50,9 @@ f32vec2 complex_mul(f32vec2 a, f32vec2 b) } #if Demodulate -SAMPLE_TYPE rotate_iq(SAMPLE_TYPE iq, uint index) +SAMPLE_TYPE rotate_iq(SAMPLE_TYPE iq, u32 index) { - float arg = radians(360) * DemodulationFrequency * index / SamplingFrequency; + f32 arg = radians(360) * index / SamplesPerWave; SAMPLE_TYPE result = SAMPLE_TYPE(complex_mul(iq, f32vec2(cos(arg), -sin(arg)))); return result; }