0001 function img= inv_solve_gn( inv_model, data1, data2);
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0031
0032
0033
0034 if isfield(inv_model, 'inv_solve_gn')
0035 if isfield(inv_model, 'inv_solve_core')
0036 error('inv_model.inv_solve_gn replaces inv_model.inv_solve_core, parameters will be lost');
0037 end
0038 inv_model.inv_solve_core = inv_model.inv_solve_gn;
0039 inv_model = rmfield(inv_model, 'inv_solve_gn');
0040 end
0041
0042 if nargin > 2
0043 img = inv_solve_core(inv_model, data1, data2);
0044 else
0045 img = inv_solve_core(inv_model, data1);
0046 end
0047
0048 if isfield(img, 'inv_solve_core')
0049 img.inv_solve_gn = img.inv_solve_core;
0050 img=rmfield(img, 'inv_solve_core');
0051 end
0052
0053 function do_unit_test()
0054 imdl=mk_common_model('a2c');
0055 imdl.reconst_type='absolute';
0056 imdl.solve=@inv_solve_gn;
0057 fimg=mk_image(imdl,1);
0058 fimg.elem_data(5:10)=1.1;
0059 vi=fwd_solve(fimg);
0060 img=inv_solve(imdl,vi);
0061 clf; subplot(121); show_fem(fimg,1); title('forward model');
0062 subplot(122); show_fem(img,1); title('reconstruction');
0063 unit_test_cmp('fwd vs. reconst', fimg.elem_data, img.elem_data, 0.08);
0064
0065 do_unit_test_abs_diff();
0066
0067 function do_unit_test_abs_diff()
0068 imdl = mk_geophysics_model('h2c',32);
0069 imdl.solve = 'inv_solve_gn';
0070 imdl.RtR_prior = 'prior_tikhonov';
0071 imdl.inv_solve_gn.verbose = 0;
0072 imdl.inv_solve_gn.elem_working = 'log_conductivity';
0073
0074
0075 ctrs = interp_mesh(imdl.fwd_model);
0076 x = ctrs(:,1);
0077 y = ctrs(:,2);
0078 r1 = sqrt((x-20).^2 + (y+25).^2);
0079 r2 = sqrt((x-125).^2 + (y+45).^2);
0080 r3 = sqrt((x-50).^2 + (y+35).^2);
0081 imga = mk_image(imdl,1);
0082 imga.elem_data(r1<15)= 0.5;
0083 imga.elem_data(r2<20)= 2;
0084 va = fwd_solve(imga);
0085 imgb = imga;
0086 imgb.elem_data(r3<15)= 1.3;
0087 vb = fwd_solve(imgb);
0088 clf; subplot(121); show_fem(imga,1); subplot(122); show_fem(imgb,1); drawnow;
0089
0090 disp('*** INVERSE CRIME ALERT **** INVERSE CRIME ALERT ***');
0091 disp(' 1. This reconstruction has no measurement noise.');
0092 disp(' 2. This reconstruction is being performed on the same mesh the');
0093 disp(' measurements were simulated from: no discretization errors have');
0094 disp(' been introduced.');
0095 imdl.reconst_type = 'absolute';
0096 imgaa = inv_solve(imdl,va);
0097
0098 imdl.reconst_type = 'difference';
0099 imdl.inv_solve_gn.prior_data = imgaa.elem_data;
0100 imgab = inv_solve(imdl,va,vb);
0101 clf; subplot(221); show_fem(imga,1); title('A'); subplot(222); show_fem(imgb,1); title('B');
0102 subplot(223); show_fem(imgaa,1); title('rec A'); subplot(224); show_fem(imgab,1); title('rec B-A');
0103 unit_test_cmp('abs ', imga.elem_data, imgaa.elem_data, max(imga.elem_data)/2);
0104 imgba = imga; imgba.elem_data = imgb.elem_data - imga.elem_data;
0105 imgba.elem_data = imgba.elem_data/max(imgba.elem_data);
0106 imgab.elem_data = imgab.elem_data/max(imgab.elem_data);
0107 unit_test_cmp('diff', norm(imgba.elem_data - imgab.elem_data), 0, 20);