Commit: e8e91267bc6f03da7de75be806398303b151eff5
Parent: 170c7446582f4cacf1761e0ec8c5e2924884bd81
Author: Randy Palamar
Date: Thu, 23 Jan 2025 10:16:54 -0700
core/das: use xdc_transform matrices sent in UBO
Diffstat:
4 files changed, 32 insertions(+), 85 deletions(-)
diff --git a/beamformer.c b/beamformer.c
@@ -290,39 +290,6 @@ beamform_work_queue_push(BeamformerCtx *ctx, Arena *a, enum beamform_work work_t
return result;
}
-static m4
-v3_to_xdc_space(v3 direction, v3 origin, v3 corner1, v3 corner2)
-{
- v3 edge1 = sub_v3(corner1, origin);
- v3 edge2 = sub_v3(corner2, origin);
- v3 xdc_normal = cross(edge1, edge2);
- if (xdc_normal.z < 0)
- xdc_normal = cross(edge2, edge1);
- ASSERT(xdc_normal.z >= 0);
-
- v3 e1 = normalize_v3(sub_v3(direction, xdc_normal));
- v3 e2 = {.y = 1};
- v3 e3 = normalize_v3(cross(e2, e1));
- v4 e4 = {.x = -origin.x, .y = -origin.y, .z = -origin.z, .w = 1};
-
- m4 result = {
- .c[0] = (v4){.x = e3.x, .y = e2.x, .z = e1.x, .w = 0},
- .c[1] = (v4){.x = e3.y, .y = e2.y, .z = e1.y, .w = 0},
- .c[2] = (v4){.x = e3.z, .y = e2.z, .z = e1.z, .w = 0},
- .c[3] = e4,
- };
-
- return result;
-}
-
-static v4
-f32_4_to_v4(f32 *in)
-{
- v4 result;
- store_f32x4(load_f32x4(in), result.E);
- return result;
-}
-
static void
export_frame(BeamformerCtx *ctx, iptr handle, BeamformFrame *frame)
{
@@ -365,13 +332,8 @@ do_beamform_shader(ComputeShaderCtx *cs, BeamformerParameters *bp, BeamformFrame
for (u32 i = 0; i < frame->dim.w; i++) {
u32 texture = frame->textures[i];
- m4 xdc_transform = v3_to_xdc_space((v3){.z = 1},
- f32_4_to_v4(bp->xdc_origin + (4 * i)).xyz,
- f32_4_to_v4(bp->xdc_corner1 + (4 * i)).xyz,
- f32_4_to_v4(bp->xdc_corner2 + (4 * i)).xyz);
glBindImageTexture(0, texture, 0, GL_TRUE, 0, GL_WRITE_ONLY, GL_RG32F);
glUniform1i(cs->xdc_index_id, i);
- glUniformMatrix4fv(cs->xdc_transform_id, 1, GL_FALSE, xdc_transform.E);
glDispatchCompute(ORONE(dispatch_dim.x / 32),
ORONE(dispatch_dim.y),
ORONE(dispatch_dim.z / 32));
diff --git a/beamformer.h b/beamformer.h
@@ -170,7 +170,6 @@ typedef struct {
X(CS_DAS, volume_export_dim_offset) \
X(CS_DAS, volume_export_pass) \
X(CS_DAS, xdc_index) \
- X(CS_DAS, xdc_transform) \
X(CS_DAS, cycle_t) \
X(CS_MIN_MAX, mips_level) \
X(CS_SUM, sum_prescale)
diff --git a/beamformer_parameters.h b/beamformer_parameters.h
@@ -30,13 +30,12 @@ typedef struct {
u32 uforces_channels[128]; /* Channels used for virtual UFORCES elements */
f32 focal_depths[128]; /* [m] Focal Depths for each transmit of a RCA imaging scheme*/
f32 transmit_angles[128]; /* [radians] Transmit Angles for each transmit of a RCA imaging scheme*/
- f32 xdc_origin[4 * MAX_MULTI_XDC_COUNT]; /* [m] Corner of transducer being treated as origin */
- f32 xdc_corner1[4 * MAX_MULTI_XDC_COUNT]; /* [m] Corner of transducer along first axis */
- f32 xdc_corner2[4 * MAX_MULTI_XDC_COUNT]; /* [m] Corner of transducer along second axis */
+ f32 xdc_transforms[16 * MAX_MULTI_XDC_COUNT]; /* IMPORTANT: column major order */
uv4 dec_data_dim; /* Samples * Channels * Acquisitions; last element ignored */
uv4 output_points; /* Width * Height * Depth * (Frame Average Count) */
v4 output_min_coordinate; /* [m] Back-Top-Left corner of output region (w ignored) */
v4 output_max_coordinate; /* [m] Front-Bottom-Right corner of output region (w ignored)*/
+ v2 xdc_element_pitch; /* [m] Transducer Element Pitch */
uv2 rf_raw_dim; /* Raw Data Dimensions */
u32 decode; /* Decode or just reshape data */
u32 xdc_count; /* Number of Transducer Arrays (4 max) */
@@ -49,7 +48,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[3];
+ f32 _pad[1];
} BeamformerParameters;
/* NOTE: garbage to get the prepocessor to properly stringize the value of a macro */
@@ -64,13 +63,12 @@ layout(std140, binding = 0) uniform parameters {\n\
uvec4 uforces_channels[32]; /* Channels used for virtual UFORCES elements */\n\
vec4 focal_depths[32]; /* [m] Focal Depths for each transmit of a RCA imaging scheme*/\n\
vec4 transmit_angles[32]; /* [radians] Transmit Angles for each transmit of a RCA imaging scheme*/\n\
- vec4 xdc_origin[" str(MAX_MULTI_XDC_COUNT) "]; /* [m] Corner of transducer being treated as origin */\n\
- vec4 xdc_corner1[" str(MAX_MULTI_XDC_COUNT) "]; /* [m] Corner of transducer along first axis (arbitrary) */\n\
- vec4 xdc_corner2[" str(MAX_MULTI_XDC_COUNT) "]; /* [m] Corner of transducer along second axis (arbitrary) */\n\
+ mat4 xdc_transforms[" str(MAX_MULTI_XDC_COUNT) "]; /* IMPORTANT: column major order */\n\
uvec4 dec_data_dim; /* Samples * Channels * Acquisitions; last element ignored */\n\
uvec4 output_points; /* Width * Height * Depth * (Frame Average Count) */\n\
vec4 output_min_coord; /* [m] Top left corner of output region */\n\
vec4 output_max_coord; /* [m] Bottom right corner of output region */\n\
+ vec2 xdc_element_pitch; /* [m] Transducer Element Pitch */\n\
uvec2 rf_raw_dim; /* Raw Data Dimensions */\n\
uint decode; /* Decode or just reshape data */\n\
uint xdc_count; /* Number of Transducer Arrays (4 max) */\n\
diff --git a/shaders/das.glsl b/shaders/das.glsl
@@ -9,9 +9,8 @@ layout(rg32f, binding = 0) writeonly uniform image3D u_out_data_tex;
layout(location = 2) uniform int u_volume_export_pass;
layout(location = 3) uniform ivec3 u_volume_export_dim_offset;
-layout(location = 4) uniform mat4 u_xdc_transform;
-layout(location = 5) uniform int u_xdc_index;
-layout(location = 6) uniform float u_cycle_t;
+layout(location = 4) uniform int u_xdc_index;
+layout(location = 5) uniform float u_cycle_t;
#define C_SPLINE 0.5
@@ -78,14 +77,16 @@ vec2 apodize(vec2 value, float apodization_arg, float distance)
return value * a * a;
}
-vec3 orientation_project(bool tx_rows)
+vec3 orientation_projection(bool tx_rows)
{
return vec3(float(!tx_rows), float(tx_rows), 1);
}
vec3 world_space_to_rca_space(vec3 image_point, int transmit_orientation)
{
- return (u_xdc_transform * vec4(image_point, 1)).xyz * orientation_project(transmit_orientation != TX_ROWS);
+ vec3 result = (xdc_transforms[u_xdc_index] * vec4(image_point, 1)).xyz;
+ result *= orientation_projection(transmit_orientation != TX_ROWS);
+ return result;
}
float sample_index(float distance)
@@ -96,14 +97,14 @@ float sample_index(float distance)
float planewave_transmit_distance(vec3 point, float transmit_angle, int transmit_orientation)
{
- return dot(point * orientation_project(transmit_orientation == TX_ROWS),
+ return dot(point * orientation_projection(transmit_orientation == 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)
{
vec3 f = focal_depth * vec3(sin(transmit_angle), sin(transmit_angle), cos(transmit_angle));
- return length((point - f) * orientation_project(transmit_orientation == TX_ROWS));
+ return length((point - f) * orientation_projection(transmit_orientation == TX_ROWS));
}
vec2 RCA(vec3 image_point, vec3 delta, uint starting_offset, float apodization_arg)
@@ -162,7 +163,10 @@ vec2 HERCULES(vec3 image_point, vec3 delta, uint starting_offset, float apodizat
float focal_depth = focal_depths[0][0];
float transmit_angle = transmit_angles[0][0];
- vec3 recieve_point = (u_xdc_transform * vec4(image_point, 1)).xyz;
+ vec3 recieve_point = (xdc_transforms[u_xdc_index] * vec4(image_point, 1)).xyz;
+
+ if (transmit_orientation == TX_ROWS)
+ delta = delta.yxz;
float transmit_distance;
if (isinf(focal_depth)) {
@@ -180,7 +184,7 @@ vec2 HERCULES(vec3 image_point, vec3 delta, uint starting_offset, float apodizat
/* 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.yxz;
+ if (transmit_orientation == TX_ROWS) 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));
@@ -196,13 +200,21 @@ vec2 HERCULES(vec3 image_point, vec3 delta, uint starting_offset, float apodizat
return sum;
}
-vec2 uFORCES(vec3 image_point, vec3 delta, float y_off, uint starting_offset, float apodization_arg)
+vec2 uFORCES(vec3 image_point, vec3 delta, uint starting_offset, float apodization_arg)
{
/* NOTE: skip first acquisition in uforces since its garbage */
uint uforces = uint(dec_data_dim.y != dec_data_dim.z);
uint ridx = starting_offset + dec_data_dim.y * dec_data_dim.x * uforces;
- image_point = (u_xdc_transform * vec4(image_point, 1)).xyz;
+ image_point = (xdc_transforms[u_xdc_index] * 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;
+
+ vec3 focal_point_offset = vec3(0, delta.y * floor(dec_data_dim.y / 2), 0);
+ delta.y = 0;
vec2 sum = vec2(0);
for (uint i = uforces; i < dec_data_dim.z; i++) {
@@ -210,7 +222,7 @@ vec2 uFORCES(vec3 image_point, vec3 delta, float y_off, uint starting_offset, fl
uint sub_idx = (i - uforces) % 4;
vec2 rdist = vec2(image_point.x, image_point.z);
- vec3 focal_point = uforces_channels[base_idx][sub_idx] * delta + vec3(0, y_off, 0);
+ vec3 focal_point = uforces_channels[base_idx][sub_idx] * delta + focal_point_offset;
float transmit_dist = distance(image_point, focal_point);
for (uint j = 0; j < dec_data_dim.y; j++) {
@@ -244,40 +256,16 @@ void main()
/* NOTE: skip over channels corresponding to other arrays */
uint starting_offset = u_xdc_index * (dec_data_dim.y / xdc_count) * dec_data_dim.x * dec_data_dim.z;
- /* NOTE: array edge vectors for calculating element step delta */
- vec3 edge1 = xdc_corner1[u_xdc_index].xyz - xdc_origin[u_xdc_index].xyz;
- vec3 edge2 = xdc_corner2[u_xdc_index].xyz - xdc_origin[u_xdc_index].xyz;
- vec3 delta;
-
vec2 sum;
switch (das_shader_id) {
case DAS_ID_UFORCES:
- /* TODO: there should be a smarter way of detecting this.
- * (hack since we assume we beamform in x-z even when we are actually doing
- * y-z. this will be cleaned up by using an initial transform to shift the
- * voxel into the correct space) */
- float y_off;
- if (edge2.x != 0) {
- delta = vec3(edge2.x, 0, 0) / float(dec_data_dim.y - 1);
- y_off = edge1.y / 2;
- } else {
- delta = vec3(edge1.x, 0, 0) / float(dec_data_dim.y - 1);
- y_off = edge2.y / 2;
- }
-
- sum = uFORCES(image_point, delta, y_off, starting_offset, apod_arg);
+ sum = uFORCES(image_point, vec3(xdc_element_pitch, 0), starting_offset, apod_arg);
break;
case DAS_ID_HERCULES:
- /* TODO: there should be a smarter way of detecting this */
- if (edge2.x != 0) delta = vec3(edge2.x, edge1.y, 0) / float(dec_data_dim.y - 1);
- else delta = vec3(edge1.x, edge2.y, 0) / float(dec_data_dim.y - 1);
- sum = HERCULES(image_point, delta, starting_offset, apod_arg);
+ sum = HERCULES(image_point, vec3(xdc_element_pitch, 0), starting_offset, apod_arg);
break;
case DAS_ID_RCA:
- /* TODO: there should be a smarter way of detecting this */
- if (edge2.x != 0) delta = vec3(edge2.x, edge1.y, 0) / float(dec_data_dim.y - 1);
- else delta = vec3(edge1.x, edge2.y, 0) / float(dec_data_dim.y - 1);
- sum = RCA(image_point, delta, starting_offset, apod_arg);
+ sum = RCA(image_point, vec3(xdc_element_pitch, 0), starting_offset, apod_arg);
break;
}