clc
clear all
%% načtení informací z DCRAW
[~,r] = system('dcraw -v -w -T 71516.CR2');
fotak=[];
for i=9:length(r)
    fotak=[fotak r(i)];
    if r(i)==' ' && r(i+1)=='i' && r(i+2)=='m' && r(i+3)=='a' && r(i+4)=='g'
        break
    end
end
f=1;
cerna=[];
for i=1:length(r)
    if r(i)=='n' && r(i+1)=='e' && r(i+2)=='s' && r(i+3)=='s' && r(i+4)==' '
        while r(f+i+4)~=','
            cerna=[cerna r(i+4+f)];
            f=f+1;
        end
    end
end
cerna = str2double(cerna);
f=1;
saturace=[];
for i=1:length(r)
    if r(i)=='t' && r(i+1)=='i' && r(i+2)=='o' && r(i+3)=='n' && r(i+4)==' '
        while r(f+i+4)~=','
            saturace=[saturace r(i+4+f)];
            f=f+1;
        end
    end
end
saturace = str2double(saturace);
f=1;
for i=1:length(r)
    if r(i)=='i' && r(i+1)=='e' && r(i+2)=='r' && r(i+3)=='s' && r(i+4)==' '
        skala_cervena=r(i+5:i+12);
        skala_zelena1=r(i+14:i+21);
        skala_modra=r(i+23:i+30);
        skala_zelena2=r(i+32:i+39);
    end
end
skala_cervena = str2double(skala_cervena);
skala_zelena1 = str2double(skala_zelena1);
skala_modra = str2double(skala_modra);
skala_zelena2 = str2double(skala_zelena2);
%% Normalizace senzorových dat
system('dcraw -4 -D -T 71516.CR2');
raw = double(imread('71516.tiff'));
raw=((raw-cerna)/(saturace-cerna));
raw_u=raw/2;
raw_o=raw*2;
%% normální snímek
R=raw(1:2:end, 1:2:end);
B=raw(2:2:end, 2:2:end);
G1=raw(1:2:end, 2:2:end);
G2=raw(2:2:end, 1:2:end);
G=(G1+G2)/2;
loc= G>1;
G(loc)=1;
loc= G<0;
G(loc)=0;

loc= R>1;
R(loc)=1;
loc= R<0;
R(loc)=0;

loc= B>1;
B(loc)=1;
loc= B<0;
B(loc)=0;

norm_G=G*(skala_zelena1); % Zde se mění množství zelené
norm_R=R*skala_cervena; % Zde se mění množství červené
norm_B=B*skala_modra; % Zde se mění množství modré
loc= norm_G>1;
norm_G(loc)=1;

loc= norm_R>1;
norm_R(loc)=1;

loc= norm_B>1;
norm_B(loc)=1;
%% podexponovaný snímek
R_u=raw_u(1:2:end, 1:2:end);
B_u=raw_u(2:2:end, 2:2:end);
G1_u=raw_u(1:2:end, 2:2:end);
G2_u=raw_u(2:2:end, 1:2:end);
G_u=(G1_u+G2_u)/2;
loc= G_u>1;
G_u(loc)=1;
loc= G_u<0;
G_u(loc)=0;

loc= R_u>1;
R_u(loc)=1;
loc= R_u<0;
R_u(loc)=0;

loc= B_u>1;
B_u(loc)=1;
loc= B_u<0;
B_u(loc)=0;

norm_G_u=G_u*(skala_zelena1); % Zde se mění množství zelené
norm_R_u=R_u*skala_cervena; % Zde se mění množství červené
norm_B_u=B_u*skala_modra; % Zde se mění množství modré
loc= norm_G_u>1;
norm_G_u(loc)=1;

loc= norm_R_u>1;
norm_R_u(loc)=1;

loc= norm_B_u>1;
norm_B_u(loc)=1;
%% přeexponoovaný snímek
R_o=raw_o(1:2:end, 1:2:end);
B_o=raw_o(2:2:end, 2:2:end);
G1_o=raw_o(1:2:end, 2:2:end);
G2_o=raw_o(2:2:end, 1:2:end);
G_o=(G1_o+G2_o)/2;
loc= G_o>1;
G_o(loc)=1;
loc= G_o<0;
G_o(loc)=0;

loc= R_o>1;
R_o(loc)=1;
loc= R_o<0;
R_o(loc)=0;

loc= B_o>1;
B_o(loc)=1;
loc= B_o<0;
B_o(loc)=0;

norm_G_o=G_o*(skala_zelena1); % Zde se mění množství zelené
norm_R_o=R_o*skala_cervena; % Zde se mění množství červené
norm_B_o=B_o*skala_modra; % Zde se mění množství modré
loc= norm_G_o>1;
norm_G_o(loc)=1;

