Untitled

%% isValid.m

function result = isValid (H, W, border, Y)
result = 0;
quants = size(border);
shots = zeros (size(border, 1)-1, 1);
for i = 1:H
    for j = 1:W
        for k = 1:quants
            if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                shots(k) = 1;
            end
        end
    end
end
if min(shots) == 1
    result = 1;
end
end


%% newBorder.m

function border = newBorder(H, W, border, Y)
quants = size(border, 1) - 1;
shots = zeros (size(border, 1)-1, 1);
plain = zeros (H*W, 1);
for i = 1:H
    for j = 1:W
        for k = 1:quants
            plain(i*W + j) = Y(i, j);
            if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                shots(k) = 1;
            end
        end
    end
end
[flag, pos] = min(shots);
m = mean(plain);
M = fix(quants/2); %интенсивность сдвига
while flag == 0
    if m > border(pos)
        M = quants - pos - 2;
        for i = 1:M
            border(pos+i) = border(pos+i)+(M-i)/M;
        end
    else
        M = pos-1;
        for i = 1:M
            border(pos+1-i) = border(pos+1-i)-(M-i)/M*(border(pos+1-i)-border(pos-i));
        end
    end
    shots = zeros (size(border, 1)-1, 1);
    for i = 1:H
        for j = 1:W
            for k = 1:quants
                if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                    shots(k) = 1;
                end
            end
        end
    end
    for i = 1:quants
        if border(i) > border(i + 1)
            disp('ERROR!');
        end
    end
    [flag, pos] = min(shots);
end
end


%% Основной код, собсна

clear;
close all;
MAXR = 7;
picture = imread('C:\Users\Tuhlomon\Desktop\Bright_2.bmp');
H = size(picture, 1);
W = size(picture, 2);
Y = zeros(H, W);
for i = 1:H
    for j = 1:W
        Y(i, j) = picture(i, j, 1)*0.299 + picture(i, j, 2)*0.587 + picture(i, j, 3)*0.114;
    end
end

U = zeros(H, W, MAXR);
for R = 1:MAXR
    quants = 2^R;
    minimum = min(min(Y));
    maximum = max(max(Y));
    border = zeros (quants+1, 1);
    border(1) = -0.000001;
    approx = zeros (quants, 1);
    delta = (maximum-minimum)/quants;
    for i = 1:quants
        border(i+1) = minimum + i*delta;
        approx(i) = (border(i)+border(i+1))/2;
    end
    border(quants+1) = 255.000001;
    valid = isValid(H, W, border, Y);
    disp (valid);
    if (valid == 0)
        border = newBorder(H, W, border, Y);
        valid = isValid(H, W, border, Y);
        disp ('new valid = ');
        disp (valid);
    end
    Teps = 0.0001;
    DELTAeps = 0.0001;
    epsilonPred = 0;
    pmax = 20;
    newU = zeros(H, W);
    for i = 1:H
        for j = 1:W
            newU = 255;
        end
    end
    %Алгоритм Макса-Ллойда
    for P = 1:pmax
        Q = zeros (quants, 1);
        average = zeros (quants, 1);
        for i = 1:H
            for j = 1:W
                index = quants;
                for k = 1:quants
                    if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                        index = k;
                        newU(i, j) = approx(k);
                        break
                    end
                end
                Q(index) = Q(index)+1;
                average(index) = average(index)+Y(i,j);
            end
        end
        average = average ./ Q;
        approx = average;
        epsilon = 0;
        for i = 1:H
            for j = 1:W
                epsilon = epsilon + (newU(i, j) - Y(i, j))^2;
            end
        end
        epsilon = epsilon / (H * W);
        if epsilon < Teps
            break
        end
        if abs(epsilon - epsilonPred) < DELTAeps
            break
        end
        for i = 2:quants
            border(i) = (approx(i-1)+approx(i))/2;
        end
        epsilonPred = epsilon;
    end
    %%% квантование
    for i = 1:H
        for j = 1:W
            for k = 1:quants
               if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                   U(i, j, R) = approx(k);
                   break;
               end
            end
        end
    end
end

figure(1);
approx(quants+1) = approx(quants);
stairs(border, approx);

SNR = zeros (MAXR, 1);
A = Y.^2;
for R = 1:MAXR
    B = (Y-U(:, :, R)).^2;
    SNR(R) = 10*log10(sum(A)/sum(B));
end
t = 1:1:MAXR;
figure(2);
plot(t, SNR), legend('SNR(R) bright');
hold on;

PSNR = zeros (MAXR, 1);
for R = 1:MAXR
    epsilon = 0;
    for i = 1:H
        for j = 1:W
            epsilon = epsilon + (U(i, j, R) - Y(i, j))^2;
        end
    end
    PSNR(R) = 10*log10(W*H*(255^2)/epsilon);
end
figure(3);
plot(t, PSNR), legend('PSNR(R) bright');
hold on;

