% Demonstration of Electrical Impedance Tomography using a 2D diffusion
% model in a rectangular grid. Here we show how EIT is nonlinear. 
%
% We construct three conductivities c1, c2 and c3 for which 
%
%    c3 = c1+c2.
%
% We simulate EIT measurements for all of them and see that they do not add
% up linearly. 
%
% Samuli Siltanen March 2025

%% Preliminaries

% Set the dimension of the model
N = 3;

% Number of electrodes
Ne = 4*(N-1);

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

% Font size in plots
fsize = 12;


%% Simulations with first conductivity 

% 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.
condmatH1 = ones(N,N-1);
condmatV1 = ones(N-1,N);
% Low conductivity inclusion
condmatH1(1,:) = .5;
condmatV1(1,:) = .5;
% High conductivity inclusion
condmatH1(end,:) = 1.5;
condmatV1(end,:) = 1.5;


% Construct system matrix
A1 = FD_EIT_2D_system_matrix(N,condmatH1,condmatV1);

% Construct right-hand side. First the part related to boundary conditions.
% We list them in a vector in counterclockwise order, starting from the
% middle entry on the right ("Eastern") edge. Positive value means current
% coming out, negative value current going into the domain.
rhs = zeros(Ne+N^2,1);
index_curr_out = 1; % Electrode 1
%index_curr_in  = Ne/2+1; % Opposite to electrode 1
index_curr_in  = 3; 
rhs(index_curr_in) = I;
rhs(index_curr_out) = -I;

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

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

% Save image of the solution including inside
disp('saving solution image')
fignum = 1;
plot_curr = 1;
show_inside = 1;
potmat1 = reshape(u1((Ne+1):end),[N,N]);
bdata_vol1 = u1(1:Ne);
[currmatH1,currmatV1,bdata_cur1] = FD_EIT_2D_currents_comp(N,potmat1,bdata_vol1,condmatH1,condmatV1);
EIT_FD_demo2D_potential_plot(potmat1,currmatH1,currmatV1,fignum,bdata_vol1,bdata_cur1, plot_curr,show_inside)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 1120;
im2 = im1(850+[1:imN],1460+[1:imN],:);
% Save
imwrite(im2,'images/EIT_FD_3x3_case3_currim_inside1.png')
disp('solution image saved')


% Record image of the solution excluding inside
fignum = 2;
show_inside = 0;
EIT_FD_demo2D_potential_plot(potmat1,currmatH1,currmatV1,fignum,bdata_vol1,bdata_cur1, plot_curr,show_inside)
eval(['print -dpng images/EIT_FD_3x3_case3_measim1.png'])

% Plot potentials as numbers
disp('saving potential numbers image')
fignum = 4;
plot_x = 0;
show_inside = 1;
EIT_FD_demo2D_text_plot(potmat1,fignum,bdata_vol1,plot_x,show_inside,fsize,3)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 1800;
im2 = im1(580+[1:imN],1140+[1:imN],:);
% Save
imwrite(im2,'images/EIT_FD_3x3_case3_potential_numbers1.png')
disp('potential numbers image saved')




% Plot the conductivity
disp('saving conductivity image')
fignum = 5;
EIT_FD_demo2D_conductivity_plot(N,condmatH1,condmatV1,fignum)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 1000;
im2 = im1(900+[1:imN],1510+[1:imN],:);
% Save
imwrite(im2,'images/EIT_FD_3x3_case3_conductivity1.png')
disp('conductivity image saved')



%% Simulations with second conductivity 

% Create conductivity 2 by flipping conductivity 1
condmatH2 = flipud(condmatH1);
condmatV2 = flipud(condmatV1);

% Construct system matrix
A2 = FD_EIT_2D_system_matrix(N,condmatH2,condmatV2);

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

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

