% Demonstration of Electrical Impedance Tomography using a 2D diffusion
% model in a rectangular grid. We apply a maximal linearly independent set
% of current patterns and save the resulting voltage vectors as columns in
% a voltage matrix. 
%
% 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);

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

% Loop over trig 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);
theta = 2*pi*[0:(Ne-1)]/Ne; % "Angle" for trig patterns
for curind = 1:(Ne-1)

    % Initialize right-hand side. 
    rhs = zeros(Ne+N^2,1);

    % 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.
    % Here we use the trigonometric current patterns although the domain is
    % a square instead of a disc. 
    if mod(curind,2)==1
        trig_curr_pattern = I*cos((curind+1)*theta/2);
    else
        trig_curr_pattern = I*sin(curind*theta/2);
    end
    rhs(1:Ne) = trig_curr_pattern(:);

    % 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_voltmat_trig1_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_voltmat_trig1_measim_',num2str(curind),'.png'])
        end
    end

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

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

    pause
end

% Save the results to file
save data/FD_EIT_2D_demo_voltmat_trig1 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_voltmat_trig1_voltmat.png

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

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

