ogl_beamforming

Ultrasound Beamforming Implemented with OpenGL
git clone anongit@rnpnr.xyz:ogl_beamforming.git
Log | Files | Refs | Feed | Submodules | LICENSE

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 }