loc= norm_R_o>1;
norm_R_o(loc)=1;

loc= norm_B_o>1;
norm_B_o(loc)=1;

%% Převedení do sRGB - normální snímek
xyz2cam= [6071 -747 -856; -7653 15365 2441; -2025 2553 7315]/10000;
srgb2xyz = [0.4124564 0.3575761 0.1804375; 0.2126729 0.7151522 0.0721750; 0.0193339 0.1191920 0.9503041];
rgb2cam = xyz2cam * srgb2xyz;
rgb2cam = rgb2cam ./ repmat(sum(rgb2cam,2),1,3);
cam2rgb = rgb2cam^-1;
r = cam2rgb(1,1)*norm_R(:,:)+cam2rgb(1,2)*norm_G(:,:)+cam2rgb(1,3)*norm_B(:,:);
g = cam2rgb(2,1)*norm_R(:,:)+cam2rgb(2,2)*norm_G(:,:)+cam2rgb(2,3)*norm_B(:,:);
b = cam2rgb(3,1)*norm_R(:,:)+cam2rgb(3,2)*norm_G(:,:)+cam2rgb(3,3)*norm_B(:,:);

rgb=cat(3,r,g,b);

xyz=applycform(rgb,makecform('srgb2xyz','AdaptedWhitePoint', whitepoint ('D55'))); %nastaveni whitepoint nic nedela
rgb=applycform(xyz,makecform('xyz2srgb','AdaptedWhitePoint', whitepoint ('D55'))); %nastaveni whitepoint nic nedela

%D55 - použití blesku z fotoaparátu jako jediné světlo
%D65 - denní světlo
%D50 - zapadající světlo
loc= rgb>1;
rgb(loc)=1;
loc= rgb<0;
rgb(loc)=0;
%% podexponovaný snímek
r_u = cam2rgb(1,1)*norm_R_u(:,:)+cam2rgb(1,2)*norm_G_u(:,:)+cam2rgb(1,3)*norm_B_u(:,:);
g_u = cam2rgb(2,1)*norm_R_u(:,:)+cam2rgb(2,2)*norm_G_u(:,:)+cam2rgb(2,3)*norm_B_u(:,:);
b_u = cam2rgb(3,1)*norm_R_u(:,:)+cam2rgb(3,2)*norm_G_u(:,:)+cam2rgb(3,3)*norm_B_u(:,:);

rgb_u=cat(3,r_u,g_u,b_u);

xyz_u=applycform(rgb_u,makecform('srgb2xyz','AdaptedWhitePoint', whitepoint ('D55'))); %nastaveni whitepoint nic nedela
rgb_u=applycform(xyz_u,makecform('xyz2srgb','AdaptedWhitePoint', whitepoint ('D55'))); %nastaveni whitepoint nic nedela

%D55 - použití blesku z fotoaparátu jako jediné světlo
%D65 - denní světlo
%D50 - zapadající světlo
loc= rgb_u>1;
rgb_u(loc)=1;
loc= rgb_u<0;
rgb_u(loc)=0;
%% přeexponovaný snímek
r_o = cam2rgb(1,1)*norm_R_o(:,:)+cam2rgb(1,2)*norm_G_o(:,:)+cam2rgb(1,3)*norm_B_o(:,:);
g_o = cam2rgb(2,1)*norm_R_o(:,:)+cam2rgb(2,2)*norm_G_o(:,:)+cam2rgb(2,3)*norm_B_o(:,:);
b_o = cam2rgb(3,1)*norm_R_o(:,:)+cam2rgb(3,2)*norm_G_o(:,:)+cam2rgb(3,3)*norm_B_o(:,:);

rgb_o=cat(3,r_o,g_o,b_o);

xyz_o=applycform(rgb_o,makecform('srgb2xyz','AdaptedWhitePoint', whitepoint ('D55'))); %nastaveni whitepoint nic nedela
rgb_o=applycform(xyz_o,makecform('xyz2srgb','AdaptedWhitePoint', whitepoint ('D55'))); %nastaveni whitepoint nic nedela

%D55 - použití blesku z fotoaparátu jako jediné světlo
%D65 - denní světlo
%D50 - zapadající světlo
loc= rgb_o>1;
rgb_o(loc)=1;
loc= rgb_o<0;
rgb_o(loc)=0;

%% filtrace normalu
R=rgb(:,:,1); G=rgb(:,:,2); B=rgb(:,:,3);

[~,loc]=findpeaks(imhist(R),'SortStr','descend','NPeaks',2,'MINPEAKHEIGHT', 10000,'MINPEAKDISTANCE',26);
if size(loc)==1
    loc=[1 loc];
