EIDORS: Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software

Hosted by

Modeling EIT current flow in a human thorax model

Load thorax shape and identify contours

% get contours
thorax = shape_library('get','adult_male','boundary');
rlung  = shape_library('get','adult_male','right_lung');
llung  = shape_library('get','adult_male','left_lung');
% one could also run:
% shape_library('get','adult_male');
% to get all the info at once in a struct

% show the library image
shape_library('show','adult_male');
print_convert thoraxmdl01a.jpg '-density 100'

Use Netgen to build a FEM model of the thorax

shape = { 1,                      % height
          {thorax, rlung, llung}, % contours
          [4,50],                 % perform smoothing with 50 points
          0.04};                  % small maxh (fine mesh)

elec_pos = [ 16,                  % number of elecs per plane
             1,                   % equidistant spacing
             0.5]';               % a single z-plane
         
elec_shape = [0.05,               % radius
              0,                  % circular electrode
              0.01 ]';             % maxh (electrode refinement) 

fmdl = ng_mk_extruded_model(shape, elec_pos, elec_shape);
% this similar model is also available as:
% fmdl = mk_library_model('adult_male_16el_lungs');

[stim,meas_sel] = mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current'}, 1);
fmdl.stimulation = stim;

img=mk_image(fmdl,1);
img.elem_data(fmdl.mat_idx{2})= 0.3; % rlung
img.elem_data(fmdl.mat_idx{3})= 0.3; % llung

clf; show_fem(img); view(0,70);
print_convert thoraxmdl02a.jpg '-density 100'

Voltage and Current Distribution in a slice

img_v = img;
% Stimulate between elecs 16 and 5 to get more interesting pattern
img_v.fwd_model.stimulation(1).stim_pattern = sparse([16;5],1,[1,-1],16,1);
img_v.fwd_solve.get_all_meas = 1;
vh = fwd_solve(img_v);

img_v = rmfield(img, 'elem_data');
img_v.node_data = vh.volt(:,1);
img_v.calc_colours.npoints = 128;

PLANE= [inf,inf,0.35]; % show voltages on this slice

subplot(221);
show_slices(img_v,PLANE); axis off; axis equal
print_convert thoraxmdl03a.jpg
% 

img_v = img;
img_v.fwd_model.mdl_slice_mapper.npx = 64;
img_v.fwd_model.mdl_slice_mapper.npy = 64;
img_v.fwd_model.mdl_slice_mapper.level = PLANE;
q = show_current(img_v, vh.volt(:,1));
quiver(q.xp,q.yp, q.xc,q.yc,10,'b');
axis tight; axis image; ylim([-1 1]);axis off
print_convert thoraxmdl04a.jpg

Left Voltage distribution and Right Current distribution in slices at z= 0.35.

Current streamlines and the original image

img_v.fwd_model.mdl_slice_mapper.npx = 1000;
img_v.fwd_model.mdl_slice_mapper.npy = 1000;
img_v.fwd_model.mdl_slice_mapper.level = PLANE;

% Calculate at high resolution
q = show_current(img_v, vh.volt(:,1));

pic = shape_library('get','adult_male','pic');
imagesc(pic.X, pic.Y, pic.img);
% imgt= flipdim(imread('thorax-mdl.jpg'),1); imagesc(imgt);
colormap(gray(256)); set(gca,'YDir','normal');
hold on

sx = linspace(-.5,.5,15)';
sy = 0.05 + linspace(-.5,.5,15)';
hh=streamline(q.xp,q.yp, q.xc, q.yc,sx,sy); set(hh,'Linewidth',2, 'color','b');
hh=streamline(q.xp,q.yp,-q.xc,-q.yc,sx,sy); set(hh,'Linewidth',2, 'color','b');

axis equal; axis tight; axis off; print_convert thoraxmdl05a.jpg

Stream lines through z= 0.35. The path of the streamslines is limited by the density of the underlying FEM. With a finer mesh FEM, it would be possible to calculate streamlines that do not display path irregularities.

Add Equipotential lines to the original image

img_v = img;
% Stimulate between elecs 16 and 5 to get more interesting pattern
img_v.fwd_model.stimulation(1).stim_pattern = sparse([16;5],1,[1,-1],16,1);
img_v.fwd_solve.get_all_meas = 1;
vh = fwd_solve(img_v);

img_v = rmfield(img, 'elem_data');
img_v.node_data = vh.volt(:,1);
img_v.calc_colours.npoints = 256;
imgs = calc_slices(img_v,PLANE);


clf
imagesc(pic.X, pic.Y, pic.img); colormap(gray(256)); set(gca,'YDir','normal');
hh=streamline(q.xp,q.yp, q.xc, q.yc,sx,sy); set(hh,'Linewidth',2);
hh=streamline(q.xp,q.yp,-q.xc,-q.yc,sx,sy); set(hh,'Linewidth',2);

[x y] = meshgrid( linspace(pic.X(1), pic.X(2),size(imgs,1)), ...
                  linspace(pic.Y(2), pic.Y(1),size(imgs,2)));
hold on;
contour(x,y,imgs,31);
hh= findobj('Type','patch'); set(hh,'LineWidth',2)

hold off; axis off; axis equal; %ylim([50,450]);
print_convert thoraxmdl06a.jpg

Stream lines and equipotiential through z= 0.35.