From 8900fb11b1c151dc82eed59eaf64335331821299 Mon Sep 17 00:00:00 2001 From: "ray@omni" Date: Wed, 12 Jun 2019 18:02:43 +0800 Subject: [PATCH] (M for modified) M mri/mri_sensemap_denoise.m Forcing parfor to avoid queueing sMaps estimation M penalty/Cdiff.m ndims()'s return always >= 2 M systems/@fatrix2/cat.m Avoid nesting @fatrix2 when numel(varargin)==1 M systems/@fatrix2/mtimes.m Enabling A*c, where A is a @fatrix2, c is a scalar M systems/@fatrix2/private/fatrix2_subsref_colon.m Use parfor for fatrix2-matrix product, when parpool is already init'ed, (init parpool can be expensive by itself, hence parfor only when parpool is already available) M systems/Gmri.m Enabling default nufft_args to support 3D; Substituting repmat w/ bsxfun, RAM friendly. --- mri/mri_sensemap_denoise.m | 20 +++++++---- penalty/Cdiff.m | 4 +-- systems/@fatrix2/cat.m | 4 +-- systems/@fatrix2/mtimes.m | 6 +++- .../@fatrix2/private/fatrix2_subsref_colon.m | 16 ++++++--- systems/Gmri.m | 36 +++++++++---------- 6 files changed, 49 insertions(+), 37 deletions(-) diff --git a/mri/mri_sensemap_denoise.m b/mri/mri_sensemap_denoise.m index 0725234..08a4797 100644 --- a/mri/mri_sensemap_denoise.m +++ b/mri/mri_sensemap_denoise.m @@ -131,17 +131,20 @@ % initial maps if isempty(arg.init) - arg.init = zeros([prod(NN) ncoil]); + init = zeros([prod(NN) ncoil]); good = abs(arg.bodycoil) > arg.thresh * max(abs(arg.bodycoil(:))); % 2010-06-24 - for ic = 1:ncoil + + bodycoil = arg.bodycoil; % extract from struct and cache to accelerate loop + parfor ic = 1:ncoil zj = stackpick(ykj,ic); - tmp = zj ./ arg.bodycoil; % usual ratio + tmp = zj ./ bodycoil; % usual ratio if 1 % set all uncertain map values to median of good ones % good = abs(zj) > arg.thresh * max(abs(zj(:))); % prior to 2010-06-24 tmp(~good) = median(abs(tmp(good))); end - arg.init(:,ic) = tmp(:); + init(:,ic) = tmp(:); end + arg.init = init; sinit = reshape(arg.init, [NN ncoil]); % return to caller if wanted end @@ -197,6 +200,7 @@ L = ichol(H, struct('michol', 'on', 'diagcomp', alpha)); end else + H = []; % any variable presents in a parfor-loop is required to be init'ed A = Gdiag(arg.bodycoil(mask(:)), 'mask', mask); end @@ -233,7 +237,9 @@ fail('bad arg.precon "%s"', arg.precon) end -for ic = 1:ncoil +C = R.C; +nargo = nargout; % cache nargout for parfor +parfor ic = 1:ncoil ytmp = stackpick(ykj,ic); if arg.chol % numerical inverse via backslash for small problems ytmp = double(ytmp); @@ -243,10 +249,10 @@ else % run qpwls algorithm for regularized fitting init = stackpick(arg.init,ic); tmp = qpwls_pcg1(init(mask), A, 1, ... - ytmp(mask), R.C, ... + ytmp(mask), C, ... 'precon', arg.precon, ... 'niter', arg.niter, 'isave', arg.isave); - if nargout > 2 % cost + if nargo > 2 % cost cost(:,ic) = pwls_cost(tmp, A, 1, ytmp(mask(:)), R); end end diff --git a/penalty/Cdiff.m b/penalty/Cdiff.m index eb12da3..cec21e3 100644 --- a/penalty/Cdiff.m +++ b/penalty/Cdiff.m @@ -84,9 +84,7 @@ end if isempty(arg.offsets) - if ndims(kappa) == 1 - arg.offsets = [1]; - elseif ndims(kappa) == 2 + if ndims(kappa) == 2 nx = arg.dim_i(1); arg.offsets = [1 nx nx+1 nx-1]; % default 2D elseif ndims(kappa) == 3 diff --git a/systems/@fatrix2/cat.m b/systems/@fatrix2/cat.m index de41c7e..4843808 100644 --- a/systems/@fatrix2/cat.m +++ b/systems/@fatrix2/cat.m @@ -11,11 +11,11 @@ % purge any empty cells so [[]; ob] returns ob is_empty = cellfun(@isempty, varargin); if any(is_empty) -% varargin = {varargin{~is_empty}}; - varargin = {varargin(~is_empty)}; % 2016-04-19 + varargin = varargin(~is_empty); % 2016-04-19 end if numel(varargin) == 1 ob = varargin{1}; + return; end switch dim diff --git a/systems/@fatrix2/mtimes.m b/systems/@fatrix2/mtimes.m index 16bc6f0..d6b4484 100644 --- a/systems/@fatrix2/mtimes.m +++ b/systems/@fatrix2/mtimes.m @@ -28,7 +28,11 @@ end end -y = fatrix2_mtimes_vector(ob, x); % "ordinary" fatrix * vector multiplication +if isscalar(x) + y = fatrix2_scalar_times(x, ob); +else + y = fatrix2_mtimes_vector(ob, x); % "ordinary" fatrix * vector multiplication +end end % mtimes() diff --git a/systems/@fatrix2/private/fatrix2_subsref_colon.m b/systems/@fatrix2/private/fatrix2_subsref_colon.m index 77fe18c..43621b2 100644 --- a/systems/@fatrix2/private/fatrix2_subsref_colon.m +++ b/systems/@fatrix2/private/fatrix2_subsref_colon.m @@ -82,10 +82,18 @@ out = zeros(size(ob,1), length(jj), class(tmp)); out(:,1) = tmp; -for nn=2:length(jj) - x = z; - x(jj(nn)) = 1; - out(:,nn) = ob * x; +if ~isempty(gcp('nocreate')) + parfor nn=2:length(jj) + x = z; + x(jj(nn)) = 1; + out(:,nn) = ob * x; + end +else + for nn=2:length(jj) + x = z; + x(jj(nn)) = 1; + out(:,nn) = ob * x; + end end diff --git a/systems/Gmri.m b/systems/Gmri.m index e51ca51..8358757 100644 --- a/systems/Gmri.m +++ b/systems/Gmri.m @@ -79,7 +79,9 @@ arg.L = []; arg.aL = []; % for autocorrelation zmap. default is arg.L N = size(mask); -arg.nufft_args = {N, [6 6], 2*N, N/2, 'table', 2^10, 'minmax:kb'}; +if iscolumn(mask), N = N(1); end +% arg.nufft_args = {N, [6 6], 2*N, N/2, 'table', 2^10, 'minmax:kb'}; +arg.nufft_args = {N, 6 * ones(1, numel(N)), 2*N, N/2, 'table', 2^10, 'minmax:kb'}; arg.n_shift = []; arg.exp_approx_args = {}; % arguments for mri_exp_approx() @@ -104,18 +106,14 @@ dd = size(kspace,2); arg.kspace = kspace; -if length(arg.fov) == 1 - arg.fov = arg.fov * ones(1,dd); -end +if isscalar(arg.fov) == 1, arg.fov = arg.fov * ones(1,dd); end arg.mask = mask; arg.Nd = size(mask); arg.dim = [size(kspace,1) sum(mask(:))]; -omega = zeros(size(kspace), class(kspace)); -for id=1:dd - omega(:,id) = 2*pi*kspace(:,id) * arg.fov(id) / arg.Nd(id); -end -if max(abs(omega(:))) > pi+1e-6 +omega = 2*pi * bsxfun(@times, kspace, arg.fov./arg.Nd); + +if any(abs(omega(:)) > pi+1e-6) warn 'omega exceeds pi. are you sure you intended this?' end @@ -384,7 +382,7 @@ end y = exp_xform_mex(complexify(single(x)), ... complexify(arg.u), complexify(arg.v)); -y = y .* repmat(arg.basis.transform, [1 ncol(y)]); % [nd *nc] +y = bsxfun(@times, y, arg.basis.transform); % Gmri_forw_exact(): y = A * x @@ -392,7 +390,7 @@ x = reshapee(x, prod(arg.Nd), []); % [(N) (nc)] to [*N *nc] x = x(arg.mask,:); % [np *nc] y = exp_xform_mex(complexify(single(x)), arg.u, arg.v); % [nd *nc] -y = y .* repmat(arg.basis.transform, [1 ncol(y)]); % [nd *nc] +y = bsxfun(@times, y, arg.basis.transform); % Gmri_forw(): y = A * x @@ -402,26 +400,25 @@ x = reshape(x, prod(arg.Nd), []); % [(N) (nc)] to [*N *nc] x = x(arg.mask,:); % [np *nc] end -nc = ncol(x); if isempty(arg.zmap) y = arg.Gnufft * x; else % approximation y = 0; for ll=1:arg.L - tmp = repmat(arg.Ct(:,ll), [1 nc]) .* x; + tmp = bsxfun(@times, arg.Ct(:,ll), x); tmp = arg.Gnufft * tmp; - tmp = repmat(arg.B(:,ll), [1 nc]) .* tmp; + tmp = bsxfun(@times, arg.B(:,ll), tmp); y = y + tmp; % y = y + arg.B(:,ll) .* (arg.Gnufft * (arg.Ct(:,ll) .* x)); end end -y = y .* repmat(arg.basis.transform, [1 nc]); % [nd *nc] +y = squeeze(bsxfun(@times, y, arg.basis.transform)); % Gmri_back_exact_Fatrix(): x = A' * y function x = Gmri_back_exact_Fatrix(arg, y) -y = y .* repmat(conj(arg.basis.transform), [1 ncol(y)]); % [nd nc] +y = bsxfun(@times, y, conj(arg.basis.transform)); % trick: conj(exp(-uv)) = exp(-conj(u) conj(v)) vc = complexify(conj(arg.v)); uc = complexify(conj(arg.u)); @@ -437,17 +434,16 @@ % Gmri_back_Fatrix(): x = A' * y % full adjoint ("back-projection") function x = Gmri_back_Fatrix(arg, y) -nc = ncol(y); -y = y .* repmat(conj(arg.basis.transform), [1 nc]); +y = bsxfun(@times, y, conj(arg.basis.transform)); if isempty(arg.zmap) x = arg.Gnufft' * y; else % approximation x = 0; for ll=1:arg.L - tmp = repmat(conj(arg.B(:,ll)), [1 nc]) .* y; + tmp = bsxfun(@times, conj(arg.B(:,ll)), y); tmp = arg.Gnufft' * tmp; - tmp = repmat(conj(arg.Ct(:,ll)), [1 nc]) .* tmp; + tmp = bsxfun(@times, conj(arg.Ct(:,ll)), tmp); x = x + tmp; end end