0001 function [NF,SE] = calc_noise_figure( inv_model, hp, iterations)
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 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'), do_unit_test, return, end
0048
0049 if nargin>=2 && numel(hp) == 1
0050 inv_model.hyperparameter.value= hp;
0051
0052 try; inv_model.hyperparameter = rmfield(inv_model.hyperparameter,'func'); end
0053 end
0054 if nargin == 3 && numel(hp) > 1
0055 h_data = hp;
0056 c_data = iterations;
0057 else
0058 [inv_model, h_data, c_data] = process_parameters( inv_model );
0059 end
0060
0061
0062 inv_model.inv_solve.calc_solution_error = false;
0063
0064
0065
0066 if nargin<3; iterations= 10; end
0067 solver = inv_model.solve;
0068 if ischar(solver) && strcmp(solver, 'eidors_default')
0069 solver = eidors_default('get','inv_solve');
0070 end
0071 if isa(solver,'function_handle')
0072 solver = func2str(solver);
0073 end
0074 switch solver
0075 case {'inv_solve_backproj'
0076 'inv_solve_conj_grad'
0077 'inv_solve_diff_GN_one_step'
0078 'inv_solve_trunc_iterative'
0079 'inv_solve_TSVD'
0080 'solve_use_matrix'}
0081 NF = nf_calc_linear( inv_model, h_data, c_data);
0082 SE = 0;
0083 otherwise
0084 [NF,SE]= nf_calc_random( inv_model, h_data, c_data, iterations);
0085 end
0086 eidors_msg('@@ NF= %f', NF, 2);
0087
0088 function [inv_model, h_data, c_data] = process_parameters( inv_model );
0089
0090 if isfield(inv_model.hyperparameter,'tgt_elems')
0091 [h_data, c_data]= simulate_targets( inv_model.fwd_model, ...
0092 inv_model.hyperparameter.tgt_elems);
0093 elseif isfield(inv_model.hyperparameter,'tgt_data')
0094 tgt_data= inv_model.hyperparameter.tgt_data;
0095 h_data= tgt_data.meas_t1;
0096 c_data= tgt_data.meas_t2;
0097 else
0098 error('unsure how to get data to measure signal');
0099 end
0100
0101
0102
0103 if isfield(inv_model.hyperparameter,'func')
0104 funcname= inv_model.hyperparameter.func;
0105 if strcmp( class(funcname), 'function_handle')
0106 funcname= func2str(funcname);
0107 end
0108
0109 if strcmp(funcname, 'choose_noise_figure')
0110 error('specifying inv_model.hp.func = choose_noise_figure will recurse');
0111 end
0112 end
0113
0114 function params = nf_calc_linear(imdl, vh, vi )
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125 if 0
0126
0127 Nnoise = 1000;
0128 noise = 0.01*std(vh)*randn(size(vh,1),Nnoise);
0129 vhn= vh*ones(1,Nnoise) + noise;
0130 else
0131 noise = 0.01*std(vh)*speye(size(vh,1));
0132 vhn= vh*ones(1,size(vh,1)) + noise;
0133 end
0134
0135 signal_y = calc_difference_data( vh, vi, imdl.fwd_model);
0136 noise_y = calc_difference_data( vh, vhn, imdl.fwd_model);
0137
0138 signal_x = inv_solve(imdl, vh, vi);
0139 signal_x = data_mapper(signal_x); signal_x = signal_x.elem_data;
0140 noise_x = inv_solve(imdl, vh, vhn);
0141 noise_x = data_mapper(noise_x); noise_x = noise_x.elem_data;
0142
0143 use_rec = 1;
0144 try
0145 use_rec = ~imdl.prior_use_fwd_not_rec;
0146 end
0147 if use_rec
0148 try
0149 VOL = get_elem_volume(imdl.rec_model);
0150 catch
0151 VOL = get_elem_volume(imdl.fwd_model);
0152 end
0153 else
0154 VOL = get_elem_volume(imdl.fwd_model);
0155 end
0156 VOL = spdiags(VOL,0, length(VOL), length(VOL));
0157
0158 signal_x = VOL*signal_x;
0159 noise_x = VOL*noise_x;
0160
0161 signal_x = mean(abs(signal_x),1);
0162 noise_x = mean(std(noise_x));
0163 snr_x = signal_x / noise_x;
0164
0165 signal_y = mean(abs(signal_y),1);
0166 noise_y = mean(std(noise_y));
0167 snr_y = signal_y / noise_y;
0168
0169 params= [snr_y(:)./snr_x(:)]';
0170
0171 try
0172 eidors_msg('@@ NF= %f (hp=%e)', NF, imdl.hyperparameter.value, 2);
0173 end
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193 function NF= nf_calc_use_matrix( inv_model, h_data, c_data)
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211 try
0212 VOL = get_elem_volume(inv_model.rec_model)';
0213 catch
0214 VOL = get_elem_volume(inv_model.fwd_model)';
0215 end
0216
0217
0218 d_len = size(h_data,1);
0219 delta = 1e-2* mean(h_data);
0220 c_noise = c_data*ones(1,d_len) + eye(d_len);
0221 h_full = h_data*ones(1,d_len);
0222
0223 sig_data = mean(abs( ...
0224 calc_difference_data( h_data, c_data , inv_model.fwd_model ) ...
0225 ));
0226 var_data = mean(sum( ...
0227 calc_difference_data( h_full, c_noise, inv_model.fwd_model ) ...
0228 .^2, 2));
0229
0230
0231
0232
0233 [img0, img0n] = get_images( inv_model, h_data, c_data, ...
0234 h_full, c_noise);
0235
0236 i_len = length(img0);
0237 sig_img= VOL*abs(img0) / i_len;;
0238 var_img= VOL.^2*sum(img0n.^2 ,2) / i_len;
0239
0240 NF = ( sig_data/ sqrt(var_data) ) / ( sig_img / sqrt(var_img) );
0241
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251 function [h_data, c_data]= simulate_targets( fwd_model, ctr_elems)
0252
0253 homg= 1;
0254
0255
0256 sigma= homg*ones( size(fwd_model.elems,1) ,1);
0257
0258 img= eidors_obj('image', 'homogeneous image', ...
0259 'elem_data', sigma, ...
0260 'fwd_model', fwd_model );
0261 h_data=fwd_solve( img );
0262 h_data= h_data.meas;
0263
0264
0265 delta = 1e-2;
0266 sigma(ctr_elems) = homg*(1 + delta);
0267 img.elem_data = sigma;
0268 c_data=fwd_solve( img );
0269 c_data= c_data.meas;
0270
0271 function [img0, img0n] = get_images( inv_model, h_data, c_data, ...
0272 h_full, c_noise);
0273 if isa(inv_model.solve,'function_handle')
0274 solve= func2str(inv_model.solve);
0275 else
0276 solve= inv_model.solve;
0277 end
0278
0279
0280 switch solve
0281 case 'ab_tv_diff_solve'
0282 error('Dont know how to calculate TV noise figure')
0283
0284 case 'inv_solve_diff_kalman'
0285 inv_model.inv_solve_diff_kalman.keep_K_k1= 1;
0286 stablize = 6;
0287 img0 = inv_solve( inv_model, h_data, ...
0288 c_data*ones(1,stablize) );
0289 K= img0.inv_solve_diff_kalman.K_k1;
0290 img0.elem_data = K*calc_difference_data( h_data , c_data , inv_model.fwd_model);
0291 img0n.elem_data= K*calc_difference_data( h_full , c_noise, inv_model.fwd_model);
0292
0293 otherwise
0294 img0 = inv_solve( inv_model, h_data, c_data);
0295 if nargin>4
0296 img0n= inv_solve( inv_model, h_full, c_noise);
0297 end
0298 end
0299
0300
0301 if isfield(img0,'node_data');
0302 img0 = img0.node_data;
0303 else
0304 img0 = img0.elem_data;
0305 end
0306
0307 if isfield(img0n,'node_data');
0308 img0n = img0n.node_data;
0309 else
0310 img0n = img0n.elem_data;
0311 end
0312
0313
0314 function NF= nf_calc_iterate( inv_model, h_data, c_data);
0315 try
0316 VOL = get_elem_volume(inv_model.rec_model)';
0317 catch
0318 VOL = get_elem_volume(inv_model.fwd_model)';
0319 end
0320
0321 d_len = size(h_data,1);
0322 delta = 1e-2* mean(h_data);
0323 sig_data = mean(abs( ...
0324 calc_difference_data( h_data, c_data , inv_model.fwd_model ) ...
0325 ));
0326
0327
0328
0329 [img0] = get_images( inv_model, h_data, c_data);
0330
0331 sig_img= VOL*abs(img0) / length(img0);
0332
0333
0334 var_data= 0;
0335 var_img = 0;
0336 for i=1:d_len
0337 this_noise = -ones(d_len, size(c_data,2))/(d_len-1);
0338 this_noise(i,:) = 1;
0339 c_noise = c_data + this_noise;
0340 [imgn] = get_images( inv_model, h_data, c_noise);
0341 if 1
0342 var_data = var_data + mean(sum( ...
0343 calc_difference_data( h_data, c_noise, inv_model.fwd_model ) ...
0344 .^2, 2));
0345
0346 var_img= var_img + (VOL.^2)*sum(imgn.elem_data.^2,2 ) / length(imgn.elem_data);
0347 else
0348
0349 var_data = var_data + var( ...
0350 calc_difference_data( h_data, c_noise, inv_model.fwd_model ) ...
0351 );
0352 var_img= var_img + var( VOL'.*imgn.elem_data );
0353 end
0354 end
0355 var_data = var_data / d_len;
0356 var_img = var_img / d_len;
0357 NF = ( sig_data/ sqrt(var_data) ) / ( sig_img / sqrt(var_img) );
0358
0359 function [NF,SE]= nf_calc_random( rec, vh, vi, N_RUNS);
0360 eidors_cache('boost_priority',-2);
0361
0362 imgr= inv_solve(rec, vh, vi);
0363
0364 if isfield(imgr,'node_data');
0365 img0 = imgr.node_data;
0366 try
0367 VOL = get_elem_volume(rec.rec_model, 1);
0368 catch
0369 VOL = get_elem_volume(rec.fwd_model, 1);
0370 end
0371 else
0372 n_els = length(rec.fwd_model);
0373 img0 = imgr.elem_data(1:n_els);
0374 try
0375 VOL = get_elem_volume(rec.rec_model, 0);
0376 catch
0377 VOL = get_elem_volume(rec.fwd_model, 0);
0378 end
0379 end
0380
0381 sig_ampl = mean( abs( VOL .* img0 )) / ...
0382 mean( abs( calc_difference_data( vh, vi, rec.fwd_model )));
0383
0384
0385 for i=1:N_RUNS
0386 vn= addnoise(vh, vi, 1.0);
0387
0388 imgr= inv_solve(rec, vh, vn);
0389
0390 if isfield(imgr,'node_data'); img0 = imgr.node_data;
0391 else; img0 = imgr.elem_data(1:n_els);
0392 end
0393
0394 noi_imag(i) = std( VOL .* img0 );
0395 noi_sgnl(i) = std( calc_difference_data( vh, vn, rec.fwd_model ));
0396 end
0397
0398 noi_ampl = noi_imag./noi_sgnl;
0399 NF = mean(noi_ampl/sig_ampl);
0400 SE = std(noi_ampl/sig_ampl)/sqrt(N_RUNS);
0401 eidors_msg('NF= %f+/-%f', NF, SE, 3);
0402
0403 eidors_cache('boost_priority',2);
0404
0405 function noise= addnoise( vh, vi, SNR);
0406 if isstruct(vh); vh= vh.meas; end
0407 if isstruct(vi); vi= vi.meas; end
0408 noise = randn(size(vh));
0409 noise = noise*std(vh-vi)/std(noise);
0410 noise = vh + SNR*noise;
0411
0412 function do_unit_test
0413 ll = eidors_msg('log_level',1);
0414 test1;
0415 imdl = mk_common_model('a2t2',16); test2(imdl);
0416 imdl = mk_common_model('d2t2',16); test2(imdl);
0417 imdl = mk_common_model('a2c2',16); test2(imdl);
0418 imdl = mk_common_model('d2c2',16); test2(imdl);
0419 ll = eidors_msg('log_level',ll);
0420
0421 function test1
0422
0423 fmdl = mk_library_model('pig_23kg_16el');
0424 [fmdl.stimulation fmdl.meas_select] = mk_stim_patterns(16,1,'{ad}','{ad}');
0425 fmdl = mdl_normalize(fmdl, 1);
0426 opt.noise_figure = 0.5; opt.imgsz = [64 64];
0427 imdl = mk_GREIT_model(fmdl, 0.25, [], opt);
0428
0429 img = mk_image(fmdl,1);
0430 vh = fwd_solve(img);
0431
0432 select_fcn = inline('(x-0).^2+(y-0).^2+(z-0.5).^2<0.1^2','x','y','z');
0433 mfrac = elem_select(fmdl, select_fcn);
0434 img.elem_data = img.elem_data + mfrac*0.1;
0435 vi = fwd_solve(img);
0436
0437 nf1 = calc_noise_params(imdl, vh.meas, vi.meas);
0438
0439
0440 imdl.hyperparameter.tgt_data.meas_t1 = vh.meas;
0441 imdl.hyperparameter.tgt_data.meas_t2 = vi.meas;
0442 try
0443
0444 nf2 = calc_noise_figure(imdl);
0445 catch
0446 nf2 = 0;
0447 end
0448 unit_test_cmp('Noise fig implementations',nf1, nf2, 1e-2);
0449
0450 function test2(imdl)
0451 fmdl = imdl.fwd_model;
0452
0453 img = mk_image(fmdl,1);
0454 vh = fwd_solve(img);
0455
0456 select_fcn = inline('(x-0).^2+(y-0).^2.^2<15^2','x','y','z');
0457 mfrac = elem_select(fmdl, select_fcn);
0458 img.elem_data = img.elem_data + mfrac*0.1;
0459 vi = fwd_solve(img);
0460
0461 nf1 = calc_noise_params(imdl, vh.meas, vi.meas);
0462 eidors_msg(nf1,0);
0463
0464 imdl.hyperparameter.tgt_data.meas_t1 = vh.meas;
0465 imdl.hyperparameter.tgt_data.meas_t2 = vi.meas;
0466
0467 nf2 = calc_noise_figure(imdl,[],1000);
0468 unit_test_cmp('Noise fig implementations',nf1, nf2, 1e-2);
0469