uforces.glsl (5548B)
1 /* See LICENSE for license details. */ 2 #version 460 core 3 layout(local_size_x = 32, local_size_y = 1, local_size_z = 32) in; 4 5 layout(std430, binding = 1) readonly restrict buffer buffer_1 { 6 vec2 rf_data[]; 7 }; 8 9 layout(std140, binding = 0) uniform parameters { 10 uvec4 channel_mapping[64]; /* Transducer Channel to Verasonics Channel */ 11 uvec4 uforces_channels[32]; /* Channels used for virtual UFORCES elements */ 12 vec4 lpf_coefficients[16]; /* Low Pass Filter Cofficients */ 13 vec4 xdc_origin; /* [m] Corner of transducer being treated as origin */ 14 vec4 xdc_corner1; /* [m] Corner of transducer along first axis (arbitrary) */ 15 vec4 xdc_corner2; /* [m] Corner of transducer along second axis (arbitrary) */ 16 uvec4 dec_data_dim; /* Samples * Channels * Acquisitions; last element ignored */ 17 uvec4 output_points; /* Width * Height * Depth; last element ignored */ 18 vec4 output_min_coord; /* [m] Top left corner of output region */ 19 vec4 output_max_coord; /* [m] Bottom right corner of output region */ 20 uvec2 rf_raw_dim; /* Raw Data Dimensions */ 21 uint channel_offset; /* Offset into channel_mapping: 0 or 128 (rows or columns) */ 22 uint lpf_order; /* Order of Low Pass Filter */ 23 float speed_of_sound; /* [m/s] */ 24 float sampling_frequency; /* [Hz] */ 25 float center_frequency; /* [Hz] */ 26 float focal_depth; /* [m] */ 27 float time_offset; /* pulse length correction time [s] */ 28 uint uforces; /* mode is UFORCES (1) or FORCES (0) */ 29 float off_axis_pos; /* [m] Position on screen normal to beamform in 2D HERCULES */ 30 int beamform_plane; /* Plane to Beamform in 2D HERCULES */ 31 }; 32 33 layout(rg32f, location = 1) writeonly uniform image3D u_out_data_tex; 34 layout(r32f, location = 2) uniform writeonly image3D u_out_volume_tex; 35 36 layout(location = 3) uniform int u_volume_export_pass; 37 layout(location = 4) uniform ivec3 u_volume_export_dim_offset; 38 39 #define C_SPLINE 0.5 40 41 #if 1 42 /* NOTE: interpolation is unnecessary if the data has been demodulated and not decimated */ 43 vec2 cubic(uint ridx, float t) 44 { 45 return rf_data[ridx + uint(floor(t))]; 46 } 47 #else 48 /* NOTE: See: https://cubic.org/docs/hermite.htm */ 49 vec2 cubic(uint ridx, float x) 50 { 51 mat4 h = mat4( 52 2, -3, 0, 1, 53 -2, 3, 0, 0, 54 1, -2, 1, 0, 55 1, -1, 0, 0 56 ); 57 58 uint xk = uint(floor(x)); 59 float t = (x - float(xk)); 60 vec4 S = vec4(t * t * t, t * t, t, 1); 61 62 vec2 P1 = rf_data[ridx + xk]; 63 vec2 P2 = rf_data[ridx + xk + 1]; 64 vec2 T1 = C_SPLINE * (P2 - rf_data[ridx + xk - 1]); 65 vec2 T2 = C_SPLINE * (rf_data[ridx + xk + 2] - P1); 66 67 vec4 C1 = vec4(P1.x, P2.x, T1.x, T2.x); 68 vec4 C2 = vec4(P1.y, P2.y, T1.y, T2.y); 69 return vec2(dot(S, h * C1), dot(S, h * C2)); 70 } 71 #endif 72 73 void main() 74 { 75 vec3 voxel = vec3(gl_GlobalInvocationID); 76 ivec3 out_coord = ivec3(gl_GlobalInvocationID); 77 ivec3 out_data_dim = imageSize(u_out_data_tex); 78 79 /* NOTE: Convert voxel to physical coordinates */ 80 vec4 xdc_size = xdc_corner1 + xdc_corner2 - xdc_origin; 81 vec3 edge1 = xdc_corner1.xyz - xdc_origin.xyz; 82 vec3 edge2 = xdc_corner2.xyz - xdc_origin.xyz; 83 vec3 xdc_normal = cross(edge2, edge1); 84 xdc_normal /= length(xdc_normal); 85 vec4 output_size = abs(output_max_coord - output_min_coord); 86 vec3 image_point = output_min_coord.xyz + voxel * output_size.xyz / out_data_dim.xyz; 87 88 /* TODO: fix the math so that the image plane can be aritrary */ 89 image_point.y = 0; 90 91 /* NOTE: used for constant F# dynamic receive apodization. This is implemented as: 92 * 93 * / |x_e - x_i|\ 94 * a(x, z) = cos(F# * π * ----------- ) ^ 2 95 * \ |z_e - z_i|/ 96 * 97 * where x,z_e are transducer element positions and x,z_i are image positions. */ 98 float f_num = 0.5; //output_size.z / output_size.x; 99 float apod_arg = f_num * 0.5 * radians(360) / abs(image_point.z); 100 101 /* NOTE: for I-Q data phase correction */ 102 float iq_time_scale = (lpf_order > 0)? radians(360) * center_frequency : 0; 103 104 /* NOTE: lerp along a line from one edge of the xdc to the other in the imaging plane */ 105 vec3 delta = edge1 / float(dec_data_dim.y); 106 vec3 xdc_start = xdc_origin.xyz; 107 xdc_start += edge2 / 2; 108 109 vec3 starting_point = image_point - xdc_start; 110 111 /* NOTE: offset correcting for both pulse length and low pass filtering */ 112 float time_correction = time_offset + lpf_order / sampling_frequency; 113 114 vec2 sum = vec2(0); 115 /* NOTE: skip first acquisition in uforces since its garbage */ 116 uint ridx = dec_data_dim.y * dec_data_dim.x * uforces; 117 for (uint i = uforces; i < dec_data_dim.z; i++) { 118 uint base_idx = (i - uforces) / 4; 119 uint sub_idx = (i - uforces) % 4; 120 121 vec3 focal_point = uforces_channels[base_idx][sub_idx] * delta + xdc_start; 122 float transmit_dist = distance(image_point, focal_point); 123 vec3 rdist = starting_point; 124 for (uint j = 0; j < dec_data_dim.y; j++) { 125 float dist = transmit_dist + length(rdist); 126 float time = dist / speed_of_sound + time_correction; 127 128 /* NOTE: apodization value for this transducer element */ 129 float a = cos(clamp(abs(apod_arg * rdist.x), 0, 0.25 * radians(360))); 130 a = a * a; 131 132 vec2 p = cubic(ridx, time * sampling_frequency); 133 p *= vec2(cos(iq_time_scale * time), sin(iq_time_scale * time)); 134 sum += p * a; 135 rdist -= delta; 136 ridx += dec_data_dim.x; 137 } 138 } 139 float val = length(sum); 140 imageStore(u_out_data_tex, out_coord, vec4(val, val, 0, 0)); 141 }