/
spm_mb_merge.m
140 lines (127 loc) · 3.97 KB
/
spm_mb_merge.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
function out = spm_mb_merge(cfg)
% Combine tissue maps together
%__________________________________________________________________________
% Copyright (C) 2020 Wellcome Centre for Human Neuroimaging
% $Id: spm_mb_merge.m 7892 2020-07-10 16:39:18Z john $
out = struct('mu','priors');
odir = cfg.odir{1};
onam = cfg.onam;
res_file = cfg.result{1};
ix = cfg.ix;
bb = cfg.bb;
[~,~,ext] = fileparts(res_file);
if strcmp(ext,'.mat')
r = load(res_file);
if ~isfield(r,'sett'), error('Incorrect format.'); end
if isfield(r.sett.mu,'create')
mu_name = r.sett.mu.create.mu;
else
mu_name = r.sett.mu.exist.mu;
end
if ~exist(mu_name,'file')
[~,nam,ext] = fileparts(mu_name);
[pth,~,~] = fileparts(res_file);
if ~exist(fullfile(pth,[nam ext]),'file')
error('Cannot find file %s or %s.',mu_name, fullfile(pth,[nam ext]));
else
mu_name = fullfile(pth,[nam ext]);
end
end
% Save the reordered priors for each population
out.priors = cell(numel(r.sett.gmm),1);
for p=1:numel(r.sett.gmm)
matname = fullfile(odir,sprintf('prior_%s_%d.mat',onam,p));
if max(r.sett.gmm(p).mg_ix) ~= numel(ix)
error('The indices are incorrectly specified.');
end
mg_ix = ix(r.sett.gmm(p).mg_ix);
[mg_ix,si] = sort(mg_ix);
pr = r.sett.gmm(p).pr;
pr{1} = pr{1}(:,si); % m
pr{2} = pr{2}(:,si); % b
pr{3} = pr{3}(:,:,si); % V
pr{4} = pr{4}(:,si); % n
save(matname,'mg_ix','pr');
out.priors{p} = matname;
end
else
mu_name = res_file;
end
Nmu = nifti(mu_name);
[ind,mat] = SubVol(Nmu,bb);
mu = single(Nmu.dat(ind{:},:));
mu1 = merge_mu(mu,ix);
Nmu.dat.fname = fullfile(odir,['mu_' onam '.nii']);
Nmu.dat.dim = size(mu1);
Nmu.mat = mat;
Nmu.mat0 = mat;
create(Nmu);
Nmu.dat(:,:,:,:) = mu1;
out.mu = {Nmu.dat.fname};
%==========================================================================
%==========================================================================
function mu1 = merge_mu(mu,ix)
K = size(mu,4)+1;
if numel(ix)~=K
error('Incorrect index dimension.\n The mean image suggests K=%d, whereas the indices assume K=%d.',K, numel(ix));
end
if any(ix>K) || any(ix<1)
error('Index out of bounds.');
end
K1 = max(ix);
if numel(unique(ix))~=K1, error('Something wrong with indices.'); end
d = [size(mu) 1 1];
d = d(1:3);
mu1 = zeros([d K1-1],'single');
mK = combine(mu,find(ix==K1));
for k1=1:(K1-1)
mk = combine(mu,find(ix==k1));
mu1(:,:,:,k1) = mk - mK;
end
%==========================================================================
%==========================================================================
function mk = combine(mu,ind)
K = size(mu,4)+1;
d = [size(mu,1) size(mu,2) size(mu,3)];
if numel(ind)==1
if ind<K
mk = mu(:,:,:,ind);
else
mk = zeros(d,'single');
end
else
mx = 0;
for k=ind
if k<K
mx = max(mu(:,:,:,k),mx);
end
end
mk = 0;
for k=ind
if k<K
mk = mk + exp(mu(:,:,:,k)-mx);
else
mk = mk + exp(-mx);
end
end
mk = log(mk)+mx;
end
%==========================================================================
%==========================================================================
function [ind,mat] = SubVol(Nii,bb)
% Information for extracting a subvolume
% FORMAT [ind,mat] = SubVol(Nii,bb)
% V - mapped file
% bb - bounding box
dm0 = [size(Nii.dat) 1];
dm0 = dm0(1:3);
msk = ~isfinite(bb(1,:));
bb(1,msk) = 1;
bb(1,:) = round(max(bb(1,:),1));
msk = ~isfinite(bb(2,:));
bb(2,msk) = dm0(msk);
bb(2,:) = round(min(bb(2,:),dm0));
bb = sort(bb);
ind = {bb(1,1):bb(2,1),bb(1,2):bb(2,2),bb(1,3):bb(2,3)};
mat = Nii.mat*[eye(4,3) [bb(1,:)'-1; 1]];
%==========================================================================