% sum of all values in normalized color space matrix
S = sum(h(:));
mi_r = 0;
mi_g = 0;
si_rr = 0.0;
si_gg = 0.0;

for i = 1:256
for j = 1:256
mi_r = mi_r + h(i,j) * i;
mi_g = mi_g + h(i,j) * j;
end
end
mi_r = mi_r / S;
mi_g = mi_g / S;
% requires higher precision for si calculations
rd = double(r);
gd = double(g);
for i = 1:256
for j = 1:256
si_rr = si_rr + h(i,j)*(rd(i,j) - mi_r).^2;
si_gg = si_gg + h(i,j)*(rd(i,j) - mi_g).^2;
end
end
si_gr = 0;
si_rg = 0;
for i = 1:256
for j = 1:256
si_gr = si_gr + h(i,j)*(rd(i,j) - mi_r)*(gd(i,j) - mi_g);
end
end
si_rg = si_gr/S;
si_gr = si_rg;
si = [ rd; gd];
mi = [ mi_r; mi_g];
C = [ si_rr, si_rg; si_gr, si_gg];
G = zeros(256, 256);
for i = 1:256
for j = 1:256
G(i,j) = calculateGauss(rd(i,j), gd(i,j), mi, C);
end
end
drugie = double(imread('reka.jpg'));
R2 = double(drugie(:,:,1));
G2 = double(drugie(:,:,2));
B2 = double(drugie(:,:,3));
L2 = R2 + G2 + B2;
R2 = double(drugie(:,:,1))./double(L2);
G2 = double(drugie(:,:,2))./double(L2);
r2 = uint16(R2 * 255) + 1;
g2 = uint16(G2 * 255) + 1;
[H2, W2] = size(r2);
h2 = zeros(H2, W2);


for i=1:H2
for j=1:W2
G2 = calculateGauss(double(r2(i,j)), double(g2(i,j)), mi, C);
h2(i,j) = G2;
end
end
figure();
imshow(h2inv);