page_screen_output(false);



# Figura 3.3
close all;
N=256;
t = linspace(0,1,256);
f = (t-t.^2).^2;
X = fft(f);
figure(1);
plot(t,f);
figure(2);
plot(-N/2:N/2-1,abs(fftshift(X)).^2/N);



# Tabela 3.1
M = max(abs(X));
# C representa o conjunto de limiares
# para eliminação de coeficientes
# "pequenos" do espectro
C = [0.5 0.1 0.03 0.02 0.01 0.005];
L=length(C);
P=[];
D=[];
for c=C
   Y=X.*(abs(X)>c*M);
   y = real(ifft(Y));
   P=[P sum(abs(X)>c*M)/N];
   D=[D norm(f-y)^2/norm(f)^2];
endfor
printf("c\tP(c)\t\tD(c)\n");
printf("%g\t%g\t%g\n",[C;P;D]);
# Recria tabela em forma de gráfico
C = 10.^linspace(-0.3,-8,50);
L=length(C);
P=[];
D=[];
for c=C
   Y=X.*(abs(X)>c*M);
   y = real(ifft(Y));
   P=[P sum(abs(X)>c*M)/N];
   D=[D norm(f-y)^2/norm(f)^2];
endfor
figure(3);
plot(1:L,C,";c;",1:L,P,";P(c);",1:L,D,";D(c);");



# Figura 3.4
close all;
N=256;
t = linspace(0,1,256);
f = zeros(1,N);
f(1:0.2*N+1) += 1;
X = fft(f);
figure(1);
plot(t,f);
figure(2);
plot(-N/2:N/2-1,abs(fftshift(X)).^2/N);



# Tabela 3.2
M = max(abs(X));
C = [0.5 0.1 0.03 0.02 0.01 0.005];
L=length(C);
P=[];
D=[];
for c=C
   Y=X.*(abs(X)>c*M);
   y = real(ifft(Y));
   P=[P sum(abs(X)>c*M)/N];
   D=[D norm(f-y)^2/norm(f)^2];
endfor
printf("c\tP(c)\t\tD(c)\n");
printf("%g\t%g\t%g\n",[C;P;D]);
# Recria tabela em forma de gráfico
C = 10.^linspace(-0.3,-8,50);
L=length(C);
P=[];
D=[];
for c=C
   Y=X.*(abs(X)>c*M);
   y = real(ifft(Y));
   P=[P sum(abs(X)>c*M)/N];
   D=[D norm(f-y)^2/norm(f)^2];
endfor
figure(4);
plot(1:L,C,";c;",1:L,P,";P(c);",1:L,D,";D(c);");



# Figura 3.5
# Reconstrução da função usando compressão de 90%
# o valor do limiar abaixo foi escolhido "a olho"
# para produzir o valor de P(c) aproximadamente
# 0.1 (10% do tamanho do arquivo original)
c=7.5e-2;
Y=X.*(abs(X)>c*M);
y = real(ifft(Y));
printf("P(%g) = %g\n",c,sum(abs(X)>c*M)/N);
printf("D(%g) = %g\n",c,norm(f-y)^2/norm(f)^2);
figure(5);
plot(t,f,t,y);



# Exemplo de compressão 3
close all;
N=256;
t = linspace(0,1,256);
f = t;
X = fft(f);
figure(1);
plot(t,f);
figure(2);
plot(-N/2:N/2-1,abs(fftshift(X)).^2/N);



# Tabela 3.3
M = max(abs(X));
C = [0.5 0.1 0.03 0.02 0.01 0.005];
L=length(C);
P=[];
D=[];
for c=C
   Y=X.*(abs(X)>c*M);
   y = real(ifft(Y));
   P=[P sum(abs(X)>c*M)/N];
   D=[D norm(f-y)^2/norm(f)^2];
endfor
printf("c\tP(c)\t\tD(c)\n");
printf("%g\t%g\t%g\n",[C;P;D]);
# Recria tabela na forma de gráfico
C = 10.^linspace(-0.3,-8,50);
L=length(C);
P=[];
D=[];
for c=C
   Y=X.*(abs(X)>c*M);
   y = real(ifft(Y));
   P=[P sum(abs(X)>c*M)/N];
   D=[D norm(f-y)^2/norm(f)^2];
endfor
# Guarda os vetores P e D para comparar
# com DCT (Figura 3.10)
PFFT=P;DFFT=D;
figure(6);
plot(1:L,C,";c;",1:L,P,";P(c);",1:L,D,";D(c);");



# Figura 3.6
# Reconstrução da função usando compressão de 90%
# o valor do limiar abaixo foi escolhido
# "a olho" para produzir o valor de P(c)
# aproximadamente 0.1 (10% do tamanho do
# arquivo original)
c=2.5e-2;
Y=X.*(abs(X)>c*M);
y = real(ifft(Y));
printf("P(%g) = %g\n",c,sum(abs(X)>c*M)/N);
printf("D(%g) = %g\n",c,norm(f-y)^2/norm(f)^2);
figure(7);
plot(t,f,t,y);



