%% Program pro výpočet parametrů pomocí cirkulárně symetrického Markovského okolí
%  Daniel Daněk
%  Bakalářská práce, 2012

%% upozornění - funkce počítá jeden vzorek až 60sekund, podle zvoleného okolí!

%% Program
% textura - n obrázků textury oka skádaných každý do jedné z-roviny
% pocet - pocet cirkulárně symetrických okolí kolem každého pixelu (def. pocet=3)
% mezera - vzdálenost jednotlivých okolí od sebe (def. mezera=1)
function [A,sigma]=GMRF_cirkular(textura,pocet,mezera)

if exist('textura','var')==0
    display('Chyba, vstupem není žádná textura!');
end
if exist('pocet','var')==0
    pocet=3;
end
if exist('mezera','var')==0
    mezera=1;
end

textura=double(textura);

okoli=pocet*mezera; % nejvíce vzdálené okolí bude právě v této vzdálenosti
okraj=ceil(okoli); % celé číslo, udává potřebný počet pixelů kolem odhadovaného pixelu

[~,~,z]=size(textura); % počet snímků

% předdefinování proměnných
A=zeros(z,pocet);
sigma=zeros(z,1);

for size_z=1:z
R=textura(:,:,size_z); % postupné načítání vzorků
[N,M]=size(R); % velikost vzorku
num=(N-(okraj*2))*(M-(okraj*2)); % počet polí, na kterých je možné provádět odhad
G=zeros(num,1); % matice hodnot pixelů ve středu
S=zeros(pocet,num); % skórovací matice

k=1;
for i=okraj+1:1:N-okraj
    for j=okraj+1:1:M-okraj
        G(k,1)=R(i,j);
        
        for r=mezera:mezera:okoli %prvně vnitřní kruh, potom vnější
            bodu=r*8; % počet sběrných bodů závisí na vzdálenosti okolí
            
            for alfa=0:(360/bodu):(360-(360/bodu)) % pro každý bod
                souradnice=[i-(r*sind(alfa)) j+(r*cosd(alfa))]; % souřadnice sběrného bodu

                if souradnice(1) == floor(souradnice(1)) &&...
                        souradnice(2) == floor(souradnice(2)) % pokud se bod nepotřebuje aproximaci
                    % Skórovací matice z pixelů z okolí
                    S(r/mezera,k)=S(r/mezera,k)+R(i-(r*sind(alfa)),j+(r*cosd(alfa)));
                elseif souradnice(1) > floor(souradnice(1)) || souradnice(2) > floor(souradnice(2)) % pokud se bod aproximuje
                    % poměrová vzdálenost pro váhování 4 nejbližších pixelů
                    vzdalenost(1)=sqrt((floor(souradnice(1))-souradnice(1))^2+(floor(souradnice(2))-souradnice(2))^2);
                    vzdalenost(2)=sqrt((floor(souradnice(1))-souradnice(1))^2+(ceil(souradnice(2))-souradnice(2))^2);
                    vzdalenost(3)=sqrt((ceil(souradnice(1))-souradnice(1))^2+(floor(souradnice(2))-souradnice(2))^2);
                    vzdalenost(4)=sqrt((ceil(souradnice(1))-souradnice(1))^2+(ceil(souradnice(2))-souradnice(2))^2);
                    vzdalenost=vzdalenost.^(-1);
                    vzdalenost=vzdalenost./sum(vzdalenost);

                    % Skórovací matice z pixelů z okolí
                    S(r/mezera,k)=S(r/mezera,k)+(vzdalenost(1)*R(floor(souradnice(1)),floor(souradnice(2)))+...
                        vzdalenost(2)*R(floor(souradnice(1)),ceil(souradnice(2)))+...
                        vzdalenost(3)*R(ceil(souradnice(1)),floor(souradnice(2)))+...
                        vzdalenost(4)*R(ceil(souradnice(1)),ceil(souradnice(2))));
                end
            end
            S(r/mezera,k)=S(r/mezera,k)/bodu; % pro každé okolí poděleno počtem bodů
        end
        k=k+1;
    end
end

A(size_z,:)=((S*S')^(-1))*(S*G); % Odhadovaný parametr A
sigma(size_z,:)=sqrt(sum((G'-A(size_z,:)*S).^2)/num); % Odhadovaný parametr sigma

end
end
