Commit: b3a86a2e37131f3aad0007b3cdd589d245f9ff6c
Parent: 309a9a4de94cac0cdbfa0e0b8964b2b24df78936
Author: zemplab
Date: Fri, 6 Dec 2024 16:17:43 -0700
add a couple normal RCA beamforming methods to das shader
Diffstat:
3 files changed, 92 insertions(+), 21 deletions(-)
diff --git a/beamformer_parameters.h b/beamformer_parameters.h
@@ -19,6 +19,7 @@ enum compute_shaders {
#define DAS_ID_UFORCES 0
#define DAS_ID_HERCULES 1
+#define DAS_ID_RCA 2
#define MAX_BEAMFORMED_SAVED_FRAMES 16
#define MAX_MULTI_XDC_COUNT 4
@@ -27,6 +28,8 @@ enum compute_shaders {
typedef struct {
u16 channel_mapping[512]; /* Transducer Channel to Verasonics Channel */
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 */
@@ -35,12 +38,12 @@ typedef struct {
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)*/
uv2 rf_raw_dim; /* Raw Data Dimensions */
+ u32 decode; /* Decode or just reshape data */
u32 xdc_count; /* Number of Transducer Arrays (4 max) */
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] */
- f32 focal_depth; /* [m] */
f32 time_offset; /* pulse length correction time [s] */
f32 off_axis_pos; /* [m] Position on screen normal to beamform in 2D HERCULES */
i32 beamform_plane; /* Plane to Beamform in 2D HERCULES */
@@ -59,6 +62,8 @@ typedef struct {
layout(std140, binding = 0) uniform parameters {\n\
uvec4 channel_mapping[64]; /* Transducer Channel to Verasonics Channel */\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\
@@ -67,12 +72,12 @@ layout(std140, binding = 0) uniform parameters {\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\
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\
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\
- float focal_depth; /* [m] */\n\
float time_offset; /* pulse length correction time [s] */\n\
float off_axis_pos; /* [m] Position on screen normal to beamform in 2D HERCULES */\n\
int beamform_plane; /* Plane to Beamform in 2D HERCULES */\n\
@@ -82,4 +87,5 @@ layout(std140, binding = 0) uniform parameters {\n\
\n\
#define DAS_ID_UFORCES " str(DAS_ID_UFORCES) "\n\
#define DAS_ID_HERCULES " str(DAS_ID_HERCULES) "\n\
-\n"
+#define DAS_ID_RCA " str(DAS_ID_RCA) "\n\n\
+#line 0\n"
diff --git a/shaders/das.glsl b/shaders/das.glsl
@@ -62,6 +62,7 @@ vec3 calc_image_point(vec3 voxel)
image_point.y = 0;
break;
case DAS_ID_HERCULES:
+ case DAS_ID_RCA:
if (u_volume_export_pass == 0)
image_point.y = off_axis_pos;
break;
@@ -82,30 +83,85 @@ vec3 row_column_point_scale(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 * row_column_point_scale(transmit_orientation != TX_ROWS);
+}
+
float sample_index(float distance)
{
float time = distance / speed_of_sound + time_offset;
return time * sampling_frequency;
}
+float planewave_transmit_distance(vec3 point, float transmit_angle)
+{
+ return dot(point, vec3(sin(transmit_angle), sin(transmit_angle), cos(transmit_angle)));
+}
+
+vec2 RCA(vec3 image_point, vec3 delta, uint starting_offset, 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). */
+ int transmit_orientation = TX_ROWS;
+ uint ridx = starting_offset;
+ int direction = beamform_plane * (u_volume_export_pass ^ 1);
+ if (direction == TX_COLS) image_point = image_point.yxz;
+
+ vec3 transmit_point = image_point * row_column_point_scale(transmit_orientation == TX_ROWS);
+ vec3 recieve_point = world_space_to_rca_space(image_point, transmit_orientation);
+ // vec3 recieve_point = (u_xdc_transform * vec4(image_point, 1)).xyz;
+
+ vec2 sum = vec2(0);
+ /* NOTE: For Each Acquistion in Raw Data */
+ // uint i = (dec_data_dim.z - 1) * uint(clamp(u_cycle_t, 0, 1)); {
+ for (uint i = 0; i < dec_data_dim.z; i++) {
+ uint base_idx = i / 4;
+ uint sub_idx = i % 4;
+
+ float focal_depth = focal_depths[base_idx][sub_idx];
+ float transmit_angle = transmit_angles[base_idx][sub_idx];
+ float tdist;
+ if (isinf(focal_depth)) {
+ tdist = planewave_transmit_distance(transmit_point, transmit_angle);
+ } else {
+ vec3 f = vec3(sin(transmit_angle), sin(transmit_angle), cos(transmit_angle));
+ f *= focal_depth * row_column_point_scale(transmit_orientation == TX_ROWS);
+ tdist = distance(transmit_point, f);
+ }
+
+ vec3 rdist = recieve_point;
+ /* NOTE: For Each Receiver */
+ // uint j = (dec_data_dim.z - 1) * uint(clamp(u_cycle_t, 0, 1)); {
+ for (uint j = 0; j < dec_data_dim.y; j++) {
+ float sidx = sample_index(tdist + length(rdist));
+ vec2 valid = vec2(sidx < dec_data_dim.x);
+ sum += apodize(cubic(ridx, sidx), apodization_arg, rdist[0]) * valid;
+ rdist[0] -= delta[0];
+ ridx += dec_data_dim.x;
+ }
+ }
+ return sum;
+}
+
vec2 HERCULES(vec3 image_point, vec3 delta, uint starting_offset, 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). */
int transmit_orientation = TX_ROWS;
- float transmit_dist;
+ float focal_depth = focal_depths[0][0];
+ float transmit_angle = transmit_angles[0][0];
vec3 transmit_point = image_point * row_column_point_scale(transmit_orientation == TX_ROWS);
+ //vec3 recieve_point = world_space_to_rca_space(image_point, transmit_orientation);
vec3 recieve_point = (u_xdc_transform * vec4(image_point, 1)).xyz;
+ float tdist;
if (isinf(focal_depth)) {
- /* NOTE: plane wave */
- transmit_dist = image_point.z;
+ tdist = planewave_transmit_distance(transmit_point, transmit_angle);
} else {
- /* NOTE: cylindrical diverging wave */
- if (transmit_orientation == TX_ROWS)
- transmit_dist = length(vec2(image_point.y, image_point.z - focal_depth));
- else
- transmit_dist = length(vec2(image_point.x, image_point.z - focal_depth));
+ vec3 f = vec3(sin(transmit_angle), sin(transmit_angle), cos(transmit_angle));
+ f *= focal_depth * row_column_point_scale(transmit_orientation == TX_ROWS);
+ tdist = distance(transmit_point, f);
}
uint ridx = starting_offset;
@@ -117,10 +173,8 @@ vec2 HERCULES(vec3 image_point, vec3 delta, uint starting_offset, float apodizat
for (uint i = 0; i < dec_data_dim.z; i++) {
/* NOTE: For Each Virtual Source */
for (uint j = 0; j < dec_data_dim.y; j++) {
- float sidx = sample_index(transmit_dist + length(rdist));
+ float sidx = sample_index(tdist + length(rdist));
vec2 valid = vec2(sidx < dec_data_dim.x);
- /* NOTE: tribal knowledge; this is a problem with the imaging sequence */
- if (i == 0) valid *= inversesqrt(128);
sum += apodize(cubic(ridx, sidx), apodization_arg, rdist.x) * valid;
@@ -202,6 +256,12 @@ void main()
else delta = vec3(edge1.x, edge2.y, 0) / float(dec_data_dim.y);
sum = HERCULES(image_point, delta, 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);
+ else delta = vec3(edge1.x, edge2.y, 0) / float(dec_data_dim.y);
+ sum = RCA(image_point, delta, starting_offset, apod_arg);
+ break;
}
imageStore(u_out_data_tex, out_coord, vec4(sum.x, sum.y, 0, 0));
diff --git a/shaders/hadamard.glsl b/shaders/hadamard.glsl
@@ -24,9 +24,6 @@ void main()
uint channel = gl_GlobalInvocationID.y;
uint acq = gl_GlobalInvocationID.z;
- /* NOTE: offset to get the correct column in hadamard matrix */
- uint hoff = dec_data_dim.z * acq;
-
/* 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;
@@ -53,10 +50,18 @@ void main()
/* NOTE: Compute N-D dot product */
int sum = 0;
- for (int i = 0; i < dec_data_dim.z; i++) {
- int data = (rf_data[ridx] << lfs) >> 16;
- sum += hadamard[hoff + i] * data;
- ridx += ridx_delta;
+ if (decode == 1) {
+ /* NOTE: offset to get the correct column in hadamard matrix */
+ uint hoff = dec_data_dim.z * acq;
+
+ for (int i = 0; i < dec_data_dim.z; i++) {
+ int data = (rf_data[ridx] << lfs) >> 16;
+ sum += hadamard[hoff + i] * data;
+ ridx += ridx_delta;
+ }
+ } else {
+ ridx += ridx_delta * acq;
+ sum = (rf_data[ridx] << lfs) >> 16;
}
out_data[out_off] = vec2(float(sum), 0);
}