# Figura 3.7
# Usa o espectro original para sintetizar
# 3 períodos da forma de onda:
# basta criar um vetor 3 vezes maior
# e espaçar os coeficientes do
# espectro original de 3 em 3, o que
# corresponde a multiplicar a
# frequência fundamental por 3.
Y3=zeros(1,3*N);
Y3(1:3:3*N)=Y;
y3=real(ifft(Y3));
figure(8);
plot(linspace(-1,2,3*N),y3);



# Motivação para DCT:
# espelhamento da função original para eliminar
# a descontinuidade na fronteira da janela.
close all;
tt = linspace(0,2,2*N);
ff = [f f(N:-1:1)];
XX = fft(ff);
M=max(abs(XX));
# Reconstrução usando o mesmo limiar c anterior,
# porém com muito mais compressão e muito mais
# qualidade (menor distorção)
c=2.5e-2;
YY=XX.*(abs(XX)>c*M);
yy = real(ifft(YY));
printf("P(%g) = %g\n",c,sum(abs(XX)>c*M)/(2*N));
printf("D(%g) = %g\n",c,norm(ff-yy)^2/norm(ff)^2);
figure(9);
plot(tt,ff,tt,yy);
# Reconstrução usando um limiar c 100x menor
c=2.5e-4;
YY=XX.*(abs(XX)>c*M);
yy = real(ifft(YY));
printf("P(%g) = %g\n",c,sum(abs(XX)>c*M)/(2*N));
printf("D(%g) = %g\n",c,norm(ff-yy)^2/norm(ff)^2);
figure(10);
plot(tt,ff,tt,yy);



# Figura 3.8
# Funções básicas usadas na DCT
close all;
N=256;
t = 0:N-1;
hold on;
for k=0:10
   rgbcolor = int8([bitget(k+1,1) bitget(k+1,2) bitget(k+1,3)]);
   plot(t,sqrt(2/N)*cos(pi*k*t/N),
        "color",rgbcolor);
   drawnow;
   pause(1);
endfor
hold off;



# Repetição dos exemplos anteriores de
# compressão usando DCT



# Função quádrica f(t) = (t-t^2)^2
close all;
N=256;
t = linspace(0,1,256);
f = (t-t.^2).^2;
X = dct(f);
figure(11);
plot(t,f);
figure(12);
plot(0:N-1,abs(X).^2/N);



# Tabela 3.4
M = max(abs(X));
C = 10.^linspace(-0.3,-8,50);
L=length(C);
P=[];
D=[];
for c=C
   Y=X.*(abs(X)>c*M);
   y = real(idct(Y));
   P=[P sum(abs(X)>c*M)/N];
   D=[D norm(f-y)^2/norm(f)^2];
endfor
figure(13);
plot(1:L,C,";c;",1:L,P,";P(c);",1:L,D,";D(c);");



# Função "degrau"
close all;
N=256;
t = linspace(0,1,256);
f = zeros(1,N);
f(1:0.2*N+1) += 1;
X = dct(f);
figure(11);
plot(t,f);
figure(12);
plot(0:N-1,abs(X).^2/N);



# Tabela 3.5
M = max(abs(X));
C = 10.^linspace(-0.3,-8,50);
L=length(C);
P=[];
D=[];
for c=C
   Y=X.*(abs(X)>c*M);
   y = real(idct(Y));
   P=[P sum(abs(X)>c*M)/N];
   D=[D norm(f-y)^2/norm(f)^2];
endfor
figure(14);
plot(1:L,C,";c;",1:L,P,";P(c);",1:L,D,";D(c);");



# Reconstrução da função degrau usando DCT
# com compressão de 90%
# o valor do limiar abaixo foi escolhido
# "a olho" para produzir o valor de P(c)
# aproximadamente 0.1 (10% do tamanho do
# arquivo original)
c=4.1e-2;
Y=X.*(abs(X)>c*M);
y = real(idct(Y));
printf("P(%g) = %g\n",c,sum(abs(X)>c*M)/N);
printf("D(%g) = %g\n",c,norm(f-y)^2/norm(f)^2);
figure(15);
plot(t,f,t,y);



# Função f(t) = t
close all;
N=256;
t = linspace(0,1,256);
f = t;
X = dct(f);
figure(11);
plot(t,f);
figure(12);
plot(0:N-1,abs(fftshift(X)).^2/N);



# Tabela 3.6
M = max(abs(X));
C = 10.^linspace(-0.3,-8,50);
L=length(C);
P=[];
D=[];
for c=C
   Y=X.*(abs(X)>c*M);
   y = real(idct(Y));
   P=[P sum(abs(X)>c*M)/N];
   D=[D norm(f-y)^2/norm(f)^2];
