
%
% EWA Splatting
% Entrada: Caminho = String entre aspas simples que representa o caminho de
% uma imagem RBG ou Gray scale, Exxemplo: 'xadrez.tif', 'xadrez.bmp'; 
% FTM = string entre aspas que representa o
% formato da imagem de entreda, Exemplo: 'tif', 'bmp',etc..
%
% SAIDA: Imagem em tons de cinza supostamente sem Aliasing
%
% Autor: Algemiro Augusto da Silva Neto
% Data: 12/Fevereiro/2007
% Pontifica Universidade Católica do Rio de Janeiro, INF2064
% 
% Para imagens maiores do que 256 x 256, dependendo da máquina
% pode demorar vários minutos
%
%


function Im2 = EWASplatting(Caminho,FMT)

% fecha todas as imagens
close all;

% Matriz de Projeção
H = [0.4 0.3 -0.3;
    -1.1 1.0 -0.6;
    -1.11 0.002 0.4];   

% Le uma imagem de arquivo
Im = imread(Caminho,FMT);

% transoforma para GrayScale
if ~isgray(Im)
  Im = rgb2gray(Im);
end

% exibe a imagem de entrada
imshow(Im);
Im = round(double(Im));

%% construir a gaussiana
[Lin, Col] = size(Im)
N = 2; % N*2 + 1 = 5 é a dimensão do suporte da Gaussiana do EWA
V = [0 0; 0 0]; % matriz para variancia
P = zeros(1,N*N); % vetor para desvio padrao da vizinhanca do ponto-base

%% aqui começa o EWA
for i=N+1:Lin-N-1
   for j=N+1:Col-N-1                   

      % calculo do desvio-padrao da vizinhança
      % que é tomado em torno do ponto-base
      k = 0;
      for p=i-N:i+N
         for q=j-N:j+N 
            k = k + 1;
            P(k) = sqrt(p^2 + q^2);
         end
      end
      
      %% calcula a matriz de variancia
      dp = std(P); % desvio padrao-da vizinhança
      var = dp^2;  % variancia
      V = zeros(2,2); % matriz de variancias
      V(1,1) = var;   % a matriz Diagonal contem as variancias
      V(2,2) = var;   % na diagonal, de acordo com a Tese do
                      % Zwicker, Seção 3.1.2
     
      % este trecho calcula a Jacobiana 
      % x0_d é um ponto posterior horizontal ao ponto base x0,
      % x0_a é um ponto anterior horizontal ao ponto-base x0,
      % d_x0 é o diferencial horizontal em torno de x0,
      % d_x1 é o diferencial vertical em torno de x0
      % d_u0 é o diferencial horizontal em torno de u0
      % d_u1 é o diferencial vertical em torno de u1
      % A matriz Jacobiana J é montada de acordo com
      % a equação 2.36 da Tese do Zwicker
      x0_d = H*[i j+1 0]'; x0_a = H*[i j-1 0]';
      d_x0 = x0_d(1) - x0_a(1);
      x1_d = H*[i+1 j 0]'; x1_a = H*[i-1 j 0]';
      d_x1 = x1_d(1) - x1_a(1);
      
      d_u0 = 2; d_u1 = 2;
      
      % A matriz Jacobiana J é montada de acordo com
      % a equação 2.36 da Tese do Zwicker      
      J = [d_x0/d_u0 d_x0/d_u1; d_x1/d_u0 d_x1/d_u1];
     
      % A matriz do EWA é montada de acordo com a Equação 
      % 2.51 da Tese do Zwicker
      V = J*V*J' + [1 0; 0 1];
      
      %% Cálculo da inversa da matriz V. Essa inversa, VI
      % será usada no expoente da gaussiana do EWA, de acordo
      % com a Equação 2.51 da Tese do Zwicker, bem como
      % o seu determinante DV
      VI = inv(V*10e+7);
      DV = det(V);
 
      %% Aqui começa a filtragem EWA 
      S = 0;
      n = 0;
      cont = 0;
      for p=i-N:i+N
          for q=j-N:j+N                                  
            if var > 0  
              % k é o espoente da Gaussiana que aqui é calculado 
              % separado
              k = (-1/2)*([p q]*VI*[p q]');  
              
              %% Esta é a Gaussiana              
              g = exp(k)/(2*pi*DV^(1/2)); 
              
              %% Este trecho corresponde 
              %% aos passos 1, 2 e 3 do EWA
              G = g * Im(p,q);
              %% contribuição da vizinhança
              %% de acordo com a equação
              %% 2.32 da Tese do Zwicker
              S = S + G;                                        
            end             
          end
      end
      % Passo 4 do EWA, reamostragem
      Im2(i,j) = Im(i,j) * S;         
                              
   end
end

%% O resultado final é normalizado para 
%% exibição no Matlab que precisa que os
%% valores de intensidade das imagens
%% estejam entre 0 e 1, por isso divido
%% tudo pelo valor máximo.
Im2 = Im2./max(max(Im2));

%% salva o resultado em disco
imwrite(Im2,'sem_aliasing.tif','tif');

%% exibe o resultado
figure;
imshow(Im2);