He = zeros (MAXR, 1);
for R = 1:MAXR
    p = zeros(2^R, 1);
    for i = 1:H
        for j = 1:W
            for k = 1:2^R
               if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                   p(k) = p(k) + 1;
                   break;
               end
            end
        end
    end
    for i = 1:2^R
        if p(i) ~= 0
            p(i) = p(i)/(H*W);
            He(R) = He(R) + (p(i)*log2(p(i)));
        end
    end
    He(R) = -He(R);
end
figure(4);
plot(He, PSNR), legend('PSNR(H) bright');
hold on;

imwrite(U(:,:,1)./255, 'C:\Users\Tuhlomon\Desktop\учеба\цос\tests2\R1.bmp');
imwrite(U(:,:,2)./255, 'C:\Users\Tuhlomon\Desktop\учеба\цос\tests2\R2.bmp');
imwrite(U(:,:,3)./255, 'C:\Users\Tuhlomon\Desktop\учеба\цос\tests2\R3.bmp');
imwrite(U(:,:,4)./255, 'C:\Users\Tuhlomon\Desktop\учеба\цос\tests2\R4.bmp');
imwrite(U(:,:,5)./255, 'C:\Users\Tuhlomon\Desktop\учеба\цос\tests2\R5.bmp');
imwrite(U(:,:,6)./255, 'C:\Users\Tuhlomon\Desktop\учеба\цос\tests2\R6.bmp');
imwrite(U(:,:,7)./255, 'C:\Users\Tuhlomon\Desktop\учеба\цос\tests2\R7.bmp');


picture = imread('C:\Users\Tuhlomon\Desktop\dark_2.bmp');
H = size(picture, 1);
W = size(picture, 2);
Y = zeros(H, W);
for i = 1:H
    for j = 1:W
        Y(i, j) = picture(i, j, 1)*0.299 + picture(i, j, 2)*0.587 + picture(i, j, 3)*0.114;
    end
end

p = zeros(256, 1);
for i = 1:H
    for j = 1:W
        for k = 1:256
            if Y(i, j) == k
                p(k) = 1;
                break;
            end
        end
    end
end
sump = sum(p);
MAXR = 6;

U = zeros(H, W, MAXR);
for R = 1:MAXR
    quants = 2^R;
    minimum = min(min(Y));
    maximum = max(max(Y));
    border = zeros (quants+1, 1);
    border(1) = -0.000001;
    approx = zeros (quants, 1);
    delta = (maximum-minimum)/quants;
    for i = 1:quants
        border(i+1) = minimum + i*delta;
        approx(i) = (border(i)+border(i+1))/2;
    end
    border(quants+1) = 255.000001;
    valid = isValid(H, W, border, Y);
    disp (valid);
    if (valid == 0)
        border = newBorder(H, W, border, Y);
        valid = isValid(H, W, border, Y);
        disp ('new valid = ');
        disp (valid);
    end
    Teps = 0.0001;
    DELTAeps = 0.0001;
    epsilonPred = 0;
    pmax = 20;
    newU = zeros(H, W);
    for i = 1:H
        for j = 1:W
            newU = 255;
        end
    end
    %Алгоритм Макса-Ллойда
    for P = 1:pmax
        Q = zeros (quants, 1);
        average = zeros (quants, 1);
        for i = 1:H
            for j = 1:W
                index = quants;
                for k = 1:quants
                    if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                        index = k;
                        newU(i, j) = approx(k);
                        break
                    end
                end
                Q(index) = Q(index)+1;
                average(index) = average(index)+Y(i,j);
            end
        end
        flag = 0;
        for i = 1:quants
            if average(index) == 0
                disp('ERROR! NO HIT IN QUANT!');
                flag = 1;
                break;
            end
        end
        if flag == 1
            break;
        end
        average = average ./ Q;
        approx = average;
        epsilon = 0;
        for i = 1:H
            for j = 1:W
                epsilon = epsilon + (newU(i, j) - Y(i, j))^2;
            end
        end
        epsilon = epsilon / (H * W);
        if epsilon < Teps
            break
        end
        if abs(epsilon - epsilonPred) < DELTAeps
            break
        end
        for i = 2:quants
            border(i) = (approx(i-1)+approx(i))/2;
        end
        epsilonPred = epsilon;
    end
    %%% квантование
    for i = 1:H
        for j = 1:W
            for k = 1:quants
               if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                   U(i, j, R) = approx(k);
                   break;
               end
            end
        end
    end
end

figure(1);
approx(quants+1) = approx(quants);
stairs(border, approx);

SNR = zeros (MAXR, 1);
A = Y.^2;
for R = 1:MAXR
    B = (Y-U(:, :, R)).^2;
    SNR(R) = 10*log10(sum(A)/sum(B));
end
t = 1:1:MAXR;
figure(2);
plot(t, SNR), legend('SNR(R) bright');
hold on;

