From d482af4a400ed3fef9b03dfa6b4e985b0ef240ff Mon Sep 17 00:00:00 2001 From: Samuel Le Meur-Diebolt Date: Wed, 7 Oct 2026 18:16:33 +0100 Subject: [PATCH 01/11] fix(matlab): correct BS67BW Fourier frequencies --- .../core/imaging/echoframe_setup_fourier.m | 3 +- .../core/verasonics/get_system_parameters.m | 3 +- echoframe/matlab/tests/verify_core_headless.m | 54 +++++++++++++++++++ 3 files changed, 58 insertions(+), 2 deletions(-) diff --git a/echoframe/matlab/core/imaging/echoframe_setup_fourier.m b/echoframe/matlab/core/imaging/echoframe_setup_fourier.m index f460fd3..03f3551 100644 --- a/echoframe/matlab/core/imaging/echoframe_setup_fourier.m +++ b/echoframe/matlab/core/imaging/echoframe_setup_fourier.m @@ -134,7 +134,8 @@ case 'BS67BW' Nz = single(ReceiveSpec.nSamples * 2); freq_scale = 2; - lower = Nz/2; + % The negative-frequency half starts after the positive half in FFT order. + lower = Nz/2 + 1; freq_window = lower:(lower + Nz/4 - 1); freq_mapping = 1:Nz/4; freq_axis_shift = false; diff --git a/echoframe/matlab/core/verasonics/get_system_parameters.m b/echoframe/matlab/core/verasonics/get_system_parameters.m index 3c1f12d..64566e9 100644 --- a/echoframe/matlab/core/verasonics/get_system_parameters.m +++ b/echoframe/matlab/core/verasonics/get_system_parameters.m @@ -14,7 +14,7 @@ % % nChannels 128 or 256, from Resource.Parameters.numRcvChannels % and Trans.numelements -% Fs Trans.frequency * 4 (Hz) +% Fs pre-quadrature sample rate: 4*Fc, or 4*Fc/3 for BS67BW (Hz) % samplesPerWavelength from samplingMode: 1 (BS50BW), 4/3 (BS67BW), % 2 (BS100BW), 4 (NS200BW) % @@ -52,6 +52,7 @@ case 'BS50BW' ReceiveSpec.samplesPerWavelength = 1; case 'BS67BW' + ReceiveSpec.Fs = Trans.frequency * (4/3) * 1e6; ReceiveSpec.samplesPerWavelength = 4/3; case 'BS100BW' ReceiveSpec.samplesPerWavelength = 2; diff --git a/echoframe/matlab/tests/verify_core_headless.m b/echoframe/matlab/tests/verify_core_headless.m index f2210ce..cbaeaf0 100644 --- a/echoframe/matlab/tests/verify_core_headless.m +++ b/echoframe/matlab/tests/verify_core_headless.m @@ -156,6 +156,8 @@ assert(threw, 'validate_structs accepted a ReceiveSpec missing samplingMode'); fprintf(' missing required field rejected OK\n'); +local_verify_fourier_sampling(); + %% ------------------------------------------------------- 5. cropped frames fprintf('\n=== stored_frame_size (cropped storage geometry) ===\n'); @@ -490,3 +492,55 @@ PDISpec.cropPDI = logical(false); PDISpec.svdMethod = 'Covariance'; end + +function local_verify_fourier_sampling() +% A known transmit delay must be cancelled at the dominant receive-tone bin. +modes = {'BS50BW', 'BS67BW', 'BS100BW', 'NS200BW'}; +rateFactors = [4, 4/3, 4, 4]; +P = struct('pitch', 300e-6, 'nElements', 8); +T = struct('c0', 1540, 'steer', 0, 'apodization', ones(8,1)); +Resource = struct('Parameters', struct('numRcvChannels', 128)); +for fc = [1e6, 31.25e6] + Trans = struct('frequency', fc/1e6, 'numelements', 8); + for n = [128, 512] + rec = struct('c0', 1540, 'nx', 128, 'nz', n/2); + for m = 1:numel(modes) + r = struct('samplingMode', modes{m}, 'nBuffers', 1, ... + 'nSamples', n, 'nSamplesIQ', n/2, 'nTransmissions', 1); + [~, ~, ~, r] = get_system_parameters(Resource, Trans, T, r); + assert(abs(r.Fs/fc - rateFactors(m)) < 1e-6, ... + 'Wrong pre-quadrature sample rate for %s', modes{m}); + switch modes{m} + case 'BS50BW' + starts = 2*(0:n/2-1)/fc; + times = reshape([starts; starts + 1/(4*fc)], [], 1); + case 'BS100BW' + starts = (0:n/2-1)/fc; + times = reshape([starts; starts + 1/(4*fc)], [], 1); + case 'BS67BW' + times = (0:n-1)'*3/(4*fc); + case 'NS200BW' + times = (0:n-1)'/(4*fc); + end + for shift = [0, 1/(8*fc)] + T.transmitDelays = shift * ones(8,1); + s = echoframe_setup_fourier(P, T, r, rec); + assert(numel(s.frequencyAxis) == n/2); + assert(all(s.delayIndices(:) >= 0 & s.delayIndices(:) < n/2)); + for tone = fc*[0.875, 1, 1.125] + a = int16(round(16000*cos(2*pi*tone*times))); + b = int16(round(16000*cos(2*pi*tone*(times-shift)))); + a = fft(complex(single(a(1:2:end)), -single(a(2:2:end)))); + b = fft(complex(single(b(1:2:end)), -single(b(2:2:end)))); + [~, k] = max(abs(a)); + rotation = exp(1i*double(s.frequencyAxis(k))*double(s.planewaveDelays(1,2))); + assert(abs(double(b(k))/double(a(k))*rotation - 1) < 2e-4, ... + 'Fourier delay correction failed: %s, Fc=%g, n=%d, tone=%g', ... + modes{m}, fc, n, tone); + end + end + end + end +end +fprintf(' all four sampling modes cancel known tone delays OK\n'); +end From 5af8db1a6bb355fe5f0d878ef3124299a74f1196 Mon Sep 17 00:00:00 2001 From: Samuel Le Meur-Diebolt Date: Wed, 30 Sep 2026 12:37:28 +0100 Subject: [PATCH 02/11] feat(das): add ffdas-backed delay-and-sum beamformer Add a DAS beamformer alongside the existing Fourier imaging path, backed by the ffdas submodule (pinned to ea3a25e). FFDASBeamformer owns the ffdas handle and tensor descriptors, applies the TGC vector to the IQ data, and calls ffdas_das with the channel and voxel geometry tables plus optional source directivity. DASReconSpec carries that geometry: channel and voxel positions, per-tx offsets, weights, the TGC vector, the wavenum, and the algorithm and compute-type selectors. Thread the new spec through the resource bundle, the core, the MEX and Python conversions, the Python wrapper, and the CLI MAT loader, so DAS is selectable everywhere the beamformer already was. On the MATLAB side add echoframe_setup_das to build the geometry tables, relax echoframe_validate_structs to require the Fourier tables only for Fourier, and let initialize_image_reconstruction dispatch on method. Add compare_beamformers.py for Fourier-versus-DAS comparison on a stored RF frame, and generate_synthetic_session.py to produce a session to compare against. Resolve and link ffdas from CMake, and bump the dev image to CMake 3.26 or newer so the submodule can be built with add_subdirectory. Also fall back to a conservative 4096-byte O_DIRECT alignment in LinuxFileIO when the kernel does not report STATX_DIOALIGN (Linux older than 6.1) instead of throwing. --- .gitmodules | 3 + docker/ubuntu22.04-cuda12.8.Dockerfile | 2 +- echoframe/cpp/libs/README.md | 13 +- echoframe/cpp/libs/ffdas | 1 + echoframe/cpp/src/CMakeLists.txt | 50 +++- .../cpp/src/beamformer/das/das_recon_spec.h | 32 +++ .../cpp/src/beamformer/das/ffdas_beamformer.h | 53 ++++ .../src/beamformer/das/ffdas_beamformer.t.hpp | 173 +++++++++++++ .../beamformer/echoframe_resources_bundle.h | 4 +- echoframe/cpp/src/beamformer/resources.h | 3 + .../src/echoframe_cli/load_specs_from_mat.cpp | 76 ++++-- echoframe/cpp/src/efcore/echoframe_core.cu | 98 ++++--- echoframe/cpp/src/efcore/echoframe_core.h | 34 +-- .../cpp/src/mex/mex_resources_conversions.cpp | 89 ++++++- .../cpp/src/mex/mex_resources_conversions.h | 11 +- .../src/python/echoframe_py_conversions.cpp | 50 +++- .../cpp/src/python/echoframe_py_wrapper.cpp | 16 +- .../matlab/core/imaging/echoframe_setup_das.m | 59 +++++ .../core/imaging/echoframe_validate_structs.m | 28 +- .../imaging/initialize_image_reconstruction.m | 11 +- .../python/core/validate_structs_py_to_ef.py | 26 +- .../python/examples/compare_beamformers.py | 240 ++++++++++++++++++ .../examples/generate_synthetic_session.py | 166 ++++++++++++ justfile | 3 + 24 files changed, 1139 insertions(+), 102 deletions(-) create mode 160000 echoframe/cpp/libs/ffdas create mode 100644 echoframe/cpp/src/beamformer/das/das_recon_spec.h create mode 100644 echoframe/cpp/src/beamformer/das/ffdas_beamformer.h create mode 100644 echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp create mode 100644 echoframe/matlab/core/imaging/echoframe_setup_das.m create mode 100644 echoframe/python/examples/compare_beamformers.py create mode 100644 echoframe/python/examples/generate_synthetic_session.py diff --git a/.gitmodules b/.gitmodules index 9390625..243aa31 100644 --- a/.gitmodules +++ b/.gitmodules @@ -2,3 +2,6 @@ path = echoframe/cpp/libs/GSL url = https://github.com/microsoft/GSL.git +[submodule "echoframe/cpp/libs/ffdas"] + path = echoframe/cpp/libs/ffdas + url = https://github.com/BrainEchoLab/ffdas.git diff --git a/docker/ubuntu22.04-cuda12.8.Dockerfile b/docker/ubuntu22.04-cuda12.8.Dockerfile index 708e4b2..4e10ff6 100644 --- a/docker/ubuntu22.04-cuda12.8.Dockerfile +++ b/docker/ubuntu22.04-cuda12.8.Dockerfile @@ -29,7 +29,6 @@ ENV PATH="${VCPKG_ROOT}:${CARGO_HOME}/bin:${PATH}" RUN apt-get update && apt-get install -y --no-install-recommends \ build-essential \ ca-certificates \ - cmake \ curl \ fontconfig \ git \ @@ -74,6 +73,7 @@ RUN apt-get update && apt-get install -y --no-install-recommends \ unzip \ xauth \ zip \ + && python3 -m pip install --no-cache-dir "cmake>=3.26,<4" \ && rm -rf /var/lib/apt/lists/* # Shell tooling for distrobox sessions that inherit host dotfiles. fzf comes from diff --git a/echoframe/cpp/libs/README.md b/echoframe/cpp/libs/README.md index 40bd8d7..7302e24 100644 --- a/echoframe/cpp/libs/README.md +++ b/echoframe/cpp/libs/README.md @@ -11,16 +11,21 @@ across the core for safer array views and pointer contracts. Header-only. Note: this is Microsoft's **Guidelines** Support Library, not the GNU Scientific Library. +## ffdas + +A submodule (`github.com/BrainEchoLab/ffdas`) providing the optional CUDA +DAS beamformer backend. EchoFrame only builds it when configured with +`-DEF_USE_FFDAS=ON`. + ## Storage An in-tree library (`Storage/`, not a submodule) that manages writing and reading real-time RF, BF and PDI data to disk. -## Updating the GSL submodule +## Updating submodules -GSL is the only submodule here (Storage is a plain in-tree directory). To -initialise or update it from the repository root: +To initialise or update external submodules from the repository root: ```bash -git submodule update --init echoframe/cpp/libs/GSL +git submodule update --init echoframe/cpp/libs/GSL echoframe/cpp/libs/ffdas ``` diff --git a/echoframe/cpp/libs/ffdas b/echoframe/cpp/libs/ffdas new file mode 160000 index 0000000..6466de3 --- /dev/null +++ b/echoframe/cpp/libs/ffdas @@ -0,0 +1 @@ +Subproject commit 6466de37dc3cb2f18946a2ef4543b788af783e95 diff --git a/echoframe/cpp/src/CMakeLists.txt b/echoframe/cpp/src/CMakeLists.txt index fb9d167..e2cb6c8 100644 --- a/echoframe/cpp/src/CMakeLists.txt +++ b/echoframe/cpp/src/CMakeLists.txt @@ -116,6 +116,50 @@ option(EF_BUILD_MEX "Build the MATLAB MEX interface" ON) option(EF_BUILD_PYTHON "Build the Python extension module" ON) option(ENABLE_CUDA_TIMING "Enable CUDA kernel timing instrumentation" OFF) +option(EF_USE_FFDAS "Build the ffdas-backed DAS beamformer" OFF) +set(EF_FFDAS_CUDA_ARCHITECTURES "" CACHE STRING + "CUDA architectures for the ffdas submodule fallback build") + +if(EF_USE_FFDAS) + find_package(ffdas CONFIG QUIET) + if(NOT TARGET ffdas::ffdas AND NOT TARGET ffdas) + set(EF_FFDAS_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/../libs/ffdas") + if(EXISTS "${EF_FFDAS_SOURCE}/CMakeLists.txt") + if(CMAKE_VERSION VERSION_LESS 3.26) + message(FATAL_ERROR + "EF_USE_FFDAS=ON needs CMake >= 3.26 to build the ffdas submodule. Update the dev image or install CMake from PyPI/Kitware.") + endif() + set(EF_SAVED_CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") + set(EF_SAVED_CUDA_SEPARABLE_COMPILATION + "${CMAKE_CUDA_SEPARABLE_COMPILATION}") + if(NOT EF_FFDAS_CUDA_ARCHITECTURES) + if(CUDAToolkit_VERSION VERSION_GREATER_EQUAL 13) + set(EF_FFDAS_CUDA_ARCHITECTURES "80;86;89;90;100;103;121") + else() + set(EF_FFDAS_CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") + endif() + endif() + set(CMAKE_CUDA_ARCHITECTURES "${EF_FFDAS_CUDA_ARCHITECTURES}") + set(CMAKE_CUDA_SEPARABLE_COMPILATION OFF) + add_subdirectory("${EF_FFDAS_SOURCE}" "${CMAKE_BINARY_DIR}/ffdas") + if(TARGET ffdas) + set_property(TARGET ffdas PROPERTY + CUDA_SEPARABLE_COMPILATION OFF) + endif() + set(CMAKE_CUDA_ARCHITECTURES "${EF_SAVED_CUDA_ARCHITECTURES}") + set(CMAKE_CUDA_SEPARABLE_COMPILATION + "${EF_SAVED_CUDA_SEPARABLE_COMPILATION}") + else() + message(FATAL_ERROR + "EF_USE_FFDAS=ON requires an installed ffdas package or the submodule at ${EF_FFDAS_SOURCE}") + endif() + endif() + if(TARGET ffdas::ffdas) + set(EF_FFDAS_TARGET ffdas::ffdas) + else() + set(EF_FFDAS_TARGET ffdas) + endif() +endif() # -------------------------------------------------------------------------- # Common include directories # -------------------------------------------------------------------------- @@ -177,8 +221,12 @@ target_include_directories(echoframe_core target_link_libraries(echoframe_core PRIVATE - CUDA::cudart + CUDA::cudart ) +if(EF_USE_FFDAS) + target_link_libraries(echoframe_core PRIVATE ${EF_FFDAS_TARGET}) + target_compile_definitions(echoframe_core PRIVATE EF_USE_FFDAS=1) +endif() target_compile_definitions(echoframe_core PRIVATE NOMINMAX diff --git a/echoframe/cpp/src/beamformer/das/das_recon_spec.h b/echoframe/cpp/src/beamformer/das/das_recon_spec.h new file mode 100644 index 0000000..f1e0962 --- /dev/null +++ b/echoframe/cpp/src/beamformer/das/das_recon_spec.h @@ -0,0 +1,32 @@ +/** + * @file das_recon_spec.h + * @brief Delay-and-sum reconstruction parameters. + */ +#pragma once + +#include + +namespace Beamform { + +struct DASReconSpec { + float *channelPositions{}; ///< Host pointer, (nActiveChannels, 3). + float *d_channelPositions{}; ///< Device pointer, (nActiveChannels, 3). + float *voxelPositions{}; ///< Host pointer, (totalSize, 3), z fastest. + float *d_voxelPositions{}; ///< Device pointer, (totalSize, 3). + float *offsets{}; ///< Host pointer, (nTX, totalSize). + float *d_offsets{}; ///< Device pointer, (nTX, totalSize). + float *weights{}; ///< Host pointer, (nTX, totalSize). + float *d_weights{}; ///< Device pointer, (nTX, totalSize). + float *tgcVector{}; ///< Host pointer, (nSamplesIQ). + float *d_tgcVector{}; ///< Device pointer, (nSamplesIQ). + float + *sourceDirections{}; ///< Optional host pointer, (nActiveChannels, 4). + float *d_sourceDirections{}; ///< Optional device pointer, + ///< (nActiveChannels, 4). + float wavenum{0.0f}; ///< IQ phase rotation, -2*pi*Fc/Fs. + int32_t algorithm{1}; ///< ffdas algorithm: 1=ALG1, 2=ALG2, 4=ALG4. + int32_t computeType{0}; ///< ffdas compute type: 0=default. + bool useDirectivity{false}; +}; + +} // namespace Beamform diff --git a/echoframe/cpp/src/beamformer/das/ffdas_beamformer.h b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.h new file mode 100644 index 0000000..49f5195 --- /dev/null +++ b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.h @@ -0,0 +1,53 @@ +/** + * @file ffdas_beamformer.h + * @brief ffdas-backed delay-and-sum beamformer. + */ +#pragma once + +#include + +#include "../beamformer.h" +#include "das_recon_spec.h" + +#ifdef EF_USE_FFDAS +#include +#endif + +namespace Beamform { + +template +class FFDASBeamformer : public Beamformer { + private: + DASReconSpec dasReconSpec; + +#ifdef EF_USE_FFDAS + ffdas_handle_t handle{}; + ffdas_tensor_desc_t xDesc{}; + ffdas_tensor_desc_t outDesc{}; +#endif + + bool mInitialized{false}; + + public: + FFDASBeamformer() = default; + + FFDASBeamformer(ReceiveSpec pReceiveSpec, ReconSpec pReconSpec, + DASReconSpec pDASReconSpec, bfType_t *d_pRF) + : Beamformer(std::move(pReceiveSpec), pReconSpec, d_pRF), + dasReconSpec(pDASReconSpec) { + initialize(); + mInitialized = true; + } + + ~FFDASBeamformer(); + + FFDASBeamformer(FFDASBeamformer &&x) noexcept = delete; + FFDASBeamformer &operator=(FFDASBeamformer &&x) noexcept = delete; + + void initialize(); + void initGPU(); + void clearGPU(); + void process(); +}; + +} // namespace Beamform diff --git a/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp new file mode 100644 index 0000000..9cf59cc --- /dev/null +++ b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp @@ -0,0 +1,173 @@ +/** + * @file ffdas_beamformer.t.hpp + * @brief ffdas-backed delay-and-sum beamformer implementation. + */ +#pragma once + +#include "ffdas_beamformer.h" + +#include +#include +#include +#include + +#include "../../cuda/cuda_error.h" +#include "../fourier_imaging/fourier_imaging_kernels.cuh" + +namespace Beamform { + +#ifdef EF_USE_FFDAS +namespace { +inline void ffdasErrchk(ffdas_error_t err) { + if (err != FFDAS_SUCCESS) { + throw std::runtime_error(std::string("ffdas error: ") + + ffdas_error_string(err)); + } +} +} // namespace +#endif + +template +FFDASBeamformer::~FFDASBeamformer() { + if (mInitialized) clearGPU(); +} + +template +void FFDASBeamformer::initialize() { +#ifndef EF_USE_FFDAS + throw std::runtime_error( + "DAS beamforming requires configuring with -DEF_USE_FFDAS=ON."); +#else + static_assert( + std::is_same_v, + "ffdas beamformer currently supports complex single IQ only."); + initGPU(); +#endif +} + +template +void FFDASBeamformer::initGPU() { +#ifndef EF_USE_FFDAS + throw std::runtime_error( + "DAS beamforming requires configuring with -DEF_USE_FFDAS=ON."); +#else + const size_t nChannels = this->receiveSpec.nActiveChannels; + const size_t nVoxels = this->reconSpec.totalSize; + const size_t nTx = this->receiveSpec.nTX; + + if (!dasReconSpec.channelPositions || !dasReconSpec.voxelPositions || + !dasReconSpec.offsets || !dasReconSpec.weights || + !dasReconSpec.tgcVector) { + throw std::runtime_error( + "DASReconSpec requires channelPositions, voxelPositions, offsets, " + "weights, and tgcVector."); + } + + gpuErrchk(cudaMalloc(&dasReconSpec.d_channelPositions, + nChannels * 3 * sizeof(float))); + gpuErrchk(cudaMemcpy( + dasReconSpec.d_channelPositions, dasReconSpec.channelPositions, + nChannels * 3 * sizeof(float), cudaMemcpyHostToDevice)); + + gpuErrchk(cudaMalloc(&dasReconSpec.d_voxelPositions, + nVoxels * 3 * sizeof(float))); + gpuErrchk(cudaMemcpy(dasReconSpec.d_voxelPositions, + dasReconSpec.voxelPositions, + nVoxels * 3 * sizeof(float), cudaMemcpyHostToDevice)); + + gpuErrchk( + cudaMalloc(&dasReconSpec.d_offsets, nTx * nVoxels * sizeof(float))); + gpuErrchk(cudaMemcpy(dasReconSpec.d_offsets, dasReconSpec.offsets, + nTx * nVoxels * sizeof(float), + cudaMemcpyHostToDevice)); + + gpuErrchk( + cudaMalloc(&dasReconSpec.d_weights, nTx * nVoxels * sizeof(float))); + gpuErrchk(cudaMemcpy(dasReconSpec.d_weights, dasReconSpec.weights, + nTx * nVoxels * sizeof(float), + cudaMemcpyHostToDevice)); + + gpuErrchk(cudaMalloc(&dasReconSpec.d_tgcVector, + this->receiveSpec.nSamplesIQ * sizeof(float))); + gpuErrchk(cudaMemcpy(dasReconSpec.d_tgcVector, dasReconSpec.tgcVector, + this->receiveSpec.nSamplesIQ * sizeof(float), + cudaMemcpyHostToDevice)); + + if (dasReconSpec.useDirectivity) { + if (!dasReconSpec.sourceDirections) + throw std::runtime_error( + "DASReconSpec.useDirectivity requires sourceDirections."); + gpuErrchk(cudaMalloc(&dasReconSpec.d_sourceDirections, + nChannels * 4 * sizeof(float))); + gpuErrchk(cudaMemcpy( + dasReconSpec.d_sourceDirections, dasReconSpec.sourceDirections, + nChannels * 4 * sizeof(float), cudaMemcpyHostToDevice)); + } + + ffdasErrchk(ffdas_create(&handle)); + + const int64_t xDims[] = { + this->receiveSpec.nRepeats, this->receiveSpec.nActiveChannels, + this->receiveSpec.nTX, this->receiveSpec.nSamplesIQ}; + const int64_t xStrides[] = { + static_cast(this->receiveSpec.nTX) * + this->receiveSpec.nSamplesIQ, + static_cast(this->receiveSpec.nRepeats) * + this->receiveSpec.nTX * this->receiveSpec.nSamplesIQ, + this->receiveSpec.nSamplesIQ, 1}; + ffdasErrchk( + ffdas_create_tensor_desc(&xDesc, 4, xDims, xStrides, FFDAS_C_32F)); + + const int64_t outDims[] = {this->reconSpec.ensembleSize, this->reconSpec.nx, + this->reconSpec.nz}; + const int64_t outStrides[] = {this->reconSpec.totalSize, this->reconSpec.nz, + 1}; + ffdasErrchk(ffdas_create_tensor_desc(&outDesc, 3, outDims, outStrides, + FFDAS_C_32F)); +#endif +} + +template +void FFDASBeamformer::clearGPU() { +#ifdef EF_USE_FFDAS + if (xDesc) ffdas_destroy_tensor_desc(xDesc); + if (outDesc) ffdas_destroy_tensor_desc(outDesc); + if (handle) ffdas_destroy(handle); + if (dasReconSpec.d_channelPositions) + cudaFree(dasReconSpec.d_channelPositions); + if (dasReconSpec.d_voxelPositions) cudaFree(dasReconSpec.d_voxelPositions); + if (dasReconSpec.d_offsets) cudaFree(dasReconSpec.d_offsets); + if (dasReconSpec.d_weights) cudaFree(dasReconSpec.d_weights); + if (dasReconSpec.d_tgcVector) cudaFree(dasReconSpec.d_tgcVector); + if (dasReconSpec.d_sourceDirections) + cudaFree(dasReconSpec.d_sourceDirections); +#endif +} + +template +void FFDASBeamformer::process() { +#ifndef EF_USE_FFDAS + throw std::runtime_error( + "DAS beamforming requires configuring with -DEF_USE_FFDAS=ON."); +#else + constexpr int threadsPerBlock = 512; + const int nRF = this->receiveSpec.nFastTimeSamples * this->receiveSpec.nTX * + this->receiveSpec.nSlowTimeSamples * + this->receiveSpec.nActiveChannels; + TGCWeighting_kernel<<>>( + this->d_RF, dasReconSpec.d_tgcVector, this->receiveSpec); + gpuErrchk(cudaPeekAtLastError()); + + const float2 beta{0.0f, 0.0f}; + ffdasErrchk(ffdas_das( + handle, dasReconSpec.d_channelPositions, + dasReconSpec.useDirectivity ? dasReconSpec.d_sourceDirections : nullptr, + dasReconSpec.wavenum, xDesc, this->d_RF, dasReconSpec.d_voxelPositions, + dasReconSpec.d_offsets, dasReconSpec.d_weights, &beta, outDesc, + this->d_BF, static_cast(dasReconSpec.computeType), + static_cast(dasReconSpec.algorithm))); + gpuErrchk(cudaDeviceSynchronize()); +#endif +} + +} // namespace Beamform diff --git a/echoframe/cpp/src/beamformer/echoframe_resources_bundle.h b/echoframe/cpp/src/beamformer/echoframe_resources_bundle.h index ac755f3..ba852a5 100644 --- a/echoframe/cpp/src/beamformer/echoframe_resources_bundle.h +++ b/echoframe/cpp/src/beamformer/echoframe_resources_bundle.h @@ -21,6 +21,7 @@ #include "../../libs/Storage/src/storage/storage_spec.h" #include "../pdi/pdi_spec.h" +#include "./das/das_recon_spec.h" #include "./fourier_imaging/fourier_recon_spec.h" #include "resources.h" @@ -32,7 +33,8 @@ struct EchoframeResources { PDI::PDISpec pdiSpec; ///< Power Doppler Imaging specification. Beamform::FourierReconSpec fourierReconSpec; ///< Fourier reconstruction specification. - Storage::StorageSpec storageSpec; ///< General storage specification. + Beamform::DASReconSpec dasReconSpec; ///< DAS reconstruction specification. + Storage::StorageSpec storageSpec; ///< General storage specification. Storage::StorageSpec pdiStorageSpec; ///< Storage specification for PDI data. Storage::StorageSpec diff --git a/echoframe/cpp/src/beamformer/resources.h b/echoframe/cpp/src/beamformer/resources.h index 9ebb032..7e09714 100644 --- a/echoframe/cpp/src/beamformer/resources.h +++ b/echoframe/cpp/src/beamformer/resources.h @@ -21,6 +21,8 @@ namespace Beamform { +enum class BeamformerType : int32_t { Fourier = 0, DAS = 1 }; + /** * @brief Structure describing RF data acquisition and transducer properties. */ @@ -74,6 +76,7 @@ struct ReconSpec { bool getBF{false}; ///< Output beamformed data. bool getPDI{false}; ///< Output PDI data. bool cropBF{false}; ///< Crop beamformed data. + BeamformerType beamformerType{BeamformerType::Fourier}; bool initialized{false}; ///< Initialization flag. }; diff --git a/echoframe/cpp/src/echoframe_cli/load_specs_from_mat.cpp b/echoframe/cpp/src/echoframe_cli/load_specs_from_mat.cpp index 6d79098..b12ba87 100644 --- a/echoframe/cpp/src/echoframe_cli/load_specs_from_mat.cpp +++ b/echoframe/cpp/src/echoframe_cli/load_specs_from_mat.cpp @@ -77,6 +77,19 @@ static void validateFourierMatStruct(matvar_t *s) { validateField(s, "tgcVector"); } +static void validateDASMatStruct(matvar_t *s) { + validateField(s, "dasChannelPositions"); + validateField(s, "dasVoxelPositions"); + validateField(s, "dasOffsets"); + validateField(s, "dasWeights"); + validateField(s, "dasWavenum"); +} + +static bool hasField(matvar_t *s, const char *fld) { + matvar_t *f = Mat_VarGetStructFieldByName(s, fld, 0); + return f && f->data; +} + template static T scalar(matvar_t *s, const char *fld) { matvar_t *f = Mat_VarGetStructFieldByName(s, fld, 0); @@ -260,25 +273,52 @@ void loadSpecsFromMat(const std::string &matFile, EchoframeResources &res) { P.getPDI = scalar(Pm.get(), "getPDI"); P.cropBF = scalar(Pm.get(), "cropBF"); P.ensembleSize = R.nRepeats; + if (hasField(Pm.get(), "beamformerType")) { + std::string bf = str(Pm.get(), "beamformerType"); + std::transform(bf.begin(), bf.end(), bf.begin(), ::tolower); + if (bf == "das") P.beamformerType = Beamform::BeamformerType::DAS; + } - validateFourierMatStruct(Pm.get()); - - auto &F = res.fourierReconSpec; - - static auto delayIdx = vecRowMajor( - Mat_VarGetStructFieldByName(Pm.get(), "delayIndices", 0)); - - static auto interpW = vecCplxRowMajor( - Mat_VarGetStructFieldByName(Pm.get(), "interpolationWeights", 0)); - static auto freqAxis = vec(Pm.get(), "frequencyAxis"); - static auto pwDelays = vec(Pm.get(), "planewaveDelays"); - static auto tgc = vec(Pm.get(), "tgcVector"); - - F.delayIndices = delayIdx.data(); - F.interpolationWeights = interpW.data(); - F.frequencyAxis = freqAxis.data(); - F.planewaveDelays = pwDelays.data(); - F.tgcVector = tgc.data(); + if (P.beamformerType == Beamform::BeamformerType::DAS) { + validateDASMatStruct(Pm.get()); + auto &DAS = res.dasReconSpec; + static auto channelPos = vec(Pm.get(), "dasChannelPositions"); + static auto voxelPos = vec(Pm.get(), "dasVoxelPositions"); + static auto offsets = vec(Pm.get(), "dasOffsets"); + static auto weights = vec(Pm.get(), "dasWeights"); + static auto tgc = vec(Pm.get(), "tgcVector"); + DAS.channelPositions = channelPos.data(); + DAS.voxelPositions = voxelPos.data(); + DAS.offsets = offsets.data(); + DAS.weights = weights.data(); + DAS.tgcVector = tgc.data(); + DAS.wavenum = scalar(Pm.get(), "dasWavenum"); + DAS.algorithm = hasField(Pm.get(), "dasAlgorithm") + ? scalar(Pm.get(), "dasAlgorithm") + : 1; + DAS.computeType = hasField(Pm.get(), "dasComputeType") + ? scalar(Pm.get(), "dasComputeType") + : 0; + } else { + validateFourierMatStruct(Pm.get()); + + auto &F = res.fourierReconSpec; + + static auto delayIdx = vecRowMajor( + Mat_VarGetStructFieldByName(Pm.get(), "delayIndices", 0)); + + static auto interpW = vecCplxRowMajor( + Mat_VarGetStructFieldByName(Pm.get(), "interpolationWeights", 0)); + static auto freqAxis = vec(Pm.get(), "frequencyAxis"); + static auto pwDelays = vec(Pm.get(), "planewaveDelays"); + static auto tgc = vec(Pm.get(), "tgcVector"); + + F.delayIndices = delayIdx.data(); + F.interpolationWeights = interpW.data(); + F.frequencyAxis = freqAxis.data(); + F.planewaveDelays = pwDelays.data(); + F.tgcVector = tgc.data(); + } /* ------------ PDISpec ------------ */ auto &D = res.pdiSpec; diff --git a/echoframe/cpp/src/efcore/echoframe_core.cu b/echoframe/cpp/src/efcore/echoframe_core.cu index a955c5f..72a8596 100644 --- a/echoframe/cpp/src/efcore/echoframe_core.cu +++ b/echoframe/cpp/src/efcore/echoframe_core.cu @@ -42,6 +42,7 @@ #include "../../libs/Storage/src/storage/Handler.t.hpp" #include "../beamformer/beamformer.t.hpp" +#include "../beamformer/das/ffdas_beamformer.t.hpp" #include "../beamformer/fourier_imaging/fourier_imaging.t.hpp" #include "../cuda/cuda_error.h" @@ -92,8 +93,8 @@ void refreshStorageSlotRings() { // during a disk stall that the ringed streams came through clean. // // RF has no ring. Staging it costs a whole frame per slot -- about 50 ms on - // the clinic protocol, which put the loop over its acquisition period -- and - // the acquisition ring already does the job for free: the hardware only + // the clinic protocol, which put the loop over its acquisition period -- + // and the acquisition ring already does the job for free: the hardware only // reuses a receive frame after cycling through the others, so a deep enough // ring outlasts the write. init_storage warns when the depths make that // impossible. @@ -104,7 +105,8 @@ void refreshStorageSlotRings() { const bool globalSet = (global != nullptr && global[0] != '\0'); const bool globalOn = globalSet && global[0] == '1'; - g_ringBF = envFlag("EF_STORAGE_SLOT_RINGS_BF", globalSet ? globalOn : false); + g_ringBF = + envFlag("EF_STORAGE_SLOT_RINGS_BF", globalSet ? globalOn : false); g_ringPDI = envFlag("EF_STORAGE_SLOT_RINGS_PDI", globalSet ? globalOn : true); g_ringTag = @@ -130,9 +132,12 @@ void refreshStorageSlotRings() { bool storageSlotRingsEnabled(StorageStream stream) { if (!g_ringsRead) refreshStorageSlotRings(); switch (stream) { - case StorageStream::BF: return g_ringBF; - case StorageStream::PDI: return g_ringPDI; - case StorageStream::RFTimeTag: return g_ringTag; + case StorageStream::BF: + return g_ringBF; + case StorageStream::PDI: + return g_ringPDI; + case StorageStream::RFTimeTag: + return g_ringTag; } return false; } @@ -159,7 +164,8 @@ EchoFrameCore::EchoFrameCore(const EchoframeResources &res, bool useStorage) // let a violation surface later as an aio_error two layers away. mTimeTagsCount = static_cast(res.receiveSpec.nRepeats) * res.receiveSpec.nTX; - resizeTimeTagRing(storageSlotsFor(res.rfTimeTagStorageSpec, StorageStream::RFTimeTag)); + resizeTimeTagRing( + storageSlotsFor(res.rfTimeTagStorageSpec, StorageStream::RFTimeTag)); // Initialization banner banner() << "+--------------------------- INITIALIZING ECHOFRAME " @@ -187,7 +193,8 @@ EchoFrameCore::EchoFrameCore(const EchoframeResources &res, bool useStorage) } pdiObject = new PDI::PDI(resources.pdiSpec, resources.receiveSpec, resources.reconSpec, false); - pdiObject->setStorageSlots(storageSlotsFor(resources.pdiStorageSpec, StorageStream::PDI)); + pdiObject->setStorageSlots( + storageSlotsFor(resources.pdiStorageSpec, StorageStream::PDI)); banner() << "+----------------------------- ECHOFRAME INITIALIZED " "-----------------------------+" @@ -228,7 +235,8 @@ EchoFrameCore::EchoFrameCore(const PDI::PDISpec &pdiSpec, } pdiObject = new PDI::PDI(resources.pdiSpec, resources.receiveSpec, resources.reconSpec, true); - pdiObject->setStorageSlots(storageSlotsFor(resources.pdiStorageSpec, StorageStream::PDI)); + pdiObject->setStorageSlots( + storageSlotsFor(resources.pdiStorageSpec, StorageStream::PDI)); if (storageInitialized && resources.pdiStorageSpec.save) { pdiStorageHandler = Storage::Handler( @@ -266,9 +274,9 @@ void EchoFrameCore::resizeTimeTagRing(int nSlots) { mTimeTagsBase = nullptr; mTimeTags = nullptr; } - gpuErrchk(cudaMallocHost(&mTimeTagsBase, - static_cast(mTimeTagsSlots) * - mTimeTagsStride * sizeof(double))); + gpuErrchk( + cudaMallocHost(&mTimeTagsBase, static_cast(mTimeTagsSlots) * + mTimeTagsStride * sizeof(double))); mTimeTags = mTimeTagsBase; if (reinterpret_cast(mTimeTags) % kAssumedPageAlignment != 0) { @@ -437,9 +445,11 @@ void EchoFrameCore::reinitStorage(const EchoframeResources &newRes) { // Saving may have just been switched on; resize before updateCropSpecs() // reallocates the cropped ring. - bfFormatter.setStorageSlots(storageSlotsFor(resources.bfStorageSpec, StorageStream::BF), - resources.bfStorageSpec.crop); - resizeTimeTagRing(storageSlotsFor(resources.rfTimeTagStorageSpec, StorageStream::RFTimeTag)); + bfFormatter.setStorageSlots( + storageSlotsFor(resources.bfStorageSpec, StorageStream::BF), + resources.bfStorageSpec.crop); + resizeTimeTagRing(storageSlotsFor(resources.rfTimeTagStorageSpec, + StorageStream::RFTimeTag)); if (resources.bfStorageSpec.crop) { bfFormatter.updateCropSpecs(resources.reconSpec); @@ -455,7 +465,8 @@ void EchoFrameCore::reinitStorage(const EchoframeResources &newRes) { resources.reconSpec); } if (pdiObject != nullptr) - pdiObject->setStorageSlots(storageSlotsFor(resources.pdiStorageSpec, StorageStream::PDI)); + pdiObject->setStorageSlots( + storageSlotsFor(resources.pdiStorageSpec, StorageStream::PDI)); } void EchoFrameCore::reinitExperiment(const EchoframeResources &newRes) { @@ -478,11 +489,14 @@ void EchoFrameCore::reinitExperiment(const EchoframeResources &newRes) { rfTimeTagStorageHandler = std::move(newRFTimeTagStorageHandler); // The new specs may enable or disable saving; re-size the rings to match. - bfFormatter.setStorageSlots(storageSlotsFor(resources.bfStorageSpec, StorageStream::BF), - resources.bfStorageSpec.crop); + bfFormatter.setStorageSlots( + storageSlotsFor(resources.bfStorageSpec, StorageStream::BF), + resources.bfStorageSpec.crop); if (pdiObject != nullptr) - pdiObject->setStorageSlots(storageSlotsFor(resources.pdiStorageSpec, StorageStream::PDI)); - resizeTimeTagRing(storageSlotsFor(resources.rfTimeTagStorageSpec, StorageStream::RFTimeTag)); + pdiObject->setStorageSlots( + storageSlotsFor(resources.pdiStorageSpec, StorageStream::PDI)); + resizeTimeTagRing(storageSlotsFor(resources.rfTimeTagStorageSpec, + StorageStream::RFTimeTag)); } void EchoFrameCore::updatePDIThreshold(const float newThreshold) { @@ -497,15 +511,16 @@ void EchoFrameCore::updatePDILowerThreshold(const float newLowerThreshold) { EchoFrameTimings EchoFrameCore::getLastTimings() { EchoFrameTimings t{}; - t.rf_transfer = eventTimer.elapsedTime("RFTransfer"); - t.rf_formatting = eventTimer.elapsedTime("RFFormatting"); - t.beamforming = eventTimer.elapsedTime("Beamforming"); - t.bf_formatting = eventTimer.elapsedTime("BFFormatting"); + t.rf_transfer = eventTimer.elapsedTime("RFTransfer"); + t.rf_formatting = eventTimer.elapsedTime("RFFormatting"); + t.beamforming = eventTimer.elapsedTime("Beamforming"); + t.bf_formatting = eventTimer.elapsedTime("BFFormatting"); t.pdi_processing = eventTimer.elapsedTime("PDIProcessing"); - t.pdi_transfer = eventTimer.elapsedTime("PDITransfer"); - t.bf_storage = eventTimer.elapsedTime("BFStorage"); // queueing, not the write - t.pdi_storage = eventTimer.elapsedTime("PDIStorage"); - t.total = eventTimer.elapsedTime("TotalTime"); + t.pdi_transfer = eventTimer.elapsedTime("PDITransfer"); + t.bf_storage = + eventTimer.elapsedTime("BFStorage"); // queueing, not the write + t.pdi_storage = eventTimer.elapsedTime("PDIStorage"); + t.total = eventTimer.elapsedTime("TotalTime"); return t; } @@ -517,10 +532,10 @@ EchoFrameStorageStats EchoFrameCore::getStorageStats() const { s.pdi = makeStreamStats( pdiStorageHandler, resources.pdiStorageSpec, storageSlotsFor(resources.pdiStorageSpec, StorageStream::PDI)); - s.timetag = makeStreamStats( - rfTimeTagStorageHandler, resources.rfTimeTagStorageSpec, - storageSlotsFor(resources.rfTimeTagStorageSpec, - StorageStream::RFTimeTag)); + s.timetag = + makeStreamStats(rfTimeTagStorageHandler, resources.rfTimeTagStorageSpec, + storageSlotsFor(resources.rfTimeTagStorageSpec, + StorageStream::RFTimeTag)); return s; } @@ -530,23 +545,32 @@ void EchoFrameCore::prepare_beamform( Beamform::RFFormatter &rfFormatter, Beamform::BFFormatter &bfFormatter, Beamform::Beamformer **beamformerPtr) { - rfFormatter = std::move(Beamform::RFFormatter( - resources.receiveSpec)); + rfFormatter = std::move( + Beamform::RFFormatter(resources.receiveSpec)); bfFormatter = std::move(Beamform::BFFormatter( resources.receiveSpec, resources.reconSpec, resources.reconSpec.totalSize * resources.reconSpec.ensembleSize, - resources.reconSpec.totalSizeCropped * resources.reconSpec.ensembleSize)); + resources.reconSpec.totalSizeCropped * + resources.reconSpec.ensembleSize)); // Size the BF ring before initializeResources() allocates the cropped one. - bfFormatter.setStorageSlots(storageSlotsFor(resources.bfStorageSpec, StorageStream::BF), - resources.bfStorageSpec.crop); + bfFormatter.setStorageSlots( + storageSlotsFor(resources.bfStorageSpec, StorageStream::BF), + resources.bfStorageSpec.crop); bfFormatter.initializeResources(); // Heap-allocate so each init gets a fresh instance sized for the current // specs; a static instance could only be sized once per MEX load. Owned by // EchoFrameCore::beamformerPtr and deleted in the destructor. + if (resources.reconSpec.beamformerType == Beamform::BeamformerType::DAS) { + *beamformerPtr = new Beamform::FFDASBeamformer( + std::move(resources.receiveSpec), resources.reconSpec, + std::move(resources.dasReconSpec), rfFormatter.getFormattedRF()); + return; + } + *beamformerPtr = new Beamform::FourierImaging( std::move(resources.receiveSpec), resources.reconSpec, std::move(resources.fourierReconSpec), rfFormatter.getFormattedRF()); diff --git a/echoframe/cpp/src/efcore/echoframe_core.h b/echoframe/cpp/src/efcore/echoframe_core.h index 154b5f7..a9b682c 100644 --- a/echoframe/cpp/src/efcore/echoframe_core.h +++ b/echoframe/cpp/src/efcore/echoframe_core.h @@ -7,7 +7,7 @@ * The EchoFrameCore class is responsible for managing the processing pipeline * of the EchoFrame system, including beamforming, PDI processing, and storage * management. It provides methods for initializing the system, processing RF - * data, and retrieving outputs. + * data, and retrieving outputs. * @version 0.1 * @date 2025-06-20 * @@ -26,6 +26,7 @@ #include "../../libs/Storage/src/storage/Handler.h" #include "../beamformer/BF_formatter.h" #include "../beamformer/RF_formatter.h" +#include "../beamformer/das/ffdas_beamformer.h" #include "../beamformer/echoframe_resources_bundle.h" #include "../beamformer/fourier_imaging/fourier_imaging.h" #include "../cuda/cuda_event_timer.hpp" @@ -99,19 +100,20 @@ struct EchoFrameTimings { * the disk held the source buffer, and how close the queue came to full. */ struct EchoFrameStreamStats { - bool saving; ///< Whether this stream is being written. - unsigned long long writes; ///< Completed writes. + bool saving; ///< Whether this stream is being written. + unsigned long long writes; ///< Completed writes. unsigned long long buffersQueued; ///< storeBuffer calls accepted. - double latencyMeanMs; ///< Mean queue-to-completion time. - double latencyMaxMs; ///< Worst queue-to-completion time. - int peakInFlight; ///< Highest concurrent outstanding writes. - int queueCapacity; ///< nBuffers-1: most that may be outstanding. - int slotRingDepth; ///< Producer slot-ring depth: nBuffers - ///< with the ring on, 1 with it off. - double blockedTotalMs; ///< Time storeBuffer waited for a slot. - double blockedMaxMs; ///< Worst single wait. - unsigned long long verified; ///< Writes checked by EF_STORAGE_VERIFY. - unsigned long long corrupted; ///< Of those, sources that changed in flight. + double latencyMeanMs; ///< Mean queue-to-completion time. + double latencyMaxMs; ///< Worst queue-to-completion time. + int peakInFlight; ///< Highest concurrent outstanding writes. + int queueCapacity; ///< nBuffers-1: most that may be outstanding. + int slotRingDepth; ///< Producer slot-ring depth: nBuffers + ///< with the ring on, 1 with it off. + double blockedTotalMs; ///< Time storeBuffer waited for a slot. + double blockedMaxMs; ///< Worst single wait. + unsigned long long verified; ///< Writes checked by EF_STORAGE_VERIFY. + unsigned long long + corrupted; ///< Of those, sources that changed in flight. }; /** @@ -136,7 +138,8 @@ inline EchoFrameStreamStats makeStreamStats( EchoFrameStreamStats s{}; s.saving = spec.save; s.writes = w.completed; - s.buffersQueued = static_cast(handler.getBuffersQueued()); + s.buffersQueued = + static_cast(handler.getBuffersQueued()); s.latencyMeanMs = w.latencyMeanMs(); s.latencyMaxMs = w.latencyMaxMs; s.peakInFlight = w.outstandingHighWater; @@ -267,7 +270,8 @@ class EchoFrameCore { // (aio_write) never races a deallocation. cudaMallocHost rather than // std::vector because O_DIRECT needs a page-aligned address, which // alignof(double) does not give; the constructor verifies the alignment. - // mTimeTags is a view into the ring, not the allocation: free mTimeTagsBase. + // mTimeTags is a view into the ring, not the allocation: free + // mTimeTagsBase. // PDI slot process() stored, so processAndGetOutputs returns that buffer // rather than calling getResults() again -- each call advances the PDI diff --git a/echoframe/cpp/src/mex/mex_resources_conversions.cpp b/echoframe/cpp/src/mex/mex_resources_conversions.cpp index 8e37bf9..2146553 100644 --- a/echoframe/cpp/src/mex/mex_resources_conversions.cpp +++ b/echoframe/cpp/src/mex/mex_resources_conversions.cpp @@ -18,6 +18,8 @@ #include +#include +#include #include #include #include @@ -50,6 +52,18 @@ bool isComplexSingle(const mxArray *array) { return mxIsComplex(array) && mxIsSingle(array); } +Beamform::BeamformerType getBeamformerType(const mxArray *reconStruct) { + const mxArray *field = + mxGetField(const_cast(reconStruct), 0, "beamformerType"); + if (!field) return Beamform::BeamformerType::Fourier; + std::string value(mxArrayToString(field)); + std::transform(value.begin(), value.end(), value.begin(), ::tolower); + if (value == "fourier") return Beamform::BeamformerType::Fourier; + if (value == "das") return Beamform::BeamformerType::DAS; + throw std::invalid_argument( + "ReconSpec.beamformerType must be 'Fourier' or 'DAS'."); +} + // Validate ReceiveSpec struct void validateReceiveStruct(const mxArray *s) { checkField(s, "nSamples", mxIsInt32, "int32"); @@ -70,8 +84,13 @@ void convertMexResources(const mxArray *prhs[], EchoframeResources &res) { convertPDISpec(res.pdiSpec, res.receiveSpec, res.reconSpec, prhs[INIT_INPUT_PDI_POS]); - convertFourierReconSpecStructs(res.fourierReconSpec, + if (res.reconSpec.beamformerType == Beamform::BeamformerType::DAS) { + convertDASReconSpecStructs(res.dasReconSpec, prhs[INIT_INPUT_RECON_POS]); + } else { + convertFourierReconSpecStructs(res.fourierReconSpec, + prhs[INIT_INPUT_RECON_POS]); + } disableStorageSpec(res.storageSpec); } @@ -84,8 +103,13 @@ void convertMexResourcesWithStorage(const mxArray *prhs[], int nrhs, convertPDISpec(res.pdiSpec, res.receiveSpec, res.reconSpec, prhs[INIT_INPUT_PDI_POS]); - convertFourierReconSpecStructs(res.fourierReconSpec, + if (res.reconSpec.beamformerType == Beamform::BeamformerType::DAS) { + convertDASReconSpecStructs(res.dasReconSpec, prhs[INIT_INPUT_RECON_POS]); + } else { + convertFourierReconSpecStructs(res.fourierReconSpec, + prhs[INIT_INPUT_RECON_POS]); + } // Storage initialization convertStorageStruct(res.bfStorageSpec, prhs[INIT_INPUT_BFSTORE_POS]); convertStorageStruct(res.pdiStorageSpec, prhs[INIT_INPUT_PDISTORE_POS]); @@ -98,7 +122,8 @@ void convertMexResourcesWithStorage(const mxArray *prhs[], int nrhs, } } -void convertReinitStorage(const mxArray *prhs[], int nrhs, EchoframeResources &res) { +void convertReinitStorage(const mxArray *prhs[], int nrhs, + EchoframeResources &res) { // Initialization code for Receive, Recon, and PDI structs convertStorageStruct(res.bfStorageSpec, prhs[REINIT_INPUT_BFSTORE_POS]); convertStorageStruct(res.pdiStorageSpec, prhs[REINIT_INPUT_PDISTORE_POS]); @@ -114,8 +139,13 @@ void convertReinitStorage(const mxArray *prhs[], int nrhs, EchoframeResources &r if (res.bfStorageSpec.crop) { convertReconStruct(res.reconSpec, res.receiveSpec, prhs[REINIT_INPUT_RECON_POS]); - convertFourierReconSpecStructs(res.fourierReconSpec, + if (res.reconSpec.beamformerType == Beamform::BeamformerType::DAS) { + convertDASReconSpecStructs(res.dasReconSpec, prhs[REINIT_INPUT_RECON_POS]); + } else { + convertFourierReconSpecStructs(res.fourierReconSpec, + prhs[REINIT_INPUT_RECON_POS]); + } } if (res.pdiStorageSpec.crop) { convertReconStruct(res.reconSpec, res.receiveSpec, @@ -125,7 +155,8 @@ void convertReinitStorage(const mxArray *prhs[], int nrhs, EchoframeResources &r } } -void convertReinitExperiment(const mxArray *prhs[], int nrhs, EchoframeResources &res) { +void convertReinitExperiment(const mxArray *prhs[], int nrhs, + EchoframeResources &res) { // Initialization code for Receive, Recon, and PDI structs convertStorageStruct(res.bfStorageSpec, prhs[REINIT_INPUT_BFSTORE_POS]); convertStorageStruct(res.pdiStorageSpec, prhs[REINIT_INPUT_PDISTORE_POS]); @@ -217,6 +248,7 @@ void convertReconStruct(Beamform::ReconSpec &reconSpec, reconSpec.getBF = mxGetLogicals(mxGetField(reconStruct, 0, "getBF"))[0]; reconSpec.getPDI = mxGetLogicals(mxGetField(reconStruct, 0, "getPDI"))[0]; reconSpec.cropBF = mxGetLogicals(mxGetField(reconStruct, 0, "cropBF"))[0]; + reconSpec.beamformerType = getBeamformerType(reconStruct); reconSpec.ensembleSize = receiveSpec.nRepeats; } @@ -241,9 +273,9 @@ void convertPDISpec(PDI::PDISpec &pdiSpec, // Optional lower (noise) threshold: default 0 when the field is absent, so // specs / ScanParameters.mat written before this field still init. const mxArray *lowerField = mxGetField(pdiStruct, 0, "lowerThreshold"); - pdiSpec.lowerThreshold = - (lowerField && mxIsSingle(lowerField)) ? mxGetSingles(lowerField)[0] - : 0.0f; + pdiSpec.lowerThreshold = (lowerField && mxIsSingle(lowerField)) + ? mxGetSingles(lowerField)[0] + : 0.0f; pdiSpec.shiftSize = mxGetInt32s(mxGetField(pdiStruct, 0, "shiftSize"))[0]; pdiSpec.cropPDI = mxGetLogicals(mxGetField(pdiStruct, 0, "cropPDI"))[0]; @@ -288,6 +320,47 @@ void convertFourierReconSpecStructs( mxGetSingles(mxGetField(reconStruct, 0, "tgcVector")); } +void validateDASReconStruct(const mxArray *s) { + checkField(s, "dasChannelPositions", mxIsSingle, "single array"); + checkField(s, "dasVoxelPositions", mxIsSingle, "single array"); + checkField(s, "dasOffsets", mxIsSingle, "single array"); + checkField(s, "dasWeights", mxIsSingle, "single array"); + checkField(s, "dasWavenum", mxIsSingle, "single scalar"); + checkField(s, "tgcVector", mxIsSingle, "single array"); +} + +void convertDASReconSpecStructs(Beamform::DASReconSpec &dasReconSpec, + const mxArray *reconStruct) { + validateDASReconStruct(reconStruct); + dasReconSpec.channelPositions = + mxGetSingles(mxGetField(reconStruct, 0, "dasChannelPositions")); + dasReconSpec.voxelPositions = + mxGetSingles(mxGetField(reconStruct, 0, "dasVoxelPositions")); + dasReconSpec.offsets = + mxGetSingles(mxGetField(reconStruct, 0, "dasOffsets")); + dasReconSpec.weights = + mxGetSingles(mxGetField(reconStruct, 0, "dasWeights")); + dasReconSpec.tgcVector = + mxGetSingles(mxGetField(reconStruct, 0, "tgcVector")); + dasReconSpec.wavenum = + mxGetSingles(mxGetField(reconStruct, 0, "dasWavenum"))[0]; + + const mxArray *alg = + mxGetField(const_cast(reconStruct), 0, "dasAlgorithm"); + dasReconSpec.algorithm = (alg && mxIsInt32(alg)) ? mxGetInt32s(alg)[0] : 1; + + const mxArray *compute = + mxGetField(const_cast(reconStruct), 0, "dasComputeType"); + dasReconSpec.computeType = + (compute && mxIsInt32(compute)) ? mxGetInt32s(compute)[0] : 0; + + const mxArray *dirs = mxGetField(const_cast(reconStruct), 0, + "dasSourceDirections"); + dasReconSpec.useDirectivity = dirs && !mxIsEmpty(dirs); + dasReconSpec.sourceDirections = + dasReconSpec.useDirectivity ? mxGetSingles(dirs) : nullptr; +} + // Validate StorageSpec struct void validateStorageStruct(const mxArray *s) { checkField(s, "save", mxIsLogical, "logical"); diff --git a/echoframe/cpp/src/mex/mex_resources_conversions.h b/echoframe/cpp/src/mex/mex_resources_conversions.h index 4a4f40f..1919fb8 100644 --- a/echoframe/cpp/src/mex/mex_resources_conversions.h +++ b/echoframe/cpp/src/mex/mex_resources_conversions.h @@ -93,7 +93,8 @@ void convertMexResourcesWithStorage(const mxArray *prhs[], int nrhs, * REINIT_INPUT_RFSTORE_POS, and RF storage is disabled without it. * @param res Reference to the native EchoframeResources structure to update. */ -void convertReinitStorage(const mxArray *prhs[], int nrhs, EchoframeResources &res); +void convertReinitStorage(const mxArray *prhs[], int nrhs, + EchoframeResources &res); /** * @brief Converts MATLAB input arrays for re-initializing experiment resources. @@ -102,7 +103,8 @@ void convertReinitStorage(const mxArray *prhs[], int nrhs, EchoframeResources &r * REINIT_INPUT_RFSTORE_POS, and RF storage is disabled without it. * @param res Reference to the native EchoframeResources structure to update. */ -void convertReinitExperiment(const mxArray *prhs[], int nrhs, EchoframeResources &res); +void convertReinitExperiment(const mxArray *prhs[], int nrhs, + EchoframeResources &res); /** * @brief Validates the fields and types of a MATLAB ReceiveSpec struct. @@ -176,6 +178,11 @@ void validateFourierReconStruct(const mxArray *s); void convertFourierReconSpecStructs( Beamform::FourierReconSpec &fourierReconSpec, const mxArray *reconStruct); +void validateDASReconStruct(const mxArray *s); + +void convertDASReconSpecStructs(Beamform::DASReconSpec &dasReconSpec, + const mxArray *reconStruct); + /** * @brief Validates the fields and types of a MATLAB StorageSpec struct. * @param s Pointer to the MATLAB struct mxArray. diff --git a/echoframe/cpp/src/python/echoframe_py_conversions.cpp b/echoframe/cpp/src/python/echoframe_py_conversions.cpp index 504a357..97699a8 100644 --- a/echoframe/cpp/src/python/echoframe_py_conversions.cpp +++ b/echoframe/cpp/src/python/echoframe_py_conversions.cpp @@ -19,6 +19,8 @@ #include #include +#include +#include #include #include "../../libs/Storage/src/storage/storage_spec.h" @@ -73,6 +75,16 @@ static ReceiveSpec to_receive(const py::dict &d) { return r; } +static BeamformerType to_beamformer_type(const py::dict &d) { + if (!d.contains("beamformerType")) return BeamformerType::Fourier; + std::string value = d["beamformerType"].cast(); + std::transform(value.begin(), value.end(), value.begin(), ::tolower); + if (value == "fourier") return BeamformerType::Fourier; + if (value == "das") return BeamformerType::DAS; + throw std::runtime_error( + "ReconSpec.beamformerType must be 'Fourier' or 'DAS'"); +} + static ReconSpec to_recon(const py::dict &d, const ReceiveSpec &rs) { ReconSpec r; r.nz = scalar_or_first(d["nz"], "nz"); @@ -103,6 +115,7 @@ static ReconSpec to_recon(const py::dict &d, const ReceiveSpec &rs) { r.getBF = scalar_or_first(d["getBF"], "getBF"); r.getPDI = scalar_or_first(d["getPDI"], "getPDI"); r.cropBF = scalar_or_first(d["cropBF"], "cropBF"); + r.beamformerType = to_beamformer_type(d); r.initialized = true; return r; } @@ -137,6 +150,37 @@ static Beamform::FourierReconSpec to_fourier_spec(const py::dict &src) { return f; } +static Beamform::DASReconSpec to_das_spec(const py::dict &src) { + Beamform::DASReconSpec d; + d.channelPositions = const_cast( + src["dasChannelPositions"].cast>().data()); + d.voxelPositions = const_cast( + src["dasVoxelPositions"].cast>().data()); + d.offsets = const_cast( + src["dasOffsets"].cast>().data()); + d.weights = const_cast( + src["dasWeights"].cast>().data()); + d.tgcVector = + const_cast(src["tgcVector"].cast>().data()); + d.wavenum = scalar_or_first(src["dasWavenum"], "dasWavenum"); + d.algorithm = + src.contains("dasAlgorithm") + ? scalar_or_first(src["dasAlgorithm"], "dasAlgorithm") + : 1; + d.computeType = + src.contains("dasComputeType") + ? scalar_or_first(src["dasComputeType"], "dasComputeType") + : 0; + if (src.contains("dasSourceDirections")) { + auto dirs = src["dasSourceDirections"].cast>(); + if (dirs.size() > 0) { + d.sourceDirections = const_cast(dirs.data()); + d.useDirectivity = true; + } + } + return d; +} + static PDI::PDISpec to_pdi(const py::dict &d, const ReconSpec &rec) { PDI::PDISpec p; p.ensemble_size = @@ -163,7 +207,11 @@ static EchoframeResources make_resources(const py::dict &recv, res.reconSpec = to_recon(recon, res.receiveSpec); res.pdiSpec = to_pdi(pdi, res.reconSpec); - res.fourierReconSpec = to_fourier_spec(recon); + if (res.reconSpec.beamformerType == BeamformerType::DAS) { + res.dasReconSpec = to_das_spec(recon); + } else { + res.fourierReconSpec = to_fourier_spec(recon); + } /* disable storage for offline */ res.bfStorageSpec.save = false; diff --git a/echoframe/cpp/src/python/echoframe_py_wrapper.cpp b/echoframe/cpp/src/python/echoframe_py_wrapper.cpp index b256d89..8733c9a 100644 --- a/echoframe/cpp/src/python/echoframe_py_wrapper.cpp +++ b/echoframe/cpp/src/python/echoframe_py_wrapper.cpp @@ -57,6 +57,18 @@ static void bind_specs(py::module_ &m) { RW(Beamform::ReceiveSpec, nActiveChannels) .def_readwrite("initialized", &Beamform::ReceiveSpec::initialized); + py::enum_(m, "BeamformerType") + .value("Fourier", Beamform::BeamformerType::Fourier) + .value("DAS", Beamform::BeamformerType::DAS); + + py::class_(m, "FourierReconSpec") + .def(py::init<>()); + py::class_(m, "DASReconSpec") + .def(py::init<>()) RW(Beamform::DASReconSpec, wavenum) + RW(Beamform::DASReconSpec, algorithm) + RW(Beamform::DASReconSpec, computeType) + RW(Beamform::DASReconSpec, useDirectivity); + /* ---- ReconSpec ------------------------------------------------------ */ py::class_(m, "ReconSpec") .def(py::init<>()) @@ -78,6 +90,7 @@ static void bind_specs(py::module_ &m) { /* behaviour flags */ RW(Beamform::ReconSpec, filterFrequencies) RW(Beamform::ReconSpec, getBF) RW(Beamform::ReconSpec, getPDI) + RW(Beamform::ReconSpec, beamformerType) .def_readwrite("initialized", &Beamform::ReconSpec::initialized); } @@ -156,7 +169,8 @@ PYBIND11_MODULE(echoframe, m) { .def_readwrite("rfTimeTagStorageSpec", &EchoframeResources::rfTimeTagStorageSpec) .def_readwrite("fourierReconSpec", - &EchoframeResources::fourierReconSpec); + &EchoframeResources::fourierReconSpec) + .def_readwrite("dasReconSpec", &EchoframeResources::dasReconSpec); py::class_(m, "EchoFrame") .def(py::init(), py::arg("resources"), diff --git a/echoframe/matlab/core/imaging/echoframe_setup_das.m b/echoframe/matlab/core/imaging/echoframe_setup_das.m new file mode 100644 index 0000000..7ed78e5 --- /dev/null +++ b/echoframe/matlab/core/imaging/echoframe_setup_das.m @@ -0,0 +1,59 @@ +function ReconSpec = echoframe_setup_das(ProbeSpec, TransmitSpec, ReceiveSpec, ReconSpec) +%ECHOFRAME_SETUP_DAS Build ffdas delay-and-sum geometry tables. +% EchoFrame IQ samples are depth-like: one complex sample advances the image +% depth by c/Fs. Geometric DAS sums transmit and receive paths, so scale +% positions by Fs/(2c) to map the two-way path onto EchoFrame's sample axis. + +if ~strcmpi(TransmitSpec.type, 'planewave') + error('echoframe_setup_das:unsupportedTransmit', ... + 'DAS setup currently supports TransmitSpec.type = ''planewave'' only.'); +end + +nz = double(ReconSpec.nz); +nx = double(ReconSpec.nx); +nTx = double(ReceiveSpec.nTransmissions); +fsOverTwoC = single(ReceiveSpec.Fs) / (2 * single(ReconSpec.c0)); +startSamples = single(ReconSpec.zAxis(1) * 1e-3) * 2 * fsOverTwoC; + +nActive = double(ProbeSpec.nElements); +channelX = ((0:nActive-1)' - (nActive-1)/2) * double(ProbeSpec.pitch); +ReconSpec.dasChannelPositions = single([channelX'; zeros(2, nActive)]) * fsOverTwoC; + +[xGridMm, zGridMm] = meshgrid(double(ReconSpec.xAxis), double(ReconSpec.zAxis)); +xGrid = xGridMm * 1e-3; +zGrid = zGridMm * 1e-3; +ReconSpec.dasVoxelPositions = single([xGrid(:)'; zeros(1, nz * nx); zGrid(:)']) * fsOverTwoC; + +receiveSamples = sqrt((xGrid(:) - channelX').^2 + zGrid(:).^2) * fsOverTwoC; +minOffset = reshape(1 - min(receiveSamples, [], 2), nz, nx); +maxOffset = reshape((double(ReceiveSpec.nSamplesIQ) - 2) - max(receiveSamples, [], 2), nz, nx); + +useTransmitDelays = isfield(TransmitSpec, 'transmitDelays') && any(TransmitSpec.transmitDelays(:) ~= 0); +txDelaySamples = double(TransmitSpec.transmitDelays) * double(ReceiveSpec.Fs) / 2; +if ndims(txDelaySamples) == 3 + txDelaySamples = squeeze(txDelaySamples(:, 1, :)); +end + +ReconSpec.dasOffsets = zeros(nz, nx, nTx, 'single'); +for iTx = 1:nTx + if useTransmitDelays + active = TransmitSpec.apodization(:) ~= 0; + txX = channelX(active)'; + txDelay = txDelaySamples(active, iTx)'; + txSamples = min(sqrt((xGrid(:) - txX).^2 + zGrid(:).^2) * fsOverTwoC + txDelay, [], 2); + offsets = reshape(txSamples, nz, nx) - startSamples; + else + theta = double(TransmitSpec.steer(iTx)); + offsets = (zGrid * cosd(theta) + xGrid * sind(theta)) * fsOverTwoC - startSamples; + end + ReconSpec.dasOffsets(:, :, iTx) = single(min(max(offsets, minOffset), maxOffset)); +end + +ReconSpec.dasWeights = ones(nz, nx, nTx, 'single') ./ single(max(1, nTx * nActive)); +% RFFormatter stores IQ as I - iQ, so ffdas' usual negative phase rotation +% is conjugated here. +ReconSpec.dasWavenum = single(4 * pi * double(ProbeSpec.Fc) / double(ReceiveSpec.Fs)); +ReconSpec.dasAlgorithm = int32(1); +ReconSpec.dasComputeType = int32(0); +ReconSpec.dasSourceDirections = single(repmat([0; 0; 1; cosd(35)], 1, nActive)); +end diff --git a/echoframe/matlab/core/imaging/echoframe_validate_structs.m b/echoframe/matlab/core/imaging/echoframe_validate_structs.m index f40a28e..737cfb8 100644 --- a/echoframe/matlab/core/imaging/echoframe_validate_structs.m +++ b/echoframe/matlab/core/imaging/echoframe_validate_structs.m @@ -71,14 +71,34 @@ 'extraVoxelsX', 'int32'; 'c0', 'single'; 'tgcVector', 'single'; - 'delayIndices', 'int32'; - 'interpolationWeights', 'single'; - 'frequencyAxis', 'single'; - 'planewaveDelays', 'single'; % Size = [TransmitSpec.nTransmissions,2] 'xAxis', 'double'; 'zAxis', 'double' }; +if ~isfield(ReconSpec, 'beamformerType') || isempty(ReconSpec.beamformerType) + ReconSpec.beamformerType = 'Fourier'; +end +expectedReconSpecFields(end+1, :) = {'beamformerType', 'char'}; + +if strcmpi(ReconSpec.beamformerType, 'DAS') + expectedReconSpecFields = [expectedReconSpecFields; { + 'dasChannelPositions', 'single'; + 'dasVoxelPositions', 'single'; + 'dasOffsets', 'single'; + 'dasWeights', 'single'; + 'dasWavenum', 'single'; + 'dasAlgorithm', 'int32'; + 'dasComputeType', 'int32' + }]; +else + expectedReconSpecFields = [expectedReconSpecFields; { + 'delayIndices', 'int32'; + 'interpolationWeights', 'single'; + 'frequencyAxis', 'single'; + 'planewaveDelays', 'single' % Size = [TransmitSpec.nTransmissions,2] + }]; +end + %% Expected fields for PDISpec expectedPDISpecFields = { 'ensembleSize', 'int32'; diff --git a/echoframe/matlab/core/imaging/initialize_image_reconstruction.m b/echoframe/matlab/core/imaging/initialize_image_reconstruction.m index f4ba459..58b4db8 100644 --- a/echoframe/matlab/core/imaging/initialize_image_reconstruction.m +++ b/echoframe/matlab/core/imaging/initialize_image_reconstruction.m @@ -67,10 +67,17 @@ tgc = ones(1,ReceiveSpec.nSamplesIQ - length(zeros_samples_start) - length(zeros_samples_end)); ReconSpec.tgcVector = convn([zeros_samples_start tgc zeros_samples_end],temp_kernel,'same'); -% Get the Fourier interpolation lookup tables -ReconSpec = echoframe_setup_fourier(ProbeSpec,TransmitSpec,ReceiveSpec,ReconSpec); +if ~isfield(ReconSpec, 'beamformerType') || isempty(ReconSpec.beamformerType) + ReconSpec.beamformerType = 'Fourier'; +end nElementRf = double(ReceiveSpec.nChannels); ReconSpec.xAxis = linspace(-(nElementRf/2)*ProbeSpec.pitch, (nElementRf/2)*ProbeSpec.pitch, Nx) * 1e3; ReconSpec.zAxis = linspace(ReceiveSpec.startDepthMm,ReceiveSpec.startDepthMm + ReceiveSpec.actualEndDepthMm, Nz); ReconSpec.imageSize = [Nz Nx]; + +if strcmpi(ReconSpec.beamformerType, 'DAS') + ReconSpec = echoframe_setup_das(ProbeSpec,TransmitSpec,ReceiveSpec,ReconSpec); +else + ReconSpec = echoframe_setup_fourier(ProbeSpec,TransmitSpec,ReceiveSpec,ReconSpec); +end diff --git a/echoframe/python/core/validate_structs_py_to_ef.py b/echoframe/python/core/validate_structs_py_to_ef.py index d962d10..5533d81 100644 --- a/echoframe/python/core/validate_structs_py_to_ef.py +++ b/echoframe/python/core/validate_structs_py_to_ef.py @@ -51,12 +51,9 @@ ("extraVoxelsX", INT32), ("c0", FLOAT), ("tgcVector", FLOAT), - ("delayIndices", INT32), - ("interpolationWeights",CFLOAT), - ("frequencyAxis", FLOAT), - ("planewaveDelays", FLOAT), ("xAxis", DOUBLE), ("zAxis", DOUBLE), + ("beamformerType", str), ], "PDISpec": [ ("ensembleSize", INT32), @@ -86,7 +83,7 @@ def _cast(value, dtype): def _validate_and_cast(name: str, spec: dict) -> dict: """Check mandatory fields & cast each to the expected dtype.""" - out = {} + out = dict(spec) for field, dtype in _EXPECTED[name]: if field not in spec: raise ValueError(f"{name}: missing field '{field}'") @@ -107,6 +104,7 @@ def validate_specs( PDISpec: dict, ): ReceiveSpec["nElements"] = ProbeSpec["nElements"] + ReconSpec.setdefault("beamformerType", "Fourier") """Replicates MATLAB echoframe_validate_structs in Python.""" @@ -117,8 +115,22 @@ def validate_specs( PDISpec = _validate_and_cast("PDISpec", PDISpec) ReconSpec["croppingROI"] = ReconSpec["croppingROI"].astype(np.int32).ravel(order='C') - ReconSpec["interpolationWeights"] = ReconSpec["interpolationWeights"].ravel(order='C') - ReconSpec["delayIndices"] = ReconSpec["delayIndices"].ravel(order='C') + if ReconSpec["beamformerType"].lower() == "das": + for key in ("dasChannelPositions", "dasVoxelPositions", "dasOffsets", "dasWeights"): + if key not in ReconSpec: + raise ValueError(f"ReconSpec: missing field '{key}'") + ReconSpec[key] = np.asarray(ReconSpec[key], dtype=FLOAT).ravel(order='C') + ReconSpec["dasWavenum"] = FLOAT(ReconSpec["dasWavenum"]) + ReconSpec["dasAlgorithm"] = INT32(ReconSpec.get("dasAlgorithm", 1)) + ReconSpec["dasComputeType"] = INT32(ReconSpec.get("dasComputeType", 0)) + if "dasSourceDirections" in ReconSpec: + ReconSpec["dasSourceDirections"] = np.asarray(ReconSpec["dasSourceDirections"], dtype=FLOAT).ravel(order='C') + else: + for key, dtype in (("interpolationWeights", CFLOAT), ("delayIndices", INT32), + ("frequencyAxis", FLOAT), ("planewaveDelays", FLOAT)): + if key not in ReconSpec: + raise ValueError(f"ReconSpec: missing field '{key}'") + ReconSpec[key] = np.asarray(ReconSpec[key], dtype=dtype).ravel(order='C') ReconSpec["xAxis"] = ReconSpec["xAxis"].ravel(order='C') ReconSpec["zAxis"] = ReconSpec["zAxis"].ravel(order='C') diff --git a/echoframe/python/examples/compare_beamformers.py b/echoframe/python/examples/compare_beamformers.py new file mode 100644 index 0000000..8f4a57f --- /dev/null +++ b/echoframe/python/examples/compare_beamformers.py @@ -0,0 +1,240 @@ +""" +Compare EchoFrame Fourier and ffdas DAS beamformers on one stored RF frame. + +Requires an EchoFrame Python build configured with ``-DEF_USE_FFDAS=ON``. +Edit ``load_path`` below or pass it on the command line. +""" + +from __future__ import annotations + +import argparse +import copy +import sys +import time +from pathlib import Path + +import h5py + +# import matplotlib.pyplot as plt +import numpy as np + +import echoframe as ef + +_HERE = Path(__file__).resolve().parent +sys.path.insert(1, str((_HERE / "../core").resolve())) +import rf_io # noqa: E402 +from validate_structs_py_to_ef import validate_specs # noqa: E402 + +load_path = Path("") + + +def h5struct_to_dict(group: h5py.Group) -> dict: + """Convert a MATLAB v7.3 struct group to a dict.""" + out = {} + for k, item in group.items(): + name = k.rstrip("\x00") + if isinstance(item, h5py.Dataset): + data = item[()] + if isinstance(data, np.ndarray) and data.dtype.names == ("real", "imag"): + data = data["real"] + 1j * data["imag"] + elif item.attrs.get("MATLAB_complex", None) == np.bytes_("1"): + real, imag = np.split(data, 2, axis=-1) + data = real.squeeze(-1) + 1j * imag.squeeze(-1) + elif data.dtype.kind in ("u", "i") and data.dtype.itemsize == 2: + data = "".join(chr(x) for x in data.squeeze()) + if isinstance(data, np.ndarray) and data.shape == (): + data = data.item() + out[name] = data + elif isinstance(item, h5py.Group): + out[name] = h5struct_to_dict(item) + return out + + +def normalize_specs(probe: dict, transmit: dict, receive: dict, recon: dict, pdi: dict) -> tuple[dict, dict, dict, dict, dict]: + """Map Effusive field names onto EchoFrame Python example names.""" + probe = dict(probe) + transmit = dict(transmit) + receive = dict(receive) + recon = dict(recon) + pdi = dict(pdi) + + if "pitch" not in probe and "pitchX" in probe: + probe["pitch"] = probe["pitchX"] + if "nElements" not in probe and "nElementsX" in probe: + probe["nElements"] = probe["nElementsX"] + if "elementPosition" not in probe and "element_position" in probe: + probe["elementPosition"] = probe["element_position"] + + if "steer" not in transmit and "steerX" in transmit: + transmit["steer"] = transmit["steerX"] + + if "samplingMode" not in receive and "sampling_mode" in receive: + receive["samplingMode"] = receive["sampling_mode"] + if "samplesPerWavelength" not in receive and "samples_per_wavelength" in receive: + receive["samplesPerWavelength"] = receive["samples_per_wavelength"] + + if "extraVoxelsZ" not in recon and "extra_voxels_z" in recon: + recon["extraVoxelsZ"] = recon["extra_voxels_z"] + if "extraVoxelsX" not in recon and "extra_voxels_x" in recon: + recon["extraVoxelsX"] = recon["extra_voxels_x"] + + return probe, transmit, receive, recon, pdi + + +def find_session_files(path: Path) -> tuple[Path, Path]: + """Accept either example folder or an Effusive *_rf.dat file.""" + if path.is_file(): + if path.name.endswith("_rf.dat"): + return path.with_name(path.name.replace("_rf.dat", "_seq.mat")), path + return path, path.with_name(path.name.replace("_seq.mat", "_rf.dat")) + return path / "ScanParameters.mat", path / "rf_acq.dat" + + +def add_das_spec(probe: dict, transmit: dict, receive: dict, recon: dict) -> dict: + """Return a ReconSpec copy with the DAS fields EchoFrame expects.""" + out = copy.deepcopy(recon) + out["beamformerType"] = "DAS" + + nz = int(np.asarray(out["nz"]).ravel()[0]) + nx = int(np.asarray(out["nx"]).ravel()[0]) + ntx = int(np.asarray(receive["nTransmissions"]).ravel()[0]) + nelem = int(np.asarray(probe["nElements"]).ravel()[0]) + pitch = float(np.asarray(probe["pitch"]).ravel()[0]) + fs = float(np.asarray(receive["Fs"]).ravel()[0]) + c0 = float(np.asarray(out["c0"]).ravel()[0]) + # EchoFrame IQ samples are depth-like: one complex sample advances the image + # depth by c/Fs. Geometric DAS sums transmit and receive paths, so scale + # positions by Fs/(2c) to map the two-way path onto EchoFrame's sample axis. + fs_over_two_c = fs / (2 * c0) + + channel_x = (np.arange(nelem, dtype=np.float32) - (nelem - 1) / 2) * pitch + out["dasChannelPositions"] = ( + np.column_stack([channel_x, np.zeros((nelem, 2), dtype=np.float32)]) + * fs_over_two_c + ).astype(np.float32).ravel() + + x_axis = np.asarray(out["xAxis"], dtype=np.float64).ravel() * 1e-3 + z_axis = np.asarray(out["zAxis"], dtype=np.float64).ravel() * 1e-3 + start_samples = z_axis[0] * 2 * fs_over_two_c + x_grid, z_grid = np.meshgrid(x_axis, z_axis) + voxel_x = x_grid.ravel(order="F") + voxel_z = z_grid.ravel(order="F") + out["dasVoxelPositions"] = ( + np.column_stack([voxel_x, np.zeros(nz * nx), voxel_z]) * fs_over_two_c + ).astype(np.float32).ravel() + + receive_samples = ( + np.sqrt((voxel_x[:, None] - channel_x[None, :]) ** 2 + voxel_z[:, None] ** 2) + * fs_over_two_c + ) + min_offset = 1 - receive_samples.min(axis=1) + max_offset = (int(np.asarray(receive["nSamplesIQ"]).ravel()[0]) - 2) - receive_samples.max(axis=1) + + steer = np.asarray(transmit["steer"], dtype=np.float64).ravel() + transmit_delays = np.asarray(transmit.get("transmitDelays", 0), dtype=np.float64) + use_transmit_delays = transmit_delays.size > 1 and np.any(transmit_delays != 0) + if use_transmit_delays: + transmit_delays = np.squeeze(transmit_delays) + if transmit_delays.shape[0] == ntx and transmit_delays.shape[-1] == nelem: + transmit_delays = transmit_delays.T + else: + transmit_delays = transmit_delays.reshape(nelem, -1) + transmit_delays = transmit_delays[:, :ntx] * fs / 2 + apod = np.asarray(transmit.get("apodization", np.ones(nelem)), dtype=bool).ravel() + tx_x = channel_x[apod] + + offsets = np.empty((ntx, nz * nx), dtype=np.float32) + for i, theta in enumerate(steer[:ntx]): + if use_transmit_delays: + tx_samples = ( + np.sqrt((voxel_x[:, None] - tx_x[None, :]) ** 2 + voxel_z[:, None] ** 2) + * fs_over_two_c + + transmit_delays[apod, i][None, :] + ).min(axis=1) + raw_offsets = tx_samples - start_samples + else: + raw_offsets = ( + ( + z_grid * np.cos(np.deg2rad(theta)) + + x_grid * np.sin(np.deg2rad(theta)) + ).ravel(order="F") + * fs_over_two_c + - start_samples + ) + offsets[i] = np.clip(raw_offsets, min_offset, max_offset) + out["dasOffsets"] = offsets.ravel() + out["dasWeights"] = np.full( + ntx * nz * nx, 1.0 / max(1, ntx * nelem), dtype=np.float32 + ) + # RFFormatter stores IQ as I - iQ, so ffdas' usual negative phase rotation + # is conjugated here. + out["dasWavenum"] = np.float32( + 4 * np.pi * float(np.asarray(probe["Fc"]).ravel()[0]) / fs + ) + out["dasAlgorithm"] = np.int32(1) + out["dasComputeType"] = np.int32(0) + out["dasSourceDirections"] = np.tile( + np.array([0, 0, 1, np.cos(np.deg2rad(35))], dtype=np.float32), nelem + ) + return out + + +def run_once(receive: dict, recon: dict, pdi: dict, rf: np.ndarray) -> tuple[np.ndarray, float]: + """Process one RF buffer and return B-mode plus elapsed seconds.""" + res = ef.make_resources(receive, recon, pdi) + core = ef.EchoFrame(res, use_storage=False) + t0 = time.perf_counter() + _, bmode, _ = core.process(rf, start_storage=False) + return bmode, time.perf_counter() - t0 + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("path", nargs="?", type=Path, default=load_path) + parser.add_argument("--no-show", action="store_true") + args = parser.parse_args() + if not args.path or not args.path.exists(): + raise SystemExit("Pass a folder containing ScanParameters.mat/rf_acq.dat or an Effusive *_rf.dat file.") + + seq_path, rf_path = find_session_files(args.path) + with h5py.File(seq_path, "r") as f: + probe = h5struct_to_dict(f["ProbeSpec"]) + transmit = h5struct_to_dict(f["TransmitSpec"]) + receive = h5struct_to_dict(f["ReceiveSpec"]) + recon = h5struct_to_dict(f["ReconSpec"]) + pdi = h5struct_to_dict(f["PDISpec"]) + probe, transmit, receive, recon, pdi = normalize_specs(probe, transmit, receive, recon, pdi) + + _, _, receive_f, recon_f, pdi_f = validate_specs(probe, transmit, receive, recon, pdi) + rf = rf_io.read_stored_rf(rf_path, receive_f).ravel(order="F") + fourier_bmode, fourier_s = run_once(receive_f, recon_f, pdi_f, rf) + + recon_d = add_das_spec(probe, transmit, receive, recon) + _, _, receive_d, recon_d, pdi_d = validate_specs(probe, transmit, receive, recon_d, pdi) + das_bmode, das_s = run_once(receive_d, recon_d, pdi_d, rf) + + print(f"Fourier: {fourier_s:.3f} s") + print(f"DAS: {das_s:.3f} s") + print(f"Mean abs B-mode diff: {np.mean(np.abs(fourier_bmode - das_bmode)):.6g}") + + if args.no_show: + return + vmax = max(float(fourier_bmode.max()), float(das_bmode.max()), np.finfo(np.float32).eps) + imgs = [fourier_bmode, das_bmode] + titles = ["Fourier", "DAS"] + # fig, axes = plt.subplots(1, 2, figsize=(8, 4), constrained_layout=True, dpi=300) + # for ax, img, title in zip(axes, imgs, titles): + # im = ax.imshow( + # 20 * np.log10(img.T / vmax + np.finfo(np.float32).eps), + # cmap="gray", + # vmin=-60, + # vmax=0, + # ) + # ax.set_title(title) + # ax.axis("image") + # fig.colorbar(im, ax=axes, label="dB", shrink=0.8) + # plt.show() + + +if __name__ == "__main__": + main() diff --git a/echoframe/python/examples/generate_synthetic_session.py b/echoframe/python/examples/generate_synthetic_session.py new file mode 100644 index 0000000..9e0383f --- /dev/null +++ b/echoframe/python/examples/generate_synthetic_session.py @@ -0,0 +1,166 @@ +""" +Generate a tiny synthetic EchoFrame session for offline examples. + +Writes ``ScanParameters.mat`` and ``rf_acq.dat`` with the minimum fields needed by +``compare_beamformers.py``. The data are deterministic and intentionally small. +""" + +from __future__ import annotations + +import argparse +import struct +from pathlib import Path + +import h5py +import numpy as np + + +def _char(s: str) -> np.ndarray: + return np.array([ord(c) for c in s], dtype=np.uint16) + + +def _dataset(group: h5py.Group, name: str, value) -> None: + group.create_dataset(name, data=value) + + +def _write_struct(root: h5py.File, name: str, fields: dict) -> None: + group = root.create_group(name) + for key, value in fields.items(): + _dataset(group, key, value) + + +def make_specs() -> tuple[dict, dict, dict, dict, dict]: + """Build deterministic small EchoFrame structs.""" + n_iq = 64 + n_samples = 2 * n_iq + n_channels = 16 + n_tx = 3 + n_repeats = 8 + pitch = 300e-6 + fs = 20e6 + fc = 5e6 + c0 = 1540.0 + + # DAS delay is transmit distance + receive distance in samples. Keep this + # tiny grid inside the 64 IQ samples so the smoke test exercises in-bounds + # interpolation instead of mostly returning zeros. + z_axis = np.linspace(0.2, 2.2, n_iq, dtype=np.float64) + x_axis = np.linspace(-(n_channels / 2) * pitch, (n_channels / 2) * pitch, n_channels) * 1e3 + + z_idx = np.arange(n_iq, dtype=np.int32)[:, None, None] + delay_indices = np.broadcast_to(z_idx, (n_iq, n_channels, n_tx)).copy() + interpolation_weights = np.ones( + (n_iq, n_channels, n_tx), dtype=np.complex64 + ) / np.float32(n_tx * n_channels * n_iq * n_channels * 2) + + probe = { + "pitch": np.array(pitch, dtype=np.float64), + "Fc": np.array(fc, dtype=np.float32), + "nElements": np.array(n_channels, dtype=np.int32), + "elementPosition": np.column_stack([ + (np.arange(n_channels) - (n_channels - 1) / 2) * pitch, + np.zeros((n_channels, 4)), + ]).astype(np.float32), + } + transmit = { + "c0": np.array(c0, dtype=np.float64), + "type": _char("planewave"), + "steer": np.array([-8.0, 0.0, 8.0], dtype=np.float32), + "apodization": np.ones(n_channels, dtype=np.float64), + "transmitDelays": np.zeros((n_channels, 1, n_tx), dtype=np.float64), + } + receive = { + "nSamples": np.array(n_samples, dtype=np.int32), + "nSamplesIQ": np.array(n_iq, dtype=np.int32), + "nTransmissions": np.array(n_tx, dtype=np.int32), + "nRepeats": np.array(n_repeats, dtype=np.int32), + "nChannels": np.array(n_channels, dtype=np.int32), + "channel2ElementMap": np.arange(n_channels, dtype=np.int32), + "nElements": np.array(n_channels, dtype=np.int32), + "Fs": np.array(fs, dtype=np.float32), + "samplingMode": _char("NS200BW"), + "samplesPerWavelength": np.array(4, dtype=np.int32), + } + recon = { + "bfDataType": _char("complex single"), + "getBF": np.array(True, dtype=np.bool_), + "getPDI": np.array(True, dtype=np.bool_), + "nz": np.array(n_iq, dtype=np.int32), + "nx": np.array(n_channels, dtype=np.int32), + "filterFrequencies": np.array(False, dtype=np.bool_), + "cropBF": np.array(False, dtype=np.bool_), + "croppingROI": np.array([0, n_iq - 1, 0, n_channels - 1], dtype=np.int32), + "extraVoxelsZ": np.array(0, dtype=np.int32), + "extraVoxelsX": np.array(0, dtype=np.int32), + "c0": np.array(c0, dtype=np.float32), + "tgcVector": np.ones(n_iq, dtype=np.float32), + "delayIndices": delay_indices, + "interpolationWeights": interpolation_weights, + "frequencyAxis": np.linspace(-fs / 4, fs / 4, n_iq, dtype=np.float32), + "planewaveDelays": np.zeros((n_tx, 2), dtype=np.float32), + "xAxis": x_axis.astype(np.float64), + "zAxis": z_axis, + } + pdi = { + "ensembleSize": np.array(n_repeats, dtype=np.int32), + "threshold": np.array(0.4, dtype=np.float32), + "shiftSize": np.array(n_repeats, dtype=np.int32), + "cropPDI": np.array(False, dtype=np.bool_), + "svdMethod": _char("Covariance"), + } + return probe, transmit, receive, recon, pdi + + +def make_rf(receive: dict) -> np.ndarray: + """Create interleaved int16 IQ RF data as EchoFrame stores it.""" + n_iq = int(receive["nSamplesIQ"]) + n_samples = int(receive["nSamples"]) + n_tx = int(receive["nTransmissions"]) + n_repeats = int(receive["nRepeats"]) + n_channels = int(receive["nChannels"]) + rows_iq = n_iq * n_tx * n_repeats + rf = np.zeros((n_samples * n_tx * n_repeats, n_channels), dtype=np.int16, order="F") + + t = np.arange(n_iq, dtype=np.float32) + for ch in range(n_channels): + center = 12 + ch * 2.2 + echo = 9000 * np.exp(-0.5 * ((t - center) / 2.5) ** 2) + phase = ch * 0.25 + iq_real = echo * np.cos(phase) + iq_imag = echo * np.sin(phase) + for tx in range(n_tx): + for rep in range(n_repeats): + base = (rep * n_tx + tx) * n_iq + amp = 1.0 + 0.02 * rep + 0.05 * tx + rf[2 * (base + np.arange(n_iq)), ch] = np.round(amp * iq_real).astype(np.int16) + rf[2 * (base + np.arange(n_iq)) + 1, ch] = np.round(amp * iq_imag).astype(np.int16) + return rf + + +def write_rf(path: Path, rf: np.ndarray) -> None: + """Write EchoFrame RF storage header plus one RF buffer.""" + payload = rf.ravel(order="F") + header = struct.pack("<6Q", 1, 48, 1, payload.nbytes, 0, 0) + path.write_bytes(header + payload.tobytes()) + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("out", type=Path, help="Output session directory") + args = parser.parse_args() + args.out.mkdir(parents=True, exist_ok=True) + + probe, transmit, receive, recon, pdi = make_specs() + with h5py.File(args.out / "ScanParameters.mat", "w") as f: + _write_struct(f, "ProbeSpec", probe) + _write_struct(f, "TransmitSpec", transmit) + _write_struct(f, "ReceiveSpec", receive) + _write_struct(f, "ReconSpec", recon) + _write_struct(f, "PDISpec", pdi) + write_rf(args.out / "rf_acq.dat", make_rf(receive)) + print(f"Wrote {args.out / 'ScanParameters.mat'}") + print(f"Wrote {args.out / 'rf_acq.dat'}") + + +if __name__ == "__main__": + main() diff --git a/justfile b/justfile index 967fde5..b7f32db 100644 --- a/justfile +++ b/justfile @@ -9,6 +9,7 @@ matlab_root := if os() == "windows" { env_var_or_default("Matlab_ROOT_DIR", "C:/ cuda_host_compiler := if os() == "windows" { "" } else { env_var_or_default("CMAKE_CUDA_HOST_COMPILER", "/usr/bin/g++") } default_cuda_architectures := if os() == "windows" { "" } else { `if nvcc --version 2>/dev/null | grep -Eq 'release 1[3-9]\.'; then printf '75;80;86;89;90;100;103;121'; else printf '61;75;86;89;90'; fi` } cuda_architectures := env_var_or_default("CMAKE_CUDA_ARCHITECTURES", default_cuda_architectures) +use_ffdas := env_var_or_default("EF_USE_FFDAS", "OFF") vcpkg_root := env_var_or_default("VCPKG_ROOT", "C:/vcpkg") vcpkg_triplet := "x64-windows-static" test_img := env_var_or_default("HOME", "") / "ef_test.img" @@ -34,6 +35,7 @@ configure: -DEF_BUILD_CLI=ON \ -DEF_BUILD_MEX=ON \ -DEF_BUILD_PYTHON=ON \ + -DEF_USE_FFDAS={{use_ffdas}} \ -DMatlab_ROOT_DIR="{{matlab_root}}" [windows] @@ -42,6 +44,7 @@ configure: # with $VCPKG_ROOT if it's not at C:/vcpkg. cmake -S {{source_dir}} -B {{build_dir}} -G "Visual Studio 17 2022" -A x64 \ -DCMAKE_BUILD_TYPE=Release \ + -DEF_USE_FFDAS={{use_ffdas}} \ -DCMAKE_TOOLCHAIN_FILE="{{vcpkg_root}}/scripts/buildsystems/vcpkg.cmake" \ -DCMAKE_PREFIX_PATH="{{vcpkg_root}}/installed/{{vcpkg_triplet}}" From adeeb876ae340b85d94046835a2458c09f6c99c0 Mon Sep 17 00:00:00 2001 From: Samuel Le Meur-Diebolt Date: Wed, 30 Sep 2026 13:35:44 +0100 Subject: [PATCH 03/11] fix(das): planewave TX path with Verasonics bulk delay Rebuild DAS offsets from steer angles plus the per-tx bulk delay that Verasonics stores by clamping negative TX.Delay to zero. Omitting that bulk desynchronizes steered txs and blurs mid-depth. Drop the in-bounds max_offset clip so deep voxels keep correct depth with partial aperture (ffdas zeros OOB reads). Fix zAxis to treat actualEndDepthMm as absolute. Point the ffdas submodule at the OOB IQ-load guard on sdiebolt/ffdas:fix/das-alg1-oob-iq-load. --- echoframe/cpp/libs/ffdas | 2 +- .../matlab/core/imaging/echoframe_setup_das.m | 36 ++++++---- .../imaging/initialize_image_reconstruction.m | 3 +- .../python/examples/compare_beamformers.py | 67 ++++++++++--------- 4 files changed, 62 insertions(+), 46 deletions(-) diff --git a/echoframe/cpp/libs/ffdas b/echoframe/cpp/libs/ffdas index 6466de3..ea3a25e 160000 --- a/echoframe/cpp/libs/ffdas +++ b/echoframe/cpp/libs/ffdas @@ -1 +1 @@ -Subproject commit 6466de37dc3cb2f18946a2ef4543b788af783e95 +Subproject commit ea3a25e4d7d8c330fa8c8ee61cdb9a52143fd87f diff --git a/echoframe/matlab/core/imaging/echoframe_setup_das.m b/echoframe/matlab/core/imaging/echoframe_setup_das.m index 7ed78e5..d428b23 100644 --- a/echoframe/matlab/core/imaging/echoframe_setup_das.m +++ b/echoframe/matlab/core/imaging/echoframe_setup_das.m @@ -26,27 +26,39 @@ receiveSamples = sqrt((xGrid(:) - channelX').^2 + zGrid(:).^2) * fsOverTwoC; minOffset = reshape(1 - min(receiveSamples, [], 2), nz, nx); -maxOffset = reshape((double(ReceiveSpec.nSamplesIQ) - 2) - max(receiveSamples, [], 2), nz, nx); +% Planewave TX path is z*cos(theta)+x*sin(theta) plus the per-tx bulk delay +% Verasonics stores by clamping negative TX.Delay to zero. That bulk is in +% the RF; omitting it desynchronizes steered txs and blurs mid-depth. Do not +% rebuild the wavefront via min-over-elements or clip maxOffset: the +% geometric path plus bulk matches min-over(active) to <0.01 sample, and +% ffdas already zeros channel reads past the IQ buffer so deep rows keep +% correct depth with partial aperture. useTransmitDelays = isfield(TransmitSpec, 'transmitDelays') && any(TransmitSpec.transmitDelays(:) ~= 0); -txDelaySamples = double(TransmitSpec.transmitDelays) * double(ReceiveSpec.Fs) / 2; -if ndims(txDelaySamples) == 3 - txDelaySamples = squeeze(txDelaySamples(:, 1, :)); +if useTransmitDelays + txDelays = double(TransmitSpec.transmitDelays); + if ndims(txDelays) == 3 + txDelays = squeeze(txDelays(:, 1, :)); + end + if size(txDelays, 1) == nTx + txDelays = txDelays.'; + end + active = double(TransmitSpec.apodization(:)) ~= 0; end ReconSpec.dasOffsets = zeros(nz, nx, nTx, 'single'); for iTx = 1:nTx + theta = double(TransmitSpec.steer(iTx)); + sinTheta = sind(theta); if useTransmitDelays - active = TransmitSpec.apodization(:) ~= 0; - txX = channelX(active)'; - txDelay = txDelaySamples(active, iTx)'; - txSamples = min(sqrt((xGrid(:) - txX).^2 + zGrid(:).^2) * fsOverTwoC + txDelay, [], 2); - offsets = reshape(txSamples, nz, nx) - startSamples; + residual = txDelays(:, iTx) - channelX * sinTheta / double(ReconSpec.c0); + bulk = median(residual(active)); else - theta = double(TransmitSpec.steer(iTx)); - offsets = (zGrid * cosd(theta) + xGrid * sind(theta)) * fsOverTwoC - startSamples; + bulk = 0; end - ReconSpec.dasOffsets(:, :, iTx) = single(min(max(offsets, minOffset), maxOffset)); + offsets = (zGrid * cosd(theta) + xGrid * sinTheta) * fsOverTwoC ... + + bulk * double(ReceiveSpec.Fs) / 2 - startSamples; + ReconSpec.dasOffsets(:, :, iTx) = single(max(offsets, minOffset)); end ReconSpec.dasWeights = ones(nz, nx, nTx, 'single') ./ single(max(1, nTx * nActive)); diff --git a/echoframe/matlab/core/imaging/initialize_image_reconstruction.m b/echoframe/matlab/core/imaging/initialize_image_reconstruction.m index 58b4db8..9fab3f5 100644 --- a/echoframe/matlab/core/imaging/initialize_image_reconstruction.m +++ b/echoframe/matlab/core/imaging/initialize_image_reconstruction.m @@ -73,7 +73,8 @@ nElementRf = double(ReceiveSpec.nChannels); ReconSpec.xAxis = linspace(-(nElementRf/2)*ProbeSpec.pitch, (nElementRf/2)*ProbeSpec.pitch, Nx) * 1e3; -ReconSpec.zAxis = linspace(ReceiveSpec.startDepthMm,ReceiveSpec.startDepthMm + ReceiveSpec.actualEndDepthMm, Nz); +% actualEndDepthMm is absolute (includes startDepthMm), not a span. +ReconSpec.zAxis = linspace(ReceiveSpec.startDepthMm, ReceiveSpec.actualEndDepthMm, Nz); ReconSpec.imageSize = [Nz Nx]; if strcmpi(ReconSpec.beamformerType, 'DAS') diff --git a/echoframe/python/examples/compare_beamformers.py b/echoframe/python/examples/compare_beamformers.py index 8f4a57f..02c9b02 100644 --- a/echoframe/python/examples/compare_beamformers.py +++ b/echoframe/python/examples/compare_beamformers.py @@ -14,8 +14,7 @@ from pathlib import Path import h5py - -# import matplotlib.pyplot as plt +import matplotlib.pyplot as plt import numpy as np import echoframe as ef @@ -128,8 +127,14 @@ def add_das_spec(probe: dict, transmit: dict, receive: dict, recon: dict) -> dic * fs_over_two_c ) min_offset = 1 - receive_samples.min(axis=1) - max_offset = (int(np.asarray(receive["nSamplesIQ"]).ravel()[0]) - 2) - receive_samples.max(axis=1) + # Planewave TX path is z*cos(theta)+x*sin(theta) plus the per-tx bulk delay + # Verasonics stores by clamping negative TX.Delay to zero. That bulk is in + # the RF; omitting it desynchronizes steered txs and blurs mid-depth. Do not + # rebuild the wavefront via min-over-elements or clip max_offset: the + # geometric path plus bulk matches min-over(active) to <0.01 sample, and + # ffdas already zeros channel reads past the IQ buffer so deep rows keep + # correct depth with partial aperture. steer = np.asarray(transmit["steer"], dtype=np.float64).ravel() transmit_delays = np.asarray(transmit.get("transmitDelays", 0), dtype=np.float64) use_transmit_delays = transmit_delays.size > 1 and np.any(transmit_delays != 0) @@ -139,29 +144,27 @@ def add_das_spec(probe: dict, transmit: dict, receive: dict, recon: dict) -> dic transmit_delays = transmit_delays.T else: transmit_delays = transmit_delays.reshape(nelem, -1) - transmit_delays = transmit_delays[:, :ntx] * fs / 2 + transmit_delays = transmit_delays[:, :ntx] apod = np.asarray(transmit.get("apodization", np.ones(nelem)), dtype=bool).ravel() - tx_x = channel_x[apod] offsets = np.empty((ntx, nz * nx), dtype=np.float32) for i, theta in enumerate(steer[:ntx]): + sin_t = np.sin(np.deg2rad(theta)) if use_transmit_delays: - tx_samples = ( - np.sqrt((voxel_x[:, None] - tx_x[None, :]) ** 2 + voxel_z[:, None] ** 2) - * fs_over_two_c - + transmit_delays[apod, i][None, :] - ).min(axis=1) - raw_offsets = tx_samples - start_samples + residual = transmit_delays[:, i] - channel_x * sin_t / c0 + bulk = float(np.median(residual[apod])) if np.any(apod) else 0.0 else: - raw_offsets = ( - ( - z_grid * np.cos(np.deg2rad(theta)) - + x_grid * np.sin(np.deg2rad(theta)) - ).ravel(order="F") - * fs_over_two_c - - start_samples - ) - offsets[i] = np.clip(raw_offsets, min_offset, max_offset) + bulk = 0.0 + raw_offsets = ( + ( + z_grid * np.cos(np.deg2rad(theta)) + + x_grid * sin_t + ).ravel(order="F") + * fs_over_two_c + + bulk * fs / 2 + - start_samples + ) + offsets[i] = np.maximum(raw_offsets, min_offset) out["dasOffsets"] = offsets.ravel() out["dasWeights"] = np.full( ntx * nz * nx, 1.0 / max(1, ntx * nelem), dtype=np.float32 @@ -222,18 +225,18 @@ def main() -> None: vmax = max(float(fourier_bmode.max()), float(das_bmode.max()), np.finfo(np.float32).eps) imgs = [fourier_bmode, das_bmode] titles = ["Fourier", "DAS"] - # fig, axes = plt.subplots(1, 2, figsize=(8, 4), constrained_layout=True, dpi=300) - # for ax, img, title in zip(axes, imgs, titles): - # im = ax.imshow( - # 20 * np.log10(img.T / vmax + np.finfo(np.float32).eps), - # cmap="gray", - # vmin=-60, - # vmax=0, - # ) - # ax.set_title(title) - # ax.axis("image") - # fig.colorbar(im, ax=axes, label="dB", shrink=0.8) - # plt.show() + fig, axes = plt.subplots(1, 2, figsize=(8, 4), constrained_layout=True, dpi=300) + for ax, img, title in zip(axes, imgs, titles): + im = ax.imshow( + 20 * np.log10(img.T / vmax + np.finfo(np.float32).eps), + cmap="gray", + vmin=-60, + vmax=0, + ) + ax.set_title(title) + ax.axis("image") + fig.colorbar(im, ax=axes, label="dB", shrink=0.8) + plt.show() if __name__ == "__main__": From 9693e8ae144ad7dd089478725c5f6b36e6c7f235 Mon Sep 17 00:00:00 2001 From: Samuel Le Meur-Diebolt Date: Wed, 30 Sep 2026 15:58:26 +0100 Subject: [PATCH 04/11] build: always build the ffdas DAS beamformer Drop the EF_USE_FFDAS toggle. The justfile read it from the environment and defaulted it to OFF, and build_echoframe.sh never passed it at all, so the DAS backend was opt-in and easy to lose: a build could succeed and then throw "DAS beamforming requires configuring with -DEF_USE_FFDAS=ON" at runtime, long after the mistake was made. Remove the option() call, the EF_USE_FFDAS compile definition, and the conditional compilation in the FFDAS beamformer, so ffdas is resolved, linked, and compiled unconditionally. A missing ffdas submodule is now a configure-time error naming the git submodule command, rather than a silently reduced build. Document the requirement in both READMEs: the ffdas submodule must be initialised before configuring, and building it from source needs CMake 3.26 or newer. --- README.md | 6 +- echoframe/cpp/libs/README.md | 7 +- echoframe/cpp/src/CMakeLists.txt | 74 +++++++++---------- .../cpp/src/beamformer/das/ffdas_beamformer.h | 8 +- .../src/beamformer/das/ffdas_beamformer.t.hpp | 20 ----- .../python/examples/compare_beamformers.py | 3 +- justfile | 3 - 7 files changed, 47 insertions(+), 74 deletions(-) diff --git a/README.md b/README.md index ae60047..cc1aab0 100644 --- a/README.md +++ b/README.md @@ -183,13 +183,17 @@ or: ## Quick Start -### 🔧 Submodules +### 🔧 Submodules Before continuing, make sure that the submodules included in the project are properly initialized. ```bash # make sure you are in the EchoFrame directory git submodule init git submodule update ``` +The [`ffdas`](echoframe/cpp/libs/README.md) submodule is **required**: it provides the +delay-and-sum (DAS) beamformer, which is built by default. If it is missing, CMake +stops with an error telling you to run the command above. Only GSL and ffdas are +needed; the rest of the build has no submodule dependencies. You can continue with the rest of the build instructions if you have no errors. Below you can find a set of commands to setup the project for both Linux and Windows. diff --git a/echoframe/cpp/libs/README.md b/echoframe/cpp/libs/README.md index 7302e24..89cf631 100644 --- a/echoframe/cpp/libs/README.md +++ b/echoframe/cpp/libs/README.md @@ -13,9 +13,10 @@ Library. ## ffdas -A submodule (`github.com/BrainEchoLab/ffdas`) providing the optional CUDA -DAS beamformer backend. EchoFrame only builds it when configured with -`-DEF_USE_FFDAS=ON`. +A submodule (`github.com/BrainEchoLab/ffdas`) providing the CUDA +DAS beamformer backend. EchoFrame always builds it, so the submodule must +be initialised before configuring; a missing submodule is a hard CMake +error rather than a silently reduced build. ## Storage diff --git a/echoframe/cpp/src/CMakeLists.txt b/echoframe/cpp/src/CMakeLists.txt index e2cb6c8..528c45a 100644 --- a/echoframe/cpp/src/CMakeLists.txt +++ b/echoframe/cpp/src/CMakeLists.txt @@ -116,50 +116,47 @@ option(EF_BUILD_MEX "Build the MATLAB MEX interface" ON) option(EF_BUILD_PYTHON "Build the Python extension module" ON) option(ENABLE_CUDA_TIMING "Enable CUDA kernel timing instrumentation" OFF) -option(EF_USE_FFDAS "Build the ffdas-backed DAS beamformer" OFF) set(EF_FFDAS_CUDA_ARCHITECTURES "" CACHE STRING "CUDA architectures for the ffdas submodule fallback build") -if(EF_USE_FFDAS) - find_package(ffdas CONFIG QUIET) - if(NOT TARGET ffdas::ffdas AND NOT TARGET ffdas) - set(EF_FFDAS_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/../libs/ffdas") - if(EXISTS "${EF_FFDAS_SOURCE}/CMakeLists.txt") - if(CMAKE_VERSION VERSION_LESS 3.26) - message(FATAL_ERROR - "EF_USE_FFDAS=ON needs CMake >= 3.26 to build the ffdas submodule. Update the dev image or install CMake from PyPI/Kitware.") - endif() - set(EF_SAVED_CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") - set(EF_SAVED_CUDA_SEPARABLE_COMPILATION - "${CMAKE_CUDA_SEPARABLE_COMPILATION}") - if(NOT EF_FFDAS_CUDA_ARCHITECTURES) - if(CUDAToolkit_VERSION VERSION_GREATER_EQUAL 13) - set(EF_FFDAS_CUDA_ARCHITECTURES "80;86;89;90;100;103;121") - else() - set(EF_FFDAS_CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") - endif() - endif() - set(CMAKE_CUDA_ARCHITECTURES "${EF_FFDAS_CUDA_ARCHITECTURES}") - set(CMAKE_CUDA_SEPARABLE_COMPILATION OFF) - add_subdirectory("${EF_FFDAS_SOURCE}" "${CMAKE_BINARY_DIR}/ffdas") - if(TARGET ffdas) - set_property(TARGET ffdas PROPERTY - CUDA_SEPARABLE_COMPILATION OFF) - endif() - set(CMAKE_CUDA_ARCHITECTURES "${EF_SAVED_CUDA_ARCHITECTURES}") - set(CMAKE_CUDA_SEPARABLE_COMPILATION - "${EF_SAVED_CUDA_SEPARABLE_COMPILATION}") - else() +find_package(ffdas CONFIG QUIET) +if(NOT TARGET ffdas::ffdas AND NOT TARGET ffdas) + set(EF_FFDAS_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/../libs/ffdas") + if(EXISTS "${EF_FFDAS_SOURCE}/CMakeLists.txt") + if(CMAKE_VERSION VERSION_LESS 3.26) message(FATAL_ERROR - "EF_USE_FFDAS=ON requires an installed ffdas package or the submodule at ${EF_FFDAS_SOURCE}") + "Building the ffdas submodule needs CMake >= 3.26. Update the dev image or install CMake from PyPI/Kitware.") endif() - endif() - if(TARGET ffdas::ffdas) - set(EF_FFDAS_TARGET ffdas::ffdas) + set(EF_SAVED_CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") + set(EF_SAVED_CUDA_SEPARABLE_COMPILATION + "${CMAKE_CUDA_SEPARABLE_COMPILATION}") + if(NOT EF_FFDAS_CUDA_ARCHITECTURES) + if(CUDAToolkit_VERSION VERSION_GREATER_EQUAL 13) + set(EF_FFDAS_CUDA_ARCHITECTURES "80;86;89;90;100;103;121") + else() + set(EF_FFDAS_CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") + endif() + endif() + set(CMAKE_CUDA_ARCHITECTURES "${EF_FFDAS_CUDA_ARCHITECTURES}") + set(CMAKE_CUDA_SEPARABLE_COMPILATION OFF) + add_subdirectory("${EF_FFDAS_SOURCE}" "${CMAKE_BINARY_DIR}/ffdas") + if(TARGET ffdas) + set_property(TARGET ffdas PROPERTY + CUDA_SEPARABLE_COMPILATION OFF) + endif() + set(CMAKE_CUDA_ARCHITECTURES "${EF_SAVED_CUDA_ARCHITECTURES}") + set(CMAKE_CUDA_SEPARABLE_COMPILATION + "${EF_SAVED_CUDA_SEPARABLE_COMPILATION}") else() - set(EF_FFDAS_TARGET ffdas) + message(FATAL_ERROR + "The ffdas-backed DAS beamformer needs an installed ffdas package or the submodule at ${EF_FFDAS_SOURCE}. Run: git submodule update --init echoframe/cpp/libs/ffdas") endif() endif() +if(TARGET ffdas::ffdas) + set(EF_FFDAS_TARGET ffdas::ffdas) +else() + set(EF_FFDAS_TARGET ffdas) +endif() # -------------------------------------------------------------------------- # Common include directories # -------------------------------------------------------------------------- @@ -223,10 +220,7 @@ target_link_libraries(echoframe_core PRIVATE CUDA::cudart ) -if(EF_USE_FFDAS) - target_link_libraries(echoframe_core PRIVATE ${EF_FFDAS_TARGET}) - target_compile_definitions(echoframe_core PRIVATE EF_USE_FFDAS=1) -endif() +target_link_libraries(echoframe_core PRIVATE ${EF_FFDAS_TARGET}) target_compile_definitions(echoframe_core PRIVATE NOMINMAX diff --git a/echoframe/cpp/src/beamformer/das/ffdas_beamformer.h b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.h index 49f5195..8b5183a 100644 --- a/echoframe/cpp/src/beamformer/das/ffdas_beamformer.h +++ b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.h @@ -6,13 +6,11 @@ #include +#include + #include "../beamformer.h" #include "das_recon_spec.h" -#ifdef EF_USE_FFDAS -#include -#endif - namespace Beamform { template @@ -20,11 +18,9 @@ class FFDASBeamformer : public Beamformer { private: DASReconSpec dasReconSpec; -#ifdef EF_USE_FFDAS ffdas_handle_t handle{}; ffdas_tensor_desc_t xDesc{}; ffdas_tensor_desc_t outDesc{}; -#endif bool mInitialized{false}; diff --git a/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp index 9cf59cc..f254663 100644 --- a/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp +++ b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp @@ -16,8 +16,6 @@ namespace Beamform { -#ifdef EF_USE_FFDAS -namespace { inline void ffdasErrchk(ffdas_error_t err) { if (err != FFDAS_SUCCESS) { throw std::runtime_error(std::string("ffdas error: ") + @@ -25,7 +23,6 @@ inline void ffdasErrchk(ffdas_error_t err) { } } } // namespace -#endif template FFDASBeamformer::~FFDASBeamformer() { @@ -34,23 +31,14 @@ FFDASBeamformer::~FFDASBeamformer() { template void FFDASBeamformer::initialize() { -#ifndef EF_USE_FFDAS - throw std::runtime_error( - "DAS beamforming requires configuring with -DEF_USE_FFDAS=ON."); -#else static_assert( std::is_same_v, "ffdas beamformer currently supports complex single IQ only."); initGPU(); -#endif } template void FFDASBeamformer::initGPU() { -#ifndef EF_USE_FFDAS - throw std::runtime_error( - "DAS beamforming requires configuring with -DEF_USE_FFDAS=ON."); -#else const size_t nChannels = this->receiveSpec.nActiveChannels; const size_t nVoxels = this->reconSpec.totalSize; const size_t nTx = this->receiveSpec.nTX; @@ -124,12 +112,10 @@ void FFDASBeamformer::initGPU() { 1}; ffdasErrchk(ffdas_create_tensor_desc(&outDesc, 3, outDims, outStrides, FFDAS_C_32F)); -#endif } template void FFDASBeamformer::clearGPU() { -#ifdef EF_USE_FFDAS if (xDesc) ffdas_destroy_tensor_desc(xDesc); if (outDesc) ffdas_destroy_tensor_desc(outDesc); if (handle) ffdas_destroy(handle); @@ -141,15 +127,10 @@ void FFDASBeamformer::clearGPU() { if (dasReconSpec.d_tgcVector) cudaFree(dasReconSpec.d_tgcVector); if (dasReconSpec.d_sourceDirections) cudaFree(dasReconSpec.d_sourceDirections); -#endif } template void FFDASBeamformer::process() { -#ifndef EF_USE_FFDAS - throw std::runtime_error( - "DAS beamforming requires configuring with -DEF_USE_FFDAS=ON."); -#else constexpr int threadsPerBlock = 512; const int nRF = this->receiveSpec.nFastTimeSamples * this->receiveSpec.nTX * this->receiveSpec.nSlowTimeSamples * @@ -167,7 +148,6 @@ void FFDASBeamformer::process() { this->d_BF, static_cast(dasReconSpec.computeType), static_cast(dasReconSpec.algorithm))); gpuErrchk(cudaDeviceSynchronize()); -#endif } } // namespace Beamform diff --git a/echoframe/python/examples/compare_beamformers.py b/echoframe/python/examples/compare_beamformers.py index 02c9b02..570e202 100644 --- a/echoframe/python/examples/compare_beamformers.py +++ b/echoframe/python/examples/compare_beamformers.py @@ -1,7 +1,8 @@ """ Compare EchoFrame Fourier and ffdas DAS beamformers on one stored RF frame. -Requires an EchoFrame Python build configured with ``-DEF_USE_FFDAS=ON``. +Requires an EchoFrame Python build. The DAS beamformer needs ffdas, which is +built by default. Edit ``load_path`` below or pass it on the command line. """ diff --git a/justfile b/justfile index b7f32db..967fde5 100644 --- a/justfile +++ b/justfile @@ -9,7 +9,6 @@ matlab_root := if os() == "windows" { env_var_or_default("Matlab_ROOT_DIR", "C:/ cuda_host_compiler := if os() == "windows" { "" } else { env_var_or_default("CMAKE_CUDA_HOST_COMPILER", "/usr/bin/g++") } default_cuda_architectures := if os() == "windows" { "" } else { `if nvcc --version 2>/dev/null | grep -Eq 'release 1[3-9]\.'; then printf '75;80;86;89;90;100;103;121'; else printf '61;75;86;89;90'; fi` } cuda_architectures := env_var_or_default("CMAKE_CUDA_ARCHITECTURES", default_cuda_architectures) -use_ffdas := env_var_or_default("EF_USE_FFDAS", "OFF") vcpkg_root := env_var_or_default("VCPKG_ROOT", "C:/vcpkg") vcpkg_triplet := "x64-windows-static" test_img := env_var_or_default("HOME", "") / "ef_test.img" @@ -35,7 +34,6 @@ configure: -DEF_BUILD_CLI=ON \ -DEF_BUILD_MEX=ON \ -DEF_BUILD_PYTHON=ON \ - -DEF_USE_FFDAS={{use_ffdas}} \ -DMatlab_ROOT_DIR="{{matlab_root}}" [windows] @@ -44,7 +42,6 @@ configure: # with $VCPKG_ROOT if it's not at C:/vcpkg. cmake -S {{source_dir}} -B {{build_dir}} -G "Visual Studio 17 2022" -A x64 \ -DCMAKE_BUILD_TYPE=Release \ - -DEF_USE_FFDAS={{use_ffdas}} \ -DCMAKE_TOOLCHAIN_FILE="{{vcpkg_root}}/scripts/buildsystems/vcpkg.cmake" \ -DCMAKE_PREFIX_PATH="{{vcpkg_root}}/installed/{{vcpkg_triplet}}" From 09ae207e411221a4d1d0796302abaff6df881b05 Mon Sep 17 00:00:00 2001 From: Samuel Le Meur-Diebolt Date: Wed, 30 Sep 2026 16:10:07 +0100 Subject: [PATCH 05/11] build(cuda): add sm_120 for consumer Blackwell Neither architecture list included sm_120, so RTX 50 series GPUs ran every kernel by JIT-compiling the sm_90 PTX at launch rather than using native SASS. That is correct but not free, and it matters for a real-time acquisition loop. sm_120 requires CUDA 12.8 or newer, so it cannot simply be appended to the existing 12.4 list: on 12.4 through 12.7 nvcc does not know the architecture and the build would fail. Split that case out into its own 12.8 branch instead, and add 120 to the CUDA 13 list and to the ffdas fallback list. The dev image is CUDA 12.8, so it now picks up sm_120 natively. --- README.md | 2 +- echoframe/cpp/src/CMakeLists.txt | 11 +++++++++-- 2 files changed, 10 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index cc1aab0..e019e55 100644 --- a/README.md +++ b/README.md @@ -294,7 +294,7 @@ This section is meant to clarify what changes you can make to the CMake build co | `EF_BUILD_PYTHON` | **ON** | build the Python `echoframe` module | | `EF_BUILD_CLI` | **ON** | build the EchoFrame `cli` (command line) module | | `ENABLE_CUDA_TIMING` | **OFF** | print per-frame CUDA kernel timings from `process()` to the console | -| `CMAKE_CUDA_ARCHITECTURES` | set from the CUDA version: `61;75;86;89;90` on CUDA 12.4+, `75;80;86;89;90;100;103;121` on CUDA 13+ | compute capability list. `CMakeLists.txt` sets it after `project()`, so a value passed on the command line is replaced — edit the list there to narrow it | +| `CMAKE_CUDA_ARCHITECTURES` | set from the CUDA version: `61;75;86;89;90` on CUDA 12.4–12.7, `61;75;86;89;90;120` on CUDA 12.8–12.9, `75;80;86;89;90;100;103;120;121` on CUDA 13+ | compute capability list. `sm_120` covers consumer Blackwell (RTX 50 series) and needs CUDA 12.8 or newer; without it those GPUs JIT the `sm_90` PTX on every kernel launch. `CMakeLists.txt` sets it after `project()`, so a value passed on the command line is replaced — edit the list there to narrow it | Example: disabling Python in Linux: diff --git a/echoframe/cpp/src/CMakeLists.txt b/echoframe/cpp/src/CMakeLists.txt index 528c45a..29833a5 100644 --- a/echoframe/cpp/src/CMakeLists.txt +++ b/echoframe/cpp/src/CMakeLists.txt @@ -102,7 +102,14 @@ endif() if(CUDAToolkit_VERSION VERSION_GREATER_EQUAL 13) message(STATUS "CUDA Toolkit ${CUDAToolkit_VERSION} - using the CUDA 13 architecture list") add_compile_definitions(GGML_CUDA_USE_GRAPHS) # try to enable cuda graphs on cu13 build - set(CMAKE_CUDA_ARCHITECTURES "75;80;86;89;90;100;103;121") # lowest CUDA 13 standard + set(CMAKE_CUDA_ARCHITECTURES "75;80;86;89;90;100;103;120;121") # lowest CUDA 13 standard +elseif(CUDAToolkit_VERSION VERSION_GREATER_EQUAL 12.8) + # sm_120 (consumer Blackwell, e.g. RTX 50 series) needs CUDA 12.8 or newer, + # so it cannot go in the 12.4 list below. Without it those GPUs fall back to + # JIT-compiling the sm_90 PTX on every kernel launch. + message(STATUS "CUDA Toolkit ${CUDAToolkit_VERSION} - using the CUDA 12.8 architecture list") + add_compile_definitions(GGML_CUDA_USE_GRAPHS) # try to enable cuda graphs on cu12 build + set(CMAKE_CUDA_ARCHITECTURES "61;75;86;89;90;120") # lowest CUDA 12 standard + consumer Blackwell elseif(CUDAToolkit_VERSION VERSION_GREATER_EQUAL 12.4) message(STATUS "CUDA Toolkit ${CUDAToolkit_VERSION} - using the CUDA 12 architecture list") add_compile_definitions(GGML_CUDA_USE_GRAPHS) # try to enable cuda graphs on cu12 build @@ -132,7 +139,7 @@ if(NOT TARGET ffdas::ffdas AND NOT TARGET ffdas) "${CMAKE_CUDA_SEPARABLE_COMPILATION}") if(NOT EF_FFDAS_CUDA_ARCHITECTURES) if(CUDAToolkit_VERSION VERSION_GREATER_EQUAL 13) - set(EF_FFDAS_CUDA_ARCHITECTURES "80;86;89;90;100;103;121") + set(EF_FFDAS_CUDA_ARCHITECTURES "80;86;89;90;100;103;120;121") else() set(EF_FFDAS_CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") endif() From 3f548ae1a9342fd7c00543254b8c73248c93e3e9 Mon Sep 17 00:00:00 2001 From: Samuel Le Meur-Diebolt Date: Wed, 30 Sep 2026 16:25:41 +0100 Subject: [PATCH 06/11] fix(das): close the Beamform namespace after removing the anonymous one Dropping the EF_USE_FFDAS guards also dropped the `namespace {` that wrapped ffdasErrchk, but its closing brace stayed. That brace closed namespace Beamform early, so every FFDASBeamformer definition after it sat outside the namespace and the file no longer compiled: ffdas_beamformer.t.hpp(28): error: FFDASBeamformer is not a template The errors that followed all pointed at fourier_imaging.t.hpp, which is included later and was an innocent bystander. --- build_scripts/build_echoframe.sh | 26 +++++++++++++++---- .../src/beamformer/das/ffdas_beamformer.t.hpp | 1 - 2 files changed, 21 insertions(+), 6 deletions(-) mode change 100644 => 100755 build_scripts/build_echoframe.sh diff --git a/build_scripts/build_echoframe.sh b/build_scripts/build_echoframe.sh old mode 100644 new mode 100755 index fcb2ddc..c6cd8f7 --- a/build_scripts/build_echoframe.sh +++ b/build_scripts/build_echoframe.sh @@ -13,7 +13,8 @@ # $HOME/MATLAB/R2024a) # --no-mex skip the MATLAB MEX layer # --no-python skip the Python extension module -# --no-cli skip the CLI (this is what needs vcpkg) +# --cli also build the CLI (off by default: needs vcpkg) +# --configure-only configure the tree, then stop without building # --clean delete the build tree first # -j N parallel jobs (default: nproc) # -h, --help this message @@ -31,7 +32,8 @@ BUILD_TYPE="Release" MATLAB_ROOT="${Matlab_ROOT_DIR:-${HOME}/MATLAB/R2024a}" BUILD_MEX="ON" BUILD_PYTHON="ON" -BUILD_CLI="ON" +BUILD_CLI="OFF" +CONFIGURE_ONLY=0 CLEAN=0 JOBS="$(nproc)" @@ -54,7 +56,8 @@ while [ $# -gt 0 ]; do --matlab-root) MATLAB_ROOT="$2"; shift 2 ;; --no-mex) BUILD_MEX="OFF"; shift ;; --no-python) BUILD_PYTHON="OFF"; shift ;; - --no-cli) BUILD_CLI="OFF"; shift ;; + --cli) BUILD_CLI="ON"; shift ;; + --configure-only) CONFIGURE_ONLY=1; shift ;; --clean) CLEAN=1; shift ;; -j) JOBS="$2"; shift 2 ;; -h|--help) sed -n '2,24p' "${BASH_SOURCE[0]}" | sed 's/^# \{0,1\}//'; exit 0 ;; @@ -78,8 +81,14 @@ if [ ! -e "${REPO_ROOT}/echoframe/cpp/libs/GSL/include/gsl/gsl" ]; then die "GSL submodule is empty. Run: git submodule update --init --recursive" fi +# ffdas is a gitlink too, and it is not optional: CMake configures it in +# unconditionally now, so an empty checkout stops the configure step. +if [ ! -e "${REPO_ROOT}/echoframe/cpp/libs/ffdas/CMakeLists.txt" ]; then + die "ffdas submodule is empty. Run: git submodule update --init --recursive" +fi + if [ "${BUILD_CLI}" = "ON" ] && [ -z "${VCPKG_ROOT:-}" ]; then - die "EF_BUILD_CLI is ON but VCPKG_ROOT is unset. Set it, or pass --no-cli. + die "EF_BUILD_CLI is ON but VCPKG_ROOT is unset. Set it, or drop --cli. Dependencies come from echoframe/cpp/src/vcpkg.json in manifest mode; no manual 'vcpkg install' is needed." fi @@ -95,7 +104,10 @@ CUDA_VER="$(nvcc --version | sed -n 's/.*release \([0-9][0-9]*\.[0-9][0-9]*\).*/ CUDA_MAJOR="${CUDA_VER%%.*}" if [ -z "${CMAKE_CUDA_ARCHITECTURES:-}" ]; then if [ "${CUDA_MAJOR}" -ge 13 ]; then - CMAKE_CUDA_ARCHITECTURES="75;80;86;89;90;100;103;121" + CMAKE_CUDA_ARCHITECTURES="75;80;86;89;90;100;103;120;121" + elif [ "${CUDA_VER}" = "12.8" ] || [ "${CUDA_VER}" = "12.9" ]; then + # sm_120 (consumer Blackwell) needs CUDA 12.8 or newer. + CMAKE_CUDA_ARCHITECTURES="61;75;86;89;90;120" else CMAKE_CUDA_ARCHITECTURES="61;75;86;89;90" fi @@ -123,6 +135,10 @@ CONFIGURE_ARGS=( command -v ninja >/dev/null && CONFIGURE_ARGS+=(-G Ninja) cmake "${CONFIGURE_ARGS[@]}" +if [ "${CONFIGURE_ONLY}" = "1" ]; then + note "configured only; not building (--configure-only)" + exit 0 +fi cmake --build "${BUILD_DIR}" -j "${JOBS}" # --- report --------------------------------------------------------------- diff --git a/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp index f254663..b96be7e 100644 --- a/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp +++ b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp @@ -22,7 +22,6 @@ inline void ffdasErrchk(ffdas_error_t err) { ffdas_error_string(err)); } } -} // namespace template FFDASBeamformer::~FFDASBeamformer() { From b24eb340b7eb5dacc9ca42749faa3ba67d398a07 Mon Sep 17 00:00:00 2001 From: Samuel Le Meur-Diebolt Date: Wed, 30 Sep 2026 16:25:44 +0100 Subject: [PATCH 07/11] fix(build): link ffdas publicly so MEX and Python targets see its headers ffdas_beamformer.h now includes unconditionally, but the ffdas link was PRIVATE to echoframe_core, so the MEX and Python targets that include that header failed: ffdas_beamformer.h:9:10: fatal error: ffdas.h: No such file or directory They also link libffdas, since the out-of-line definitions call ffdas_*. This was hidden before because those targets never saw the EF_USE_FFDAS define, so the include was skipped rather than satisfied. --- echoframe/cpp/src/CMakeLists.txt | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/echoframe/cpp/src/CMakeLists.txt b/echoframe/cpp/src/CMakeLists.txt index 29833a5..0b32972 100644 --- a/echoframe/cpp/src/CMakeLists.txt +++ b/echoframe/cpp/src/CMakeLists.txt @@ -227,7 +227,10 @@ target_link_libraries(echoframe_core PRIVATE CUDA::cudart ) -target_link_libraries(echoframe_core PRIVATE ${EF_FFDAS_TARGET}) +# PUBLIC, not PRIVATE: ffdas_beamformer.h includes and the out-of-line +# definitions call ffdas_*, so every target that includes that header needs the +# headers and the library, not just echoframe_core. +target_link_libraries(echoframe_core PUBLIC ${EF_FFDAS_TARGET}) target_compile_definitions(echoframe_core PRIVATE NOMINMAX From bcb1051c25d20218b616f85c39762148bada749e Mon Sep 17 00:00:00 2001 From: Samuel Le Meur-Diebolt Date: Wed, 30 Sep 2026 16:25:49 +0100 Subject: [PATCH 08/11] build: make build_echoframe.sh the single Linux build implementation The justfile and build_echoframe.sh each carried their own copy of the CMake configure line, and they drifted: EF_USE_FFDAS reached the justfile but not the script, so a build through the script silently produced a DAS-less binary. Two flag lists will drift again. Keep both entry points -- build_echoframe.bat already has a .sh counterpart, and the script is usable without just -- but make the justfile's configure recipe call the script with --configure-only, so the flags and the preflight checks exist once. Add a just report recipe for the artefact listing the script ended with. While consolidating, fix three things in the script: - It was committed without the executable bit, so the invocation the README documents failed with "Permission denied". - The CLI defaulted to ON, which needs vcpkg. It is now opt-in via --cli, so a host without vcpkg can build; --no-cli is gone with it. - Add the ffdas submodule to the preflight checks. It is a gitlink like GSL and is no longer optional, so an empty checkout should fail with the submodule command rather than deep inside CMake. Also add sm_120 to the script's CUDA 12.8/12.9 branch, matching the CMakeLists change. --- build_scripts/README.md | 23 +++++++++++++++---- justfile | 49 +++++++++++++++++++++++++++++------------ 2 files changed, 54 insertions(+), 18 deletions(-) diff --git a/build_scripts/README.md b/build_scripts/README.md index 11b7a4f..53f099a 100644 --- a/build_scripts/README.md +++ b/build_scripts/README.md @@ -6,7 +6,7 @@ it, set up its dependencies, or check that a build works. | Script | What it does | |---|---| | `build_echoframe.bat` | Automated Windows build. Walks the CUDA / MATLAB / Visual Studio combinations and reports what succeeded. The clearest reference for how a working Windows build is configured. Builds into `build_/` at the repository root; set `EF_BUILD_ROOT` to put the build trees somewhere else. | -| `build_echoframe.sh` | Single-configuration Linux build — the counterpart of the `.bat`, minus the matrix: one configure, one build. Flags: `--build-dir`, `--build-type`, `--matlab-root`, `--no-mex`, `--no-python`, `--no-cli`, `--clean`, `-j N`. Checks for the GSL submodule and `VCPKG_ROOT` up front, and picks `CMAKE_CUDA_ARCHITECTURES` from the detected toolkit the same way `CMakeLists.txt` does. | +| `build_echoframe.sh` | Single-configuration Linux/macOS build — the counterpart of the `.bat`, minus the matrix: one configure, one build. This is the Linux build implementation: the `justfile`'s `configure` recipe delegates here, so the flag list and the preflight checks live in one place. Flags: `--build-dir`, `--build-type`, `--matlab-root`, `--no-mex`, `--no-python`, `--cli`, `--configure-only`, `--clean`, `-j N`. Checks the GSL and ffdas submodules, MATLAB and vcpkg up front, and picks `CMAKE_CUDA_ARCHITECTURES` from the detected toolkit the same way `CMakeLists.txt` does. | | `pack_build.ps1` | Packages a `Release/` output directory as `binaries/build_/build.zip`, with the `build/Release/` layout the README documents. Called by `build_echoframe.bat`; can be run on its own to repackage without rebuilding. | | `setup_vcpkg.ps1` | Clones and bootstraps vcpkg, then installs the ports EchoFrame needs. Defaults to `C:\vcpkg`; override with `-VcpkgRoot`. | | `verify_builds.ps1` | Build verification harness. Copies the tree to a short path to dodge Windows `MAX_PATH`, then runs the MEX, CLI and Python builds. Non-zero exit on failure. Flags: `-SkipMex`, `-SkipCli`, `-Keep`, `-Python`, `-WorkDir`. | @@ -19,11 +19,26 @@ Run them from the repository root: pwsh -File build_scripts/verify_builds.ps1 ``` -On Linux: +On Linux and macOS, use the `justfile` at the repository root. `just configure` calls +`build_echoframe.sh`, so the flags below reach the same script the justfile uses: ```bash -build_scripts/build_echoframe.sh # everything -build_scripts/build_echoframe.sh --no-mex # no MATLAB installed +just configure # preflight, then configure +just build # build every target +just report # list artefacts with sizes +just build-target echoframe # build a single target + +EF_BUILD_MEX=OFF just configure # no MATLAB on this host +EF_BUILD_PYTHON=OFF just configure # skip the Python extension +EF_BUILD_CLI=ON just configure # opt in to the CLI (needs vcpkg) +``` + +The CLI is off by default because it pulls Matio, HDF5 and ZLIB through vcpkg. The +script can also be run directly, which is what the justfile does: + +```bash +build_scripts/build_echoframe.sh # everything +build_scripts/build_echoframe.sh --no-mex --configure-only ``` Each resolves the repository root from its own location, so they also work when diff --git a/justfile b/justfile index 967fde5..d0492b3 100644 --- a/justfile +++ b/justfile @@ -7,8 +7,14 @@ jobs := if os() == "windows" { env_var_or_default("NUMBER_OF_PROCESSORS", "4") } make_program := if os() == "windows" { "" } else { `command -v make` } matlab_root := if os() == "windows" { env_var_or_default("Matlab_ROOT_DIR", "C:/Program Files/MATLAB/R2024a") } else { `printf '%s' "${Matlab_ROOT_DIR:-$HOME/MATLAB/R2024a}"` } cuda_host_compiler := if os() == "windows" { "" } else { env_var_or_default("CMAKE_CUDA_HOST_COMPILER", "/usr/bin/g++") } -default_cuda_architectures := if os() == "windows" { "" } else { `if nvcc --version 2>/dev/null | grep -Eq 'release 1[3-9]\.'; then printf '75;80;86;89;90;100;103;121'; else printf '61;75;86;89;90'; fi` } +default_cuda_architectures := if os() == "windows" { "" } else { `if nvcc --version 2>/dev/null | grep -Eq 'release 1[3-9]\.'; then printf '75;80;86;89;90;100;103;120;121'; elif nvcc --version 2>/dev/null | grep -Eq 'release 12\.[89]'; then printf '61;75;86;89;90;120'; else printf '61;75;86;89;90'; fi` } cuda_architectures := env_var_or_default("CMAKE_CUDA_ARCHITECTURES", default_cuda_architectures) +# Component toggles, named after the CMake options they forward. The CLI is +# opt-in because it pulls Matio, HDF5 and ZLIB through vcpkg, which not every +# build host has. +build_cli := env_var_or_default("EF_BUILD_CLI", "OFF") +build_mex := env_var_or_default("EF_BUILD_MEX", "ON") +build_python := env_var_or_default("EF_BUILD_PYTHON", "ON") vcpkg_root := env_var_or_default("VCPKG_ROOT", "C:/vcpkg") vcpkg_triplet := "x64-windows-static" test_img := env_var_or_default("HOME", "") / "ef_test.img" @@ -21,20 +27,16 @@ help: # Configure and build. all: configure build -# Configure the CMake build (Unix Makefiles on Linux/macOS, Visual Studio 2022 on Windows). +# Configure the build (delegates to build_echoframe.sh on Linux/macOS). [unix] configure: - # Set CUDA archs up front; EchoFrame's CMakeLists sets them too late for - # CMake's initial CUDA compiler checks. - cmake -S {{source_dir}} -B {{build_dir}} \ - -DCMAKE_BUILD_TYPE=Release \ - -DCMAKE_MAKE_PROGRAM={{make_program}} \ - -DCMAKE_CUDA_HOST_COMPILER={{cuda_host_compiler}} \ - -DCMAKE_CUDA_ARCHITECTURES="{{cuda_architectures}}" \ - -DEF_BUILD_CLI=ON \ - -DEF_BUILD_MEX=ON \ - -DEF_BUILD_PYTHON=ON \ - -DMatlab_ROOT_DIR="{{matlab_root}}" + #!/usr/bin/env bash + set -euo pipefail + flags=(--configure-only --build-dir "{{build_dir}}") + [ "{{build_mex}}" = "ON" ] || flags+=(--no-mex) + [ "{{build_python}}" = "ON" ] || flags+=(--no-python) + [ "{{build_cli}}" = "OFF" ] || flags+=(--cli) + "{{justfile_directory()}}/build_scripts/build_echoframe.sh" "${flags[@]}" [windows] configure: @@ -42,8 +44,27 @@ configure: # with $VCPKG_ROOT if it's not at C:/vcpkg. cmake -S {{source_dir}} -B {{build_dir}} -G "Visual Studio 17 2022" -A x64 \ -DCMAKE_BUILD_TYPE=Release \ + -DEF_BUILD_CLI={{build_cli}} \ + -DEF_BUILD_MEX={{build_mex}} \ + -DEF_BUILD_PYTHON={{build_python}} \ -DCMAKE_TOOLCHAIN_FILE="{{vcpkg_root}}/scripts/buildsystems/vcpkg.cmake" \ - -DCMAKE_PREFIX_PATH="{{vcpkg_root}}/installed/{{vcpkg_triplet}}" + -DCMAKE_PREFIX_PATH="{{vcpkg_root}}/installed/{{vcpkg_triplet}}" \ + -DMatlab_ROOT_DIR="{{matlab_root}}" + +# List the build artefacts, with sizes. +[unix] +report: + #!/usr/bin/env bash + set -euo pipefail + found=0 + while IFS= read -r artefact; do + printf '%s (%s)\n' "$(basename "$artefact")" "$(du -h "$artefact" | cut -f1)" + found=1 + done < <(find "{{build_dir}}" -maxdepth 2 \ + \( -name 'echoframe_mex.mexa64' -o -name 'storage.mexa64' \ + -o -name 'echoframe*.so' -o -name 'echoframe*.pyd' \ + -o -name 'echoframe_cli' \) -type f 2>/dev/null) + [ "$found" = "1" ] || echo "(no artefacts matched -- check the build log)" # Build all targets. [unix] From 365021f8199ca7702f688c4f9ac5bc4cc93b3340 Mon Sep 17 00:00:00 2001 From: Samuel Le Meur-Diebolt Date: Wed, 30 Sep 2026 14:56:36 +0100 Subject: [PATCH 09/11] fix(cmake): stop the ffdas build from exporting its C++ runtime ffdas links the C++ runtime statically so it cannot clash with the runtime its host already has, but static does not mean private: without a version script the linker also puts those symbols in libffdas.so's dynamic symbol table, and the loader then binds the whole process to a second copy of the standard library. The two copies keep separate std::locale::id tables, so a std::use_facet<> lookup keyed on the host's id misses the facet the host registered. It throws std::bad_cast, and with exceptions off on std::cout the stream is marked bad -- silently, since nothing else reads rdstate. Every later write to it is discarded, which cut the banners off mid-row: | Device PY: after EchoFrame() Confirmed with a probe on printGPUInfo: rdstate goes 0 -> 5 (badbit | failbit) across the first std::cout << dev and never recovers. Only the first insertion of an arithmetic type fails, because the string insertions before it go straight to the streambuf; every banner after it is dropped. Why it needs matplotlib to show up is load order -- importing it perturbs the lookup enough for libffdas' copy to win -- which is why echoframe_cli on the same machine looks fine and it reads as a matplotlib bug. Restrict the fallback ffdas build to its ffdas_* C API with cmake/ffdas_exports.map, and check the linker's actual output rather than trusting the flag: a build that exports anything else now fails, as does a prebuilt ffdas found through find_package(), which cannot be fixed from here but should not fail silently either. Verified on Linux / CUDA 12.8 / RTX 5080, Python 3.13 and 3.14: banners are complete with and without matplotlib imported, a DAS frame still runs through ffdas, and pytest passes. Not verified: the MEX layer (no MATLAB here), and the Windows build, where a static-CRT DLL exports nothing beyond __declspec(dllexport) and the module is a no-op. --- echoframe/cpp/src/CMakeLists.txt | 60 ++++--- .../cpp/src/cmake/CheckFfdasExports.cmake | 15 ++ .../cpp/src/cmake/Modules/ffdas_exports.cmake | 161 ++++++++++++++++++ echoframe/cpp/src/cmake/ffdas_exports.map | 30 ++++ 4 files changed, 240 insertions(+), 26 deletions(-) create mode 100644 echoframe/cpp/src/cmake/CheckFfdasExports.cmake create mode 100644 echoframe/cpp/src/cmake/Modules/ffdas_exports.cmake create mode 100644 echoframe/cpp/src/cmake/ffdas_exports.map diff --git a/echoframe/cpp/src/CMakeLists.txt b/echoframe/cpp/src/CMakeLists.txt index 0b32972..bdfb8a4 100644 --- a/echoframe/cpp/src/CMakeLists.txt +++ b/echoframe/cpp/src/CMakeLists.txt @@ -65,6 +65,10 @@ project(EchoFrame VERSION 1.0 LANGUAGES CXX CUDA) # find_package(matio MODULE ...) below. list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules") +# Keeps the fallback ffdas build from exporting its statically linked C++ runtime +# (see cmake/ffdas_exports.map for what goes wrong when it does). +include("${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules/ffdas_exports.cmake") + # -------------------------------------------------------------------------- # Global compile/architectures options @@ -129,38 +133,42 @@ set(EF_FFDAS_CUDA_ARCHITECTURES "" CACHE STRING find_package(ffdas CONFIG QUIET) if(NOT TARGET ffdas::ffdas AND NOT TARGET ffdas) set(EF_FFDAS_SOURCE "${CMAKE_CURRENT_SOURCE_DIR}/../libs/ffdas") - if(EXISTS "${EF_FFDAS_SOURCE}/CMakeLists.txt") - if(CMAKE_VERSION VERSION_LESS 3.26) - message(FATAL_ERROR - "Building the ffdas submodule needs CMake >= 3.26. Update the dev image or install CMake from PyPI/Kitware.") - endif() - set(EF_SAVED_CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") - set(EF_SAVED_CUDA_SEPARABLE_COMPILATION - "${CMAKE_CUDA_SEPARABLE_COMPILATION}") - if(NOT EF_FFDAS_CUDA_ARCHITECTURES) - if(CUDAToolkit_VERSION VERSION_GREATER_EQUAL 13) - set(EF_FFDAS_CUDA_ARCHITECTURES "80;86;89;90;100;103;120;121") - else() - set(EF_FFDAS_CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") - endif() - endif() - set(CMAKE_CUDA_ARCHITECTURES "${EF_FFDAS_CUDA_ARCHITECTURES}") - set(CMAKE_CUDA_SEPARABLE_COMPILATION OFF) - add_subdirectory("${EF_FFDAS_SOURCE}" "${CMAKE_BINARY_DIR}/ffdas") - if(TARGET ffdas) - set_property(TARGET ffdas PROPERTY - CUDA_SEPARABLE_COMPILATION OFF) - endif() - set(CMAKE_CUDA_ARCHITECTURES "${EF_SAVED_CUDA_ARCHITECTURES}") - set(CMAKE_CUDA_SEPARABLE_COMPILATION - "${EF_SAVED_CUDA_SEPARABLE_COMPILATION}") - else() + if(NOT EXISTS "${EF_FFDAS_SOURCE}/CMakeLists.txt") message(FATAL_ERROR "The ffdas-backed DAS beamformer needs an installed ffdas package or the submodule at ${EF_FFDAS_SOURCE}. Run: git submodule update --init echoframe/cpp/libs/ffdas") endif() + if(CMAKE_VERSION VERSION_LESS 3.26) + message(FATAL_ERROR + "Building the ffdas submodule needs CMake >= 3.26. Update the dev image or install CMake from PyPI/Kitware.") + endif() + set(EF_SAVED_CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") + set(EF_SAVED_CUDA_SEPARABLE_COMPILATION + "${CMAKE_CUDA_SEPARABLE_COMPILATION}") + if(NOT EF_FFDAS_CUDA_ARCHITECTURES) + if(CUDAToolkit_VERSION VERSION_GREATER_EQUAL 13) + set(EF_FFDAS_CUDA_ARCHITECTURES "80;86;89;90;100;103;120;121") + else() + set(EF_FFDAS_CUDA_ARCHITECTURES "${CMAKE_CUDA_ARCHITECTURES}") + endif() + endif() + set(CMAKE_CUDA_ARCHITECTURES "${EF_FFDAS_CUDA_ARCHITECTURES}") + set(CMAKE_CUDA_SEPARABLE_COMPILATION OFF) + add_subdirectory("${EF_FFDAS_SOURCE}" "${CMAKE_BINARY_DIR}/ffdas") + if(TARGET ffdas) + set_property(TARGET ffdas PROPERTY + CUDA_SEPARABLE_COMPILATION OFF) + efconfine_ffdas_exports(ffdas) + endif() + set(CMAKE_CUDA_ARCHITECTURES "${EF_SAVED_CUDA_ARCHITECTURES}") + set(CMAKE_CUDA_SEPARABLE_COMPILATION + "${EF_SAVED_CUDA_SEPARABLE_COMPILATION}") endif() if(TARGET ffdas::ffdas) set(EF_FFDAS_TARGET ffdas::ffdas) + # A prebuilt ffdas was linked before this project could reach it, so the + # version script cannot be added here. Check what it exports instead of + # leaving the banner truncation to be diagnosed from a Python session. + efreport_ffdas_exports(ffdas::ffdas) else() set(EF_FFDAS_TARGET ffdas) endif() diff --git a/echoframe/cpp/src/cmake/CheckFfdasExports.cmake b/echoframe/cpp/src/cmake/CheckFfdasExports.cmake new file mode 100644 index 0000000..217d850 --- /dev/null +++ b/echoframe/cpp/src/cmake/CheckFfdasExports.cmake @@ -0,0 +1,15 @@ +# Post-link driver for the ffdas export check. Run as: +# +# cmake -DFFDAS_LIBRARY= -P cmake/CheckFfdasExports.cmake +# +# Exits non-zero, with the reason on stderr, if the library exports anything +# outside its C API. The implementation is in Modules/ffdas_exports.cmake; it is +# a module rather than a plain script so CMakeLists.txt can call the same check on +# a prebuilt ffdas at configure time. +include("${CMAKE_CURRENT_LIST_DIR}/Modules/ffdas_exports.cmake") + +if(NOT DEFINED FFDAS_LIBRARY) + message(FATAL_ERROR "FFDAS_LIBRARY must be set to the ffdas library to check.") +endif() + +efcheck_ffdas_exports("${FFDAS_LIBRARY}") diff --git a/echoframe/cpp/src/cmake/Modules/ffdas_exports.cmake b/echoframe/cpp/src/cmake/Modules/ffdas_exports.cmake new file mode 100644 index 0000000..3c0a547 --- /dev/null +++ b/echoframe/cpp/src/cmake/Modules/ffdas_exports.cmake @@ -0,0 +1,161 @@ +# Confine a fallback ffdas build to its own C API. +# +# ffdas links the C++ runtime statically so it cannot clash with the runtime its +# host already has, but a static runtime in a shared library is still exported by +# default, and the loader then binds the whole process to that second copy of the +# standard library. The two disagree about std::locale::id, so a std::use_facet<> +# lookup misses, throws std::bad_cast, and the stream it was called on goes into +# an error state where every later write is silently discarded -- the Python +# banners stop partway through the first GPU-info row. cmake/ffdas_exports.map +# has the full account. +# +# Applied from the parent project rather than ffdas' own CMakeLists.txt so the +# submodule stays untouched. A prebuilt ffdas found through find_package() is +# outside this module's reach; efreport_ffdas_exports() covers that case. + +# efconfine_ffdas_exports() +# +# Restrict 's dynamic symbol table to the ffdas_* C API. +function(efconfine_ffdas_exports target) + # Apple's linker spells this -exported_symbols_list and takes a different file + # format, and MSVC exports nothing but __declspec(dllexport) symbols from a + # static-CRT library. Neither has a runtime to confine. + if(CMAKE_HOST_WIN32) + return() + endif() + if(CMAKE_HOST_APPLE) + message(WARNING + "No --version-script support on Apple linkers, so ffdas will export " + "its statically linked C++ runtime and the host's std::cout may stop " + "working partway through a banner. Add an exported-symbols list to the " + "ffdas build to fix it.") + return() + endif() + + set(map "${CMAKE_CURRENT_SOURCE_DIR}/cmake/ffdas_exports.map") + if(NOT EXISTS "${map}") + message(FATAL_ERROR + "ffdas version script not found at ${map}; without it the ffdas " + "build exports a private copy of the C++ runtime and silently " + "truncates the Python banner output.") + endif() + + # --version-script, unlike -Bsymbolic, leaves the undefined CUDA runtime + # entry points alone, so ffdas still resolves them against libcudart. + # set_property rather than target_link_options: was created in the + # ffdas subdirectory, not this one. LINK_DEPENDS is what makes an edit to the + # map re-link the library, and through it re-run the check below. + set_property(TARGET ${target} APPEND PROPERTY + LINK_OPTIONS "-Wl,--version-script=${map}") + set_property(TARGET ${target} APPEND PROPERTY LINK_DEPENDS "${map}") + message(STATUS "ffdas exports: restricted to the ffdas_* C API") + + # Verify what the linker actually did, not what the flag asked for: a linker + # that ignores --version-script leaves the runtime exported and the banners + # broken, and this is the only place that says so. A stamped custom target + # rather than add_custom_command(TARGET ... POST_BUILD), which CMake refuses + # for a target from another directory. Depending on the library file as well + # as the target is what makes the stamp go stale when ffdas relinks; + # add_dependencies() only orders the two, it does not invalidate anything. + set(stamp "${CMAKE_CURRENT_BINARY_DIR}/ffdas_exports_checked.stamp") + add_custom_command( + OUTPUT "${stamp}" + COMMAND ${CMAKE_COMMAND} + -DFFDAS_LIBRARY=$ + -P "${CMAKE_CURRENT_SOURCE_DIR}/cmake/CheckFfdasExports.cmake" + COMMAND ${CMAKE_COMMAND} -E touch "${stamp}" + DEPENDS $ + COMMENT "Checking that ffdas exports only its C API" + VERBATIM) + add_custom_target(ef_ffdas_exports_check ALL DEPENDS "${stamp}") + add_dependencies(ef_ffdas_exports_check ${target}) +endfunction() + +# efreport_ffdas_exports() +# +# Check the library an imported ffdas target points at. An install(EXPORT) target +# only carries the per-configuration IMPORTED_LOCATION_, so fall back to +# those when LOCATION comes back empty. +function(efreport_ffdas_exports target) + get_target_property(location ${target} LOCATION) + if(NOT location) + foreach(config IN LISTS CMAKE_CONFIGURATION_TYPES + ITEMS Release Debug RelWithDebInfo MinSizeRel) + get_target_property(location ${target} IMPORTED_LOCATION_${config}) + if(location) + break() + endif() + endforeach() + endif() + + if(NOT location) + message(WARNING + "Could not locate the ffdas library behind ${target}, so its exports " + "went unchecked. If it exports C++ symbols the host's std::cout will " + "silently stop working partway through a banner.") + return() + endif() + efcheck_ffdas_exports("${location}") +endfunction() + +# efcheck_ffdas_exports() +# +# Fail if exports anything outside the ffdas_* C API. Defined symbols only +# -- `local: *` in a version script cannot affect the undefined imports. A +# prebuilt ffdas that trips this needs its own CMakeLists.txt to add +# cmake/ffdas_exports.map; EchoFrame cannot link it away after the fact. +function(efcheck_ffdas_exports library) + if(CMAKE_HOST_WIN32 OR NOT EXISTS "${library}") + return() + endif() + find_program(EF_NM NAMES nm llvm-nm) + if(NOT EF_NM) + message(WARNING + "nm not found, so the exported symbols of ${library} went unchecked. " + "If it exports C++ symbols the host's std::cout will silently stop " + "working partway through a banner.") + return() + endif() + + execute_process( + COMMAND "${EF_NM}" --dynamic --defined-only --format=posix "${library}" + OUTPUT_VARIABLE symbols + RESULT_VARIABLE status) + if(NOT status EQUAL 0) + message(WARNING "Could not read the dynamic symbols of ${library}.") + return() + endif() + + # Each line is " ". Absolute (A) entries are + # linker-defined, not code: the version node the version script itself + # creates. Everything else has to be part of the C API. + string(REPLACE "\n" ";" lines "${symbols}") + set(leaked "") + foreach(line IN LISTS lines) + if(line MATCHES "^([^ \t]+) ([A-Za-z]) ") + set(name "${CMAKE_MATCH_1}") + set(kind "${CMAKE_MATCH_2}") + else() + continue() + endif() + if(kind STREQUAL "A") + continue() + endif() + string(REGEX REPLACE "@.*$" "" name "${name}") + if(NOT name MATCHES "^ffdas_") + list(APPEND leaked "${name}") + endif() + endforeach() + + if(leaked) + list(LENGTH leaked count) + list(GET leaked 0 first) + message(FATAL_ERROR + "${library} exports ${count} symbol(s) outside the ffdas_* C API " + "(${first} and others). It carries a private copy of the C++ runtime, " + "which the loader will interpose on the whole process. Fix the ffdas " + "build to link with cmake/ffdas_exports.map; see the comment there for " + "the failure this causes.") + endif() + message(STATUS "ffdas exports: only the ffdas_* C API (${library})") +endfunction() diff --git a/echoframe/cpp/src/cmake/ffdas_exports.map b/echoframe/cpp/src/cmake/ffdas_exports.map new file mode 100644 index 0000000..a927e8c --- /dev/null +++ b/echoframe/cpp/src/cmake/ffdas_exports.map @@ -0,0 +1,30 @@ +/* Only the ffdas C API leaves this library. + * + * ffdas links the C++ runtime statically (see its CMakeLists.txt) so it stays + * self-contained and cannot clash with the runtime the host already has -- under + * MATLAB, which ships its own, or under CPython. Static, though, does not mean + * private: without a version script the linker puts every one of those runtime + * symbols in the dynamic symbol table, and the loader then binds the *whole + * process* -- the echoframe extension, the CLI, anything else -- to this second, + * private copy of the standard library. + * + * Two copies do not agree. Each has its own std::locale::id table, so a + * std::use_facet<> lookup keyed on the host's id misses the facet the host + * registered, throws std::bad_cast, and the stream it was called on is marked + * bad. That is silent by default: with exceptions off on std::cout the first + * std::cout << leaves the stream in an error state and every later + * write to it is discarded. The banners stop mid-line, partway through the first + * GPU-info row, and only once something else in the process has perturbed the + * lookup order -- importing matplotlib from Python is enough, which is why it + * reads as a matplotlib bug and why echoframe_cli on the same machine looks fine. + * + * `local: *` is what makes the static link private. It applies to defined + * symbols only, so the undefined CUDA runtime entry points ffdas still imports + * resolve exactly as before. + */ +FFDAS_1.0 { + global: + ffdas_*; + local: + *; +}; From 1a04c2756684bf237d2a8bb9048f904b2d9e0866 Mon Sep 17 00:00:00 2001 From: Samuel Le Meur-Diebolt Date: Thu, 1 Oct 2026 15:45:37 +0100 Subject: [PATCH 10/11] feat(das): compare regenerated Fourier and DAS --- .gitignore | 1 + .../src/echoframe_cli/load_specs_from_mat.cpp | 5 + .../matlab/core/imaging/echoframe_setup_das.m | 48 +++- .../core/verasonics/vsx_to_ef_structs.m | 5 + .../python/examples/benchmark_beamformers.py | 92 +++++++ .../python/examples/compare_beamformers.py | 225 ++++++++++++++++-- 6 files changed, 357 insertions(+), 19 deletions(-) create mode 100644 echoframe/python/examples/benchmark_beamformers.py diff --git a/.gitignore b/.gitignore index 33e99e9..32e75a0 100644 --- a/.gitignore +++ b/.gitignore @@ -19,6 +19,7 @@ echoframe/cpp/out/ *.h5 *.mat +/Perrot et al. - 2021 - So you think you can DAS A viewpoint on delay-and-sum beamforming.pdf __pycache__/ diff --git a/echoframe/cpp/src/echoframe_cli/load_specs_from_mat.cpp b/echoframe/cpp/src/echoframe_cli/load_specs_from_mat.cpp index b12ba87..99d966e 100644 --- a/echoframe/cpp/src/echoframe_cli/load_specs_from_mat.cpp +++ b/echoframe/cpp/src/echoframe_cli/load_specs_from_mat.cpp @@ -299,6 +299,11 @@ void loadSpecsFromMat(const std::string &matFile, EchoframeResources &res) { DAS.computeType = hasField(Pm.get(), "dasComputeType") ? scalar(Pm.get(), "dasComputeType") : 0; + if (hasField(Pm.get(), "dasSourceDirections")) { + static auto sourceDirs = vec(Pm.get(), "dasSourceDirections"); + DAS.sourceDirections = sourceDirs.data(); + DAS.useDirectivity = true; + } } else { validateFourierMatStruct(Pm.get()); diff --git a/echoframe/matlab/core/imaging/echoframe_setup_das.m b/echoframe/matlab/core/imaging/echoframe_setup_das.m index d428b23..2c220a5 100644 --- a/echoframe/matlab/core/imaging/echoframe_setup_das.m +++ b/echoframe/matlab/core/imaging/echoframe_setup_das.m @@ -61,11 +61,55 @@ ReconSpec.dasOffsets(:, :, iTx) = single(max(offsets, minOffset)); end -ReconSpec.dasWeights = ones(nz, nx, nTx, 'single') ./ single(max(1, nTx * nActive)); +fNumber = directivity_f_number(ProbeSpec, ReconSpec); +cosAlpha = single(cos(atan(1 / (2 * fNumber)))); + +channelXRow = reshape(channelX, 1, 1, []); +directivityMask = single(zGrid ./ sqrt((xGrid - channelXRow).^2 + zGrid.^2) > double(cosAlpha)); +ReconSpec.dasWeights = zeros(nz, nx, nTx, 'single'); +for iTx = 1:nTx + phase = double(ReconSpec.dasOffsets(:, :, iTx)) + reshape(receiveSamples, nz, nx, []); + valid = directivityMask & single(phase >= 0 & phase < double(ReceiveSpec.nSamplesIQ - 1)); + counts = sum(valid, 3); + ReconSpec.dasWeights(:, :, iTx) = single(counts > 0) ./ single(max(1, nTx * counts)); +end + % RFFormatter stores IQ as I - iQ, so ffdas' usual negative phase rotation % is conjugated here. ReconSpec.dasWavenum = single(4 * pi * double(ProbeSpec.Fc) / double(ReceiveSpec.Fs)); ReconSpec.dasAlgorithm = int32(1); ReconSpec.dasComputeType = int32(0); -ReconSpec.dasSourceDirections = single(repmat([0; 0; 1; cosd(35)], 1, nActive)); +ReconSpec.dasSourceDirections = single(repmat([0; 0; 1; cosAlpha], 1, nActive)); +end + +function fNumber = directivity_f_number(ProbeSpec, ReconSpec) +% -3 dB piston-element directivity cutoff from Perrot et al. Defaults are +% deliberately boring: pitch as element width, Fc as highest useful frequency. +width = getfield_default(ProbeSpec, {'elementWidth', 'ElementWidth', 'width'}, ProbeSpec.pitch); +fmax = getfield_default(ProbeSpec, {'fmax', 'Fmax'}, ProbeSpec.Fc); +if isfield(ProbeSpec, 'bandwidth') + fmax = double(ProbeSpec.Fc) + double(ProbeSpec.bandwidth) / 2; +elseif isfield(ProbeSpec, 'bandwidthFraction') + fmax = double(ProbeSpec.Fc) * (1 + double(ProbeSpec.bandwidthFraction) / 2); +end +lambda = double(ReconSpec.c0) / double(fmax); +theta = linspace(0, pi/2 - 1e-3, 4096); +directivity = abs(cos(theta) .* sinc(double(width) / lambda * sin(theta))); +idx = find(directivity <= 0.71, 1, 'first'); +if isempty(idx) + alpha = theta(end); +else + alpha = theta(idx); +end +fNumber = 1 / (2 * tan(alpha)); +end + +function value = getfield_default(s, names, defaultValue) +value = defaultValue; +for i = 1:numel(names) + if isfield(s, names{i}) && ~isempty(s.(names{i})) + value = s.(names{i}); + return; + end +end end diff --git a/echoframe/matlab/core/verasonics/vsx_to_ef_structs.m b/echoframe/matlab/core/verasonics/vsx_to_ef_structs.m index 37cfb0f..ea2340c 100644 --- a/echoframe/matlab/core/verasonics/vsx_to_ef_structs.m +++ b/echoframe/matlab/core/verasonics/vsx_to_ef_structs.m @@ -28,6 +28,11 @@ ProbeSpec.pitch = Trans.spacingMm / 1e3; % distance between the transducers in x-direction [m] ProbeSpec.Fc = Trans.frequency * 1e6; % center frequency of the transducers [Hz] ProbeSpec.elementPosition = Trans.ElementPos; % in mm +if isfield(Trans, 'ElementWidth') + ProbeSpec.elementWidth = Trans.ElementWidth / 1e3; +elseif isfield(Trans, 'elementWidth') + ProbeSpec.elementWidth = Trans.elementWidth / 1e3; +end %% TransmitSpec steer = [TX.Steer] * 180 / pi; diff --git a/echoframe/python/examples/benchmark_beamformers.py b/echoframe/python/examples/benchmark_beamformers.py new file mode 100644 index 0000000..a3baf95 --- /dev/null +++ b/echoframe/python/examples/benchmark_beamformers.py @@ -0,0 +1,92 @@ +""" +Benchmark EchoFrame Fourier and ffdas DAS beamformers on one stored RF frame. + +Requires an EchoFrame Python build. Pass a folder, *_rf.dat, or *_seq.mat file. +""" + +from __future__ import annotations + +import argparse +import statistics +import time +from pathlib import Path + +import numpy as np + +from compare_beamformers import ( + add_das_spec, + add_fourier_spec, + ef, + find_session_files, + load_path, + load_specs, + normalize_specs, + rf_io, + validate_specs, +) + + +def _disable_pdi(recon: dict) -> dict: + out = dict(recon) + out["getBF"] = True + out["getPDI"] = False + return out + + +def _make_case(name: str, probe: dict, transmit: dict, receive: dict, recon: dict, pdi: dict): + recon = add_fourier_spec(probe, transmit, receive, recon) if name == "Fourier" else add_das_spec(probe, transmit, receive, recon) + recon = _disable_pdi(recon) + _, _, receive, recon, pdi = validate_specs(probe, transmit, receive, recon, pdi) + t0 = time.perf_counter() + core = ef.EchoFrame(ef.make_resources(receive, recon, pdi), use_storage=False) + return receive, core, time.perf_counter() - t0 + + +def _time_process(core, rf: np.ndarray, warmup: int, iterations: int) -> list[float]: + for _ in range(warmup): + core.process(rf, start_storage=False) + times = [] + for _ in range(iterations): + t0 = time.perf_counter() + core.process(rf, start_storage=False) + times.append(time.perf_counter() - t0) + return times + + +def _summarize(times: list[float]) -> str: + return ( + f"median {statistics.median(times) * 1e3:.2f} ms, " + f"mean {statistics.mean(times) * 1e3:.2f} ms, " + f"min {min(times) * 1e3:.2f} ms, " + f"std {(statistics.stdev(times) if len(times) > 1 else 0.0) * 1e3:.2f} ms" + ) + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("path", nargs="?", type=Path, default=load_path) + parser.add_argument("--warmup", type=int, default=3) + parser.add_argument("--iterations", "-n", type=int, default=20) + args = parser.parse_args() + if not args.path or not args.path.exists(): + raise SystemExit("Pass a folder containing ScanParameters.mat/rf_acq.dat or an Effusive *_rf.dat file.") + + seq_path, rf_path = find_session_files(args.path) + probe, transmit, receive, recon, pdi = load_specs(seq_path) + probe, transmit, receive, recon, pdi = normalize_specs(probe, transmit, receive, recon, pdi) + + cases = {} + for name in ("Fourier", "DAS"): + receive_i, core, init_s = _make_case(name, probe, transmit, receive, recon, pdi) + cases[name] = (receive_i, core, init_s) + + rf = rf_io.read_stored_rf(rf_path, cases["Fourier"][0]).ravel(order="F") + print(f"RF: {rf_path}") + print(f"Iterations: {args.iterations} (+{args.warmup} warmup), PDI disabled") + for name, (_, core, init_s) in cases.items(): + times = _time_process(core, rf, args.warmup, args.iterations) + print(f"{name:7} init {init_s * 1e3:.2f} ms, process {_summarize(times)}") + + +if __name__ == "__main__": + main() diff --git a/echoframe/python/examples/compare_beamformers.py b/echoframe/python/examples/compare_beamformers.py index 570e202..bad1940 100644 --- a/echoframe/python/examples/compare_beamformers.py +++ b/echoframe/python/examples/compare_beamformers.py @@ -28,6 +28,20 @@ load_path = Path("") +def _tukey_window(n: int, alpha: float) -> np.ndarray: + if alpha <= 0: + return np.ones(n) + if alpha >= 1: + return np.hanning(n) + x = np.linspace(0, 1, n) + w = np.ones(n) + left = x < alpha / 2 + right = x >= 1 - alpha / 2 + w[left] = 0.5 * (1 + np.cos(2 * np.pi / alpha * (x[left] - alpha / 2))) + w[right] = 0.5 * (1 + np.cos(2 * np.pi / alpha * (x[right] - 1 + alpha / 2))) + return w + + def h5struct_to_dict(group: h5py.Group) -> dict: """Convert a MATLAB v7.3 struct group to a dict.""" out = {} @@ -50,6 +64,38 @@ def h5struct_to_dict(group: h5py.Group) -> dict: return out +def matstruct_to_dict(obj) -> dict: + """Convert a scipy-loaded MATLAB v5 struct to a dict.""" + fields = getattr(obj, "_fieldnames", None) + if fields is None: + return obj + return {name: matstruct_to_dict(getattr(obj, name)) for name in fields} + + +def load_specs(path: Path) -> tuple[dict, dict, dict, dict, dict]: + try: + with h5py.File(path, "r") as f: + return tuple(h5struct_to_dict(f[k]) for k in ("ProbeSpec", "TransmitSpec", "ReceiveSpec", "ReconSpec", "PDISpec")) + except OSError: + from scipy.io import loadmat + + data = loadmat(path, squeeze_me=True, struct_as_record=False) + specs = tuple(matstruct_to_dict(data[k]) for k in ("ProbeSpec", "TransmitSpec", "ReceiveSpec", "ReconSpec", "PDISpec")) + recon = specs[3] + # scipy preserves MATLAB column-major shapes; Python bindings expect C-style + # arrays with coordinate/tx dimensions leading. + if str(recon.get("beamformerType", "")).lower() == "das": + recon = dict(recon) + recon["dasChannelPositions"] = np.asarray(recon["dasChannelPositions"]).T + recon["dasVoxelPositions"] = np.asarray(recon["dasVoxelPositions"]).T + recon["dasOffsets"] = np.transpose(np.asarray(recon["dasOffsets"]), (2, 1, 0)) + recon["dasWeights"] = np.transpose(np.asarray(recon["dasWeights"]), (2, 1, 0)) + if "dasSourceDirections" in recon: + recon["dasSourceDirections"] = np.asarray(recon["dasSourceDirections"]).T + specs = (specs[0], specs[1], specs[2], recon, specs[4]) + return specs + + def normalize_specs(probe: dict, transmit: dict, receive: dict, recon: dict, pdi: dict) -> tuple[dict, dict, dict, dict, dict]: """Map Effusive field names onto EchoFrame Python example names.""" probe = dict(probe) @@ -82,12 +128,140 @@ def normalize_specs(probe: dict, transmit: dict, receive: dict, recon: dict, pdi def find_session_files(path: Path) -> tuple[Path, Path]: - """Accept either example folder or an Effusive *_rf.dat file.""" + """Accept a folder, an RF .dat, or a parameter .mat file.""" if path.is_file(): - if path.name.endswith("_rf.dat"): - return path.with_name(path.name.replace("_rf.dat", "_seq.mat")), path - return path, path.with_name(path.name.replace("_seq.mat", "_rf.dat")) - return path / "ScanParameters.mat", path / "rf_acq.dat" + if path.suffix == ".dat": + mats = sorted(path.parent.glob("*.mat")) + preferred = path.with_name(path.name.replace("_rf.dat", "_seq.mat")) + return (preferred if preferred.exists() else mats[0]), path + dats = sorted(path.parent.glob("*rf*.dat")) + return path, dats[0] + mat = path / "ScanParameters.mat" + if not mat.exists(): + mats = sorted(path.glob("*.mat")) + mat = mats[0] + rf = path / "rf_acq.dat" + if not rf.exists(): + rfs = sorted(path.glob("*rf*.dat")) + rf = rfs[0] + return mat, rf + + +def _mode_setup(receive: dict) -> tuple[int, float, np.ndarray, np.ndarray, bool]: + n_samples = int(np.asarray(receive["nSamples"]).ravel()[0]) + samples_per_wavelength = int(np.asarray(receive["samplesPerWavelength"]).ravel()[0]) + n_factor = round(4 / samples_per_wavelength) + mode = str(np.asarray(receive["samplingMode"]).item()) + if mode == "BS50BW": + nz = n_samples * n_factor + freq_window = np.arange(int(np.ceil(nz / 4 - nz / 16)), int(np.floor(nz / 4 + nz / 16))) + freq_mapping = np.fft.fftshift(np.arange(nz // 8)) + return nz, 1.0, freq_window, freq_mapping, True + if mode == "BS67BW": + nz = n_samples * 2 + lower = nz // 2 - 1 + freq_window = np.arange(lower, lower + nz // 4) + return nz, 2.0, freq_window, np.arange(nz // 4), False + if mode == "BS100BW": + nz = n_samples * n_factor + freq_window = np.arange(int(np.ceil(nz / 4 - nz / 8)), int(np.floor(nz / 4 + nz / 8))) + freq_mapping = np.fft.fftshift(np.arange(nz // 4)) + return nz, 1.0, freq_window, freq_mapping, True + if mode == "NS200BW": + nz = n_samples * n_factor + return nz, 1.0, np.arange(nz // 2), np.arange(nz // 2), False + raise ValueError(f"Unsupported samplingMode {mode!r}") + + +def add_fourier_spec(probe: dict, transmit: dict, receive: dict, recon: dict) -> dict: + """Return a ReconSpec copy with the Fourier fields EchoFrame expects.""" + out = copy.deepcopy(recon) + out["beamformerType"] = "Fourier" + + nz, freq_scale, freq_window, freq_mapping, freq_axis_shift = _mode_setup(receive) + fs = float(np.asarray(receive.get("Fs_base", receive["Fs"])).ravel()[0]) + c0 = float(np.asarray(out["c0"]).ravel()[0]) + nx = int(np.asarray(receive["nChannels"]).ravel()[0]) + ntx = int(np.asarray(receive["nTransmissions"]).ravel()[0]) + pitch = float(np.asarray(probe["pitch"]).ravel()[0]) + dz = c0 / (fs * freq_scale) + + frequency_axis = np.fft.fftshift(np.arange(-0.5, 0.5, 1 / nz)) * fs * freq_scale + frequency_axis = frequency_axis[freq_window] + if freq_axis_shift: + frequency_axis = np.fft.fftshift(frequency_axis) + + kx_vector = np.arange(-0.5, 0.5, 1 / nx) * 2 * np.pi / pitch + kz_vector = np.arange(-0.5, 0.5, 1 / nz) * 4 * np.pi / dz + max_k = np.max(np.abs(np.concatenate([kx_vector, kz_vector]))) + kx_vector = np.fft.fftshift(kx_vector / max_k).astype(np.float32) + kz_vector = np.fft.fftshift(kz_vector / max_k).astype(np.float32) + gam = 2 / nz + + kz_vector = np.fft.fftshift(kz_vector[freq_window]) + kz, kx = np.meshgrid(kz_vector, kx_vector, indexing="ij") + spectrum_weighting = np.fft.fftshift( + np.outer(_tukey_window(len(kz_vector), 0.2), _tukey_window(len(kx_vector), 0.2)) + ) + + theta = np.asarray(transmit["steer"], dtype=np.float32).ravel()[:ntx] + cos_t = np.cos(np.deg2rad(theta))[None, None, :] + sin_t = np.sin(np.deg2rad(theta))[None, None, :] + kz3 = kz[:, :, None] + kx3 = kx[:, :, None] + with np.errstate(divide="ignore", invalid="ignore"): + k = (kz3**2 + kx3**2) / (2 * kz3 * cos_t + 2 * kx3 * sin_t) + valid = np.abs(k) < np.abs(kz3) * 2 + k = 2 * k * valid + k[np.isnan(k)] = 0 + + phase_window = spectrum_weighting[:, :, None] * (((2 * kz3 * cos_t) * (kx3 * sin_t)) < 0.01) + kbin = k / gam + koff = np.floor(kbin - 0.5) + kk = np.mod(koff + 1, nz).astype(np.int64) + idx = kk - int(freq_window[0]) + arg = kbin - koff + phase_w = np.exp(-1j * np.pi * arg) * phase_window + + n_iq = int(np.asarray(receive["nSamplesIQ"]).ravel()[0]) + idx[(idx < 0) | (idx >= n_iq)] = 0 + idx = freq_mapping[idx] + scaling = float(ntx * int(np.asarray(out["nz"]).ravel()[0]) * int(np.asarray(out["nx"]).ravel()[0]) * 2) + + out["delayIndices"] = np.transpose(idx.astype(np.int32), (2, 1, 0)) + out["interpolationWeights"] = np.transpose((phase_w / scaling).astype(np.complex64), (2, 1, 0)) + out["frequencyAxis"] = frequency_axis.astype(np.float32) + + transmit_delays = np.asarray(transmit["transmitDelays"], dtype=np.float64) + transmit_delays = np.squeeze(transmit_delays) + if transmit_delays.shape[0] == ntx: + transmit_delays = transmit_delays.T + transmit_delays = transmit_delays[:, :ntx] * 2 * np.pi + active = np.flatnonzero(np.asarray(transmit["apodization"]).ravel() != 0) + delay_ax = np.median(np.diff(transmit_delays[active, :], axis=0), axis=0) + mid_element = int(np.asarray(probe["nElements"]).ravel()[0]) // 2 - 1 + mid_channel = nx // 2 + delay_b = transmit_delays[mid_element, :] - mid_channel * delay_ax + out["planewaveDelays"] = np.concatenate([delay_ax, delay_b]).astype(np.float32) + return out + + +def _directivity_f_number(probe: dict, recon: dict) -> float: + width = float(np.asarray( + probe.get("elementWidth", probe.get("ElementWidth", probe.get("width", probe["pitch"]))) + ).ravel()[0]) + fc = float(np.asarray(probe["Fc"]).ravel()[0]) + fmax = float(np.asarray(probe.get("fmax", probe.get("Fmax", fc))).ravel()[0]) + if "bandwidth" in probe: + fmax = fc + float(np.asarray(probe["bandwidth"]).ravel()[0]) / 2 + elif "bandwidthFraction" in probe: + fmax = fc * (1 + float(np.asarray(probe["bandwidthFraction"]).ravel()[0]) / 2) + wavelength = float(np.asarray(recon["c0"]).ravel()[0]) / fmax + theta = np.linspace(0, np.pi / 2 - 1e-3, 4096) + directivity = np.abs(np.cos(theta) * np.sinc(width / wavelength * np.sin(theta))) + idx = np.flatnonzero(directivity <= 0.71) + alpha = theta[idx[0]] if idx.size else theta[-1] + return 1 / (2 * np.tan(alpha)) def add_das_spec(probe: dict, transmit: dict, receive: dict, recon: dict) -> dict: @@ -167,9 +341,25 @@ def add_das_spec(probe: dict, transmit: dict, receive: dict, recon: dict) -> dic ) offsets[i] = np.maximum(raw_offsets, min_offset) out["dasOffsets"] = offsets.ravel() - out["dasWeights"] = np.full( - ntx * nz * nx, 1.0 / max(1, ntx * nelem), dtype=np.float32 + f_number = _directivity_f_number(probe, out) + cos_alpha = np.float32(np.cos(np.arctan(1 / (2 * f_number)))) + directivity_mask = ( + voxel_z[:, None] + / np.sqrt((voxel_x[:, None] - channel_x[None, :]) ** 2 + voxel_z[:, None] ** 2) + > float(cos_alpha) ) + weights = np.zeros((ntx, nz * nx), dtype=np.float32) + for i in range(ntx): + phase = offsets[i, :, None] + receive_samples + valid = directivity_mask & (phase >= 0) & (phase < float(np.asarray(receive["nSamplesIQ"]).ravel()[0] - 1)) + counts = valid.sum(axis=1) + weights[i] = np.divide( + 1.0, + ntx * counts, + out=np.zeros_like(counts, dtype=np.float32), + where=counts > 0, + ) + out["dasWeights"] = weights.ravel() # RFFormatter stores IQ as I - iQ, so ffdas' usual negative phase rotation # is conjugated here. out["dasWavenum"] = np.float32( @@ -178,18 +368,23 @@ def add_das_spec(probe: dict, transmit: dict, receive: dict, recon: dict) -> dic out["dasAlgorithm"] = np.int32(1) out["dasComputeType"] = np.int32(0) out["dasSourceDirections"] = np.tile( - np.array([0, 0, 1, np.cos(np.deg2rad(35))], dtype=np.float32), nelem + np.array([0, 0, 1, cos_alpha], dtype=np.float32), nelem ) return out +def _ef_image_layout(image: np.ndarray) -> np.ndarray: + """Return EchoFrame image data with Python axis layout.""" + return image.ravel(order="C").reshape(image.shape, order="F") + + def run_once(receive: dict, recon: dict, pdi: dict, rf: np.ndarray) -> tuple[np.ndarray, float]: """Process one RF buffer and return B-mode plus elapsed seconds.""" res = ef.make_resources(receive, recon, pdi) core = ef.EchoFrame(res, use_storage=False) t0 = time.perf_counter() _, bmode, _ = core.process(rf, start_storage=False) - return bmode, time.perf_counter() - t0 + return _ef_image_layout(bmode), time.perf_counter() - t0 def main() -> None: @@ -201,15 +396,11 @@ def main() -> None: raise SystemExit("Pass a folder containing ScanParameters.mat/rf_acq.dat or an Effusive *_rf.dat file.") seq_path, rf_path = find_session_files(args.path) - with h5py.File(seq_path, "r") as f: - probe = h5struct_to_dict(f["ProbeSpec"]) - transmit = h5struct_to_dict(f["TransmitSpec"]) - receive = h5struct_to_dict(f["ReceiveSpec"]) - recon = h5struct_to_dict(f["ReconSpec"]) - pdi = h5struct_to_dict(f["PDISpec"]) + probe, transmit, receive, recon, pdi = load_specs(seq_path) probe, transmit, receive, recon, pdi = normalize_specs(probe, transmit, receive, recon, pdi) - _, _, receive_f, recon_f, pdi_f = validate_specs(probe, transmit, receive, recon, pdi) + recon_f = add_fourier_spec(probe, transmit, receive, recon) + _, _, receive_f, recon_f, pdi_f = validate_specs(probe, transmit, receive, recon_f, pdi) rf = rf_io.read_stored_rf(rf_path, receive_f).ravel(order="F") fourier_bmode, fourier_s = run_once(receive_f, recon_f, pdi_f, rf) @@ -229,7 +420,7 @@ def main() -> None: fig, axes = plt.subplots(1, 2, figsize=(8, 4), constrained_layout=True, dpi=300) for ax, img, title in zip(axes, imgs, titles): im = ax.imshow( - 20 * np.log10(img.T / vmax + np.finfo(np.float32).eps), + 20 * np.log10(img / vmax + np.finfo(np.float32).eps), cmap="gray", vmin=-60, vmax=0, From 22e404f8e40a08ea77df7007ea84530c0216c953 Mon Sep 17 00:00:00 2001 From: Samuel Le Meur-Diebolt Date: Thu, 1 Oct 2026 18:08:00 +0100 Subject: [PATCH 11/11] feat: set DR to -40 in beamformer comparison script --- echoframe/python/examples/compare_beamformers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/echoframe/python/examples/compare_beamformers.py b/echoframe/python/examples/compare_beamformers.py index bad1940..fcfb2d9 100644 --- a/echoframe/python/examples/compare_beamformers.py +++ b/echoframe/python/examples/compare_beamformers.py @@ -422,7 +422,7 @@ def main() -> None: im = ax.imshow( 20 * np.log10(img / vmax + np.finfo(np.float32).eps), cmap="gray", - vmin=-60, + vmin=-40, vmax=0, ) ax.set_title(title)