EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

Hosted by

2D Geophysical models using square mesh elements

Create fine mesh

% $Id: square_mesh01.m 1535 2008-07-26 15:36:27Z aadler $

z_contact= 0.01;
n_elec= 17;
nodes_per_elec= 5;
elec_width= 0.2;
elec_spacing= 1.0;

xllim=-12; xrlim= 12; ydepth=-15;
[x,y] = meshgrid( linspace(xllim,xrlim,49), linspace(ydepth,0,31) );
vtx= [x(:),y(:)];
% Refine points close to electrodes - don't worry if points overlap
[x,y] = meshgrid( -9:.25:9, -3:.25:0 );
vtx= [vtx; x(:),y(:)];

xgrid=  linspace(-elec_width/2, +elec_width/2, nodes_per_elec)';
x2grid= elec_width* [-5,-4,-3,3,4,5]'/4;
for i=1:n_elec
% Electrode centre
  x0= (i-1-(n_elec-1)/2)*elec_spacing;
  y0=0;
  elec_nodes{i}= [x0+ xgrid, y0+0*xgrid];
  vtx= [ vtx; ...
        [x0 + x2grid   ,  y0             + 0*x2grid];
        [x0 + xgrid*1.5,  y0-elec_width/2+ 0*xgrid];
        [x0 + x2grid*1.5, y0-elec_width/2+ 0*x2grid];
        [x0 + xgrid*2   , y0-elec_width  + 0*xgrid];
        [x0 + xgrid*2   , y0-elec_width*2+ 0*xgrid]];
end

fmdl= mk_fmdl_from_nodes( vtx, elec_nodes, z_contact, 'sq_m1');


subplot(121)
show_fem(fmdl); axis image
subplot(122)
show_fem(fmdl); axis image; axis([-2 2 -2.5 0.5]);

print -dpng -r150 square_mesh01a.png

Figure: Left: Fine mesh model with electrodes at surface Right: close up view of mesh near electrodes

Create Dual Mesh

% Create and show square model $Id: square_mesh02.m 1535 2008-07-26 15:36:27Z aadler $

% Create square mesh model
[cmdl,c2f]= mk_grid_model(fmdl, linspace(-8,8,17), linspace(-11.5,-0.5,13) );

clf;
show_fem(fmdl);
hold on;
h= trimesh(cmdl.elems,cmdl.nodes(:,1),cmdl.nodes(:,2));
set(h,'Color',[0,0,1],'LineWidth',2);
hold off
axis image

print -dpng -r125 square_mesh02a.png

Figure: Left: Uniform mesh density Right: Mesh density refined going from surface to bottom

Create Simulation Pattern

% simulate targets $Id: square_mesh03.m 1887 2009-06-23 20:05:12Z aadler $

fmdl.stimulation= mk_stim_patterns(length(elec_nodes), 1, '{ad}','{ad}', {}, 1);

img= eidors_obj('image','fmdl','fwd_model',fmdl, ...
                'elem_data', ones(size(fmdl.elems,1),1) );
vh= fwd_solve(img);

% interpolate onto mesh
xym= interp_mesh( fmdl, 3);
x_xym= xym(:,1,:); y_xym= xym(:,2,:);
% non-conductive target
ff  = (x_xym>-3) & (x_xym<-2) & (y_xym<-4) & (y_xym>-7);
img.elem_data= img.elem_data - 0.1*mean(ff,3);
% conductive target
ff  = (x_xym> 2) & (x_xym< 4) & (y_xym<-5) & (y_xym>-7);
img.elem_data= img.elem_data + 0.1*mean(ff,3);

% inhomogeneous image
vi= fwd_solve(img);

show_fem(img); axis image;
print -dpng -r125 square_mesh03a.png

Figure: Left: Uniform mesh density Right: Mesh density refined going from surface to bottom

Inverse Solution

% 2D solver $Id: square_mesh04.m 1535 2008-07-26 15:36:27Z aadler $

% Create a new inverse model, and set
% reconstruction model and fwd_model
imdl= mk_common_model('c2c2',16);
imdl.rec_model= cmdl;
imdl.fwd_model= fmdl;

