2d_hercules.glsl (5458B)
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 uvec4 dec_data_dim; /* Samples * Channels * Acquisitions; last element ignored */ 14 uvec4 output_points; /* Width * Height * Depth; last element ignored */ 15 vec4 output_min_coord; /* [m] Top left corner of output region */ 16 vec4 output_max_coord; /* [m] Bottom right corner of output region */ 17 uvec2 rf_raw_dim; /* Raw Data Dimensions */ 18 vec2 xdc_min_xy; /* [m] Min center of transducer elements */ 19 vec2 xdc_max_xy; /* [m] Max center of transducer elements */ 20 uint channel_offset; /* Offset into channel_mapping: 0 or 128 (rows or columns) */ 21 uint lpf_order; /* Order of Low Pass Filter */ 22 float speed_of_sound; /* [m/s] */ 23 float sampling_frequency; /* [Hz] */ 24 float center_frequency; /* [Hz] */ 25 float focal_depth; /* [m] */ 26 float time_offset; /* pulse length correction time [s] */ 27 uint uforces; /* mode is UFORCES (1) or FORCES (0) */ 28 }; 29 30 layout(rg32f, location = 1) uniform writeonly image3D u_out_data_tex; 31 layout(r32f, location = 2) uniform writeonly image3D u_out_volume_tex; 32 33 layout(location = 3) uniform int u_volume_export_pass; 34 35 #define C_SPLINE 0.5 36 37 #if 0 38 /* NOTE: interpolation is unnecessary if the data has been demodulated and not decimated */ 39 vec2 cubic(uint ridx, float t) 40 { 41 return rf_data[ridx + uint(floor(t))]; 42 } 43 #else 44 /* NOTE: See: https://cubic.org/docs/hermite.htm */ 45 vec2 cubic(uint ridx, float x) 46 { 47 mat4 h = mat4( 48 2, -3, 0, 1, 49 -2, 3, 0, 0, 50 1, -2, 1, 0, 51 1, -1, 0, 0 52 ); 53 54 uint xk = uint(floor(x)); 55 float t = (x - float(xk)); 56 vec4 S = vec4(t * t * t, t * t, t, 1); 57 58 vec2 P1 = rf_data[ridx + xk]; 59 vec2 P2 = rf_data[ridx + xk + 1]; 60 vec2 T1 = C_SPLINE * (P2 - rf_data[ridx + xk - 1]); 61 vec2 T2 = C_SPLINE * (rf_data[ridx + xk + 2] - P1); 62 63 vec4 C1 = vec4(P1.x, P2.x, T1.x, T2.x); 64 vec4 C2 = vec4(P1.y, P2.y, T1.y, T2.y); 65 return vec2(dot(S, h * C1), dot(S, h * C2)); 66 } 67 #endif 68 69 void main() 70 { 71 vec3 voxel = vec3(gl_GlobalInvocationID.xyz); 72 ivec3 out_coord = ivec3(gl_GlobalInvocationID.xyz); 73 ivec3 out_data_dim; 74 if (u_volume_export_pass == 0) out_data_dim = imageSize(u_out_data_tex); 75 else out_data_dim = imageSize(u_out_volume_tex); 76 77 /* NOTE: Convert pixel to physical coordinates */ 78 vec2 xdc_size = abs(xdc_max_xy - xdc_min_xy); 79 vec4 output_size = abs(output_max_coord - output_min_coord); 80 vec3 image_point = vec3( 81 output_min_coord.x + voxel.x * output_size.x / out_data_dim.x, 82 output_min_coord.y + voxel.y * output_size.y / out_data_dim.y, 83 output_min_coord.z + voxel.z * output_size.z / out_data_dim.z 84 ); 85 86 /* NOTE: used for constant F# dynamic receive apodization. This is implemented as: 87 * 88 * / |x_e - x_i|\ 89 * a(x, z) = cos(F# * π * ----------- ) ^ 2 90 * \ |z_e - z_i|/ 91 * 92 * where x,z_e are transducer element positions and x,z_i are image positions. */ 93 float f_num = output_size.z / output_size.x; 94 float apod_arg = f_num * 0.5 * radians(360) / abs(image_point.z); 95 96 /* NOTE: for I-Q data phase correction */ 97 float iq_time_scale = (lpf_order > 0)? radians(360) * center_frequency : 0; 98 99 vec3 starting_dist = vec3(image_point.x - xdc_min_xy.x, image_point.y - xdc_min_xy.y, image_point.z); 100 vec3 delta = vec3(xdc_size.x / float(dec_data_dim.y), xdc_size.y / float(dec_data_dim.y), 0); 101 float dzsign = sign(image_point.z - focal_depth); 102 103 /* NOTE: offset correcting for both pulse length and low pass filtering */ 104 float time_correction = time_offset + (lpf_order + 1) / sampling_frequency; 105 106 vec2 sum = vec2(0); 107 vec3 rdist = starting_dist; 108 109 int direction = 1 * (u_volume_export_pass ^ 1); 110 uint ridx = 0; 111 /* NOTE: For Each Acquistion in Raw Data */ 112 for (uint i = 0; i < dec_data_dim.z; i++) { 113 uint base_idx = (i - uforces) / 4; 114 uint sub_idx = (i - uforces) % 4; 115 116 float transmit_dist = image_point.z; 117 118 /* NOTE: For Each Virtual Source */ 119 for (uint j = 0; j < dec_data_dim.y; j++) { 120 float dist = transmit_dist + length(rdist); 121 float time = dist / speed_of_sound + time_correction; 122 123 /* NOTE: apodization value for this transducer element */ 124 float a = cos(clamp(abs(apod_arg * rdist.x), 0, 0.25 * radians(360))); 125 a = a * a; 126 127 vec2 p = cubic(ridx, time * sampling_frequency); 128 /* NOTE: tribal knowledge; this is a problem with the imaging sequence */ 129 if (i == 0) p *= inversesqrt(128); 130 //p *= vec2(cos(iq_time_scale * time), sin(iq_time_scale * time)); 131 sum += p; 132 133 rdist[direction] -= delta[direction]; 134 ridx += dec_data_dim.x; 135 } 136 137 rdist[direction] = starting_dist[direction]; 138 rdist[direction ^ 1] -= delta[direction ^ 1]; 139 } 140 float val = length(sum); 141 if (u_volume_export_pass == 0) imageStore(u_out_data_tex, out_coord, vec4(val)); 142 else imageStore(u_out_volume_tex, out_coord, vec4(val)); 143 }