demo_3d_simdata

 How to make simulation data using EIDORS3D

0001 % How to make simulation data using EIDORS3D
0002 
0003 % (C) 2005 Nick Polydorides + Andy Adler. License: GPL version 2 or version 3
0004 % $Id: demo_3d_simdata.m 3776 2013-03-18 18:10:07Z aadler $
0005 
0006 % STIMULATION PATTERN
0007 n_elec= 16;
0008 n_rings= 1;
0009 %options = {'no_meas_current','rotate_meas'};
0010  options = {'no_meas_current','no_rotate_meas'};
0011 stimulation= mk_stim_patterns(n_elec, n_rings, '{ad}','{ad}', ...
0012                             options, 10);
0013 %
0014 % Example 1: Create 16 electrode 3D model
0015 %
0016 disp('STEP 1: Model simultion 3D');
0017 
0018 % get parameters for model from mk_circ_tank
0019 % param= mk_circ_tank(rings, levels, n_elec, n_planes )
0020 levels= [-.5:.1:.5];
0021 e_levels= [4, 8];
0022 % params= mk_circ_tank( 8, levels, n_elec );
0023   params= mk_circ_tank( 8, levels, { 'zigzag', n_elec, e_levels } );
0024 % params= mk_circ_tank(12, levels, { 'zigzag', n_elec, [3,5,7] , ...
0025 %                                    'planes', n_elec, 2} );
0026 
0027 params.stimulation= stimulation;
0028 params.solve=      @fwd_solve_1st_order;
0029 %params.solve=      @eidors_default;
0030 params.system_mat= @system_mat_1st_order;
0031 params.jacobian=   @jacobian_adjoint;
0032 mdl_3d = eidors_obj('fwd_model', params);
0033 
0034 
0035 disp('STEP 1A: simultion 3D - homogeneous');
0036 % create homogeneous image + simulate data
0037 cond= ones( size(mdl_3d.elems,1) ,1);
0038 homg_img= eidors_obj('image', 'homogeneous image', ...
0039                      'elem_data', cond, ...
0040                      'fwd_model', mdl_3d );
0041 homg_data=fwd_solve( homg_img);
0042 
0043 disp('STEP 1B: simultion 3D - inhomogeneous');
0044 
0045 % create inhomogeneous image + simulate data
0046 inhv= [38,50,51,66,67,83];
0047 for inhlev= (e_levels(1)-1)*3 + [-3:2];
0048      cond(256*inhlev+inhv) =2;
0049 %    cond(64*inhlev+[5,9,10,17,18,26]) =2;
0050 end
0051 
0052 inh_img= eidors_obj('image', 'inhomogeneous image', ...
0053                      'elem_data', cond, ...
0054                      'fwd_model', mdl_3d );
0055 inh_data=fwd_solve( inh_img);
0056 show_fem( inh_img);
0057 disp([inh_img.name, '. Press a key']); pause;
0058 
0059 % Add 10% noise
0060 sig= std( inh_data.meas - homg_data.meas );
0061 inh_data.meas= inh_data.meas + 0.10 * sig* randn( size(inh_data.meas) );
0062 
0063 %
0064 % Step 2: Reconstruction in 2D
0065 %
0066 params= mk_circ_tank(8, [], n_elec);
0067 
0068 params.stimulation= stimulation;
0069 params.solve=      'fwd_solve_1st_order';
0070 params.system_mat= 'system_mat_1st_order';
0071 params.jacobian=   'jacobian_adjoint';
0072 mdl_2d_2 = eidors_obj('fwd_model', params);
0073 
0074 inv2d.name= 'EIT inverse';
0075 inv2d.solve=       'inv_solve_diff_GN_one_step';
0076 %inv2d.hyperparameter.func = 'select_noise_figure';
0077 %inv2d.hyperparameter.noise_figure= 2;
0078 %inv2d.hyperparameter.tgt_elems= 1:4;
0079  inv2d.hyperparameter.value = 1e-1;
0080  inv2d.RtR_prior= 'prior_laplace';
0081 %inv2d.RtR_prior= 'prior_gaussian_HPF';
0082 inv2d.jacobian_bkgnd.value= 1;
0083 inv2d.reconst_type= 'difference';
0084 inv2d.fwd_model= mdl_2d_2;
0085 inv2d= eidors_obj('inv_model', inv2d);
0086 
0087 img2= inv_solve( inv2d, inh_data, homg_data);
0088 img2.name= '2D inverse solution';
0089 if ~exist('OCTAVE_VERSION')
0090    show_fem(img2);
0091 else
0092    show_slices(img2);
0093 end
0094 disp([img2.name, '. Press a key']); pause;
0095 
0096 %
0097 % Step 2: Reconstruction in 3D (using np_2003 code) and point
0098 %          electrode models with 'zigzag' electrodes
0099 %
0100 disp('STEP 2: Reconstruction 3D');
0101 clear inv3d;
0102 
0103  levels= [-.4:.2:.4];
0104  params= mk_circ_tank( 8, levels, { 'zigzag', n_elec, [2,4] } );
0105 %params= mk_circ_tank( 8, levels, { 'zigzag', n_elec, e_levels } );
0106 %params= mk_circ_tank( 4, levels, { 'zigzag', n_elec, e_levels } );
0107 %params= mk_circ_tank( 4, levels, n_elec );
0108 params.stimulation= stimulation;
0109 params.solve=      @fwd_solve_1st_order;
0110 params.system_mat= @system_mat_1st_order;
0111 params.jacobian=   @jacobian_adjoint;
0112 params.misc.perm_sym= '{n}';
0113 fm3d = eidors_obj('fwd_model', params);
0114 
0115 inv3d.name=  'EIT inverse: 3D';
0116 inv3d.solve= @inv_solve_diff_GN_one_step;
0117 inv3d.hyperparameter.value = 1e-2;
0118 inv3d.jacobian_bkgnd.value= 1;
0119 inv3d.RtR_prior= 'prior_laplace';
0120 inv3d.reconst_type= 'difference';
0121 inv3d.fwd_model= fm3d;
0122 inv3d= eidors_obj('inv_model', inv3d);
0123 
0124  img3= inv_solve( inv3d, inh_data, homg_data);
0125  img3.name= '3D inverse solution';
0126 
0127 if ~exist('OCTAVE_VERSION')
0128    show_fem(img3);
0129 else
0130    show_slices(img3, [-.35:.2:.4]'*[inf,inf,1])
0131 end
0132