ogl_beamforming

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

filter.glsl (3374B)


      1 /* See LICENSE for license details. */
      2 #if DataKind == DataKind_Float32
      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 (ShaderFlags & ShaderFlags_ComplexFilter)
     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 const bool map_channels = (ShaderFlags & ShaderFlags_MapChannels) != 0;
     31 
     32 vec2 complex_mul(vec2 a, vec2 b)
     33 {
     34 	mat2 m = mat2(b.x, b.y, -b.y, b.x);
     35 	vec2 result = m * a;
     36 	return result;
     37 }
     38 
     39 #if (ShaderFlags & ShaderFlags_Demodulate)
     40 vec2 rotate_iq(vec2 iq, int index)
     41 {
     42 	vec2 result;
     43 	switch (SamplingMode) {
     44 	case SamplingMode_4X:{
     45 		// fs = 2 * fd
     46 		// arg = PI * index
     47 		// cos -> 1 -1  1 -1
     48 		// sin -> 0  0  0  0
     49 		/* NOTE(rnp): faster than taking iq or -iq, good job shader compiler */
     50 		if (bool(index & 1)) result = mat2(-1, 0, 0, -1) * iq;
     51 		else                 result = mat2( 1, 0, 0,  1) * iq;
     52 	}break;
     53 	case SamplingMode_2X:{
     54 		// fs  = fd
     55 		// arg = 2 * PI * index
     56 		// cos -> 1 1 1 1
     57 		// sin -> 0 0 0 0
     58 		result = iq;
     59 	}break;
     60 	default:{
     61 		float arg    = radians(360) * demodulation_frequency * index / sampling_frequency;
     62 		mat2  phasor = mat2(cos(arg), -sin(arg),
     63 		                    sin(arg),  cos(arg));
     64 		result = phasor * iq;
     65 	}break;
     66 	}
     67 	return result;
     68 }
     69 #endif
     70 
     71 vec2 sample_rf(uint index)
     72 {
     73 	vec2 result = SAMPLE_TYPE_CAST(in_data[index]);
     74 	return result;
     75 }
     76 
     77 void main()
     78 {
     79 	uint in_sample  = gl_GlobalInvocationID.x * decimation_rate;
     80 	uint out_sample = gl_GlobalInvocationID.x;
     81 	uint channel    = gl_GlobalInvocationID.y;
     82 	uint transmit   = gl_GlobalInvocationID.z;
     83 
     84 	uint in_channel = map_channels ? imageLoad(channel_mapping, int(channel)).x : channel;
     85 	uint in_offset  = input_channel_stride * in_channel + input_transmit_stride * transmit;
     86 	uint out_offset = output_channel_stride  * channel +
     87 	                  output_transmit_stride * transmit +
     88 	                  output_sample_stride   * out_sample;
     89 
     90 	int target;
     91 	if (map_channels) {
     92 		target = int(output_channel_stride / output_sample_stride);
     93 	} else {
     94 		target = int(output_transmit_stride);
     95 	}
     96 
     97 	if (out_sample < target) {
     98 		target *= int(decimation_rate);
     99 
    100 		vec2 result  = vec2(0);
    101 		int a_length = target;
    102 		int b_length = imageSize(filter_coefficients).x;
    103 		int index    = int(in_sample);
    104 
    105 		const float scale = bool(ShaderFlags & ShaderFlags_ComplexFilter) ? 1 : sqrt(2);
    106 
    107 		for (int j = max(0, index - b_length); j < min(index, a_length); j++) {
    108 			vec2 iq  = sample_rf(in_offset + j);
    109 			vec4 h   = imageLoad(filter_coefficients, index - j);
    110 		#if (ShaderFlags & ShaderFlags_Demodulate)
    111 			result  += scale * apply_filter(rotate_iq(iq * vec2(1, -1), -j), h);
    112 		#else
    113 			result  += apply_filter(iq, h);
    114 		#endif
    115 		}
    116 
    117 		out_data[out_offset] = RESULT_TYPE_CAST(result);
    118 	}
    119 }