v1.1.1 Changes

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.0
+version: v1.1.1
 license: Apache-2.0
 date-released: 2024-04-07
 repository-code: "https://github.com/thorstone25/qups"

buildfile.m

@@ -389,6 +389,20 @@ function patch_kWaveTask(context)
     writelines(txt, fl); % save modified file
 end
 
+% if recording pressure vs. time on gpu, keep storage on CPU when
+% allocating memory
+fl = dir(fullfile(fld, '**', "kspaceFirstOrder_createStorageVariables.m"));
+fl = string(fullfile({fl.folder}, {fl.name}));
+pat = "        sensor_data.p = castZeros([num_sensor_points, num_recorded_time_points]);";
+rep = join([
+    "        sensor_data.p = castZeros([num_sensor_points, 1]);", ...
+    "        if isa(sensor_data.p, 'gpuArray'), sensor_data.p = gather(sensor_data.p); end % always retain on CPU", ...
+    "        sensor_data.p = repmat(sensor_data.p, [1, num_recorded_time_points]);", ...
+    ], newline);
+txt = string(readlines(fl));
+txt = replace(txt, pat, rep); % find & replace
+writelines(txt, fl); % save modified file
+
 end
 
 function benchmarkTask(context)

kern/convd.m

@@ -1,14 +1,15 @@
 function [C, lags] = convd(x, y, dim, shape, kwargs)
 % CONVD - GPU-enabled Convolution in one dimension
 %
-% C = CONVD(A, B) computes the convolution of A with B in dimension 1. A 
-% and B may differ in size in dimension 1 only. All other dimensions must
-% be of equal size.
+% C = CONVD(A, B) computes the convolution of A with B in the first
+% non-singleton dimension. A and B must have compatible dimensions.
 % 
-% C = CONVD(A, B, dim) executes in dimension "dim" instead of dimension 1
+% C = CONVD(A) computes the auto-correlation i.e. B = conj(flip(A));
+% 
+% C = CONVD(A, B, dim) executes in dimension dim instead of dimension 1
 % 
 % C = CONVD(A, B, dim, shape) selects the shape of the output. Must be one 
-% of {'full'*|'same'|'valid'}. The default is 'full'.
+% of {'full', 'same', 'valid'}. The default is 'full'.
 % 
 % C = CONVD(..., 'gpu', true) selects whether to use a gpu. A ptx-file will 
 % be used if compiled. The default is true if x or y is a gpuArray.
@@ -52,7 +53,7 @@
 
 arguments
     x {mustBeFloat}
-    y {mustBeFloat}
+    y {mustBeFloat} = conj(flip(x))
     dim (1,1) {mustBePositive, mustBeInteger} = findSingletonDim(x, y)
     shape (1,1) string {mustBeMember(shape, ["full", "same", "valid"])} = 'full'
     kwargs.gpu (1,1) logical = isa(x, 'gpuArray') || isa(y, 'gpuArray')
@@ -69,9 +70,19 @@
 % check data sizes
 sz_x = size(x, 1:D);
 sz_y = size(y, 1:D);
