0001 function img = inv_solve( 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
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0055
0056 inv_model = prepare_model( inv_model );
0057 opts = parse_parameters( inv_model );
0058
0059 print_hello(inv_model.solve);
0060 try inv_model.solve = str2func(inv_model.solve); end
0061
0062 if opts.abs_solve
0063 if nargin~=2;
0064 error('only one data set is allowed for a static reconstruction');
0065 end
0066
0067 fdata = filt_data(inv_model,data1);
0068
0069 imgc= eidors_cache( inv_model.solve, {inv_model, fdata}, 'inv_solve');
0070
0071 else
0072 if nargin~=3;
0073 error('two data sets are required for a difference reconstruction');
0074 end
0075
0076
0077 data_width= max(num_frames(data1), num_frames(data2));
0078
0079 fdata1 = filt_data( inv_model, data1, data_width );
0080 fdata2 = filt_data( inv_model, data2, data_width );
0081
0082 imgc= eidors_cache( inv_model.solve, {inv_model, fdata1, fdata2},'inv_solve');
0083
0084
0085 end
0086
0087 check_parametrization_handling(inv_model,imgc);
0088
0089 img = eidors_obj('image', imgc );
0090
0091 if ~isfield(img,'current_params')
0092 img.current_params = [];
0093 end
0094
0095
0096
0097 if isfield(inv_model,'rec_model')
0098 img.fwd_model= inv_model.rec_model;
0099 end
0100
0101
0102 if opts.select_parameters;
0103 img.elem_data = img.elem_data( opts.select_parameters, :);
0104 end
0105
0106 if ~opts.reconst_to_elems
0107 img.node_data= img.elem_data;
0108 img = rmfield(img,'elem_data');
0109 end
0110
0111
0112 try; img.elem_data = opts.offset + opts.scale * img.elem_data; end
0113 try; img.node_data = opts.offset + opts.scale * img.node_data; end
0114
0115
0116
0117
0118
0119 img.info.error = NaN;
0120 img = orderfields(img);
0121
0122
0123 try
0124 do_calc = inv_model.inv_solve.calc_solution_error;
0125 catch
0126 eidors_msg('inv_solve: not calculting solution residual',3);
0127 do_calc = false;
0128 end
0129 if ~do_calc;
0130 return;
0131 end
0132 eidors_msg('inv_solve: Calculating solution residual (inv_model.inv_solve.calc_solution_error = 0 to disable)',2);
0133 try
0134 if opts.abs_solve
0135 img.info.error = calc_solution_error( imgc, inv_model, fdata);
0136 else
0137 img.info.error = calc_solution_error( imgc, inv_model, fdata1, fdata2 );
0138 end
0139 eidors_msg('inv_solve: Solution Error: %f', img.info.error, 2);
0140 catch
0141 eidors_msg('inv_solve: Solution Error calculation failed.',2);
0142 end
0143
0144 function print_hello(solver)
0145 if isa(solver,'function handle')
0146 solver = func2str(solver);
0147 end
0148 if strcmp(solver, 'eidors_default');
0149 solver = eidors_default('get','inv_solve');
0150 end
0151 eidors_msg('inv_solve: %s', solver,3);
0152
0153
0154 function opts = parse_parameters( imdl )
0155 if strcmp(imdl.reconst_type,'static') || ...
0156 strcmp(imdl.reconst_type,'absolute')
0157 opts.abs_solve = 1;
0158 elseif strcmp(imdl.reconst_type,'difference')
0159 opts.abs_solve = 0;
0160 else
0161 error('inv_model.reconst_type (%s) not understood', imdl.reconst_type);
0162 end
0163
0164 opts.select_parameters = [];
0165 try
0166 opts.select_parameters = imdl.inv_solve.select_parameters;
0167 end
0168
0169 opts.reconst_to_elems = 1;
0170 try if strcmp( imdl.reconst_to, 'nodes' )
0171 opts.reconst_to_elems = 0;
0172 end; end
0173
0174 opts.scale = 1;
0175 try opts.scale = imdl.inv_solve.scale_solution.scale; end
0176
0177 opts.offset = 0;
0178 try opts.offset = imdl.inv_solve.scale_solution.offset; end
0179
0180 function imdl = deprecate_imdl_parameters(imdl)
0181 if isfield(imdl, 'parameters')
0182 if isstruct(imdl.parameters)
0183
0184
0185
0186
0187 if ~isfield(imdl, 'inv_solve')
0188 imdl.inv_solve = imdl.parameters;
0189 else
0190
0191
0192 A = imdl.parameters;
0193 B = imdl.inv_solve;
0194 M = [fieldnames(A)' fieldnames(B)'; struct2cell(A)' struct2cell(B)'];
0195 try
0196 imdl.inv_solve=struct(M{:});
0197 catch
0198 [tmp, rows] = unique(M(1,:), 'last');
0199 M=M(:,rows);
0200 imdl.inv_solve=struct(M{:});
0201 end
0202 end
0203 imdl = rmfield(imdl, 'parameters');
0204 imdl.parameters = imdl.inv_solve;
0205 else
0206 error('unexpected inv_model.parameters where parameters is not a struct... i do not know what to do');
0207 end
0208 end
0209
0210 function mdl = prepare_model( mdl )
0211 fmdl = mdl.fwd_model;
0212 fmdl = mdl_normalize(fmdl,mdl_normalize(fmdl));
0213 if ~isfield(fmdl,'elems');
0214 return;
0215 end
0216
0217 fmdl.elems = double(fmdl.elems);
0218 fmdl.n_elem = size(fmdl.elems,1);
0219 fmdl.n_node = size(fmdl.nodes,1);
0220 if isfield(fmdl,'electrode');
0221 fmdl.n_elec = length(fmdl.electrode);
0222 else
0223 fmdl.n_elec = 0;
0224 end
0225
0226 mdl.fwd_model= fmdl;
0227 if ~isfield(mdl,'reconst_type');
0228 mdl.reconst_type= 'difference';
0229 end
0230
0231
0232 mdl = deprecate_imdl_parameters(mdl);
0233
0234 function check_parametrization_handling(inv_model,imgc)
0235 if isfield(inv_model, 'jacobian_bkgnd') && ...
0236 has_params(inv_model.jacobian_bkgnd) && ~has_params(imgc)
0237 if isa(inv_model.solve,'function_handle')
0238 solver = func2str(inv_model.solve);
0239 else
0240 solver = inv_model.solve;
0241 end
0242 if strcmp(solver,'eidors_default');
0243 solver = eidors_default('get','inv_solve');
0244 end
0245 warning('EIDORS:ParamsObliviousSolver',...
0246 ['The solver %s did not handle the parametrization data properly.\n'...
0247 'The results may be incorrect. Please check the code to verify.'], ...
0248 solver);
0249 end
0250
0251
0252 function b = has_params(s)
0253 b = false;
0254 if isstruct(s)
0255 b = any(ismember(fieldnames(s),supported_params));
0256 end
0257
0258
0259
0260 function nf= num_frames(d0)
0261 if isnumeric( d0 )
0262 nf= size(d0,2);
0263 elseif d0(1).type == 'data';
0264 nf= size( horzcat( d0(:).meas ), 2);
0265 else
0266 error('Problem calculating number of frames. Expecting numeric or data object');
0267 end
0268
0269
0270 function d2= filt_data(inv_model, d0, data_width )
0271 if ~isnumeric( d0 )
0272
0273
0274 d1 = [];
0275 for i=1:length(d0)
0276 if strcmp( d0(i).type, 'data' )
0277 d1 = [d1, d0(i).meas];
0278 else
0279 error('expecting an object of type data');
0280 end
0281 end
0282
0283 else
0284
0285 d1 = d0;
0286 end
0287
0288 d1= double(d1);
0289
0290 if isfield(inv_model.fwd_model,'meas_select') && ...
0291 ~isempty(inv_model.fwd_model.meas_select)
0292
0293
0294 meas_select= inv_model.fwd_model.meas_select;
0295 if size(d1,1) == length(meas_select)
0296 d2= d1(meas_select,:);
0297 elseif size(d1,1) == sum(meas_select==1)
0298 d2= d1;
0299 else
0300 error('inconsistent difference data: (%d ~= %d). Maybe check fwd_model.meas_select', ...
0301 size(d1,1), length(meas_select));
0302 end
0303 else
0304 d2= d1;
0305 end
0306
0307 if nargin==3
0308 d2_width= size(d2,2);
0309 if d2_width == data_width
0310
0311 elseif d2_width == 1
0312 d2= d2(:,ones(1,data_width));
0313 else
0314 error('inconsistent difference data: (%d ~= %d)', ...
0315 d2_width, data_width);
0316 end
0317 end
0318
0319
0320 function do_unit_test
0321 k=0; N=5; nd = 5;
0322
0323 imdl = mk_common_model('d2c2',16);
0324 imdl = select_imdl( imdl, {'Choose NF=2.5'});
0325 mvx = linspace(-0.8,0.8,nd);
0326 [vh,vi] = simulate_movement(mk_image(imdl), [mvx;0*mvx;0.05+0*mvx]);
0327
0328 img= inv_solve(imdl,vh,vi); img.show_slices.img_cols = 5;
0329 k=k+1; subplot(N,1,k); show_slices(img);
0330 unit_test_cmp('inv_solve: 1a', mean(img.elem_data), 3e-3, 1e-3);
0331 unit_test_cmp('inv_solve: 1b', std(img.elem_data), 1.5e-2, 1e-2);
0332
0333 vhm= eidors_obj('data','nn','meas',vh);
0334 img= inv_solve(imdl,vhm,vi);
0335 unit_test_cmp('inv_solve: 2a', mean(img.elem_data), 3e-3, 1e-3);
0336 unit_test_cmp('inv_solve: 2b', std(img.elem_data), 1.5e-2, 1e-2);
0337
0338 img= inv_solve(imdl,vh*ones(1,nd),vi);
0339 unit_test_cmp('inv_solve: 3a', mean(img.elem_data), 3e-3, 1e-3);
0340 unit_test_cmp('inv_solve: 3b', std(img.elem_data), 1.5e-2, 1e-2);
0341
0342 vim= eidors_obj('data','nn','meas',vi);
0343 img= inv_solve(imdl,vhm,vim);
0344 unit_test_cmp('inv_solve: 4a', mean(img.elem_data), 3e-3, 1e-3);
0345 unit_test_cmp('inv_solve: 4b', std(img.elem_data), 1.5e-2, 1e-2);
0346
0347 vhm= eidors_obj('data','nn','meas',vh*ones(1,nd));
0348 img= inv_solve(imdl,vhm,vi);
0349 unit_test_cmp('inv_solve: 5a', mean(img.elem_data), 3e-3, 1e-3);
0350 unit_test_cmp('inv_solve: 5b', std(img.elem_data), 1.5e-2, 1e-2);
0351
0352 vhm(1)= eidors_obj('data','nn','meas',vh*ones(1,2));
0353 vhm(2)= eidors_obj('data','nn','meas',vh*ones(1,nd-2));
0354 img= inv_solve(imdl,vhm,vi);
0355 unit_test_cmp('inv_solve: 6a', mean(img.elem_data), 3e-3, 1e-3);
0356 unit_test_cmp('inv_solve: 6b', std(img.elem_data), 1.5e-2, 1e-2);
0357
0358 vim(1)= eidors_obj('data','nn','meas',vi(:,1:3));
0359 vim(2)= eidors_obj('data','nn','meas',vi(:,4:end));
0360 img= inv_solve(imdl,vhm,vim);
0361 unit_test_cmp('inv_solve: 7a', mean(img.elem_data), 3e-3, 1e-3);
0362 unit_test_cmp('inv_solve: 7b', std(img.elem_data), 1.5e-2, 1e-2);
0363
0364 im2 = imdl; im2.inv_solve.select_parameters = 1:5;
0365 img= inv_solve(im2,vh,vi);
0366 unit_test_cmp('inv_solve: 10', size(img.elem_data), [5,nd]);
0367
0368
0369 im2 = imdl;
0370 im2.inv_solve.scale_solution.offset = 1;
0371 im2.inv_solve.scale_solution.scale = 2;
0372 img= inv_solve(im2,vh,vi); img.show_slices.img_cols = 5;
0373 unit_test_cmp('inv_solve: 20a', mean(img.elem_data), 1.006, 2e-3);
0374 unit_test_cmp('inv_solve: 20b', std(img.elem_data), 3e-2, 1e-2);
0375 im2.inv_solve.calc_solution_error = 1;
0376 img= inv_solve(im2,vh,vi);
0377 unit_test_cmp('inv_solve: 20e', mean(img.elem_data), 1.006, 2e-3);
0378
0379 im2.inv_solve.scale_solution.offset = 0;
0380 d = interp_mesh( imdl.fwd_model); d= sqrt(sum(d.^2,2));
0381 im2.inv_solve.scale_solution.scale = spdiags(1-d,0,length(d),length(d));
0382 img= inv_solve(imdl,vh,vi); img.show_slices.img_cols = 5;
0383 k=k+1; subplot(N,1,k); show_slices(img);
0384
0385 im2 = select_imdl(imdl, {'Nodal GN dif'} );
0386 img= inv_solve(im2,vh,vi);
0387 unit_test_cmp('inv_solve: 30a', mean(img.node_data), 3e-3, 1e-3);
0388 unit_test_cmp('inv_solve: 30b', std(img.node_data), 1.5e-2, 1e-2);
0389
0390 im2.inv_solve.scale_solution.offset = 1;
0391 im2.inv_solve.scale_solution.scale = 2;
0392 img= inv_solve(im2,vh,vi); img.show_slices.img_cols = 5;
0393 unit_test_cmp('inv_solve: 31a', mean(img.node_data), 1.006, 2e-3);
0394 unit_test_cmp('inv_solve: 31b', std(img.node_data), 3e-2, 1.5e-2);
0395
0396 im2 = select_imdl(imdl, {'Basic GN abs'} );
0397 im2.inv_solve.calc_solution_error = 0;
0398 img= inv_solve(im2,vi(:,1));
0399 unit_test_cmp('inv_solve: 40a', mean(img.elem_data), 1.004, 1e-3);
0400 unit_test_cmp('inv_solve: 40b', std(img.elem_data), 1.5e-2, 1e-2);
0401 im2.inv_solve.calc_solution_error = 1;
0402 img= inv_solve(im2,vi(:,1));
0403 unit_test_cmp('inv_solve: 40e', mean(img.elem_data), 1.004, 1e-3);
0404
0405