diff --git a/CITATION.cff b/CITATION.cff index 7cdffcd..584ae38 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -6,7 +6,7 @@ authors: orcid: "https://orcid.org/0000-0002-0720-6955" title: "Quick Ultrasound Processing & Simulation" type: software -version: v1.1.1 +version: v1.2.1 license: Apache-2.0 -date-released: 2024-04-07 +date-released: 2024-08-22 repository-code: "https://github.com/thorstone25/qups" diff --git a/README.md b/README.md index c88a705..b799256 100644 --- a/README.md +++ b/README.md @@ -43,7 +43,10 @@ This package can readily be used to develop new transducer array designs by spec - Documentation via `help` and `doc` -## Installation +## Installation +### Prerequisites +QUPS requires the [Signal Processing Toolbox](https://www.mathworks.com/products/signal.html) and the [Parallel Computing Toolbox](https://www.mathworks.com/products/parallel-computing.html) to be installed. + ### MATLAB R2023b+ & git Starting in MATLAB R2023b+, QUPS and most of it's extension packages can be installed from within MATLAB via [buildtool](https://www.mathworks.com/help/matlab/ref/buildtool.html) if you have setup [git for MATLAB](https://www.mathworks.com/help/matlab/matlab_prog/set-up-git-source-control.html). 1. Install qups @@ -61,7 +64,7 @@ buildtool test ``` ### Legacy Installation -If the above procedure does not work for you, you can manually download and install each [extension](##Extensions). +You can manually download and install each [extension](##Extensions) separately. 1. Download the desired extension packages into a folder adjacent to the "qups" folder e.g. if qups is located at `/path/to/my/qups`, kWave should be downloaded to an adjacent folder `/path/to/my/kWave`. 2. Create a MATLAB [Project](https://www.mathworks.com/help/matlab/matlab_prog/create-projects.html) and add the root folder of the extension to the path e.g. `/path/to/my/kWave`. @@ -218,3 +221,29 @@ First, be sure you can run `nvcc` from a terminal or command-line interface per #### Windows First, setup your system for CUDA per [CUDA installation instructions](https://docs.nvidia.com/cuda/cuda-installation-guide-microsoft-windows/index.html). On Windows you must set the path for both CUDA and the _correct_ MSVC compiler for C/C++. Start a PowerShell terminal within Visual Studio. Run `echo %CUDA_PATH%` to find the base CUDA_PATH and run `echo %VCToolsInstallDir%` to find the MSVC path. Then, in MATLAB, set these paths with `setenv('MW_NVCC_PATH', YOUR_CUDA_BIN_PATH); setenv('VCToolsInstallDir', YOUR_MSVC_PATH);`, where `YOUR_CUDA_BIN_PATH` is the path to the `bin` folder in the `CUDA_PATH` folder. Finally, run `setup CUDA`. From here the proper paths should be added. +## Troubleshooting +### Parpool management +Most functions will perform best with an active `parallel.ThreadPool` available. This can be started with `setup parallel` or `parpool Threads`. Alternatively, a `parallel.ProcessPool` can be started with `parpool Processes` (or `parpool local` on earlier MATLAB releases). Functions that allow specifying a parallel environment can avoid using the active parallel pool by using an argument of `0`, or can limit the number of . + +### GPU management +Some QUPS functions use the currently selected `parallel.gpu.CUDADevice`, or select a `parallel.gpu.CUDADevice` by default. Use [`gpuDevice`](https://www.mathworks.com/help/parallel-computing/parallel.gpu.gpudevice.html) to manually select the gpu. Within a parallel pool, each worker can have a unique selection. By default, GPUs are spread evenly across all workers. +If [CUDA support](#cuda-support) is enabled, [ptx-files](https://www.mathworks.com/help/parallel-computing/parallel.gpu.gpudevice.html) will be compiled to target the currently selected device. If the currently selected device changes, you may need to recompile binares using `UltrasoundSystem.recompileCUDA` or `setup cache`, particularly if the computer contains GPUs from different [virtual architectures](https://docs.nvidia.com/cuda/cuda-compiler-driver-nvcc/index.html#virtual-architecture-feature-list) or have different compute capabilities. + +### OOM (Out of memory) errors +Beamforming and simulation routines tend to require a trade-off between performance and memory usage. While QUPS attempts to balance this, you can still run into OOM errors if the GPU is almost full or if the problem is too large for either the GPU or CPU. There are several options to mitigate this such as: +* (GPU OOM) use `x = gather(x);` to move some variables from the GPU (device) to the CPU (host) - this works on `ChannelData` objects too. +* for beamformers, set the `bsize` optional keyword argument to a smaller value +* (GPU OOM) use `gpuDevice().reset()` to reset the currently selected GPU - NOTE: this will clear any variables that were not moved from the GPU to the CPU with `gather` first +* consider reducing the problem and using a for-loop e.g. by processing an image per frame or per depth, or simulating groups of 1000 scatterers at a time +* (GPU OOM) for `UltrasoundSystem.DAS` and `UltrasoundSystem.greens`, set the `device` argument to `0` to avoid using a GPU device. +* if a `parallel.Pool` is active and a function accepts an optional parallel environment keyword argument `penv`, use a small number e.g. `4` to limit the number of workers, or `0` to avoid using the active `parallel.Pool` +* if a `parallel.Pool` is active, shut it down with `delete(gcp('nocreate'))` and [disable automatic parallel pool creation](https://www.mathworks.com/help/parallel-computing/parallel-preferences.html). + + + + + + + + + diff --git a/build/coverage.xml b/build/coverage.xml index 0977ff0..7e2f223 100644 --- a/build/coverage.xml +++ b/build/coverage.xml @@ -1,10 +1,10 @@ - + /home/thurston/sandbox-2/qups/ - + @@ -14,7 +14,7 @@ - + @@ -140,15 +140,15 @@ - - - + + + - + @@ -277,7 +277,7 @@ - + @@ -350,26 +350,26 @@ - + - - - - - + + + + + - + - + @@ -385,110 +385,110 @@ - + - - - - - - + + + + + + - - + + - - + + - + - - + + - + - - - + + + - - - + + + - - - + + + - - - + + + - - - + + + - - + + - - + + - + @@ -510,47 +510,47 @@ - + - - - - - - - + + + + + + + - + - + - - + + - - - - + + + + @@ -829,92 +829,92 @@ - - - - - - - - + + + + + + + + - - - - - - - - - - - + + + + + + + + + + + - + - - - - - - + + + + + + - + - - - - - - - + + + + + + + - + - - - - - + + + + + - - + + - - - - + + + + @@ -982,19 +982,19 @@ - - + + - - - - - + + + + + @@ -1008,7 +1008,7 @@ - + @@ -1034,68 +1034,68 @@ - - + + - - - - - - - - - + + + + + + + + + - + - + - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + - + - + @@ -1443,7 +1443,7 @@ - + @@ -1499,44 +1499,44 @@ - - - + + + - + - + - + - + - + - - + + - + - - + + @@ -1565,9 +1565,9 @@ - - - + + + @@ -1582,281 +1582,281 @@ - - + + - + - - + + - - + + - + - - - - + + + + - + - - + + - - - + + + - + - + - + - + - + - - + + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - + - - + + - - + + - - + + - - - + + + - - + + - - + + - - + + - - - + + + - - + + - - - + + + - - + + @@ -1890,165 +1890,199 @@ - + - - - - - - - - - + + + + + + - + + + + + + + - - - - - + + + + + - + - - - - - - + + + + + + + + + + + + + - + - - - + + + - - - - - - - - - - - + + + + + + + + + + + - - - - + + + + - - - - - + + + + + - - - - - - + + + + + + + + + + + + - - - + + + - + - + - + + + + + + + - - - + + + + - + - - - + + + - - - - - - + + + + + + + + + + + + + + + + + @@ -2906,28 +2940,28 @@ - - - - - - - - - + + + + + + + + + - - + + - + @@ -2937,24 +2971,24 @@ - - - - - - + + + + + + - + - + @@ -2967,20 +3001,20 @@ - + - + - - + + @@ -2988,9 +3022,9 @@ - - - + + + @@ -3003,107 +3037,107 @@ - + - - - + + + - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + - + - - + + - + - - - - - - - - + + + + + + + + - - + + - + - + - - - - - - - - + + + + + + + + - + @@ -3114,139 +3148,139 @@ - - + + - + - + - + - + - + - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + - - - - - - - - - - + + + + + + + + + + - + - - - - + + + + - - - + + + - + - - - - - + + + + + - - - - - - - + + + + + + + - - + + - + @@ -3256,24 +3290,24 @@ - - - - - - + + + + + + - + - + @@ -3286,20 +3320,20 @@ - + - + - - + + @@ -3307,148 +3341,148 @@ - - - + + + - - - - + + + + - + - - - - - - + + + + + + - + - - + + - - - - - - - + + + + + + + - - - - - - - - - - - + + + + + + + + + + + - + - - + + - + - - - - - - - - - + + + + + + + + + - + - + - - - - - - - - - + + + + + + + + + - - + + - + @@ -3544,7 +3578,7 @@ - + @@ -3554,7 +3588,7 @@ - + @@ -3579,16 +3613,16 @@ - - + + - - - + + + @@ -3840,7 +3874,7 @@ - + @@ -3860,7 +3894,7 @@ - + @@ -3890,7 +3924,7 @@ - + @@ -4134,8 +4168,8 @@ - - + + @@ -4457,31 +4491,31 @@ - + - + - - - + + + - - - - - + + + + + @@ -4598,76 +4632,76 @@ - + - + - + - + - + - - - + + + - + - - - - - - - + + + + + + + - - - - - - - + + + + + + + - + @@ -4798,38 +4832,38 @@ - + - + - + - + - + - - - - - - + + + + + + @@ -4839,120 +4873,120 @@ - - - - - - - - - - - - + + + + + + + + + + + + - + - - - - - + + + + + - - + + - + - - - - - - - - - + + + + + + + + + - - - + + + - - - - + + + + - - - - - - + + + + + + - - + + - - - + + + - - - - - - - + + + + + + + - - - - + + + + @@ -5011,17 +5045,17 @@ - + - + - + @@ -5036,13 +5070,13 @@ - + - + @@ -5059,16 +5093,16 @@ - - + + - - - + + + @@ -5311,7 +5345,7 @@ - + @@ -5323,13 +5357,13 @@ - + - + @@ -5441,9 +5475,9 @@ - - - + + + @@ -5674,31 +5708,31 @@ - + - + - - - + + + - - - - - + + + + + @@ -5799,76 +5833,76 @@ - + - + - + - + - + - - - + + + - + - - - - - - - + + + + + + + - - - - - - - + + + + + + + - + @@ -5976,127 +6010,127 @@ - - - - - + + + + + - - - - - - + + + + + + - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + - - - - - + + + + + - - + + - + - - - - - - - - - + + + + + + + + + - - - + + + - - - - - - - - - - - - + + + + + + + + + + + + - - - - - - - - - - + + + + + + + + + + - - - - + + + + @@ -6135,9 +6169,9 @@ - - - + + + @@ -6146,13 +6180,13 @@ - + - + @@ -6529,14 +6563,14 @@ - - + + - + @@ -6884,7 +6918,7 @@ - + @@ -6989,7 +7023,7 @@ - + @@ -6999,23 +7033,23 @@ - + - - - - - - - - - - - + + + + + + + + + + + @@ -7026,22 +7060,22 @@ - + - + - + - + - + @@ -7243,13 +7277,13 @@ - + - - + + @@ -7347,35 +7381,35 @@ - + - + - - - - - - - - - - - + + + + + + + + + + + - - - - + + + + @@ -7558,9 +7592,9 @@ - - - + + + @@ -7573,26 +7607,26 @@ - + - + - - + + - - + + @@ -7631,70 +7665,70 @@ - - - - - - + + + + + + - + - + - + - + - + - + - + - + - + - + - + @@ -7714,7 +7748,7 @@ - + @@ -7724,7 +7758,7 @@ - + @@ -7734,7 +7768,7 @@ - + @@ -7744,7 +7778,7 @@ - + @@ -7754,7 +7788,7 @@ - + @@ -7764,7 +7798,7 @@ - + @@ -7774,7 +7808,7 @@ - + @@ -7789,29 +7823,29 @@ - - - + + + - + - - - + + + - - + + @@ -7834,109 +7868,109 @@ - - - - - - - + + + + + + + - + - + - - - - - - - - - + + + + + + + + + - + - + - + - + - + - + - + - - - + + + - + - + - - - - - + + + + + - + @@ -7959,19 +7993,19 @@ - + - - - - - - + + + + + + @@ -7989,12 +8023,12 @@ - + - + @@ -8009,46 +8043,46 @@ - + - - - - + + + + - + - + - + - + - + - + @@ -8093,7 +8127,7 @@ - + @@ -8158,33 +8192,33 @@ - - - - + + + + - + - - - - - - + + + + + + - + @@ -8195,20 +8229,20 @@ - - + + - - - - - - + + + + + + @@ -8222,30 +8256,30 @@ - - + + - + - - - - - - - - - - + + + + + + + + + + @@ -8258,7 +8292,7 @@ - + @@ -8291,37 +8325,37 @@ - - - - + + + + - + - + - + - - + + - - - + + + @@ -8345,34 +8379,34 @@ - - - - - + + + + + - + - + - - - - - - + + + + + + @@ -8380,21 +8414,21 @@ - + - + - + - - - - - - - + + + + + + + @@ -8418,32 +8452,32 @@ - + - + - + - + - + - + @@ -8553,31 +8587,31 @@ - - - - + + + + - + - - - + + + - - - + + + @@ -8589,49 +8623,49 @@ - - - - - - + + + + + + - + - + - - - - - - + + + + + + - + - + - - - - - - - + + + + + + + @@ -8643,12 +8677,12 @@ - - - - - - + + + + + + @@ -8694,23 +8728,23 @@ - - - - + + + + - + - + @@ -8729,34 +8763,34 @@ - - - - - + + + + + - + - + - - - - - - + + + + + + @@ -8801,32 +8835,32 @@ - + - + - + - + - + - + @@ -8936,21 +8970,21 @@ - - - - + + + + - + - + @@ -8962,30 +8996,30 @@ - - - - - + + + + + - + - + - - - - - - + + + + + + @@ -9014,12 +9048,12 @@ - - - - - - + + + + + + @@ -9065,25 +9099,25 @@ - - - - + + + + - + - - - + + + @@ -9108,13 +9142,13 @@ - + - + @@ -9128,25 +9162,25 @@ - - - + + + - - - + + + - - - + + + @@ -9186,8 +9220,8 @@ - - + + @@ -9227,13 +9261,13 @@ - - - - - - - + + + + + + + @@ -9273,41 +9307,41 @@ - - + + - + - + - + - + - + - - - + + + @@ -9332,13 +9366,13 @@ - + - + @@ -9348,25 +9382,25 @@ - - - + + + - - - + + + - - - + + + @@ -9395,8 +9429,8 @@ - - + + @@ -9427,13 +9461,13 @@ - - - - - - - + + + + + + + @@ -9465,19 +9499,19 @@ - - - - - - + + + + + + - + @@ -9487,7 +9521,7 @@ - + @@ -9584,31 +9618,31 @@ - - - + + + - + - - - + + + - - - + + + @@ -9621,7 +9655,7 @@ - + @@ -10041,37 +10075,37 @@ - - - + + + - - + + - - + + - - - - + + + + - + @@ -10083,23 +10117,23 @@ - - + + - - + + - + - + @@ -10121,35 +10155,35 @@ - + - - + + - - - + + + - + - + @@ -10183,14 +10217,14 @@ - - + + - + @@ -10210,55 +10244,55 @@ - - + + - - + + - + - + - + - + - + - + - + @@ -10275,19 +10309,19 @@ - + - + - + @@ -10348,13 +10382,13 @@ - + - + @@ -10373,17 +10407,17 @@ - + - + - + @@ -10403,11 +10437,11 @@ - - - + + + - + @@ -10417,7 +10451,7 @@ - + @@ -10504,31 +10538,31 @@ - - - + + + - + - - - + + + - - - + + + @@ -10541,7 +10575,7 @@ - + @@ -10952,36 +10986,36 @@ - - - - - - - - + + + + + + + + - - + + - - - - + + + + - + @@ -10993,23 +11027,23 @@ - - + + - - + + - + - + @@ -11027,35 +11061,35 @@ - + - - + + - - - + + + - + - + @@ -11081,14 +11115,14 @@ - - + + - + @@ -11108,44 +11142,44 @@ - - + + - - + + - - - - + + + + - - + + - + - + @@ -11155,19 +11189,19 @@ - + - + - + @@ -11205,13 +11239,13 @@ - + - + @@ -11226,9 +11260,9 @@ - - - + + + @@ -11303,56 +11337,56 @@ - - - - + + + + - + - + - + - + - - - + + + - + - - - - + + + + @@ -11362,7 +11396,7 @@ - + @@ -11388,7 +11422,7 @@ - + @@ -11411,63 +11445,63 @@ - - - - - - - - + + + + + + + + - + - + - + - - - - + + + + - - - - - - - - + + + + + + + + - + - + @@ -11489,39 +11523,39 @@ - + - + - + - - + + - + - + @@ -11540,24 +11574,24 @@ - - + + - - - + + + - - - + + + @@ -11568,40 +11602,40 @@ - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + - - - + + + - - - + + + @@ -11614,8 +11648,8 @@ - - + + @@ -11623,42 +11657,42 @@ - + - - - - + + + + - - + + - - - - - - + + + + + + - - - - + + + + @@ -11945,13 +11979,13 @@ - - - - - - - + + + + + + + @@ -11967,12 +12001,12 @@ - + - + @@ -12043,8 +12077,8 @@ - - + + @@ -12055,39 +12089,39 @@ - + - + - + - - + + - + - + @@ -12099,63 +12133,63 @@ - - + + - - - + + + - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + - - - + + + - - - + + + @@ -12163,40 +12197,40 @@ - + - - - - - - - - + + + + + + + + - - - - - - + + + + + + - - - - + + + + @@ -12448,20 +12482,20 @@ - - - - - - - + + + + + + + - - + + @@ -12515,8 +12549,8 @@ - - + + @@ -12525,23 +12559,23 @@ - + - - - + + + - - + + @@ -12551,48 +12585,48 @@ - - - + + + - - + + - - - + + + - - - - - - + + + + + + - - - - - - - - - - - + + + + + + + + + + + @@ -12732,8 +12766,8 @@ - - + + @@ -12749,15 +12783,15 @@ - - - - - - - - - + + + + + + + + + @@ -12801,15 +12835,15 @@ - - - - - - - - - + + + + + + + + + @@ -12874,60 +12908,60 @@ - + - + - - - + + + - - + + - - - + + + - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + @@ -13049,8 +13083,8 @@ - - + + @@ -13058,15 +13092,15 @@ - - - - - - - - - + + + + + + + + + @@ -13094,15 +13128,15 @@ - - - - - - - - - + + + + + + + + + @@ -13162,35 +13196,35 @@ - + - - - + + + - + - - + + - + @@ -13200,37 +13234,37 @@ - - - + + + - - + + - - - - - - + + + + + + - - - - - - - + + + + + + + @@ -13385,7 +13419,7 @@ - + @@ -13417,26 +13451,26 @@ - + - + - - - - - - - - - - + + + + + + + + + + @@ -13518,65 +13552,65 @@ - - - + + + - - - + + + - + - - + + - + - - - + + + - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + @@ -13714,7 +13748,7 @@ - + @@ -13726,18 +13760,18 @@ - - - - - - - - - - - - + + + + + + + + + + + + @@ -13812,66 +13846,66 @@ - - - + + + - - - + + + - - + + - - - + + + - + - - - + + + - - + + - + - - - - - + + + + + @@ -13924,56 +13958,56 @@ - - - - - + + + + + - - - + + + - - + + - - - + + + - - - - + + + + - - - - - - - - + + + + + + + + @@ -14020,45 +14054,45 @@ - + - + - + - + - + - + - + @@ -14068,11 +14102,11 @@ - - - - - + + + + + @@ -14080,21 +14114,21 @@ - - + + - + - + @@ -14106,37 +14140,37 @@ - - - + + + - - + + - - - - - - - + + + + + + + - - - - - - + + + + + + @@ -14203,45 +14237,45 @@ - - + + - - - - - - - + + + + + + + - - - - - - - - - + + + + + + + + + - - - - - - - - - + + + + + + + + + @@ -14307,45 +14341,45 @@ - + - + - + - + - + - + - + @@ -14355,11 +14389,11 @@ - - - - - + + + + + @@ -14367,21 +14401,21 @@ - - + + - + - + @@ -14389,29 +14423,29 @@ - - - + + + - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + @@ -14464,34 +14498,34 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + @@ -14548,50 +14582,50 @@ - + - - + + - + - + - - - - + + + + - - - - - - + + + + + + - + @@ -14601,7 +14635,7 @@ - + @@ -14610,12 +14644,12 @@ - - - - - - + + + + + + @@ -14628,7 +14662,7 @@ - + @@ -14676,27 +14710,27 @@ - - + + - - - - - - - - - - - - - - + + + + + + + + + + + + + + @@ -14736,17 +14770,17 @@ - - - + + + - + - + @@ -14807,168 +14841,168 @@ - - - + + + - - + + - - + + - + - + - + - - + + - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + - + - + - + - + - - + + - - - - - + + + + + - - - - - - - - - - - - - + + + + + + + + + + + + + - - + + - + - + - + - - - - + + + + - + - + @@ -14987,14 +15021,14 @@ - - + + - + @@ -15013,23 +15047,23 @@ - - - - - - + + + + + + - - - - - - + + + + + + @@ -15042,22 +15076,22 @@ - - - - - - - - - - + + + + + + + + + + - + @@ -15156,20 +15190,20 @@ - + - - - - - - + + + + + + - + @@ -15427,19 +15461,19 @@ - - + + - - - - - - + + + + + + @@ -15455,58 +15489,58 @@ - + - - - - - - - - - - - + + + + + + + + + + + - + - - - + + + - + - - + + - + @@ -15516,47 +15550,47 @@ - - - + + + - + - + - - + + - - - + + + - - - - - - + + + + + + @@ -16243,7 +16277,7 @@ - + @@ -16267,9 +16301,9 @@ - - - + + + @@ -16283,7 +16317,7 @@ - + @@ -16297,13 +16331,13 @@ - + - + @@ -16330,8 +16364,8 @@ - - + + @@ -16350,22 +16384,22 @@ - - - - - - + + + + + + - - - - - + + + + + @@ -16389,81 +16423,81 @@ - + - - + + - + - - + + - - - + + + - - - - + + + + - - - + + + - - - + + + - + - - - - + + + + - + @@ -16473,29 +16507,29 @@ - - - - + + + + - + - - - - - - - - + + + + + + + + @@ -16533,30 +16567,30 @@ - - - - - - - - + + + + + + + + - + - - - - + + + + @@ -16857,71 +16891,71 @@ - - - - + + + + - + - + - - - - - - + + + + + + - + - + - - - - - - + + + + + + - - - + + + - - - + + + @@ -16932,7 +16966,7 @@ - + @@ -16954,8 +16988,8 @@ - - + + @@ -16966,7 +17000,7 @@ - + @@ -16988,58 +17022,58 @@ - - + + - + - - + + - + - - - - + + + + - - - + + + - - + + - - + + @@ -17047,11 +17081,11 @@ - - - - - + + + + + @@ -17181,25 +17215,25 @@ - - - + + + - - - + + + - - + + @@ -17208,26 +17242,26 @@ - - + + - - - - + + + + - - - - + + + + @@ -17236,42 +17270,42 @@ - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + - + - - + + - - - - + + + + @@ -17288,29 +17322,29 @@ - + - + - + - - - - - - - - - - - + + + + + + + + + + + @@ -17321,8 +17355,8 @@ - - + + @@ -17366,69 +17400,80 @@ - - - - - - - - - + + + + + + + + + + + + + + + + + + + + - + - - + + - + - + - - - - - + + + + + - + - + - + - + - + - + @@ -17436,364 +17481,386 @@ - + - - - - - - - - - - - - + + + + + + + + + + + + - - + + - + - + + + + + + + - - - - - - - - - - - - - - + + - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - + + + + + + - - - - - - - - - + + + + + + + + + - + - + - + - + - - - - - - - - - - - - - - + + + + + + + + + + + + + + - - + + - - + + - - - + + + - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + - - - + + + - + - + - - + + - - + + - - - - - - - - - - - - - + + + + + + + + + + + + + - - - - - - - - - - - - - + + + + + + + + + + + + + - - - - - - - - - - - - + + + + + + + + + + + + - - - - - - - - - - - - - + + + + + + + + + + + + + - - - - - - - - - - - - - + + + + + + + + + + + + + - + - - + + - - + + - - - - - - - - - + + + + + + + + + - - - - + + + + - + - + - - - - + + + + - - - - - - + + + + + + - + @@ -17803,7 +17870,7 @@ - + @@ -17812,12 +17879,12 @@ - - - - - - + + + + + + @@ -17830,7 +17897,7 @@ - + @@ -17874,27 +17941,27 @@ - - + + - - - - - - - - - - - - - - + + + + + + + + + + + + + + @@ -17926,17 +17993,17 @@ - - - + + + - + - + @@ -17986,173 +18053,173 @@ - - - + + + - - + + - - + + - + - - - - - - - - - - - + + + + + + + + + + + - - + + - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + - + - + - + - + - - + + - - - - - + + + + + - - - - - - - - - - - - - + + + + + + + + + + + + + - - + + - + - + - + - - - - + + + + - + - + @@ -18171,14 +18238,14 @@ - - + + - + @@ -18197,23 +18264,23 @@ - - - - - - + + + + + + - - - - - - + + + + + + @@ -18226,22 +18293,22 @@ - - - - - - - - - - + + + + + + + + + + - + @@ -18340,20 +18407,20 @@ - + - - - - - - + + + + + + - + @@ -18599,28 +18666,28 @@ - - - - - - + + + + + + - - - + + + - - - - - - + + + + + + @@ -18636,58 +18703,58 @@ - + - - - - - - - - - - - + + + + + + + + + + + - + - - - + + + - + - - + + - + @@ -18697,47 +18764,47 @@ - - - + + + - + - + - - + + - - - + + + - - - - - - + + + + + + @@ -19445,7 +19512,7 @@ - + @@ -19469,9 +19536,9 @@ - - - + + + @@ -19485,7 +19552,7 @@ - + @@ -19499,13 +19566,13 @@ - + - + @@ -19532,8 +19599,8 @@ - - + + @@ -19544,31 +19611,31 @@ - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + - - - - - + + + + + @@ -19592,115 +19659,115 @@ - + - - + + - + - - + + - - - + + + - - - - + + + + - - - + + + - - - + + + - + - - - - + + + + - + - - - - - - - - - + + + + + + + + + - + - - - - - - - - + + + + + + + + @@ -19738,30 +19805,30 @@ - - - - - - - - + + + + + + + + - + - - - - + + + + @@ -20074,85 +20141,85 @@ - - - - - - - - + + + + + + + + - - - - + + + + - + - + - - - - - - + + + + + + - + - + - - - - - - + + + + + + - - - - - - - - - + + + + + + + + + - - - - + + + + @@ -20163,7 +20230,7 @@ - + @@ -20185,8 +20252,8 @@ - - + + @@ -20197,7 +20264,7 @@ - + @@ -20219,58 +20286,58 @@ - - + + - + - - + + - + - - - - + + + + - - - + + + - - + + - - + + @@ -20278,11 +20345,11 @@ - - - - - + + + + + @@ -20412,49 +20479,49 @@ - - - + + + - - - + + + - - + + - - + + - - - - + + + + - - - - + + + + @@ -20463,39 +20530,39 @@ - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + - - - - + + + + @@ -20512,29 +20579,29 @@ - + - + - + - - - - - - - - - - - + + + + + + + + + + + @@ -20543,8 +20610,8 @@ - - + + @@ -20582,406 +20649,441 @@ - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + - - - - + + + + - - - - - - - - - + + + + + + + + + - + - - + + - - - - - - - - - - - - + + + + + + + + + + + + - - - - - + + + + + - + - - - - - - - - - - - - - - + - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - + - + - + - - - - - - - - - - - - - - + + + + + + + + + + + + + + - - + + - - + + - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + - - - - - - - - + + + + + + + + - - + + - - + + - - - - - - - - - - - - - + + + + + + + + + + + + + - - - - - - - - - - - - - + + + + + + + + + + + + + - - - - - - - - - - - - + + + + + + + + + + + + - - - - - - - - - - - - - + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + - - + + - - + + - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + - + - + - + - + - + - + - - + + - + @@ -20992,7 +21094,7 @@ - + @@ -21115,46 +21217,46 @@ - - - + + + - + - - - - - - + + + + + + - + - + - - - - + + + + @@ -21171,7 +21273,7 @@ - + @@ -21181,28 +21283,28 @@ - + - + - + - + - + @@ -21257,27 +21359,27 @@ - - - + + + - - + + - - + + - + @@ -21288,7 +21390,7 @@ - + @@ -21401,57 +21503,57 @@ - - - + + + - - - - - - - - + + + + + + + + - + - + - - - - + + + + - + - - + + - - - + + + @@ -21796,25 +21898,25 @@ - + - + - + - + @@ -21822,9 +21924,9 @@ - - - + + + @@ -21901,21 +22003,21 @@ - + - - + + - + @@ -21968,13 +22070,13 @@ - + - + @@ -21982,16 +22084,16 @@ - - - - + + + + - + @@ -22002,20 +22104,20 @@ - + - + - - + + @@ -22040,85 +22142,85 @@ - - - - - - + + + + + + - - - - - - - - - + + + + + + + + + - - - + + + - + - + - + - + - + - - + + - - + + - - + + @@ -22129,38 +22231,40 @@ - - - - - - - + + + + + + + - - - - - - + + + + + + - - + + - - + + + + diff --git a/buildfile.m b/buildfile.m index b3007ae..82aa1a6 100644 --- a/buildfile.m +++ b/buildfile.m @@ -49,7 +49,6 @@ "Dependencies", "patch_"+[ext_nms([2 4 3])] ... ); - %% Testing and Artifacts % source files (folders) src = ["kern","src","utils"] + filesep; @@ -72,6 +71,8 @@ "Tag","syntax","SourceFiles",src ... ,CodeCoverageResults=fullfile("build","code-coverage-brief",["coverage.xml" "html/index.html"]) ... ,TestResults="build/test-results-brief/report.html" ... + ,Dependencies="compatability" ... + ,OutputDetail="Concise" ... ); %% Add compilation dependencies @@ -81,7 +82,7 @@ mfls = dir(fullfile(base, "src", "**", "msfm*.c")); [mfls, nms] = deal(string(fullfile({mfls.folder}, {mfls.name})), string({mfls.name})); -ofls = replace(fullfile(base, "bin", nms), '.c', "."+mexext()); +ofls = fullfile(base, "bin", replace(nms, ".c" + lineBoundary("end"), "."+mexext())); tnm = "compile_mex_"+extractBefore(nms,'.c'); for i = 1:numel(mfls) plan(tnm(i)) = matlab.buildtool.tasks.MexTask(mfls(i), fileparts(ofls(i)), ... @@ -417,3 +418,18 @@ function benchmarkTask(context) writelines(txt, ofl); end + +function compatabilityTask(context) +% Check for required toolboxes +req = ["Signal Processing", "Parallel Computing"]; +rec = ["Image Processing"]; +sug = reshape(string.empty,1,0); +rep = "The following packages are " + ["required", "recommended", "suggested"] ... + + ": " + sprintf('\t') + cellfun(@(s)join(s,", ",2),{req,rec,sug}); +arrayfun(@disp, rep(~ismissing(rep))); +vn = string({ver().Name}); +i = ismember(req + " Toolbox", vn); +assert(all(i), "The following toolboxes are required, but not installed:"+newline+join(req(~i)+ " Toolbox", newline)); +disp("All required toolboxes are installed."); + +end diff --git a/examples/beamforming/beamforming_walkthrough.mlx b/examples/beamforming/beamforming_walkthrough.mlx index 97c4643..90f8fd0 100644 Binary files a/examples/beamforming/beamforming_walkthrough.mlx and b/examples/beamforming/beamforming_walkthrough.mlx differ diff --git a/kern/pwznxcorr.m b/kern/pwznxcorr.m index 7a88c7f..4176f72 100644 --- a/kern/pwznxcorr.m +++ b/kern/pwznxcorr.m @@ -41,6 +41,18 @@ % y = PWZNXCORR(..., 'norm', true) normalizes the windowed vectors u and v % such that |u||_2 == ||v||_2 == 1. The default is true. % +% y = PWZNXCORR(..., 'tdim', tdim) uses dimension tdim as the time +% dimension. The default is 1. +% +% y = PWZNXCORR(..., 'ndim', ndim) uses dimension ndim as the channel +% dimension. The default is 2. +% +% y = PWZNXCORR(..., 'ldim', ldim) uses dimension ldim as the (output) lag +% dimension. The default is ndims(x) + 1. +% +% y = PWZNXCORR(..., 'lvec', true) vectorizes the lag dimension rather than +% computing each lag iteratively. The default is false. +% % y = PWZNXCORR(..., 'ref', 'neighbor') compares each channel n to the % neighboring channel n + 1. This is the default. % @@ -60,6 +72,10 @@ % end, end, end, end, % ``` % +% y = PWZNXCORR(..., 'ref', 'center', 'multi', true) compares each channel n +% to the central M == size(w, ndim) channels rather than the median channel. +% The default is false. +% % y = PWZNXCORR(..., 'ref', 'x0', 'x0', signal) compares each channel to % the given signal x0. The dimensions of x0 must be compatible with x. This % is roughly equivalent to modifying the psuedo-code above to contain the @@ -72,15 +88,6 @@ % end, end, end, end, % ``` % -% y = PWZNXCORR(..., 'tdim', tdim) uses dimension tdim as the time -% dimension. The default is 1. -% -% y = PWZNXCORR(..., 'ndim', ndim) uses dimension ndim as the channel -% dimension. The default is 2. -% -% y = PWZNXCORR(..., 'ldim', ldim) uses dimension ldim as the (output) lag -% dimension. The default is ndims(x) + 1. -% % y = PWZNXCORR(..., 'pad', false) disables 0-padding prior to computing % the cross-correlation. If `imfilter` (image processing toolbox) is not % available, this introduces wrap-around artefacts, but prevents the need @@ -130,13 +137,16 @@ kwargs.tdim (1,1) int64 {mustBeNumeric, mustBeInteger, mustBePositive} = 1 kwargs.ndim (1,1) int64 {mustBeNumeric, mustBeInteger, mustBePositive} = 2 kwargs.ldim (1,1) int64 {mustBeNumeric, mustBeInteger, mustBePositive} = ndims(x) + 1; + kwargs.multi (1,1) logical = false; % whether to reference multiple channels + kwargs.lvec (1,1) logical = true; % whether to vectorize lags end % alias x0 = kwargs.x0; % verify that w/W operates in the correct dimensions wsz = size(W,1:max(kwargs.tdim, kwargs.ndim)); -if any(wsz(setdiff(1:numel(wsz),[kwargs.tdim, kwargs.ndim])) ~= 1) +wsz_chk = wsz; wsz_chk([kwargs.tdim, kwargs.ndim]) = []; % size of W except in tdim,ndim +if any(wsz_chk ~= 1) error("QUPS:pwznxcorr:incompatibleWeightSize", ... "The filter weights w must be scalar in all dimensions except time (" ... +kwargs.tdim+") and channel ("+kwargs.ndim+")." ... @@ -155,6 +165,10 @@ [w, W] = deal(W, numel(W)); %#ok (for debug) end +% window sizing +[C, Wn] = deal(1, wsz(kwargs.ndim)); % number of channels: to correlate | normalization +if kwargs.multi, [C, Wn] = deal(Wn, C); end + % vector accumulation function isiflt = exist('imfilter', 'file'); % requires image processing toolbox if isiflt % whether to use imfilter @@ -187,11 +201,12 @@ case "center" % find the middle channel - mid = (N + 1) / 2; + mid = (N + 1 - C + 1)/2 + (0:C-1); n = unique([floor(mid), ceil(mid)]); % mean signal approach % n = floor(mid); % pick one signal, although spatially biased xl = x; % (T x N x ...) - xr = mean(sub(x, n, kwargs.ndim), kwargs.ndim); % (T x 1 x ...) + xr = sub(x, n, kwargs.ndim); + if ~kwargs.multi, xr = mean(xr, kwargs.ndim); end % (T x 1 x ...) case "x0" % TODO: validate data size / cast type etc. @@ -217,11 +232,13 @@ % windowed inner product for each lag % function y = windowed_inner_lag(lag) for i = numel(lags) : -1 : 1 - % get the lag for this iteration - lag = lags(i); - % delay right channel by lag - xr_l = conj(circshift(xr, -lag, tdim)); + if kwargs.lvec + xr_l = arrayfun(@(l){conj(circshift(xr, -l, tdim))}, lags); + xr_l = cat(kwargs.ldim, xr_l{:}); + else + xr_l = conj(circshift(xr, -lags(i), tdim)); + end % select portion if U > 1, xr_l = sub(xr_l, 1:U:T, tdim); end @@ -240,7 +257,7 @@ xrn_l = kernfun(real(xrz_l .* conj(xrz_l))); % normalization denominator - complex for (garbage) negative amplitudes - r = sqrt(complex(xln)) .* sqrt(complex(xrn_l)) .* sqrt(wsz(kwargs.ndim)) ; % denom + r = sqrt(complex(xln)) .* sqrt(complex(xrn_l)) .* sqrt(Wn); % denom % normalize y = y ./ r; @@ -248,11 +265,14 @@ % store yi{i} = y; + + % short-circuit + if kwargs.lvec, break; end end % compute for all lags and unpack % y = arrayfun(@windowed_inner_lag, lags, 'UniformOutput', false); -y = cast(cat(kwargs.ldim, yi{:}), 'like', x); % (T x [N|N-S] x ... x lags) +if ~kwargs.lvec, y = cast(cat(kwargs.ldim, yi{:}), 'like', x); end % (T x [N|N-S] x ... x lags) % undo manual 0-padding if kwargs.pad && ~isiflt, y = sub(y, 1:size(y,kwargs.tdim)-P, kwargs.tdim); end diff --git a/src/UltrasoundSystem.m b/src/UltrasoundSystem.m index 3664eb2..58f4250 100644 --- a/src/UltrasoundSystem.m +++ b/src/UltrasoundSystem.m @@ -698,7 +698,7 @@ function delete(us) k.GridSize = [ceil(QS ./ k.ThreadBlockSize(1)), N, M]; % get the computational bounds - sblk = (0 : k.GridSize - 1)' * k.ThreadBlockSize(1); % block starting time index + sblk = (0 : k.GridSize(1) - 1)' * k.ThreadBlockSize(1); % block starting time index blocks = sblk <= [-inf, sb(2,:), inf]; % point at which we are in bounds blocks = blocks(:,2:end) - blocks(:,1:end-1); % detect change sbk = cellfun(@(x) find(x,1,'first'), num2cell(blocks,2)); % get transition point @@ -4688,7 +4688,7 @@ function delete(us) % colormap gray; colorbar; title("Aperture growth apodization"); % % - % See also ULTRASOUNDSYSTEM/APACCEPTANCEANGLE + % See also APACCEPTANCEANGLE APCOSINEANGLE % defaults arguments @@ -4759,7 +4759,7 @@ function delete(us) % figure; imagesc(us.scan, ap); % animate(ap,'fs',1,'loop',false,"title","Angle: "+seq.angles); % - % See also apAcceptanceAngle + % See also APACCEPTANCEANGLE APCOSINEANGLE arguments us (1,1) UltrasoundSystem theta (1,:) = atan2d(us.seq.focus(1,:), us.seq.focus(3,:)); @@ -4817,7 +4817,7 @@ function delete(us) % nexttile(); imagesc(us.scan, bima, [-80 0] + max(bima(:))); % colormap gray; colorbar; title("Acceptance angle apodization"); % - % See also ULTRASOUNDSYSTEM/APAPERTUREGROWTH + % See also APAPERTUREGROWTH APCOSINEANGLE % defaults arguments @@ -4846,6 +4846,44 @@ function delete(us) % accept if greater than the cutoff angle apod = r >= cosd(theta); end + + + function apod = apCosineAngle(us, theta) + % APCOSINEANGLE - Create an cosine-weighted apodization array + % + % apod = APCOSINEANGLE(us) creates an ND-array to weight + % delayed data from the UltrasoundSystem us by the angle + % between the element normal and the element to pixel vector. + % The weighting is + % + % `cosd(min(90, (phi / theta)))` + % + % where phi is the angle and theta is the maximum angle. + % + % apod = APCOSINEANGLE(us, theta) uses an acceptance angle + % of theta in degrees. The default is 45. + % + % The output apod has dimensions I1 x I2 x I3 x N x 1 where + % I1 x I2 x I3 are the dimensions of the scan, N is the number + % of receive elements. + % + % See also APAPERTUREGROWTH APACCEPTANCEANGLE + + % defaults + arguments + us (1,1) UltrasoundSystem + theta (1,1) {mustBePositive} = 45 + end + + % cosine apodization + [pg, pn] = deal(us.scan.positions(), us.xdc.positions()); % pixels | elems + [~,~,nn] = us.xdc.orientations(); % elem normals + [pn, nn] = deal(swapdim(pn,2,5), swapdim(nn,2,5)); % match dims + r = (pg - pn); % elem -> pixel vector + r = r ./ vecnorm(r,2,1); % normalized + r = swapdim(pagemtimes(nn, 'transpose', r, 'none'), 2:6); % normalized inner product + apod = cosd(min(90,(90/theta)*acosd(r))); % cosine gradient with (scaled) angle + end end % dependent methods @@ -4943,10 +4981,10 @@ function delete(us) end end - function com = recompileCUDA(us, defs, arch) + function com = recompileCUDA(us, defs, arch, kwargs) % RECOMPILECUDA - Recompile CUDA ptx files % - % RECOMPILECUDA(us) recompiles all CUDA files and stores + % RECOMPILECUDA(us) recompiles all CUDA files to ptx and stores % them in us.tmp_folder. % % RECOMPILECUDA(us, defs) compiles for the compiler @@ -4963,6 +5001,12 @@ function delete(us) % If arch is an empty string, no architecture argument is used. % The default is the architecture of the current GPU device. % + % RECOMPILECUDA(..., 'mex', true) uses `mexcuda` rather than + % `nvcc` to compile ptx files (requires R2023a or later). + % Architecture selection, warning selection and suppression, + % and other compiler options are not available via mexcuda. The + % default is true if `nvcc` is not found on the path. + % % nvcom = RECOMPILECUDA(...) returns the command sent to nvcc. % This command can be tweaked and resent to debug or fix % compilation errors. @@ -4975,7 +5019,8 @@ function delete(us) % % change the source code ... % % % recompile the 3rd CUDA file manually - % system(nvcom(3)); + % system(nvcom(3)); % via nvcc + % % mexcuda(nvcom{3}{:}); % via mexcuda % % See also ULTRASOUNDSYSTEM.RECOMPILE ULTRASOUNDSYSTEM.RECOMPILEMEX % ULTRASOUNDSYSTEM.GENCUDADEFS @@ -4984,12 +5029,17 @@ function delete(us) us (1,1) UltrasoundSystem defs (1,:) struct = UltrasoundSystem.genCUDAdefs(); arch (:,1) string {mustBeScalarOrEmpty, mustBeArch(arch)} = nvarch() + kwargs.mex (1,1) logical = isempty(argn(2, @system, 'which nvcc')) && ~isMATLABReleaseOlderThan("R2023a") end + + % ensure that the output folder exists + if ~exist(us.tmp_folder,'dir'), warning("QUPS:UltrasoundSystem:ReinitializingTmpFolder","Recreating a tmp_folder that was deleted (somehow)."); mkdir(us.tmp_folder); addpath(us.tmp_folder); end % compile each for i = numel(defs):-1:1 d = defs(i); % make full command + if ~kwargs.mex % via nvcc com(i) = join(cat(1,... "nvcc ", ... "--ptx " + which(string(d.Source)), ... @@ -5001,8 +5051,22 @@ function delete(us) join("-W" + d.Warnings), ... join("-D" + d.DefinedMacros)... )); - - try s = system(com(i)); + comp = @system; + else % via mexcuda + com{i} = cellstr(cat(1,... + "-ptx", ... + "-output", fullfile(us.tmp_folder, argn(2, @fileparts, d.Source) + ".ptx"), ... + ... ("--" + d.CompileOptions),... + ("-I" + d.IncludePath), ... + ("-L" + d.Libraries), ... + ... ("-W" + d.Warnings), ... + ("-D" + d.DefinedMacros),... + which(string(d.Source)) ... + ... "-arch=" + arch + " ", ... compile for active gpu + )); + comp = @(s) mexcuda(s{1}{:}); + end + try s = comp(com(i)); if s, warning( ... "QUPS:recompile:UnableToRecompile", ... "Unable to recompile " + d.Source ... diff --git a/src/interpd.cu b/src/interpd.cu index 3bdc0c3..07be7b6 100644 --- a/src/interpd.cu +++ b/src/interpd.cu @@ -93,17 +93,17 @@ __device__ T2 cubic(const T2 * x, U tau, T2 no_v) { if (!(0 <= (ti - 1) && (ti + 2) < QUPS_T)) return no_v; - T2 s0 = x[ti - 1]; - T2 s1 = x[ti + 0]; - T2 s2 = x[ti + 1]; - T2 s3 = x[ti + 2]; + const T2 s0 = x[ti - 1]; + const T2 s1 = x[ti + 0]; + const T2 s2 = x[ti + 1]; + const T2 s3 = x[ti + 2]; // Cubic Hermite interpolation (increased precision using fused multiply-adds) // (Catmull-Rom) - U a0 = 0 + u * (-1 + u * (+2 * u - 1)); - U a1 = 2 + u * (+0 + u * (-5 * u + 3)); - U a2 = 0 + u * (+1 + u * (+4 * u - 3)); - U a3 = 0 + u * (+0 + u * (-1 * u + 1)); + const U a0 = 0.f + u * (-1.f + u * (+2.f * u - 1.f)); + const U a1 = 2.f + u * (+0.f + u * (-5.f * u + 3.f)); + const U a2 = 0.f + u * (+1.f + u * (+4.f * u - 3.f)); + const U a3 = 0.f + u * (+0.f + u * (-1.f * u + 1.f)); // // Cubic Hermite interpolation (naive, less precise implementation) // float a0 = -1 * u * u * u + 2 * u * u - 1 * u + 0; // float a1 = +3 * u * u * u - 5 * u * u + 0 * u + 2; diff --git a/test/BFTest.m b/test/BFTest.m index 0541dcb..d8672b4 100644 --- a/test/BFTest.m +++ b/test/BFTest.m @@ -209,7 +209,7 @@ function reselectgpu() % Github test routine - methods(Test, ParameterCombination = 'sequential', TestTags={'Github'}) + methods(Test, ParameterCombination = 'sequential', TestTags={'Github', 'build'}) function github_psf(test, bf_name)%, prec, terp) if(any(bf_name == ["Eikonal","Adjoint"])), return; end % Too much compute @@ -247,9 +247,9 @@ function psf(test, gdev, bf_name, prec, terp, apodization) % for the Eikonal beamformer, pass if not given FSA delays if(bf_name == "Eikonal" && us.seq.type ~= "FSA"), return; end % if using a frequency domain method, skip half precision - the phase errors are too large - if(ismember(bf_name, ["Adjoint"]) && prec == "halfT"), return; end + if(ismember(bf_name, "Adjoint") && prec == "halfT"), return; end % weird phase error, but it's not important right now - skip it - if(ismember(bf_name, ["Adjoint"]) && us.seq.type == "VS") return; end % && isa(us.xdc, 'TransducerConvex')); + if(ismember(bf_name, "Adjoint") && us.seq.type == "VS"), return; end % && isa(us.xdc, 'TransducerConvex')); % is using the adjoint method, pagemtimes,pagetranspose must be supported test.assumeTrue( bf_name ~= "Adjoint" || (... logical(exist('pagemtimes' , 'builtin')) ... @@ -264,7 +264,7 @@ function psf(test, gdev, bf_name, prec, terp, apodization) % for frequency domain methods, the time-axis must be extended % so that the replications of the miage are outside of the % imaging range - if ismember(bf_name, ["Adjoint"]) + if ismember(bf_name, "Adjoint") dr = hypot(range(scanc.xb), range(scanc.zb)) / 2; % half the largest range across the image dT = round(2 * dr / scat.c0 * chd.fs); % temporal buffer in indices chd = zeropad(chd, dT, dT); % add buffer on both sides diff --git a/test/ChdTest.m b/test/ChdTest.m index 2a0fef0..5fa8c0e 100644 --- a/test/ChdTest.m +++ b/test/ChdTest.m @@ -39,7 +39,7 @@ function typeCheck(tst) % tall, sparse chdo = tall(chd); - sparse(ChannelData('data', randi([0 1], [16 32]))) + sparse(ChannelData('data', randi([0 1], [16 32]))); % unit conversion check chdo = angle(complex(chd)); @@ -75,9 +75,9 @@ function freqDomainCheck(tst) D3 = chd.getLowpassFilter(0.2); D4 = chd.getPassbandFilter([0.1 0.4]); for D = [D1 D2] - filter( chd, D) - filtfilt(chd, D) - fftfilt( chd, D) + filter( chd, D); + filtfilt(chd, D); + fftfilt( chd, D); end % sampling @@ -129,13 +129,13 @@ function transforming(tst) tst.assertEqual(-wvr.tend, wv.t0); tst.assertEqual(-wv.tend, wvr.t0); tst.assertEqual(wv.samples, reverse(reverse(wv)).samples); - for s = ["full", "same", "valid"], convt(chd, wv, s), end + for s = ["full", "same", "valid"], convt(chd, wv, s); end % join/splice - ChannelData.empty().join(4) - join(chd.splice(),4) - chd.splice(4) - chd.splice(1,4) + ChannelData.empty().join(4); + join(chd.splice(),4); + chd.splice(4); + chd.splice(1,4); % data order ord = [1,4,2,3]; diff --git a/test/ExampleTest.m b/test/ExampleTest.m index 2c662e3..be7610a 100644 --- a/test/ExampleTest.m +++ b/test/ExampleTest.m @@ -7,6 +7,9 @@ scat_lim (1,1) double % limit of scatterers pulse_lim (1,2) double % max number of pulses to sim in k-Wave, greens end + properties(Constant) + debug (1,1) logical = false; % whether to print debug warnings + end % blacklists properties @@ -56,9 +59,9 @@ 'FieldII', ["calc_scat"+["","_all", "_multi"], "getFieldII"+["Aperture","Patches","Positions"]], ... FieldII 'MUST', "simus", ... MUST 'fullwave', ["getFullwaveTransducer", "fullwaveSim", "fullwaveConf", "fullwaveJob", "mapToCoords"], ... % fullwave - 'mex', ["bfEikonal", "msfm"+["","2d","3d"]], ... mex binaries required ... compilation (optional, requires CUDA) + 'mex', "recompile" + ["","Mex"], ... mex binaries required ... compilation (optional, requires CUDA) 'gpu', ["wbilerpg","feval"], ... CUDAKernel or oclKernel support required - 'comp', ("recompile" + ["","Mex","CUDA"]), ... compilations setup required + 'comp', ("recompile" + ["","CUDA"]), ... compilations setup required 'ext', "vol3d", ... external functions }, 2, [])'; % + "("; bl_fcn(:,2) = cellfun(@(s){s+"("}, bl_fcn(:,2)); % append "(" @@ -93,12 +96,14 @@ if all(arrayfun(@(k) exist(fullfile(prj_rt, "bin", k), 'file'), kerns)) bl_fcn(bl_fcn(:,1) == "gpu",:) = []; end - kerns = ["msfm2d", "msfm3d"]+"."+mexext; % require mex binaries - if 1||all(arrayfun(@(k) exist(fullfile(prj_rt, "bin", k), 'file'), kerns)) + kerns = ["msfm2d", "msfm3d"]+".c"; % require mex source files + lst = arrayfun(@(k) dir(fullfile(prj_rt, "**", k)), kerns, 'uni', 0); + if all(ismember(kerns, {vertcat(lst{:}).name})) + % if all(arrayfun(@(k) exist(fullfile(prj_rt, "bin", k), 'file'), kerns)) bl_fcn(bl_fcn(:,1) == "mex",:) = []; end % fullwave executable must exist in 'bin' folder under this name - kerns = ["fullwave2_executable"]; + kerns = "fullwave2_executable"; if all(arrayfun(@(k) exist(fullfile(prj_rt, "bin", k), 'file'), kerns)) ... && exist("mapToCoords", "file") bl_fcn(bl_fcn(:,1) == "fullwave",:) = []; @@ -136,22 +141,22 @@ function filterBlacklist(test, code, blk, fls) ,"Matching packages: " + join(string(pkg(1,p)), ", ") ... ]), newline)); end - function silenceAcceptableWarnings(testCase) + function silenceAcceptableWarnings(test) lids = [... "QUPS:kspaceFirstOrder:upsampling", ... % in US "MATLAB:ver:ProductNameDeprecated", ... % in k-Wave "QUPS:recompile:UnableToRecompile" ... % in US ]; W = warning(); % get current state - arrayfun(@(l) warning( 'off', l), lids); % silence - testCase.addTeardown(@() warning(W)); % restore on exit + arrayfun(@(l) warning('off', l), lids); % silence + test.addTeardown(@() warning(W)); % restore on exit if ~isempty(gcp('nocreate')) % any pool - execute on each worker ws = parfevalOnAll(@warning, 1); wait(ws);% current state ws = fetchOutputs(ws); % retrieve [~, i] = unique(string({ws.identifier}), 'stable'); ws = ws(i); % select first unique set (assumer identical) wait(arrayfun(@(l) parfevalOnAll(@warning, 0, 'off', l), lids)); % silence - testCase.addTeardown(@() parfevalOnAll(@warning, 0, ws)); % restore on exit + test.addTeardown(@() parfevalOnAll(@warning, 0, ws)); % restore on exit end end function txt = truncateJobs(test, txt, kwargs) @@ -286,14 +291,14 @@ function cleanup_test(test) if isempty(bdln) bdln = nan; elseif ~isscalar(bdln) - warning("Multiple lines match " + string(pt(j)) + " in file " + fls_(m) + ": choosing the "+ord(j)+"."); + if ExampleTest.debug, warning("Multiple lines match " + string(pt(j)) + " in file " + fls_(m) + ": choosing the "+ord(j)+"."); end switch j, case 1, bdln = bdln(1); case 2, bdln = bdln(end); end end k(j) = bdln; %#ok end k = k + [1, -1]; % go after/before matching start/end if isnan(k(1)) - % warning("No example found for lines "+join(string(blk{:}([1 end])),"-")+" in file "+f+"."); + % if ExampleTest.debug, warning("No example found for lines "+join(string(blk{:}([1 end])),"-")+" in file "+f+"."); end continue; end if isnan(k(2)), k(2) = length(code_); end % can't find an ending - assume it's okay @@ -330,7 +335,7 @@ function run_main_example(test) test.assumeTrue(logical(exist("kspaceFirstOrder3D", "file")), "Missing package kWave" ); % kWave for f = fnms copyfile(which(f+".m"), "./"); % copy to here (temp folder) - eval(f); % run + evalc(f); % run end end end @@ -343,7 +348,7 @@ function run_kernel_examples(test, fls, lns) else base_dir = prj.RootFolder; end - flds = fullfile(base_dir, ["examples"]); % blacklist classes and examples + flds = fullfile(base_dir, "examples"); % blacklist classes and examples run_examples(test, fls, lns, flds); end end @@ -418,9 +423,9 @@ function run_examples(test, fls, lns, flt_fld) % assert a clean run test.silenceAcceptableWarnings(); if any(contains(code, "compile")) - f = str2func("@"+fnm); f(); % run directly + evalc(fnm); % run directly else - test.assertWarningFree(str2func("@"+fnm), "Example "+fnm+" did not complete without a warning!"); + test.assertWarningFree(@()evalc(fnm), "Example "+fnm+" did not complete without a warning!"); end end end @@ -437,7 +442,7 @@ function run_examples(test, fls, lns, flt_fld) i = find(contains(code,pat)); % lines with call to kspaceFirstOrder if isempty(i), V = 0; return; end % 0 if greens not found -if ~isscalar(i), i = i(1); warning("Multiple calls to "+fnm+"(): choosing the first."); end +if ~isscalar(i), i = i(1); if ExampleTest.debug, warning("Multiple calls to "+fnm+"(): choosing the first."); end, end v = strip(extractBetween(code(i), pat, ("," | ")"))); % us variable % write code to a file @@ -467,3 +472,4 @@ function run_examples(test, fls, lns, flt_fld) fclose(fid); end end + diff --git a/test/InitTest.m b/test/InitTest.m index 617d9db..f84c6f0 100644 --- a/test/InitTest.m +++ b/test/InitTest.m @@ -16,14 +16,6 @@ ]); end - properties - end - - methods(TestClassSetup, ParameterCombination = 'exhaustive') - end - methods(TestClassTeardown) - end - methods(TestMethodSetup) % Setup for each test function fig(~), figure; end @@ -31,6 +23,24 @@ methods(TestMethodTeardown) function cls(~), close; end end + methods(TestClassSetup) + function silenceAcceptableWarnings(tst) + lids = ... + "MATLAB:structOnObject" ... calling struct directly + ; + W = warning(); % get current state + arrayfun(@(l) warning('off', l), lids); % silence + tst.addTeardown(@() warning(W)); % restore on exit + if ~isempty(gcp('nocreate')) % any pool - execute on each worker + ws = parfevalOnAll(@warning, 1); wait(ws);% current state + ws = fetchOutputs(ws); % retrieve + [~, i] = unique(string({ws.identifier}), 'stable'); + ws = ws(i); % select first unique set (assumer identical) + wait(arrayfun(@(l) parfevalOnAll(@warning, 0, 'off', l), lids)); % silence + tst.addTeardown(@() parfevalOnAll(@warning, 0, ws)); % restore on exit + end + end + end methods(Test) function initxdc(test) % INITXDC - Assert that Transducer constructors initialize @@ -43,10 +53,10 @@ function initxdc(test) arrayfun(@plot, xdcs); % supports plotting arrayfun(@patch, xdcs); % supports patch arrayfun(@(x)patch(x, nexttile(), 'el_sub_div', [2 2]), xdcs); % supports args - xdcs(end+2) = xdcs(1), %#ok implicit empty value, display - arrayfun(@disp, xdcs) % display scalar - arrayfun(@(x)disp(x([])), xdcs) % display empty - x = copy(xdcs); x.delete(); x, arrayfun(@disp, x); % display deleted + xdcs(end+2) = xdcs(1); % implicit empty value + arrayfun(@edisp, xdcs) % display scalar + arrayfun(@(x)edisp(x([])), xdcs) % display empty + x = copy(xdcs); x.delete(); arrayfun(@edisp, x); % display deleted test.assertWarning(@()xdcs(1).ultrasoundTransducerImpulse(), "QUPS:Transducer:DeprecatedMethod"); [xdcs.origin] = deal([0 0 -10e-3]); % offset (to be deprecated?) [xdcs.rot] = deal([20 -10]); % offset (to be deprecated?) @@ -65,15 +75,15 @@ function initseq(test) [seqs([seqs.type] == "FSA").numPulse] = deal(1); arrayfun(@(seq) test.assertThat(seq, IsInstanceOf('Sequence')), seqs); arrayfun(@(scn) scale(scn, 'dist', 1e3), seqs, "UniformOutput",false); % can scale - seqs(end+2) = seqs(1), %#ok implicit empty value, display + seqs(end+2) = seqs(1); % implicit empty value seqs(end-1).numPulse = 1; % fix implicit seq s = arrayfun(@obj2struct, seqs, 'UniformOutput', false); % supports specialized struct conversion arrayfun(@plot, seqs, 'UniformOutput',false); % supports plotting cellfun(@(s) test.assertThat(s.pulse, IsInstanceOf('struct')), s); % recursive check - arrayfun(@disp, seqs) % display scalar - arrayfun(@(x)disp(x([])), seqs) % display empty + arrayfun(@edisp, seqs) % display scalar + arrayfun(@(x)edisp(x([])), seqs) % display empty arrayfun(@splice, seqs, 'UniformOutput', false); % can splice - x = copy(seqs); x.delete(); x, arrayfun(@disp, x); % display deleted + x = copy(seqs); x.delete(); arrayfun(@edisp, x); % display deleted % polar: manipulate range/angles seq = SequenceRadial("focus",[0 0 1]'); % fine @@ -82,7 +92,7 @@ function initseq(test) seq = SequenceRadial("ranges",[1 1 1]); % fine seq = SequenceRadial("angles",-1:1, "ranges",[1 1 1]); % fine seq = SequenceRadial("angles", 0, "ranges",[1 1 1]); % fine - seq = SequenceRadial("angles",-1:1, "ranges",[1]); % fine + seq = SequenceRadial("angles",-1:1, "ranges",1); % fine seq.angles = [-5 0 5]; % fine seq.ranges = [2 3 2]; % fine test.assertError(@()SequenceRadial("focus",[0 0 1]',"angles",0),""); % bad @@ -114,10 +124,10 @@ function initscan(test) arrayfun(@obj2struct, scns, 'UniformOutput', false); % supports specialized struct conversion arrayfun(@plot, scns); % supports plotting arrayfun(@(s) imagesc(s,randn(s.size)), scns); % supports imagesc - scns(end+2) = scns(1),%#ok implicit empty value, display - arrayfun(@disp, scns) % display scalar - arrayfun(@(x)disp(x([])), scns) % display empty - x = copy(scns); x.delete(); x, arrayfun(@disp, x); % display deleted + scns(end+2) = scns(1); % implicit empty value + arrayfun(@edisp, scns) % display scalar + arrayfun(@(x)edisp(x([])), scns) % display empty + x = copy(scns); x.delete(); arrayfun(@edisp, x); % display deleted % supports req'd overload methods [~,~,~] = arrayfun(@getImagingGrid, scns, "UniformOutput",false); @@ -144,8 +154,8 @@ function initscan(test) end % convert polar to cart - ScanCartesian(ScanPolar()), - ScanCartesian(ScanSpherical()), + ScanCartesian(ScanPolar()); + ScanCartesian(ScanSpherical()); % deprecations test.assertWarning(@()scanCartesian(ScanPolar() ), "QUPS:ScanPolar:syntaxDeprecated" ); @@ -155,7 +165,7 @@ function initscan(test) function initwv(test) % INITWV - Assert that the Waveform constructor initializes import matlab.unittest.constraints.IsInstanceOf; - wv = Waveform(); % init + wv = Waveform(); edisp(wv); % construct & display test.assertThat(wv, IsInstanceOf('Waveform')); % multiple ways to init @@ -186,7 +196,7 @@ function initchd(test) % INITCHD - Assert that a ChannelData constructor initializes % without arguments import matlab.unittest.constraints.IsInstanceOf; - chd = ChannelData('data',0), %#ok implicit display + chd = ChannelData('data',0); test.assertThat(chd, IsInstanceOf('ChannelData')); scale(chd, 'time', 1e6); arrayfun(@obj2struct, chd, 'UniformOutput', false); % supports specialized struct conversion @@ -200,19 +210,20 @@ function initus(test) % INITUS - Assert that an UltrasoundSystem constructor % initializes without arguments import matlab.unittest.constraints.IsInstanceOf; - us = UltrasoundSystem(), %#ok implicit display + us = UltrasoundSystem(); edisp(us); % construct & display test.assertThat(us, IsInstanceOf('UltrasoundSystem')); fld = us.tmp_folder; test.assertTrue(logical(exist(fld,"dir"))); UltrasoundSystem('copybin', true); - UltrasoundSystem('recompile',true); + evalc("UltrasoundSystem('recompile',true);"); us2 = copy(us); fld2 = us2.tmp_folder; % new instance test.assertFalse(fld == fld2); % different folders clear us2; test.assertFalse(logical(exist(fld2,'dir'))); % deleted on clear us2 = copy(us); fld2 = us2.tmp_folder; % new instance copyfile(which("Qups.prj"), fld2); % dirty test.addTeardown(@()delete(fullfile(fld2,"Qups.prj"))); % cleanup + try clear us2 % should error [~, lid] = lastwarn; test.assertTrue(lid == "MATLAB:class:DestructorError", ... @@ -242,8 +253,8 @@ function initus(test) uss = copy([us us us]); % array uss(end+2) = copy(us); % implicit creation - arrayfun(@disp, uss) % display scalar - arrayfun(@(x)disp(x([])), uss) % display empty + arrayfun(@edisp, uss) % display scalar + arrayfun(@(x)edisp(x([])), uss) % display empty % deprecation test.assertWarning(@()setfield(us,'sequence',Sequence()),"QUPS:UltrasoundSystem:syntaxDeprecated") @@ -257,8 +268,8 @@ function initmedscat(test) % initialize without arguments import matlab.unittest.constraints.IsInstanceOf; - sct = Scatterers(), %#ok implicit display - med = Medium(), %#ok implicit display + sct = Scatterers(); edisp(sct); % construct & display + med = Medium(); edisp(med); % construct & display test.assertThat(sct, IsInstanceOf('Scatterers')); test.assertThat(med, IsInstanceOf('Medium' )); scale(sct, "dist", 1e3, "time", 1e6); @@ -270,9 +281,9 @@ function initmedscat(test) % special constructors grd = ScanCartesian(); - Medium.Sampled(grd) - Scatterers.Diffuse(grd) - Scatterers.Grid() + Medium.Sampled(grd); + Scatterers.Diffuse(grd); + Scatterers.Grid(); % deprecated properties test.assertWarning(@()setfield(med,'alphap0',1.2), "QUPS:Medium:syntaxDeprecated") @@ -337,4 +348,6 @@ function propTest(test, clss) end end -end \ No newline at end of file +end + +function edisp(x), evalc("disp(x)"); end %#ok diff --git a/test/InteropTest.m b/test/InteropTest.m index cf6c759..bce14be 100644 --- a/test/InteropTest.m +++ b/test/InteropTest.m @@ -74,13 +74,32 @@ if isempty(vsx_files), vsx_files = {''}; end end end + methods(TestClassSetup) + function silenceAcceptableWarnings(tst) + lids = [... + "QUPS:Sequence:DeprecatedValue" ... in Sequence.UFF + "QUPS:QUPS2USTB:ambiguousSequence" ... in Sequence.QUPS2USTB + ]; + W = warning(); % get current state + arrayfun(@(l) warning('off', l), lids); % silence + tst.addTeardown(@() warning(W)); % restore on exit + if ~isempty(gcp('nocreate')) % any pool - execute on each worker + ws = parfevalOnAll(@warning, 1); wait(ws);% current state + ws = fetchOutputs(ws); % retrieve + [~, i] = unique(string({ws.identifier}), 'stable'); + ws = ws(i); % select first unique set (assumer identical) + wait(arrayfun(@(l) parfevalOnAll(@warning, 0, 'off', l), lids)); % silence + tst.addTeardown(@() parfevalOnAll(@warning, 0, ws)); % restore on exit + end + end + end methods(Test) function fieldII_xdc(tst, xdcs) tst.assumeTrue(logical(exist('field_init','file'))) % need FieldII for this import matlab.unittest.constraints.IsEqualTo; import matlab.unittest.constraints.AbsoluteTolerance; - xdcs.patches() % assumes defaults + xdcs.patches(); % assumes defaults pchq = xdcs.patches([4 2]); pchf = xdcs.getFieldIIPatches([4 2]); [pchq, pchf] = dealfun(@(x) cell2mat(swapdim(cellfun(@(x) {cat(3, x{:})}, x),1:2,4:5)), pchq, pchf); diff --git a/test/KernTest.m b/test/KernTest.m index 1bcd5ce..29206f5 100644 --- a/test/KernTest.m +++ b/test/KernTest.m @@ -62,6 +62,24 @@ function set_data(test) b = complex(randn(sz), randn(sz)); test.dargs{3} = {us, b}; end + + function silenceAcceptableWarnings(test) + lids = ... + "QUPS:animate:MatchingPermutation" ... in animate + ; + W = warning(); % get current state + arrayfun(@(l) warning('off', l), lids); % silence + test.addTeardown(@() warning(W)); % restore on exit + if ~isempty(gcp('nocreate')) % any pool - execute on each worker + ws = parfevalOnAll(@warning, 1); wait(ws);% current state + ws = fetchOutputs(ws); % retrieve + [~, i] = unique(string({ws.identifier}), 'stable'); + ws = ws(i); % select first unique set (assumer identical) + wait(arrayfun(@(l) parfevalOnAll(@warning, 0, 'off', l), lids)); % silence + test.addTeardown(@() parfevalOnAll(@warning, 0, ws)); % restore on exit + end + end + end methods function set_device(test, dev), test.setDev(dev); end end diff --git a/test/interpTest.m b/test/interpTest.m index a55053c..27e3172 100644 --- a/test/interpTest.m +++ b/test/interpTest.m @@ -47,7 +47,7 @@ function setupData(test, dsize,type, dev) tau_ = cfun(tau_); % gpu/cpu - switch dev, + switch dev case "CPU", x0_ = gather( x0_); case "GPU", x0_ = gpuArray(x0_); end @@ -90,7 +90,7 @@ function setDev(test, dev) % all permutations to test % ord = {[1,2,3,4], [1,3,2,4], [1,4,2,3]} % permute the data? terp = {'cubic', 'nearest', 'linear'}; % overlapping methods - dsum = {[], [2], [3,4], 6} % summation dimensions + dsum = {[], 2, [3,4], 6} % summation dimensions wvecd = {[], 3, [3,4]} % different outer-product weights dimensions end methods(Test, ParameterCombination = 'exhaustive') diff --git a/utils/animate.m b/utils/animate.m index 8e75866..d4310df 100644 --- a/utils/animate.m +++ b/utils/animate.m @@ -111,7 +111,7 @@ if r(1) == c, r = r(2); else, r = r(1); end % extract row dim if all([c r] == [1 2]), continue; end % already in order ord = [c r, setdiff(dvec, [c,r])]; % permutation order - put '[c r]' first - warning("Permuting "+ith(i)+" argument to match the image (ord = ["+join(string(ord),",")+"])."); + warning("QUPS:animate:MatchingPermutation","Permuting "+ith(i)+" argument to match the image (ord = ["+join(string(ord),",")+"])."); x{i} = permute(x{i}, ord); % permute end diff --git a/utils/swapdim.m b/utils/swapdim.m index f975f18..0d5c498 100644 --- a/utils/swapdim.m +++ b/utils/swapdim.m @@ -31,12 +31,12 @@ ord0 = ord; % get permutation -if isempty(intersect(i, o)) % can be handled separately +if (noIntersect(i, o)) % can be handled separately ord([i o]) = [o i]; % swap dimensions else % must fill in missing indices l = min(min(i),min(o)) : max(max(i),max(o)); % all indices within swap - i = [i, setdiff(l, i)]; % expanded input indices - o = [o, setdiff(l, o)]; % expanded output indices + i = appendDiff(i,l); % expanded input indices + o = appendDiff(o,l); % expanded output indices ord(o) = i; % full permutation ordering end @@ -53,3 +53,8 @@ x = permute(x, ord); end +function tf = noIntersect(i,o), tf = all(i'~=o,'all'); +% replaces `isempty(intersect(i, o))` + +function i = appendDiff(i, l), i = [i,l(all(l~=i',1))]; +% replaces `i = [i, setdiff(l, i)];`