Daniel@0: seed = 0;
Daniel@0: rand('state', seed);
Daniel@0: randn('state', seed);
Daniel@0: 
Daniel@0: nrows = 3;
Daniel@0: ncols = 3;
Daniel@0: npixels = nrows*ncols;
Daniel@0: 
Daniel@0: % we number pixels in transposed raster scan order (top to bottom, left to right)
Daniel@0: 
Daniel@0: % hidden var
Daniel@0: HV = reshape(1:npixels, nrows, ncols);
Daniel@0: % observed var
Daniel@0: OV = reshape(1:npixels, nrows, ncols) + length(HV(:));
Daniel@0: 
Daniel@0: % observed factor
Daniel@0: OF = reshape(1:npixels, nrows, ncols);
Daniel@0: % vertical edge factor VEF(i,j) is the factor for edge HV(i,j) - HV(i+1,j)
Daniel@0: VEF = reshape((1:(nrows-1)*ncols), nrows-1, ncols) + length(OF(:));
Daniel@0: % horizontal edge factor HEF(i,j) is the factor for edge HV(i,j) - HV(i,j+1)
Daniel@0: HEF = reshape((1:nrows*(ncols-1)), nrows, ncols-1) + length(OF(:)) + length(VEF(:));
Daniel@0: 
Daniel@0: nvars = length(HV(:))+length(OV(:));
Daniel@0: assert(nvars == 2*npixels);
Daniel@0: nfac = length(OF(:)) + length(VEF(:)) + length(HEF(:));
Daniel@0: 
Daniel@0: K = 2; % number of discrete values for the hidden vars
Daniel@0: %O = 1; % each observed pixel is a scalar
Daniel@0: O = 2; % each observed pixel is binary
Daniel@0: 
Daniel@0: factors = cell(1,3);
Daniel@0: 
Daniel@0: % hidden states generate observed 0 or 1 plus noise
Daniel@0: %factors{2} = cond_gauss1_kernel(K, O, 'mean', [0 1], 'cov', [0.1 0.1]);
Daniel@0: pnoise = 0.2;
Daniel@0: factors{1} = tabular_kernel([K O], [1-pnoise pnoise; pnoise 1-pnoise]);
Daniel@0: ofactor = 1;
Daniel@0: 
Daniel@0: % encourage compatibility between neighboring vertical pixels
Daniel@0: factors{2} = tabular_kernel([K K], [0.8 0.2; 0.2 0.8]);
Daniel@0: vedge_factor = 2;
Daniel@0: 
Daniel@0: %% no constraint between neighboring horizontal pixels
Daniel@0: %factors{3} = tabular_kernel([K K], [0.5 0.5; 0.5 0.5]);
Daniel@0: 
Daniel@0: factors{3} = tabular_kernel([K K], [0.8 0.2; 0.2 0.8]);
Daniel@0: hedge_factor = 3;
Daniel@0: 
Daniel@0: 
Daniel@0: 
Daniel@0: factor_ndx = zeros(1, 3);
Daniel@0: G = zeros(nvars, nfac);
Daniel@0: ns = [K*ones(1,length(HV(:))) O*ones(1,length(OV(:)))];
Daniel@0: 
Daniel@0: N = length(ns);
Daniel@0: %cnodes = OV(:);
Daniel@0: cnodes = [];
Daniel@0: dnodes = 1:N;
Daniel@0: 
Daniel@0: for i=1:nrows
Daniel@0:   for j=1:ncols
Daniel@0:     G([HV(i,j), OV(i,j)], OF(i,j)) = 1;
Daniel@0:     factor_ndx(OF(i,j)) = ofactor;
Daniel@0: 
Daniel@0:     if i < nrows
Daniel@0:       G(HV(i:i+1,j), VEF(i,j)) = 1;
Daniel@0:       factor_ndx(VEF(i,j)) = vedge_factor;
Daniel@0:     end
Daniel@0: 
Daniel@0:     if j < ncols
Daniel@0:       G(HV(i,j:j+1), HEF(i,j)) = 1;
Daniel@0:       factor_ndx(HEF(i,j)) = hedge_factor;
Daniel@0:     end
Daniel@0: 
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: 
Daniel@0: fg = mk_fgraph(G, ns, factors, 'discrete', dnodes, 'equiv_class', factor_ndx);
Daniel@0: 
Daniel@0: if 1
Daniel@0:   % make image with vertical stripes
Daniel@0:   I = zeros(nrows, ncols);
Daniel@0:   for j=1:2:ncols
Daniel@0:     I(:,j) = 1;
Daniel@0:   end
Daniel@0: else
Daniel@0:   % make image with square in middle
Daniel@0:   I = zeros(nrows, ncols);
Daniel@0:   I(3:6,3:6) = 1;
Daniel@0: end
Daniel@0: 
Daniel@0:   
Daniel@0: % corrupt image
Daniel@0: O = mod(I + (rand(nrows,ncols)> (1-pnoise)), 2);
Daniel@0: 
Daniel@0: maximize = 1;
Daniel@0: engine = belprop_fg_inf_engine(fg, 'maximize', maximize, 'max_iter', npixels*5);
Daniel@0: 
Daniel@0: evidence = cell(1, nvars);
Daniel@0: onodes = OV(:);
Daniel@0: evidence(onodes) = num2cell(O+1); % values must be in range {1,2}
Daniel@0: 
Daniel@0: engine = enter_evidence(engine, evidence);
Daniel@0: 
Daniel@0: for i=1:nrows
Daniel@0:   for j=1:ncols
Daniel@0:     m = marginal_nodes(engine, HV(i,j));
Daniel@0:     Ihat(i,j) = argmax(m.T)-1;
Daniel@0:   end
Daniel@0: end
Daniel@0: 
Daniel@0: Ihat