text:
function VO = spm_segment(VF,PG,flags)
% Segment an MR image into Gray, White & CSF.
%
% FORMAT VO = spm_segment(PF,PG,flags)
% PF    - name(s) of image(s) to segment (must have same dimensions).
% PG    - name(s) of template image(s) for realignment.
%       - or a 4x4 transformation matrix which maps from the image to
%         the set of templates.
% flags - a structure normally based on defaults.segment
% VO    - optional output volume
%
%                      The algorithm is four step:
%
% 1) Determine the affine transform which best matches the image with a
%    template image. If the name of more than one image is passed, then
%    the first image is used in this step. This step is not performed if
%    no template images are specified.
%
% 2) Perform Cluster Analysis with a modified Mixture Model and a-priori
%    information about the likelihoods of each voxel being one of a
%    number of different tissue types. If more than one image is passed,
%    then they they are all assumed to be in register, and the voxel
%    values are fitted to multi-normal distributions.
%
% 3) Perform morphometric operations on the grey and white partitions
%    in order to more accurately identify brain tissue. This is then used
%    to clean up the grey and white matter segments.
%
% 4) If no output argument is specified, then the segmented images are
%    written to disk. The names of these images have "_seg1", "_seg2"
%    & "_seg3" appended to the name of the first image passed.
%
%_______________________________________________________________________
% Refs:
%
% Ashburner J & Friston KJ (1997) Multimodal Image Coregistration and
% Partitioning - a Unified Framework. NeuroImage 6:209-217
%
%_______________________________________________________________________
%
% The template image, and a-priori likelihood images are modified
% versions of those kindly supplied by Alan Evans, MNI, Canada
% (ICBM, NIH P-20 project, Principal Investigator John Mazziotta).
%_______________________________________________________________________
% @(#)spm_segment.m	2.28a John Ashburner 04/04/01

% Create some suitable default values
%-----------------------------------------------------------------------

def_flags.estimate.priors = str2mat(...
        fullfile(spm('Dir'),'apriori','gray.mnc'),...
        fullfile(spm('Dir'),'apriori','white.mnc'),...
        fullfile(spm('Dir'),'apriori','csf.mnc'));
def_flags.estimate.reg    = 0.01;
def_flags.estimate.cutoff = 30;
def_flags.estimate.samp   = 3;
def_flags.estimate.bb     =  [[-88 88]' [-122 86]' [-60 95]'];
def_flags.estimate.affreg.smosrc = 8;
def_flags.estimate.affreg.regtype = 'mni';
def_flags.estimate.affreg.weight = '';
def_flags.write.cleanup   = 1;
def_flags.write.wrt_cor   = 1;
def_flags.graphics        = 1;

if nargin<3, flags = def_flags; end;
if ~isfield(flags,'estimate'),        flags.estimate        = def_flags.estimate;        end;
if ~isfield(flags.estimate,'priors'), flags.estimate.priors = def_flags.estimate.priors; end;
if ~isfield(flags.estimate,'reg'),    flags.estimate.reg    = def_flags.estimate.reg;    end;
if ~isfield(flags.estimate,'cutoff'), flags.estimate.cutoff = def_flags.estimate.cutoff; end;
if ~isfield(flags.estimate,'samp'),   flags.estimate.samp   = def_flags.estimate.samp;   end;
if ~isfield(flags.estimate,'bb'),     flags.estimate.bb     = def_flags.estimate.bb;     end;
if ~isfield(flags.estimate,'affreg'), flags.estimate.affreg = def_flags.estimate.affreg; end;
if ~isfield(flags.estimate.affreg,'smosrc'),
	flags.estimate.affreg.smosrc = def_flags.estimate.affreg.smosrc;
end;
if ~isfield(flags.estimate.affreg,'regtype'),
	flags.estimate.affreg.regtype = def_flags.estimate.affreg.regtype;
end;
if ~isfield(flags.estimate.affreg,'weight'),
	flags.estimate.affreg.weight = def_flags.estimate.affreg.weight;
end;
if ~isfield(flags,'write'),         flags.write         = def_flags.write;         end;
if ~isfield(flags.write,'cleanup'), flags.write.cleanup = def_flags.write.cleanup; end;
if ~isfield(flags.write,'wrt_cor'), flags.write.wrt_cor = def_flags.write.wrt_cor; end;
if ~isfield(flags,'graphics'),      flags.graphics      = def_flags.graphics;      end;

%-----------------------------------------------------------------------

if ischar(VF), VF= spm_vol(VF); end;

SP         = init_sp(flags.estimate,VF,PG);
[x1,x2,x3] = get_sampling(SP.MM,VF,flags.estimate.samp,flags.estimate.bb);
BP         = init_bp(VF, flags.estimate.cutoff, flags.estimate.reg);
CP         = init_cp(VF,x3);
sums       = zeros(8,1);

for pp=1:length(x3),
	[raw,msk] = get_raw(VF,x1,x2,x3(pp));
	s         = get_sp(SP,x1,x2,x3(pp));
	CP        = update_cp_est(CP,s,raw,msk,pp);
	sums      = sums + reshape(sum(sum(s,1),2),8,1);
end;
sums = sums/sum(sums);
CP   = shake_cp(CP);

[CP,BP,SP] = run_segment(CP,BP,SP,VF,sums,x1,x2,x3);

%save segmentation_results.mat CP BP SP VF sums

[g,w,c] = get_gwc(VF,BP,SP,CP,sums,flags.write.wrt_cor);
if flags.write.cleanup, [g,w,c] = clean_gwc(g,w,c); end;

% Create the segmented images.
%-----------------------------------------------------------------------
%offs  = cumsum(repmat(prod(VF(1).dim(1:2)),1,VF(1).dim(3)))-prod(VF(1).dim(1:2));
%pinfo = [repmat([1/255 0]',1,VF(1).dim(3)) ; offs];
[pth,nm,xt] = fileparts(deblank(VF(1).fname));
for j=1:3,
	tmp   = fullfile(pth,[nm '_seg' num2str(j) xt]);
	VO(j) = struct(...
		'fname',tmp,...
		'dim',    [VF(1).dim(1:3) 2],...
		'mat',    VF(1).mat,...
		'pinfo',  [1/255 0 0]',...
		'descrip','Segmented image');
end;

if nargout==0,
	VO = spm_create_vol(VO);

	spm_progress_bar('Init',VF(1).dim(3),'Writing Segmented','planes completed');
	for pp=1:VF(1).dim(3),
		VO(1) = spm_write_plane(VO(1),double(g(:,:,pp))/255,pp);
		VO(2) = spm_write_plane(VO(2),double(w(:,:,pp))/255,pp);
		VO(3) = spm_write_plane(VO(3),double(c(:,:,pp))/255,pp);
		spm_progress_bar('Set',pp);
	end;
	VO = spm_close_vol(VO);
	spm_progress_bar('Clear');
end;

VO(1).dat = g; VO(1).pinfo = VO(1).pinfo(1:2,:);
VO(2).dat = w; VO(2).pinfo = VO(2).pinfo(1:2,:);
VO(3).dat = c; VO(3).pinfo = VO(3).pinfo(1:2,:);

if flags.graphics, display_graphics(VF,VO,CP.mn,CP.cv,CP.mg); end;

return;
%=======================================================================

%=======================================================================
function [y1,y2,y3] = affine_transform(x1,x2,x3,M)
	y1 = M(1,1)*x1 + M(1,2)*x2 + M(1,3)*x3 + M(1,4);
	y2 = M(2,1)*x1 + M(2,2)*x2 + M(2,3)*x3 + M(2,4);
	y3 = M(3,1)*x1 + M(3,2)*x2 + M(3,3)*x3 + M(3,4);
return;
%=======================================================================

%=======================================================================
function display_graphics(VF,VS,mn,cv,mg)
% Do the graphics
nb = 3;
spm_figure('Clear','Graphics');
fg = spm_figure('FindWin','Graphics');
if ~isempty(fg),
	% Show some text
	%-----------------------------------------------------------------------
	ax = axes('Position',[0.05 0.8 0.9 0.2],'Visible','off','Parent',fg);
	text(0.5,0.80, 'Segmentation','FontSize',16,'FontWeight','Bold',...
		'HorizontalAlignment','center','Parent',ax);

	text(0,0.65, ['Image:  ' spm_str_manip(VF(1).fname,'k50d')],...
		'FontSize',14,'FontWeight','Bold','Parent',ax);

	text(0,0.40, 'Means:','FontSize',12,'FontWeight','Bold','Parent',ax);
	text(0,0.30, 'Std devs:' ,'FontSize',12,'FontWeight','Bold','Parent',ax);
	text(0,0.20, 'N vox:','FontSize',12,'FontWeight','Bold','Parent',ax);
	for j=1:nb,
		text((j+0.5)/(nb+1),0.40, num2str(mn(1,j)),...
			'FontSize',12,'FontWeight','Bold',...
			'HorizontalAlignment','center','Parent',ax);
		text((j+0.5)/(nb+1),0.30, num2str(sqrt(cv(1,1,j))),...
			'FontSize',12,'FontWeight','Bold',...
			'HorizontalAlignment','center','Parent',ax);
		text((j+0.5)/(nb+1),0.20, num2str(mg(1,j)/sum(mg(1,:))),...
			'FontSize',12,'FontWeight','Bold',...
			'HorizontalAlignment','center','Parent',ax);
	end;
	if length(VF) > 1,
		text(0,0.10,...
		'Note: only means and variances for the first image are shown',...
		'Parent',ax,'FontSize',12);
	end;

	M1 = VS(1).mat;
	M2 = VF(1).mat;
	for i=1:5,
		M   = spm_matrix([0 0 i*VF(1).dim(3)/6]);
		img = spm_slice_vol(VF(1),M,VF(1).dim(1:2),1);
		img(1,1) = eps;
		ax = axes('Position',...
			[0.05 0.75*(1-i/5)+0.05 0.9/(nb+1) 0.75/5],...
			'Visible','off','Parent',fg);
		imagesc(rot90(img), 'Parent', ax);
		set(ax,'Visible','off','DataAspectRatio',[1 1 1]);

		for j=1:3,
			img = spm_slice_vol(VS(j),M2\M1*M,VF(1).dim(1:2),1);
			ax  = axes('Position',...
				[0.05+j*0.9/(nb+1) 0.75*(1-i/5)+0.05 0.9/(nb+1) 0.75/5],...
				'Visible','off','Parent',fg);
			image(rot90(img*64), 'Parent', ax);
			set(ax,'Visible','off','DataAspectRatio',[1 1 1]);
		end;
	end;

	spm_print;
	drawnow;
end;
return;
%=======================================================================

%=======================================================================
function M = get_affine_mapping(VF,VG,aflags)

if ~isempty(VG) & ischar(VG), VG = spm_vol(VG); end;

if ~isempty(VG) & isstruct(VG),
	% Affine registration so that a priori images match the image to
	% be segmented.
	%-----------------------------------------------------------------------

	VFS = spm_smoothto8bit(VF(1),aflags.smosrc);

	% Scale all images approximately equally
	% ---------------------------------------------------------------
	for i=1:length(VG),
		VG(i).pinfo(1:2,:) = VG(i).pinfo(1:2,:)/spm_global(VG(i));
	end;
	VFS(1).pinfo(1:2,:) = VFS(1).pinfo(1:2,:)/spm_global(VFS(1));

	spm_chi2_plot('Init','Affine Registration','Mean squared difference','Iteration');
	flags     = struct('sep',aflags.smosrc, 'regtype',aflags.regtype,'WG',[],'globnorm',0,'debug',0);
	M         = eye(4);
	[M,scal]  = spm_affreg(VG, VFS, flags, M);

	if ~isempty(aflags.weight), flags.WG = spm_vol(aflags.weight); end;

	flags.sep = aflags.smosrc/2;
	M         = spm_affreg(VG, VFS, flags, M,scal);
	spm_chi2_plot('Clear');

elseif all(size(VG) == [4 4])
	% Assume that second argument is a matrix that will do the job
	%-----------------------------------------------------------------------
	M = VG;
else
	% Assume that image is normalized
	%-----------------------------------------------------------------------
	M = eye(4);
end
return;
%=======================================================================

%=======================================================================
function [x1,x2,x3] = get_sampling(MM,VF,samp,bb1)
% Voxels to sample during the cluster analysis
%-----------------------------------------------------------------------

% A bounding box for the brain in Talairach space.
%bb = [ [-88 88]' [-122 86]' [-60 95]'];
%c = [bb(1,1) bb(1,2) bb(1,3) 1
%     bb(1,1) bb(1,2) bb(2,3) 1
%     bb(1,1) bb(2,2) bb(1,3) 1
%     bb(1,1) bb(2,2) bb(2,3) 1
%     bb(2,1) bb(1,2) bb(1,3) 1
%     bb(2,1) bb(1,2) bb(2,3) 1
%     bb(2,1) bb(2,2) bb(1,3) 1
%     bb(2,1) bb(2,2) bb(2,3) 1]';
%tc = MM\c;
%tc = tc(1:3,:)';
%mx = max(tc);
%mn = min(tc);
%bb = [mn ; mx];
%vx = sqrt(sum(VF(1).mat(1:3,1:3).^2));
%samp = round(max(abs([4 4 4]./vx), [1 1 1]));
%x1 = bb(1,1):samp(1):bb(2,1);
%x2 = bb(1,2):samp(2):bb(2,2);
%x3 = bb(1,3):samp(3):bb(2,3);
%return;

% A bounding box for the brain in Talairach space.
if nargin<4, bb1 = [ [-88 88]' [-122 86]' [-60 95]']; end;

% A mapping from a unit radius sphere to a hyper-ellipse
% that is just enclosed by the bounding box in Talairach
% space.
M0 = [diag(diff(bb1)/2) mean(bb1)';[0 0 0 1]];

% The mapping from voxels to Talairach space is MM,
% so the ellipse in the space of the image becomes:
M0 = MM\M0;

% So to work out the bounding box in the space of the
% image that just encloses the hyper-ellipse.
tmp = M0(1:3,1:3);
tmp = diag(tmp*tmp'/diag(sqrt(diag(tmp*tmp'))));
bb  = round([M0(1:3,4)-tmp M0(1:3,4)+tmp])';
bb  = min(max(bb,[1 1 1 ; 1 1 1]),[VF(1).dim(1:3) ; VF(1).dim(1:3)]);

% Want to sample about every 3mm
tmp  = sqrt(sum(VF(1).mat(1:3,1:3).^2))';
samp = round(max(abs(tmp.^(-1)*samp), [1 1 1]'));

x1 = bb(1,1):samp(1):bb(2,1);
x2 = bb(1,2):samp(2):bb(2,2);
x3 = bb(1,3):samp(3):bb(2,3);

return;
%=======================================================================

%=======================================================================
function [CP,BP,SP] = run_segment(CP,BP,SP,VF,sums,x1,x2,x3)
oll = -Inf;
spm_chi2_plot('Init','Segmenting','Log-likelihood','Iteration #');

for iter = 1:64,
	ll= 0;
	for pp = 1:length(x3), % Loop over planes
		bf        = get_bp(BP,x1,x2,x3(pp));
		[raw,msk] = get_raw(VF,x1,x2,x3(pp));
		s         = get_sp(SP,x1,x2,x3(pp));
		cor       = bf.*raw;
		[P,ll0]   = get_p(cor,msk,s,sums,CP,bf);
		ll        = ll + ll0;
		CP        = update_cp_est(CP,P,cor,msk,pp);
		BP        = update_bp_est(BP,P,cor,CP,msk,x1,x2,x3(pp));
	end;

	BP = update_bp(BP);
	if iter>1, spm_chi2_plot('Set',ll); end;
	%fprintf('\t%g\n', ll);

	% Stopping criterion
	%-----------------------------------------------------------------------
	if iter == 2,
		ll2 = ll;
	elseif iter > 2 & abs((ll-oll)/(ll-ll2)) < 0.0001
		break;
	end;
	oll = ll;
end;
spm_chi2_plot('Clear');
return;
%=======================================================================

%=======================================================================
function BP = init_bp(VF,co,reg)
m        = length(VF);
tmp      = sqrt(sum(VF(1).mat(1:3,1:3).^2));
BP.nbas  = max(round((VF(1).dim(1:3).*tmp)/co),[1 1 1]);
BP.B1    = spm_dctmtx(VF(1).dim(1),BP.nbas(1));
BP.B2    = spm_dctmtx(VF(1).dim(2),BP.nbas(2));
BP.B3    = spm_dctmtx(VF(1).dim(3),BP.nbas(3));

nbas     = BP.nbas;
if prod(BP.nbas)>1,
	% Set up a priori covariance matrix
	vx = sqrt(sum(VF(1).mat(1:3,1:3).^2));
	kx=(pi*((1:nbas(1))'-1)*pi/vx(1)/VF(1).dim(1)*10).^2;
	ky=(pi*((1:nbas(2))'-1)*pi/vx(2)/VF(1).dim(2)*10).^2;
	kz=(pi*((1:nbas(3))'-1)*pi/vx(3)/VF(1).dim(3)*10).^2;

	% Cost function based on sum of squares of 4th derivatives
	IC0 =  (1*kron(kz.^4,kron(ky.^0,kx.^0)) +...
	        1*kron(kz.^0,kron(ky.^4,kx.^0)) +...
	        1*kron(kz.^0,kron(ky.^0,kx.^4)) +...
	        4*kron(kz.^3,kron(ky.^1,kx.^0)) +...
	        4*kron(kz.^3,kron(ky.^0,kx.^1)) +...
	        4*kron(kz.^1,kron(ky.^3,kx.^0)) +...
	        4*kron(kz.^0,kron(ky.^3,kx.^1)) +...
	        4*kron(kz.^1,kron(ky.^0,kx.^3)) +...
	        4*kron(kz.^0,kron(ky.^1,kx.^3)) +...
	        6*kron(kz.^2,kron(ky.^2,kx.^0)) +...
	        6*kron(kz.^2,kron(ky.^0,kx.^2)) +...
	        6*kron(kz.^0,kron(ky.^2,kx.^2)) +...
	       12*kron(kz.^2,kron(ky.^1,kx.^1)) +...
	       12*kron(kz.^1,kron(ky.^2,kx.^1)) +...
	       12*kron(kz.^1,kron(ky.^1,kx.^2)) )*reg;

	%IC0(1) = max(IC0);
	BP.IC0 = diag(IC0(2:end));

	% Initial estimate for intensity modulation field
	BP.T   = zeros(nbas(1),nbas(2),nbas(3),length(VF));
	%-----------------------------------------------------------------------
else
	BP.T   = zeros([1 1 1 length(VF)]);
	BP.IC0 = [];
end;
BP.Alpha = zeros(prod(BP.nbas(1:3)),prod(BP.nbas(1:3)),m);
BP.Beta  = zeros(prod(BP.nbas(1:3)),m);
return;
%=======================================================================

%=======================================================================
function BP = update_bp_est(BP,p,cor,CP,msk,x1,x2,x3)
if prod(BP.nbas)<=1, return; end;
B1 = BP.B1(x1,:);
B2 = BP.B2(x2,:);
B3 = BP.B3(x3,:);
for j=1:size(BP.Alpha,3),
	cr  = cor(:,:,j);
	w1 = zeros(size(cr));
	w2 = zeros(size(cr));
	for i=[1 2 3 4 5 6 7 8],
		tmp = p(:,:,i)*CP.cv(j,j,i)^(-1);
		w1  = w1 + tmp.*(CP.mn(j,i) - cr);
		w2  = w2 + tmp;
	end;
	wt1       = 1 + cr.*w1;
	wt2       = cr.*(cr.*w2 - w1);
	wt1(~msk) = 0;
	wt2(~msk) = 0;

	BP.Beta(:,j)    = BP.Beta(:,j)    + kron(B3',spm_krutil(wt1,B1,B2,0));
	BP.Alpha(:,:,j) = BP.Alpha(:,:,j) + kron(B3'*B3,spm_krutil(wt2,B1,B2,1));
end;
return;
%=======================================================================

%=======================================================================
function BP = update_bp(BP)
if prod(BP.nbas)<=1, return; end;
for j=1:size(BP.Alpha,3),
	x     = BP.T(:,:,:,j);
	x     = x(:);
	x     = x(2:end);
	Alpha = BP.Alpha(2:end,2:end,j);
	Beta  = BP.Beta(2:end,j);
	x     = (Alpha + BP.IC0)\(Alpha*x + Beta);

	BP.T(:,:,:,j) = reshape([0 ; x],BP.nbas(1:3));
	BP.Alpha      = zeros(size(BP.Alpha));
	BP.Beta       = zeros(size(BP.Beta));
end;
return;
%=======================================================================

%=======================================================================
function bf = get_bp(BP,x1,x2,x3)
bf = ones(length(x1),length(x2),size(BP.Alpha,3));
if prod(BP.nbas)<=1, return; end;
B1 = BP.B1(x1,:);
B2 = BP.B2(x2,:);
B3 = BP.B3(x3,:);
for i=1:size(BP.Alpha,3),
	t = reshape(reshape(BP.T(:,:,:,i),...
		BP.nbas(1)*BP.nbas(2),BP.nbas(3))*B3', BP.nbas(1), BP.nbas(2));
	bf(:,:,i) = exp(B1*t*B2');
end;
return;
%=======================================================================

%=======================================================================
function [dat,msk] = get_raw(VF,x1,x2,x3)
[X1,X2,X3] = ndgrid(x1,x2,x3);
for i=1:length(VF),
	[Y1,Y2,Y3] = affine_transform(X1,X2,X3,VF(i).mat\VF(1).mat);
	dat(:,:,i) = spm_sample_vol(VF(i),Y1,Y2,Y3,1);
end;
msk = all(dat,3) & all(isfinite(double(dat)),3);
return;
%=======================================================================

%=======================================================================
function CP = init_cp(VF,x3)
n = 8;
m = length(VF);
p = length(x3);
CP.mom0 = zeros(1,n,p)+eps;
CP.mom1 = zeros(m,n,p);
CP.mom2 = zeros(m,m,n,p)+eps;

% Occasionally the dynamic range of the images is such that many voxels
% all have the same intensity.  Adding cv0 is an attempt to improve the
% stability of the algorithm if this occurs. The value 0.083 was obtained
% from var(rand(1000000,1)).  It prbably isn't the best way of doing
% things, but it appears to work.
CP.cv0 = zeros(m,m);
for i=1:m,
	if spm_type(VF(i).dim(4),'intt'),
		CP.cv0(i,i)=0.083*mean(VF(i).pinfo(1,:));
	end;
end;
return;
%=======================================================================

%=======================================================================
function CP = shake_cp(CP)
CP.mom0(:,5,:)   = CP.mom0(:,1,:);
CP.mom0(:,6,:)   = CP.mom0(:,2,:);
CP.mom0(:,7,:)   = CP.mom0(:,3,:);
CP.mom1(:,5,:)   = CP.mom1(:,1,:);
CP.mom1(:,6,:)   = CP.mom1(:,2,:);
CP.mom1(:,7,:)   = CP.mom1(:,3,:);
CP.mom1(:,8,:)   = 0;
CP.mom2(:,:,5,:) = CP.mom2(:,:,1,:);
CP.mom2(:,:,6,:) = CP.mom2(:,:,2,:);
CP.mom2(:,:,7,:) = CP.mom2(:,:,3,:);
return;
%=======================================================================

%=======================================================================
function CP = update_cp_est(CP,P,dat,msk,p)
m   = size(dat,3);
d   = size(P);
P   = reshape(P,[d(1)*d(2),d(3)]);
dat = reshape(dat,[d(1)*d(2),m]);
P(~msk(:),:)   = [];
dat(~msk(:),:) = [];
for i=1:size(CP.mom0,2),
	CP.mom0(1,i,p)   = sum(P(:,i));
	CP.mom1(:,i,p)   = sum((P(:,i)*ones(1,m)).*dat)';
	CP.mom2(:,:,i,p) = ((P(:,i)*ones(1,m)).*dat)'*dat;
end;

for i=1:size(CP.mom0,2),
	CP.mg(1,i)   = sum(CP.mom0(1,i,:),3);
	CP.mn(:,i)   = sum(CP.mom1(:,i,:),3)/CP.mg(1,i);

	tmp          = (CP.mg(1,i).*CP.mn(:,i))*CP.mn(:,i)';
	tmp          = tmp-eye(size(tmp))*eps*10000;
	CP.cv(:,:,i) = (sum(CP.mom2(:,:,i,:),4) - tmp)/CP.mg(1,i) + CP.cv0;
end;
CP.mg   = CP.mg/sum(CP.mg);
return;
%=======================================================================

%=======================================================================
function [p,ll] = get_p(cor,msk,s,sums,CP,bf)
d   = [size(cor) 1 1];
n   = size(CP.mg,2);
cor = reshape(cor,d(1)*d(2),d(3));
cor = cor(msk,:);
p   = zeros(d(1)*d(2),n);
if ~any(msk), p  = reshape(p,d(1),d(2),n); ll=0; return; end;

for i=1:n,
	amp       = 1/sqrt((2*pi)^d(3) * det(CP.cv(:,:,i)));
	dst       = (cor-ones(size(cor,1),1)*CP.mn(:,i)')/sqrtm(CP.cv(:,:,i));
	dst       = sum(dst.*dst,2);
	tmp       = s(:,:,i);
	p(msk,i)  = (amp*CP.mg(1,i)/sums(i))*exp(-0.5*dst).*tmp(msk) +eps;
end;
sp = sum(p,2);
ll = sum(log(sp(msk).*bf(msk)+eps));
sp(~msk) = Inf;
for i=1:n, p(:,i) = p(:,i)./sp; end;
p  = reshape(p,d(1),d(2),n);
return;
%=======================================================================

%=======================================================================
function SP = init_sp(flags,VF,PG)
SP.VB       = spm_vol(flags.priors);
MM          = get_affine_mapping(VF,PG,flags.affreg);
%VF          = spm_vol(PF);
SP.MM       = MM*VF(1).mat;
SP.w        = 0.98;
return;
%=======================================================================

%=======================================================================
function s = get_sp(SP,x1,x2,x3)
[X1,X2,X3] = ndgrid(x1,x2,x3);
[Y1,Y2,Y3] = affine_transform(X1,X2,X3,SP.VB(1).mat\SP.MM);
w1  = SP.w;
w2  = (1-w1)/2;
s   = zeros([size(Y1),4]);
for i=1:3,
	s(:,:,i) = spm_sample_vol(SP.VB(i),Y1,Y2,Y3,1)*w1+w2;
end;
s(:,:,4:8)   = repmat(abs(1-sum(s(:,:,1:3),3))/5,[1 1 5]);
return;
%=======================================================================

%=======================================================================
function [g,w,c] = get_gwc(VF,BP,SP,CP,sums,wc)

if wc,
	VC = VF;
	for j=1:length(VF),
		[pth,nm,xt]    = fileparts(deblank(VF(j).fname));
		VC(j).fname    = fullfile(pth,['m' nm xt]);
		VC(j).descrip  = 'Bias corrected image';
	end;
	VC = spm_create_vol(VC);
end;

spm_progress_bar('Init',VF(1).dim(3),'Creating Segmented','planes completed');
x1 = 1:VF(1).dim(1);
x2 = 1:VF(1).dim(2);
x3 = 1:VF(1).dim(3);

g = uint8(0); g(VF(1).dim(1),VF(1).dim(2),VF(1).dim(3)) = 0;
w = uint8(0); w(VF(1).dim(1),VF(1).dim(2),VF(1).dim(3)) = 0;
c = uint8(0); c(VF(1).dim(1),VF(1).dim(2),VF(1).dim(3)) = 0;

for pp=1:length(x3),
	bf        = get_bp(BP,x1,x2,x3(pp));
	[raw,msk] = get_raw(VF,x1,x2,x3(pp));
	cor       = raw.*bf;
	if wc,
		for j=1:length(VC),
			VC(j) = spm_write_plane(VC(j),cor(:,:,j),pp);
		end;
	end;
	s         = get_sp(SP,x1,x2,x3(pp));
	p         = get_p(cor,msk,s,sums,CP,bf);
	g(:,:,pp) = uint8(round(p(:,:,1)*255));
	w(:,:,pp) = uint8(round(p(:,:,2)*255));
	c(:,:,pp) = uint8(round(p(:,:,3)*255));

	spm_progress_bar('Set',pp);
end;
spm_progress_bar('Clear');

if wc, spm_close_vol(VC); end;
return;
%=======================================================================

%=======================================================================
function [g,w,c] = clean_gwc(g,w,c)
b    = w;
b(1) = w(1);

% Build a 3x3x3 seperable smoothing kernel
%-----------------------------------------------------------------------
kx=[0.75 1 0.75];
ky=[0.75 1 0.75];
kz=[0.75 1 0.75];
sm=sum(kron(kron(kz,ky),kx))^(1/3);
kx=kx/sm; ky=ky/sm; kz=kz/sm;

% Erosions and conditional dilations
%-----------------------------------------------------------------------
niter = 32;
spm_progress_bar('Init',niter,'Extracting Brain','Iterations completed');
for j=1:niter,
	if j>2, th=0.15; else th=0.6; end; % Dilate after two its of erosion.
	for i=1:size(b,3),
		gp = double(g(:,:,i));
		wp = double(w(:,:,i));
		bp = double(b(:,:,i))/255;
		bp = (bp>th).*(wp+gp);
		b(:,:,i) = uint8(round(bp));
	end;
	spm_conv_vol(b,b,kx,ky,kz,-[1 1 1]);
	spm_progress_bar('Set',j);
end;
th = 0.05;
for i=1:size(b,3),
	gp       = double(g(:,:,i))/255;
	wp       = double(w(:,:,i))/255;
	cp       = double(c(:,:,i))/255;
	bp       = double(b(:,:,i))/255;
	bp       = ((bp>th).*(wp+gp))>th;
	g(:,:,i) = uint8(round(255*gp.*bp./(gp+wp+cp+eps)));
	w(:,:,i) = uint8(round(255*wp.*bp./(gp+wp+cp+eps)));
	c(:,:,i) = uint8(round(255*(cp.*bp./(gp+wp+cp+eps)+cp.*(1-bp))));
end;
spm_progress_bar('Clear');
return;


text:
function fx = p06_fx ( x )

%*****************************************************************************80
%
%% P06_FX evaluates exp ( x ) - 2 - 1 / ( 10 * x )^2 + 2 / ( 100 * x )^3.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    07 July 2011
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(*), the point at which F is to be evaluated.
%
%    Output, real FX(*), the value of the function at X.
%
  fx = exp ( x ) - 2.0 - 1.0 ./ ( 10.0 * x ).^2 + 2.0 ./ ( 100.0 * x ).^3;

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function [features, num_of_success] = extract_features(db, func)
	% db - a database to extract the features from
	% func - a function that extracts the features for a single database entry

    TRIM_SILENCE = true;
    GYRO = true;

    NUM_OF_ENTRIES = length(db);
	GYRO_DIM = 1;
    FS = 8000;

    features = {};
	num_of_success = 0;

	for k = 1:NUM_OF_ENTRIES
		try
			[wavdata, samp_rate] = read(db, k);
			wavdata = (wavdata{1});
            if GYRO
                wavdata = wavdata(:, GYRO_DIM);
            end
            wavdata = resample(wavdata, FS, samp_rate);
            fs = FS;

            if TRIM_SILENCE
                if GYRO
                    % cut 0.5 second at the beginning and a second from the end
                    wavdata = wavdata(4000:end-8000);
                end
                [wavdata, nseg] = get_voiced_segments(wavdata, fs);
            else
                nseg = 1;
            end

            if nseg > 0
                num_of_success = num_of_success + 1;
                features{num_of_success} = func(wavdata, fs);
            end
		catch ME
			% print error source
			ME.stack(1)
		end
	end
end

author: yueyuzhao

repo: gyrophone

save_path: github-repos/MATLAB/yueyuzhao-gyrophone


text:
function clqs = compute_minimal_interface(intra, inter)

int = compute_fwd_interface(intra, inter);
ss  = length(intra);
Z = zeros(ss);
dag = [intra inter;
       Z     intra];
G = moralize(dag);
intra2 = G(1:ss,1:ss);
inter2 = G(1:ss,(1:ss)+ss);
G = unroll_dbn_topology(intra2, inter2, ss);
T = ss;
last_slice = (1:ss) + (T-1)*ss;
G = (G + G')/2; % mk symmetric
G2 = (expm(full(G)) > 0); % closure of graph
G3 = G2(last_slice, last_slice);
[c,v] = scc(G3); % connected components
ncomp = size(v,1);
clqs = cell(1,ncomp);
for i=1:ncomp
  ndx = find(v(i,:)>0);
  clqs{i} = v(i,ndx);
  clqs{i} = myintersect(clqs{i}, int);
end

author: bayesnet

repo: bnt

save_path: github-repos/MATLAB/bayesnet-bnt


text:
function [volToEval,maskBoundingBox3M,gridS,paramS,diagS] = ...
    preProcessForRadiomics(scanNum,structNum, paramS, planC)
% preProcessForRadiomics.m
% Pre-process image for radiomics feature extraction. Includes
% perturbation, resampling, cropping, and intensity-thresholding.
%
% AI 3/28/19

%% Get scan array, mask and x,y,z grid for reslicing
if numel(scanNum) == 1

    if ~exist('planC','var')
        global planC
    end

    indexS = planC{end};

    % Get scan
    scanArray3M = getScanArray(planC{indexS.scan}(scanNum));
    scanArray3M = double(scanArray3M) - ...
        planC{indexS.scan}(scanNum).scanInfo(1).CTOffset;

    scanSiz = size(scanArray3M);

    if numel(structNum) == 1
        % Get structure mask
        [rasterSegments] = getRasterSegments(structNum,planC);
        [maskUniq3M, uniqueSlices] = rasterToMask(rasterSegments, scanNum, planC);
        mask3M = false(scanSiz);
        mask3M(:,:,uniqueSlices) = maskUniq3M;
    else
        mask3M = structNum;
    end

    % Get x,y,z grid for reslicing and calculating the shape features
    [xValsV, yValsV, zValsV] = getScanXYZVals(planC{indexS.scan}(scanNum));
    if yValsV(1) > yValsV(2)
        yValsV = fliplr(yValsV);
    end
    %zValsV = zValsV(uniqueSlices);

%     % Get image types with various parameters
%     fieldNamC = fieldnames(paramS.imageType);
%     for iImg = 1:length(fieldNamC)
%         for iFilt = 1:length(paramS.imageType.(fieldNamC{iImg}))
%             imageType = fieldNamC{iImg};
%             if strcmpi(imageType,'SUV')
%                 scanInfoS = planC{indexS.scan}(scanNum).scanInfo;
%                 paramS.imageType.(fieldNamC{iImg})(iFilt).scanInfoS = scanInfoS;
%             end
%         end
%     end

else

    scanArray3M = scanNum;
    mask3M = structNum;

    % Assume xValsV,yValsV,zValsV increase monotonically.
    xValsV = planC{1};
    yValsV = planC{2};
    zValsV = planC{3};

end

if isempty(scanArray3M)
    volToEval = [];
    maskBoundingBox3M = [];
    gridS = [];
    return;
end

% Pixelspacing (dx,dy,dz)
PixelSpacingX = abs(xValsV(1) - xValsV(2));
PixelSpacingY = abs(yValsV(1) - yValsV(2));
PixelSpacingZ = abs(zValsV(1) - zValsV(2));

diagS.numVoxelsOrig = sum(mask3M(:));

%% Apply global settings
whichFeatS = paramS.whichFeatS;

perturbX = 0;
perturbY = 0;
perturbZ = 0;

%--- 1. Perturbation ---
if whichFeatS.perturbation.flag

    [scanArray3M,mask3M] = perturbImageAndSeg(scanArray3M,mask3M,planC,...
        scanNum,whichFeatS.perturbation.sequence,...
        whichFeatS.perturbation.angle2DStdDeg,...
        whichFeatS.perturbation.volScaleStdFraction,...
        whichFeatS.perturbation.superPixVol);

    % Get grid perturbation deltas
    if ismember('T',whichFeatS.perturbation.sequence)
        perturbFractionV = [-0.75,-0.25,0.25,0.75];
        perturbFractionV(randsample(4,1));
        perturbX = PixelSpacingX*perturbFractionV(randsample(4,1));
        perturbY = PixelSpacingY*perturbFractionV(randsample(4,1));
        perturbZ = PixelSpacingZ*perturbFractionV(randsample(4,1));
    end

end

%--- 2. Crop scan around mask and pad as required ---
padScaleX = 1;
padScaleY = 1;
padScaleZ = 1;
if whichFeatS.resample.flag
    dx = median(diff(xValsV));
    dy = median(diff(yValsV));
    dz = median(diff(zValsV));
    padScaleX = ceil(whichFeatS.resample.resolutionXCm/dx);
    padScaleY = ceil(whichFeatS.resample.resolutionYCm/dy);
    padScaleZ = ceil(whichFeatS.resample.resolutionZCm/dz);
end

%---For IBSI2 ------
cropForResampling = 0;
padSizV = [];
padMethod = 'none';
%------------------
if cropForResampling
    %Default:Pad by 5 voxels (original image intensities) before resampling
    padMethod = 'expand';
    padSizV = [5,5,5];
    if whichFeatS.padding.flag
        if isfield(whichFeatS.padding,'size')
            filtPadSizV = whichFeatS.padding.size;
            filtPadSizV = reshape(filtPadSizV,1,[]);
            if length(filtPadSizV)==2
                filtPadSizV = [filtPadSizV,0];
            end
            filtPadSizV = filtPadSizV.*[padScaleX padScaleY padScaleZ];
            repIdxV = filtPadSizV > padSizV;
            padSizV(repIdxV) = filtPadSizV(repIdxV);
            padMethod = whichFeatS.padding.method;
        end
    end
end

% Crop to ROI and pad
scanArray3M = double(scanArray3M);
[padScanBoundsForResamp3M,padMaskBoundsForResamp3M,outLimitsV] = ...
    padScan(scanArray3M,mask3M,padMethod,padSizV,cropForResampling);
xValsV = xValsV(outLimitsV(3):outLimitsV(4));
yValsV = yValsV(outLimitsV(1):outLimitsV(2));
zValsV = zValsV(outLimitsV(5):outLimitsV(6));
slcIndV = outLimitsV(5):outLimitsV(6);

%--- 3. Resampling ---

%Get resampling method
if whichFeatS.resample.flag
    % Pixelspacing (dx,dy,dz) after resampling
    if ~isempty(whichFeatS.resample.resolutionXCm)
        PixelSpacingX = whichFeatS.resample.resolutionXCm;
    end
    if ~isempty(whichFeatS.resample.resolutionYCm)
        PixelSpacingY = whichFeatS.resample.resolutionYCm;
    end
    if ~isempty(whichFeatS.resample.resolutionZCm)
        PixelSpacingZ = whichFeatS.resample.resolutionZCm;
    end

    % Interpolate using the method defined in settings file
    roiInterpMethod = 'linear';
    scanInterpMethod = whichFeatS.resample.interpMethod;
    extrapVal = 0;
    gridResampleMethod = 'center';

    % Interpolate using the method defined in settings file
    origVolToEval = padScanBoundsForResamp3M;
    origMask = padMaskBoundsForResamp3M;
    outputResV = [PixelSpacingX,PixelSpacingY,PixelSpacingZ];
    originV = [xValsV(1),yValsV(1),zValsV(end)];

    %Get resampling grid
    [xResampleV,yResampleV,zResampleV] = getResampledGrid(outputResV,...
        xValsV,yValsV,zValsV,originV,gridResampleMethod,...
        [perturbX,perturbY,perturbZ]);

    %Resample scan
    resampScanBounds3M = imgResample3d(origVolToEval,xValsV,yValsV,zValsV,...
        xResampleV,yResampleV,zResampleV,scanInterpMethod,extrapVal);
    %Option to round
    if isfield(whichFeatS.resample,'intensityRounding') && ...
            strcmpi(whichFeatS.resample.intensityRounding,'on')
        resampScanBounds3M = round(resampScanBounds3M);
    end

    %Resample mask
    resampMaskBounds3M = imgResample3d(single(origMask),xValsV,yValsV,...
        zValsV,xResampleV,yResampleV,zResampleV,roiInterpMethod) >= 0.5;

    newSlcIndV = zeros(1,length(zResampleV));
    for iSlc = 1:length(zResampleV)
        newSlcIndV(iSlc) = findnearest(zValsV, zResampleV(iSlc));
    end

else
    resampScanBounds3M = padScanBoundsForResamp3M;
    resampMaskBounds3M = padMaskBoundsForResamp3M;
    xResampleV = xValsV;
    yResampleV = yValsV;
    zResampleV = zValsV;
    newSlcIndV = 1:length(slcIndV);
    dx = abs(median(diff(xResampleV)));
    dy = abs(median(diff(yResampleV)));
    dz = abs(median(diff(zResampleV)));
    outputResV = [dx dy dz];
end

%Apply padding as required for convolutional filtering
filtPadMethod = 'none';
filtPadSizeV = [0 0 0];
cropFlag = 1;  %Default
if ~cropForResampling
    if whichFeatS.padding.flag
        if isfield(whichFeatS.padding,'cropToMaskBounds')
            cropFlag = strcmpi(whichFeatS.padding.cropToMaskBounds,'yes');
        end
        if isfield(whichFeatS.padding,'method') && ...
                ~strcmpi(whichFeatS.padding.method,'none')
            filtPadMethod = whichFeatS.padding.method;
            filtPadSizeV = reshape(whichFeatS.padding.size,1,[]);
            if length(filtPadSizeV)==2
                filtPadSizeV = [filtPadSizeV,0];
            end
        end
    end
    [volToEval,maskBoundingBox3M,outLimitsV] = padScan(resampScanBounds3M,...
        resampMaskBounds3M,filtPadMethod,filtPadSizeV,cropFlag);

    %Extend resampling grid if padding original image (cropFlag=0)
    if outLimitsV(1)<1
        numPad = (1-outLimitsV(1));
        yExtendV = yResampleV(1)-(1:numPad)*outputResV(2);
        yResampleV = [yExtendV,yResampleV];
        outLimitsV(1) = outLimitsV(1)+numPad;
        outLimitsV(2) = outLimitsV(2)+numPad;
    end
    if outLimitsV(3)<1
        numPad = (1-outLimitsV(3));
        xExtendV = xResampleV(1)-(1:numPad)*outputResV(1);
        xResampleV = [xExtendV,xResampleV];
        outLimitsV(3) = outLimitsV(3)+numPad;
        outLimitsV(4) = outLimitsV(4)+numPad;
    end
    if outLimitsV(5)<1
        numPad = (1-outLimitsV(5));
        zExtendV = zResampleV(1)-(1:numPad)*outputResV(3);
        zResampleV = [zExtendV,zResampleV];
        outLimitsV(5) = outLimitsV(5)+numPad;
        outLimitsV(6) = outLimitsV(6)+numPad;
    end
    if outLimitsV(2)>length(yResampleV)
        numPad = (outLimitsV(2)-length(yResampleV));
        yExtendV = yResampleV(end)+(1:numPad)*outputResV(2);
        yResampleV = [yResampleV,yExtendV];
    end
    if outLimitsV(4)>length(xResampleV)
        numPad = (outLimitsV(4)-length(xResampleV));
        xExtendV = xResampleV(end)+(1:numPad)*outputResV(1);
        xResampleV = [xResampleV,xExtendV];
    end
    if outLimitsV(6)>length(zResampleV)
        numPad = (outLimitsV(6)-length(zResampleV));
        zExtendV = zResampleV(end)+(1:numPad)*outputResV(3);
        zResampleV = [zResampleV,zExtendV];
    end
    xResampleV = xResampleV(outLimitsV(3):outLimitsV(4));
    yResampleV = yResampleV(outLimitsV(1):outLimitsV(2));
    zResampleV = zResampleV(outLimitsV(5):outLimitsV(6));
else
    %     resampSizeV = size(resampScanBounds3M);
    %     volToEval = resampScanBounds3M(padSizV(1)+1:resampSizeV(1)-padSizV(1),...
    %         padSizV(2)+1:resampSizeV(2)-padSizV(2),...
    %         padSizV(3)+1:resampSizeV(3)-padSizV(3));
    %     maskBoundingBox3M = resampMaskBounds3M(padSizV(1)+1:resampSizeV(1)-padSizV(1),...
    %         padSizV(2)+1:resampSizeV(2)-padSizV(2),...
    %         padSizV(3)+1:resampSizeV(3)-padSizV(3));
    %    resampSizeV = size(resampScanBounds3M);
    volToEval = resampScanBounds3M;
    maskBoundingBox3M = resampMaskBounds3M;
end

%--- 4. Ignore voxels below and above cutoffs, if defined ----
minSegThreshold = [];
maxSegThreshold = [];
if isfield(paramS,'textureParamS')
    if isfield(paramS.textureParamS,'minSegThreshold')
        minSegThreshold = paramS.textureParamS.minSegThreshold;
    end
    if isfield(paramS.textureParamS,'maxSegThreshold')
        maxSegThreshold = paramS.textureParamS.maxSegThreshold;
    end
    if ~isempty(minSegThreshold)
        maskBoundingBox3M(volToEval < minSegThreshold) = 0;
    end
    if ~isempty(maxSegThreshold)
        maskBoundingBox3M(volToEval > maxSegThreshold) = 0;
    end
end
%volToEval(~maskBoundingBox3M) = NaN;
volV = volToEval(maskBoundingBox3M);
diagS.numVoxelsInterpReseg = sum(maskBoundingBox3M(:));
diagS.MeanIntensityInterpReseg = mean(volV);
diagS.MaxIntensityInterpReseg = max(volV);
diagS.MinIntensityInterpReseg = min(volV);

% Return grid and pixel-spacing
gridS.xValsV = xResampleV;
gridS.yValsV = yResampleV;
gridS.zValsV = zResampleV;
gridS.PixelSpacingV = [PixelSpacingX,PixelSpacingY,PixelSpacingZ];

% Pass scanInfo as an additional parameter for imageType = SUV
if numel(scanNum) == 1
    % Get image types with various parameters
    fieldNamC = fieldnames(paramS.imageType);
    for iImg = 1:length(fieldNamC)
        for iFilt = 1:length(paramS.imageType.(fieldNamC{iImg}))
            imageType = fieldNamC{iImg};
            if strcmpi(imageType,'SUV')
                scanInfoS = planC{indexS.scan}(scanNum).scanInfo;
                scanInfoS = scanInfoS(slcIndV); % slices for the cropped volume
                scanInfoS = scanInfoS(newSlcIndV); % resampled slices
                paramS.imageType.(fieldNamC{iImg})(iFilt).scanInfoS = scanInfoS;
            end
        end
    end
end

end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function value = p12_f ( dim_num, point_num, x )

%*****************************************************************************80
%
%% P12_F evaluates the integrand for problem 12.
%
%  Discussion:
%
%    The highly oscillatory nature of the integrand makes this
%    a difficult and perhaps even dubious test.
%
%  Dimension:
%
%    DIM_NUM arbitrary.
%
%  Region:
%
%    0 <= X(1:DIM_NUM) <= 1
%
%  Integrand:
%
%    product ( 1 <= i <= dim_num ) ( i * cos ( i * x(i) ) )
%
%  Exact Integral:
%
%    product ( 1 <= I <= DIM_NUM ) sin ( i )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    02 June 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Paul Bratley, Bennett Fox, Harald Niederreiter,
%    Implementation and Tests of Low-Discrepancy Sequences,
%    ACM Transactions on Modeling and Computer Simulation,
%    Volume 2, Number 3, July 1992, pages 195-213.
%
%  Parameters:
%
%    Input, integer DIM_NUM, the dimension of the argument.
%
%    Input, integer POINT_NUM, the number of points.
%
%    Input, real X(DIM_NUM,POINT_NUM), the evaluation points.
%
%    Output, real VALUE(POINT_NUM), the integrand values.
%
  value(1:point_num) = 1.0;

  for point = 1 : point_num
    for dim = 1 : dim_num
      value(point) = value(point) * dim * cos ( dim * x(dim,point) );
    end
  end

  p12_i4 ( 'I', '#', point_num );

  return
end


text:
%kvcentfilt 'Center Filter'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vcentfilt.pane file
%
% Parameters:
% InputFile: i 'Input Image', required: 'input image'
% InputFile: str 'Struct Element', required: 'Structuring Element'
% OutputFile: o 'Output Image ', required: 'resulting output image'
%
% Example: o = kvcentfilt({i, str}, {'i','';'str','';'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vcentfilt - Center Filter
%
%  DESCRIPTION
% .I vcentfilt
% will do the center filter of an input image by a set structuring element.
%
% The output image data type will be the same as that of the highest data type
% of the input image.
%
%
%
%  EXAMPLES
% vcentfilt -i ball.xv -str b272.str -o outimage.xv
%
% Will do the center filter of the image "ball.xv" by the structuring set b272.str, with the resulting image written to "outimage.xv".
%
%  "SEE ALSO"
%
%  RESTRICTIONS
% .I vcentfilt
% can be defined for all data types supported by Khoros, but at the moment it has been implemented just for the bit and unsigned char types.
% The structuring elements are subsets of the 3x3 matrix and the origin is always at the center of this matrix.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993-1997 Junior Barrera, Roberto Lotufo.  All rights reserved.
%

function varargout = kvcentfilt(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvcentfilt(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'str', '__input';'o', '__output'};
maxval={0,0,0};
minval={0,0,0};
istoggle=[0,0,0];
was_set=istoggle * 0;
paramtype={'InputFile','InputFile','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vcentfilt"  '],Inputs,narglist);

author: aludnam

repo: MATLAB

save_path: github-repos/MATLAB/aludnam-MATLAB


text:
%kkbitwise 'Dual Operand Bitwise Operations'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros kbitwise.pane file
%
% Parameters:
% InputFile: i1 'Input 1', required: 'first input data object'
% OutputFile: o 'Output', required: 'resulting output data object'
% InputFile: i2 'Input 2', optional: 'second input data object'
% Integer: real 'Constant ', default: 1: 'constant value'
% Toggle: and 'AND', default: 0: 'Bitwise AND: o = i1 & (i2 or val)'
% Toggle: ls 'Left Shift', default: 0: 'Bitwise Left Shift: o = i1 << by (i2 or val)'
% Toggle: or 'OR', default: 0: 'Bitwise OR: o = i1 | (i2 or val)'
% Toggle: rs 'Unsigned Right Shift', default: 0: 'Bitwise Right Shift: o = i1 >> by (i2 or val)'
% Toggle: xor 'XOR', default: 0: 'Bitwise XOR: o = i1 ^ (i2 or val)'
% Toggle: srs 'Signed Right Shift', default: 0: 'Bitwise Signed Right Shift by (i2 or val)'
% Toggle: nand 'NAND', default: 0: 'Bitwise NAND: o = (~i1) | (~(i2 or val))'
% Toggle: rr 'Right Rotate', default: 0: 'Bitwise Right Rotate by (i2 or val)'
% Toggle: nor 'NOR', default: 0: 'Bitwise NOR: o = (~i1) & (~(i2 or val))'
% Toggle: lr 'Left Rotate', default: 0: 'Bitwise Left Rotate by (i2 or val)'
% Toggle: andrev '1 AND (NOT 2)', default: 0: 'Bitwise AND Reverse: o = i1 & (~(i2 or val))'
% Toggle: andinv '(NOT 1) AND 2', default: 0: 'Bitwise AND Inverted: o = (~i1) & (i2 or val)'
% Toggle: orrev '1 OR (NOT 2)', default: 0: 'Bitwise OR Reverse: o = i1 | (~(i2 or val))'
% Toggle: orinv '(NOT 1) OR 2', default: 0: 'Bitwise OR Inverted: o = (~i1) | i2'
%
% Example: o = kkbitwise({i1, i2}, {'i1','';'o','';'i2','';'real',1;'and',0;'ls',0;'or',0;'rs',0;'xor',0;'srs',0;'nand',0;'rr',0;'nor',0;'lr',0;'andrev',0;'andinv',0;'orrev',0;'orinv',0})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% kbitwise - Dual Operand Bitwise Operations
%
%  DESCRIPTION
% The bitwise operator, kbitwise, performs a the specified bitwise operation
% between each point in \fBInput 1" (i1) and each corresponding point in
% \fBInput 2" (i2), if the second input is specified.  If Input 2 is not
% specified, the operation will be performed using the \fBConstant" value
% (real).
%
% Therefore, if kbitwise were called with two input objects:
%
% 	kbitwise -i1 image:ball -i2 image:moon -o out -and
%
% then the following operation would be performed
%
% 	out = i1 & i2
%
% Whereas if kbitwise were called with a single data object and a constant:
%
% 	kbitwise -i1 image:ball -real 10 -o out -and
%
% then the following would be performed:
%
% 	o = i1 & 10
%
%
% The supported functions are:
%
%
% 	AND (and)                   o = i1 & (i2 or val)
% 	OR   (or)                   o = i1 | (i2 or val)
% 	Exclusive Or (xor)          o = i1 ^ (i2 or val)
% 	NAND (nand)                 o = (~i1) | (~(i2 or val))
% 	NOR  (nor)                  o = (~i1) & (~(i2 or val))
% 	AND Reverse  (andrev)       o = i1 & (~(i2 or val))
% 	AND Inverted (andinv)       o = (~i1) & (i2 or val)
% 	OR Reverse  (orrev)         o = i1 | (~(i2 or val))
% 	OR Inverted (orinv)         o = (~i1) | i2
% 	Unsigned Right Shift (rs)   o = i1 >> by (i2 or val)
% 	Unsigned Left  Shift (ls)   o = i1 << by (i2 or val)
%
%
%
%  "Data Type - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_1input
%
%  "Data Type - Two Input Objects" 5
% .cI $DATAMANIP/repos/shared/man/sections/value_type_2input
% The bitwise operators will not operate on float or complex data types.
%
%  "Map Data" 5
% The bitwise operations have not been written to be fully polymorphic yet.
% They does not
% check for map data, and will therefore always operate on the value data,
% even if a map exists.  This will most likely corrupt indexing into the map.
% In the case of a single input, it is recommended to use the "Copy to
% Value" (kcptoval) segment operator to temporarily move the map data into
% the value segment, run kbitwise on the data, then move it back to the
% map with the "Copy from Value" (kcpfromval) operator.
% When operating with two input objects, where at least one of the objects
% contains map data, data should be mapped prior to running the bitwise
% operator.  The "Map Data" operator (kmapdata) can be used to perform
% the mapping operation.
%
%  "Validity Mask - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/mask_1input
%
%  "Validity Mask - Two Input Objects" 5
% Masking has not been implemented yet for the bitwise operators.  Therefore,
% only the mask from the first input object will be transferred to the output.
%
%  "Input Objects of Different Sizes" 5
% .cI $DATAMANIP/repos/shared/man/sections/resize_2input
% The value used to pad the data when the input files are not the same size is
% zero.
%
%  "Explicit Location and Time Data - Single Input" 5
% .cI $DATAMANIP/repos/shared/man/sections/loc_and_time_1input
%
%  "Explicit Location and Time Data - Two Input Objects" 5
% The bitwise operations have not been extended to understand location and
% time data.  Therefore, only location and time data present in the first input
% object will be transferred to the output.
%
%
%
%  EXAMPLES
%
%  "SEE ALSO"
%
%  RESTRICTIONS
%
% THIS ROUTINE is still under development, and has not been modified to fully
% support the polymorphic data model.  See paragraphs above for discussions
% concerning the polymorphic data segments.
%
% All processing is currently performed as UNSIGNED long, so the operations
% may not work properly on negative values.
%
% The bitwise operators will not operate on float or complex data types.
%
%  REFERENCES
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
%

function varargout = kkbitwise(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kkbitwise(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i1', '__input';'o', '__output';'i2', '__input';'real', 1;'and', 0;'ls', 0;'or', 0;'rs', 0;'xor', 0;'srs', 0;'nand', 0;'rr', 0;'nor', 0;'lr', 0;'andrev', 0;'andinv', 0;'orrev', 0;'orinv', 0};
maxval={0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
minval={0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
istoggle=[0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1];
was_set=istoggle * 0;
paramtype={'InputFile','OutputFile','InputFile','Integer','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle','Toggle'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs))
        error(['Input ' narglist{jj,1} ' has no corresponding input!']);
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout)
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']);
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in')
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = '';
        end;
        narglist{jj,2} = '';
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = '';
          narglist{jj,2} = '';
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs))
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']);
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs))
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1)
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']);
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs)
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'kbitwise"  '],Inputs,narglist);


text:
function T_new = ExplicitScheme(T_old,alpha,dx,dy,dt,x_intervals,y_intervals)

T_new = T_old;

for x_index = 2:1:(x_intervals-1)
    for y_index = 2:1:(y_intervals-1)

         T_new(x_index,y_index) = T_old(x_index,y_index) + (alpha*dt)*((T_old(x_index+1,y_index) - 2*T_old(x_index,y_index)+ T_old(x_index-1,y_index))/dx^2 + (T_old(x_index,y_index+1) - 2*T_old(x_index,y_index)+ T_old(x_index,y_index-1))/dy^2);

    end
end

end

author: Sable

repo: mcbench-benchmarks


text:
function perm = shuffle_text_file(filename,dest_file,perm)

the_lines = fileread(filename);
the_lines = regexp(the_lines,'(.*?)\n','tokens');

if(~exist('perm','var'))
    perm = randperm(length(the_lines));
end

dest_lines = the_lines(perm);

fp = fopen(dest_file,'wt');
for i = 1 : length(dest_lines)
    fprintf(fp,'%s\n',dest_lines{i}{1});
end
fclose(fp);

author: lmb-freiburg

repo: orion

save_path: github-repos/MATLAB/lmb-freiburg-orion


text:
function log_evidence = ldae_chibms(words, topics, topic_prior, ms_iters)
%LDAE_CHIBMS Approximate evidence for LDA using Murray & Salakhutdinov's Chib-style method
%
% log_evidence = ldae_chibms(words, topics, topic_prior);
%
% Inputs:
%             words 1xNd
%            topics TxV each row is a distribution over a vocabulary of size V
%       topic_prior 1xT parameters of Dirichlet from which document topic vector is drawn
%          ms_iters 1x1 Default: 1000
%
% Outputs:
%     log_evidence  1x1

% Iain Murray, January 2009

BURN_ITERS = 3;

if ~exist('ms_iters', 'var')
    ms_iters = 1000;
end

[T, V] = size(topics);
Nd = length(words);

% Sanity checking input sizes
assert(isvector(topic_prior));
assert(T == length(topic_prior));
assert(isvector(words));

topic_prior = topic_prior(:)';
topic_alpha = sum(topic_prior);
%topic_mean = topic_prior / topic_alpha;

% Assign latents to words in isolation as a simple initialization
Nz = zeros(1, T);
for t = 1:Nd
    pz = topics(:, words(t))' .* topic_prior;
    zz(t) = discreternd(1, pz);
    Nz(zz(t)) = Nz(zz(t)) + 1;
end

% Run some sweeps of Gibbs sampling
for sweeps = 1:BURN_ITERS
    for t = 1:Nd
        Nz(zz(t)) = Nz(zz(t)) - 1;
        pz = topics(:, words(t))' .* (Nz + topic_prior);
        zz(t) = discreternd(1, pz);
        Nz(zz(t)) = Nz(zz(t)) + 1;
    end
end

% Find local optimim to use as z^*, "iterative conditional modes"
% But don't spend forever on this, bail out if necessary
for i = 1:12
    old_zz = zz;
    for t = 1:Nd
        Nz(zz(t)) = Nz(zz(t)) - 1;
        pz = topics(:, words(t))' .* (Nz + topic_prior);
        [dummy, zz(t)] = max(pz);
        Nz(zz(t)) = Nz(zz(t)) + 1;
    end
    if ~isequal(old_zz, zz)
        break;
    end
end

% Run Murray & Salakhutdinov algorithm
zstar = zz;
log_Tvals = zeros(ms_iters, 1);
log_Tprob = @(zto, zfrom, Nzfrom) log_Tprob_base(zto, zfrom, Nzfrom, words, topics, topic_prior);
% draw starting position
ss = ceil(rand() * ms_iters);
% Draw z^(s)
for t = Nd:-1:1
    Nz(zz(t)) = Nz(zz(t)) - 1;
    pz = topics(:, words(t))' .* (Nz + topic_prior);
    zz(t) = discreternd(1, pz);
    Nz(zz(t)) = Nz(zz(t)) + 1;
end
zs = zz;
log_Tvals(ss) = log_Tprob(zstar, zz, Nz);
% Draw forward stuff
for sprime = (ss+1):ms_iters
    for t = 1:Nd
        Nz(zz(t)) = Nz(zz(t)) - 1;
        pz = topics(:, words(t))' .* (Nz + topic_prior);
        zz(t) = discreternd(1, pz);
        Nz(zz(t)) = Nz(zz(t)) + 1;
    end
    log_Tvals(sprime) = log_Tprob(zstar, zz, Nz);
end
% Draw backward stuff
for sprime = (ss-1):-1:1
    for t = Nd:-1:1
        Nz(zz(t)) = Nz(zz(t)) - 1;
        pz = topics(:, words(t))' .* (Nz + topic_prior);
        zz(t) = discreternd(1, pz);
        Nz(zz(t)) = Nz(zz(t)) + 1;
    end
    log_Tvals(sprime) = log_Tprob(zstar, zz, Nz);
end
% Final estimate
Nkstar = histc(zstar, 1:T); Nkstar = Nkstar(:)'; % 1xT
log_pz = sum(gammaln(Nkstar + topic_prior)) + gammaln(topic_alpha) ...
        - sum(gammaln(topic_prior)) - gammaln(Nd + topic_alpha);
log_w_given_z = 0;
for t = 1:Nd
    log_w_given_z = log_w_given_z + log(topics(zstar(t), words(t)));
end
log_joint = log_pz + log_w_given_z;
log_evidence = log_joint - (logsumexp(log_Tvals) - log(ms_iters));

function lp = log_Tprob_base(zto, zfrom, Nz, words, topics, topic_prior)
Nd = length(words);
lp = 0;
for t = 1:Nd
    Nz(zfrom(t)) = Nz(zfrom(t)) - 1;
    pz = topics(:, words(t))' .* (Nz + topic_prior);
    pz = pz/sum(pz);
    lp = lp + log(pz(zto(t)));
    Nz(zto(t)) = Nz(zto(t)) + 1;
end

author: jacobeisenstein

repo: SAGE


text:
function scales = hsvargplvmShowScales(model, displ, varargin)
if nargin < 2 || isempty(displ)
    displ = true;
end

% Layers to visualise
layers = model.H:-1:1;

% Do not show individual scales in mult-output
skipComp = false;

displayAsMatrix = false;

if nargin > 2
    if ~isempty(varargin{1}), layers = varargin{1}; end
    if length(varargin) > 1 && ~isempty(varargin{2})
        skipComp = varargin{2};
    end
    if length(varargin) > 2 && ~isempty(varargin{3})
        displayAsMatrix = varargin{3};
    end
end

%{
if ~model.multOutput
    for h=1:model.H
        if displ
            subplot(model.H,1,h)
        end
        vargplvmShowScales(model.layer{h}.comp{1},displ); title(num2str(h));
    end
    return
end
%}

for hCount=1:length(layers)
    h = layers(hCount);
    scalesAll{h} = zeros(1, model.layer{h}.q);
    if model.layer{h}.M > 10 && displ && ~displayAsMatrix
        for i=1:model.layer{h}.M
            sc = vargplvmShowScales(model.layer{h}.comp{i}, ~skipComp);
            scalesAll{h} = scalesAll{h} + sc;
            if ~skipComp
                title(['Scales for layer ' num2str(h) ', model ' num2str(i)])
                pause
            end
        end
        scalesAll{h} = scalesAll{h} ./ max(scalesAll{h});
    else
        if displ && ~displayAsMatrix
            if ~model.multOutput
                subplot(model.H,1,hCount)
            else
                figure
            end
        end
        scales{h} = svargplvmShowScales(model.layer{h}, (displ && ~displayAsMatrix));
        if displ && ~displayAsMatrix
            title(['Layer ' num2str(h)]);
        end
        if model.layer{h}.M < 2
            scalesAll{h} = scales{h}{1};
            scalesAll{h} = scalesAll{h} ./ max(scalesAll{h});
        end
    end

    if model.layer{h}.M > 10 && displ && ~displayAsMatrix
        bar(scalesAll{h}); title(['Normalised sum of scales for layer ' num2str(h)])
    end
end

if displayAsMatrix
    maxQ = length(scalesAll{1});
    for h = 2:model.H
        if length(scalesAll{h}) > maxQ
            maxQ = length(scalesAll{h});
        end
    end

    scalesAllMat = zeros(model.H, maxQ);
    for hh = 1:model.H
        %h = hh;             % This will put layer 1 on top
        h = model.H - hh +1; % This will put layer H on top
        for q = 1:maxQ
            if q <= length(scalesAll{hh})
                scalesAllMat(h,q) = scalesAll{hh}(q);
            else
                scalesAllMat(h,q) = NaN;
            end
        end
    end
    h=imagesc(scalesAllMat);
    set(h,'alphadata',~isnan(scalesAllMat))
    colorbar
end

author: SheffieldML

repo: deepGP

save_path: github-repos/MATLAB/SheffieldML-deepGP


text:
function [A_or,C_or,S_or,P_or,srt,srt_val] = order_ROIs(A,C,S,P,srt)

% ordering of the found components based on their maximum temporal
% activation and their size (through their l_inf norm)
% you can also pre-specify the ordering sequence

nA = full(sqrt(sum(A.^2)));
nr = length(nA);
A = A/spdiags(nA(:),0,nr,nr);
C = spdiags(nA(:),0,nr,nr)*C;
mA = sum(A.^4).^(1/4);
%sA = sum(A);
mC = max(C,[],2);
if ~exist('srt', 'var')||isempty(srt)
    [srt_val,srt] = sort(mC.*mA','descend');
end
A_or = A(:,srt);
C_or = C(srt,:);

if nargin < 4
    P_or = [];
else
    P_or = P;
    if isfield(P,'gn'); P_or.gn=P.gn(srt); end
    if isfield(P,'b'); P_or.b = cellfun(@times,P.b(srt),num2cell(nA(srt)'),'UniformOutput',false); end
    if isfield(P,'c1'); P_or.c1 = cellfun(@times,P.c1(srt),num2cell(nA(srt)'),'UniformOutput',false); end
    if isfield(P,'neuron_sn'); P_or.neuron_sn=num2cell(nA(srt)'.*cell2mat(P.neuron_sn(srt))); end
end

if nargin < 3 || isempty(S)
    S_or = [];
else
    S = spdiags(nA(:),0,nr,nr)*S;
    S_or = S(srt,:);
end

author: flatironinstitute

repo: CaImAn-MATLAB


text:
function outC = anonymize(inC,sensitiveString,replacement)
%outC = anonymize(inC,sensitiveString,replacement)
%Replace all field values in the cell array inC (usually a CERR planC archive)
%which contains the strings in 'sensitiveString' with the replacement string 'replacement'.
%This anonymizer function is useful for removing protected health information from CERR archives.
%Note:  *No warranty is expressed or implied*.  Use at your own risk.
%Created Dec 02, JOD.
%Example of use:
%str = '12383_4dec02'
%
%planC = anonymize(planC,'studyNumberOfOrigin',str)
%Normally one should loop over all cells with potential protected health information.
%The following fields should always be viewed suspiciously:
%PHI = {'studyNumberOfOrigin','patientName','caseNumber','archive'}.
%Last modified:  20 May 03, JOD, corrected loss of CERROptions cell.
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
% CERR is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
% without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
% See the GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with CERR.  If not, see <http://www.gnu.org/licenses/>.

%Loop over all the cells

for i = 1 : length(inC)

    if i ~= inC{end}.CERROptions
        if length(inC{i})~= 0
            entry = inC{i}(1);

            if ~isempty(entry)
                %Get all struct field names
                fields = subfieldnames(entry);
                %which have the string sensitiveString

                for j = 1 : length(fields)
                    name = fields{j};
                    indV = strfind(lower(name),lower(sensitiveString));
                    if ~isempty(indV)
                        if length(indV) ~= 1
                            error(['On field name:  ' name])
                        end
                        %anonymize!
                        disp(['Found Patient Name in ' num2str(name)])
                        %eval([ ['[' name ']'] '= deal(replacement);']);
                        for k = 1:2
                            [nameCell{k} name] = strtok(name,'.');
                        end
                        for k = 1 : length(inC{i})
                            if isempty(name(2:end))
                                eval(['inC{i}' '(',num2str(k),').' nameCell{2} '= deal(replacement);' ])
                            else
                                for l = 1 : length(eval([nameCell{1} '.' nameCell{2}]))
                                    eval(['inC{i}' '(',num2str(k),').' nameCell{2} '(',num2str(l),').' name(2:end) '= deal(replacement);' ])
                                end

                            end
                        end
                    end

                end
            end
        end
        outC{i} =inC{i};
    end
end
outC{inC{end}.CERROptions} = inC{inC{end}.CERROptions};  %Put in options


text:
%  Infrared small target detection via line-based reconstruction and entropy-induced suppression
clearvars;
close all;
clc;
flag = 1;
wflag = 1;
channel = 1;
readimg;
for kk = 1
    fold = '.\data\';% 27 images
    try
        img = imread([fold, num2str(kk), '.jpg']);
    catch
        img = imread([fold, num2str(kk), '.bmp']);
    end
    img = img(:,:,1);
    img = double(img);
    img = ones(500, 500);
    img(:, 1:250) = 100;
    img(:, 251:end) = 180;
    img(100:104, 100:104) = 100 +100* fspecial('gaussian', [5, 5]);
    img = img(:,:,1);
    [ re, re1, re2] = lrfunc(img);
    bw = bwfunc(re);
    if flag
        figure; imshow(uint8(img));
        figure; mesh(img);
        figure; imshow(re, []);
        figure; imshow(re1, []);
        figure; imshow(re2, []);
    end
    if wflag
        refold = '.\'
        if ~exist(refold, 'dir')
            mkdir(refold);
        end
        imwrite(uint8(img), [refold, num2str(kk), '1.png']);
        imwrite(uint8( mat2gray(re) * 255), [refold, num2str(kk), '2.png']);
        imwrite(bw, [refold, num2str(kk), '3.png']);
    end
end

author: daxjuanxiong

repo: infrared-small-target-detection


text:
function lagrange_nd_test10 ( )

%*****************************************************************************80
%
%% LAGRANGE_ND_TEST10 tests LAGRANGE_PARTIAL2 in 2D.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    29 January 2014
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'LAGRANGE_ND_TEST10\n' );
  fprintf ( 1, '  LAGRANGE_PARTIAL2 determines\n' );
  fprintf ( 1, '  the Lagrange interpolating polynomials L(x)\n' );
  fprintf ( 1, '  for ND points in D dimensions, assuming that\n' );
  fprintf ( 1, '  the number of points is less than or equal to\n' );
  fprintf ( 1, '  R = Pi(D,N), the number of monomials of degree N or less\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  For this example, the data points are the same as those\n' );
  fprintf ( 1, '  used by the level 2 Clenshaw Curtis sparse grid in 2D.\n' );

  d = 2;
  n = 4;
  r = mono_upto_enum ( d, n );
  nd = 13;
  xd = [ 0.0,  0.0;
        -1.0,  0.0;
         1.0,  0.0,
         0.0, -1.0,
         0.0,  1.0,
        -1.0,  1.0,
         1.0,  1.0,
        -1.0, -1.0,
         1.0, -1.0,
        -0.5,  0.0,
         0.0, -0.5,
         0.0, +0.5,
         0.5,  0.0 ]';

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Spatial dimension D = %d\n', d );
  fprintf ( 1, '  Maximum degree N = %d\n', n );
  fprintf ( 1, '  Number of monomials R = %d\n', r );
  fprintf ( 1, '  Number of data points ND = %d\n', nd );

  r8mat_transpose_print ( d, nd, xd, '  Data points XD:' );

  [ po, pc, pe ] = lagrange_partial2 ( d, n, r, nd, xd );
%
%  Print the polynomials.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Lagrange polynomials for XD data points:\n' );
  fprintf ( 1, '\n' );
  for i = 1 : nd
    o = po(i);
    label = sprintf ( '  P(%d)(x) =', i );
    polynomial_print ( d, o, pc(i,1:o), pe(i,1:o), label );
  end
%
%  Evaluate the polynomials at XD.
%
  value = zeros ( nd, nd );

  for j = 1 : nd
    o = po(j);
    label = sprintf ( '  P(%d)(x) =', j );
    value(1:nd,j) = polynomial_value ( d, o, pc(j,1:o), pe(j,1:o), nd, xd );
  end

  err = r8mat_is_identity ( nd, value );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Frobenius norm of Lagrange matrix error = %g\n', err );
%
%  Now evaluate a function at the data points.
%
  pd = zeros ( nd, 1 );
  for i = 1 : nd
    pd(i) = sin ( xd(1,i) ) * cos ( xd(2,i) );
  end
%
%  Compare exact function and interpolant at a grid of points.
%
  ni = 11 * 11;
  xyi = zeros ( 2, ni );

  k = 0;
  for j = 1 : 11
    for i = 1 : 11
      k = k + 1;
      xyi(1,k) = ( ( 11 - i     ) * ( - 1.0 )   ...
                 + (      i - 1 ) * ( + 1.0 ) ) ...
                 / ( 11     - 1 );
      xyi(2,k) = ( ( 11 - j     ) * ( - 1.0 )   ...
                 + (      j - 1 ) * ( + 1.0 ) ) ...
                 / ( 11     - 1 );
    end
  end

  pn = nd;
  zi = interpolant_value ( d, r, pn, po, pc, pe, pd, ni, xyi );

  err = 0.0;
  for k = 1 : ni
    f = sin ( xyi(1,k) ) * cos ( xyi(2,k) );
    err = max ( err, abs ( zi(k) - f ) );
  end
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Maximum absolute interpolant error on 11x11 grid = %g\n', err );

%
%  Plot the interpolant at a grid of points.
%
%  Reshuffling the shapes of MATLAB vectors is my definition of hell.
%
  x1 = linspace ( -1, +1, 11 );
  y1 = linspace ( -1, +1, 11 );
  [ x2, y2 ] = meshgrid ( x1, y1 );
  x3 = reshape ( x2, 1, 11 * 11 );
  y3 = reshape ( y2, 1, 11 * 11 );
  xyi = [ x3; y3 ];
  ni = 11 * 11;
  zi = interpolant_value ( d, r, pn, po, pc, pe, pd, ni, xyi );
  z2 = reshape ( zi, 11, 11 );

  surf ( x2, y2, z2 );
  title ( 'Lagrange interpolant to sin(x)*cos(y)', 'FontSize', 16 )
  xlabel ( '<---X--->', 'FontSize', 16 )
  ylabel ( '<---Y--->', 'FontSize', 16 )
  zlabel ( '<---F(X,Y)--->', 'FontSize', 16 )

  filename = 'test10_interpolant.png';
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Saving plot file in "%s"\n', filename );

  return
end


text:
function aic = ml_gmm_aic(Data, Priors,Mu,Sigma,cov_type)
%ML_GMM_AIC Akaike Information criterion
%
%  input ------------------------------------------------------------------
%
%   o Data:     D x N array representing N datapoints of D dimensions.
%
%   o Priors:   1 x K array representing the prior probabilities of the
%               K GMM components.
%   o Mu:       D x K array representing the centers of the K GMM components.
%
%   o Sigma:    D x D x K array representing the covariance matrices of the
%               K GMM components.
%
%   o cov_type: string, covariance type = 'full','diag' or 'iso'
%
%  output -----------------------------------------------------------------
%
%   o aic   : (1 x 1)
%
%
[D,~]       = size(Data);
K           = length(Priors);
num_param   = K-1 + K * D;

if strcmp(cov_type,'full') == true
    num_param = num_param + K * ( D * ( D - 1)/2 );
elseif strcmp(cov_type,'diag') == true
    num_param = num_param + K * D;
elseif strcmp(cov_type,'iso') == true
    num_param = num_param + K * 1;
else
   error(['no such covariance type: ' cov_type '  only full | diag | isot ']);
end

% compute the loglikelihood of the data given the model
% loglik: (N x 1)
loglik = ml_LogLikelihood_gmm(Data,Priors,Mu,Sigma);

aic = - 2 * loglik + 2 * num_param;

end

author: epfl-lasa


text:
function [d,fp,dt,tc,t]=readhtk(file)
%READHTK  read an HTK parameter file [D,FP,DT,TC,T]=(FILE)
%
% d is data, fp is frame period in seconds
% dt is data type, tc is full type code, t is a text version of the full typecode
% tc is the sum of the following values:
%			0		WAVEFORM
%			1		LPC
%			2		LPREFC
%			3		LPCEPSTRA
%			4		LPDELCEP
%			5		IREFC
%			6		MFCC
%			7		FBANK
%			8		MELSPEC
%			9		USER
%			10		DISCRETE
%                       11              PLP
%			64		-E		Includes energy terms
%			128	_N		Suppress absolute energy
%			256	_D		Include delta coefs
%			512	_A		Include acceleration coefs
%			1024	_C		Compressed
%			2048	_Z		Zero mean static coefs
%			4096	_K		CRC checksum (not implemented yet)
%			8192	_0		Include 0'th cepstral coef

%      Copyright (C) Mike Brookes 1997
%
%      This version modified to read HTK's compressed feature files
%      2005-05-18 dpwe@ee.columbia.edu
%
%   VOICEBOX home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This program is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation; either version 2 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You can obtain a copy of the GNU General Public License from
%   ftp://prep.ai.mit.edu/pub/gnu/COPYING-2.0 or by writing to
%   Free Software Foundation, Inc.,675 Mass Ave, Cambridge, MA 02139, USA.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fid=fopen(file,'r','b');
if fid < 0
   error(sprintf('Cannot read from file %s',file));
end
nf=fread(fid,1,'long');
fp=fread(fid,1,'long')*1.E-7;
by=fread(fid,1,'short');
tc=fread(fid,1,'short');
hb=floor(tc*pow2(-14:-6));
hd=hb(9:-1:2)-2*hb(8:-1:1);
dt=tc-64*hb(9);

% hd(7)=1 CRC check
% hd(5)=1 compressed data

if ( dt == 0 ),
  d=fread(fid,Inf,'short');
else
  if (hd(5) == 1)
    % compressed data - first read scales
    nf = nf - 4;
    ncol = by / 2;
    scales = fread(fid, ncol, 'float');
    biases = fread(fid, ncol, 'float');
    d = fread(fid,[ncol, nf], 'short');
    d = repmat(1./scales,1,nf).*(d+repmat(biases,1,nf));
    d = d.';
  else
    d=fread(fid,[by/4,nf],'float').';
  end
end;
fclose(fid);
if nargout > 2
   hd(7)=0;
   hd(5)=0;
   ns=sum(hd);
   kinds=['WAVEFORM  ';'LPC       ';'LPREFC    ';'LPCEPSTRA ';'LPDELCEP  ';'IREFC     ';'MFCC      ';'FBANK     ';'MELSPEC   ';'USER      ';'DISCRETE  ';'PLP       ';'???       '];
   kind=kinds(min(dt+1,size(kinds,1)),:);
   cc='ENDACZK0';
   t=[kind(1:min(find(kind==' '))-1) reshape(['_'*ones(1,ns);cc(hd>0)],1,2*ns)];
end

author: stephencwelch

repo: Perceptual-Coding-In-Python


text:
% Copyright (C) 1993-2017, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

%%begin

% Our robot will operate within a grid world that contains free
% space where it can drive and obstacles.  We load a map of the world
load map1
% which loads a variable called map
about map
% and the cells contain zero if it is free space (driveable) and one if it is an
% obstacle.

% Now we create an instance of a robot with the PRM navigation algorithm
prm = PRM(map);

% and because PRM is a probabilistic method we will reset the random number
% generator to a known state
randinit

% Now we define the goal and start coordinates
goal = [50,30];
start = [20, 10];

% then ask the robot to plan a path to goal (it will take few seconds)
prm.plan();

% The roadmap planner does not need to know, yet, the goal or start positions, it
% chooses random points in the world and tries to find obstacle-free paths between
% them (like railway lines or freeways)

% Now we can display the obstacles and the cost to reach the goal from every
% point in the world
prm.plot();

% Now we can find a path from start to goal using the roadmap
p = prm.query(start, goal);

% and now display it
prm.plot(p);

author: petercorke

repo: robotics-toolbox-matlab


text:
% Surround matrix with NaNs
function [R]=NanMat(Data,GridSize,varargin)
%Input = Data: a m*n matrix
% GridSize: Optional, only needed if data is entered as a list
%           a 1x2 vector of size of original grid. (npx,npy)
%output = R: a m+2*n+2 matrix in which Data is preserved and surrounded by
%NaN's

%Reshape the data list to a grid
if (length(varargin) == 2)
L = reshape (Data,GridSize(1),GridSize(2));
else
    L = Data;
end
% create a matrix sames size as L but with 2 extra rows and columns, and
% fill it with NaN
R=zeros((size(L,1)+2),(size(L,2)+2))*NaN;
%Insert L into middle of matrix
R(2:end-1,2:end-1)= L;
end


text:
function spy( prob, reduce )
if nargin < 2 || ~reduce,
    A = extract( prob );
else
    A = eliminate( prob );
end
spy( A' );

% Copyright 2005-2014 CVX Research, Inc.
% See the file LICENSE.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

author: yu-jiang

repo: radpbook

save_path: github-repos/MATLAB/yu-jiang-radpbook


text:
function showTrendCorrelation(fitVals, correlations, titleString)

figure('Name', titleString, 'Color', [1, 1, 1]);
hold on
line([-1, 1], [0, 0], 'LineWidth', 2', 'Color', [0.9, 0.9, 0.9])
plot(correlations, fitVals(:, 1), 'k', 'Marker', '+', 'MarkerSize', 12', ...
    'LineWidth', 2, 'LineStyle', 'None')
hold off
set(gca, 'XLim', [-1, 1], 'XLimMode', 'manual')
xlabel('Channel-trend correlation')
ylabel('Slope of trend')
title(titleString, 'Interpreter', 'none')

figure('Name', titleString, 'Color', [1, 1, 1]);
hold on
line( [0, 0], [-1, 1], 'LineWidth', 2', 'Color', [0.9, 0.9, 0.9])
plot(fitVals(:, 1), correlations, 'k', 'Marker', '+', 'MarkerSize', 12', ...
    'LineWidth', 2, 'LineStyle', 'None')
hold off
set(gca, 'YLim', [-1, 1], 'YLimMode', 'manual')
ylabel('Channel-trend correlation')
xlabel('Slope of trend')
title(titleString, 'Interpreter', 'none')

author: VisLab

repo: EEG-Clean-Tools


text:
function [tt] = move_tt_block(tt, spos, epos, eps)
%Performs a bubble movement of a block inside a train
%   [TT] = MOVE_TT_BLOCK(TT, SPOS, EPOS, EPOS) Performs the bubble movement
%   of the SPOS-th block to the position EPOS, the intermediate blocks
%   shift to the vacant places, i.e. by (-1) if spos<eps, and +1,
%   otherwise. Requires truncation with the L2-norm accuracy EPSs
%
%
% TT-Toolbox 2.2, 2009-2012
%
%This is TT Toolbox, written by Ivan Oseledets et. al.
%Institute of Numerical Mathematics, Moscow, Russia
%webpage: http://spring.inm.ras.ru/osel
%
%For all questions, bugs and suggestions please mail
%ivan.oseledets@gmail.com
%---------------------------

if (spos==epos)
    return;
end;

d = tt.d;
n = tt.n;
r = tt.r;

% QR to spos
for i=1:spos-1
    cr = reshape(tt{i}, r(i)*n(i), r(i+1));
    [cr, rv]=qr(cr, 0);
    cr2 = reshape(tt{i+1}, r(i+1), n(i+1)*r(i+2));
    cr2 = rv*cr2;
    r(i+1) = size(cr, 2);
    tt{i} = reshape(cr, r(i), n(i), r(i+1));
    tt{i+1} = reshape(cr2, r(i+1), n(i+1), r(i+2));
end;
for i=d:-1:spos+1
    cr = reshape(tt{i}, r(i), n(i)*r(i+1));
    [cr, rv]=qr(cr.', 0);
    cr2 = reshape(tt{i-1}, r(i-1)*n(i-1), r(i));
    cr2 = cr2*(rv.');
    r(i) = size(cr, 2);
    tt{i} = reshape(cr.', r(i), n(i), r(i+1));
    tt{i-1} = reshape(cr2, r(i-1), n(i-1), r(i));
end;

% Now, start permutation
if (spos<epos) % From left to right
    for i=spos:epos-1
        cr1 = reshape(tt{i}, r(i)*n(i), r(i+1));
        cr2 = reshape(tt{i+1}, r(i+1), n(i+1)*r(i+2));
        cr = cr1*cr2;
        cr = reshape(cr, r(i), n(i), n(i+1), r(i+2));
        cr = permute(cr, [1, 3, 2, 4]);
        tempvar = n(i); n(i) = n(i+1); n(i+1) = tempvar;
        cr = reshape(cr, r(i)*n(i), n(i+1)*r(i+2));
        [u,s,v]=svd(cr, 'econ');
        s = diag(s);
        nrm = norm(s);
        r(i+1) = my_chop2(s, eps*nrm/sqrt(d-1));
        tt{i} = reshape(u(:,1:r(i+1)), r(i), n(i), r(i+1));
        v = diag(s(1:r(i+1)))*(v(:,1:r(i+1))');
        tt{i+1} = reshape(v, r(i+1), n(i+1), r(i+2));
    end;
end;

if (epos<spos) % From right  to left
    for i=spos:-1:epos+1
        cr1 = reshape(tt{i}, r(i), n(i)*r(i+1));
        cr2 = reshape(tt{i-1}, r(i-1)*n(i-1), r(i));
        cr = cr2*cr1;
        cr = reshape(cr, r(i-1), n(i-1), n(i), r(i+1));
        cr = permute(cr, [1, 3, 2, 4]);
        tempvar = n(i); n(i) = n(i-1); n(i-1) = tempvar;
        cr = reshape(cr, r(i-1)*n(i-1), n(i)*r(i+1));
        [u,s,v]=svd(cr, 'econ');
        s = diag(s);
        nrm = norm(s);
        r(i) = my_chop2(s, eps*nrm/sqrt(d-1));
        tt{i} = reshape(v(:,1:r(i))', r(i), n(i), r(i+1));
        u = u(:,1:r(i))*diag(s(1:r(i)));
        tt{i-1} = reshape(u, r(i-1), n(i-1), r(i));
    end;
end;

end

author: oseledets

repo: TT-Toolbox


text:
function test_suite=SimTest_qmt_bssfp_test
try % assignment of 'localfunctions' is necessary in Matlab >= 2016
test_functions=localfunctions();
catch % no problem; early Matlab versions can use initTestSuite fine
end
initTestSuite;

function TestSetup
setenv('ISDISPLAY','0') % go faster! Fit only 2 voxels in FitData.m
setenv('ISCITEST','1')

function test_Sim
disp('===========================================')
disp('Running simulation test for qmt_bssfp');
disp('testing Simulation Single Voxel Curve...');

Model = str2func('qmt_bssfp'); Model = Model();
savedModel_fname = fullfile(fileparts(which('qMRLab')),'Test','MoxUnitCompatible','static_savedModelsforRetrocompatibility',['qmt_bssfp.qmrlab.mat']);
if ~exist(savedModel_fname,'file')
Model.saveObj(savedModel_fname);
else
Model = Model.loadObj(savedModel_fname);
end

disp(class(Model))
try Opt = button2opts(Model.Sim_Single_Voxel_Curve_buttons,1); end
try st = Model.st; catch, try st = mean([Model.lb(:),Model.ub(:)],2); catch, st = ones(length(Model.xnames),1); end; end
if exist('Opt','var') && length(Opt)>1
[Opt(:).SNR] = deal(1000);
else
Opt.SNR=1000;
end
for iopt=1:length(Opt) % Test all simulation options
disp(['Testing ' class(Model) ' simulation option:'])
disp(Opt(iopt))
FitResults = Model.Sim_Single_Voxel_Curve(st,Opt(iopt));
% Compare inputs and outputs
fnm=fieldnames(FitResults);
FitResults = rmfield(FitResults,fnm(~ismember(fnm,Model.xnames))); fnm=fieldnames(FitResults);
[~,FitResults,GroundTruth]=comp_struct(FitResults,mat2struct(st,Model.xnames),[],[],.30);
assertTrue(isempty(FitResults) & isempty(GroundTruth),evalc('FitResults, GroundTruth'))
end
disp ..ok

function TestTeardown
setenv('ISDISPLAY','') % go faster! Fit only 2 voxels in FitData.m
setenv('ISCITEST','')

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function net = cnn_mnist_init(varargin)
% CNN_MNIST_LENET Initialize a CNN similar for MNIST
opts.batchNormalization = true ;
opts.networkType = 'simplenn' ;
opts = vl_argparse(opts, varargin) ;

rng('default');
rng(0) ;

f=1/100 ;
net.layers = {} ;
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(5,5,1,20, 'single'), zeros(1, 20, 'single')}}, ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(5,5,20,50, 'single'),zeros(1,50,'single')}}, ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(4,4,50,500, 'single'),  zeros(1,500,'single')}}, ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(1,1,500,10, 'single'), zeros(1,10,'single')}}, ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'softmaxloss') ;

% optionally switch to batch normalization
if opts.batchNormalization
  net = insertBnorm(net, 1) ;
  net = insertBnorm(net, 4) ;
  net = insertBnorm(net, 7) ;
end

% Meta parameters
net.meta.inputSize = [27 27 1] ;
net.meta.trainOpts.learningRate = 0.001 ;
net.meta.trainOpts.numEpochs = 20 ;
net.meta.trainOpts.batchSize = 100 ;

% Fill in defaul values
net = vl_simplenn_tidy(net) ;

% Switch to DagNN if requested
switch lower(opts.networkType)
  case 'simplenn'
    % done
  case 'dagnn'
    net = dagnn.DagNN.fromSimpleNN(net, 'canonicalNames', true) ;
    net.addLayer('error', dagnn.Loss('loss', 'classerror'), ...
             {'prediction','label'}, 'error') ;
  otherwise
    assert(false) ;
end

% --------------------------------------------------------------------
function net = insertBnorm(net, l)
% --------------------------------------------------------------------
assert(isfield(net.layers{l}, 'weights'));
ndim = size(net.layers{l}.weights{1}, 4);
layer = struct('type', 'bnorm', ...
               'weights', {{ones(ndim, 1, 'single'), zeros(ndim, 1, 'single')}}, ...
               'learningRate', [1 1 0.05], ...
               'weightDecay', [0 0]) ;
net.layers{l}.biases = [] ;
net.layers = horzcat(net.layers(1:l), layer, net.layers(l+1:end)) ;

author: willard-yuan

repo: cnn-for-image-retrieval


text:
% In this version of this function we try to return a (NxNxN) cube, while still
% keeping the aspect ratio

function [voxel_model] = pointcloud_to_voxels_3(model,N,bbox,voxelization_type,binaryMax)

if(isempty(model.x))
    warning('empty model passed to this function!');
    voxel_model = zeros([N N N]);

    return;
end

if(~exist('bbox','var'))
    bbox = [];
end

if(~exist('voxelization_type','var'))
    voxelization_type = 'binary';
end

if(~exist('binaryMax','var'))
    binaryMax = 1; %This is only for backward compatibilty. Otherwise it should be 255 to be compatible with 'intensity' type, in addition to other benefits!
end

[X,Y,Z] = model_to_components(model);

%--- Obtain the boundaries of the object, taking into account the bbox if
%needed
if(isempty(bbox))
    mX = min(X);
    MX = max(X);
    mY = min(Y);
    MY = max(Y);
    mZ = min(Z);
    MZ = max(Z);
else
    %-- in case of a bounding box, we only need to check the 4 corners
    x1 = bbox(1);
    y1 = bbox(2);
    z1 = bbox(3);
    w = bbox(4);
    h = bbox(5);
    d = bbox(6);
    x2 = x1+w;
    y2 = y1+h;
    z2 = z1+d;
    phi = bbox(7);

    cornersX = [x1 x1 x2 x2];
    cornersY = [y1 y2 y1 y2];
    xc = mean([x1,x2]);
    yc = mean([y1,y2]);

    % rotate the coordinates around the center
    cornersX = cornersX-repmat(xc,1,size(cornersX,2));
    cornersY = cornersY-repmat(yc,1,size(cornersX,2));
    rotated = [cosd(phi) -sind(phi); sind(phi) cosd(phi)]*[cornersX;cornersY];
    cornersX = rotated(1,:)+repmat(xc,1,size(cornersX,2));
    cornersY = rotated(2,:)+repmat(yc,1,size(cornersY,2));

    %-- now pick the limits of the corner coordinates
    mX = min(cornersX);
    MX = max(cornersX);
    mY = min(cornersY);
    MY = max(cornersY);
    mZ = min([z1,z2]); %z is interpreted differently, as the rotation is only around z
    MZ = max([z1,z2]);

end

maxmax = max([MX-mX,MY-mY,MZ-mZ]);
if(maxmax == 0) maxmax = 1; end %This happens only when there's a single point in the scene (the empty scene case has been taken care of before);

N1 = round(N*(MX-mX)/maxmax);
N2 = round(N*(MY-mY)/maxmax);
N3 = round(N*(MZ-mZ)/maxmax);

% This happens when the input point cloud is in in fact less than 3D
if(~N1), N1 = 1;end;
if(~N2), N2 = 1;end;
if(~N3), N3 = 1;end;

voxel_model = zeros(N1,N2,N3);

ax = 1/(MX-mX)*(N1-1);
bx = -mX;
ay = 1/(MY-mY)*(N2-1);
by = -mY;
az = 1/(MZ-mZ)*(N3-1);
bz = -mZ;

% This happens when the input point cloud is in fact less than 3D
if isnan(ax), ax = 0; end;
if isnan(ay), ay = 0; end;
if isnan(az), az = 0; end;

X2 = round(ax*(X+bx))+1;
Y2 = round(ay*(Y+by))+1;
Z2 = round(az*(Z+bz))+1;

switch(voxelization_type )
    case 'binary'
        voxel_model(sub2ind(size(voxel_model),X2,Y2,Z2)) = binaryMax ;
    case 'density'
        [C,ia,ic] = unique([X2,Y2,Z2],'rows');
        h = hist(ic,1:max(ic));
        ind = sub2ind(size(voxel_model),C(:,1),C(:,2),C(:,3));
        voxel_model(ind) = h;

        %max normalization
        voxel_model = voxel_model ./ max(voxel_model(:)) * 255;
    otherwise
        error('voxelization type not supported');
end

% Add the required symmetric zero padding to make it a cube
Nzx = floor( (N - N1)/2 );
Nzy = floor( (N - N2)/2 );
Nzz = floor( (N - N3)/2 );
final = zeros(N,N,N);
final(Nzx+1:Nzx+N1 , Nzy+1:Nzy+N2 ,Nzz+1:Nzz+N3) = voxel_model;
voxel_model = final;