ogl_beamforming

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

filter.glsl (3381B)


      1 /* See LICENSE for license details. */
      2 #if   defined(INPUT_DATA_TYPE_FLOAT)
      3   #define DATA_TYPE           vec2
      4   #define RESULT_TYPE_CAST(v) (v)
      5   #define SAMPLE_TYPE_CAST(v) (v)
      6 #else
      7   #define DATA_TYPE           uint
      8   #define RESULT_TYPE_CAST(v) packSnorm2x16(v)
      9   #define SAMPLE_TYPE_CAST(v) unpackSnorm2x16(v)
     10 #endif
     11 
     12 layout(std430, binding = 1) readonly restrict buffer buffer_1 {
     13 	DATA_TYPE in_data[];
     14 };
     15 
     16 layout(std430, binding = 2) writeonly restrict buffer buffer_2 {
     17 	DATA_TYPE out_data[];
     18 };
     19 
     20 layout(r16i, binding = 1) readonly restrict uniform iimage1D channel_mapping;
     21 
     22 #if COMPLEX_FILTER
     23 	layout(rg32f, binding = 0) readonly restrict uniform  image1D filter_coefficients;
     24 	#define apply_filter(iq, h) complex_mul((iq), (h).xy)
     25 #else
     26 	layout(r32f, binding = 0) readonly restrict uniform  image1D filter_coefficients;
     27 	#define apply_filter(iq, h) ((iq) * (h).x)
     28 #endif
     29 
     30 vec2 complex_mul(vec2 a, vec2 b)
     31 {
     32 	mat2 m = mat2(b.x, b.y, -b.y, b.x);
     33 	vec2 result = m * a;
     34 	return result;
     35 }
     36 
     37 vec2 rotate_iq(vec2 iq, int index)
     38 {
     39 	vec2 result;
     40 	/* TODO(rnp): this doesn't give us the same performance boost as hardcoding the mode */
     41 	switch (shader_flags & ShaderFlags_SamplingModeMask) {
     42 	case SamplingMode_NS200BW:{
     43 		// fs = 2 * fd
     44 		// arg = PI * index
     45 		// cos -> 1 -1  1 -1
     46 		// sin -> 0  0  0  0
     47 		/* NOTE(rnp): faster than taking iq or -iq, good job shader compiler */
     48 		if (bool(index & 1)) result = mat2(-1, 0, 0, -1) * iq;
     49 		else                 result = mat2( 1, 0, 0,  1) * iq;
     50 	}break;
     51 	case SamplingMode_BS100BW:{
     52 		// fs  = fd
     53 		// arg = 2 * PI * index
     54 		// cos -> 1 1 1 1
     55 		// sin -> 0 0 0 0
     56 		result = iq;
     57 	}break;
     58 	default:{
     59 		float arg    = radians(360) * demodulation_frequency * index / sampling_frequency;
     60 		mat2  phasor = mat2(cos(arg), -sin(arg),
     61 		                    sin(arg),  cos(arg));
     62 		result = phasor * iq;
     63 	}break;
     64 	}
     65 	return result;
     66 }
     67 
     68 vec2 sample_rf(uint index)
     69 {
     70 	vec2 result = SAMPLE_TYPE_CAST(in_data[index]);
     71 	return result;
     72 }
     73 
     74 void main()
     75 {
     76 	uint in_sample  = gl_GlobalInvocationID.x * decimation_rate;
     77 	uint out_sample = gl_GlobalInvocationID.x;
     78 	uint channel    = gl_GlobalInvocationID.y;
     79 	uint transmit   = gl_GlobalInvocationID.z;
     80 
     81 	bool map_channels = bool(shader_flags & ShaderFlags_MapChannels);
     82 	uint in_channel = map_channels ? imageLoad(channel_mapping, int(channel)).x : channel;
     83 	uint in_offset  = input_channel_stride * in_channel + input_transmit_stride * transmit;
     84 	uint out_offset = output_channel_stride  * channel +
     85 	                  output_transmit_stride * transmit +
     86 	                  output_sample_stride   * out_sample;
     87 
     88 	int target;
     89 	if (map_channels) {
     90 		target = int(output_channel_stride / output_sample_stride);
     91 	} else {
     92 		target = int(output_transmit_stride);
     93 	}
     94 
     95 	if (out_sample < target) {
     96 		target *= int(decimation_rate);
     97 
     98 		vec2 result  = vec2(0);
     99 		int a_length = target;
    100 		int b_length = imageSize(filter_coefficients).x;
    101 		int index    = int(in_sample);
    102 
    103 		const float scale = bool(COMPLEX_FILTER) ? 1 : sqrt(2);
    104 
    105 		for (int j = max(0, index - b_length); j < min(index, a_length); j++) {
    106 			vec2 iq  = sample_rf(in_offset + j);
    107 			vec4 h   = imageLoad(filter_coefficients, index - j);
    108 		#if defined(DEMODULATE)
    109 			result  += scale * apply_filter(rotate_iq(iq * vec2(1, -1), -j), h);
    110 		#else
    111 			result  += apply_filter(iq, h);
    112 		#endif
    113 		}
    114 
    115 		out_data[out_offset] = RESULT_TYPE_CAST(result);
    116 	}
    117 }