ogl_beamforming

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

filter.glsl (3003B)


      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 #if ComplexFilter
     13 	#define FILTER_TYPE         vec2
     14 	#define apply_filter(iq, h) complex_mul((iq), (h))
     15 #else
     16 	#define FILTER_TYPE         float
     17 	#define apply_filter(iq, h) ((iq) * (h))
     18 #endif
     19 
     20 layout(std430, binding = 1) readonly restrict buffer buffer_1 {
     21 	DATA_TYPE in_data[];
     22 };
     23 
     24 layout(std430, binding = 2) writeonly restrict buffer buffer_2 {
     25 	DATA_TYPE out_data[];
     26 };
     27 
     28 layout(std430, binding = 3) readonly restrict buffer buffer_3 {
     29 	FILTER_TYPE filter_coefficients[];
     30 };
     31 
     32 layout(r16i, binding = 1) readonly restrict uniform iimage1D channel_mapping;
     33 
     34 vec2 complex_mul(vec2 a, vec2 b)
     35 {
     36 	mat2 m = mat2(b.x, b.y, -b.y, b.x);
     37 	vec2 result = m * a;
     38 	return result;
     39 }
     40 
     41 #if Demodulate
     42 vec2 rotate_iq(vec2 iq, int index)
     43 {
     44 	vec2 result;
     45 	switch (SamplingMode) {
     46 	case SamplingMode_4X:{
     47 		// fs = 2 * fd
     48 		// arg = PI * index
     49 		// cos -> 1 -1  1 -1
     50 		// sin -> 0  0  0  0
     51 		const float scales[2] = {1, -1};
     52 		result = scales[index & 1] * iq;
     53 	}break;
     54 	case SamplingMode_2X:{
     55 		// fs  = fd
     56 		// arg = 2 * PI * index
     57 		// cos -> 1 1 1 1
     58 		// sin -> 0 0 0 0
     59 		result = iq;
     60 	}break;
     61 	default:{
     62 		float arg    = radians(360) * DemodulationFrequency * index / SamplingFrequency;
     63 		mat2  phasor = mat2(cos(arg), -sin(arg),
     64 		                    sin(arg),  cos(arg));
     65 		result = phasor * iq;
     66 	}break;
     67 	}
     68 	return result;
     69 }
     70 #endif
     71 
     72 vec2 sample_rf(uint index)
     73 {
     74 	vec2 result = SAMPLE_TYPE_CAST(in_data[index]);
     75 	return result;
     76 }
     77 
     78 void main()
     79 {
     80 	uint in_sample  = gl_GlobalInvocationID.x * DecimationRate;
     81 	uint out_sample = gl_GlobalInvocationID.x;
     82 	uint channel    = gl_GlobalInvocationID.y;
     83 	uint transmit   = gl_GlobalInvocationID.z;
     84 
     85 	uint in_channel = bool(MapChannels) ? imageLoad(channel_mapping, int(channel)).x : channel;
     86 	uint in_offset  = InputChannelStride * in_channel + InputTransmitStride * transmit;
     87 	uint out_offset = OutputChannelStride  * channel +
     88 	                  OutputTransmitStride * transmit +
     89 	                  OutputSampleStride   * out_sample;
     90 
     91 	int target;
     92 	if (bool(MapChannels)) {
     93 		target = OutputChannelStride / OutputSampleStride;
     94 	} else {
     95 		target = OutputTransmitStride;
     96 	}
     97 
     98 	if (out_sample < target) {
     99 		target *= DecimationRate;
    100 
    101 		vec2 result  = vec2(0);
    102 		int a_length = target;
    103 		int index    = int(in_sample);
    104 
    105 		const float scale = bool(ComplexFilter) ? 1 : sqrt(2);
    106 
    107 		for (int j = max(0, index - FilterLength); j < min(index, a_length); j++) {
    108 			vec2        iq = sample_rf(in_offset + j);
    109 			FILTER_TYPE h  = filter_coefficients[index - j];
    110 		#if 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 }