Instructions

%% Aluno: Thiago Carvalho Miranda
%%Basic Matlab/Scilab Instructions

% O primeiro tutorial mostra apenas a introdução a comandos do
% Matlab/Scilab para o processamento de sinais.

a = 1; a = 2+1i; % real and complex numbers
b = [1 2 3 4]; % row vector
c = [1; 2; 3; 4]; % column vector
d = 1:2:7; % here one has d=[1 3 5 7]

size(d) %tamanho de d
d(1) %mostra a primeira entrada em d
d(1:2) %mostra o sub-array com os dois primeiros valores de d
% Os códigos abaixo servem para criar arrays pre definidos, neste caso
% utilizando conceitos de identidade, 1's e numeros aleatorios
A = eye(4,4); 
B = ones(4,4);
C = rand(4,4);
C = b'; %C é a transposta de b

D = B*A %Multiplicação de matrizes
D = B.*A %Multiplica cada entrada do vetor

%Aplicando funções em matrizes diretamente
E = A./B; % division
E = sin(A); % sinus is applied to each entry
E = abs(A + 1i*B); % modulus of each entry

%Diferentes métodos de modificar matrizes e arrays
b = sort(b); % sort values
b = b .* (b>2); % set to zeros (threshold) the values below 2
b(3) = []; % suppress the 3rd entry of a vector
B = [b; b]; % create a matrix of size 2x4
c = B(:,2); % to access 2nd column

b(end-2:end) = 1; % to access the last entries
b = b(end:-1:1); % reverse a vector

%Instruções avançadas
disp('Hello'); % display a text
x = 1.23456;
disp( sprintf('Value of x=%.2f', x) ); % print a values with 2 digits
A(A==Inf) = 3; % replace Inf values by 3
A(:); % flatten a matrix into a column vector
max(A(:)); % max of a matrix

%Transforma em 0 valores menores que 3.
C = C .* (abs(C)>.3);

%Construções básicas de programas usando for e while
for i=1:3 % repeat the loop for i=1, i=2, i=3
    disp(i); % make here something
end
i = 3;
while i>0 % while syntax
    disp(i); % do smth
    i = i-1;
end

%Introdução a Processamento de Imagens
n = 256; % size of the image
M = load_image('lena', n);
clf;
imageplot(M);
%Manipulação da imagem como um array arbitrário
clf;
imageplot(M(1:50,1:50), 'Zoom', 1,2,1);
imageplot(-M, 'Reversed contrast', 1,2,2);


%% Introduction to Signal Processing


% Loading and displaying signals
n = 512;
f = load_signal('Piece-Regular', n);
f = f(:); %Making sure it's a vector
f = rescale(f); % Rescaling the entries of the signal

% Displaying the signal - basic style
clf;
plot(1:n, f);
axis('tight');
title('My title');
set_label('variable x', 'variable y');

% Different ways to plot using |subplot| 
clf;
subplot(2, 2, 1);
plot(f); axis('tight');
clf;
subplot(2, 2, 4);
plot(f.^2); axis('tight');

