function ps3p1
% Edit this file. 
% Modify the two functions inside marked with "IMPLEMENT THIS"
% NAME: write your name here

% read in the image
imagename = 'Cactus.tiff';  im =imread(imagename);

% make the relaxation labeling structure
rls=makeRLStruct(im, 2);

% initialize the network: this is a function that you need to implement,
% see below to find its definition
rls=init_p(rls, im);

% display the segmentaion in fig. 1 using the initial values
display_labels(rls, 1);

% update the probabilities using the rules from the handout
% you need to implement the function compatibility() in order to do this
% you may want to call update more than once
rls = update(rls);

% display the segmenetation in fig. 2 after running the relaxation labeling
% network
display_labels(rls, 2);

 
%% IMPLEMENT THIS: compatibility()
% This corresponds to r_ij(lambda, mu) from the problem set. Node i is
% represented as (xi,yi) and node j=(xj,yj)
    function r=compatibility(xi, yi,xj, yj, lambda, mu)
        % implement your own rule(s)
        r=0; % set r to the value that your principle suggests instead of 0
    end

%% IMPLEMENT THIS: init_p()
% You need to change the two lines below 
    function rls=init_p(rls, im)
        mxim = max(max(im));
        for x=1:size(rls.p,1)
            for y=1:size(rls.p,2)
                % initial confidence in label 1 (object)
                rls.p(x,y,1)=0.50;                % IMPLEMENT YOUR OWN
                % initial confidence in label 2 (background)
                rls.p(x,y,2)=0.50;                % IMPLEMENT YOUR OWN

                sm = sum(rls.p(x,y,:));
                rls.p(x,y,:) = rls.p(x,y,:) / sm;
            end
        end
    end

%% display_labels()
    function display_labels(rls, figid)
        iml = zeros(size(rls.p,1), size(rls.p,2));
        
        for x=1:size(rls.p,1)
           for y=1:size(rls.p,2)
               iml(x,y) = ( rls.p(x,y,1) == max(rls.p(x,y,:)));
           end
        end
        
        figure(figid)
        imshow(iml)
    end
%% makeRLStruct()
    function rls=makeRLStruct(im, nlabels)
       rls.p = ones(size(im,1), size(im,2), nlabels); 
    end

%% compute_support()
    function s=compute_support(rls, xi, yi, lambda, nbx, nby)
        s=0;
        %nbx=2; nby=2;
        if xi < nbx || xi > size(rls.p,1)-nbx ...
                || yi < nby|| yi > size(rls.p,2)
           return 
        end
        
        for nx=xi-nbx+1:xi+nbx-1
            for ny=yi-nby+1:yi+nby-1
                for mu=1:size(rls.p,3)
                    r=compatibility(xi, yi,nx, ny, lambda, mu);

                    p = rls.p(nx,ny,mu);

                    s = s+r*p;
                end

            end
        end
    end

%% update()
    function rls = update(rls)
        nbx=2; nby=2;
        newp = rls.p*0;
        for x=nbx:size(rls.p,1)-nbx
            for y=nby:size(rls.p,2)-nby
                for lambda=1:size(rls.p,3)
                    den = 0;
                    for mu=1:size(rls.p,3)
                        simu=compute_support(rls, x,y, mu, nbx, nby);
                        den = den + rls.p(x,y,mu)*simu;
                    end
                    silambda=compute_support(rls, x,y, lambda, nbx, nby);
                    newp(x,y,lambda) = rls.p(x,y,lambda)*silambda / den;
                end
            end
        end
        
        rls.p = newp;
    end
end