cnsslis2.m

% Complex Networks Semi-Supervised Learning Image Segmentation v2
% Trabalha na primeira fase com imagem redimensionada para um 1/9 do tamanho
% original. Inclui ExR, ExB, e ExG. Exclui desvios padrões (v.2)
% Usage: [owner, pot] = cnsslis(img, imgslab, k, fw, disttype, valpha, maxiter)
% INPUT:
% img       - Image to be segmented
% imgslab   - Image with labeled/unlabeled pixel information
% k         - each node is connected to its k-neirest neighbors
% fw        - vector of feature weights
% disttype  - use 'euclidean', etc.
% valpha    - Default: 20 (lower it to stop earlier, accuracy may be lower)
% maxiter   - maximum amount of iterations
% OUTPUT:
% owner     - vector of classes assigned to each data item
% pot     

function [owner, pot] = cnsslis2(img, imgslab, fw, k, disttype, valpha, maxiter)
if (nargin < 7) || isempty(maxiter)
    maxiter = 500000; % número de iterações
end
if (nargin < 6) || isempty(valpha)
    valpha = 20;
end
if (nargin < 5) || isempty(disttype)
    disttype = 'euclidean'; % distância euclidiana não normalizada
end
if (nargin < 4) || isempty(k)
    k = 8; % quantidade de vizinhos mais próximos
end
if (nargin < 3) || isempty(fw)
    fw = ones(1,20);
end
% tratamento da entrada
k = uint16(k);

