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