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
140
141
142
143
144
145
146
147
148
149
150
|
seed = 0;
rand('state', seed);
randn('state', seed);
nrows = 5;
ncols = 5;
npixels = nrows*ncols;
% we number pixels in transposed raster scan order (top to bottom, left to right)
% H(i,j) is the number of the hidden node at (i,j)
H = reshape(1:npixels, nrows, ncols);
% O(i,j) is the number of the obsevred node at (i,j)
O = reshape(1:npixels, nrows, ncols) + length(H(:));
% Make a Bayes net where each hidden pixel generates an observed pixel
% but there are no connections between the hidden pixels.
% We use this just to generate noisy versions of known images.
N = 2*npixels;
dag = zeros(N);
for i=1:nrows
for j=1:ncols
dag(H(i,j), O(i,j)) = 1;
end
end
K = 2; % number of discrete values for the hidden vars
ns = ones(N,1);
ns(H(:)) = K;
ns(O(:)) = 1;
% make image with vertical stripes
I = zeros(nrows, ncols);
for j=1:2:ncols
I(:,j) = 1;
end
% each "hidden" node will be instantiated to the pixel in the known image
% each observed node has conditional Gaussian distribution
eclass = ones(1,N);
%eclass(H(:)) = 1;
%eclass(O(:)) = 2;
eclass(H(:)) = 1:npixels;
eclass(O(:)) = npixels+1;
bnet = mk_bnet(dag, ns, 'discrete', H(:), 'equiv_class', eclass);
%bnet.CPD{1} = tabular_CPD(bnet, H(1), 'CPT', normalise(ones(1,K)));
for i=1:nrows
for j=1:ncols
bnet.CPD{H(i,j)} = root_CPD(bnet, H(i,j), I(i,j) + 1);
end
end
% If H(i,j)=1, O(i,j)=+1 plus noise
% If H(i,j)=2, O(i,j)=-1 plus noise
sigma = 0.5;
bnet.CPD{eclass(O(1,1))} = gaussian_CPD(bnet, O(1,1), 'mean', [1 -1], 'cov', reshape(sigma*ones(1,K), [1 1 K]));
ofactor = bnet.CPD{eclass(O(1,1))};
%ofactor = gaussian_CPD('self', 2, 'dps', 1, 'cps', [], 'sz', [K O], 'mean', [1 -1], 'cov', reshape(sigma*ones(1,K), [1 1 K)));
data = sample_bnet(bnet);
img = reshape(data(O(:)), nrows, ncols)
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now create MRF represented as a factor graph to try and recover the scene
% VEF(i,j) is the number of the factor for the vertical edge between HV(i,j) - HV(i+1,j)
VEF = reshape((1:(nrows-1)*ncols), nrows-1, ncols);
% HEF(i,j) is the number of the factor for the horizontal edge between HV(i,j) - HV(i,j+1)
HEF = reshape((1:nrows*(ncols-1)), nrows, ncols-1) + length(VEF(:));
nvars = npixels;
nfac = length(VEF(:)) + length(HEF(:));
G = zeros(nvars, nfac);
N = length(ns);
eclass = zeros(1, nfac); % eclass(i)=j means factor i gets its params from factors{j}
vfactor_ndx = 1; % all vertcial edges get their params from factors{1}
hfactor_ndx = 2; % all vertcial edges get their params from factors{2}
for i=1:nrows
for j=1:ncols
if i < nrows
G(H(i:i+1,j), VEF(i,j)) = 1;
eclass(VEF(i,j)) = vfactor_ndx;
end
if j < ncols
G(H(i,j:j+1), HEF(i,j)) = 1;
eclass(HEF(i,j)) = hfactor_ndx;
end
end
end
% "kitten raised in cage" prior - more likely to see continguous vertical lines
vfactor = tabular_kernel([K K], softeye(K, 0.9));
hfactor = tabular_kernel([K K], softeye(K, 0.5));
factors = cell(1,2);
factors{vfactor_ndx} = vfactor;
factors{hfactor_ndx} = hfactor;
ev_eclass = ones(1,N); % every observation factor gets is params from ofactor
ns = K*ones(1,nvars);
%fg = mk_fgraph_given_ev(G, ns, factors, {ofactor}, num2cell(img), 'equiv_class', eclass, 'ev_equiv_class', ev_eclass);
fg = mk_fgraph_given_ev(G, ns, factors, {ofactor}, img, 'equiv_class', eclass, 'ev_equiv_class', ev_eclass);
bnet2 = fgraph_to_bnet(fg);
% inference
maximize = 1;
engine = {};
engine{1} = belprop_fg_inf_engine(fg, 'max_iter', npixels*2);
engine{2} = jtree_inf_engine(bnet2);
nengines = length(engine);
% on fg, we have already included the evidence
evidence = cell(1,npixels);
tic; [engine{1}, ll(1)] = enter_evidence(engine{1}, evidence, 'maximize', maximize); toc
% on bnet2, we must add evidence to the dummy nodes
V = fg.nvars;
dummy = V+1:V+fg.nfactors;
N = max(dummy);
evidence = cell(1, N);
evidence(dummy) = {1};
tic; [engine{2}, ll(2)] = enter_evidence(engine{2}, evidence); toc
Ihat = zeros(nrows, ncols, nengines);
for e=1:nengines
for i=1:nrows
for j=1:ncols
m = marginal_nodes(engine{e}, H(i,j));
Ihat(i,j,e) = argmax(m.T)-1;
end
end
end
Ihat
|
