das.glsl (8760B)
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 float u_cycle_t; 13 14 #define C_SPLINE 0.5 15 16 #define TX_ROWS 0 17 #define TX_COLS 1 18 19 #define TX_MODE_TX_COLS(a) (((a) & 2) != 0) 20 #define TX_MODE_RX_COLS(a) (((a) & 1) != 0) 21 22 #if 1 23 /* NOTE: interpolation is unnecessary if the data has been demodulated and not decimated */ 24 vec2 cubic(uint ridx, float t) 25 { 26 return rf_data[ridx + uint(floor(t))]; 27 } 28 #else 29 /* NOTE: See: https://cubic.org/docs/hermite.htm */ 30 vec2 cubic(uint ridx, float x) 31 { 32 mat4 h = mat4( 33 2, -3, 0, 1, 34 -2, 3, 0, 0, 35 1, -2, 1, 0, 36 1, -1, 0, 0 37 ); 38 39 uint xk = uint(floor(x)); 40 float t = (x - float(xk)); 41 vec4 S = vec4(t * t * t, t * t, t, 1); 42 43 vec2 P1 = rf_data[ridx + xk]; 44 vec2 P2 = rf_data[ridx + xk + 1]; 45 vec2 T1 = C_SPLINE * (P2 - rf_data[ridx + xk - 1]); 46 vec2 T2 = C_SPLINE * (rf_data[ridx + xk + 2] - P1); 47 48 vec4 C1 = vec4(P1.x, P2.x, T1.x, T2.x); 49 vec4 C2 = vec4(P1.y, P2.y, T1.y, T2.y); 50 return vec2(dot(S, h * C1), dot(S, h * C2)); 51 } 52 #endif 53 54 vec3 calc_image_point(vec3 voxel) 55 { 56 ivec3 out_data_dim = imageSize(u_out_data_tex); 57 vec4 output_size = abs(output_max_coord - output_min_coord); 58 vec3 image_point = output_min_coord.xyz + voxel * output_size.xyz / out_data_dim; 59 60 switch (das_shader_id) { 61 case DAS_ID_FORCES: 62 case DAS_ID_UFORCES: 63 /* TODO: fix the math so that the image plane can be aritrary */ 64 image_point.y = 0; 65 break; 66 case DAS_ID_HERCULES: 67 case DAS_ID_RCA_TPW: 68 case DAS_ID_RCA_VLS: 69 if (u_volume_export_pass == 0) 70 image_point.y = off_axis_pos; 71 break; 72 } 73 74 return image_point; 75 } 76 77 vec2 apodize(vec2 value, float apodization_arg, float distance) 78 { 79 /* NOTE: apodization value for this transducer element */ 80 float a = cos(clamp(abs(apodization_arg * distance), 0, 0.25 * radians(360))); 81 return value * a * a; 82 } 83 84 vec3 orientation_projection(vec3 point, bool rows) 85 { 86 return point * vec3(!rows, rows, 1); 87 } 88 89 vec3 world_space_to_rca_space(vec3 image_point, bool rx_rows) 90 { 91 return orientation_projection((xdc_transform * vec4(image_point, 1)).xyz, rx_rows); 92 } 93 94 float sample_index(float distance) 95 { 96 float time = distance / speed_of_sound + time_offset; 97 return time * sampling_frequency; 98 } 99 100 float planewave_transmit_distance(vec3 point, float transmit_angle, bool tx_rows) 101 { 102 return dot(orientation_projection(point, tx_rows), 103 vec3(sin(transmit_angle), sin(transmit_angle), cos(transmit_angle))); 104 } 105 106 float cylindricalwave_transmit_distance(vec3 point, float focal_depth, float transmit_angle, bool tx_rows) 107 { 108 vec3 f = focal_depth * vec3(sin(transmit_angle), sin(transmit_angle), cos(transmit_angle)); 109 return length(orientation_projection(point - f, tx_rows)); 110 } 111 112 vec2 RCA(vec3 image_point, vec3 delta, float apodization_arg) 113 { 114 uint ridx = 0; 115 int direction = beamform_plane * (u_volume_export_pass ^ 1); 116 if (direction != TX_ROWS) image_point = image_point.yxz; 117 118 bool tx_col = TX_MODE_TX_COLS(transmit_mode); 119 bool rx_col = TX_MODE_RX_COLS(transmit_mode); 120 121 vec3 receive_point = world_space_to_rca_space(image_point, !rx_col); 122 delta = orientation_projection(delta, !rx_col); 123 124 vec2 sum = vec2(0); 125 /* NOTE: For Each Acquistion in Raw Data */ 126 // uint i = (dec_data_dim.z - 1) * uint(clamp(u_cycle_t, 0, 1)); { 127 for (uint i = 0; i < dec_data_dim.z; i++) { 128 uint base_idx = i / 4; 129 uint sub_idx = i % 4; 130 131 float focal_depth = focal_depths[base_idx][sub_idx]; 132 float transmit_angle = transmit_angles[base_idx][sub_idx]; 133 134 float transmit_distance; 135 if (isinf(focal_depth)) { 136 transmit_distance = planewave_transmit_distance(image_point, transmit_angle, 137 !tx_col); 138 } else { 139 transmit_distance = cylindricalwave_transmit_distance(image_point, focal_depth, 140 transmit_angle, 141 !tx_col); 142 } 143 144 vec3 receive_distance = receive_point; 145 /* NOTE: For Each Receiver */ 146 // uint j = (dec_data_dim.z - 1) * uint(clamp(u_cycle_t, 0, 1)); { 147 for (uint j = 0; j < dec_data_dim.y; j++) { 148 float sidx = sample_index(transmit_distance + length(receive_distance)); 149 vec2 valid = vec2(sidx >= 0) * vec2(sidx < dec_data_dim.x); 150 sum += apodize(cubic(ridx, sidx), apodization_arg, length(receive_distance.xy)) * valid; 151 receive_distance -= delta; 152 ridx += dec_data_dim.x; 153 } 154 } 155 return sum; 156 } 157 158 vec2 HERCULES(vec3 image_point, vec3 delta, float apodization_arg) 159 { 160 uint ridx = 0; 161 int direction = beamform_plane * (u_volume_export_pass ^ 1); 162 if (direction != TX_ROWS) image_point = image_point.yxz; 163 164 bool tx_col = TX_MODE_TX_COLS(transmit_mode); 165 bool rx_col = TX_MODE_RX_COLS(transmit_mode); 166 float focal_depth = focal_depths[0][0]; 167 float transmit_angle = transmit_angles[0][0]; 168 169 vec3 receive_point = (xdc_transform * vec4(image_point, 1)).xyz; 170 171 float transmit_distance; 172 if (isinf(focal_depth)) { 173 transmit_distance = planewave_transmit_distance(image_point, transmit_angle, 174 !tx_col); 175 } else { 176 transmit_distance = cylindricalwave_transmit_distance(image_point, focal_depth, 177 transmit_angle, 178 !tx_col); 179 } 180 181 vec2 sum = vec2(0); 182 /* NOTE: For Each Acquistion in Raw Data */ 183 for (uint i = 0; i < dec_data_dim.z; i++) { 184 /* NOTE: For Each Virtual Source */ 185 for (uint j = 0; j < dec_data_dim.y; j++) { 186 vec3 element_position; 187 if (rx_col) element_position = vec3(j, i, 0) * delta; 188 else element_position = vec3(i, j, 0) * delta; 189 vec3 receive_distance = receive_point - element_position; 190 float sidx = sample_index(transmit_distance + length(receive_distance)); 191 vec2 valid = vec2(sidx >= 0) * vec2(sidx < dec_data_dim.x); 192 193 /* NOTE: tribal knowledge */ 194 if (i == 0) valid *= inversesqrt(dec_data_dim.z); 195 196 sum += apodize(cubic(ridx, sidx), apodization_arg, 197 length(receive_distance.xy)) * valid; 198 ridx += dec_data_dim.x; 199 } 200 } 201 return sum; 202 } 203 204 vec2 uFORCES(vec3 image_point, vec3 delta, float apodization_arg) 205 { 206 /* NOTE: skip first acquisition in uforces since its garbage */ 207 uint uforces = uint(das_shader_id == DAS_ID_UFORCES); 208 uint ridx = dec_data_dim.y * dec_data_dim.x * uforces; 209 210 image_point = (xdc_transform * vec4(image_point, 1)).xyz; 211 212 vec3 focal_point_offset = vec3(0, delta.y * floor(dec_data_dim.y / 2), 0); 213 delta.y = 0; 214 215 vec2 sum = vec2(0); 216 for (uint i = uforces; i < dec_data_dim.z; i++) { 217 uint base_idx = ((i - uforces) / 8); 218 uint sub_idx = ((i - uforces) % 8) / 2; 219 uint shift = (~(i - uforces) & 1u) * 16u; 220 uint channel = (uforces_channels[base_idx][sub_idx] << shift) >> 16u; 221 222 vec2 rdist = vec2(image_point.x, image_point.z); 223 vec3 focal_point = channel * delta + focal_point_offset; 224 float transmit_dist = distance(image_point, focal_point); 225 226 for (uint j = 0; j < dec_data_dim.y; j++) { 227 float sidx = sample_index(transmit_dist + length(rdist)); 228 vec2 valid = vec2(sidx >= 0) * vec2(sidx < dec_data_dim.x); 229 sum += apodize(cubic(ridx, sidx), apodization_arg, rdist.x) * valid; 230 rdist.x -= delta.x; 231 ridx += dec_data_dim.x; 232 } 233 } 234 return sum; 235 } 236 237 void main() 238 { 239 /* NOTE: Convert voxel to physical coordinates */ 240 ivec3 out_coord = ivec3(gl_GlobalInvocationID) + u_volume_export_dim_offset; 241 vec3 image_point = calc_image_point(vec3(gl_GlobalInvocationID) 242 + vec3(u_volume_export_dim_offset)); 243 244 /* NOTE: used for constant F# dynamic receive apodization. This is implemented as: 245 * 246 * / |x_e - x_i|\ 247 * a(x, z) = cos(F# * π * ----------- ) ^ 2 248 * \ |z_e - z_i|/ 249 * 250 * where x,z_e are transducer element positions and x,z_i are image positions. */ 251 float apod_arg = f_number * radians(180) / abs(image_point.z); 252 253 /* NOTE: skip over channels corresponding to other arrays */ 254 vec2 sum; 255 switch (das_shader_id) { 256 case DAS_ID_FORCES: 257 case DAS_ID_UFORCES: 258 sum = uFORCES(image_point, vec3(xdc_element_pitch, 0), apod_arg); 259 break; 260 case DAS_ID_HERCULES: 261 sum = HERCULES(image_point, vec3(xdc_element_pitch, 0), apod_arg); 262 break; 263 case DAS_ID_RCA_TPW: 264 case DAS_ID_RCA_VLS: 265 sum = RCA(image_point, vec3(xdc_element_pitch, 0), apod_arg); 266 break; 267 } 268 269 imageStore(u_out_data_tex, out_coord, vec4(sum, 0, 0)); 270 }