ogl_beamforming

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

decode.glsl (3061B)


      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 layout(r16i, binding = 1) readonly restrict uniform iimage1D channel_mapping;
     28 
     29 INPUT_DATA_TYPE sample_rf_data(int index, uint lfs)
     30 {
     31 	INPUT_DATA_TYPE result;
     32 #if   defined(INPUT_DATA_TYPE_FLOAT)
     33 	result = rf_data[index];
     34 #elif defined(INPUT_DATA_TYPE_FLOAT_COMPLEX)
     35 	result = rf_data[index];
     36 #else
     37 	result = (rf_data[index] << lfs) >> 16;
     38 #endif
     39 	return result;
     40 }
     41 
     42 void main()
     43 {
     44 	/* NOTE: each invocation takes a time sample and a receive channel.
     45 	 * It first maps the column to the correct column in the rf data then
     46 	 * does the dot product with the equivalent row of the hadamard matrix.
     47 	 * The result is stored to the equivalent row, column index of the output.
     48 	 */
     49 	int time_sample = int(gl_GlobalInvocationID.x);
     50 	int channel     = int(gl_GlobalInvocationID.y);
     51 	int acq         = int(gl_GlobalInvocationID.z);
     52 
     53 	/* NOTE: offsets for storing the results in the output data */
     54 	uint out_off = dec_data_dim.x * dec_data_dim.y * acq + dec_data_dim.x * channel + time_sample;
     55 
     56 	int rf_channel = imageLoad(channel_mapping, channel).x;
     57 
     58 	/* NOTE: stride is the number of samples between acquistions; off is the
     59 	 * index of the first acquisition for this channel and time sample  */
     60 	int rf_stride = int(dec_data_dim.x);
     61 	int rf_off    = int(rf_raw_dim.x) * rf_channel + time_sample;
     62 
     63 	/* NOTE: rf_data index and stride considering the data is i16 not i32 */
     64 	int ridx       = rf_off    / RF_SAMPLES_PER_INDEX;
     65 	int ridx_delta = rf_stride / RF_SAMPLES_PER_INDEX;
     66 
     67 	/* NOTE: rf_data is i16 so each access grabs two time samples at time.
     68 	 * We need to shift arithmetically (maintaining the sign) to get the
     69 	 * desired element. If the time sample is even we take the upper half
     70 	 * and if its odd we take the lower half. */
     71 	uint lfs = ((~time_sample) & 1u) * 16;
     72 
     73 	/* NOTE: Compute N-D dot product */
     74 	vec2 result = vec2(0);
     75 	switch (decode) {
     76 	case DECODE_MODE_NONE: {
     77 		result = RESULT_TYPE_CAST(sample_rf_data(ridx + ridx_delta * acq, lfs));
     78 	} break;
     79 	case DECODE_MODE_HADAMARD: {
     80 		INPUT_DATA_TYPE sum = INPUT_DATA_TYPE(0);
     81 		for (int i = 0; i < dec_data_dim.z; i++) {
     82 			sum  += imageLoad(hadamard, ivec2(i, acq)).x * sample_rf_data(ridx, lfs);
     83 			ridx += ridx_delta;
     84 		}
     85 		result = RESULT_TYPE_CAST(sum) / float(dec_data_dim.z);
     86 	} break;
     87 	}
     88 	out_data[out_off] = result;
     89 }