hercules.glsl (6461B)
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) uniform writeonly 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 #define TX_ROWS 0 42 #define TX_COLS 1 43 44 #if 1 45 /* NOTE: interpolation is unnecessary if the data has been demodulated and not decimated */ 46 vec2 cubic(uint ridx, float t) 47 { 48 return rf_data[ridx + uint(floor(t))]; 49 } 50 #else 51 /* NOTE: See: https://cubic.org/docs/hermite.htm */ 52 vec2 cubic(uint ridx, float x) 53 { 54 mat4 h = mat4( 55 2, -3, 0, 1, 56 -2, 3, 0, 0, 57 1, -2, 1, 0, 58 1, -1, 0, 0 59 ); 60 61 uint xk = uint(floor(x)); 62 float t = (x - float(xk)); 63 vec4 S = vec4(t * t * t, t * t, t, 1); 64 65 vec2 P1 = rf_data[ridx + xk]; 66 vec2 P2 = rf_data[ridx + xk + 1]; 67 vec2 T1 = C_SPLINE * (P2 - rf_data[ridx + xk - 1]); 68 vec2 T2 = C_SPLINE * (rf_data[ridx + xk + 2] - P1); 69 70 vec4 C1 = vec4(P1.x, P2.x, T1.x, T2.x); 71 vec4 C2 = vec4(P1.y, P2.y, T1.y, T2.y); 72 return vec2(dot(S, h * C1), dot(S, h * C2)); 73 } 74 #endif 75 76 void main() 77 { 78 vec3 voxel = vec3(gl_GlobalInvocationID.xyz) + vec3(u_volume_export_dim_offset); 79 ivec3 out_coord = ivec3(gl_GlobalInvocationID.xyz) + u_volume_export_dim_offset; 80 ivec3 out_data_dim; 81 82 if (u_volume_export_pass == 0) out_data_dim = imageSize(u_out_data_tex); 83 else out_data_dim = imageSize(u_out_volume_tex); 84 85 /* NOTE: Convert pixel to physical coordinates */ 86 vec4 xdc_size = xdc_corner1 + xdc_corner2 - xdc_origin; 87 vec3 edge1 = xdc_corner1.xyz - xdc_origin.xyz; 88 vec3 edge2 = xdc_corner2.xyz - xdc_origin.xyz; 89 vec3 xdc_normal = cross(edge2, edge1); 90 xdc_normal /= length(xdc_normal); 91 vec4 output_size = abs(output_max_coord - output_min_coord); 92 vec3 image_point = output_min_coord.xyz + voxel * output_size.xyz / out_data_dim.xyz; 93 94 if (u_volume_export_pass == 0) 95 image_point.y = off_axis_pos; 96 97 /* NOTE: used for constant F# dynamic receive apodization. This is implemented as: 98 * 99 * / |x_e - x_i|\ 100 * a(x, z) = cos(F# * π * ----------- ) ^ 2 101 * \ |z_e - z_i|/ 102 * 103 * where x,z_e are transducer element positions and x,z_i are image positions. */ 104 float f_num = output_size.z / output_size.x; 105 float apod_arg = f_num * 0.5 * radians(360) / abs(image_point.z); 106 107 /* NOTE: for I-Q data phase correction */ 108 float iq_time_scale = (lpf_order > 0)? radians(360) * center_frequency : 0; 109 110 /* NOTE: lerp along a line from one edge of the xdc to the other in the imaging plane */ 111 vec3 delta = edge1 / float(dec_data_dim.y); 112 vec3 xdc_start = xdc_origin.xyz; 113 xdc_start += edge2 / 2; 114 115 vec3 starting_point = image_point - xdc_start; 116 117 /* NOTE: offset correcting for both pulse length and low pass filtering */ 118 float time_correction = time_offset + lpf_order / sampling_frequency; 119 120 vec2 sum = vec2(0); 121 vec3 rdist = starting_point; 122 123 /* TODO: pass this in (there is a problem in that it depends on the orientation 124 * of the array relative to the target/subject). */ 125 int transmit_orientation = TX_ROWS; 126 float transmit_dist; 127 if (isinf(focal_depth)) { 128 /* NOTE: plane wave */ 129 transmit_dist = image_point.z; 130 } else { 131 /* NOTE: cylindrical diverging wave */ 132 if (transmit_orientation == TX_ROWS) 133 transmit_dist = length(vec2(image_point.y, image_point.z - focal_depth)); 134 else 135 transmit_dist = length(vec2(image_point.x, image_point.z - focal_depth)); 136 } 137 138 int direction = beamform_plane * (u_volume_export_pass ^ 1); 139 uint ridx = 0; 140 /* NOTE: For Each Acquistion in Raw Data */ 141 for (uint i = 0; i < dec_data_dim.z; i++) { 142 /* NOTE: For Each Virtual Source */ 143 for (uint j = 0; j < dec_data_dim.y; j++) { 144 float dist = transmit_dist + length(rdist); 145 float time = dist / speed_of_sound + time_correction; 146 147 /* NOTE: apodization value for this transducer element */ 148 float a = cos(clamp(abs(apod_arg * rdist.x), 0, 0.25 * radians(360))); 149 a = a * a; 150 151 vec2 p = cubic(ridx, time * sampling_frequency); 152 /* NOTE: tribal knowledge; this is a problem with the imaging sequence */ 153 if (i == 0) p *= inversesqrt(128); 154 //p *= vec2(cos(iq_time_scale * time), sin(iq_time_scale * time)); 155 sum += p; 156 157 rdist[direction] -= delta[direction]; 158 ridx += dec_data_dim.x; 159 } 160 161 rdist[direction] = starting_point[direction]; 162 rdist[direction ^ 1] -= delta[direction ^ 1]; 163 } 164 float val = length(sum); 165 if (u_volume_export_pass == 0) imageStore(u_out_data_tex, out_coord, vec4(val)); 166 else imageStore(u_out_volume_tex, out_coord, vec4(val)); 167 }