endfor
PDCT=P;DDCT=D;
figure(16);
plot(1:L,C,";c;",1:L,P,";P(c);",1:L,D,";D(c);");



# Reconstrução da função f(t)=t usando DCT
# com compressão de 90%
# o valor do limiar abaixo foi escolhido
# "a olho" para produzir o valor de P(c)
# aproximadamente 0.1 (10% do tamanho do
# arquivo original)
c=2e-4;
Y=X.*(abs(X)>c*M);
y = real(idct(Y));
printf("P(%g) = %g\n",c,sum(abs(X)>c*M)/N);
printf("D(%g) = %g\n",c,norm(f-y)^2/norm(f)^2);
figure(17);
plot(t,f,t,y);



# Figura 3.10
# Recria gráfico da FFT
X = fft(f);
M = max(abs(X));
C = 10.^linspace(-0.3,-8,50);
L=length(C);
P=[];
D=[];
for c=C
   Y=X.*(abs(X)>c*M);
   y = real(ifft(Y));
   P=[P sum(abs(X)>c*M)/N];
   D=[D norm(f-y)^2/norm(f)^2];
endfor
PFFT=P;DFFT=D;
figure(18);
plot(PFFT(1:14),log10(DFFT(1:14)),";FFT;",
     PDCT(1:42),log10(DDCT(1:42)),";DCT;");
title("Compressao usando FFT e DCT");
xlabel("P(c)");
ylabel("log10(D(c))");



# Figura 3.11
close all;
M=rgb2gray(imread("double_ferris.jpg"));
S=size(M);
MM=[M M(1:S(1),S(2):-1:1); M(S(1):-1:1,1:S(2)) M(S(1):-1:1,S(2):-1:1)];
figure(1);
imshow(MM);



# Figura 3.12
# Recorta a figura um pouco para obter divisibilidade por unidades de 8x8 pixels
M=M(1:S(1)-mod(S(1),8),1:S(2)-mod(S(2),8),:);
S=size(M);
figure(1);
imshow(M);
N=dct2(M);
figure(2);
Nlog = log(1+abs(N));
maximo = max(max(Nlog));
imshow(Nlog/maximo);
MBAK = M;
NBAK=N;
for j=1:8:S(1)-1
   for k=1:8:S(2)-1
      M=MBAK(j:j+7,k:k+7);
      N = dct2(M);
      NBAK(j:j+7,k:k+7) = N;
   endfor
endfor
figure(3);
imshow(log(1+abs(NBAK))/maximo);
# detalhe da roda gigante superior
#figure(4);
#imshow(log(1+abs(NBAK(1:400,300:600)))/maximo);



# Figura 3.13
# Reordena a DCT em blocos para juntar os coeficientes correspondentes a uma mesma frequência (j,k)
NOVO = NBAK;
for j=0:7
   for k=0:7
      NOVO(j*S(1)/8+1:(j+1)*S(1)/8,k*S(2)/8+1:(k+1)*S(2)/8)  = NBAK(j+1:8:S(1),k+1:8:S(2));
   endfor
endfor
figure(4);
imshow(log(1+abs(NOVO))/maximo);
# Detalhe da parte correspondente aos coeficientes DC de cada bloco
figure(5);
imshow(log(1+abs(NOVO(1:S(1)/8,1:S(2)/8,:)))/maximo);



# Matriz de quantização usada em cada bloco da DCT
# durante a compressão JPG
close all;
E = [  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   ] ;
figure(1);
imshow(E/max(max(E)));



# Exemplo 3.1, figuras 3.14 e 3.15
B = [  47  32  75 148 192 215 216 207 ;
       36  82 161 196 205 207 190 140 ;
       86 154 200 203 213 181 143  82 ;
      154 202 209 203 159 145 147 127 ;
      184 207 199 147 134 127 137 138 ;
      205 203 125  72 123 129 150 115 ;
      209 167 126 107 111  94 105 107 ;
      191 129 126 136 106  54  99 165   ] ;
figure(1);
maxB = max(max(B));
imshow(B/maxB);
# Calcula a DCT da imagem com valores centralizados (entre -127 e +128)
Bhat = dct2(B-127);
# Os valores de R abaixo correspondem ao fator de escala associado
# à matriz de quantização. Quanto maior este fator, maior a distância
# entre os níveis de quantização, e consequentemente serão maiores
# tanto a compressão quanto a perda (ou ruído de quantização).
R = [ 1 0.1 0.5 2 ];
for k=1:4
  # Calcula valores quantizados
  Bquant = round(Bhat./(R(k)*E));
  # Retorna à escala original (este é o passo da de(s)quantização)
  Btil = Bquant.*(R(k)*E);
  # Aplica a DCT, voltando os valores à faixa 0...255
  Bn = round(idct2(Btil)+127);
  figure(k+1);
  imshow(Bn/maxB);
  # Calcula a distorção da imagem codificada em relação à original
  printf("Distorção(r=%g) = %g\n",R(k),norm(B-Bn)/norm(B));
  # Estimador da compressão a partir do número de zeros da DCT
  printf("Compressão(r=%g) = %g%%\n",R(k),100*sum(sum(Bquant==0))/64);
