% function par = imncut; or
%          [X,V,S,emag,ephase,W] = imncut(I,par)
% Input:
%   I = image
%   par = parameters
%       par.nb_r = neighbour hood radius
%       par.sample_rate = sampling rate
%       par.mask = mask for bias, label matrix of image size
%       par.sig_ic = sigma for computing affinity from intervening contour
%       par.rep = offset for repulsion
%       par.reg = offset for regularization
%       par.nv = number of eigenvectors
% Output:
%   X = K-way partition matrix
%   [V,S] = eigenvector and eigenvalues from constrained ncuts
%   [emag,ephase] = edge detection magnitudes and phases
%   W = affinity matrix
%   figure display of K-way discrete segmentation
% Reference:
%   www.cs.berkeley.edu/~stellayu/thesis.html
% Examples:
%   par = imncut; % return the default parameters
%   imncut(I); % ncuts using the default parameters

% Stella X. Yu, July 13 2003.

function [X,V,S,emag,ephase,W] = imncut(I,par)

if nargin==0,
    par.nb_r = 8;
    par.sample_rate = 0.15;
    par.mask = [];
    par.sig_ic = 0.15;
    par.rep = 0;
    par.reg = 0;
    par.nv = 3;  
    X = par;
    return;
end

if nargin<2,
    par = imncut;
end
nb_r = par.nb_r;
sample_rate = par.sample_rate;
mask = par.mask;
sig_ic = par.sig_ic;
rep = par.rep;
reg = par.reg;
nv = par.nv;


% sampling
sz = size(I); np = prod(sz);
if length(sample_rate)==1,
    samp = double(rand(sz)<sample_rate);
else
    samp = sample_rate;
end

% pixel affinity
[emag,ephase] = quadeg(I);
[ci,cj] = imnb(nb_r,samp);
cw = ic(emag,ephase,ci,cj,mask); 
if abs(rep)<eps, % attraction
    offset = 0.005;
else % repulsion
    offset = -rep;
end
a = 0.5 / (sig_ic * sig_ic);
W = exp(-cw * a) + offset;
W = csparse(ci,cj,W,np);            
clear ci cj cw

% constraints
p = find(samp);
if isempty(mask),
    U = [];
else
    U = pargrp(double(mask & samp));
    U = U(p,:);
end

% cuts on samples; and then interpolate to the rest pixels
[Vp,S] = cncut(W(p,p),U,nv);
q = find(not(samp));
if isempty(q),
    Vq = [];
else
    C = W(q,p);
    C = spmd(C,1./double(sum(C,2)+eps),[]);
    Vq = C * Vp;
end
V = zeros(np,nv);
V([p;q],:) = [Vp; Vq];
for j=1:nv,
    z = reshape(V(:,j),sz);
    z = medfilt2(z,[7,7],'symmetric');
    V(:,j) = z(:);
end

% discretization
X = getbinsol(V); 

% visualization
figure;
range = [min(V(:)),max(V(:))];
for j=1:nv,
    z = reshape(V(:,j),sz);
    subplot(nv,2,j+j-1); imagesc(z,range); axis image; axis off;
    colormap(gray);
    title(sprintf('%5.4f',S(j)));
    z = reshape(X(:,j),sz);
    subplot(nv,2,j+j); imagesc(I.*z+not(z),[0,1]); axis image; axis off;
    colormap(gray);
end