% Save image of the solution including inside
disp('saving solution image')
fignum = 12;
plot_curr = 1;
show_inside = 1;
potmat2 = reshape(u2((Ne+1):end),[N,N]);
bdata_vol2 = u2(1:Ne);
[currmatH2,currmatV2,bdata_cur2] = FD_EIT_2D_currents_comp(N,potmat2,bdata_vol2,condmatH2,condmatV2);
EIT_FD_demo2D_potential_plot(potmat2,currmatH2,currmatV2,fignum,bdata_vol2,bdata_cur2, plot_curr,show_inside)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 1120;
im2 = im1(850+[1:imN],1460+[1:imN],:);
% Save
imwrite(im2,'images/EIT_FD_3x3_case3_currim_inside2.png')
disp('solution image saved')

% Record image of the solution excluding inside
fignum = 22;
show_inside = 0;
EIT_FD_demo2D_potential_plot(potmat2,currmatH2,currmatV2,fignum,bdata_vol2,bdata_cur2, plot_curr,show_inside)
eval(['print -dpng images/EIT_FD_3x3_case3_measim2.png'])

% Plot potentials as numbers
disp('saving potential numbers image')
fignum = 42;
plot_x = 0;
show_inside = 1;
EIT_FD_demo2D_text_plot(potmat2,fignum,bdata_vol2,plot_x,show_inside,fsize,3)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 1800;
im2 = im1(580+[1:imN],1140+[1:imN],:);
% Save
imwrite(im2,'images/EIT_FD_3x3_case3_potential_numbers2.png')
disp('potential numbers image saved')



% Plot the conductivity
disp('saving conductivity image')
fignum = 52;
EIT_FD_demo2D_conductivity_plot(N,condmatH2,condmatV2,fignum)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 1000;
im2 = im1(900+[1:imN],1510+[1:imN],:);
% Save
imwrite(im2,'images/EIT_FD_3x3_case3_conductivity2.png')
disp('conductivity image saved')



%% Simulations with constant conductivity with value 2 everywhere

% Here the idea is to take conductivity3 = conductivity1 + conductivity2,
% which is equal to 2 everywhere. Note that usually we assume that the
% conductivity equals 1 at the boundary edges between pixels and
% electrodes. But now we take the boundary conductivity to be 2. We need to
% use a special routine for constructing the system matrix. 

condmatH3 = 2*ones(N,N-1);
condmatV3 = 2*ones(N-1,N);

% Construct system matrix
A3 = FD_EIT_2D_system_matrix_boundarycond2(N,condmatH3,condmatV3);

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

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


% Save image of the solution including inside
disp('saving solution image')
fignum = 12;
plot_curr = 1;
show_inside = 1;
potmat3 = reshape(u3((Ne+1):end),[N,N]);
bdata_vol3 = u3(1:Ne);
[currmatH3,currmatV3,bdata_cur3] = FD_EIT_2D_currents_comp(N,potmat3,bdata_vol3,condmatH3,condmatV3);
EIT_FD_demo2D_potential_plot(potmat3,currmatH3,currmatV3,fignum,bdata_vol3,bdata_cur3, plot_curr,show_inside)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 1120;
im2 = im1(850+[1:imN],1460+[1:imN],:);
% Save
imwrite(im2,'images/EIT_FD_3x3_case3_currim_inside3.png')
disp('solution image saved')

% Record image of the solution excluding inside
fignum = 22;
show_inside = 0;
EIT_FD_demo2D_potential_plot(potmat3,currmatH3,currmatV3,fignum,bdata_vol3,bdata_cur3, plot_curr,show_inside)
eval(['print -dpng images/EIT_FD_3x3_case3_measim3.png'])

% Plot potentials as numbers
disp('saving potential numbers image')
fignum = 42;
plot_x = 0;
show_inside = 1;
EIT_FD_demo2D_text_plot(potmat3,fignum,bdata_vol3,plot_x,show_inside,fsize,3)
% Initial image
im1 = print('-r500','-RGBImage');
[row1,col1] = size(im1(:,:,1));
% Crop
imN = 1800;
im2 = im1(580+[1:imN],1140+[1:imN],:);
% Save
imwrite(im2,'images/EIT_FD_3x3_case3_potential_numbers3.png')
disp('potential numbers image saved')