endfor



# Figura 3.16, usando a imagem da roda gigante
close all;
M=rgb2gray(imread("double_ferris.jpg"));
S=size(M);
# Recorta a imagem para que largura e altura sejam múltiplos de 8
M=M(1:S(1)-mod(S(1),8),1:S(2)-mod(S(2),8),:);
S=size(M);
figure(1);
imshow(M);
# Calcula a DCT em blocos de 8x8, e mostra o resultado
N=blkdct2(M,8,8);
Nlog=log(1+abs(N));
maxN = max(max(Nlog));
figure(2);
imshow(Nlog/maxN);
# Codifica cada bloco 8x8 usando o mesmo esquema de
# quantização/dequantização da DCT do exemplo anterior.
Nquant = N;
Ntil = N;
# Parâmetro de granularidade da quantização
# (faz pouca diferença nesse exemplo; experimente R=0.1)
R=1;
# Processa todos os blocos
for j=1:8:S(1)
  for k=1:8:S(2)
    Nquant(j:j+7,k:k+7) = round(N(j:j+7,k:k+7)./(R*E));
    Ntil(j:j+7,k:k+7) = (R*E).*Nquant(j:j+7,k:k+7);
  endfor;
endfor;
# Faz a iDCT em blocos, restaurando a imagem original
# (uint8 é apenas uma conversão de tipo para inteiros de 8 bits
#  sem sinal, para compatibilidade com a função imshow)
Mn = uint8(round(blkidct2(Ntil,8,8)));
figure(3);
imshow(Mn);
# Mostra imagem "diferença"
figure(4);
imshow(abs(M-Mn));
# Mostra imagem "diferença" normalizada
figure(5);
imshow(256*(abs(M-Mn)/14));
# Mede níveis de distorção e compressão (estimada)
printf("Distorção(r=%g) = %g\n",R,norm(double(M-Mn))/norm(double(M)));
printf("Compressão(r=%g) = %g%%\n",R,100*sum(sum(Nquant==0))/(S(1)*S(2)));


# Exemplos de reconstrução progressiva

# Reconstrução progressiva usando sub-blocos de lxl (l=1...8)
# (este esquema não aparece no livro)
close all;
figure(2);
for l=1:8
  Nl=[]; Ml=[];
  for j=1:8:S(1)
    for k=1:8:S(2)
      # encontra a posição do início (j,k)-ésimo bloco
      # dentro da DCT reconstruída com sub-blocos lxl
      fj = floor(j/8)*l; fk = floor(k/8)*l;
      # copia o bloco lxl, ajustando a escala da DCT
      Nl(fj+1:fj+l,fk+1:fk+l) = Ntil(j:j+l-1,k:k+l-1)*(sqrt(l*l)/sqrt(64));
      # Obs: o fator (sqrt(l*l)/sqrt(64)) serve para compensar as constantes
      # diferentes na DCT de tamanho 8x8 (dos dados originais) em relação às
      # DCTs de lxl que serão usadas na reconstrução abaixo.
    endfor;
  endfor;
  # Reconstrói a imagem pela iDCT com blocos lxl
  Ml = uint8(round(blkidct2(Nl,l,l)));
  imshow(Ml);
  title(sprintf("Reconstrução com sub-blocos de %dx%d, distorção estimada=%g",l,l,norm(double(M(floor(1:(8/l):S(1)),floor(1:(8/l):S(2)))-Ml))/norm(double(M(floor(1:(8/l):S(1)),floor(1:(8/l):S(2)))))));
  drawnow;
  pause(5);
endfor



# Exemplo 3.2 e figura 3.17
# Reconstrução progressiva usando coeficientes DCT até j+k=l (para l=0...14)
close all;
Nl=0*Ntil;
figure(2);
for l = 0:14
  for j=0:l
    k = l-j;
    Nl(j+1:8:S(1),k+1:8:S(2)) = Ntil(j+1:8:S(1),k+1:8:S(2));
  endfor
  Ml = uint8(round(blkidct2(Nl,8,8)));
  imshow(Ml);
  title(sprintf("Reconstrução progressiva até j+k=%d, distorção estimada=%g",l,norm(double(M-Ml))/norm(double(M))));
  drawnow;
  pause(8/3); # calculado para durar +/- o mesmo que a outra animação
endfor