object_in_tank_2d

 2D demo example for reconstruction of object floating inside tank with
 ambient fluid. Only internat structure of the object is reconstruted.
 External shape od object, and parameters of fluid is treated as a known.
 Dual meshing for inverse problem has been deployed.
 
 (C) 2010 Bartosz Sawicki. License: GPL v.3

0001 % 2D demo example for reconstruction of object floating inside tank with
0002 % ambient fluid. Only internat structure of the object is reconstruted.
0003 % External shape od object, and parameters of fluid is treated as a known.
0004 % Dual meshing for inverse problem has been deployed.
0005 %
0006 % (C) 2010 Bartosz Sawicki. License: GPL v.3
0007 
0008 %% Create object model for reconstruction
0009 
0010 object_center = [0.0, 0.1];
0011 object_radius = 0.6;
0012 
0013 mdl = create_circle_mesh_gmsh('circle', object_center, object_radius, 0.04 );
0014 object_mdl = eidors_obj('fwd_model', mdl);
0015 
0016 %% Create fine 'tank with object' model for forward problem
0017 
0018 electrodes_per_plane = 16;
0019 tank_radius = 1.0;
0020 
0021 extra={'object','solid object = cylinder(0.0,0.,0;0,0,1;0.6) and plane(0,0,0;0,0,-1) and plane(0,0,0.2;0,0,1) -maxh=0.05;'};
0022 [tank_mdl, mat_indices] = ng_mk_cyl_models([0,tank_radius,0.03],[electrodes_per_plane],[0.1,0,0.01], extra);
0023 
0024 fluid_indices=mat_indices{1};
0025 object_indices=mat_indices{2};
0026 
0027 options = {'no_meas_current','no_rotate_meas'};
0028 [st, els]= mk_stim_patterns(electrodes_per_plane, 1, '{ad}','{ad}', options, 10);
0029 
0030 tank_mdl.stimulation= st;
0031 tank_mdl.meas_select= els;
0032 tank_mdl.solve=      'fwd_solve_1st_order';
0033 tank_mdl.system_mat= 'bs_calc_system_mat';
0034 tank_mdl.jacobian=   'jacobian_adjoint';
0035 tank_mdl.normalize_measurements= 0;
0036 tank_mdl.np_fwd_solve.perm_sym= '{n}';
0037 
0038 %% Solve model without inclusion
0039 disp('Model without inclusion');
0040 
0041 %Initialize
0042 mat = ones(size(tank_mdl.elems,1),1);
0043 %Fluid conductivity
0044 mat(fluid_indices) = 0.75;
0045 %Object avarage conductivity
0046 mat(object_indices) = 1.0;
0047 
0048 homg_img= eidors_obj('image', 'homogeneous image', ...
0049                      'elem_data', mat, ...
0050                      'fwd_model', tank_mdl );
0051 
0052 homg_data=fwd_solve( tank_mdl, homg_img);
0053 
0054 %% Create inclusion and solve model
0055 disp('Model with inclusion');
0056 
0057 inclusion_center = [0.3, 0.3];
0058 inclusion_radius = 0.05;
0059 inclusion_conductivity = 1.1;
0060 
0061 inhomg_img = create_inclusion(homg_img, inclusion_center, inclusion_radius, ...
0062     inclusion_conductivity);
0063 
0064 inhomg_img.name = 'Inhomogeneous model';
0065 
0066 figure; show_fem( inhomg_img );
0067 
0068 inhomg_data = fwd_solve( tank_mdl, inhomg_img);
0069 
0070 %% Add some noise
0071 disp('Add noise');
0072 
0073 SNR = -30; %dB
0074 
0075 hdata = homg_data.meas;
0076 idata = inhomg_data.meas;
0077 noise = 10^(SNR/20)*std(idata-hdata)*randn(size(idata));
0078 idata = idata + noise;
0079 inhomg_data.meas = idata ;
0080 
0081 %% Create inverse model
0082 disp('Inverse model')
0083 
0084 inv2d.name=  'Object in tank 2D EIT inverse';
0085 inv2d.fwd_model= tank_mdl;
0086 inv2d.rec_model= object_mdl;
0087 
0088 disp('Calculating coarse2fine mapping ...');
0089 inv2d.fwd_model.coarse2fine = ...
0090        mk_coarse_fine_mapping( tank_mdl, object_mdl);
0091 disp('   ... done');
0092 
0093 % Absolute reconstruction
0094 inv2d.reconst_type= 'static';
0095 inv2d.solve= @bs_nonlinearGN;
0096 inv2d.parameters.term_tolerance= 1e-8;
0097 inv2d.parameters.max_iterations= 5;
0098 inv2d.nonlinearGN.init_backgnd= 0.00;
0099 %
0100 inv2d.RtR_prior= @prior_tikhonov;
0101 inv2d.hyperparameter.value= 1e-3;
0102 %
0103 inv2d.fwd_model.background= homg_img.elem_data;
0104 
0105 inv2d= eidors_obj('inv_model', inv2d);
0106 recon_img= inv_solve(inv2d, inhomg_data);
0107 
0108 figure; show_fem( recon_img );