ogl_beamforming

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

decode.glsl (3209B)


      1 /* See LICENSE for license details. */
      2 layout(local_size_x = 32, local_size_y = 32, local_size_z = 1) in;
      3 
      4 #if   defined(INPUT_DATA_TYPE_FLOAT)
      5 	#define INPUT_DATA_TYPE      float
      6 	#define RF_SAMPLES_PER_INDEX 1
      7 	#define RESULT_TYPE_CAST(x)  vec2(x, 0)
      8 #elif defined(INPUT_DATA_TYPE_FLOAT_COMPLEX)
      9 	#define INPUT_DATA_TYPE      vec2
     10 	#define RF_SAMPLES_PER_INDEX 1
     11 	#define RESULT_TYPE_CAST(x)  vec2(x)
     12 #else
     13 	#define INPUT_DATA_TYPE      int
     14 	#define RF_SAMPLES_PER_INDEX 2
     15 	#define RESULT_TYPE_CAST(x)  vec2(x, 0)
     16 #endif
     17 
     18 layout(std430, binding = 1) readonly restrict buffer buffer_1 {
     19 	INPUT_DATA_TYPE rf_data[];
     20 };
     21 
     22 layout(std430, binding = 2) writeonly restrict buffer buffer_2 {
     23 	vec2 out_data[];
     24 };
     25 
     26 layout(r8i, binding = 0) readonly restrict uniform iimage2D hadamard;
     27 
     28 INPUT_DATA_TYPE sample_rf_data(uint index, uint lfs)
     29 {
     30 	INPUT_DATA_TYPE result;
     31 #if   defined(INPUT_DATA_TYPE_FLOAT)
     32 	result = rf_data[index];
     33 #elif defined(INPUT_DATA_TYPE_FLOAT_COMPLEX)
     34 	result = rf_data[index];
     35 #else
     36 	result = (rf_data[index] << lfs) >> 16;
     37 #endif
     38 	return result;
     39 }
     40 
     41 void main()
     42 {
     43 	/* NOTE: each invocation takes a time sample and a receive channel.
     44 	 * It first maps the column to the correct column in the rf data then
     45 	 * does the dot product with the equivalent row of the hadamard matrix.
     46 	 * The result is stored to the equivalent row, column index of the output.
     47 	 */
     48 	uint time_sample = gl_GlobalInvocationID.x;
     49 	uint channel     = gl_GlobalInvocationID.y;
     50 	uint acq         = gl_GlobalInvocationID.z;
     51 
     52 	/* NOTE: offsets for storing the results in the output data */
     53 	uint out_off = dec_data_dim.x * dec_data_dim.y * acq + dec_data_dim.x * channel + time_sample;
     54 
     55 	/* NOTE: channel mapping is stored as u16s so we must do this to extract the final value */
     56 	uint ch_array_idx = (channel / 8);
     57 	uint ch_vec_idx   = (channel % 8) / 2;
     58 	uint ch_elem_lfs  = ((~channel) & 1u) * 16;
     59 	uint rf_channel   = (channel_mapping[ch_array_idx][ch_vec_idx] << ch_elem_lfs) >> 16;
     60 
     61 	/* NOTE: stride is the number of samples between acquistions; off is the
     62 	 * index of the first acquisition for this channel and time sample  */
     63 	uint rf_stride = dec_data_dim.x;
     64 	uint rf_off    = rf_raw_dim.x * rf_channel + time_sample;
     65 
     66 	/* NOTE: rf_data index and stride considering the data is i16 not i32 */
     67 	uint ridx       = rf_off    / RF_SAMPLES_PER_INDEX;
     68 	uint ridx_delta = rf_stride / RF_SAMPLES_PER_INDEX;
     69 
     70 	/* NOTE: rf_data is i16 so each access grabs two time samples at time.
     71 	 * We need to shift arithmetically (maintaining the sign) to get the
     72 	 * desired element. If the time sample is even we take the upper half
     73 	 * and if its odd we take the lower half. */
     74 	uint lfs = ((~time_sample) & 1u) * 16;
     75 
     76 	/* NOTE: Compute N-D dot product */
     77 	vec2 result = vec2(0);
     78 	switch (decode) {
     79 	case DECODE_MODE_NONE: {
     80 		result = RESULT_TYPE_CAST(sample_rf_data(ridx + ridx_delta * acq, lfs));
     81 	} break;
     82 	case DECODE_MODE_HADAMARD: {
     83 		INPUT_DATA_TYPE sum = INPUT_DATA_TYPE(0);
     84 		for (int i = 0; i < dec_data_dim.z; i++) {
     85 			sum  += imageLoad(hadamard, ivec2(i, acq)).x * sample_rf_data(ridx, lfs);
     86 			ridx += ridx_delta;
     87 		}
     88 		result = RESULT_TYPE_CAST(sum) / float(dec_data_dim.z);
     89 	} break;
     90 	}
     91 	out_data[out_off] = result;
     92 }