function demoing(write)
if ~exist('write', 'var'), write = false; end;
SI =5; SX = 6; r = 1.5; sNcut = 0.02;
I = imread('3.jpg');
segI = main(I, SI, SX, r, sNcut);
% show
for i=1:length(segI)
    if ~write
        figure; imshow(segI{i});
    end
end
end

function SegI = main(image, sig_i, sig_x, r, min_cut)

[num_row, num_col, components] = size(image);
N=num_row * num_col;

% create graph (vertical list)
V=reshape(image, N, components);

% Step 1. Compute weight matrix W, and D
W = computeW(image, sig_i, sig_x, r);

% nodes are labled by pixel #
Seg=(1:N)';% the first segment has whole nodes. [1 2 3 ... N]'
[Seg]=Ncuts(Seg, W, min_cut);

% Convert node ids into images
for i=1:length(Seg)
    subV = zeros(N, components); %ones(N, c) * 255;
    subV(Seg{i}, :) = V(Seg{i}, :);
    SegI{i} = uint8(reshape(subV, num_row, num_col, components));
end
end

function W = computeW(in_seg, sig_i, sig_x, r);
% input will be the image, sig_i (in the exponential weighting), sig_x
% (in the exponential weighting), and r the threshold distance
% get the parameters of the image the number of rows, columns, and the
[num_row, num_col, components] = size(in_seg);
N = num_row * num_col;
W = zeros(N,N);

% Feature Vectors
F = in_seg;
F = reshape(F, N, 1, components); % col vector

%create a matrix of pixel w/ coordinate values
xcoord=repmat((1:num_row)', 1, num_col);
ycoord=repmat((1:num_col), num_row, 1);
X=cat(3, xcoord, ycoord);
X=reshape(X, N, 1, 2); % col vector

%create symmetrical matrix W with weights
% we have two sets of coordinates here one i for the node we are
% using as the center node, j nodes which are the ones which it is
% connected to
for i_col=1:num_col
    for i_row=1:num_row
%         % compute the j nodes' coordinates using the r value
%         % j will be within i_row+r and i_row-r same applies for columns
% use floor here because discrete
        j_col = (i_col - floor(r)) : (i_col + floor(r)); % vector
        j_row = ((i_row - floor(r)) :(i_row + floor(r)))';
        j_col = j_col(j_col >= 1 & j_col <= num_col);
        j_row = j_row(j_row >= 1 & j_row <= num_row);

        % compute the locations of the i nodes
        % index in 1D array format
        % compute the locations of the j nodes connected to i
        % index in 1D array format ??
        i = i_row + (i_col - 1) * num_row;
        j = repmat(j_row, 1, length(j_col)) + repmat((j_col -1) * num_row, length(j_row), 1);
        j = reshape(j, length(j_col) * length(j_row), 1); % a col vector

        %compute X_i-X_j and F_i-F_j
        X_j = X(j, 1, :);% store coordinates of each index j [x,y] of each connected node
        X_i = repmat(X(i, 1, :), length(j), 1);
        %compute difference
        DiffX = X_i - X_j;
        %calculate ||F(i)-F(j)||2 where 2 represents L2 the euclidian distance
        DiffX = sum(DiffX .* DiffX, 3); %distance formula in 3D

        %retain all points that are relevant using the condition that they are <r
        retained = find(sqrt(DiffX) <= r);
        j = j(retained);
        DiffX = DiffX(retained);

        % feature vector disimilarity
        F_j = F(j, 1, :);
        F_i = repmat(F(i, 1, :), length(j), 1);
        DiffF = F_i - F_j;
        DiffF = sum(DiffF .* DiffF, 3); % squared euclid distance

        %calculate similarity W matrix
        W(i, j) = exp(-DiffF / (sig_i*sig_i)) .* exp(-DiffX / (sig_x*sig_x));
    end
end
end

function [seg_out ] = Ncuts(in_seg, W, min_cut)
% seg_out is the segments returned after cutting
% ncut is the normalized cut values of each seg_out part
% seg_in represents the segment fed in to be cut
% W is the weight matrix computed
% min_cut is the minimum cut value (min Ncut(A,B))
% min_area is the minimum input area (in paper they used a 5x5 grid)

% compute NxN D matrix from section 2.1
% recall d(i)=summation_{j} w(i,j) and D is NxN with d on its diag
N=length(W);
d=sum(W, 2);
D=spdiags(d,0,N,N);

% as per section section 3 we use the eig_vec with the second smallest
% eig_val given by the definition of the rayleigh quotient
[eig_vec,diag_eigval] = eigs(D-W, D, 2, 'sm');
eig_vec2=eig_vec(:, 2);% grab second smallest

% https://www.quora.com/What-is-an-eigenvector-of-a-covariance-matrix
% "...and the second eigenvector is orthogonal (perpendicular) to the
% first. (for euclidian space)" which i'm guessing means is the center
% of the data. knowing that fminsearch will find data based on its initial
% value input we can choose the center value to predict a minimum as
% follows

% HOW TO GET THE MINIMUM EQUATION THING?????
t=mean(eig_vec2);
t=fminsearch('NcutValue', t, [], eig_vec2, W, D);

% get cuts in A and B with the condition that one is less than one is
% greatre than
A=find(eig_vec2>t);
B=find(eig_vec2<=t);

% ncut computed value is larger than the condition we set break recursion and end the
% function
ncut=NcutValue(t, eig_vec2, W, D);
if (ncut>min_cut)
    seg_out{1}=in_seg;
    return;
end

% recursive call on the A part
[seg1] = Ncuts(in_seg(A), W(A, A), min_cut);
% recursive call on the B part
[seg2] = Ncuts(in_seg(B), W(B, B), min_cut);

% return complete array of segments
seg_out=[seg1 seg2];
end