demo_real

 [inhomg_img, demo_img] = demo_real;
 DEMO to show usage of EIDORS3D

0001 function [inhomg_img, demo_img] = demo_real;
0002 % [inhomg_img, demo_img] = demo_real;
0003 % DEMO to show usage of EIDORS3D
0004 
0005 % (C) 2005 Nick Polydorides + Andy Adler. License: GPL version 2 or version 3
0006 % $Id: demo_real.m 3630 2012-11-18 18:29:52Z aadler $
0007 
0008 eidors_msg('log_level',2); % most messages
0009 
0010 disp('step 1: create FEM model structure');
0011 %
0012 % create a 'fwd_model' object with name demo_mdl
0013 %
0014 
0015 [vtx,simp, bdy] = get_model_elems;
0016 
0017 demo_mdl.name = 'demo real model';
0018 demo_mdl.nodes= vtx;
0019 demo_mdl.elems= simp;
0020 demo_mdl.boundary= bdy;
0021 demo_mdl.solve=      'np_fwd_solve';
0022 demo_mdl.jacobian=   'np_calc_jacobian';
0023 demo_mdl.system_mat= 'np_calc_system_mat';
0024 
0025 disp('step 2: create FEM model electrodes definitions');
0026 
0027 [gnd_ind, electrodes, perm_sym, elec, protocol, no_pl] = get_model_elecs;
0028 demo_mdl.gnd_node=           gnd_ind;
0029 demo_mdl.electrode =         electrodes;
0030 demo_mdl.np_fwd_solve.perm_sym =          perm_sym;
0031 
0032 disp('step 3: create FEM model stimulation and measurement patterns');
0033 
0034 [stimulations ] = get_model_stim( demo_mdl );
0035 demo_mdl.stimulation= stimulations;
0036 
0037 demo_mdl= eidors_obj('fwd_model', demo_mdl); %create object
0038 
0039 disp('step 4: simulate data for homogeneous medium');
0040 
0041 %
0042 % create a homogeneous image
0043 %
0044 
0045 mat= ones( size(demo_mdl.elems,1) ,1);
0046 
0047 homg_img= eidors_obj('image', 'homogeneous image', ...
0048                      'elem_data', mat, ...
0049                      'fwd_model', demo_mdl );
0050 
0051 homg_data=fwd_solve( demo_mdl, homg_img);                    
0052 
0053 disp('step 5: simulate data for inhomogeneous medium');
0054 %
0055 % create an inhomogeneous image
0056 % A,B are Indices of the elements to represent the inhomogeneity
0057 %
0058 load( 'datacom.mat' ,'A','B')
0059 mat(A)= mat(A)+0.15;
0060 mat(B)= mat(B)-0.20;
0061 
0062 inhomg_img= eidors_obj('image', 'inhomogeneous image', ...
0063                        'elem_data', mat, ...
0064                        'fwd_model', demo_mdl );
0065 
0066 show_fem( inhomg_img );
0067 
0068 inhomg_data=fwd_solve( demo_mdl, inhomg_img);
0069 
0070 disp('step 6: add noise to simulated data');
0071 
0072 inhomg_data.meas = inhomg_data.meas + 1e-5*randn(size(inhomg_data.meas));
0073   homg_data.meas =   homg_data.meas + 1e-5*randn(size(  homg_data.meas));
0074 
0075 disp('step 7: create inverse model');
0076 
0077 % create an inv_model structure of name 'demo_inv'
0078 demo_inv.name= 'Nick Polydorides EIT inverse';
0079 demo_inv.solve= 'np_inv_solve';
0080 demo_inv.hyperparameter.value = 1e-3;
0081 demo_inv.R_prior= 'np_calc_image_prior';
0082 demo_inv.np_calc_image_prior.parameters= [3 1]; % see iso_f_smooth: deg=1, w=1
0083 demo_inv.jacobian_bkgnd.value= 1;
0084 demo_inv.reconst_type= 'difference';
0085 demo_inv.fwd_model= demo_mdl;
0086 demo_inv= eidors_obj('inv_model', demo_inv);
0087 
0088 disp('step 8: solve inverse model');
0089 
0090 demo_img= inv_solve( demo_inv, homg_data, inhomg_data);
0091 
0092 disp('step 9: display results');
0093 
0094 levels=[ 2.63 2.10 1.72 1.10 0.83 0.10 ];
0095 
0096 figure; image_levels( inhomg_img, levels );
0097 
0098 demo_img.name= 'Reconstructed conductivity distribution';
0099 demo_img.calc_colours.ref_level = 0;
0100 figure; image_levels( demo_img, levels );
0101 
0102 function [vtx,simp,bdy] = get_model_elems;
0103 %bdy : the boundary surfaces (triangles)
0104 %vtx : the vertices of the model (coordinates of the nodes)
0105 %simp: the simplices of the model (connectivity in tetrahedral)
0106 load('datareal.mat','vtx','simp');
0107 bdy= find_boundary( simp );
0108 
0109 
0110 
0111 function [gnd_ind, electrodes, perm_sym, elec, protocol, no_pl] = get_model_elecs;
0112 %elec : The electrodes matrix.
0113 %np_pl : Number of electrode planes (in planar arrangements)
0114 %protocol : Adjacent or Opposite or Customized.
0115 %zc : Contact impedances of the electrodes
0116 %perm_sym : Boolean entry for efficient forward computations
0117 %perm_sym='{n}';
0118 
0119 load('datareal.mat','gnd_ind','elec','zc','protocol','no_pl');
0120 
0121 for i=1:length(zc)
0122     electrodes(i).z_contact= zc(i);
0123     electrodes(i).nodes=     unique( elec(i,:) );
0124 end
0125 
0126 perm_sym='{n}';
0127 
0128 % get the current stimulation patterns
0129 function [stimulations] = get_model_stim( mdl );
0130 
0131 load('datareal.mat','protocol','no_pl','elec');
0132 [I,Ib] = set_3d_currents(protocol, ...
0133                          elec, ...
0134                          mdl.nodes, ...
0135                          mdl.gnd_node, ...
0136                          no_pl);
0137 
0138 % get the measurement patterns, only indH is used in this model
0139 %   here we only want to get the meas pattern from 'get_3d_meas',
0140 %   not the voltages, so we enter zeros
0141 [jnk,jnk,indH,indV,jnk] = get_3d_meas( ...
0142                   elec, mdl.nodes, ...
0143                   zeros(size(I)), ... % Vfwd
0144                   Ib, no_pl );
0145 n_elec= size(elec,1);
0146 n_meas= size(indH,1) / size(Ib,2);
0147 for i=1:size(Ib,2)
0148     stimulations(i).stimulation= 'Amp';
0149     stimulations(i).stim_pattern= Ib(:,i);
0150     idx= ( 1+ (i-1)*n_meas ):( i*n_meas );
0151     meas_pat = sparse( (1:n_meas)'*[1,1], ...
0152                        indH( idx, : ), ...
0153                        ones(n_meas,2)*[1,0;0,-1], ...
0154                        n_meas, n_elec );
0155     stimulations(i).meas_pattern= meas_pat;
0156 end
0157 
0158 
0159 
0160 
0161 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0162 % This is part of the EIDORS suite.
0163 % Copyright (c) N. Polydorides 2003
0164 % Copying permitted under terms of GNU GPL
0165 % See enclosed file gpl.html for details.
0166 % EIDORS 3D version 2.0
0167 % MATLAB version 6.1 R12
0168 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%