-assert(numel(sz_x) == numel(sz_y) && all(sz_x(idim) == sz_y(idim)),...
+assert(all(sz_x(idim) == sz_y(idim) | sz_x(idim) == 1 | sz_y(idim) == 1),...
     "Incompatible sizes [" + join(string((sz_x))+",") + "], and [" + join(string((sz_y))+",") + "].");
 
+% explicit broadcast
+i = sz_x(idim) ~= sz_y(idim) & sz_x(idim) == 1;
+j = sz_y(idim) ~= sz_x(idim) & sz_y(idim) == 1;
+[xrp, yrp] = deal(ones(1,D)); % replication sizes
+xrp(idim(i)) = sz_y(idim(i));
+yrp(idim(j)) = sz_x(idim(j));
+if any(xrp ~= 1), x = repmat(x, xrp); sz_x = size(x, 1:D); end
+if any(yrp ~= 1), y = repmat(y, yrp); sz_y = size(y, 1:D); end
+assert(isequal(size(x, idim), size(y, idim)));
+
 % get computation/output data type and precision
 complex_type = ~isreal(x) || ~isreal(y);
 gpu_type    = isa(x, 'gpuArray') || isa(y, 'gpuArray');

src/ChannelData.m

@@ -389,12 +389,14 @@
                     ysz = size(x);
                     ysz(3) = Trans.numelements;
                     y = zeros(ysz, 'like', x);
+                    N = 256; % max number of channels - modulo for matrix arrays
+                    mod1 = @(x,N) mod(x-1,N)+1; % convert 0-based to 1-based modulo
 
                     % load into output
                     for a = unique(as)' % for each aperture
                         j = a == as; % matching aperture
                         k = aps(:,a); % channel indices
-                        y(:,j,k~=0,:) = x(:,j,k(k~=0),:); % load
+                        y(:,j,k~=0,:) = x(:,j,mod1(k(k~=0),N),:); % load
                     end
                     x = y; % (time x acq x elem x frame)
                 end

src/Sequence.m

@@ -742,7 +742,7 @@
             elseif ~any(tau,'all') % no delays -> FSA
                 seq = Sequence("type","FSA", "numPulse",numel(TX));
 
-            elseif all(all(abs(pog) <= 1e-12, 2) & all(rf == 0,1) & any(ang,'all'),1) % PW (with tolerance)
+            elseif all(all(rf == 0,1) & any(ang,'all'),1) % PW (with tolerance)
                 az = rad2deg(ang(:,1)'); % azimuth
                 el = rad2deg(ang(:,2)'); % elevation
                 if any(el)

src/Transducer.m

@@ -571,25 +571,25 @@ function getSIMUSParam(xdc)
 
     % FDTD functions
     methods
-        function [mask, el_weight, el_dist, el_ind] = elem2grid(xdc, scan, el_sub_div)
+        function [mask, el_weight, el_dist, el_ind] = elem2grid(xdc, grid, el_sub_div)
             % ELEM2GRID - Transducer element to grid mapping
             %
-            % [mask, el_weight, el_dist, el_ind] = ELEM2GRID(xdc, scan)
+            % [mask, el_weight, el_dist, el_ind] = ELEM2GRID(xdc, grid)
             % returns a binary mask and a set of indices, weights, and
             % distances defined for each element give a Transducer xdc and
-            % a ScanCartesian scan. The element indicies are defined on the
+            % a ScanCartesian grid. The element indicies are defined on the
             % vectorized non-zero indices of the mask i.e. on the indices
             % of `find(mask)`.
             %
-            % [...] = ELEM2GRID(xdc, scan, el_sub_div) subdivides the
+            % [...] = ELEM2GRID(xdc, grid, el_sub_div) subdivides the
             % elements in width and height by el_sub_div. The default is
             % [1,1].
             % 
             % See also ULTRASOUNDSYSTEM/KSPACEFIRSTORDER
 
             arguments
                 xdc (1,1) Transducer
-                scan (1,1) ScanCartesian
+                grid (1,1) ScanCartesian
                 el_sub_div (1,2) {mustBeInteger, mustBePositive} = [1,1]
             end
 
@@ -599,19 +599,19 @@ function getSIMUSParam(xdc)
             vec = @(x)x(:); % define locally for compatibility
 
             % cast grid points to single type for efficiency and base grid on the origin
-            [gxv, gyv, gzv] = deal(single(scan.x), single(scan.y), single(scan.z)); % grid {dim} vector
-            [gx0, gy0, gz0] = deal(scan.x(1), scan.y(1), scan.z(1)); % grid {dim} first point
-            [gxd, gyd, gzd] = deal(single(scan.dx), single(scan.dy), single(scan.dz)); % grid {dim} delta
+            [gxv, gyv, gzv] = deal(single(grid.x), single(grid.y), single(grid.z)); % grid {dim} vector
+            [gx0, gy0, gz0] = deal(grid.x(1), grid.y(1), grid.z(1)); % grid {dim} first point
+            [gxd, gyd, gzd] = deal(single(grid.dx), single(grid.dy), single(grid.dz)); % grid {dim} delta
             % pg = scan.positions(); % 3 x Z x X x Y
-            pdims = arrayfun(@(c){find(c == scan.order)}, 'XYZ'); % dimensions mapping
+            pdims = arrayfun(@(c){find(c == grid.order)}, 'XYZ'); % dimensions mapping
             [xdim, ydim, zdim] = deal(pdims{:});
             iord = arrayfun(@(d)find([pdims{:}] == d), [1,2,3]); % inverse mapping
 
             % get array sizing - size '1' and step 'inf' for sliced dimensions
-            [Nx, Ny, Nz] = deal(scan.nx, scan.ny, scan.nz); %#ok<ASGLU>
+            [Nx, Ny, Nz] = deal(grid.nx, grid.ny, grid.nz); %#ok<ASGLU>
 
             % get local element size reference
-            [width_, height_, sz_] = deal(xdc.width, xdc.height, scan.size); %#ok<ASGLU>
+            [width_, height_, sz_] = deal(xdc.width, xdc.height, grid.size); %#ok<ASGLU>
 
             % get regions for each receive transducer element
             el_cen = xdc.positions(); % center of each element
@@ -701,7 +701,7 @@ function getSIMUSParam(xdc)
             assert(~any(cellfun(@isempty, mask_ind), 'all'), 'Unable to assign all (sub-)elements to the grid.');
 
             % reduce to make full recording sensor mask
-            mask = false(scan.size);
+            mask = false(grid.size);
             mask(unique(cat(1, mask_ind{:}))) = true;
 
             % get the translation map for grid mask index to receiver element

src/UltrasoundSystem.m

@@ -312,7 +312,11 @@ function delete(us)
             arguments, self (1,1) UltrasoundSystem, end % always scalar in            
             other = [email protected](self); % shallow copy handle
             other.tmp_folder = tempname(); % new temp dir
-            copyfile(self.tmp_folder, other.tmp_folder); % copy binaries over
+            if exist(self.tmp_folder, 'dir')
+                copyfile(self.tmp_folder, other.tmp_folder); % copy binaries over
+            else
+                mkdir(other.tmp_folder); % create empty
+            end
             addpath(other.tmp_folder); % add to path
             % disp("[DEBUG]: Adding path " + other.tmp_folder + " with " + (numel(dir(other.tmp_folder)) - 2) + " files.");
         end
@@ -920,13 +924,13 @@ function delete(us)
 
             arguments
                 uchannel_data (1,1) uff.channel_data
-                uscan (1,1) uff.scan
+                uscan uff.scan {mustBeScalarOrEmpty} = uff.scan.empty
             end
             fs = uchannel_data.sampling_frequency; % sampling frequency
             seq = Sequence.UFF(uchannel_data.sequence, uchannel_data.sound_speed); % Sequence
             xdc = Transducer.UFF(uchannel_data.probe); % Transducer
             us = UltrasoundSystem('xdc', xdc, 'seq', seq, 'fs', fs);
-            if nargin >= 2, us.scan = Scan.UFF(uscan); end
+            if ~isempty(uscan), us.scan = Scan.UFF(uscan); end
             if nargout >= 2, chd = ChannelData.UFF(uchannel_data, us.seq, us.xdc); end % data
         end
     end
@@ -2206,7 +2210,7 @@ function delete(us)
             % See also ULTRASOUNDSYSTEM/FULLWAVESIM PARALLEL.JOB/FETCHOUTPUTS
             arguments % required arguments
                 us (1,1) UltrasoundSystem
-                med Medium
+                med (1,1) Medium = Medium("c0", us.seq.c0, "rho0", us.seq.c0 / 1.5);
                 sscan (1,1) ScanCartesian = us.scan
             end
             arguments(Repeating)
@@ -2225,9 +2229,9 @@ function delete(us)
                 kwargs.gpu (1,1) logical = logical(gpuDeviceCount()) % whether to employ gpu during pre-processing
             end
             arguments % kWaveArray arguments - these are passed to kWaveArray
-                karray_args.UpsamplingRate (1,1) double =  10, ...
-                karray_args.BLITolerance (1,1) double = 0.05, ...
-                karray_args.BLIType (1,1) string {mustBeMember(karray_args.BLIType, ["sinc", "exact"])} = 'sinc', ... stencil - exact or sinc
+                karray_args.UpsamplingRate (1,1) double =  10
+                karray_args.BLITolerance (1,1) double = 0.05
+                karray_args.BLIType (1,1) string {mustBeMember(karray_args.BLIType, ["sinc", "exact"])} = 'sinc' % stencil - exact or sinc
             end
             arguments % kWave 1.1 arguments - these are passed to kWave
                 kwave_args.BinaryPath (1,1) string = getkWavePath('binaries')
@@ -2358,10 +2362,9 @@ function delete(us)
                         psigv = sample(txsig, psigv); % sample the time delays (J''' x M x T' x V)
                         psigv = wnorm .* el_weight .* swapdim(apodv,3,4) .* psigv; % per sub-element transmit waveform (J''' x M x T' x V x 3)
                         psigv = reshape(psigv, [prod(size(psigv,1:2)), size(psigv,3:4), 1, size(psigv,5)]); % per element transmit waveform (J'' x T' x V x 1 x 3)
-                        psigvproto = psigv(1);
                         psigv = double(psigv); % in-place ?
                         psigv = reshape(el_map_grd * psigv(:,:), [size(el_map_grd,1), size(psigv,2:5)]); % per grid-point transmit waveform (J' x T' x V x 1 x 3)
-                        psigv = cast(psigv, 'like', psigvproto); % in-place ?
+                        psigv = cast(psigv, 'like', psigv([])); % in-place ?
                         psig(:,:,v+b,:,:) = gather(psigv); % store
                     end
                     end
@@ -4864,18 +4867,15 @@ function delete(us)
                 arch (:,1) string {mustBeScalarOrEmpty, mustBeArch(arch)} = nvarch()
             end
 
-            % src file folder
-            src_folder = UltrasoundSystem.getSrcFolder();
-
             % compile each
             for i = numel(defs):-1:1
                 d = defs(i);
                 % make full command
                 com(i) = join(cat(1,...
                     "nvcc ", ...
-                    "--ptx " + fullfile(src_folder, d.Source), ...
+                    "--ptx " + which(string(d.Source)), ...
                     "-arch=" + arch + " ", ... compile for active gpu
-                    "-o " + fullfile(us.tmp_folder, strrep(d.Source, '.cu', '.ptx')), ...
+                    "-o " + fullfile(us.tmp_folder, argn(2, @fileparts, d.Source) + ".ptx"), ...
                     join("--" + d.CompileOptions),...
                     join("-I" + d.IncludePath), ...
                     join("-L" + d.Libraries), ...

test/KernTest.m

@@ -95,9 +95,13 @@ function runconvd(test, prec, complexity, dev)
             dims = [1,2,3]; 
 
             % check that defaults and options work
+            test.assertTrue(isalmostn(convd(A), xcorr(A)));
             convd(A, B);
             convd(B, A, 'lowmem', true);
 
+            % broadcasting
+            test.assertTrue( all(convd(cat(4,A,A), B) == convd(A, B), 'all') );
+
             % convolve and check outputs
             for s = shapes, for d = dims %#ok<ALIGN> 
                 z1 = cell2mat(cellfun(@(A,B) {gather(conv(A, B, char(s)))}, num2cell(A,1), num2cell(B,1)));