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

%% Program
% textura - n obrázků textury oka skádaných každý do jedné z-roviny
% typ - odpovídá rozložení a řádu okolí
function [A,sigma]=GMRF_normal(textura,typ)

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

textura=double(textura);

switch typ % přiřazení ke každému typu okolí a počet
    case 1
        okoli=1; % nejvíce vzdálené okolí bude právě v této vzdálenosti
        pocet=1; % pocet parametrů, které se bude odhadovat
    case 2
        okoli=1;
        pocet=2;
    case 3
        okoli=2;
        pocet=3;
    case 4
        okoli=2;
        pocet=4;
    case 5
        okoli=3;
        pocet=5;
    case 6
        okoli=3;
        pocet=6;
end

[~,~,z]=size(textura);

% 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-(okoli*2))*(M-(okoli*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=okoli+1:1:N-okoli
    for j=okoli+1:1:M-okoli
        G(k,1)=R(i,j);
 
        switch typ % podle typu se počítá skórovací matice z pixelů z okolí
            case 1
                S(1,k)=R(i-1,j)+R(i+1,j)-R(i,j-1)-R(i,j+1);
            case 2
                S(1,k)=R(i-1,j)+R(i+1,j)+R(i,j-1)+R(i,j+1);
                S(2,k)=R(i-1,j+1)+R(i+1,j-1)+R(i-1,j-1)+R(i+1,j+1);
            case 3
                S(1,k)=R(i-1,j)+R(i+1,j)+R(i,j-1)+R(i,j+1);
                S(2,k)=R(i-1,j-1)+R(i-1,j+1)+R(i+1,j-1)+R(i+1,j+1);
                S(3,k)=R(i-2,j)+R(i+2,j)+R(i,j-2)+R(i,j+2);
            case 4
                S(1,k)=R(i-1,j)+R(i+1,j)+R(i,j-1)+R(i,j+1);
                S(2,k)=R(i-1,j-1)+R(i-1,j+1)+R(i+1,j-1)+R(i+1,j+1);
                S(3,k)=R(i-2,j)+R(i+2,j)+R(i,j-2)+R(i,j+2);
                S(4,k)=R(i-1,j-2)+R(i-2,j-1)+R(i-2,j+1)+R(i-1,j+2)+...
                    R(i+1,j-2)+R(i+2,j-1)+R(i+2,j+1)+R(i+1,j+2);
            case 5
                S(1,k)=R(i-1,j)+R(i+1,j)+R(i,j-1)+R(i,j+1);
                S(2,k)=R(i-1,j-1)+R(i-1,j+1)+R(i+1,j-1)+R(i+1,j+1);
                S(3,k)=R(i-2,j)+R(i+2,j)+R(i,j-2)+R(i,j+2);
                S(4,k)=R(i-1,j-2)+R(i-2,j-1)+R(i-2,j+1)+R(i-1,j+2)+...
                    R(i+1,j-2)+R(i+2,j-1)+R(i+2,j+1)+R(i+1,j+2);
                S(5,k)=R(i-3,j)+R(i+3,j)+R(i,j-3)+R(i,j+3);
            case 6
                S(1,k)=R(i-1,j)+R(i+1,j)+R(i,j-1)+R(i,j+1);
                S(2,k)=R(i-1,j-1)+R(i-1,j+1)+R(i+1,j-1)+R(i+1,j+1);
                S(3,k)=R(i-2,j)+R(i+2,j)+R(i,j-2)+R(i,j+2);
                S(4,k)=R(i-1,j-2)+R(i-2,j-1)+R(i-2,j+1)+R(i-1,j+2)+...
                    R(i+1,j-2)+R(i+2,j-1)+R(i+2,j+1)+R(i+1,j+2);
                S(5,k)=R(i-3,j)+R(i+3,j)+R(i,j-3)+R(i,j+3);
                S(6,k)=R(i-3,j+1)+R(i-2,j+2)+R(i-1,j+3)+...
                    R(i+1,j+3)+R(i+2,j+2)+R(i+3,j+1)+...
                    R(i+3,j-1)+R(i+2,j-2)+R(i+1,j-3)+...
                    R(i-1,j-3)+R(i-2,j-2)+R(i-3,j-1);
        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