% Demonstration of Electrical Impedance Tomography using a 2D diffusion
% model in a rectangular grid.
%
% Samuli Siltanen March 2025

% Set the dimension of the model
N = 9;

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

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

% Save a picture of every measurement? Put 1 to both of these
allpics_flag = 1;
allpics_save_flag = 0;

% Define halfN so that N = 2*halfN+1
halfN = floor(N/2);

% Define the 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.
condmatH = ones(N,N-1);
condmatV = ones(N-1,N);
condmatH(1:3,1:3) = .1;
condmatV(1:3,1:3) = .1;
condmatH(end-2:end,end-2:end) = 3;
condmatV(end-2:end,end-2:end) = 3;

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

% Loop over skip current patterns and collect voltage vectors as columns
% in a voltage measurement matrix. Save images of intermediate results to
% folder images.
voltmat = zeros(Ne,Ne-1);
currmat = zeros(Ne,Ne-1);
skipn = 8;
for curind = 1:(Ne-1)

    % 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_in  = curind;
    if curind+skipn+1 <= Ne
        index_curr_out = curind+skipn+1;
    else
        index_curr_out = curind+skipn+1-Ne;
    end
    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:Ne));

    % Pick out the boundary voltages
    bdata_vol = u(1:Ne);

    % Save many images if user wants
    if allpics_flag
        % Record image of the solution including inside
        fignum = 1;
        plot_curr = 1;
        show_inside = 1;
        potmat = reshape(u((Ne+1):end),[N,N]);
        [currmatH,currmatV,bdata_cur] = FD_EIT_2D_currents_comp(N,potmat,bdata_vol,condmatH,condmatV);
        EIT_FD_demo2D_potential_plot(potmat,currmatH,currmatV,fignum,bdata_vol,bdata_cur, plot_curr,show_inside)
        if allpics_save_flag
            eval(['print -dpng images/FD_EIT_voltmat2_currim_inside_',num2str(curind),'.png'])
        end


        % Record image of the solution excluding inside
        fignum = 2;
        show_inside = 0;
        EIT_FD_demo2D_potential_plot(potmat,currmatH,currmatV,fignum,bdata_vol,bdata_cur, plot_curr,show_inside)
        if allpics_save_flag
            eval(['print -dpng images/FD_EIT_voltmat2_measim_',num2str(curind),'.png'])
        end
    end

    % Record the measurement and the current pattern
    voltmat(:,curind) = bdata_vol;
    currmat(:,curind) = rhs(1:Ne);

    % Monitor the run
    disp([curind Ne-1])

    %pause
end

% Save the results to file
save data/FD_EIT_2D_demo_voltmat2 currmat voltmat N condmatV condmatH I

% Save images of
% (i)   current patterns
% (ii)  conductivity
% (iii) voltage matrix

% Plot the voltage matrix
figure(3)
clf
imagesc(imresize(voltmat,10,'nearest'))
axis image
axis off
print -dpng images/FD_EIT_voltmat2_voltmat.png

% Plot the current patterns
figure(4)
clf
imagesc(imresize(currmat,10,'nearest'))
axis image
axis off
print -dpng images/FD_EIT_voltmat2_currmat.png

% Plot the conductivity
fignum = 5;
EIT_FD_demo2D_conductivity_plot(N,condmatH,condmatV,fignum)
print -dpng images/_voltmat2_conductivity.png

