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