% Displaying two signals at once
clf;
plot(1:n, [f f.^2]');
legend('signal', 'signal^2');

%% Introduction to Image Processing
% Carregando a imagem
name = 'lena';
n = 256;
M = load_image(name, []);
M = rescale(crop(M,n));
%Mostrando a imagem na tela, selecionando uma área de interesse para dar
%zoom
figure;imageplot(M, 'Original', 1,2,1);
imageplot(crop(M,50), 'Zoom', 1,2,2);
% As imagens podem ser modificadas como arrays arbitrários. Imagens em tons
% de cinza são arrays 2D e podem ser modificadas como tal.
figure;imageplot(-M, '-M', 1,2,1);
imageplot(M(n:-1:1,:), 'Flipped', 1,2,2);
% Exemplo de modificação de imagem chamado borramento, que consiste em
% convoluir a imagem com uma máscara
% Construindo a máscara
k = 9;
h = ones(k,k);
h = h/sum(h(:));
% Convoluindo a imagem
Mh = perform_convolution(M, h);
% Mostrando a imagem
figure;imageplot(M, 'Image', 1,2,1);
imageplot(Mh, 'Blurred', 1,2,2);
% Existem diversas implementações de convolução e diferenciação
G = grad(M);
figure;imageplot(G(:,:,1), 'd/dx', 1,2,1);
imageplot(G(:,:,2), 'd/dy', 1,2,2);
% A transformada de fourier serve para interpolação da imagem e aproximação
% em baixa frequência
Mf = fft2(M);
Lf = fftshift(log( abs(Mf)+1e-1 ));
figure;imageplot(M, 'Image', 1,2,1);
imageplot(Lf, 'Fourier transform', 1,2,2);
% Interpolação de imagens usando alta frequência
p = 64;
n = p*4;
M = load_image('boat', 2*p); M = crop(M,p);
Mf = fftshift(fft2(M));
MF = zeros(n,n);
sel = n/2-p/2+1:n/2+p/2;
sel = sel;
MF(sel, sel) = Mf;
MF = fftshift(MF);
Mpad = real(ifft2(MF));
figure;imageplot( crop(M), 'Image', 1,2,1);
imageplot( crop(Mpad), 'Interpolated', 1,2,2);
% interpolação de imagem usando cubic-splines.
Mspline = image_resize(M,n,n);
figure;imageplot( crop(Mpad), 'Fourier (sinc)', 1,2,1);
imageplot( crop(Mspline), 'Spline', 1,2,2);

%% Image Approximation with Fourier and Wavelets 

% Loading and displaying the image
name = 'lena';
n0 = 512;
f = rescale( load_image(name,n0) );
clf;
imageplot( f, 'Image f');

% Zoom
clf;
imageplot( crop(f,64), 'Zoom' );

% Flip image
clf;
imageplot(-f, '-f', 1,2,1);
imageplot(f(n0:-1:1,:), 'Flipped', 1,2,2);

% Blurrying image
k = 9; % size of the kernel
h = ones(k,k);
h = h/sum(h(:)); % normalize
fh = perform_convolution(f,h);
clf;
imageplot(fh, 'Blurred image');

% Fourier transform
F = fft2(f) / n0;
L = fftshift(log( abs(F)+1e-1 ));
clf;
imageplot(L, 'Log(Fourier transform)');

% Wavelet transform

Jmin = 0;
fw = perform_wavelet_transf(f,Jmin,+1);
clf;
plot_wavelet(fw);

%% Image Processing with Wavelets 
% Podemos fazer aproximação de imagens utilizando coeficientes wavelet
% Carregando e mostrando a imagem na tela
name = 'cortex';
n0 = 512;
f = load_image(name,n0);
f = rescale( sum(f,3) );
figure;imageplot(f);
% Setando os parâmetros para a aproximação via wavelet.
Jmin = 0;
% Utilizando uma função pronta que aplica a transformação wavelet
% utilizando a imagem e a menor escala para a transformação.
Psi = @(f)perform_wavelet_transf(f,Jmin,+1);
% Função pronta para a transformação inversa
PsiS = @(fw)perform_wavelet_transf(fw,Jmin,-1);
% Aplicando a transformação e mostrando a imagem na tela
fW = Psi(f);
figure; plot_wavelet(fW);
% Setando novos parâmetros para a transformada
T = .5;
Thresh = @(fW,T)fW .* (abs(fW)>T); % Utilizando uma função para o operador de limite
% Setando os limites do coeficiente
fWT = Thresh(fW,T);
% Mostrando os coeficientes de limite
figure;subplot(1,2,1);
plot_wavelet(fW);
title('Original coefficients');
subplot(1,2,2);
plot_wavelet(fWT);
%Reconstruindo a imagem usando a transformada inversa e mostrando a
%aproximação obtida
f1 = PsiS(fWT);
figure;imageplot(f, 'Image', 1,2,1);
imageplot(clamp(f1), strcat(['Approximation, SNR=' num2str(snr(f,f1),3) 'dB']), 1,2,2);
% Novamente setando os limites
fWT = Thresh(fW,T);
% Checando o número de coeficientes diferentes de 0 em |fWT|
disp(strcat(['      M=' num2str(M)]));
disp(strcat(['|fWT|_0=' num2str(sum(fWT(:)~=0))]));

% Remoção de ruído usando coeficientes wavelet
% Setando um parâmetro para adicionar ruído a imagem e gerando a imagem
sigma = .1;
y = f + randn(n0,n0)*sigma;
% Mostrando a imagem
imageplot(f, 'Clean image', 1,2,1);
imageplot(clamp(y), ['Noisy image, SNR=' num2str(snr(f,y),3) 'dB'], 1,2,2);
% Computando os coeficienes de ruído
fW = Psi(y);
% Computando os valores de limite
T = 3*sigma;
fWT = Thresh(fW,T);
% Plotando os valores dos coeficientes
figure;subplot(1,2,1);
plot_wavelet(fW);
title('Original coefficients');
subplot(1,2,2);
plot_wavelet(fWT);
% Reconstrução da imagem e mostrando a remoção de ruído
f1 = PsiS(fWT);
figure;imageplot(clamp(y), 'Noisy image', 1,2,1);
imageplot(clamp(f1), strcat(['Denoising, SNR=' num2str(snr(f,f1),3) 'dB']), 1,2,2);