clear, clc, close all

file = "18031998_171.fits";   % jméno souboru
r = 2; R = 1.07;   % parametry pro mediánový filtr
gama = 1.7;   % koeficient gama korekce
r1 = 10; r2 = 15; t = 1.0;   % parametry frekvenčního filtru

%% Příprava vstupních dat

% načtení surových dat ze souboru
data = fitsread(file);
data = data(3:1022, 3:1022); % odstranění krajních pixelů
[m, n] = size(data);
load("P.mat") % lokální maxima zprůměrovaných dat

% odstranění impulzního šumu
f = medianovy_filtr(data, r, R);

% gama korekce
f = gama_korekce(f, gama);

%% Fourierova transformace
F = fft2(f);
F_shift = fftshift(F);
A_shift = abs(F_shift);

%% Nalezení lokálních maxim
tmp = P;
P = [];
for l = 1:size(tmp, 1)
    sx = tmp(l, 1); sy = tmp(l, 2);
    value = A_shift(sx, sy);
    a = sx; b = sy;
    for i = 1:m
        for j = 1:n
            if ((i-sx)^2 + (j-sy)^2 < 10^2) && A_shift(i, j) > value
                a = i; b = j;
                value = A_shift(i, j);
            end
        end
    end
    P = [P; a,b];
end

%% Filtrování
A_shift_filtered = frekvencni_filtr(A_shift, P, r1, r2, t);

% stanovení váhové funkce
H = A_shift_filtered./A_shift;

% zrcadlení váhové funkce
for i = 1:m
    for j = 1:n-i
        H(m-i+1, n-j+1) = H(i+1, j+1);
    end
end

% filtrování Fourierova spektra
G_shift = F_shift.*H;

% inverzní Fourierova transformace
G = ifftshift(G_shift);
g = ifft2(G);
g = preskalovani(real(g));

%% Zobrazení výsledků
figure
subplot(1,2,1), imshow(f), title("Původní snímek")
subplot(1,2,2), imshow(g), title("Snímek po odstranění mřížky metodou FFPS")

% zobrazení vyfiltrovaného spektra
% figure, imagesc(log(abs(G_shift) + 1)), axis equal, shading flat, colormap gray

% histogram
% figure, imhist(g)