EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

Hosted by

Electrode movement reconstruction for simulated 3D data

% Simulate 3D movement $Id: move_3d01.m 1535 2008-07-26 15:36:27Z aadler $

noiselev = .1;
movement = 2;

% Generate eidors 3D finite element model
mdl3dim = mk_common_model( 'n3r2' );
mdl3dim.fwd_model.nodes(:,3) = mdl3dim.fwd_model.nodes(:,3)/3;

% Instantiate a homogeneous forward model.
% ref_level is 1 since we use ones( ).
sigma = ones( size(mdl3dim.fwd_model.elems,1) ,1);
% Create eidors_obj type image
img = eidors_obj('image', 'homogeneous image', ...
    'elem_data', sigma, ...
    'fwd_model', mdl3dim.fwd_model, ...
    'ref_level', 1);
% Solve the homogeneous forward problem
vh = fwd_solve( img );

% Set local inhomogeneites using eidors 3D model
load('datacom.mat','A','B');
sigma(A) = 1.2; 
sigma(B) = 0.8;
img.elem_data = sigma;
node0 = img.fwd_model.nodes;  % Node variable before movement
node1 = node0;                % Node variable after movement
z_axis = node1(:,3);

exaggeration = 10;
movement= 0.01*movement;
% Do a 3D twist - exaggerated for clearer illustration of distortion
node1(:,1) = node0(:,1).*(1 + exaggeration*movement*z_axis);
node1(:,2) = node0(:,2).*(1 + exaggeration*movement*(1-z_axis));

img.fwd_model.nodes = node1;
show_fem( img );
xlabel('x'); ylabel('y');
view(-44,22)

print -r75 -dpng move_3d01.png

% Do a 3D twist - we'll actually use this one
centre = 1- movement/2;
node1(:,1) = node0(:,1).*(centre + movement*z_axis);
node1(:,2) = node0(:,2).*(centre + movement*(1-z_axis));

% Solve inhomogeneous forward problem with movements and normal noise.
img.fwd_model.nodes = node1;
vi = fwd_solve( img );
noise = noiselev*std( vh.meas - vi.meas )*randn( size(vi.meas) );
vi.meas = vi.meas + noise;
move = node1 - node0;

Figure: Illustration of the deformation of the 3D model (exaggerated 10x).

% Reconstruct 3D movement $Id: move_3d02.m 1535 2008-07-26 15:36:27Z aadler $

img3dim= img;

mdl3dim.RtR_prior = 'laplace_image_prior';
mdl3dim.hyperparameter.value = 3e-3;

% Show slices of 3D model with true movement vectors
subplot(1,3,1);
img3dim.elem_data = img3dim.elem_data - 1;
img3dim.calc_colours.backgnd=[.8 .8 .9];
show_slices_move( img3dim, move );

% Inverse solution of data without movement consideration
img3dim= inv_solve(mdl3dim, vh, vi);
img3dim.calc_colours.backgnd=[.8 .8 .9];
subplot(1,3,2)
show_slices_move( img3dim );

% Inverse solution of data with movement consideration
move_vs_conduct = 20;  % Movement penalty (symbol mu in paper)
% Define a eidglobalors_obj Movement model solved by electrode movement 
% algorithms.
mdlM = mdl3dim;
mdlM.fwd_model.conductivity_jacobian = mdlM.fwd_model.jacobian;

% this is a perturbation type jacobian
%mdlM.fwd_model.jacobian = 'aa_e_move_jacobian';

mdlM.fwd_model.jacobian = @calc_move_jacobian; % faster / more accurate
mdlM.RtR_prior = 'aa_e_move_image_prior';
mdlM.aa_e_move_image_prior.parameters = move_vs_conduct;
% Solve inversglobale problem and show slices
imgM = inv_solve(mdlM, vh, vi);
imgM.calc_colours.backgnd=[.8 .8 .9];
subplot(1,3,3)
show_slices_move( imgM );
print -r100 -dpng move_3d02.png

Figure: Inverse solutions of the problem above. Each column shows three cross-sectional slices of the model. Forward problem with true electrode movement (scaled 10x) (left column). Inverse solution without movement compensation (centre column) and with movement estimation and compensation (scaled 10x) (right column).