if k>0
    % reduzindo imagem
    rs_img = imresize(img,1/3,'bilinear');
    rs_imgslab = imresize(imgslab,1/3,'nearest');
    [rs_dim,qtnode,X,slabel,nodeval,nclass] = getFeatures(rs_img,rs_imgslab,fw);
    
    % já estamos normalizando de qualquer forma
    if strcmp(disttype,'seuclidean')==1
        disttype='euclidean';
    end
    
    indval = find(nodeval);     % pega só os índices dos pixels que não são do fundo ignorado
    Xval = X(indval,:);         % cria lista de pixels válidos (que não são do fundo ignorado)
    qtnodeval = size(indval,1); % quantidade de nós válidos (pixels válidos)
    slabelval = slabel(indval); % rótulos dos pixels válidos (não são do fundo ignorado)    
    
    nnonlabeled = sum(slabelval==0); % quantidade de nós não rotulados
    
    stopmax = round((qtnodeval/nnonlabeled)*round(valpha*0.1)); % qtde de iterações para verificar convergência
    
    % encontrando k-vizinhos mais próximos
    KNN = uint32(knnsearch(Xval,Xval,'K',k+1,'NSMethod','kdtree','Distance',disttype));
    clear XVal;
    KNN = KNN(:,2:end); % eliminando o elemento como vizinho de si mesmo
    KNN(:,end+1:end+k) = 0; % adicionando mais k espaços para vizinhança recíproca
    knns = repmat(k,qtnodeval,1); % vetor com a quantidade de vizinhos de cada nó
    for i=1:qtnodeval
        %KNNR(sub2ind(size(KNNR),KNN(i,:),(knns(KNN(i,:))+1)'))=i; % adicionando i como vizinho dos vizinhos de i (criando reciprocidade)
        KNN(sub2ind(size(KNN),KNN(i,1:k),(knns(KNN(i,1:k))+1)'))=i; % adicionando i como vizinho dos vizinhos de i (criando reciprocidade)
        knns(KNN(i,1:k))=knns(KNN(i,1:k))+1; % aumentando contador de vizinhos nos nós que tiveram vizinhos adicionados
        if max(knns)==size(KNN,2) % se algum nó atingiu o limite de colunas da matriz de vizinhança recíproca teremos de aumentá-la
            % portanto vamos aumenta-la em 10% + 1 (para garantir no caso do tamanho ser menor que 10)
            % convém que o aumento seja de várias colunas de uma vez, pois
            % é uma operação custosa, mas não muitas para não ocupar muita
            % memória
            KNN(:,max(knns)+1:round(max(knns)*1.1)+1) = zeros(qtnodeval,round(max(knns)*0.1)+1);
        end
    end
    % removendo duplicatas
    for i=1:qtnodeval
        knnrow = unique(KNN(i,:),'stable'); % remove as duplicatas
        knns(i) = size(knnrow,2)-1; % atualiza quantidade de vizinhos (e descarta o zero no final)
        KNN(i,1:knns(i)) = knnrow(1:end-1); % copia para matriz KNN
    end
    clear knnrow;
    KNN = KNN(:,1:max(knns)); % eliminando colunas que não tem vizinhos válidos

    % lista de nós não rotulados
    indnonlabeled = uint32(find(slabelval==0));
    % lista de nós rotulados
    labelednodes = uint32(find(slabelval>0));
    % ajustando todas as distâncias na máxima possível
    potval = repmat(0.5,qtnodeval,nclass);   
    % zerando potenciais dos nós rotulados
    potval(labelednodes,:) = 0;
    % ajustando potencial da classe respectiva do nó rotulado para máximo
    potval(sub2ind(size(potval),labelednodes,slabelval(labelednodes))) = 1;
    % variável para guardar máximo potencial mais alto médio
    potval = cnsslis2loop(maxiter,nnonlabeled,indnonlabeled,stopmax,potval,knns,KNN);
    clear KNN slabelval;
           
    pot = repmat([1 0],qtnode,1);
    pot(indval,:)=potval;
    
    clear potval;
end

[dim,qtnode,X,slabel,~,nclass] = getFeatures(img,imgslab,fw);
% Redimensionar matriz de potenciais
% (antes de redimensionar é preciso passar para matriz de 3 dimensões e
% depois voltar para o formato anterior)
if k>0
    pot = reshape(imresize(reshape(pot,rs_dim(1),rs_dim(2),2),[dim(1) dim(2)],'bilinear'),qtnode,2);
else
    pot = repmat(0.5,qtnode,nclass);
end

% encontrando nos rotulados
labelednodes = find(slabel>0);
% zerando potenciais dos nós rotulados
pot(labelednodes,:) = 0;
% ajustando potencial da classe respectiva do nó rotulado para 1
pot(sub2ind(size(pot),labelednodes,slabel(labelednodes))) = 1;

% PARTE 2!
%disp('Parte 2: Encontrando vizinhos...');
indefnodesb = max(pot,[],2)<0.9; % vetor onde 1 é nó indefinido e 0 é definido
indefnodes = uint32(find(indefnodesb)); % lista de nós indefinidos
indefnodesc = size(indefnodes,1); % contagem de nós indefinidos
if indefnodesc>0
    
    %fprintf('Parte 2: %i nós indefinidos. Pegando colaboração de pixels vizinhos\n',size(indefnodes,1))
    
    Ndist = zeros(size(X,1),8);
    Nlist = zeros(size(X,1),8,'uint32');
    Nsize = zeros(size(X,1),1,'uint8');
    % Pesos das ligações horizontais
    for i=1:dim(1)
        for j=1:dim(2)-1
            ind1 = i+(j-1)*dim(1);
            ind2 = ind1 + dim(1);
            if indefnodesb(ind1) || indefnodesb(ind2)
                p2addNeighbor;
            end
        end
    end
    % Peso das ligações diagonais (\)
    for i=1:dim(1)-1
        for j=1:dim(2)-1
            ind1 = i+(j-1)*dim(1);
            ind2 = ind1+dim(1)+1;
            if indefnodesb(ind1) || indefnodesb(ind2)
                p2addNeighbor;
            end
        end
    end
    % Peso das ligações verticais
    for i=1:dim(1)-1
        for j=1:dim(2)
            ind1 = i+(j-1)*dim(1);
            ind2 = ind1+1;
            if indefnodesb(ind1) || indefnodesb(ind2)
                p2addNeighbor;
            end
        end
    end
    % Peso das ligações diagonais (/)
    for i=1:dim(1)-1
        for j=2:dim(2)
            ind1 = i+(j-1)*dim(1);
            ind2 = ind1-dim(1)+1;
            if indefnodesb(ind1) || indefnodesb(ind2)
                p2addNeighbor;
            end
        end
    end
    clear X;
    % Ajustando distâncias para intervalo 0 - 1 e invertendo (convertendo em peso de aresta);
    Ndist = 1 - Ndist/max(max(Ndist));
    % constantes
    npart = indefnodesc; % quantidade de nós ainda não rotulados
    stopmax = round((qtnode/npart)*round(valpha*0.1)); % qtde de iterações para verificar convergência
    % variável para guardar máximo potencial mais alto médio
    % chamando o arquivo mex do strwalk25
    %disp('Parte 2: Propagação de rótulos...');
    pot = strwalk25loop(maxiter, npart, nclass, stopmax, indefnodes, slabel, Nsize, Nlist, Ndist, pot);
end
[~,owner] = max(pot,[],2);

    function p2addNeighbor
        Nsize(ind1) = Nsize(ind1) + 1;
        Nsize(ind2) = Nsize(ind2) + 1;
        Ndist(ind1,Nsize(ind1)) = norm(X(ind1,:)-X(ind2,:));
        Ndist(ind2,Nsize(ind2)) = Ndist(ind1,Nsize(ind1));
        Nlist(ind1,Nsize(ind1)) = ind2;
        Nlist(ind2,Nsize(ind2)) = ind1;
    end

end

function [dim,qtnode,X,slabel,nodeval,nclass] = getFeatures(img,imgslab,fw)

% Atenção: Atributo Linha e HSV estão errados em todas as versões anteriores deste algoritmo!

% Dimensões da imagem
dim = size(img);
qtnode = dim(1)*dim(2);
X = zeros(qtnode,20);
% primeiro e segundo elementos são linha e coluna normalizadas no intervalo 0:1
X(:,1:2) = [repmat(((1:dim(1))/dim(1))',dim(2),1), reshape(repmat((1:dim(1))/dim(1),dim(2),1),dim(1)*dim(2),1)]; 
% depois vem os 3 elementos RGB normalizados em 0:1
imgvec = double(squeeze(reshape(img,dim(1)*dim(2),1,3)))/255;
X(:,3:5) = imgvec;
% depois vem os 3 elementos HSV
imghsv = rgb2hsv(double(img)/255);
X(:,6:8) = squeeze(reshape(imghsv,dim(1)*dim(2),1,3));
% em seguida ExR, ExG, e ExB
exr = 2.*double(img(:,:,1)) - double(img(:,:,2)) - double(img(:,:,3));
exg = 2.*double(img(:,:,2)) - double(img(:,:,1)) - double(img(:,:,3));
exb = 2.*double(img(:,:,3)) - double(img(:,:,1)) - double(img(:,:,2));
imgex = cat(3, exr, exg, exb);
clear exr exg exb;
X(:,9:11) = squeeze(reshape(imgex,dim(1)*dim(2),1,3));
% médias
h = fspecial('average', [3 3]);
g = imfilter(img, h,'replicate'); % adicionado replicate para que bordas não fiquem diferentes
X(:,12:14) = double(squeeze(reshape(g,dim(1)*dim(2),1,3)))/255;
g = imfilter(imghsv, h, 'replicate'); % adicionado replicate para que bordas não fiquem diferentes)
X(:,15:17) = double(squeeze(reshape(g,dim(1)*dim(2),1,3)));
g = imfilter(imgex, h, 'replicate'); % adicionado replicate para que bordas não fiquem diferentes)
X(:,18:20) = double(squeeze(reshape(g,dim(1)*dim(2),1,3)));
clear g imghsv imgex;
% s = stdfilt(img);
% X(:,18:20) = double(squeeze(reshape(s,dim(1)*dim(2),1,3)))/255;
% s = stdfilt(rgb2hsv(img));
% X(:,21:23) = double(squeeze(reshape(s,dim(1)*dim(2),1,3)));
% clear s;
% normalizando as colunas
X = zscore(X) .* repmat(fw,qtnode,1);
% Converter imagem com rótulos em vetor de rótulos
slabel = uint16(reshape(imgslab,dim(1)*dim(2),1));
% montar vetor onde 0 é nó do fundo não considerado e 1 é nó válido
nodeval = zeros(qtnode,1);
nodeval(slabel~=0)=1;
% ajustar vetor de rótulos
slabel(slabel==0)=1; % fundo não considerado
slabel(slabel==64)=1;  % c/ rótulo - fundo
slabel(slabel==255)=2; % c/ rótulo - objeto
slabel(slabel==128)=0; % sem rótulo
nclass = 2;
end