cheating_2d

 code to simulate inverse crimes in EIT

0001 % code to simulate inverse crimes in EIT
0002 
0003 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0004 % $Id: cheating_2d.m 5112 2015-06-14 13:00:41Z aadler $
0005 
0006 %TODO: calculate how well data matches priors
0007 function out=cheating_2d( figno, rand_seed )
0008    [vis,vhs,s_mdl]= small_2d_mdl;
0009 
0010    if nargin<2; rand_seed= []; end
0011 
0012    il_g = make_inv_model( 12 ); %large model
0013 
0014    if nargin==0; figno= {'1','2a','2b','3a','3b','4'}; end 
0015    if ischar(figno); figno= {figno}; end
0016 
0017    for idx= 1:length(figno)
0018        if idx~=1; pause; end
0019 
0020        if     strcmp( figno{idx}, '1' )
0021            approach1(vis, vhs, s_mdl, il_g, rand_seed)
0022        elseif strcmp( figno{idx}, '2a' )
0023            approach2a(vis, vhs, s_mdl, il_g)
0024        elseif strcmp( figno{idx}, '2b' )
0025            approach2b(vis, vhs, s_mdl, il_g)
0026        elseif strcmp( figno{idx}, '3a' )
0027            approach3a(vis, vhs, s_mdl, il_g)
0028        elseif strcmp( figno{idx}, '3b' )
0029            approach3b(vis, vhs, s_mdl, il_g)
0030        elseif strcmp( figno{idx}, '4' )
0031            approach4(vis, vhs, s_mdl, il_g, rand_seed)
0032        elseif strcmp( figno{idx}, '4b' )
0033        mdl= mk_common_model('b2c');
0034            if ~isempty(rand_seed); randn('state', rand_seed(i)); end
0035            def_mdl = eidors_obj('fwd_model', angl_deform(mdl.fwd_model, .02));
0036        show_fem(def_mdl);
0037        end
0038    end
0039 %
0040 % APPROACH 1
0041 %
0042 function approach1(vis, vhs, s_mdl, il_g, rand_seed)
0043    disp('Approach #1: reconstruct with noise');
0044 
0045    num_tries=6;
0046    levels= [];
0047    if ~isempty(rand_seed);
0048       num_tries= length(rand_seed);
0049       levels= [0,0,0,1,1];
0050    end
0051 
0052    vi_n(1:num_tries) = vis;
0053    for i= 1:num_tries % stupid matlab doesn't allow easy vectorization
0054       if ~isempty(rand_seed); randn('state', rand_seed(i)); end
0055       noise = .0002*randn(size(vis.meas));
0056       vi_n(i).meas = vi_n(i).meas + noise;
0057    end
0058    show_slices( inv_solve( il_g, vhs, vi_n ), levels);
0059 
0060 
0061 %
0062 % APPROACH 2A
0063 %
0064 function approach2a(vis, vhs, s_mdl, il_g)
0065    levels= [0,0,0,1,1];
0066    disp('Approach #2: reconstruct with tikhonov cheat - with inv crime');
0067    pp= small_face;
0068    % homog (normal) model
0069    RtR_prior= @cheat_tikhonov;
0070    param_vals.cheat_elements= [];
0071    param_vals.cheat_weight = 0.5;
0072    is_n = make_inv_model( 8 , RtR_prior, 'cheat_tikhonov', param_vals ); 
0073    % sad model
0074    param_vals.cheat_elements= [pp.eyes, pp.sad];
0075    is_s = make_inv_model( 8 , RtR_prior, 'cheat_tikhonov', param_vals ); 
0076    % happy model
0077    param_vals.cheat_elements= [pp.eyes, pp.smile];
0078    is_h = make_inv_model( 8 , RtR_prior, 'cheat_tikhonov', param_vals ); 
0079    % happy/sad (medium) model
0080    param_vals.cheat_elements= [pp.eyes, pp.rsmile, pp.lsad];
0081    is_m = make_inv_model( 8 , RtR_prior, 'cheat_tikhonov', param_vals ); 
0082 
0083    show_slices( [ inv_solve( is_n, vhs, vis ), ... 
0084                   inv_solve( is_s, vhs, vis ), ... 
0085                   inv_solve( is_h, vhs, vis ), ... 
0086                   inv_solve( is_m, vhs, vis ) ], levels );
0087 
0088 
0089 
0090 %
0091 % APPROACH 2B
0092 %
0093 function approach2b(vis, vhs, s_mdl, il_g)
0094    disp('Approach #2B: reconstruct with tikhonov cheat - without inv crime');
0095    levels= [0,0,0,1,1];
0096    pp= large_face;
0097    % homog (normal) model
0098    RtR_prior= @cheat_tikhonov;
0099    param_vals.cheat_elements= [];
0100    param_vals.cheat_weight = 0.5;
0101    il_n = make_inv_model(12 , RtR_prior, 'cheat_tikhonov', param_vals ); 
0102    % sad model
0103    param_vals.cheat_elements= pp.sad;
0104    il_s = make_inv_model(12 , RtR_prior, 'cheat_tikhonov', param_vals ); 
0105    % happy model
0106    param_vals.cheat_elements= pp.happy;
0107    il_h = make_inv_model(12 , RtR_prior, 'cheat_tikhonov', param_vals ); 
0108    % happy/sad (medium) model
0109    param_vals.cheat_elements= pp.halfy;
0110    il_m = make_inv_model(12 , RtR_prior, 'cheat_tikhonov', param_vals ); 
0111 
0112    show_slices( [ inv_solve( il_n, vhs, vis ), ... 
0113                   inv_solve( il_s, vhs, vis ), ... 
0114                   inv_solve( il_h, vhs, vis ), ... 
0115                   inv_solve( il_m, vhs, vis ) ], levels );
0116 
0117 %
0118 % APPROACH 3A
0119 %
0120 function approach3a(vis, vhs, s_mdl, il_g)
0121    disp('Approach #3: reconstruct with Laplace filter cheat - with inv crime');
0122    levels= [0,0,0,1,1];
0123    pp= small_face;
0124    % homog (normal) model
0125    RtR_prior= @cheat_laplace;
0126    param_vals.cheat_elements= [];
0127    param_vals.cheat_weight = 0.2;
0128    is_n = make_inv_model( 8 , RtR_prior, 'cheat_laplace', param_vals ); 
0129    % sad model
0130    param_vals.cheat_elements= [pp.eyes, pp.sad];
0131    is_s = make_inv_model( 8 , RtR_prior, 'cheat_laplace', param_vals ); 
0132    % happy model
0133    param_vals.cheat_elements= [pp.eyes, pp.smile];
0134    is_h = make_inv_model( 8 , RtR_prior, 'cheat_laplace', param_vals ); 
0135    % happy/sad (medium) model
0136    param_vals.cheat_elements= [pp.eyes, pp.rsmile, pp.lsad];
0137    is_m = make_inv_model( 8 , RtR_prior, 'cheat_laplace', param_vals ); 
0138 
0139    show_slices( [ inv_solve( is_n, vhs, vis ), ... 
0140                   inv_solve( is_s, vhs, vis ), ... 
0141                   inv_solve( is_h, vhs, vis ), ... 
0142                   inv_solve( is_m, vhs, vis ) ], levels );
0143 
0144 
0145 %
0146 % APPROACH 3B
0147 %
0148 function approach3b(vis, vhs, s_mdl, il_g)
0149    disp('Approach #3B: reconstruct with Laplace cheat - without inv crime');
0150    levels= [0,0,0,1,1];
0151    pp= large_face;
0152    % homog (normal) model
0153    RtR_prior= @cheat_laplace;
0154    param_vals.cheat_elements= [];
0155    param_vals.cheat_weight = 0.2;
0156    il_n = make_inv_model(12 , RtR_prior, 'cheat_laplace', param_vals ); 
0157    % sad model
0158    param_vals.cheat_elements= pp.sad;
0159    il_s = make_inv_model(12 , RtR_prior, 'cheat_laplace', param_vals ); 
0160    % happy model
0161    param_vals.cheat_elements= pp.happy;
0162    il_h = make_inv_model(12 , RtR_prior, 'cheat_laplace', param_vals ); 
0163    % happy/sad (medium) model
0164    param_vals.cheat_elements= pp.halfy;
0165    il_m = make_inv_model(12 , RtR_prior, 'cheat_laplace', param_vals ); 
0166 
0167    show_slices( [ inv_solve( il_n, vhs, vis ), ... 
0168                   inv_solve( il_s, vhs, vis ), ... 
0169                   inv_solve( il_h, vhs, vis ), ... 
0170                   inv_solve( il_m, vhs, vis ) ], levels );
0171 
0172 
0173 %
0174 % APPROACH 4
0175 %
0176 function approach4(vis, vhs, s_mdl, il_g, rand_seed)
0177    disp('Approach #4: deform the model');
0178    num_tries=6;
0179    levels= [];
0180    if ~isempty(rand_seed);
0181       num_tries= length(rand_seed);
0182       levels= [0,0,0,1,1];
0183    end
0184 
0185    params= mk_circ_tank(8, [], 16 ); 
0186    params.stimulation= mk_stim_patterns(16, 1, '{ad}','{ad}', ...
0187                          {'no_meas_current','no_rotate_meas'}, 1);
0188    params.solve=      'fwd_solve_1st_order';
0189    params.system_mat= 'system_mat_1st_order';
0190 
0191    mat= ones( size(params.elems,1), 1);
0192    pp= small_face;
0193    mat(pp.eyes)= 2;
0194    mat(pp.sad)=1.5;
0195 
0196    def_amount= .0035;
0197 
0198    for i= 1:num_tries % stupid matlab doesn't allow easy vectorization
0199       if ~isempty(rand_seed); randn('state', rand_seed(i)); end
0200 
0201       def_mdl = eidors_obj('fwd_model', angl_deform(params, def_amount ) );
0202       vi_m(i)= fwd_solve( eidors_obj('image','name',  ...
0203                      'elem_data', mat, 'fwd_model', def_mdl ));
0204    end
0205    show_slices( inv_solve( il_g, vhs, vi_m ), levels);
0206 
0207 
0208 
0209    
0210 return;
0211 
0212 function p= small_face;
0213    p.leye=   [78,97,98,117,118,141];
0214    p.reye=   [66,82,83,102,103,123];
0215    p.rsmile= [40:41, 53:55, 69:71, 86:88];
0216    p.lsmile= [43:44, 57:59, 73:75, 91:93];
0217    p.lsad =  [31,43,58,57,74,73,92,93,113,112,135];
0218    p.sad =   [28,31,40,43,58,57,53,54,74,73,92,93,113, ...
0219               112,135,69,87,70,107,88,108,129];
0220    p.eyes= [p.leye,p.reye];
0221    p.smile= [p.rsmile, p.lsmile];
0222 
0223 function p=large_face
0224 p.sad=  [ 53; 57; 69; 73; 86; 87; 91; 92; 106; 107; 111; 112; 127; 128;
0225          129; 133; 134; 135; 147; 151; 152; 153; 157; 158; 159; 165;
0226          171; 172; 177; 178; 179; 184; 185; 186; 192; 193; 199; 200;
0227          205; 206; 207; 212; 213; 214; 220; 221; 227; 228; 229; 235;
0228          236; 243; 244; 251; 252; 253; 259; 260; 261; 267; 268; 275;
0229          276; 283; 284; 285; 292; 293; 301; 310; 319; 320]; 
0230 
0231 p.happy= [ ...
0232         69; 70; 71; 73; 74; 75; 86; 87; 88; 89; 91; 92; 93; 94; 106;
0233        107; 108; 109; 111; 112; 113; 114; 127; 128; 129; 130; 131; 133;
0234        134; 135; 136; 137; 147; 151; 152; 153; 154; 155; 157; 158; 159;
0235        160; 161; 165; 171; 172; 176; 177; 178; 179; 180; 183; 184; 185;
0236        186; 187; 192; 193; 199; 200; 220; 221; 227; 228; 229; 251; 252;
0237        253; 259; 260; 261; 283; 284; 285; 292; 293; 319; 320]; 
0238 
0239 p.halfy= [ ...
0240         57; 69; 70; 71; 73; 86; 87; 88; 89; 91; 92; 106; 107; 108; 109;
0241        111; 112; 127; 128; 129; 130; 131; 133; 134; 135; 147; 151; 152;
0242        153; 154; 155; 157; 158; 159; 165; 171; 172; 176; 177; 178; 179;
0243        180; 184; 185; 186; 192; 193; 199; 200; 212; 213; 214; 220; 221;
0244        227; 228; 229; 243; 244; 251; 252; 253; 259; 260; 261; 275; 276;
0245        283; 284; 285; 292; 293; 310; 319; 320];
0246 
0247 % simulate 'sad' data for small model
0248 function [vis,vhs,mdl]= small_2d_mdl
0249    params= mk_circ_tank(8, [], 16 ); 
0250    params.stimulation= mk_stim_patterns(16, 1, '{ad}','{ad}', ...
0251                          {'no_meas_current','no_rotate_meas'}, 1);
0252    params.solve=      'fwd_solve_1st_order';
0253    params.system_mat= 'system_mat_1st_order';
0254    mdl= eidors_obj('fwd_model', params);
0255 
0256    mat= ones( size(mdl.elems,1), 1);
0257 
0258    % homogeneous data
0259    vhs= fwd_solve( eidors_obj('image','name',  ...
0260                      'elem_data', mat, 'fwd_model', mdl ));
0261 
0262    % inhomogeneous data for sad face model
0263    pp= small_face;
0264    mat(pp.eyes)= 2;
0265    mat(pp.sad)=1.5;
0266 
0267    vis= fwd_solve( eidors_obj('image','name',  ...
0268                      'elem_data', mat, 'fwd_model', mdl ));
0269 
0270 % create inv_model, and specify img_prior if required
0271 function i_mdl= make_inv_model( n_rings, img_prior, param_name, param_vals );
0272    params= mk_circ_tank(n_rings, [], 16 ); 
0273    params.stimulation= mk_stim_patterns(16, 1, '{ad}','{ad}', ...
0274                          {'no_meas_current','no_rotate_meas'}, 1);
0275    params.solve=      'fwd_solve_1st_order';
0276    params.system_mat= 'system_mat_1st_order';
0277    params.jacobian  = 'jacobian_adjoint';
0278 %  params.normalize_measurements  = 1; TODO: we have a bug here
0279    l_mdl= eidors_obj('fwd_model', params);
0280 
0281 % create inverse model
0282    hparam.value = 3e-4;
0283   %hparam.func = 'select_noise_figure';
0284   %hparam.noise_figure= 1;
0285   %hparam.tgt_elems= 1:4;
0286 
0287    if nargin < 2
0288       img_prior = 'prior_tikhonov';
0289       param_name= 'jnk___'; param_vals= 'jnk___';
0290    end
0291 
0292    i_mdl= eidors_obj( ...
0293           'inv_model', '2D inverse', ...
0294           'hyperparameter', hparam, 'RtR_prior', img_prior, ...
0295           param_name, param_vals, ...
0296           'reconst_type', 'difference', ...
0297           'fwd_model', l_mdl, 'solve', 'inv_solve_diff_GN_one_step' );
0298    i_mdl.jacobian_bkgnd.value= 1;
0299 
0300 
0301 function reconst_with_noise(i_mdl, vis, vhs, num_tries)
0302     noise = .0002*randn(size(vis.meas,1),num_tries);
0303 
0304     vi_n(1:num_tries) = vis;
0305     for i= 1:num_tries % stupid matlab doesn't allow easy vectorization
0306        vi_n(i).meas = vi_n(i).meas + noise(:,i);
0307     end
0308     show_slices( inv_solve( i_mdl, vhd, vi_n ));
0309 
0310 % Image prior with Tikhonov + cheating elements
0311 function Reg= cheat_tikhonov( inv_model )
0312 % Reg= cheat_tikhonov( inv_model )
0313 % Reg        => output regularization term
0314 % Parameters:
0315 %   elems    = inv_model.RtR_prior.cheat_elements;
0316 %            => elements weights to modify
0317 %   weight   = inv_model.RtR_prior.cheat_weight;
0318 %            => new weight to set elements to
0319 
0320 
0321 pp= fwd_model_parameters( inv_model.fwd_model );
0322 idx= 1:pp.n_elem;
0323 weight= ones(1,pp.n_elem);
0324 weight( inv_model.cheat_tikhonov.cheat_elements ) = ...
0325         inv_model.cheat_tikhonov.cheat_weight;
0326 
0327 Reg = sparse( idx, idx, weight );
0328 
0329 % find elems which are connected to elems ee
0330 function elems= find_adjoin(ee, ELEM)
0331    nn= ELEM(:,ee);
0332    [d,e]= size(ELEM);
0333    ss= zeros(1,e);
0334    for i=1:d
0335      ss= ss+ any(ELEM==nn(i));
0336    end
0337    elems= find(ss==d-1);
0338 
0339 % Image prior with FEM based Laplacian + cheating
0340 function Reg= cheat_laplace( inv_model )
0341 % Reg= cheat_laplace( inv_model )
0342 % Reg        => output regularization term
0343 % Parameters:
0344 %   elems    = inv_model.RtR_prior.cheat_elements;
0345 %            => elements weights to modify
0346 %   weight   = inv_model.RtR_prior.cheat_weight;
0347 %            => new weight to set elements to
0348 
0349 % make pseudo laplacian filter
0350 %   roielems = region to exclude boundary
0351 %   if both ii and jj are in or out of roielems
0352 %   then include the term
0353    pp= fwd_model_parameters( inv_model.fwd_model );
0354 
0355    ROI = zeros(1,pp.n_elem);
0356    ROI( inv_model.cheat_laplace.cheat_elements ) = 1;
0357 
0358    Iidx= [];
0359    Jidx= [];
0360    Vidx= [];
0361    for ii=1:pp.n_elem
0362      el_adj = find_adjoin( ii, pp.ELEM );
0363      for jj=el_adj(:)'
0364          if (ROI(ii) + ROI(jj)) == 1 %one only
0365             fac= inv_model.cheat_laplace.cheat_weight *.5;
0366          else 
0367             fac = .5;
0368          end
0369          Iidx= [Iidx,      ii, ii, jj, jj];
0370          Jidx= [Jidx,      ii, jj, ii, jj];
0371          Vidx= [Vidx, fac*([1, -1, -1,  1]) ];
0372      end
0373    end
0374    Reg = sparse(Iidx,Jidx, Vidx, pp.n_elem, pp.n_elem );
0375 
0376 % deform a model by Delta
0377 % A1 = A0 + k1*sin(A0) + k2*cos(A0) + k3*sin(2*A0) ...
0378 function mdl1 = angl_deform(mdl0, def_amount );
0379 
0380    node0= mdl0.nodes';
0381    deform = def_amount*(8:-1:1).*randn(1,8);
0382 
0383    A0 = atan2( node0(2,:), node0(1,:) );
0384    R0 = sqrt( sum( node0.^2 ));
0385    A1= A0;
0386    for k= 2:2:length(deform)
0387       m = k/2;
0388       A1= A1 + deform(k-1)*cos(m*A0) + deform(k)*sin(m*A0);
0389    end
0390    node1 = [R0.*cos(A1); R0.*sin(A1) ];
0391 
0392    mdl1 = mdl0;
0393    mdl1.nodes = node1';
0394