ogl_beamforming

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

decode.glsl (3026B)


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