% Supporting routine for a simple 2D finite difference (FD) solver for the 
% conductivity equation in a square domain. We divide the square into NxN 
% pixels and assume that the pixel width is one. 
%
% It is convenient from the point of view of Electrical Impedance
% Tomography (EIT) to think of electrodes around the square. The 4*N pixel 
% boundary segments located at the edge of the pixelated square are models
% of electrodes. 
%
% Further, we take N to be odd so that there is a unique center pixel on
% the right (Eastern) boundary of the grid square. That is electrode number
% one. Electrode numbering goes couterclockwise around the square, with the
% last electrode having number 4*(N-1).
%
% In some visualization tasks it is nice to have a vector of voltages
% measured at the electrodes to be presented in four chunks, one for each of
% the square boundaries. West, East, North and South. The corners are so
% placed that 
% Northeast corner belongs to Eastern edge
% Northwest corner belongs to Northern edge
% Southwest corner belongs to Western edge
% Southeast corner belongs to Southern edge
%
% See also the opposite routine: edges2vec.m
% 
% Arguments:
% N      the computational grid has NxN pixels
% bdata  voltages at electrodes listed counterclockwise, 4*(N-1) vector
%
% Returns:
% bdata_W   voltages on the left, size (N-1)x1, listed from top to bottom
% bdata_E   voltages on the right, size (N-1)x1, listed from top to bottom
% bdata_N   voltages on the top, size 1x(N-1), listed from left to right
% bdata_S   voltages on the bottom, size 1x(N-1), listed from left to right
%
% Samuli Siltanen March 2025


function [bdata_W,bdata_E,bdata_N,bdata_S] = vec2edges(N,bdata)

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

bdata_W = bdata((halfN+1)+(N-1)+[1:(N-1)]);
bdata_W = bdata_W(:); % Size (N-1)x1

bdata_E = zeros(N,1);
tmp1 = bdata(1:(halfN+1));
tmp1 = tmp1(:);
bdata_E(1:(halfN+1)) = flipud(tmp1);
tmp2 = bdata((halfN+1)+3*(N-1)+[1:(halfN-1)]);
tmp2 = tmp2(:);
bdata_E((halfN+2):(N-1)) = flipud(tmp2);

bdata_N = bdata((halfN+1)+[1:(N-1)]);
bdata_N = bdata_N(:).'; % Size 1x(N-1)
bdata_N = fliplr(bdata_N); % Correct order

bdata_S = bdata((halfN+1)+2*(N-1)+[1:(N-1)]);
bdata_S = bdata_S(:).'; % Size 1x(N-1)