clear all
close all
clc
%%                  DETEKCE VNITRNI HRANICE DUHOVKY - snimky lamp
% vstupem sedotonovy snimek
% vystupem souradnice stredu kruznice, polomer a souradnice kruznice
% cast Houghovy transformace prevzat od Ing. Podlipne
%% 1. nacteni snimku
snimky = ziskaniSouboru('\Tobiasova_DP\SNI\LAMP\original');
pocet = length(snimky);

for s = 33% 1:pocet ;
 I = imread(cell2mat(snimky(s)));
 figure (1)
 imshow(I); title ('original vyrez')

%% 2. histogram, nalezeni 1. minima
figure (2); imhist(I); title('histogram vyrezu') % histogram
minimum=[]; % pomocna promenna pro ulozeni minima
pom=1;      % pomocna promenna
[pocetPixelu bins] = imhist(I,256); % pocet pixelu v jednotlivych stupnich
%hledani prvniho lokalniho minima
%prvek mensi nez prvek pred nim a za nim

for i = 3:size(pocetPixelu,1)-2
    if pocetPixelu(i)<pocetPixelu(i-1) && pocetPixelu(i)<pocetPixelu(i+1) && pocetPixelu(i+1)<pocetPixelu(i+2)&& pocetPixelu(i-1)<pocetPixelu(i-2) && pocetPixelu(i)>20
        minimum(pom)=pocetPixelu(i); % zapis vsech lokalnich minim
        pom=pom+1;
    end
end
prahx = find(pocetPixelu(3:end)==minimum(1,1)); % 1. lokalni minimum jako PRAH
prah = 2 + prahx(1,1); % z duvodu prohledavani hodnot od 3 

%% 3. prahovani pomoci indexovani
prahovany = zeros(size(I)); % velikost vystupu
prahovany(I>prah)=0;       % nevyznamna oblast
prahovany(I<=prah)=1;       % vyznamna oblast
figure (3); imshow(prahovany); title('prahovany snimek')

%% 4. Detekce hran
Filt  = medfilt2(prahovany, [5 5]);
Edges = edge(Filt, 'sobel');
figure(4); imshow(Edges, []) ; title('detekce hran')

%% 5. Houghova transformace - detekce hranic zornice (vnitrni hranice duhovky)
Radius = 27:2:65; % poloměr
Accumulator = zeros(size(Edges,1),size(Edges,2),size(Radius,2));
% preddefinovani akumulatoru velikosti x, y, z; 1-sloupec,2-radek
[PointsX PointsY] = find(Edges); 
UnitCircleX = cos((1:0.5:360)*pi/180);
UnitCircleY = sin((1:0.5:360)*pi/180);

for r = 1:size(Radius,2)
    for j = 1:length(PointsX)                 
        a = round(PointsX(j) - Radius(r)*UnitCircleX); % Výpočet kružnic
        b = round(PointsY(j) - Radius(r)*UnitCircleY);
        for i = 1:360
            if a(i)>0 && a(i)<size(Edges,1) && b(i)>0 && b(i)<size(Edges,2) % Jestli souřadnice bodu leží uvnitř obrazu
                Accumulator(a(i),b(i),r) = Accumulator(a(i),b(i),r)+1;      % Inkrementujeme akumulátor
            end
        end    
    end
end
HoughSpace = mat2gray(Accumulator); % interval 0-1

%% 6. Vyhodnocení 
[maxHoughSpace idx] = max(HoughSpace(:)); % zjisteni max hodnoty 3D matice
[n m t] = ind2sub(size(HoughSpace),idx);
% figure (5); imshow(HoughSpace(:,:,t),[]); title('Hough. prostor') % zobrazeni hough. prostoru
% pom=HoughSpace(:,:,t);
% image2=imresize(I, 0.5); % podvzorkovani puvodniho snimku
figure(6);imshow(I,[]); title('Hough-Lamp')
hold on
x = m + Radius(t)*UnitCircleX;
y = n + Radius(t)*UnitCircleY;
plot( x, y, '.g', 'MarkerSize', 10)


% path='I:\Users\tobiasova\pokracovani\IMAGES\VNITRNI_hough_INTERVAL_podm\';
% saveas(f,[path,'vnitrni_hough',num2str(s),'.png']);
end