c2f= mk_coarse_fine_mapping( fmdl, cmdl);
imdl.fwd_model.coarse2fine = c2f;
%imdl.RtR_prior = @gaussian_HPF_prior;
imdl.RtR_prior = @noser_image_prior;
imdl.prior_use_fwd_not_rec = 1;
imdl.noser_image_prior.exponent= 0.5;
imdl.solve = @aa_inv_solve;
imdl.hyperparameter.value= 0.003;

imgs= inv_solve(imdl, vh, vi);

show_fem(fmdl); ax= axis;
hold on
show_fem(imgs);
hold on;
h= trimesh(cmdl.elems,cmdl.nodes(:,1),cmdl.nodes(:,2));
set(h,'Color',[0,0,1]);
hold off
hold off
axis(ax); axis image
print -dpng -r125 square_mesh04a.png

Figure: Left: Uniform mesh density Right: Mesh density refined going from surface to bottom

Internal Electrodes

% $Id: square_mesh05.m 1886 2009-06-23 19:21:55Z aadler $

z_contact= 0.01;
nodes_per_elec= 5;
n_elec=17;
elec_width= 0.2;
elec_spacing= 1.0;

xllim=-12; xrlim= 12; ydepth=-15;
[x,y] = meshgrid( linspace(xllim,xrlim,49), linspace(ydepth,0,31) );
vtx= [x(:),y(:)];
% Refine points close to electrodes - don't worry if points overlap
[x,y] = meshgrid( -9:.25:9, -3:.25:0 );
vtx= [vtx; x(:),y(:)];
[x,y] = meshgrid( -9:.25:9, -15:.25:-11 );
vtx= [vtx; x(:),y(:)];

xgrid=  linspace(-elec_width/2, +elec_width/2, nodes_per_elec)';
x2grid= elec_width* [-5,-4,-3,3,4,5]'/4;
for i=1:n_elec
  x0= (i-1-(n_elec-1)/2)*elec_spacing;

% Top electrode
  y0  = zeros(size(xgrid));
  y0_2= zeros(size(x2grid));
% elec_nodes{2*i-1}= [x0+ xgrid, y0];
  elec_nodes{2*i-1}= [x0, 0];
  vtx= [ vtx; ...
        [x0 + x2grid    , y0_2              ];
        [x0 + xgrid*1.5 , y0  - elec_width/2];
        [x0 + x2grid*1.5, y0_2- elec_width/2];
        [x0 + xgrid*2   , y0  - elec_width  ];
        [x0 + xgrid*2   , y0  - elec_width*2]];

% Bottom electrode
  y0  = -13*ones(size(xgrid));
  y0_2= -13*ones(size(x2grid));
  elec_nodes{2*i}= [x0,-13]; % Only point electrodes insidq
  vtx= [ vtx; ...
        [x0 + x2grid    , y0_2              ];
        [x0 + xgrid*1.5 , y0  - elec_width/2];
        [x0 + x2grid*1.5, y0_2- elec_width/2];
        [x0 + xgrid*1.5 , y0  + elec_width/2];
        [x0 + x2grid*1.5, y0_2+ elec_width/2];];
end

fmdl= mk_fmdl_from_nodes( vtx, elec_nodes, z_contact, 'sq_m1');
fmdl.solve=@aa_fwd_solve;
fmdl.system_mat=@aa_calc_system_mat;
fmdl.jacobian=@aa_calc_jacobian;


subplot(121)
show_fem(fmdl); axis image
subplot(122)
show_fem(fmdl); axis image; axis([-2 2 -14 -12]);

print -dpng -r150 square_mesh05a.png

Figure: Left: Uniform mesh density Right: Mesh density refined going from surface to bottom

Figure: Left: Uniform mesh density Right: Mesh density refined going from surface to bottom

Figure: Left: Uniform mesh density Right: Mesh density refined going from surface to bottom

Figure: Left: Uniform mesh density Right: Mesh density refined going from surface to bottom