ogl_beamforming

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

Commit: 78df1e275ea8d36581169c82dc96df3471c6a78a
Parent: dd12be90aa918af0e4f9f7de452a9f1b946928b7
Author: Randy Palamar
Date:   Mon, 27 Jan 2025 20:20:34 -0700

adapt shaders and parameters to defined ornot api

Co-Authored-By: Darath <48868688+darathdev@users.noreply.github.com>

Diffstat:
Mbeamformer.c | 3+--
Mbeamformer_parameters.h | 31++++++++++++++++++++-----------
Mshaders/das.glsl | 75++++++++++++++++++++++++++++++++++++++++-----------------------------------
Mshaders/hadamard.glsl | 12+++++++-----
4 files changed, 68 insertions(+), 53 deletions(-)

diff --git a/beamformer.c b/beamformer.c @@ -396,8 +396,7 @@ do_compute_shader(BeamformerCtx *ctx, Arena arena, BeamformFrame *frame, u32 raw csctx->last_output_ssbo_index = !csctx->last_output_ssbo_index; break; case CS_CUDA_DECODE: - ctx->cuda_lib.cuda_decode(raw_data_index * rf_raw_size, output_ssbo_idx, - ctx->params->raw.channel_offset); + ctx->cuda_lib.cuda_decode(raw_data_index * rf_raw_size, output_ssbo_idx, 0); csctx->raw_data_fences[raw_data_index] = glFenceSync(GL_SYNC_GPU_COMMANDS_COMPLETE, 0); csctx->last_output_ssbo_index = !csctx->last_output_ssbo_index; break; diff --git a/beamformer_parameters.h b/beamformer_parameters.h @@ -17,9 +17,14 @@ enum compute_shaders { CS_LAST }; -#define DAS_ID_UFORCES 0 -#define DAS_ID_HERCULES 1 -#define DAS_ID_RCA 2 +#define DECODE_MODE_NONE 0 +#define DECODE_MODE_HADAMARD 1 + +#define DAS_ID_FORCES 0 +#define DAS_ID_UFORCES 1 +#define DAS_ID_HERCULES 2 +#define DAS_ID_RCA_VLS 3 +#define DAS_ID_RCA_TPW 4 #define MAX_BEAMFORMED_SAVED_FRAMES 16 /* NOTE: This struct follows the OpenGL std140 layout. DO NOT modify unless you have @@ -36,9 +41,8 @@ typedef struct { v4 output_max_coordinate; /* [m] Front-Bottom-Right corner of output region (w ignored)*/ f32 xdc_element_pitch[2]; /* [m] Transducer Element Pitch {row, col} */ uv2 rf_raw_dim; /* Raw Data Dimensions */ - u32 transmit_mode; /* Method/Orientation of Transmit */ + i32 transmit_mode; /* Method/Orientation of Transmit */ u32 decode; /* Decode or just reshape data */ - u32 channel_offset; /* Offset into channel_mapping: 0 or 128 (rows or columns) */ f32 speed_of_sound; /* [m/s] */ f32 sampling_frequency; /* [Hz] */ f32 center_frequency; /* [Hz] */ @@ -47,7 +51,7 @@ typedef struct { i32 beamform_plane; /* Plane to Beamform in 2D HERCULES */ f32 f_number; /* F# (set to 0 to disable) */ u32 das_shader_id; - f32 _pad[1]; + f32 _pad[2]; } BeamformerParameters; /* NOTE: garbage to get the prepocessor to properly stringize the value of a macro */ @@ -69,9 +73,8 @@ layout(std140, binding = 0) uniform parameters {\n\ vec4 output_max_coord; /* [m] Bottom right corner of output region */\n\ vec2 xdc_element_pitch; /* [m] Transducer Element Pitch {row, col} */\n\ uvec2 rf_raw_dim; /* Raw Data Dimensions */\n\ - uint transmit_mode; /* Method/Orientation of Transmit */\n\ + int transmit_mode; /* Method/Orientation of Transmit */\n\ uint decode; /* Decode or just reshape data */\n\ - uint channel_offset; /* Offset into channel_mapping: 0 or 128 (rows or columns) */\n\ float speed_of_sound; /* [m/s] */\n\ float sampling_frequency; /* [Hz] */\n\ float center_frequency; /* [Hz] */\n\ @@ -82,7 +85,13 @@ layout(std140, binding = 0) uniform parameters {\n\ uint das_shader_id;\n\ };\n\ \n\ -#define DAS_ID_UFORCES " str(DAS_ID_UFORCES) "\n\ -#define DAS_ID_HERCULES " str(DAS_ID_HERCULES) "\n\ -#define DAS_ID_RCA " str(DAS_ID_RCA) "\n\n\ +#define DECODE_MODE_NONE " str(DECODE_MODE_NONE) "\n\ +#define DECODE_MODE_HADAMARD " str(DECODE_MODE_HADAMARD) "\n\ +\n\ +#define DAS_ID_FORCES " str(DAS_ID_FORCES) "\n\ +#define DAS_ID_UFORCES " str(DAS_ID_UFORCES) "\n\ +#define DAS_ID_HERCULES " str(DAS_ID_HERCULES) "\n\ +#define DAS_ID_RCA_VLS " str(DAS_ID_RCA_VLS) "\n\ +#define DAS_ID_RCA_TPW " str(DAS_ID_RCA_TPW) "\n\ +\n\ #line 0\n" diff --git a/shaders/das.glsl b/shaders/das.glsl @@ -16,6 +16,9 @@ layout(location = 4) uniform float u_cycle_t; #define TX_ROWS 0 #define TX_COLS 1 +#define TX_MODE_TX_COLS(a) (((a) & 2) != 0) +#define TX_MODE_RX_COLS(a) (((a) & 1) != 0) + #if 1 /* NOTE: interpolation is unnecessary if the data has been demodulated and not decimated */ vec2 cubic(uint ridx, float t) @@ -55,12 +58,14 @@ vec3 calc_image_point(vec3 voxel) vec3 image_point = output_min_coord.xyz + voxel * output_size.xyz / out_data_dim; switch (das_shader_id) { + case DAS_ID_FORCES: case DAS_ID_UFORCES: /* TODO: fix the math so that the image plane can be aritrary */ image_point.y = 0; break; case DAS_ID_HERCULES: - case DAS_ID_RCA: + case DAS_ID_RCA_TPW: + case DAS_ID_RCA_VLS: if (u_volume_export_pass == 0) image_point.y = off_axis_pos; break; @@ -76,15 +81,15 @@ vec2 apodize(vec2 value, float apodization_arg, float distance) return value * a * a; } -vec3 orientation_projection(bool tx_rows) +vec3 orientation_projection(bool rows) { - return vec3(float(!tx_rows), float(tx_rows), 1); + return vec3(float(!rows), float(rows), 1); } -vec3 world_space_to_rca_space(vec3 image_point, int transmit_orientation) +vec3 world_space_to_rca_space(vec3 image_point, bool rx_rows) { vec3 result = (xdc_transform * vec4(image_point, 1)).xyz; - result *= orientation_projection(transmit_orientation != TX_ROWS); + result *= orientation_projection(rx_rows); return result; } @@ -94,25 +99,29 @@ float sample_index(float distance) return time * sampling_frequency; } -float planewave_transmit_distance(vec3 point, float transmit_angle, int transmit_orientation) +float planewave_transmit_distance(vec3 point, float transmit_angle, bool tx_rows) { - return dot(point * orientation_projection(transmit_orientation == TX_ROWS), + return dot(point * orientation_projection(tx_rows), vec3(sin(transmit_angle), sin(transmit_angle), cos(transmit_angle))); } -float cylindricalwave_transmit_distance(vec3 point, float focal_depth, float transmit_angle, int transmit_orientation) +float cylindricalwave_transmit_distance(vec3 point, float focal_depth, float transmit_angle, bool tx_rows) { vec3 f = focal_depth * vec3(sin(transmit_angle), sin(transmit_angle), cos(transmit_angle)); - return length((point - f) * orientation_projection(transmit_orientation == TX_ROWS)); + return length((point - f) * orientation_projection(tx_rows)); } vec2 RCA(vec3 image_point, vec3 delta, float apodization_arg) { - /* TODO: pass this in (there is a problem in that it depends on the orientation - * of the array relative to the target/subject). */ uint ridx = 0; int direction = beamform_plane * (u_volume_export_pass ^ 1); - if (direction == TX_ROWS) image_point = image_point.yxz; + if (direction != TX_ROWS) image_point = image_point.yxz; + + bool tx_col = TX_MODE_TX_COLS(transmit_mode); + bool rx_col = TX_MODE_RX_COLS(transmit_mode); + + vec3 recieve_point = world_space_to_rca_space(image_point, !rx_col); + delta *= orientation_projection(!rx_col); vec2 sum = vec2(0); /* NOTE: For Each Acquistion in Raw Data */ @@ -121,20 +130,17 @@ vec2 RCA(vec3 image_point, vec3 delta, float apodization_arg) uint base_idx = i / 4; uint sub_idx = i % 4; - //int transmit_orientation = int(transmit_orientations[base_idx][sub_idx]); - int transmit_orientation = TX_ROWS; float focal_depth = focal_depths[base_idx][sub_idx]; float transmit_angle = transmit_angles[base_idx][sub_idx]; - vec3 recieve_point = world_space_to_rca_space(image_point, transmit_orientation); float transmit_distance; if (isinf(focal_depth)) { transmit_distance = planewave_transmit_distance(image_point, transmit_angle, - transmit_orientation); + !tx_col); } else { transmit_distance = cylindricalwave_transmit_distance(image_point, focal_depth, transmit_angle, - transmit_orientation); + !tx_col); } vec3 receive_distance = recieve_point; @@ -143,8 +149,8 @@ vec2 RCA(vec3 image_point, vec3 delta, float apodization_arg) for (uint j = 0; j < dec_data_dim.y; j++) { float sidx = sample_index(transmit_distance + length(receive_distance)); vec2 valid = vec2(sidx >= 0) * vec2(sidx < dec_data_dim.x); - sum += apodize(cubic(ridx, sidx), apodization_arg, receive_distance[0]) * valid; - receive_distance[transmit_orientation] -= delta[transmit_orientation]; + sum += apodize(cubic(ridx, sidx), apodization_arg, length(receive_distance.xy)) * valid; + receive_distance -= delta; ridx += dec_data_dim.x; } } @@ -157,24 +163,21 @@ vec2 HERCULES(vec3 image_point, vec3 delta, float apodization_arg) int direction = beamform_plane * (u_volume_export_pass ^ 1); if (direction != TX_ROWS) image_point = image_point.yxz; - //int transmit_orientation = int(transmit_orientations[0][0]); - int transmit_orientation = TX_ROWS; - float focal_depth = focal_depths[0][0]; - float transmit_angle = transmit_angles[0][0]; + bool tx_col = TX_MODE_TX_COLS(transmit_mode); + bool rx_col = TX_MODE_RX_COLS(transmit_mode); + float focal_depth = focal_depths[0][0]; + float transmit_angle = transmit_angles[0][0]; vec3 recieve_point = (xdc_transform * vec4(image_point, 1)).xyz; - if (transmit_orientation == TX_ROWS) - delta = delta.yxz; - float transmit_distance; if (isinf(focal_depth)) { transmit_distance = planewave_transmit_distance(image_point, transmit_angle, - transmit_orientation); + !tx_col); } else { transmit_distance = cylindricalwave_transmit_distance(image_point, focal_depth, transmit_angle, - transmit_orientation); + !tx_col); } vec2 sum = vec2(0); @@ -183,8 +186,8 @@ vec2 HERCULES(vec3 image_point, vec3 delta, float apodization_arg) /* NOTE: For Each Virtual Source */ for (uint j = 0; j < dec_data_dim.y; j++) { vec3 element_position; - if (transmit_orientation == TX_ROWS) element_position = vec3(j, i, 0) * delta; - else element_position = vec3(i, j, 0) * delta; + if (rx_col) element_position = vec3(j, i, 0) * delta; + else element_position = vec3(i, j, 0) * delta; vec3 receive_distance = recieve_point - element_position; float sidx = sample_index(transmit_distance + length(receive_distance)); vec2 valid = vec2(sidx >= 0) * vec2(sidx < dec_data_dim.x); @@ -192,8 +195,9 @@ vec2 HERCULES(vec3 image_point, vec3 delta, float apodization_arg) /* NOTE: tribal knowledge */ if (i == 0) valid *= inversesqrt(dec_data_dim.z); - sum += apodize(cubic(ridx, sidx), apodization_arg, length(receive_distance.xy)) * valid; - ridx += dec_data_dim.x; + sum += apodize(cubic(ridx, sidx), apodization_arg, + length(receive_distance.xy)) * valid; + ridx += dec_data_dim.x; } } return sum; @@ -202,13 +206,12 @@ vec2 HERCULES(vec3 image_point, vec3 delta, float apodization_arg) vec2 uFORCES(vec3 image_point, vec3 delta, float apodization_arg) { /* NOTE: skip first acquisition in uforces since its garbage */ - uint uforces = uint(dec_data_dim.y != dec_data_dim.z); + uint uforces = uint(das_shader_id == DAS_ID_UFORCES); uint ridx = dec_data_dim.y * dec_data_dim.x * uforces; image_point = (xdc_transform * vec4(image_point, 1)).xyz; int transmit_orientation = TX_ROWS; - //int transmit_orientation = int(transmit_orientations[0][0]); if (transmit_orientation == TX_ROWS) delta = delta.yxz; @@ -257,13 +260,15 @@ void main() /* NOTE: skip over channels corresponding to other arrays */ vec2 sum; switch (das_shader_id) { + case DAS_ID_FORCES: case DAS_ID_UFORCES: sum = uFORCES(image_point, vec3(xdc_element_pitch, 0), apod_arg); break; case DAS_ID_HERCULES: sum = HERCULES(image_point, vec3(xdc_element_pitch, 0), apod_arg); break; - case DAS_ID_RCA: + case DAS_ID_RCA_TPW: + case DAS_ID_RCA_VLS: sum = RCA(image_point, vec3(xdc_element_pitch, 0), apod_arg); break; } diff --git a/shaders/hadamard.glsl b/shaders/hadamard.glsl @@ -25,7 +25,6 @@ void main() /* NOTE: offsets for storing the results in the output data */ uint out_off = dec_data_dim.x * dec_data_dim.y * acq + dec_data_dim.x * channel + time_sample; - channel += channel_offset; /* NOTE: channel mapping is stored as u16s so we must do this to extract the final value */ uint ch_array_idx = (channel / 8); uint ch_vec_idx = (channel % 8) / 2; @@ -49,7 +48,12 @@ void main() /* NOTE: Compute N-D dot product */ int sum = 0; - if (decode == 1) { + switch (decode) { + case DECODE_MODE_NONE: { + ridx += ridx_delta * acq; + sum = (rf_data[ridx] << lfs) >> 16; + } break; + case DECODE_MODE_HADAMARD: { /* NOTE: offset to get the correct column in hadamard matrix */ uint hoff = dec_data_dim.z * acq; @@ -58,9 +62,7 @@ void main() sum += imageLoad(hadamard, ivec2(acq, i)).x * data; ridx += ridx_delta; } - } else { - ridx += ridx_delta * acq; - sum = (rf_data[ridx] << lfs) >> 16; + } break; } out_data[out_off] = vec2(float(sum), 0); }