%Skript slouží k ověření měření testovacích snímků na PC, liší se od funkce pro standardní měření s telefonem rozložením jamek 

%% LOAD RAW SOUBORŮ a proměnných
load('lab_kal_avg.mat'); %načte předem uložené proměnné transm_fault a well_empty
well_empty = well_empty_avg;
transm_fault = transm_fault_avg;

img = 'testmereni1.dng';
raw = double(imread(img));
info = imfinfo(img);

white = info(1).WhiteLevel;     %hodnota bílé
black = info(1).BlackLevel(4);  %hodnota černé
bps = info(1).BitsPerSample;    %bitová hloubka - kontrola

blueraw = raw(1:2:end,1:2:end); %výběr modrých pixelů, platné pro Bayer pattern RG/GB

%normalizace na rozsah 0-1, tento vystupuje dále jako original obraz
lin_raw = (blueraw - black) / (white - black);
lin_raw = max(0,min(lin_raw,1));
imshow(lin_raw); title('Originální obraz (RAW/DNG po extrakci modrého kanálu a normalizaci)');

%% BINARIZACE OBRAZU
%filtrace pro snažší binarizaci
kernel = ones(3, 3) / 9;                        %3x3 mean kernel
raw_filtr = medfilt2(lin_raw);                  %mediánová předfiltrace
raw_filtr2 = conv2(raw_filtr, kernel, 'same');  %průměrování konvolucí s průměrovací maskou

bw = imbinarize(raw_filtr2,0.04);   %binarizace s prahem 0.04 - empiricky určeno, pozadí vždy tmavší, prstence jamek vždy světlejší

SE = strel('disk',4);               
bw_dil = imdilate(bw,SE);           %morfo dilatace pro zlepšení tvaru binar. obrazu jamek
%imshow(bw_dil);

%% NAROVNÁNÍ SNÍMKU POMOCÍ DETEKCE HRAN A ROTACE - volání funkce rotace_hrany

[bw_rot,raw_rot] = rotace_hrany(bw_dil,lin_raw);

%% DETEKCE JAMEK
[centers, radii, metric] = imfindcircles(bw_rot,[20 40],'Sensitivity',0.75,'ObjectPolarity','bright'); %hledání kružnic
%imshow(bw_rot)
%viscircles(centers, radii,'EdgeColor','r');

%seřazení středů kružnic dle y-souřadnice a následně dle x-souřadnice tak,
%aby pořadí detekovaných jamek odpovídalo pozicím jamek v matici mikrotitrační destičky
sort_centers_y = sortrows (centers,2);  %plot(sort_centers_y(:,2))
sort_centers = [];
for k = 1:12:85
    sort_centers = [sort_centers; sortrows(sort_centers_y(k:k+11,:),1)];
end
sort_centers = round(sort_centers);

%zobrazení original obrazu se zaznačenými detekovanými jamkami a pořadím po seřazení + plnění matice středů jamek
imshow(bw_rot)
viscircles(sort_centers, repmat(mean(radii),96,1),'EdgeColor','r');
title('Originální obraz se zaznačenými detekovanými jamkami a pořadím po seřazení');
jamky = {};
for k = 1:length(sort_centers)
    centers_string = sprintf('%i',k);  
    text(sort_centers(k,1),sort_centers(k,2),centers_string,'color','b','HorizontalAlignment', 'center','VerticalAlignment', 'middle');
    jamky{k} = sort_centers(k,:); %plnění matice středů jamek
end

well_coor = reshape(jamky,12,8)';	%přetvarování a transpozice matice středů jamek


%% ČTECÍ RÁMCE A KOREKCE VZHLEDEM K ÚHLU NÁHLEDU
ST = strel('disk',7,0); %čtecí rámec
okoli = 7;
corected_centers = well_coor; %budoucí matice korigovaných středů

%rozdělení destičky na 9 segmentů, aplikace korekce na 8 krajních
%posunutí čtecího rámce víc ke středu destičky
%Levý horní segment
kompenzace_LH{1}(1) = okoli+2 ;kompenzace_LH{1}(2) = okoli+2;
for k = 1:2
    for l = 1:4
        corected_centers(k,l) = cellfun(@plus, corected_centers(k,l) , kompenzace_LH ,'Un',0);  
    end
end
%Levý segment
kompenzace_L{1}(1) = okoli ;kompenzace_L{1}(2) = 0;
for k = 3:6
    for l = 1:4
        corected_centers(k,l) = cellfun(@plus, corected_centers(k,l) , kompenzace_L ,'Un',0);  
    end
end
%Levý dolní segment
kompenzace_LD{1}(1) = okoli ;kompenzace_LD{1}(2) = -okoli;
for k = 7:8
    for l = 1:4
        corected_centers(k,l) = cellfun(@plus, corected_centers(k,l) , kompenzace_LD ,'Un',0);  
    end
