Commit: 0171f2ba6acb39c216654b82bbdf08fa04e4200e
Parent: 2317f027275a706ac57cb1dee459d4ccd997b49b
Author: Randy Palamar
Date: Thu, 14 Nov 2024 18:07:59 -0700
correct transducer coordinate transform calculation
This calculation should give the correct transform regardless of
what order the corners are specified in. It should also include
the translation because it is free to do so.
Diffstat:
4 files changed, 43 insertions(+), 43 deletions(-)
diff --git a/beamformer.c b/beamformer.c
@@ -231,23 +231,28 @@ beamform_work_queue_push(BeamformerCtx *ctx, Arena *a, enum beamform_work work_t
return result;
}
-static m3
+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);
- xdc_normal.z = ABS(xdc_normal.z);
+ if (xdc_normal.z < 0)
+ xdc_normal = cross(edge2, edge1);
+ ASSERT(xdc_normal.z >= 0);
- v3 e1 = normalize_v3(sub_v3(xdc_normal, direction));
+ 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};
- m3 result = {
- .c[0] = (v3){.x = e3.x, .y = e2.x, .z = e1.x},
- .c[1] = (v3){.x = e3.y, .y = e2.y, .z = e1.y},
- .c[2] = (v3){.x = e3.z, .y = e2.z, .z = e1.z},
+ 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;
}
@@ -301,13 +306,13 @@ do_beamform_shader(ComputeShaderCtx *cs, BeamformerParameters *bp, BeamformFrame
for (u32 i = 0; i < frame->dim.w; i++) {
u32 texture = frame->textures[i];
- m3 xdc_transform = v3_to_xdc_space((v3){.z = 1},
+ 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);
- glUniformMatrix3fv(cs->xdc_transform_id, 1, GL_FALSE, xdc_transform.E);
+ 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/shaders/hercules.glsl b/shaders/hercules.glsl
@@ -9,7 +9,7 @@ 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 mat3 u_xdc_transform;
+layout(location = 4) uniform mat4 u_xdc_transform;
layout(location = 5) uniform int u_xdc_index;
#define C_SPLINE 0.5
@@ -53,25 +53,30 @@ vec3 calc_image_point(vec3 voxel)
{
ivec3 out_data_dim = imageSize(u_out_data_tex);
vec4 output_size = abs(output_max_coord - output_min_coord);
- vec3 image_point = output_min_coord.xyz + voxel * output_size.xyz / out_data_dim.xyz;
+ vec4 image_point = vec4(output_min_coord.xyz + voxel * output_size.xyz / out_data_dim, 1);
if (u_volume_export_pass == 0)
image_point.y = off_axis_pos;
/* NOTE: move the image point into xdc space */
image_point = u_xdc_transform * image_point;
- return image_point;
+
+ return image_point.xyz;
}
void main()
{
- vec3 voxel = vec3(gl_GlobalInvocationID.xyz) + vec3(u_volume_export_dim_offset);
- ivec3 out_coord = ivec3(gl_GlobalInvocationID.xyz) + u_volume_export_dim_offset;
+ vec3 voxel = vec3(gl_GlobalInvocationID.xyz) + vec3(u_volume_export_dim_offset);
+ ivec3 out_coord = ivec3(gl_GlobalInvocationID.xyz) + u_volume_export_dim_offset;
/* NOTE: Convert voxel to physical coordinates */
- 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 image_point = calc_image_point(voxel);
+ 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 image_point = calc_image_point(voxel);
+ vec3 delta;
+ /* 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);
/* NOTE: used for constant F# dynamic receive apodization. This is implemented as:
*
@@ -83,15 +88,8 @@ void main()
float f_num = 0.5;
float apod_arg = f_num * 0.5 * radians(360) / abs(image_point.z);
- /* NOTE: lerp along a line from one edge of the xdc to the other in the imaging plane */
- vec3 delta = vec3(edge1.x, edge2.y, 0) / float(dec_data_dim.y);
- vec3 xdc_start = xdc_origin[u_xdc_index].xyz;
- xdc_start += edge2 / 2;
-
- vec3 starting_point = image_point - xdc_start;
-
vec2 sum = vec2(0);
- vec3 rdist = starting_point;
+ vec3 rdist = image_point;
/* TODO: pass this in (there is a problem in that it depends on the orientation
* of the array relative to the target/subject). */
@@ -127,13 +125,13 @@ void main()
vec2 p = cubic(ridx, sidx);
/* NOTE: tribal knowledge; this is a problem with the imaging sequence */
if (i == 0) p *= inversesqrt(128);
- sum += p;
+ sum += p * a;
rdist[direction] -= delta[direction];
ridx += dec_data_dim.x;
}
- rdist[direction] = starting_point[direction];
+ rdist[direction] = image_point[direction];
rdist[direction ^ 1] -= delta[direction ^ 1];
}
imageStore(u_out_data_tex, out_coord, vec4(sum.x, sum.y, 0, 0));
diff --git a/shaders/uforces.glsl b/shaders/uforces.glsl
@@ -9,7 +9,7 @@ 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 mat3 u_xdc_transform;
+layout(location = 4) uniform mat4 u_xdc_transform;
layout(location = 5) uniform int u_xdc_index;
#define C_SPLINE 0.5
@@ -50,14 +50,14 @@ vec3 calc_image_point(vec3 voxel)
{
ivec3 out_data_dim = imageSize(u_out_data_tex);
vec4 output_size = abs(output_max_coord - output_min_coord);
- vec3 image_point = output_min_coord.xyz + voxel * output_size.xyz / out_data_dim.xyz;
+ vec4 image_point = vec4(output_min_coord.xyz + voxel * output_size.xyz / out_data_dim, 1);
/* TODO: fix the math so that the image plane can be aritrary */
image_point.y = 0;
/* NOTE: move the image point into xdc space */
image_point = u_xdc_transform * image_point;
- return image_point;
+ return image_point.xyz;
}
void main()
@@ -69,6 +69,10 @@ void main()
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 image_point = calc_image_point(voxel);
+ vec3 delta;
+ /* TODO: there should be a smarter way of detecting this */
+ if (edge2.x != 0) delta = vec3(edge2.x, 0, 0) / float(dec_data_dim.y);
+ else delta = vec3(edge1.x, 0, 0) / float(dec_data_dim.y);
/* NOTE: used for constant F# dynamic receive apodization. This is implemented as:
*
@@ -80,13 +84,6 @@ void main()
float f_num = 0.5;
float apod_arg = f_num * 0.5 * radians(360) / abs(image_point.z);
- /* NOTE: lerp along a line from one edge of the xdc to the other in the imaging plane */
- vec3 delta = edge1 / float(dec_data_dim.y);
- vec3 xdc_start = xdc_origin[u_xdc_index].xyz;
- xdc_start += edge2 / 2;
-
- vec3 starting_point = image_point - xdc_start;
-
vec2 sum = vec2(0);
/* NOTE: skip over channels corresponding to other arrays */
uint ridx = u_xdc_index * (dec_data_dim.y / xdc_count) * dec_data_dim.x * dec_data_dim.z;
@@ -97,9 +94,9 @@ void main()
uint base_idx = (i - uforces) / 4;
uint sub_idx = (i - uforces) % 4;
- vec3 focal_point = uforces_channels[base_idx][sub_idx] * delta + xdc_start;
+ vec3 focal_point = uforces_channels[base_idx][sub_idx] * delta;
float transmit_dist = distance(image_point, focal_point);
- vec3 rdist = starting_point;
+ vec3 rdist = image_point;
for (uint j = 0; j < dec_data_dim.y; j++) {
float dist = transmit_dist + length(rdist);
float time = dist / speed_of_sound + time_offset;
diff --git a/util.h b/util.h
@@ -134,10 +134,10 @@ typedef union {
} v4;
typedef union {
- struct { v3 x, y, z; };
- v3 c[3];
- f32 E[9];
-} m3;
+ struct { v4 x, y, z, w; };
+ v4 c[4];
+ f32 E[16];
+} m4;
typedef union {
struct { v2 pos, size; };