Some tiny performance improvements

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

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

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

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()
 

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
 
 

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,15 +382,15 @@
 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
 function y = Gmri_forw_exact(arg, x)
 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