end
%Horní segment
kompenzace_H{1}(1) = 0 ;kompenzace_H{1}(2) = okoli;
for k = 1:2
    for l = 5:8
        corected_centers(k,l) = cellfun(@plus, corected_centers(k,l) , kompenzace_H ,'Un',0);  
    end
end
%Dolní segment
kompenzace_D{1}(1) = 0 ;kompenzace_D{1}(2) = -okoli;
for k = 7:8
    for l = 5:8
        corected_centers(k,l) = cellfun(@plus, corected_centers(k,l) , kompenzace_D ,'Un',0);  
    end
end
%Pravý horní segment
kompenzace_PH{1}(1) = -okoli ;kompenzace_PH{1}(2) = okoli;
for k = 1:2
    for l = 9:12
        corected_centers(k,l) = cellfun(@plus, corected_centers(k,l) , kompenzace_PH ,'Un',0);  
    end
end
%Pravý segment 
kompenzace_P{1}(1) = -okoli ;kompenzace_P{1}(2) = 0;
for k = 3:6
    for l = 9:12
        corected_centers(k,l) = cellfun(@plus, corected_centers(k,l) , kompenzace_P ,'Un',0);  
    end
end
%Pravý dolní segment
kompenzace_PD{1}(1) = -(okoli+3) ;kompenzace_PD{1}(2) = -okoli;
for k = 7:8
    for l = 9:12
        corected_centers(k,l) = cellfun(@plus, corected_centers(k,l) , kompenzace_PD ,'Un',0);  
    end
end


%pomocná matice pro zobrazeni rámců v imshow
corected_centers_mat = cell2mat(corected_centers);
cor_cen_row = [];
for i = 1:2:23
    cor_cen_row = [cor_cen_row; corected_centers_mat(:,i:i+1)];
end

%kontrolní zobrazení (pozn. přibližné vyznačení čtecích rámců, nikoliv skutečná reprezentace)
figure
imshow(raw_rot); 
hold on; viscircles(sort_centers, repmat(okoli,96,1),'EdgeColor','r','LineWidth',1); hold off;  %red - původní středy
hold on; viscircles(cor_cen_row,repmat(okoli,96,1),'EdgeColor','b','LineWidth',1); hold off;    %blue- korigované středy
title('Přibližné vyznačení čtecích rámců: R=původní; B=korigované');

%% ČTENÍ HODNOT Z JAMEK
SK = strel('disk',5,0);     %čtecí rámec
kernel = ones(5, 5) / 25;   %5x5 mean kernel

for g = 1:8
    for h = 1:12
        %netvarované okolí zkoumaného pixelu - plná matice 
        pom0 = raw_rot(corected_centers{g,h}(2)-okoli : corected_centers{g,h}(2) +okoli , corected_centers{g,h}(1)-okoli : corected_centers{g,h}(1)+okoli); 
        pom1 = medfilt2(pom0,'symmetric');      %mediánová předfiltrace
        pom2 = conv2(pom1, kernel, 'valid');    %průměrování konvolucí s 5x5 maskou pro zmírnění vlivu šumu
        flek = pom2 .* double(SK.Neighborhood); %z plné matice bereme pouze část ve tvaru disku, lépe odpovídá tvaru jamky
        well_vals(g,h) = median(median(flek(flek ~=0 ))); %mediánová hodnota z disku(méně ovlivněna neodfiltrovanými extrémy než průměr)
    end
end

%% KOREKCE NEHOMOGENITY PODSVĚTLENÍ

well_vals_cleared = well_vals ./ transm_fault;     %dělení pomocí faktorů relativího útlumu

%% STANDARDNÍ VZORKY A KALIBRAČNÍ KŘIVKA - indexace dle volby rozložení jamek
stds = [ 1.25 1 .75 .5 .1875 .125 .0625 .03125 .0156 .0078 .0039] %ml / 100 ml

rr_blank = 3;       %řádkové souřadnice blanků
cc_blank = 8:9;     %sloupcové souřadnice blanků
well_blank = well_vals_cleared(rr_blank,cc_blank); %uložení hodnot 2 blanků do proměnné

%% 2 možnosti jak získat OD - z korigovaného blanku, nebo z hodnot
well_OD_A = -log10( well_vals ./ well_empty )                % 1 měření a 1 prázdné jamky
well_OD_B = -log10( well_vals_cleared ./ mean(well_blank) ) ;% OD všech jamek

well_OD_A = round(100*well_OD_A)/100; % zaokrouhlení na 2 des. místa , nebo 1000
well_OD_B = round(100*well_OD_B)/100;
%% Určení kalibračních jamek a kalibr.křivka
rr_std1 = 4;    cc_std1 = 1:11;
rr_std2 = 5;    cc_std2 = 12:-1:2;

well_stds_OD_1_A = well_OD_A(rr_std1 , cc_std1) 
well_stds_OD_2_A = well_OD_A(rr_std2 , cc_std2)

well_stds_OD_1_B = well_OD_B(rr_std1 , cc_std1) 
well_stds_OD_2_B = well_OD_B(rr_std2 , cc_std2)