PSNR = zeros (MAXR, 1);
for R = 1:MAXR
    epsilon = 0;
    for i = 1:H
        for j = 1:W
            epsilon = epsilon + (U(i, j, R) - Y(i, j))^2;
        end
    end
    PSNR(R) = 10*log10(W*H*(255^2)/epsilon);
end
figure(3);
plot(t, PSNR), legend('PSNR(R) bright');
hold on;

He = zeros (MAXR, 1);
for R = 1:MAXR
    p = zeros(2^R, 1);
    for i = 1:H
        for j = 1:W
            for k = 1:2^R
               if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                   p(k) = p(k) + 1;
                   break;
               end
            end
        end
    end
    for i = 1:2^R
        if p(i) ~= 0
            p(i) = p(i)/(H*W);
            He(R) = He(R) + (p(i)*log2(p(i)));
        end
    end
    He(R) = -He(R);
end
figure(4);
plot(He, PSNR), legend('PSNR(H) bright');
hold on;


picture = imread('C:\Users\Tuhlomon\Desktop\lenna.ppm');
H = size(picture, 1);
W = size(picture, 2);
Y = zeros(H, W);
for i = 1:H
    for j = 1:W
        Y(i, j) = picture(i, j, 1)*0.299 + picture(i, j, 2)*0.587 + picture(i, j, 3)*0.114;
    end
end

MAXR = 7;
U = zeros(H, W, MAXR);
for R = 1:MAXR
    quants = 2^R;
    minimum = min(min(Y));
    maximum = max(max(Y));
    border = zeros (quants+1, 1);
    border(1) = -0.000001;
    approx = zeros (quants, 1);
    delta = (maximum-minimum)/quants;
    for i = 1:quants
        border(i+1) = minimum + i*delta;
        approx(i) = (border(i)+border(i+1))/2;
    end
    border(quants+1) = 255.000001;
    valid = isValid(H, W, border, Y);
    disp (valid);
    if (valid == 0)
        border = newBorder(H, W, border, Y);
        valid = isValid(H, W, border, Y);
        disp ('new valid = ');
        disp (valid);
    end
    Teps = 0.0001;
    DELTAeps = 0.0001;
    epsilonPred = 0;
    pmax = 20;
    newU = zeros(H, W);
    for i = 1:H
        for j = 1:W
            newU = 255;
        end
    end
    %Алгоритм Макса-Ллойда
    for P = 1:pmax
        Q = zeros (quants, 1);
        average = zeros (quants, 1);
        for i = 1:H
            for j = 1:W
                index = quants;
                for k = 1:quants
                    if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                        index = k;
                        newU(i, j) = approx(k);
                        break
                    end
                end
                Q(index) = Q(index)+1;
                average(index) = average(index)+Y(i,j);
            end
        end
        average = average ./ Q;
        approx = average;
        epsilon = 0;
        for i = 1:H
            for j = 1:W
                epsilon = epsilon + (newU(i, j) - Y(i, j))^2;
            end
        end
        epsilon = epsilon / (H * W);
        if epsilon < Teps
            break
        end
        if abs(epsilon - epsilonPred) < DELTAeps
            break
        end
        for i = 2:quants
            border(i) = (approx(i-1)+approx(i))/2;
        end
        epsilonPred = epsilon;
    end
    %%% квантование
    for i = 1:H
        for j = 1:W
            for k = 1:quants
               if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                   U(i, j, R) = approx(k);
                   break;
               end
            end
        end
    end
end

figure(1);
approx(quants+1) = approx(quants);
stairs(border, approx);

SNR = zeros (MAXR, 1);
A = Y.^2;
for R = 1:MAXR
    B = (Y-U(:, :, R)).^2;
    SNR(R) = 10*log10(sum(A)/sum(B));
end
t = 1:1:MAXR;
figure(2);
plot(t, SNR), legend('SNR(R) bright', 'SNR(R) dark', 'SNR(R) lenna');
hold on;

PSNR = zeros (MAXR, 1);
for R = 1:MAXR
    epsilon = 0;
    for i = 1:H
        for j = 1:W
            epsilon = epsilon + (U(i, j, R) - Y(i, j))^2;
        end
    end
    PSNR(R) = 10*log10(W*H*(255^2)/epsilon);
end
figure(3);
plot(t, PSNR), legend('PSNR(R) bright', 'PSNR(R) dark', 'PSNR(R) lenna');
hold on;

He = zeros (MAXR, 1);
for R = 1:MAXR
    p = zeros(2^R, 1);
    for i = 1:H
        for j = 1:W
            for k = 1:2^R
               if Y(i, j) >= border(k) && Y(i, j) < border(k+1)
                   p(k) = p(k) + 1;
                   break;
               end
            end
        end
    end
    for i = 1:2^R
        if p(i) ~= 0
            p(i) = p(i)/(H*W);
            He(R) = He(R) + (p(i)*log2(p(i)));
        end
    end
    He(R) = -He(R);
end
figure(4);
plot(He, PSNR), legend('PSNR(H) bright', 'PSNR(H) dark', 'PSNR(H) lenna');
hold on;