clear all
clc
tic

B = 4;                                % velikost bloku pri preskladavani pomoci chaoticke mapy
Blok = 8;                             % vel. bloku pro zpracovani  
Qfaktor = 20;
kolikrat_rotovat = 9;
stred = 1;                            % str. hodnota v bloku 8x8 a vyber hodnot podle ni  
alfa = 0.2;

obraz=imread('c:\metoda_2\patchwork1.jpg');             % 256 x 256
[m,n]=size(obraz);

pocet=(m*n)/(B^2)                     % kolik musime udelat bloku z obrazu kvuli chaot mape
N = pocet;                            % delka vektoru chaoticke mapy

vodoznak=im2double(imread('c:\metoda_2\obrazky\copyright.bmp')); 
[yv,xv]=size(vodoznak);


% definice velikosti vektoru pro chaotickou mapu
s = zeros(1,N);

% % pocatecni podminka meli oni nechame stejnou
% s(1) = 0.1664;
% 
% % chaoticke chovani co nejvetsi rozkmit hodnot v intervalu <0,1>
% R = 4;
% 
% % vypocet podle vzorce
% for cnt= 2:pocet
%   s(cnt) = R*s(cnt-1)*(1-s(cnt-1));
% end;
% 
% s = ceil(s.*pocet);              % zisk opakujici se posloupnosti zaokr. k nekonecnu
% 
% % figure(255)
% % plot(s);
% 
% % ******************************************************************
% % pokud se vyskytuje nektere cislo vicekrat dame 0
% 
% cnt=1;
% 
% for i=1:N
%     nalez = find(s==i);                     % fce pro nalezeni prvku ve vektoru
%     if isempty(nalez)                        % pokud nebyl nalezen prvek
%         chybejici(cnt)=i;                    % ulozim cislo ktere chybi
%         cnt=cnt+1;
%     else
%         for j=2:size(nalez,2)                % budeme pocitat od druheho prvku
%             s(nalez(j)) = 0;                 % dame nuly na mista ktera se opakuji
%         end
%     end     
% end
% 
% % ******************************************************************
% % nahrada opakujicich se prvku neopakujicimi
% 
% nula = find(s==0);                              % fce pro nalezeni 0 prvků ve vektoru
% 
% for j=1:size(nula,2)                            % budeme pocitat od prvniho prvku
%      s(nula(j)) = chybejici(j);                 % chybej. hodnoty na mista ktera byla 0
% end
% 
% % figure(251)
% % plot(s);

% s = randperm(pocet);                            % bude moznost volby
% 
% save s                                          % ulozeni chaoticke mapy        
% 
% % ******************************************************************
% 
% figure(1)
% colormap('Gray')
% imagesc(obraz)
% title('Puvodni obraz')
% 
% obraz_new = zeros(size(obraz));
% vysledek = zeros(size(obraz));
% 
% mat = struct('obr',{zeros(B)},'label',{});
% mat_help = struct('obr',{zeros(B)},'label',{});
% chaos = struct('obr',{zeros(B)},'label',{});
% chaos_help = struct('obr',{zeros(B)},'label',{});
% 
% % ***********************************************************************
% % pro bloky NxN, rozdeleni do bloku s labelem
% 
% cnt = 0;
% x_krat = m/B;
% y_krat = n/B;
% 
%     for j=0:(x_krat-1)
%         for k=0:(y_krat-1)
%              cnt = cnt+1;
%              mat{cnt}.obr = obraz(j*B+1:j*B+B,k*B+1:k*B+B);
%              mat{cnt}.label = cnt;
%         end 
%     end
% 
% % *******************************************************************
% % preskladavani podle chaoticke mapy
% 
% for i=1:kolikrat_rotovat
%     cnt = 1;  
%         for j=0:(x_krat-1)
%             for k=0:(y_krat-1)
%                 index = find(s==cnt);
%                 chaos{cnt}.obr = mat{index}.obr;
%                 chaos{cnt}.label = index;                  % kde se nachazelo cislo pro prehozeni
%                 cnt = cnt+1;                               % vyhledavani cisel 0.1.2.3...
%             end
%         end
%          mat = chaos;
% end
% 
% cnt = 0;
%     for j=0:(x_krat-1)
%         for k=0:(y_krat-1)
%              cnt = cnt+1;
%              obraz_new(j*B+1:j*B+B,k*B+1:k*B+B) = mat{cnt}.obr;
%         end 
%     end

