ogl_beamforming

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

filter.glsl (4104B)


      1 /* See LICENSE for license details. */
      2 #if  (InputDataKind == DataKind_Int16Complex         || \
      3      (InputDataKind == DataKind_Int16 && Demodulate) || \
      4      (InputDataKind == DataKind_Float16 && Demodulate))
      5   #define SAMPLE_TYPE f16vec2
      6 #elif InputDataKind == DataKind_Int16
      7   #define SAMPLE_TYPE f16
      8 #elif InputDataKind == DataKind_Float32 && Demodulate
      9   #define SAMPLE_TYPE f32vec2
     10 #endif
     11 
     12 #ifndef SAMPLE_TYPE
     13   #define SAMPLE_TYPE InputDataType
     14 #endif
     15 
     16 #define ComplexSampleType (InputDataKind == DataKind_Float32Complex || \
     17                            InputDataKind == DataKind_Float16Complex || \
     18                            InputDataKind == DataKind_Int16Complex   || \
     19                            Demodulate)
     20 #if ComplexSampleType
     21   #define RESULT_TYPE f32vec2
     22 #else
     23   #define RESULT_TYPE f32
     24 #endif
     25 
     26 #if ComplexFilter
     27   #define FILTER_TYPE f32vec2
     28 #else
     29   #define FILTER_TYPE f32
     30 #endif
     31 
     32 #if ComplexFilter && ComplexSampleType
     33   #define apply_filter(iq, h) complex_mul(f32vec2(iq), f32vec2(h))
     34 #else
     35   #define apply_filter(iq, h) ((iq) * (h))
     36 #endif
     37 
     38 layout(std430, buffer_reference) restrict readonly  buffer Input  { InputDataType x[];  };
     39 layout(std430, buffer_reference) restrict writeonly buffer Output { OutputDataType x[]; };
     40 
     41 layout(std430, buffer_reference, buffer_reference_align = 64) restrict readonly buffer Filter {
     42 	FILTER_TYPE values[FilterLength];
     43 };
     44 
     45 f32vec2 complex_mul(f32vec2 a, f32vec2 b)
     46 {
     47 	mat2 m = mat2(b.x, b.y, -b.y, b.x);
     48 	f32vec2 result = m * a;
     49 	return result;
     50 }
     51 
     52 #if Demodulate
     53 SAMPLE_TYPE rotate_iq(SAMPLE_TYPE iq, uint index)
     54 {
     55 	float arg          = radians(360) * DemodulationFrequency * index / SamplingFrequency;
     56 	SAMPLE_TYPE result = SAMPLE_TYPE(complex_mul(iq, f32vec2(cos(arg), -sin(arg))));
     57 	return result;
     58 }
     59 #endif
     60 
     61 shared SAMPLE_TYPE rf[DecimationRate * gl_WorkGroupSize.x + FilterLength - 1];
     62 
     63 void main()
     64 {
     65 	uint out_sample = gl_GlobalInvocationID.x;
     66 	uint channel    = gl_GlobalInvocationID.y;
     67 	uint transmit   = gl_GlobalInvocationID.z;
     68 
     69 	/////////////////////////
     70 	// NOTE: sample caching
     71 	{
     72 		bool offset_wraps = (DecimationRate * gl_WorkGroupID.x * gl_WorkGroupSize.x) < (FilterLength - 1);
     73 
     74 		u32 in_offset = InputDataKindByteSize * (InputChannelStride * channel + InputTransmitStride * transmit);
     75 		// NOTE(rnp): when demodulating we want to load 2 elements at a time but the
     76 		// input strides were specified in terms of a single element. therefore we
     77 		// must divide this by two. by doing this here we can gracefully handle
     78 		// the case where there are an odd number of samples (this drops the last one).
     79 		if (Demodulate)
     80 			in_offset /= 2;
     81 
     82 		// NOTE(rnp): broken out to avoid overflow from the subtraction
     83 		u64 input_address = input_buffer + in_offset;
     84 		input_address += InputDataKindByteSize * (DecimationRate * gl_WorkGroupID.x * gl_WorkGroupSize.x);
     85 		input_address -= InputDataKindByteSize * (FilterLength - 1);
     86 
     87 		const SAMPLE_TYPE scale = SAMPLE_TYPE(bool(ComplexFilter) ? 1 : sqrt(2.0f));
     88 		for (u32 index = gl_LocalInvocationIndex; index < rf.length(); index += gl_WorkGroupSize.x) {
     89 			SAMPLE_TYPE s = SAMPLE_TYPE(0);
     90 			if (!offset_wraps || index >= FilterLength - 1) {
     91 				s = SAMPLE_TYPE(Input(input_address).x[index]);
     92 				#if Demodulate
     93 				s = scale * rotate_iq(s * SAMPLE_TYPE(1, -1), index);
     94 				#endif
     95 			}
     96 			rf[index] = s;
     97 		}
     98 	}
     99 	barrier();
    100 
    101 	Filter f = Filter(HeapBase + FilterCoefficients);
    102 	if (out_sample < SampleCount / DecimationRate) {
    103 		RESULT_TYPE result = RESULT_TYPE(0);
    104 		u32 offset = DecimationRate * gl_LocalInvocationIndex;
    105 		for (u32 j = 0; j < FilterLength; j++)
    106 			result += apply_filter(rf[offset + j], f.values[j]);
    107 
    108 		u32 out_offset = OutputChannelStride  * channel +
    109 		                 OutputTransmitStride * transmit +
    110 		                 OutputSampleStride   * out_sample;
    111 
    112 		if (BatchSampleCount != 0) {
    113 			// NOTE(rnp): deinterleave
    114 			Output(output_buffer).x[out_offset] = OutputDataType(result.x);
    115 			out_offset += BatchSampleCount;
    116 			Output(output_buffer).x[out_offset] = OutputDataType(result.y);
    117 		} else {
    118 			Output(output_buffer).x[out_offset] = OutputDataType(result);
    119 		}
    120 	}
    121 }