% Test linearization as EIT reconstruction
%
% Samuli Siltanen May 2025



%% Preliminaries

% Pixel grid size 
n = 3;

% Number of electrodes
L = 4*(n-1);

% Injected current amplitude (milli-Amperes)
I = 1;



%% Calculation of the homogeneous case

% Construct simple conductance matrices
condmatH0 = ones(n,n-1);
condmatV0 = ones(n-1,n);

% Construct system matrix
A0 = FD_EIT_2D_system_matrix(n,condmatH0,condmatV0);

% Initialize voltage matrix
Uadj0 = zeros(L,L-1);

% Loop over adjacent current patterns
for iii = 1:(L-1)
    % Current pattern
    index_curr_in  = iii;
    index_curr_out = iii+1;

    % Construct right-hand side.
    rhs = zeros(L+n^2,1);
    rhs(index_curr_in) = I;
    rhs(index_curr_out) = -I;

    % Solve the FD equation, using pseudoinverse because dim(ker(A))=1
    u = pinv(A0)*rhs;

    % The solution is unique up to a constant. Make it unique by subtracting
    % the mean value of voltages at electrodes
    u = u-mean(u(1:L));

    % Insert solution as a column to the voltage matrix
    Uadj0(:,iii) = u(1:L);
end
% Drop the result to vertical vector
Uadj0 = Uadj0(:);


%% Computation of non-homogeneous case

% condmatH: conductivity values between horizontal neighbors, nx(n-1) matrix
% condmatV: conductivity values between vertical neighbors, (n-1)x n matrix
% The conductivity between boundary pixels and outside is assumed to be one.
condmatH = ones(n,n-1);
condmatV = ones(n-1,n);
% Low conductivity inclusion
condmatH(1,:) = .4;
condmatV(1,:) = .4;
% High conductivity inclusion
condmatH(end,:) = 9;
condmatV(end,:) = 9;

% Construct system matrix
A = FD_EIT_2D_system_matrix(n,condmatH,condmatV);

% Initialize voltage matrix
Uadj = zeros(L,L-1);

% Loop over adjacent current patterns
for iii = 1:(L-1)
    % Current pattern
    index_curr_in  = iii;
    index_curr_out = iii+1;

    % Construct right-hand side.
    rhs = zeros(L+n^2,1);
    rhs(index_curr_in) = I;
    rhs(index_curr_out) = -I;

    % Solve the FD equation, using pseudoinverse because dim(ker(A))=1
    u = pinv(A)*rhs;

    % The solution is unique up to a constant. Make it unique by subtracting
    % the mean value of voltages at electrodes
    u = u-mean(u(1:L));

    % Insert solution as a column to the voltage matrix
    Uadj(:,iii) = u(1:L);
end
% Drop the result to vertical vector
Uadj = Uadj(:);


%% Linearized reconstruction

% Calculate Jacobian
Jcbn = EIT_FD_demo2D_forward_map_Jacobian_adj(n,.01);

% Apply pseudoinverse
recon = pinv(Jcbn)*(Uadj-Uadj0);

% Collect result into appropriate matrix form
condmatH_recon = reshape(1+recon(1:n*(n-1)),[n,n-1]);
condmatV_recon = reshape(1+recon((n*(n-1)+1):end),[n-1,n]);

% Record minimum and maximum conductances
scaleMIN = min([min(condmatH(:)),min(condmatV(:)),min(condmatH_recon(:)),min(condmatV_recon(:))]);
scaleMAX = max([max(condmatH(:)),max(condmatV(:)),max(condmatH_recon(:)),max(condmatV_recon(:))]);

%disp([condmatH,condmatH_recon])
%disp([condmatV,condmatV_recon])


% Show original conductivity in colors
EIT_FD_demo2D_recon_colors_plot(n,condmatH,condmatV,1,scaleMIN,scaleMAX)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 850;
im2 = im1(970+[1:imN],1580+[1:imN],:);
% Save
imwrite(im2,'images/FD_EIT_2D_demoE_cond3x3largecontrast_orig.png')

% Show reconstructed conductivity in colors
EIT_FD_demo2D_recon_colors_plot(n,condmatH_recon,condmatV_recon,2,scaleMIN,scaleMAX)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 850;
im2 = im1(970+[1:imN],1580+[1:imN],:);
% Save
imwrite(im2,'images/FD_EIT_2D_demoE_cond3x3largecontrast_recon.png')


% Show original conductivity in numbers
fsize = 11;
dec_acc = 2;
x_displmt = -.3;
y_displmt = 0;
EIT_FD_demo2D_recon_numbers_plot(n,condmatH,condmatV,3,fsize,dec_acc,x_displmt,y_displmt)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 850;
im2 = im1(970+[1:imN],1580+[1:imN],:);
% Save
imwrite(im2,'images/FD_EIT_2D_demoE_cond3x3largecontrast_orig_numbers.png')

% Show reconstructed conductivity in numbers
EIT_FD_demo2D_recon_numbers_plot(n,condmatH_recon,condmatV_recon,4,fsize,dec_acc,x_displmt,y_displmt)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 850;
im2 = im1(970+[1:imN],1580+[1:imN],:);
% Save
imwrite(im2,'images/FD_EIT_2D_demoE_cond3x3largecontrast_recon_numbers.png')