EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

EIDORS (mirror)
− Image Reconst
− Data Structures
− Applications
− FEM Modelling
− Old tutorials
Contrib Data
Browse Docs
Browse SVN

Mailing list


Hosted by
SourceForge.net Logo


Image reconstruction of moving objects

These results are taken from the paper:
    Andy Adler, Tao Dai, William R.B. Lionheart Temporal Image Reconstruction in Electrical Impedance Tomography Physiol. Meas., 2007.
Electrical Impedance Tomography (EIT) calculates images of the body from body impedance measurements. While the spatial resolution of these images is relatively low, the temporal resolution of EIT data can be high. Most EIT reconstruction algorithms solve each data frame independently, although Kalman filter algorithms track the image changes across frames. This paper proposes a new approach which directly accounts for correlations between images in successive data frames. Image reconstruction is posed in terms of an augmented image x and measurement vector y, which concatenate the values from the d previous and future frames. Image reconstruction is then based on an augmented regularization matrix R, which accounts for a model of both the spatial and temporal correlations between image elements. Results are compared for reconstruction algorithms based on independent frames, Kalman filters, and the proposed approach. For low values of the regularization hyperparameter, the proposed approach performs similarly to independent frames, but for higher hyperparameter values, it uses adjacent frame data to reduce reconstructed image noise.

Sample Data

Simulated data are calculated using the function simulate_2d_movement, which models a small (0.05×radius) object moving in a circle in a 2D circular tank. These data are available here.

Since the temporal solver works best for noisy data, we simulate a fairly low SNR (high NSR). To understand why the temporal solver works best at for noisy data, consider that if data are clean, then there is no reason to look at data from the past/future, since the current measurement tells us everything we need to know.

% Sample Data $Id: temporal_solver01.m 3140 2012-06-09 10:55:15Z bgrychtol $

if ~exist('simulate_2d_movement_data.mat')
    save simulate_2d_movement_data vi vh xyr_pt
    load simulate_2d_movement_data

% Temporal solver works best at large Noise
nsr= 4.0; %Noise to Signal Ratio

signal = mean( abs(mean(vi,2) - vh) ); % remember to define signal
% Only add noise to vi. This is reasonable, since we have
% lots of data of vh to average to reduce noise
vi_n= vi + nsr*signal*randn( size(vi) );

Reconstruction Algorithms

Over time steps, k, a sequence of difference vectors, yk = Jxk, are measured (assuming the body and electrode geometry, and thus J, stay fixed). If the conductivity of the body under investigation doesn’t change too rapidly, then it is reasonable to expect that a certain number of measurements, d, into the past and future provide useful information about the current image. Labelling the current instant as t, we therefore seek to estimate xt from data [yt-d; ... ; yt-1; yt; yt+1; ... ; yt+d].

In the subsequent sections we consider three traditional approaches and the proposed temporal inverse; each estimates xt at frame t from a sequence of data starting at frame 0, using the indicated data:

  1. Gauss-Newton (GN) inverse, using yt only;
  2. GN with weighted data, using a weighted average of yt-d ... yt+d;
  3. Kalman filter inverse, using all previous and current data, y0 ... yt;
  4. Temporal inverse, using yt-d ... yt+d based on a temporal prior model.
The following code defines the reconstruction algorithms:
% Image reconstruction of moving objects $Id: temporal_solver02.m 4839 2015-03-30 07:44:50Z aadler $

time_steps=  3;
time_weight= .8;

base_model = mk_common_model( 'c2c2', 16 ); % 576 element
base_model.fwd_model = mdl_normalize(base_model.fwd_model, 0);
% HP value gives good solutions (calculated via Noise Figure)
hp= sqrt(0.101046); % hp is squared
base_model.hyperparameter.value= hp;

% GN Solver
imdl_GN = base_model;
imdl_GN.RtR_prior= @prior_noser;
imdl_GN.prior_noser.exponent= .5;
imdl_GN.solve= @inv_solve_diff_GN_one_step;
imdl_GN.hyperparameter.value= hp;
imdl_GN.fwd_model = mdl_normalize(imdl_GN.fwd_model, 0);

