EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

Hosted by

Image reconstruction of moving objects

Andy Adler, Tao Dai, William R.B. Lionheart Temporal Image Reconstruction in Electrical Impedance Tomography Physiol. Meas., 2007.

Sample Data

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')
    [vh,vi,xyr_pt]=simulate_2d_movement;
    save simulate_2d_movement_data vi vh xyr_pt
else
    load simulate_2d_movement_data
end

% 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, y_k = Jx_k, 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 x_t from data [y_t-d; ... ; y_t-1; y_t; y_t+1; ... ; y_t+d].

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

1. Gauss-Newton (GN) inverse, using y_t only;
2. GN with weighted data, using a weighted average of y_t-d ... y_t+d;
3. Kalman filter inverse, using all previous and current data, y₀ ... y_t;
4. Temporal inverse, using y_t-d ... y_t+d based on a temporal prior model.

% 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

% 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;

   end

   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);

   show_fem(img);
   axis equal
   axis([-1.1, 0.1, -1.1, 1.1]);
   set(gca,{'XTicklabel','YTicklabel'},{'',''});
%
% 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));
       set(hh,'LineWidth',3,'Color',[0,0,0]);
   end
   hold off;

   print_convert(sprintf('temporal_solver03%c.png',96+alg));
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

% 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
 k=1;
 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);
   k=k+1;
 end
 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);