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 6253 2022-03-26 16:30:22Z 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 homg_img= mk_image(mdl_3d, 1);
0038 homg_data=fwd_solve( homg_img);
0039 
0040 disp('STEP 1B: simultion 3D - inhomogeneous');
0041 
0042 % create inhomogeneous image + simulate data
0043 inh_img = homg_img;
0044 inhv= [38,50,51,66,67,83];
0045 for inhlev= (e_levels(1)-1)*3 + [-3:2];
0046     inh_img.elem_data(256*inhlev+inhv) =2;
0047 end
0048 inh_data=fwd_solve( inh_img);
0049 subplot(221);show_fem( inh_img);
0050 
0051 % Add noise SNR=20
0052 sig= std( inh_data.meas - homg_data.meas );
0053 inh_data = add_noise(20, inh_data, homg_data);
0054 
0055 %
0056 % Step 2: Reconstruction in 2D
0057 %
0058 params= mk_circ_tank(8, [], n_elec);
0059 
0060 params.stimulation= stimulation;
0061 params.solve=      'fwd_solve_1st_order';
0062 params.system_mat= 'system_mat_1st_order';
0063 params.jacobian=   'jacobian_adjoint';
0064 mdl_2d_2 = eidors_obj('fwd_model', params);
0065 
0066 inv2d.name= 'EIT inverse';
0067 inv2d.solve=       'inv_solve_diff_GN_one_step';
0068 %inv2d.hyperparameter.func = 'select_noise_figure';
0069 %inv2d.hyperparameter.noise_figure= 2;
0070 %inv2d.hyperparameter.tgt_elems= 1:4;
0071  inv2d.hyperparameter.value = 1e-1;
0072  %inv2d.RtR_prior= 'prior_TV';
0073 %inv2d.R_prior = 'prior_TV';
0074  inv2d.RtR_prior= 'prior_gaussian_HPF';
0075 inv2d.jacobian_bkgnd.value= 1;
0076 inv2d.reconst_type= 'difference';
0077 inv2d.fwd_model= mdl_2d_2;
0078 inv2d= eidors_obj('inv_model', inv2d);
0079 
0080 img2= inv_solve( inv2d, homg_data, inh_data);
0081 img2.name= '2D inverse solution';
0082 subplot(223);   show_slices(img2);
0083 
0084 %
0085 % Step 2: Reconstruction in 3D (using np_2003 code) and point
0086 %          electrode models with 'zigzag' electrodes
0087 %
0088 disp('STEP 2: Reconstruction 3D');
0089 clear inv3d;
0090 
0091  levels= [-.4:.2:.4];
0092  params= mk_circ_tank( 8, levels, { 'zigzag', n_elec, [2,4] } );
0093 %params= mk_circ_tank( 8, levels, { 'zigzag', n_elec, e_levels } );
0094 %params= mk_circ_tank( 4, levels, { 'zigzag', n_elec, e_levels } );
0095 %params= mk_circ_tank( 4, levels, n_elec );
0096 params.stimulation= stimulation;
0097 params.solve=      @fwd_solve_1st_order;
0098 params.system_mat= @system_mat_1st_order;
0099 params.jacobian=   @jacobian_adjoint;
0100 params.misc.perm_sym= '{n}';
0101 fm3d = eidors_obj('fwd_model', params);
0102 
0103 inv3d.name=  'EIT inverse: 3D';
0104 inv3d.solve= @inv_solve_diff_GN_one_step;
0105 inv3d.hyperparameter.value = 1e-2;
0106 inv3d.jacobian_bkgnd.value= 1;
0107 %inv3d.RtR_prior= 'prior_TV';
0108 inv3d.R_prior = 'prior_TV';
0109 inv3d.reconst_type= 'difference';
0110 inv3d.fwd_model= fm3d;
0111 inv3d= eidors_obj('inv_model', inv3d);
0112 
0113  img3= inv_solve( inv3d, homg_data, inh_data);
0114  img3.name= '3D inverse solution';
0115 
0116 subplot(122)
0117 level(:,3) = [-.35:.2:.4]'; level(:,1:2) = Inf;
0118 show_slices(img3, level);
0119