obraz_new = chaotic(obraz,kolikrat_rotovat,0);
    
% figure(2)
% colormap('Gray')
% imagesc(obraz_new)
% title(['Preskladany obraz - ', num2str(kolikrat_rotovat), ' krat']);

% *************************************************************************
% IBPC - registr a kontejner, vytvoreni, vstup je obraz_new, typ ABAB ....

RR = zeros(m,n/2);
RC = zeros(m,n/2);

for j=1:m
    cnt = 0;
    for k=1:n/2
            RC(j,k) = obraz_new(j,k+cnt);
            RR(j,k) = obraz_new(j,k+cnt+1);
            cnt = cnt + 1;
    end 
end

% figure(3)
% colormap('Gray')
% imagesc(RR)                 % RR = reference registr a RC = container 
% title('Referencni registr');

% *************************************************************************
% vytvoreni transformacni matice o velikosti B (8*8)

T=dctmtx(Blok);

%transformace do dct po blocich 8*8
RC_dct = blkproc(RC,[Blok Blok],'P1*x*P2',T,T');
RR_dct = blkproc(RR,[Blok Blok],'P1*x*P2',T,T');

% *************************************************************************
% kvantizace Q - faktor >50 zmenseni koef. <50 zvetseni 

QTable =   [16   11   10   16   24   40    51    61         
            12   12   14   19   26   58    60    55
            14   13   16   24   40   57    69    56
            14   17   22   29   51   87    80    62
            18   22   37   56   68   109   103   77
            24   35   55   64   81   104   113   92
            49   64   78   87   103  121   120   101
            72   92   95   98   112  100   103   99];   % kvantizacni jpeg matice

% *************************************************************************
% vypocet Qtable podle Q faktoru
       
if Qfaktor <= 50
    QTable = QTable.*50/Qfaktor;
else
    QTable = QTable.*(2-0.02*Qfaktor);
end

fun = @(x) x./QTable;                        % x=RR_dct, ktere pak vola blkproc

QA = blkproc(RC_dct,[Blok Blok],fun);
QB = blkproc(RR_dct,[Blok Blok],fun);        % samotna kvantizace

% *************************************************************************
% vkladani vodoznaku 
% je mozne vodoznak svou velikosti vlozit do obrazu?

[yqb,xqb]=size(QB);

save('QA.mat');
save('QB.mat');

suma=0;
jednicky=1;
nuly=1;

for j=0:(yqb/Blok-1)
    for k=0:(xqb/Blok-1)

        matQB = QB(j*Blok+1:j*Blok+Blok,k*Blok+1:k*Blok+Blok);
        matQA = QA(j*Blok+1:j*Blok+Blok,k*Blok+1:k*Blok+Blok);

        matQB(1,1) = 0;                                     % snulovani stredni hodnoty
        
        plus = find((matQB) > stred);                       % {+} > nez stred, indexy
        minus = find((matQB) < -stred);                     % {-} < nez stred, indexy
                     
        for y=1:size(plus,1)
            suma = suma + 1;
              if suma <(yv*xv)+1                                              % delej jen do rozmeru vodoznaku
                 if vodoznak(suma) == 1
                        if matQA(plus(y)) > 0
                           matQA(plus(y)) = matQA(plus(y)) + abs(alfa*matQB(plus(y)));        % v jejich navrhu je to matice E
                           jednicky = jednicky+1;                                                  % pocet jednicek
                        else
                           matQA(plus(y)) = matQA(plus(y)) - abs(alfa*matQB(plus(y)));        % v jejich navrhu je to matice E
                           jednicky = jednicky+1;                                                  % pocet jednicek
                        end                           
                 else  
                        nuly = nuly + 1;
                  end
             end
        end
            
        for y=1:size(minus,1)
            suma = suma + 1;
              if suma <(yv*xv)+1                                                                   % delej jen do rozmeru vodoznaku
                 if vodoznak(suma) == 1
                     if matQA(minus(y)) >= 0
                        matQA(minus(y)) = matQA(minus(y)) + abs(alfa*matQB(minus(y)));        % v jejich navrhu je to matice E
                        jednicky = jednicky + 1;                                                     % pocet jednicek
                     else
                        matQA(minus(y)) = matQA(minus(y)) - abs(alfa*matQB(minus(y)));             % v jejich navrhu je to matice E
                        jednicky = jednicky + 1;                                                     % pocet jednicek                       
                     end
                 else  
                        nuly = nuly + 1;
                        
                 end               
               end
        end
        QA(j*Blok+1:j*Blok+Blok,k*Blok+1:k*Blok+Blok) = matQA;
    end   
end
         
jednicky
nuly
suma

%  xlswrite('qa_uprav', QA);
%  xlswrite('qb_uprav', QB);

% *************************************************************************
% po vkladani vodoznaku nutno QA vynasobit kvantizacni matici


fun = @(x) x.*QTable;                           % x=RR_dct, ktere pak vola blkproc
MA = blkproc(QA,[Blok Blok],fun);

RC_IDCT = blkproc(MA,[Blok Blok],'P1*x*P2',T',T);     % zpetna dct transformace po blocich B*B

RC_IDCT = ceil(RC_IDCT);

% *************************************************************************
% IIBPC  - inverze, registr a kontejner ABABAB

for j=1:m
    cnt = 0;
    for k=1:n/2
            obraz_IIBPC(j,k+cnt) = RC_IDCT(j,k);
            obraz_IIBPC(j,k+cnt+1) = RR(j,k); 
            cnt = cnt + 1;
    end 
end

% figure(4)
% title('Zpetne prehozeni po vlozeni vodoznaku');
% colormap('Gray')
% imagesc(obraz_IIBPC)                 % reference registr A

% *******************************************************************
% preskladavani zpet podle chaoticke mapy


vysledek = chaotic(obraz_IIBPC,kolikrat_rotovat,1);

% cnt = 0;
%     for j=0:(x_krat-1)
%         for k=0:(y_krat-1)
%              cnt = cnt+1;
%              chaos_help{cnt}.obr = obraz_IIBPC(j*B+1:j*B+B,k*B+1:k*B+B);
%         end 
%     end
% 
% for i=1:kolikrat_rotovat
%     cnt = 1;  
%         for j=0:(x_krat-1)
%             for k=0:(y_krat-1)
%                 index = chaos{cnt}.label;
%                 mat_help{index}.obr = chaos_help{cnt}.obr;
%                 cnt = cnt+1;                               % vyhledavani cisel 0.1.2.3...
%             end
%         end
%            chaos_help = mat_help;
% end   
% 
% cnt = 0;
%     for j=0:(x_krat-1)
%         for k=0:(y_krat-1)
%              cnt = cnt+1;
%              vysledek(j*B+1:j*B+B,k*B+1:k*B+B) = chaos_help{cnt}.obr;
%         end 
%     end
%        

kapacita = suma / (m*n)

figure(5)
colormap('Gray')
imagesc(vysledek)                                % opetovne slozeni obrazu
title('Zpetne slozen puvodni obraz + vodoznak');

% vypocet PSNR
obr_original = double(obraz);
obr_vodoznaceny = double(uint8(vysledek));
obr_velikost = double(m*n);
MSE=(sum(sum((obr_original-obr_vodoznaceny).^2)))/obr_velikost;
PSNR=10*log10(255^2/(MSE))

% ulozeni vodoznaceneho obrazku na zvolene misto
imwrite(uint8(vysledek),'PatchworkWM1.bmp');                        % ulozeni kvuli checkmarku
imwrite(uint8(vysledek),'PatchworkWM1.jpg','Quality',100,'Bitdepth',8);
toc