uforces.glsl (3931B)
1 /* See LICENSE for license details. */ 2 layout(local_size_x = 32, local_size_y = 1, local_size_z = 32) in; 3 4 layout(std430, binding = 1) readonly restrict buffer buffer_1 { 5 vec2 rf_data[]; 6 }; 7 8 layout(rg32f, binding = 0) writeonly uniform image3D u_out_data_tex; 9 10 layout(location = 2) uniform int u_volume_export_pass; 11 layout(location = 3) uniform ivec3 u_volume_export_dim_offset; 12 layout(location = 4) uniform mat4 u_xdc_transform; 13 layout(location = 5) uniform int u_xdc_index; 14 15 #define C_SPLINE 0.5 16 17 #if 1 18 /* NOTE: interpolation is unnecessary if the data has been demodulated and not decimated */ 19 vec2 cubic(uint ridx, float t) 20 { 21 return rf_data[ridx + uint(floor(t))]; 22 } 23 #else 24 /* NOTE: See: https://cubic.org/docs/hermite.htm */ 25 vec2 cubic(uint ridx, float x) 26 { 27 mat4 h = mat4( 28 2, -3, 0, 1, 29 -2, 3, 0, 0, 30 1, -2, 1, 0, 31 1, -1, 0, 0 32 ); 33 34 uint xk = uint(floor(x)); 35 float t = (x - float(xk)); 36 vec4 S = vec4(t * t * t, t * t, t, 1); 37 38 vec2 P1 = rf_data[ridx + xk]; 39 vec2 P2 = rf_data[ridx + xk + 1]; 40 vec2 T1 = C_SPLINE * (P2 - rf_data[ridx + xk - 1]); 41 vec2 T2 = C_SPLINE * (rf_data[ridx + xk + 2] - P1); 42 43 vec4 C1 = vec4(P1.x, P2.x, T1.x, T2.x); 44 vec4 C2 = vec4(P1.y, P2.y, T1.y, T2.y); 45 return vec2(dot(S, h * C1), dot(S, h * C2)); 46 } 47 #endif 48 49 vec3 calc_image_point(vec3 voxel) 50 { 51 ivec3 out_data_dim = imageSize(u_out_data_tex); 52 vec4 output_size = abs(output_max_coord - output_min_coord); 53 vec4 image_point = vec4(output_min_coord.xyz + voxel * output_size.xyz / out_data_dim, 1); 54 55 /* TODO: fix the math so that the image plane can be aritrary */ 56 image_point.y = 0; 57 58 /* NOTE: move the image point into xdc space */ 59 image_point = u_xdc_transform * image_point; 60 return image_point.xyz; 61 } 62 63 void main() 64 { 65 vec3 voxel = vec3(gl_GlobalInvocationID); 66 ivec3 out_coord = ivec3(gl_GlobalInvocationID); 67 68 /* NOTE: Convert voxel to physical coordinates */ 69 vec3 edge1 = xdc_corner1[u_xdc_index].xyz - xdc_origin[u_xdc_index].xyz; 70 vec3 edge2 = xdc_corner2[u_xdc_index].xyz - xdc_origin[u_xdc_index].xyz; 71 vec3 image_point = calc_image_point(voxel); 72 vec3 delta; 73 /* TODO: there should be a smarter way of detecting this */ 74 if (edge2.x != 0) delta = vec3(edge2.x, 0, 0) / float(dec_data_dim.y); 75 else delta = vec3(edge1.x, 0, 0) / float(dec_data_dim.y); 76 77 /* NOTE: used for constant F# dynamic receive apodization. This is implemented as: 78 * 79 * / |x_e - x_i|\ 80 * a(x, z) = cos(F# * π * ----------- ) ^ 2 81 * \ |z_e - z_i|/ 82 * 83 * where x,z_e are transducer element positions and x,z_i are image positions. */ 84 float apod_arg = f_number * 0.5 * radians(360) / abs(image_point.z); 85 86 vec2 sum = vec2(0); 87 /* NOTE: skip over channels corresponding to other arrays */ 88 uint ridx = u_xdc_index * (dec_data_dim.y / xdc_count) * dec_data_dim.x * dec_data_dim.z; 89 /* NOTE: skip first acquisition in uforces since its garbage */ 90 uint uforces = uint(dec_data_dim.y != dec_data_dim.z); 91 ridx += dec_data_dim.y * dec_data_dim.x * uforces; 92 for (uint i = uforces; i < dec_data_dim.z; i++) { 93 uint base_idx = (i - uforces) / 4; 94 uint sub_idx = (i - uforces) % 4; 95 96 vec3 focal_point = uforces_channels[base_idx][sub_idx] * delta; 97 float transmit_dist = distance(image_point, focal_point); 98 vec3 rdist = image_point; 99 for (uint j = 0; j < dec_data_dim.y; j++) { 100 float dist = transmit_dist + length(rdist); 101 float time = dist / speed_of_sound + time_offset; 102 103 /* NOTE: apodization value for this transducer element */ 104 float a = cos(clamp(abs(apod_arg * rdist.x), 0, 0.25 * radians(360))); 105 a = a * a; 106 107 float sidx = time * sampling_frequency; 108 vec2 valid = vec2(sidx < dec_data_dim.x); 109 vec2 p = cubic(ridx, sidx) * valid; 110 sum += p * a; 111 rdist -= delta; 112 ridx += dec_data_dim.x; 113 } 114 } 115 imageStore(u_out_data_tex, out_coord, vec4(sum.x, sum.y, 0, 0)); 116 }