% Once the parameters of the logistic classfier have been learned, this
% function uses those parameters to assign an edge weight to each
% superpixel edge. This function shows the resulting edges on an image for
% convinience.
% Input: 
%   beta - 1xD matrix of learned parameters. Currently it is 1x561
%   iidBSDS - an ID of an image from the Berkeley dataset
%   matlabFilesDir - an optional parameter specifying the directory of the
%   .mat files storing the texture analysis. The default is the BSDS
%   directory
%
% Returns:
%   features - ExD matrix of features for each superpixel edge
%   edgeIdx  - a 2D (NxN) matrix mapping each pair of superpixels (i,j) to their
%   edge index k
%   edgeStrength - NxN matrix of edge strengths inferred with logistic
%   classifier
%   segimage - the superpixel segmentation image.
%   img - the image with the specified ID in the BSDS.

function [features, edgeIdx, edgeStrength, segimage, img] = testSuperpixelSegmentation(beta,iidBSDS, matFilesDir)

    % Read the image
    imFilename = imgFilename(iidBSDS);
    img = im2double(imread(imFilename));
    
    % Read the file with the corresponding texture segmentation
    % There is an option to provide a separate directory for the .mat files
    % in the 3rd parameter to the function
    if (nargin < 3)
        matlabSaveFile = [imFilename(1:end-3),'mat'];
    else
        matlabSaveFile = [matFilesDir,num2str(iidBSDS),'.mat'];
    end
    
    temp = load(matlabSaveFile);
    textureMap = temp.tmap;
    
    % Read the superpixel segmentation
    segFilename = [imFilename(1:end-4),'_1.seg'];
    segimage = readSeg(segFilename);
    
    % Compute the adjacency matrix for superpixels
    nseg = max(segimage(:));
    adjmat = zeros(nseg, nseg);

    [height, width] = size(segimage);

    for y = 1:height-1
        for x = 1:width-1
            s1 = segimage(y, x);
            s2 = segimage(y+1, x);
            s3 = segimage(y, x+1);
            adjmat(s1, s2) = 1;
            adjmat(s2, s1) = 1;
            adjmat(s1, s3) = 1;
            adjmat(s3, s1) = 1;
        end
    end

    % Compute the features for superpixels
    lowDiagFeaturesIdx = tril(ones(33,33)) > 0;
    adjmat(sub2ind(size(adjmat),1:nseg,1:nseg)) = 0;
    
    features = zeros(sum(adjmat(:))/2, sum(lowDiagFeaturesIdx(:)));
    edgeStrength = zeros(nseg,nseg);
    edgeIdx = zeros(nseg,nseg);
    k=1;
    for sup1=1:nseg
        for sup2=1:sup1-1
            % Check if two superpixels are adjacent
            if (adjmat(sup1,sup2))
                
                % Compute the features for the superpixels
                histogram1 = hist(textureMap(segimage(:)==sup1),1:32);
                histogram1 = histogram1 / sum(histogram1);
                histogram2 = hist(textureMap(segimage(:)==sup2),1:32);
                histogram2 = histogram2 / sum(histogram2);
                quadF1 = [histogram1';1] * [histogram2,1];
                quadF2 = [histogram2';1] * [histogram1,1];
                quadF = quadF1 + quadF2;
                features(k,:) = quadF(lowDiagFeaturesIdx);
                
                % Compute the current edge index
                edgeIdx(sup1,sup2) = k;
                edgeIdx(sup2,sup1) = k;
                
                % Infer the edge strength with the learned logistic classifier
                edgeStrength(sup1,sup2) = 1 - (1 / (1 + exp(-features(k,:)*beta)));
                edgeStrength(sup2,sup1) = edgeStrength(sup1,sup2);
                k=k+1;
            end
        end
    end
    
    % Create an edge image where superpixel edges are colored according to
    % their strength
    edgeImg = zeros(size(segimage));
    for y = 1:height-1
        for x = 1:width-1
            s1 = segimage(y, x);
            s2 = segimage(y+1, x);
            s3 = segimage(y, x+1);
            if (adjmat(s1,s2))
                edgeImg(y,x) = edgeStrength(s1,s2);
            elseif (adjmat(s1,s3))
                edgeImg(y,x) = edgeStrength(s1,s3);
            end
        end
    end
    
    imagesc(edgeImg);