end
y=imhist(R);
[~,x]=min(y(min(loc):max(loc)));
pom=(min(loc)+x-2)/256;
loc= R<pom;
R(loc)=0;

[~,loc]=findpeaks(imhist(G),'SortStr','descend','NPeaks',2,'MINPEAKHEIGHT', 10000,'MINPEAKDISTANCE',12);
if size(loc)==1
    loc=[1 loc];
end
y=imhist(G);
[~,x]=min(y(min(loc):max(loc)));
pom=(min(loc)+x-2)/256;
loc= G<pom;
G(loc)=0;

[~,loc]=findpeaks(imhist(B),'SortStr','descend','NPeaks',2,'MINPEAKHEIGHT', 10000,'MINPEAKDISTANCE',12);
if size(loc)==1
    loc=[1 loc];
end
y=imhist(B);
[~,x]=min(y(min(loc):max(loc)));
pom=(min(loc)+x-2)/256;
loc= B<pom;
B(loc)=0;
rgb=cat(3,R,G,B);
%% filtrace rgb_u
R=rgb_u(:,:,1); G=rgb_u(:,:,2); B=rgb_u(:,:,3);

[~,loc]=findpeaks(imhist(R),'SortStr','descend','NPeaks',2,'MINPEAKHEIGHT', 10000,'MINPEAKDISTANCE',26);
if size(loc)==1
    loc=[1 loc];
end
y=imhist(R);
[~,x]=min(y(min(loc):max(loc)));
pom=(min(loc)+x-2)/256;
loc= R<pom;
R(loc)=0;

[~,loc]=findpeaks(imhist(G),'SortStr','descend','NPeaks',2,'MINPEAKHEIGHT', 10000,'MINPEAKDISTANCE',12);
if size(loc)==1
    loc=[1 loc];
    y=imhist(G);
[~,x]=min(y(min(loc):max(loc)));
pom=(min(loc)+x-2)/256;
loc= G<pom;
G(loc)=0;
elseif size(loc)==2
    y=imhist(G);
[~,x]=min(y(min(loc):max(loc)));
pom=(min(loc)+x-2)/256;
loc= G<pom;
G(loc)=0;
end


[~,loc]=findpeaks(imhist(B),'SortStr','descend','NPeaks',2,'MINPEAKHEIGHT', 10000,'MINPEAKDISTANCE',12);
if size(loc)==1
    loc=[1 loc];
    y=imhist(B);
[~,x]=min(y(min(loc):max(loc)));
pom=(min(loc)+x-2)/256;
loc= B<pom;
B(loc)=0;
elseif size(loc)==2
    y=imhist(B);
[~,x]=min(y(min(loc):max(loc)));
pom=(min(loc)+x-2)/256;
loc= B<pom;
B(loc)=0;
end
rgb_u=cat(3,R,G,B);
%% filtrace rgb_o
R=rgb_o(:,:,1); G=rgb_o(:,:,2); B=rgb_o(:,:,3);

[~,loc]=findpeaks(imhist(R),'SortStr','descend','NPeaks',2,'MINPEAKHEIGHT', 10000,'MINPEAKDISTANCE',26);
if size(loc)==1
    loc=[1 loc];
end
y=imhist(R);
[~,x]=min(y(min(loc):max(loc)));
pom=(min(loc)+x-2)/256;
loc= R<pom;
R(loc)=0;

[~,loc]=findpeaks(imhist(G),'SortStr','descend','NPeaks',2,'MINPEAKHEIGHT', 10000,'MINPEAKDISTANCE',12);
if size(loc)==1
    loc=[1 loc];
end
y=imhist(G);
[~,x]=min(y(min(loc):max(loc)));
pom=(min(loc)+x-2)/256;
loc= G<pom;
G(loc)=0;

[~,loc]=findpeaks(imhist(B),'SortStr','descend','NPeaks',2,'MINPEAKHEIGHT', 10000,'MINPEAKDISTANCE',12);
if size(loc)==1
    loc=[1 loc];
end
y=imhist(B);
[~,x]=min(y(min(loc):max(loc)));
pom=(min(loc)+x-2)/256;
loc= B<pom;
B(loc)=0;
rgb_o=cat(3,R,G,B);
%% zobrazeni


imshow(rgb)
imwrite(rgb,'rgb.tiff');
imwrite(rgb_u,'rgb_u.tiff');
imwrite(rgb_o,'rgb_o.tiff');
figure(2)
imshow(rgb_u)
figure(3)
imshow(rgb_o)
files={'rgb_u.tiff','rgb.tiff','rgb_o.tiff'};
expTimes = [1/120, 1/60, 1/60];
hdr = makehdr(files, 'ExposureValues', expTimes);
hdr_rgb = tonemap(hdr);
figure; imshow(hdr_rgb)