% Temporal Solver
imdl_TS = base_model;
imdl_TS.RtR_prior= @prior_time_smooth;
imdl_TS.prior_time_smooth.space_prior= @prior_noser;
imdl_TS.prior_noser.exponent= .5;
imdl_TS.prior_time_smooth.time_weight= time_weight;
imdl_TS.prior_time_smooth.time_steps=  time_steps;
imdl_TS.solve= @inv_solve_time_prior;
imdl_TS.inv_solve_time_prior.time_steps=   time_steps;

% Kalman Solver
imdl_KS = base_model;
imdl_KS.RtR_prior= @prior_noser;
imdl_KS.prior_noser.exponent= .5;
imdl_KS.solve= @inv_solve_diff_kalman;

Image Reconstruction

We reconstruct the image at 9 O'clock using each algorithm, using the code:
% Image reconstruction of moving objects $Id: temporal_solver03.m 3346 2012-07-01 21:30:55Z bgrychtol $

image_select= length(xyr_pt)/2+1;; % this image is at 9 O'Clock
time_steps=  3; ts_expand= 5;
time_weight= .8;
ts_vec= -time_steps:time_steps;

% vi_sel is the inhomog data used by the algorithm
% sel is the image to show (ie. last for kalman, middle for temporal)
for alg=1:4
   if     alg==1; % GN Solver
      im_sel = image_select;
      vi_sel = vi_n(:,im_sel);

      sel  = 1;
      imdl= imdl_GN;

   elseif alg==2  % Weighted GN Solver
      im_sel= image_select+ ts_vec*ts_expand;
      weight= (time_weight.^abs(ts_vec));
      vi_sel= vi_n(:,im_sel) * weight(:) / sum(weight);

      sel  = 1;
      imdl= imdl_GN;

   elseif alg==3  % Temporal Solver
      im_sel= image_select+ ts_vec*ts_expand;
      vi_sel= vi_n(:,im_sel);

      sel  = 1 + time_steps; % choose the middle
      imdl= imdl_TS;

   elseif alg==4  % Kalman Solver
      im_sel= image_select+ (-12:0)*ts_expand; %let Kalman warm up
      vi_sel= vi_n(:,im_sel);

      sel  = length(im_sel); % choose the last
      imdl= imdl_KS;


   imdl.fwd_model = mdl_normalize(imdl.fwd_model, 0);
   img= inv_solve( imdl, vh, vi_sel);
   % only show image for sel (ie. last for kalman, middle for temporal)
   img.elem_data= img.elem_data(:,sel);

   axis equal
   axis([-1.1, 0.1, -1.1, 1.1]);
% Put circles where the data points used for reconstruction are
   theta= linspace(0,2*pi,length(xyr_pt)/ts_expand);
   xr=   cos(theta);       yr=   sin(theta);
   xpos= xyr_pt(1,im_sel); ypos= xyr_pt(2,im_sel);
   rad= xyr_pt(3,im_sel);
   hold on;
   for i=1:length(xpos)
       hh= plot(rad(i)*xr+ xpos(i),rad(i)*yr+ ypos(i));
   hold off;

end % for i

Figure: Image reconstructions of a moving ball in a medium. Top: using a classic one-stop solver, Bottom: using a Kalman solver. Note the warm up period of the Kalman images when the time correlations are being trained.

Image Reconstruction Movies

In order to show animated movies of temporal solvers, we can do the following
% Image reconstruction of moving objects $Id: temporal_solver04.m 1535 2008-07-26 15:36:27Z aadler $

time_steps=  3; ts_expand= 5;
time_weight= .8;
ts_vec= -time_steps:time_steps;

image_select= .25*length(xyr_pt)+1:ts_expand:.75*length(xyr_pt)+1;

% GN Solver
 vi_sel = vi(:,image_select);
 img= inv_solve( imdl_GN, vh, vi_sel);
 animate_reconstructions('temporal_solver04a', img);

% Temporal Solver
 for i= image_select
   vi_sel= vi(:,i+ts_vec);
   imgs= inv_solve( imdl_TS, vh, vi_sel);
   img.elem_data(:,k)= imgs.elem_data(:,1+time_steps);
 animate_reconstructions('temporal_solver04c', img);

% Kalman Solver
 vi_sel = vi(:,image_select);
 img= inv_solve( imdl_KS, vh, vi_sel);
 animate_reconstructions('temporal_solver04d', img);

Output images are:
Gauss-Newton solver
Temporal solver
Kalman solver

Last Modified: $Date: 2017-02-28 13:12:08 -0500 (Tue, 28 Feb 2017) $ by $Author: aadler $