well_stds_OD_A = mean( [well_stds_OD_1_A; well_stds_OD_2_A], 1 )
well_stds_OD_B = mean( [well_stds_OD_1_B; well_stds_OD_2_B], 1 )


calibr_curve_A = polyfit(stds,well_stds_OD_A,1); %kalibrační křivka - polynom prvního stupně -> vždy lineární
OD_computed_A = polyval(calibr_curve_A,stds); %vykreslení pro kontrolu

calibr_curve_B = polyfit(stds,well_stds_OD_B,1); %kalibrační křivka - polynom prvního stupně -> vždy lineární
OD_computed_B = polyval(calibr_curve_B,stds); %vykreslení pro kontrolu

%% SUMA ČTVERCŮ ODCHYLEK - naměřené OD, předpokládáné OD
RA = corrcoef(OD_computed_A,well_stds_OD_A); R2A = RA(2)^2;
RB = corrcoef(OD_computed_B,well_stds_OD_B); R2B = RB(2)^2;

%to samé ručně
resid = well_stds_OD_A - OD_computed_A;
SSresid = sum(resid.^2);
SStotal = (length(well_stds_OD_A)-1) * var(well_stds_OD_A);
rsq = 1 - SSresid/SStotal
%disp SSresid/SStotal
%disp rsq
resid = well_stds_OD_B - OD_computed_B;
SSresid = sum(resid.^2);
SStotal = (length(well_stds_OD_B)-1) * var(well_stds_OD_B);
rsq = 1 - SSresid/SStotal


figure(1)
plot(stds,OD_computed_A,'r'); hold on;
plot(stds,OD_computed_B,'-.k'); hold on;
scatter(stds,well_stds_OD_A,'r');hold on;
scatter(stds,well_stds_OD_B,'k');hold off;
str = ['A:   R^2 = ' num2str(R2A) ' B:   R^2 = ' num2str(R2B) ];
str2 =['A:   ' 'y = ' num2str(calibr_curve_A(1)) '*x + ' num2str(calibr_curve_A(2))]; 
str3 =['B:   ' 'y = ' num2str(calibr_curve_B(1)) '*x + ' num2str(calibr_curve_B(2))];
xlim=get(gca,'XLim');
ylim=get(gca,'YLim');
text(0.1*xlim(1)+0.1*xlim(2),0.3*ylim(1)+0.9*ylim(2),str);
text(0.1*xlim(1)+0.1*xlim(2),0.3*ylim(1)+0.8*ylim(2),str2);
text(0.1*xlim(1)+0.1*xlim(2),0.3*ylim(1)+0.77*ylim(2),str3);
title('Kalibrační křivky - červená = A, černá = B');
xlabel('Koncentrace [ml / 100 ml]') 
ylabel('OD') 
legend('Kalibrační křivka A','Kalibrační křivka B','Location','Southeast')


%% MATICE KONCENTRACÍ
C_complete_A = zeros(8,12); C_complete_B = zeros(8,12); 

C_complete_A = (well_OD_A - calibr_curve_A(2)) ./ calibr_curve_A(1);    % matice koncentrací vypočtených z kalibrační křivky
C_complete_B = (well_OD_B - calibr_curve_B(2)) ./ calibr_curve_B(1); 

C_complete_A(find(well_OD_A <=0)) = 0;                              %nuluje koncentraci, kde absorbance <=0
C_complete_A(C_complete_A<0) = 0;   
C_complete_B(find(well_OD_B <=0)) = 0;                              %nuluje koncentraci, kde absorbance <=0
C_complete_B(C_complete_B<0) = 0;                                   %nuluje koncentraci, koncentrace <=0
%C_completeC = round(10000*C_complete)/100;                         %případné zaokrouhlení


%% ZOBRAZENÍ VÝSLEDKŮ - zobrazuje hodnoty zaokrouhlené na 5 desetinných míst
C_complete_A = round(100000*C_complete_A)/100000;
C_complete_B = round(100000*C_complete_B)/100000;


figure(2)
subplot 211
imagesc(C_complete_A);colormap(flipud(bone))
set(gca,'XAxisLocation','top','xgrid', 'on', 'ygrid', 'on','XTick',1:12,'YTickLabel',{'A','B','C','D','E','F','G','H',})
title('Mapa koncentrací A [ml / 100 ml]')             
for i = 1:8
    for j = 1:12
        text(j,i ,num2str(C_complete_A(i,j)),'FontSize',15,'Color','r','HorizontalAlignment','center');
    end
end

subplot 212
imagesc(C_complete_B);colormap(flipud(bone))
set(gca,'XAxisLocation','top','xgrid', 'on', 'ygrid', 'on','XTick',1:12,'YTickLabel',{'A','B','C','D','E','F','G','H',})
title('Mapa koncentrací B [ml / 100 ml]')             
for i = 1:8
    for j = 1:12
        text(j,i ,num2str(C_complete_B(i,j)),'FontSize',15,'Color','r','HorizontalAlignment','center');
    end
end


