hercules.glsl (4544B)
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 #define TX_ROWS 0 18 #define TX_COLS 1 19 20 #if 1 21 /* NOTE: interpolation is unnecessary if the data has been demodulated and not decimated */ 22 vec2 cubic(uint ridx, float t) 23 { 24 return rf_data[ridx + uint(floor(t))]; 25 } 26 #else 27 /* NOTE: See: https://cubic.org/docs/hermite.htm */ 28 vec2 cubic(uint ridx, float x) 29 { 30 mat4 h = mat4( 31 2, -3, 0, 1, 32 -2, 3, 0, 0, 33 1, -2, 1, 0, 34 1, -1, 0, 0 35 ); 36 37 uint xk = uint(floor(x)); 38 float t = (x - float(xk)); 39 vec4 S = vec4(t * t * t, t * t, t, 1); 40 41 vec2 P1 = rf_data[ridx + xk]; 42 vec2 P2 = rf_data[ridx + xk + 1]; 43 vec2 T1 = C_SPLINE * (P2 - rf_data[ridx + xk - 1]); 44 vec2 T2 = C_SPLINE * (rf_data[ridx + xk + 2] - P1); 45 46 vec4 C1 = vec4(P1.x, P2.x, T1.x, T2.x); 47 vec4 C2 = vec4(P1.y, P2.y, T1.y, T2.y); 48 return vec2(dot(S, h * C1), dot(S, h * C2)); 49 } 50 #endif 51 52 vec3 calc_image_point(vec3 voxel) 53 { 54 ivec3 out_data_dim = imageSize(u_out_data_tex); 55 vec4 output_size = abs(output_max_coord - output_min_coord); 56 vec4 image_point = vec4(output_min_coord.xyz + voxel * output_size.xyz / out_data_dim, 1); 57 58 if (u_volume_export_pass == 0) 59 image_point.y = off_axis_pos; 60 61 /* NOTE: move the image point into xdc space */ 62 image_point = u_xdc_transform * image_point; 63 64 return image_point.xyz; 65 } 66 67 void main() 68 { 69 vec3 voxel = vec3(gl_GlobalInvocationID.xyz) + vec3(u_volume_export_dim_offset); 70 ivec3 out_coord = ivec3(gl_GlobalInvocationID.xyz) + u_volume_export_dim_offset; 71 72 /* NOTE: Convert voxel to physical coordinates */ 73 vec3 edge1 = xdc_corner1[u_xdc_index].xyz - xdc_origin[u_xdc_index].xyz; 74 vec3 edge2 = xdc_corner2[u_xdc_index].xyz - xdc_origin[u_xdc_index].xyz; 75 vec3 image_point = calc_image_point(voxel); 76 vec3 delta; 77 /* TODO: there should be a smarter way of detecting this */ 78 if (edge2.x != 0) delta = vec3(edge2.x, edge1.y, 0) / float(dec_data_dim.y); 79 else delta = vec3(edge1.x, edge2.y, 0) / float(dec_data_dim.y); 80 81 /* NOTE: used for constant F# dynamic receive apodization. This is implemented as: 82 * 83 * / |x_e - x_i|\ 84 * a(x, z) = cos(F# * π * ----------- ) ^ 2 85 * \ |z_e - z_i|/ 86 * 87 * where x,z_e are transducer element positions and x,z_i are image positions. */ 88 float apod_arg = f_number * 0.5 * radians(360) / abs(image_point.z); 89 90 vec2 sum = vec2(0); 91 vec3 rdist = image_point; 92 93 /* TODO: pass this in (there is a problem in that it depends on the orientation 94 * of the array relative to the target/subject). */ 95 int transmit_orientation = TX_ROWS; 96 float transmit_dist; 97 if (isinf(focal_depth)) { 98 /* NOTE: plane wave */ 99 transmit_dist = image_point.z; 100 } else { 101 /* NOTE: cylindrical diverging wave */ 102 if (transmit_orientation == TX_ROWS) 103 transmit_dist = length(vec2(image_point.y, image_point.z - focal_depth)); 104 else 105 transmit_dist = length(vec2(image_point.x, image_point.z - focal_depth)); 106 } 107 108 /* NOTE: skip over channels corresponding to other arrays */ 109 uint ridx = u_xdc_index * (dec_data_dim.y / xdc_count) * dec_data_dim.x * dec_data_dim.z; 110 int direction = beamform_plane * (u_volume_export_pass ^ 1); 111 /* NOTE: For Each Acquistion in Raw Data */ 112 for (uint i = 0; i < dec_data_dim.z; i++) { 113 /* NOTE: For Each Virtual Source */ 114 for (uint j = 0; j < dec_data_dim.y; j++) { 115 float dist = transmit_dist + length(rdist); 116 float time = dist / speed_of_sound + time_offset; 117 118 /* NOTE: apodization value for this transducer element */ 119 float a = cos(clamp(abs(apod_arg * rdist.x), 0, 0.25 * radians(360))); 120 a = a * a; 121 122 float sidx = time * sampling_frequency; 123 vec2 valid = vec2(sidx < dec_data_dim.x); 124 vec2 p = cubic(ridx, sidx); 125 /* NOTE: tribal knowledge; this is a problem with the imaging sequence */ 126 if (i == 0) p *= inversesqrt(128); 127 sum += p * a; 128 129 rdist[direction] -= delta[direction]; 130 ridx += dec_data_dim.x; 131 } 132 133 rdist[direction] = image_point[direction]; 134 rdist[direction ^ 1] -= delta[direction ^ 1]; 135 } 136 imageStore(u_out_data_tex, out_coord, vec4(sum.x, sum.y, 0, 0)); 137 }