manchester_tomography

PURPOSE ^

Example to show reconstructions from

SYNOPSIS ^

function manchester_tomography( example_no)

DESCRIPTION ^

 Example to show reconstructions from
 the Process Tomography 
 group from U.Manchester, UK

 manchester_tomography( example_no)

 example_no == 1
    show reconstructions from 2planes / 2rods
    Reconstruction = inv_solve_trunc_iterative

 example_no == 2
    show reconstructions from 2planes / 2rods
    Reconstruction = np_inv_solve

 (C) 2005 by Stephen Murphy.  License: GPL version 2 or version 3
 $Id: manchester_tomography.m 3934 2013-04-20 15:46:00Z aadler $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 % Example to show reconstructions from
0002 % the Process Tomography
0003 % group from U.Manchester, UK
0004 %
0005 % manchester_tomography( example_no)
0006 %
0007 % example_no == 1
0008 %    show reconstructions from 2planes / 2rods
0009 %    Reconstruction = inv_solve_trunc_iterative
0010 %
0011 % example_no == 2
0012 %    show reconstructions from 2planes / 2rods
0013 %    Reconstruction = np_inv_solve
0014 %
0015 % (C) 2005 by Stephen Murphy.  License: GPL version 2 or version 3
0016 % $Id: manchester_tomography.m 3934 2013-04-20 15:46:00Z aadler $
0017 function manchester_tomography( example_no)
0018 
0019 switch example_no
0020     case 1,
0021         example_diff_morozov_reconst;
0022 
0023     case 2,
0024         example_diff_np_reconst;
0025 
0026     case 3,
0027         example_diff_tv_reconst_sim
0028 
0029     case 4,
0030         example_diff_tv_reconst
0031 
0032     case 10,
0033         example2
0034 
0035 end
0036 
0037 % output fwd_model, vinh, vhom
0038 function [fwd_m1,meas_2rod,meas_ref]= twoplane_mdl;
0039     load 2planes_Op_drive;
0040     fwd_m1=eidors_obj('fwd_model',mdl_2plane);
0041     fwd_m1.np_fwd_solve.perm_sym = '{y}';
0042     fwd_m1.solve=    'np_fwd_solve';
0043     fwd_m1.jacobian= 'np_calc_jacobian';
0044     fwd_m1.system_mat= 'np_calc_system_mat';
0045     meas_ref.type='data';
0046     meas_2rod.type='data';
0047 
0048 function [DS_coarse, DS_dense]= coarse_dense_mdl
0049     load dual_mesh;
0050     DS_coarse= set_fwd_model(vtx_coarse,simp_coarse,[], ...
0051               elec_coarse,zc,gnd_ind_coarse,Ib_coarse, [], []);
0052     DS_dense= set_fwd_model(vtx_dense,simp_dense,[], ...
0053               elec_dense,zc,gnd_ind_dense,Ib_dense, [], []);
0054     stim= mk_stim_patterns(16,1,'{ad}','{ad}',[],1);
0055     M_coarse.stimulation= stim;
0056     M_dense.stimulation= stim;
0057 
0058 
0059 function example_diff_morozov_reconst
0060     [fwd_m1,vi,vh]= twoplane_mdl;
0061 
0062     imdl= eidors_obj('inv_model','Morozov mdl');
0063     imdl.solve= 'np_inv_solve';
0064     imdl.solve= 'inv_solve_trunc_iterative';
0065     imdl.R_prior= 'np_calc_image_prior';
0066     imdl.jacobian_bkgnd.value = .01;
0067     imdl.np_calc_image_prior.parameters= [3 1];
0068     imdl.reconst_type= 'difference';
0069     imdl.fwd_model= fwd_m1;
0070 
0071     imr= inv_solve(imdl, vh, vi);
0072     show_slices(imr, linspace(.01,.09,4)'*[inf,inf,1])
0073 
0074 function example_diff_np_reconst
0075     [fwd_m1,vi,vh]= twoplane_mdl;
0076 
0077     imdl= eidors_obj('inv_model','NP mdl');
0078     imdl.solve= 'np_inv_solve';
0079     imdl.R_prior= 'np_calc_image_prior';
0080     imdl.hyperparameter.value = 1e-2;
0081 
0082     imdl.jacobian_bkgnd.value = .01;
0083     imdl.np_calc_image_prior.parameters= [3 1];
0084     imdl.reconst_type= 'difference';
0085     imdl.fwd_model= fwd_m1;
0086 
0087     imr= inv_solve(imdl, vh, vi);
0088     show_slices(imr, linspace(.01,.09,4)'*[inf,inf,1])
0089 
0090 function example_diff_tv_reconst_sim
0091     [fwd_m1,vi,vh]= twoplane_mdl;
0092 
0093     imdl= eidors_obj('inv_model','TV mdl');
0094     imdl.solve= 'inv_solve_TV_pdipm';
0095     imdl.R_prior= 'prior_TV';
0096     imdl.parameters.max_iterations= 5;
0097     imdl.hyperparameter.value = 1e-4;
0098 
0099     imdl.jacobian_bkgnd.value = 1;
0100     imdl.reconst_type= 'difference';
0101     imdl.fwd_model= fwd_m1;
0102     [vi,vh] = sim_inhomg(fwd_m1);
0103 
0104     imr= inv_solve(imdl, vh, vi);
0105     show_slices(imr, linspace(.01,.09,4)'*[inf,inf,1])
0106 
0107 function example_diff_tv_reconst
0108     [fwd_m1,vi,vh]= twoplane_mdl;
0109 
0110     imdl= eidors_obj('inv_model','TV mdl');
0111     imdl.solve= 'inv_solve_TV_pdipm';
0112     imdl.R_prior= 'prior_TV';
0113     imdl.parameters.max_iterations= 3;
0114     imdl.hyperparameter.value = 1e-4;
0115 
0116     imdl.jacobian_bkgnd.value = .01;
0117     imdl.reconst_type= 'difference';
0118     imdl.fwd_model= fwd_m1;
0119     [vi,vh] = sim_inhomg(fwd_m1);
0120 
0121     imr= inv_solve(imdl, vh, vi);
0122     show_slices(imr, linspace(.01,.09,4)'*[inf,inf,1])
0123 
0124 
0125 
0126 function example_diff_sim_reconst
0127     load 2planes_Op_drive;
0128     fwd_m1=eidors_obj('fwd_model',mdl_2plane);
0129     fwd_m1.np_fwd_solve.perm_sym = '{y}';
0130     fwd_m1.solve=    'np_fwd_solve';
0131     fwd_m1.jacobian= 'np_calc_jacobian';
0132     fwd_m1.system_mat= 'np_calc_system_mat';
0133 
0134 function [vi,vh] = sim_inhomg( mdl);
0135     n= size(mdl.elems,1);
0136     img2=eidors_obj('image','homog');
0137     img2.elem_data= ones(n,1);
0138     img2.fwd_model= mdl;
0139     vh= fwd_solve(img2);
0140 
0141     cc= center_of_simps( mdl);
0142     [jnk,mat]= elems_in_cylinder(cc,[0,.03,.02],.01,ones(n,1),1.2);
0143     [jnk,mat]= elems_in_cylinder(cc,[0,-.03,-.02],.01,mat,0.8);
0144     img2.elem_data= mat;
0145 
0146 %   show_fem(img2);
0147     vi= fwd_solve(img2);
0148 
0149 function mdl_z= zig_zag_from_2rings(mdl_2r)
0150     mdl_z= mdl_2r;
0151     renumber=[1:2:16;18:2:32];
0152     renumber=renumber(:);
0153     mdl_z.electrode= mdl_2r.electrode(renumber);
0154     mdl_z.stimulation= mk_stim_patterns(16,1,'{ad}','{ad}',[],1);
0155 
0156 
0157 function example2( M_coarse, M_dense)
0158     [M_coarse, M_dense]= coarse_dense_mdl;
0159 
0160 imdl= mk_common_model('b2c',16);
0161 imdl.hyperparameter.value= 1e-2;
0162 imgr= inv_solve(imdl, vh,vi);
0163 show_fem(imgr);
0164 
0165 
0166 imdl2= eidors_obj('inv_model', 'DS dual mesh');
0167 imdl2.solve= 'inv_solve_dual_mesh';
0168 imdl2.hyperparameter.value = 1e-4;
0169 imdl2.R_prior= 'np_calc_image_prior';
0170 imdl2.np_calc_image_prior.parameters= [3 1]; % see iso_f_smooth: deg=1, w=1
0171 imdl2.jacobian_bkgnd.value= 1;
0172 imdl2.reconst_type= 'static';
0173 imdl2.fwd_model= M_dense;
0174 imdl2.parameters.max_iterations= 5;
0175 
0176 imdl2.inv_solve_dual_mesh.coarse_mdl= M_coarse;
0177 imdl2.inv_solve_dual_mesh.mapper_func= 'edge_refined_elem_mapper';
0178 
0179 imgr= inv_solve(imdl2, vi);

Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005