
    % Use the BSDS images in the train directory.
    % Currently only that directory contains the texture analized files as
    % well as the superpixel segmentation files.
    train_or_test = 'train';
    BSDSdir = [bsdsRoot(),'/images/',train_or_test,'/'];
    files = dir([BSDSdir,'*.jpg']);
    numFiles = length(files);
    features_pos = [];
    features_neg = [];
    
    % Loop over all the BSDS images
    for iFile = 1:numFiles
        disp(['Processing file ',num2str(iFile),' out of ',num2str(numFiles)]);
        
        % Read the image
        imgFilename = [bsdsRoot(),'/images/',train_or_test,'/', ...
                    files(iFile).name];
        img = im2double(imread(imgFilename));
        
        % Read the file with the texture analysis
        matlabSaveFile = [BSDSdir,files(iFile).name(1:end-3),'mat'];
        temp = load(matlabSaveFile);
        textureMap = temp.tmap;
        
        segFilename = [bsdsRoot(),'/images/',train_or_test,'/',...
                       files(iFile).name(1:end-4),'_1.seg'];
        iidBSDS = str2double(files(iFile).name(1:end-4));
        
        % Read the ground truth segmentations
        [segs,uids] = readSegs('gray',iidBSDS);
        
        % Read the superpixel segmentation image computed by NCuts
        segimage = readSeg(segFilename);
        
        % Compute the adjacency matrix
        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 superpixel edges
        for sup1=1:nseg
            for sup2=1:sup1-1
                % Check if two superpixels are adjacent
                if (adjmat(sup1,sup2))
                    % Compute a histogram of textons for each superpixel
                    histogram1 = hist(textureMap(segimage(:)==sup1),1:32);
                    histogram1 = histogram1 / sum(histogram1);
                    histogram2 = hist(textureMap(segimage(:)==sup2),1:32);
                    histogram2 = histogram2 / sum(histogram2);
                    sup1segment = mode(segs{1}(segimage(:) == sup1));
                    sup2segment = mode(segs{1}(segimage(:) == sup2));
                    
                    % Store the features in different arrays based on
                    % whether the superpixels belong to the same ground
                    % truth segment or not
                    if (sup1segment == sup2segment)
                        features_pos = [features_pos;[histogram1,histogram2]];
                    else
                        features_neg = [features_neg;[histogram1,histogram2]];
                    end
                end
            end
        end
    end
    
    % Compute the ground truth labels (0: different segment, 1: same
    % segment)
    labels = zeros(size(features_pos,1) + size(features_neg,1),1);
    labels(1:size(features_pos,1)) = 1;
    
    % Compute the superpixel edge features
    % Currently the covariance between the two texton histograms is
    % computed. To make it invariant to the superpixel order (i,j) or (j,i)
    % the two corresponding covariance matrices are added and the lower
    % diagonal of the result is taken as the feature vector.
    numFeatures = 33*34 / 2;
    features = zeros(numel(labels),numFeatures);
    all_features = [features_pos;features_neg];
    lowDiagFeaturesIdx = tril(ones(33,33)) > 0;
    for i=1:size(all_features,1)
        quadF1 = [all_features(i,1:32)';1] * [all_features(i,33:end),1];
        quadF2 = [all_features(i,33:end)';1] * [all_features(i,1:32),1];
        quadF = quadF1 + quadF2;
        features(i,:) = quadF(lowDiagFeaturesIdx);
    end
    
    % Train the logistic classifier on the data
    % beta are the learned logistic parameters such that:
    % P(superpixels are in the same segment) = 1 / (1 + exp(-x*beta)) where
    % x is a feature vector for these superpixels
    [beta,p,lli] = logist2(labels,features);
