close all
clear all
clc


% Uso delle Thin-plate smoothing spline

[A,map]= imread('lena.gif');        % caricamento immagine
[m n] = size(A);
rstep = 16;
cstep = 16;
rsudd = m/rstep;
csudd = n/cstep;

% visualizzo immagine originaria
B = im2double(A,'indexed');         % conversione in numeri reali
figure(1);
imshow(B,map);                      % visualizzazione

% immagine con rumore
Brum = B + randn(size(B))*50.0;
% Brum = B + randn(size(B))*10.0;
figure(2);
imshow(Brum,map);                    % visualizzazione

% immagine con rumore 
C = Brum;

% Ciclo sulle righe e sulle colonne della matrice partizionata in blocchi
% non sovrapposti
for i=1:rsudd,
    for j=1:csudd,
        i,j
        % estraggo parte di A
        i_index = (i-1)*rstep+1:i*rstep;
        j_index = (j-1)*cstep+1:j*cstep;
        Aproc = Brum(i_index,j_index);

        % Ricostruzione con thin-plate smoothing spline
        x = 1:size(Aproc,1);
        y = 1:size(Aproc,2);
        [Xgrid,Ygrid] = meshgrid(x,y);
        X = Xgrid(:);
        Y = Ygrid(:);
        Z = Aproc(:);
        
        % Thin-plate spline
        xy = [X,Y]';
        st = tpaps(xy,Z');

        % visualizzo soluzione
        points = fnplt(st);
        Bsmooth = griddata(points{1},points{2},points{3},Xgrid,Ygrid);
        C(i_index,j_index) = Bsmooth;
        
        % visualizzazione parziale ad ogni iterazione
        %figure(3);
        %imshow(C,map);                     
    end
end

% Visualizzo immagine elaborata
figure(3);
imshow(C,map);                      % visualizzazione
