Commit: e50628ac71c3ab8e3cb69e4070429101b45f0e47
Parent: 0afdb85ac35a8c54946ce288f9ea7e76ae0a7e29
Author: Randy Palamar
Date: Fri, 12 Jul 2024 06:32:25 -0600
first pass at envelope detection
Diffstat:
8 files changed, 205 insertions(+), 68 deletions(-)
diff --git a/beamformer.c b/beamformer.c
@@ -2,33 +2,13 @@
#include "beamformer.h"
static void
-alloc_shader_storage(BeamformerCtx *ctx, Arena a)
+alloc_output_image(BeamformerCtx *ctx)
{
BeamformerParameters *bp = &ctx->params->raw;
- uv4 rf_data_dim = bp->rf_data_dim;
- ctx->csctx.rf_data_dim = rf_data_dim;
- size rf_raw_size = bp->channel_data_stride * rf_data_dim.y * rf_data_dim.z * sizeof(i16);
- size rf_decoded_size = rf_data_dim.x * rf_data_dim.y * rf_data_dim.z * sizeof(f32);
+ ctx->out_data_dim.x = ORONE(bp->output_points.x);
+ ctx->out_data_dim.y = ORONE(bp->output_points.y);
+ ctx->out_data_dim.z = ORONE(bp->output_points.z);
- glDeleteBuffers(ARRAY_COUNT(ctx->csctx.rf_data_ssbos), ctx->csctx.rf_data_ssbos);
- glGenBuffers(ARRAY_COUNT(ctx->csctx.rf_data_ssbos), ctx->csctx.rf_data_ssbos);
- glBindBuffer(GL_SHADER_STORAGE_BUFFER, ctx->csctx.rf_data_ssbos[0]);
- glBufferStorage(GL_SHADER_STORAGE_BUFFER, rf_raw_size, 0,
- GL_DYNAMIC_STORAGE_BIT|GL_MAP_WRITE_BIT);
- glBindBuffer(GL_SHADER_STORAGE_BUFFER, ctx->csctx.rf_data_ssbos[1]);
- glBufferStorage(GL_SHADER_STORAGE_BUFFER, rf_decoded_size, 0, GL_DYNAMIC_STORAGE_BIT);
-
- /* NOTE: store hadamard in GPU once; it won't change for a particular imaging session */
- ctx->csctx.hadamard_dim = (uv2){ .x = rf_data_dim.z, .y = rf_data_dim.z };
- size hadamard_elements = ctx->csctx.hadamard_dim.x * ctx->csctx.hadamard_dim.y;
- i32 *hadamard = alloc(&a, i32, hadamard_elements);
- fill_hadamard(hadamard, ctx->csctx.hadamard_dim.x);
-
- rlUnloadShaderBuffer(ctx->csctx.hadamard_ssbo);
- ctx->csctx.hadamard_ssbo = rlLoadShaderBuffer(hadamard_elements * sizeof(i32), hadamard,
- GL_STATIC_DRAW);
-
- ctx->out_data_dim = bp->output_points;
/* NOTE: allocate storage for beamformed output data;
* this is shared between compute and fragment shaders */
uv4 odim = ctx->out_data_dim;
@@ -44,6 +24,59 @@ alloc_shader_storage(BeamformerCtx *ctx, Arena a)
UnloadRenderTexture(ctx->fsctx.output);
ctx->fsctx.output = LoadRenderTexture(odim.x, odim.y);
+
+ ctx->flags &= ~ALLOC_OUT_TEX;
+}
+
+static void
+upload_filter_coefficients(BeamformerCtx *ctx, Arena a)
+{
+ f32 lpf_coeff[] = {
+ 0.001504252781, 0.006636276841, 0.01834679954, 0.0386288017,
+ 0.06680636108, 0.09852545708, 0.1264867932, 0.1429549307,
+ 0.1429549307, 0.1264867932, 0.09852545708, 0.06680636108,
+ 0.0386288017, 0.01834679954, 0.006636276841, 0.001504252781,
+ };
+ u32 lpf_coeff_count = ARRAY_COUNT(lpf_coeff);
+ ctx->csctx.lpf_order = lpf_coeff_count - 1;
+ rlUnloadShaderBuffer(ctx->csctx.lpf_ssbo);
+ ctx->csctx.lpf_ssbo = rlLoadShaderBuffer(lpf_coeff_count * sizeof(f32), lpf_coeff, GL_STATIC_DRAW);
+
+ ctx->flags &= ~UPLOAD_FILTER;
+}
+
+static void
+alloc_shader_storage(BeamformerCtx *ctx, Arena a)
+{
+ BeamformerParameters *bp = &ctx->params->raw;
+ uv4 rf_data_dim = bp->rf_data_dim;
+ size rf_raw_size = bp->channel_data_stride * rf_data_dim.y * rf_data_dim.z * sizeof(i16);
+ size rf_decoded_size = rf_data_dim.x * rf_data_dim.y * rf_data_dim.z * sizeof(f32) * 2;
+ ctx->csctx.rf_data_dim = rf_data_dim;
+
+ glDeleteBuffers(ARRAY_COUNT(ctx->csctx.rf_data_ssbos), ctx->csctx.rf_data_ssbos);
+ glDeleteBuffers(1, &ctx->csctx.raw_data_ssbo);
+ glGenBuffers(1, &ctx->csctx.raw_data_ssbo);
+ glGenBuffers(ARRAY_COUNT(ctx->csctx.rf_data_ssbos), ctx->csctx.rf_data_ssbos);
+ glBindBuffer(GL_SHADER_STORAGE_BUFFER, ctx->csctx.raw_data_ssbo);
+ glBufferStorage(GL_SHADER_STORAGE_BUFFER, rf_raw_size, 0,
+ GL_DYNAMIC_STORAGE_BIT|GL_MAP_WRITE_BIT);
+
+ for (u32 i = 0; i < ARRAY_COUNT(ctx->csctx.rf_data_ssbos); i++) {
+ glBindBuffer(GL_SHADER_STORAGE_BUFFER, ctx->csctx.rf_data_ssbos[i]);
+ glBufferStorage(GL_SHADER_STORAGE_BUFFER, rf_decoded_size, 0, GL_DYNAMIC_STORAGE_BIT);
+ }
+
+ /* NOTE: store hadamard in GPU once; it won't change for a particular imaging session */
+ ctx->csctx.hadamard_dim = (uv2){.x = rf_data_dim.z, .y = rf_data_dim.z};
+ size hadamard_elements = rf_data_dim.z * rf_data_dim.z;
+ i32 *hadamard = alloc(&a, i32, hadamard_elements);
+ fill_hadamard(hadamard, rf_data_dim.z);
+
+ rlUnloadShaderBuffer(ctx->csctx.hadamard_ssbo);
+ ctx->csctx.hadamard_ssbo = rlLoadShaderBuffer(hadamard_elements * sizeof(i32), hadamard,
+ GL_STATIC_DRAW);
+ ctx->flags &= ~ALLOC_SSBOS;
}
static void
@@ -52,22 +85,34 @@ do_compute_shader(BeamformerCtx *ctx, enum compute_shaders shader)
ComputeShaderCtx *csctx = &ctx->csctx;
glUseProgram(csctx->programs[shader]);
- glBindImageTexture(ctx->out_texture_unit, ctx->out_texture, 0, GL_FALSE, 0,
- GL_WRITE_ONLY, GL_RG32F);
- glUniform1i(csctx->out_data_tex_id, ctx->out_texture_unit);
-
glBindBufferBase(GL_UNIFORM_BUFFER, 0, csctx->shared_ubo);
- u32 rf_ssbo_idx = 0;
- u32 decoded_ssbo_idx = 1;
+ u32 output_ssbo_idx = !csctx->last_active_ssbo_index;
+ u32 input_ssbo_idx = csctx->last_active_ssbo_index;
switch (shader) {
case CS_HADAMARD:
- glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, csctx->rf_data_ssbos[rf_ssbo_idx]);
- glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 2, csctx->rf_data_ssbos[decoded_ssbo_idx]);
+ glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, csctx->raw_data_ssbo);
+ glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 2, csctx->rf_data_ssbos[output_ssbo_idx]);
glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 3, csctx->hadamard_ssbo);
- glDispatchCompute(csctx->rf_data_dim.x / 32, csctx->rf_data_dim.y / 32, csctx->rf_data_dim.z);
+ glDispatchCompute(ORONE(csctx->rf_data_dim.x / 32),
+ ORONE(csctx->rf_data_dim.y / 32),
+ ORONE(csctx->rf_data_dim.z));
+ csctx->last_active_ssbo_index = !csctx->last_active_ssbo_index;
+ break;
+ case CS_LPF:
+ glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, csctx->rf_data_ssbos[input_ssbo_idx]);
+ glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 2, csctx->rf_data_ssbos[output_ssbo_idx]);
+ glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 3, csctx->lpf_ssbo);
+ glUniform1i(csctx->lpf_order_id, csctx->lpf_order);
+ glDispatchCompute(ORONE(csctx->rf_data_dim.x / 32),
+ ORONE(csctx->rf_data_dim.y / 32),
+ ORONE(csctx->rf_data_dim.z));
+ csctx->last_active_ssbo_index = !csctx->last_active_ssbo_index;
break;
case CS_MIN_MAX:
+ glBindImageTexture(ctx->out_texture_unit, ctx->out_texture, 0, GL_FALSE, 0,
+ GL_WRITE_ONLY, GL_RG32F);
+ glUniform1i(csctx->out_data_tex_id, ctx->out_texture_unit);
for (u32 i = 1; i < ctx->out_texture_mips; i++) {
u32 otu = ctx->out_texture_unit;
glBindImageTexture(otu + 1, ctx->out_texture, i - 1,
@@ -78,7 +123,6 @@ do_compute_shader(BeamformerCtx *ctx, enum compute_shaders shader)
glUniform1i(csctx->mip_view_tex_id, otu + 2);
glUniform1i(csctx->mips_level_id, i);
- #define ORONE(x) ((x)? (x) : 1)
u32 width = ctx->out_data_dim.x >> i;
u32 height = ctx->out_data_dim.y >> i;
u32 depth = ctx->out_data_dim.z >> i;
@@ -87,8 +131,15 @@ do_compute_shader(BeamformerCtx *ctx, enum compute_shaders shader)
}
break;
case CS_UFORCES:
- glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, csctx->rf_data_ssbos[decoded_ssbo_idx]);
- glDispatchCompute(ctx->out_data_dim.x / 32, ctx->out_data_dim.y / 32, ctx->out_data_dim.z);
+ glActiveTexture(GL_TEXTURE0 + ctx->out_texture_unit);
+ glBindTexture(GL_TEXTURE_3D, ctx->out_texture);
+ glBindImageTexture(ctx->out_texture_unit, ctx->out_texture, 0, GL_FALSE, 0,
+ GL_WRITE_ONLY, GL_RG32F);
+ glUniform1i(csctx->out_data_tex_id, ctx->out_texture_unit);
+ glBindBufferBase(GL_SHADER_STORAGE_BUFFER, 1, csctx->rf_data_ssbos[input_ssbo_idx]);
+ glDispatchCompute(ORONE(ctx->out_data_dim.x / 32),
+ ORONE(ctx->out_data_dim.y / 32),
+ ORONE(ctx->out_data_dim.z));
break;
default: ASSERT(0);
}
@@ -149,17 +200,19 @@ draw_settings_ui(BeamformerCtx *ctx, Arena arena, f32 dt, Rect r, v2 mouse)
f32 data_scale;
b32 editable;
} listings[] = {
- { "Sampling Rate:", " [MHz]", &bp->sampling_frequency, 1e-6, 0 },
- { "Speed of Sound:", " [m/s]", &bp->speed_of_sound, 1, 1 },
- { "Min X Point:", " [mm]", &bp->output_min_xz.x, 1e3, 1 },
- { "Max X Point:", " [mm]", &bp->output_max_xz.x, 1e3, 1 },
- { "Min Z Point:", " [mm]", &bp->output_min_xz.y, 1e3, 1 },
- { "Max Z Point:", " [mm]", &bp->output_max_xz.y, 1e3, 1 },
- { "Dynamic Range:", " [dB]", &ctx->fsctx.db, 1, 1 },
+ { "Sampling Rate:", " [MHz]", &bp->sampling_frequency, 1e-6, 0 },
+ { "Center Frequency:", " [MHz]", &bp->center_frequency, 1e-6, 1 },
+ { "Speed of Sound:", " [m/s]", &bp->speed_of_sound, 1, 1 },
+ { "Min X Point:", " [mm]", &bp->output_min_xz.x, 1e3, 1 },
+ { "Max X Point:", " [mm]", &bp->output_max_xz.x, 1e3, 1 },
+ { "Min Z Point:", " [mm]", &bp->output_min_xz.y, 1e3, 1 },
+ { "Max Z Point:", " [mm]", &bp->output_max_xz.y, 1e3, 1 },
+ { "Dynamic Range:", " [dB]", &ctx->fsctx.db, 1, 1 },
};
struct { f32 min, max; } limits[] = {
{0},
+ {0, 100e6},
{0, 1e6},
{-1e3, bp->output_max_xz.x - 1e-6},
{bp->output_min_xz.x + 1e-6, 1e3},
@@ -351,12 +404,12 @@ do_beamformer(BeamformerCtx *ctx, Arena arena)
BeamformerParameters *bp = &ctx->params->raw;
/* NOTE: Check for and Load RF Data into GPU */
if (os_poll_pipe(ctx->data_pipe)) {
- if (!uv4_equal(ctx->csctx.rf_data_dim, bp->rf_data_dim) ||
- !uv4_equal(ctx->out_data_dim, bp->output_points)) {
+ if (!uv4_equal(ctx->csctx.rf_data_dim, bp->rf_data_dim) || ctx->flags & ALLOC_SSBOS)
alloc_shader_storage(ctx, arena);
- }
+ if (!uv4_equal(ctx->out_data_dim, bp->output_points) || ctx->flags & ALLOC_OUT_TEX)
+ alloc_output_image(ctx);
- glBindBuffer(GL_SHADER_STORAGE_BUFFER, ctx->csctx.rf_data_ssbos[0]);
+ glBindBuffer(GL_SHADER_STORAGE_BUFFER, ctx->csctx.raw_data_ssbo);
void *rf_data_buf = glMapBuffer(GL_SHADER_STORAGE_BUFFER, GL_WRITE_ONLY);
ASSERT(rf_data_buf);
uv4 rf_data_dim = ctx->csctx.rf_data_dim;
@@ -368,6 +421,9 @@ do_beamformer(BeamformerCtx *ctx, Arena arena)
else ctx->partial_transfer_count++;
}
+ if (ctx->flags & UPLOAD_FILTER)
+ upload_filter_coefficients(ctx, arena);
+
if (ctx->flags & DO_COMPUTE) {
if (ctx->params->upload) {
glBindBuffer(GL_UNIFORM_BUFFER, ctx->csctx.shared_ubo);
@@ -378,6 +434,7 @@ do_beamformer(BeamformerCtx *ctx, Arena arena)
ctx->params->upload = 0;
}
do_compute_shader(ctx, CS_HADAMARD);
+ do_compute_shader(ctx, CS_LPF);
do_compute_shader(ctx, CS_UFORCES);
do_compute_shader(ctx, CS_MIN_MAX);
ctx->flags &= ~DO_COMPUTE;
diff --git a/beamformer.h b/beamformer.h
@@ -33,6 +33,7 @@ enum compute_shaders {
// CS_FORCES,
CS_HADAMARD,
// CS_HERCULES,
+ CS_LPF,
CS_MIN_MAX,
CS_UFORCES,
CS_LAST
@@ -40,14 +41,19 @@ enum compute_shaders {
enum program_flags {
RELOAD_SHADERS = 1 << 0,
- DO_COMPUTE = 1 << 1,
+ ALLOC_SSBOS = 1 << 1,
+ ALLOC_OUT_TEX = 1 << 2,
+ UPLOAD_FILTER = 1 << 3,
+ DO_COMPUTE = 1 << 30,
};
typedef struct {
u32 programs[CS_LAST];
- /* NOTE: One SSBO for raw data and one for decoded data */
+ /* NOTE: One SSBO for raw data and two for decoded data (swapped for chained stages)*/
+ u32 raw_data_ssbo;
u32 rf_data_ssbos[2];
+ u32 last_active_ssbo_index;
u32 hadamard_ssbo;
uv2 hadamard_dim;
@@ -57,6 +63,10 @@ typedef struct {
i32 out_data_tex_id;
i32 mip_view_tex_id;
i32 mips_level_id;
+
+ u32 lpf_ssbo;
+ u32 lpf_order;
+ i32 lpf_order_id;
} ComputeShaderCtx;
typedef struct {
diff --git a/beamformer_parameters.h b/beamformer_parameters.h
@@ -15,6 +15,7 @@ typedef struct {
u32 channel_offset; /* Offset into channel_mapping: 0 or 128 (rows or columns) */
f32 speed_of_sound; /* [m/s] */
f32 sampling_frequency; /* [Hz] */
+ f32 center_frequency; /* [Hz] */
f32 focal_depth; /* [m] */
f32 _pad[3];
} BeamformerParameters;
diff --git a/main.c b/main.c
@@ -2,8 +2,9 @@
#include "beamformer.h"
static char *compute_shader_paths[CS_LAST] = {
- [CS_MIN_MAX] = "shaders/min_max.glsl",
[CS_HADAMARD] = "shaders/hadamard.glsl",
+ [CS_LPF] = "shaders/lpf.glsl",
+ [CS_MIN_MAX] = "shaders/min_max.glsl",
[CS_UFORCES] = "shaders/uforces.glsl",
};
@@ -122,6 +123,7 @@ reload_shaders(BeamformerCtx *ctx, Arena a)
glDeleteShader(shader_id);
}
+ csctx->lpf_order_id = glGetUniformLocation(csctx->programs[CS_LPF], "u_lpf_order");
csctx->out_data_tex_id = glGetUniformLocation(csctx->programs[CS_UFORCES], "u_out_data_tex");
csctx->mip_view_tex_id = glGetUniformLocation(csctx->programs[CS_MIN_MAX], "u_mip_view_tex");
csctx->mips_level_id = glGetUniformLocation(csctx->programs[CS_MIN_MAX], "u_mip_map");
@@ -172,7 +174,7 @@ main(void)
glBindBuffer(GL_UNIFORM_BUFFER, ctx.csctx.shared_ubo);
glBufferData(GL_UNIFORM_BUFFER, sizeof(BeamformerParameters), 0, GL_STATIC_DRAW);
- ctx.flags |= RELOAD_SHADERS;
+ ctx.flags |= RELOAD_SHADERS|ALLOC_SSBOS|ALLOC_OUT_TEX|UPLOAD_FILTER;
while(!WindowShouldClose()) {
do_debug();
diff --git a/shaders/hadamard.glsl b/shaders/hadamard.glsl
@@ -7,7 +7,7 @@ layout(std430, binding = 1) readonly restrict buffer buffer_1 {
};
layout(std430, binding = 2) writeonly restrict buffer buffer_2 {
- float out_data[];
+ vec2 out_data[];
};
layout(std430, binding = 3) readonly restrict buffer buffer_3 {
@@ -27,6 +27,7 @@ layout(std140, binding = 0) uniform parameters {
uint channel_offset; /* Offset into channel_mapping: 0 or 128 (rows or columns) */
float speed_of_sound; /* [m/s] */
float sampling_frequency; /* [Hz] */
+ float center_frequency; /* [Hz] */
float focal_depth; /* [m] */
};
@@ -74,5 +75,7 @@ void main()
ridx += ridx_delta;
}
- out_data[out_off + out_stride * acq] = float(sum);
+ float arg = radians(360) * center_frequency * time_sample / sampling_frequency;
+ out_data[out_off + out_stride * acq].x = float(sum) * cos(arg);
+ out_data[out_off + out_stride * acq].y = float(sum) * sin(arg);
}
diff --git a/shaders/lpf.glsl b/shaders/lpf.glsl
@@ -0,0 +1,55 @@
+/* See LICENSE for license details. */
+#version 460 core
+layout(local_size_x = 32, local_size_y = 32, local_size_z = 1) in;
+
+layout(std430, binding = 1) readonly restrict buffer buffer_1 {
+ vec2 in_data[];
+};
+
+layout(std430, binding = 2) writeonly restrict buffer buffer_2 {
+ vec2 out_data[];
+};
+
+layout(std430, binding = 3) readonly restrict buffer buffer_3 {
+ float lpf_coeff[];
+};
+
+layout(std140, binding = 0) uniform parameters {
+ uvec4 channel_mapping[64]; /* Transducer Channel to Verasonics Channel */
+ uvec4 uforces_channels[32]; /* Channels used for virtual UFORCES elements */
+ uvec4 rf_data_dim; /* Samples * Channels * Acquisitions; last element ignored */
+ uvec4 output_points; /* Width * Height * Depth; last element ignored */
+ vec2 output_min_xz; /* [m] Top left corner of output region */
+ vec2 output_max_xz; /* [m] Bottom right corner of output region */
+ vec2 xdc_min_xy; /* [m] Min center of transducer elements */
+ vec2 xdc_max_xy; /* [m] Max center of transducer elements */
+ uint channel_data_stride; /* Data points between channels (samples * acq + padding) */
+ uint channel_offset; /* Offset into channel_mapping: 0 or 128 (rows or columns) */
+ float speed_of_sound; /* [m/s] */
+ float sampling_frequency; /* [Hz] */
+ float center_frequency; /* [Hz] */
+ float focal_depth; /* [m] */
+};
+
+layout(location = 1) uniform uint u_lpf_order;
+
+void main()
+{
+ uint time_sample = gl_GlobalInvocationID.x;
+ uint channel = gl_GlobalInvocationID.y;
+ uint acq = gl_GlobalInvocationID.z;
+
+ /* NOTE: offsets for storing the results in the output data */
+ uint stride = rf_data_dim.x * rf_data_dim.y;
+ uint off = rf_data_dim.x * channel + time_sample;
+
+ vec2 sum = vec2(0);
+ for (int i = 0; i <= u_lpf_order; i++) {
+ vec2 data;
+ /* NOTE: make sure data samples come from the same acquisition */
+ if (time_sample > i) data = in_data[off - i];
+ else data = vec2(0);
+ sum += lpf_coeff[i] * data;
+ }
+ out_data[stride * acq + off] = sum;
+}
diff --git a/shaders/uforces.glsl b/shaders/uforces.glsl
@@ -3,7 +3,7 @@
layout(local_size_x = 32, local_size_y = 32, local_size_z = 1) in;
layout(std430, binding = 1) readonly restrict buffer buffer_1 {
- float rf_data[];
+ vec2 rf_data[];
};
layout(std140, binding = 0) uniform parameters {
@@ -19,16 +19,16 @@ layout(std140, binding = 0) uniform parameters {
uint channel_offset; /* Offset into channel_mapping: 0 or 128 (rows or columns) */
float speed_of_sound; /* [m/s] */
float sampling_frequency; /* [Hz] */
+ float center_frequency; /* [Hz] */
float focal_depth; /* [m] */
};
-//layout(location = 6) uniform sampler2D u_element_positions;
layout(rg32f, location = 1) uniform image3D u_out_data_tex;
#define C_SPLINE 0.5
/* NOTE: See: https://cubic.org/docs/hermite.htm */
-float cubic(uint ridx, float x)
+vec2 cubic(uint ridx, float x)
{
mat4 h = mat4(
2, -3, 0, 1,
@@ -41,12 +41,14 @@ float cubic(uint ridx, float x)
float t = (x - float(xk));
vec4 S = vec4(t * t * t, t * t, t, 1);
- float P1 = rf_data[ridx + xk];
- float P2 = rf_data[ridx + xk + 1];
- float T1 = C_SPLINE * (P2 - rf_data[ridx + xk - 1]);
- float T2 = C_SPLINE * (rf_data[ridx + xk + 2] - P1);
- vec4 C = vec4(P1, P2, T1, T2);
- return dot(S, h * C);
+ vec2 P1 = rf_data[ridx + xk];
+ vec2 P2 = rf_data[ridx + xk + 1];
+ vec2 T1 = C_SPLINE * (P2 - rf_data[ridx + xk - 1]);
+ vec2 T2 = C_SPLINE * (rf_data[ridx + xk + 2] - P1);
+
+ vec4 C1 = vec4(P1.x, P2.x, T1.x, T2.x);
+ vec4 C2 = vec4(P1.y, P2.y, T1.y, T2.y);
+ return vec2(dot(S, h * C1), dot(S, h * C2));
}
void main()
@@ -71,7 +73,9 @@ void main()
float dx = xdc_size.x / float(rf_data_dim.y);
float dzsign = sign(image_point.z - focal_depth);
- float sum = 0;
+ vec2 sum = vec2(0);
+
+ float time_scale = radians(360) * center_frequency;
uint ridx = rf_data_dim.y * rf_data_dim.x;
/* NOTE: skip first acquisition since its garbage */
@@ -84,15 +88,19 @@ void main()
vec2 rdist = vec2(x, image_point.z);
for (uint j = 0; j < rf_data_dim.y; j++) {
- float dist = transmit_dist + length(rdist);
- float rsample = dist * sampling_frequency / speed_of_sound;
+ float dist = transmit_dist + length(rdist);
+ float time = dist / speed_of_sound;
/* NOTE: do cubic interp between adjacent time samples */
- sum += cubic(ridx, rsample);
+ vec2 p = cubic(ridx, time * sampling_frequency);
+ p.x *= cos(time_scale * time);
+ p.y *= sin(time_scale * time);
+ sum += p;
rdist.x -= dx;
ridx += rf_data_dim.x;
}
ridx += rf_data_dim.y * rf_data_dim.x;
}
- imageStore(u_out_data_tex, out_coord, vec4(sum, sum, 0, 0));
+ float val = length(sum);
+ imageStore(u_out_data_tex, out_coord, vec4(val, val, 0, 0));
}
diff --git a/util.h b/util.h
@@ -15,9 +15,10 @@
#define ARRAY_COUNT(a) (sizeof(a) / sizeof(*a))
#define ABS(x) ((x) < 0 ? (-x) : (x))
-#define MAX(a, b) ((a) > (b) ? (a) : (b))
#define CLAMP(x, a, b) ((x) = (x) < (a) ? (a) : (x) > (b) ? (b) : (x))
#define ISPOWEROF2(a) (((a) & ((a) - 1)) == 0)
+#define MAX(a, b) ((a) > (b) ? (a) : (b))
+#define ORONE(x) ((x)? (x) : 1)
#define MEGABYTE (1024ULL * 1024ULL)
#define GIGABYTE (1024ULL * 1024ULL * 1024ULL)