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/.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/README.md b/README.md index ae60047..e019e55 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. @@ -290,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/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/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/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..89cf631 100644 --- a/echoframe/cpp/libs/README.md +++ b/echoframe/cpp/libs/README.md @@ -11,16 +11,22 @@ 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 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 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..ea3a25e --- /dev/null +++ b/echoframe/cpp/libs/ffdas @@ -0,0 +1 @@ +Subproject commit ea3a25e4d7d8c330fa8c8ee61cdb9a52143fd87f diff --git a/echoframe/cpp/src/CMakeLists.txt b/echoframe/cpp/src/CMakeLists.txt index fb9d167..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 @@ -102,7 +106,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 @@ -116,6 +127,51 @@ 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) +set(EF_FFDAS_CUDA_ARCHITECTURES "" CACHE STRING + "CUDA architectures for the ffdas submodule fallback build") + +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(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() # -------------------------------------------------------------------------- # Common include directories # -------------------------------------------------------------------------- @@ -177,8 +233,12 @@ target_include_directories(echoframe_core target_link_libraries(echoframe_core PRIVATE - CUDA::cudart + CUDA::cudart ) +# 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 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..8b5183a --- /dev/null +++ b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.h @@ -0,0 +1,49 @@ +/** + * @file ffdas_beamformer.h + * @brief ffdas-backed delay-and-sum beamformer. + */ +#pragma once + +#include + +#include + +#include "../beamformer.h" +#include "das_recon_spec.h" + +namespace Beamform { + +template +class FFDASBeamformer : public Beamformer { + private: + DASReconSpec dasReconSpec; + + ffdas_handle_t handle{}; + ffdas_tensor_desc_t xDesc{}; + ffdas_tensor_desc_t outDesc{}; + + 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..b96be7e --- /dev/null +++ b/echoframe/cpp/src/beamformer/das/ffdas_beamformer.t.hpp @@ -0,0 +1,152 @@ +/** + * @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 { + +inline void ffdasErrchk(ffdas_error_t err) { + if (err != FFDAS_SUCCESS) { + throw std::runtime_error(std::string("ffdas error: ") + + ffdas_error_string(err)); + } +} + +template +FFDASBeamformer::~FFDASBeamformer() { + if (mInitialized) clearGPU(); +} + +template +void FFDASBeamformer::initialize() { + static_assert( + std::is_same_v, + "ffdas beamformer currently supports complex single IQ only."); + initGPU(); +} + +template +void FFDASBeamformer::initGPU() { + 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)); +} + +template +void FFDASBeamformer::clearGPU() { + 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); +} + +template +void FFDASBeamformer::process() { + 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()); +} + +} // 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/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: + *; +}; 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..99d966e 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,57 @@ 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; + if (hasField(Pm.get(), "dasSourceDirections")) { + static auto sourceDirs = vec(Pm.get(), "dasSourceDirections"); + DAS.sourceDirections = sourceDirs.data(); + DAS.useDirectivity = true; + } + } 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..2c220a5 --- /dev/null +++ b/echoframe/matlab/core/imaging/echoframe_setup_das.m @@ -0,0 +1,115 @@ +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); + +% 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); +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 + residual = txDelays(:, iTx) - channelX * sinTheta / double(ReconSpec.c0); + bulk = median(residual(active)); + else + bulk = 0; + end + offsets = (zGrid * cosd(theta) + xGrid * sinTheta) * fsOverTwoC ... + + bulk * double(ReceiveSpec.Fs) / 2 - startSamples; + ReconSpec.dasOffsets(:, :, iTx) = single(max(offsets, minOffset)); +end + +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; 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/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/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..9fab3f5 100644 --- a/echoframe/matlab/core/imaging/initialize_image_reconstruction.m +++ b/echoframe/matlab/core/imaging/initialize_image_reconstruction.m @@ -67,10 +67,18 @@ 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); +% 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') + ReconSpec = echoframe_setup_das(ProbeSpec,TransmitSpec,ReceiveSpec,ReconSpec); +else + ReconSpec = echoframe_setup_fourier(ProbeSpec,TransmitSpec,ReceiveSpec,ReconSpec); +end 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/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/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 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/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 new file mode 100644 index 0000000..fcfb2d9 --- /dev/null +++ b/echoframe/python/examples/compare_beamformers.py @@ -0,0 +1,435 @@ +""" +Compare EchoFrame Fourier and ffdas DAS beamformers on one stored RF frame. + +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. +""" + +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 _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 = {} + 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 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) + 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 a folder, an RF .dat, or a parameter .mat file.""" + if path.is_file(): + 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: + """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) + + # 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) + 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] + apod = np.asarray(transmit.get("apodization", np.ones(nelem)), dtype=bool).ravel() + + 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: + residual = transmit_delays[:, i] - channel_x * sin_t / c0 + bulk = float(np.median(residual[apod])) if np.any(apod) else 0.0 + else: + 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() + 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( + 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, 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 _ef_image_layout(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) + probe, transmit, receive, recon, pdi = load_specs(seq_path) + probe, transmit, receive, recon, pdi = normalize_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) + + 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 / vmax + np.finfo(np.float32).eps), + cmap="gray", + vmin=-40, + 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..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]