0001 function [imdl, weight]= mk_GREIT_model( fmdl, radius, weight, options )
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
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087 citeme(mfilename);
0088
0089 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0090
0091 if nargin < 4, options = [];end
0092 [imdl,fmdl,imgs] = parse_fmdl(fmdl);
0093 options = parse_options(options,fmdl,imdl, weight);
0094
0095 copt.cache_obj= { fmdl, imdl, imgs, radius, weight, options};
0096 copt.fstr = 'mk_GREIT_model';
0097 params = {fmdl, imdl, imgs, radius, weight, options};
0098
0099 [imdl, weight] = eidors_cache(@mk_GREIT_model_calc, params, copt);
0100
0101
0102 function [imdl, weight]= mk_GREIT_model_calc( fmdl, imdl, imgs, radius, weight, opt)
0103
0104 Nsim = opt.Nsim;
0105 [vi,vh,xyz,opt]= stim_targets(imgs, Nsim, opt );
0106
0107
0108 if ~isfield(imdl,'rec_model');
0109 opt.do_coarse2fine = 0;
0110 [imdl.rec_model imdl.fwd_model] = mk_pixel_slice(fmdl,opt.target_plane,opt);
0111 imdl.rec_model.nodes(:,3) = [];
0112
0113 imdl.rec_model.mdl_slice_mapper.y_pts = fliplr(imdl.rec_model.mdl_slice_mapper.y_pts);
0114 end
0115
0116 opt.rec_model = imdl.rec_model;
0117
0118 imdl.solve = @solve_use_matrix;
0119
0120
0121 if ~isempty(opt.noise_figure) || ~isempty(opt.image_SNR)
0122 if ~isempty(opt.noise_figure)
0123
0124 target = opt.noise_figure;
0125 NoisPerfName = 'Noise Figure';
0126 elseif ~isempty(opt.image_SNR)
0127
0128 NoisPerfName = 'Image SNR';
0129 target = opt.image_SNR;
0130 if isfield(opt, 'SigmaN')
0131 imdl.hyperparameter.SigmaN = opt.SigmaN;
0132 end
0133 if isfield(opt, 'image_SNR_targets')
0134 imdl.hyperparameter.xyzr_targets = opt.image_SNR_targets;
0135 end
0136 else
0137 error('internal bug: shouldn''t get here');
0138 end
0139
0140 if ~isempty(weight)
0141 eidors_msg('mk_GREIT_model: Using weight parameter as a guess, options.noise_figure or opt.image_SNR is non-empty');
0142 else
0143 if ~isempty(opt.noise_figure)
0144 weight = target;
0145 elseif ~isempty(opt.image_SNR)
0146 weight = 1/target;
0147 else
0148 error('internal bug: shouldn''t get here');
0149 end
0150 end
0151
0152 xyzr = opt.noise_figure_targets;
0153 [jnk,vi_NF] = simulate_movement(imgs,xyzr');
0154 vi_NF = sum(vi_NF,2);
0155 eidors_msg(['mk_GREIT_model: Finding noise weighting for given ', NoisPerfName],1);
0156 eidors_msg('mk_GREIT_model: This will take a while...',1);
0157 f = @(X) to_optimise(vh,vi,xyz, radius, X, opt, imdl, target, vi_NF);
0158 fms_opts.TolFun = 0.01*target;
0159
0160
0161 imdl.solve_use_matrix.RM = calc_GREIT_RM(vh,vi,xyz, radius, weight, opt);
0162 log_level = eidors_msg('log_level');
0163 if log_level > 1
0164 log_level = eidors_msg( 'log_level', 1);
0165 end
0166 [weight, NF] = fminsearch(f, weight,fms_opts);
0167 eidors_msg(['mk_GREIT_model: Optimal solution gives ', NoisPerfName, '=' ...
0168 num2str(NF+target) ' with weight=' num2str(weight)],1);
0169 assert((sqrt(NF) / target) < 0.01, ...
0170 ['Cannot find an accurate enough match for desired ', NoisPerfName]');
0171 eidors_msg( 'log_level', log_level);
0172 end
0173
0174 [RM, PJt, M] = calc_GREIT_RM(vh,vi, xyz, radius, weight, opt );
0175 imdl.solve_use_matrix.RM = RM;
0176 if opt.keep_intermediate_results
0177
0178 imdl.solve_use_matrix.PJt = PJt;
0179 imdl.solve_use_matrix.X = inv(M);
0180 end
0181
0182 imdl.jacobian_bkgnd = imgs;
0183
0184 imdl.hyperparameter.value = weight;
0185
0186 function out = to_optimise(vh,vi,xy,radius,weight, opt, imdl, ...
0187 target,vi_NF)
0188
0189
0190 imdl.solve_use_matrix.RM = calc_GREIT_RM(vh,vi,xy, radius, weight, opt);
0191 imdl.hyperparameter.value = weight;
0192
0193 if ~isempty(opt.noise_figure)
0194 NF = calc_noise_figure(imdl,vh, vi_NF);
0195 eidors_msg(['NF = ', num2str(NF), ' weight = ', num2str(weight)],1);
0196 elseif ~isempty(opt.image_SNR)
0197 NF = calc_image_SNR(imdl);
0198 eidors_msg(['SNR = ', num2str(NF), ' weight = ', num2str(weight)],1);
0199 else
0200 error('internal bug: shouldn''t get here');
0201 end
0202 out = (NF - target)^2;
0203
0204
0205
0206 function imgs = get_prepackaged_fmdls( fmdl );
0207 switch fmdl
0208 case 'c=1;h=2;r=.08;ce=16;bg=1;st=1;me=1;nd'
0209 fmdl = ng_mk_cyl_models([2,1,0.18],[16,1],[0.05]);
0210 fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
0211 fmdl = mdl_normalize(fmdl,1);
0212 imgs= mk_image( fmdl, 1);
0213 otherwise
0214 error('specified fmdl (%s) is not understood', fmdl);
0215 end
0216
0217 function [vi,vh,xyz,opt]= stim_targets(imgs, Nsim, opt );
0218 fmdl = imgs.fwd_model;
0219 ctr = mean(fmdl.nodes);
0220 maxx = max(abs(fmdl.nodes(:,1) - ctr(1)));
0221 maxy = max(abs(fmdl.nodes(:,2) - ctr(2)));
0222 if numel(opt.distr) > 1
0223 xyzr = opt.distr;
0224 xyzr(4,:) = calc_radius(mean([maxx,maxy]),opt,size(opt.distr,2));
0225 else
0226 switch opt.distr
0227 case 0
0228 r = linspace(0,0.9, Nsim);
0229 th = r*4321;
0230 xyzr = [maxx*r.*cos(th); maxy*r.*sin(th);
0231 opt.target_plane*ones(1,Nsim);
0232 0.05*mean([maxx,maxy])*ones(1,Nsim)];
0233
0234 case 1
0235 F = fourier_fit(opt.contour_boundary(:,1:2));
0236 v = linspace(0,1,Nsim*100+1); v(end)=[];
0237 pts = fourier_fit(F,v);
0238 idx_p = floor(rand(Nsim,1)*Nsim*100);
0239 xyzr = pts(idx_p,:)'.*repmat(rand(Nsim,1),[1 2])';
0240 xyzr(3,:) = calc_offset(opt.target_plane,opt,Nsim);
0241
0242
0243 xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,Nsim);
0244 case 2
0245
0246
0247
0248 pts = opt.contour_boundary(:,1:2);
0249
0250 pts = 0.9*( pts - repmat(ctr(1:2),length(pts),1) ) + repmat(ctr(1:2),length(pts),1);
0251
0252
0253 lim = max(maxx, maxy);
0254 x = ctr(1) + (rand(Nsim*10,1)-0.5)*2*lim;
0255 y = ctr(2) + (rand(Nsim*10,1)-0.5)*2*lim;
0256 IN = inpolygon(x,y,pts(:,1),pts(:,2));
0257 xyzr(1,:) = x(find(IN,Nsim));
0258 xyzr(2,:) = y(find(IN,Nsim));
0259 xyzr(3,:) = calc_offset(opt.target_plane,opt,Nsim);
0260
0261 xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,Nsim);
0262 case 3
0263
0264
0265
0266 pts = opt.contour_boundary(:,1:2);
0267 lim = max(maxx, maxy);
0268 frac = polyarea(pts(:,1),pts(:,2)) / (2*lim)^2;
0269 [x,y] = ndgrid( linspace(-lim,lim,ceil(sqrt(Nsim/frac))), ...
0270 linspace(-lim,lim,ceil(sqrt(Nsim/frac))));
0271
0272 x = x+ctr(1); y = y + ctr(2);
0273 IN = inpolygon(x,y,pts(:,1),pts(:,2));
0274 xyzr(1,:) = x(find(IN));
0275 xyzr(2,:) = y(find(IN));
0276 xyzr(3,:) = calc_offset(opt.target_plane,opt,size(xyzr,2));
0277
0278 xyzr(4,:) = calc_radius(mean([maxx maxy]),opt,size(xyzr,2));
0279 eidors_msg(['mk_GREIT_model: Using ' num2str(size(xyzr,2)) ' points']);
0280 otherwise; error('GREIT: opt.distr no such case=%d',opt.distr);
0281 end
0282 end
0283 before = size(xyzr,2);
0284 [vh,vi,xyzr] = simulate_movement(imgs, xyzr);
0285 after = size(xyzr,2);
0286 if(after~=before)
0287 eidors_msg(['mk_GREIT_model: Now using ' num2str(after) ' points']);
0288 end
0289 xyz = xyzr(1:3,:);
0290
0291 function z = calc_offset(z0,opt,Nsim)
0292 if opt.random_offset
0293 l_bnd = opt.target_offset(1);
0294 width = sum(opt.target_offset(1:2));
0295 z = z0 - l_bnd + rand(Nsim,1)*width;
0296 else
0297 z = z0*ones(Nsim,1);
0298 end
0299
0300 function r = calc_radius(R,opt,Nsim)
0301 if opt.random_size
0302 min_sz = opt.target_size(1);
0303 max_sz = opt.target_size(2);
0304 range = max_sz - min_sz;
0305 r = (min_sz + rand(Nsim,1)*range)*R;
0306 else
0307 r = opt.target_size(1)*ones(Nsim,1)*R;
0308 end
0309
0310
0311
0312 function RM = resize_if_reqd(RM,inside,rmdl)
0313 szRM = size(RM,1);
0314 if sum(inside) == szRM || ...
0315 szRM == size(rmdl.elems,1) || ...
0316 (isfield(rmdl,'coarse2fine') && szRM == size(rmdl.coarse2fine,2))
0317
0318 elseif any(size(inside)==szRM) && any(size(inside) == 1)
0319 RM = RM(inside,:);
0320 else
0321 error('mismatch in size of provided RecMatrix');
0322 end
0323
0324
0325 function [imdl,fmdl,imgs] = parse_fmdl(fmdl);
0326 imdl = [];
0327 if ischar(fmdl)
0328 imgs = get_prepackaged_fmdls( fmdl );
0329 fmdl = imgs.fwd_model;
0330 elseif isfield(fmdl,'type');
0331 switch fmdl.type
0332
0333 case 'fwd_model'; imgs = mk_image( fmdl, 1);
0334
0335 case 'image'; imgs = fmdl;
0336 fmdl = imgs.fwd_model;
0337 case 'inv_model'; imdl = fmdl;
0338 fmdl = imdl.fwd_model;
0339 imgs = calc_jacobian_bkgnd(imdl);
0340 otherwise; error('unrecognized eidors object');
0341 end
0342 else
0343 error('specified parameter must be an object or a string');
0344 end
0345
0346 if isempty(imdl)
0347 imdl = select_imdl( fmdl,{'Basic GN dif'});
0348 end
0349
0350
0351 function opt = parse_options(opt,fmdl,imdl, weight)
0352
0353 if ~isfield(opt, 'imgsz'), opt.imgsz = [32 32]; end
0354 if ~isfield(opt, 'square_pixels')
0355 opt.square_pixels = 0;
0356 end
0357
0358 if isfield(imdl,'rec_model') && ~isempty(imdl.rec_model)
0359
0360 opt.imgsz(1) = numel(unique(imdl.rec_model.nodes(:,1)))-1;
0361 opt.imgsz(2) = numel(unique(imdl.rec_model.nodes(:,2)))-1;
0362 try
0363 opt.imgsz(3) = numel(unique(imdl.rec_model.nodes(:,3)))-1;
0364 end
0365 end
0366
0367 if ~isfield(opt, 'distr'), opt.distr = 3; end
0368 if ~isfield(opt, 'Nsim' ), opt.Nsim = 1000; end
0369 if ~isfield(opt, 'noise_figure'), opt.noise_figure = []; end
0370 if ~isfield(opt, 'image_SNR'), opt.image_SNR = []; end
0371 if isempty(opt.noise_figure) && isempty(opt.image_SNR) && isempty(weight)
0372 error('EIDORS:WrongInput', ...
0373 'The weight parameter must be specified if opt.noise_figure or opt.image_SNR are empty or absent');
0374 end
0375 if isfield(opt,'extra_noise')
0376 error('mk_GREIT_model: doesn''t currently support extra_noise');
0377 end
0378 if ~isfield(opt, 'target_size')
0379 opt.target_size = 0.05;
0380 end
0381 if sum(size(opt.target_size)) > 2
0382 if opt.target_size(1) == opt.target_size(2);
0383 opt.random_size = false;
0384 else
0385 opt.random_size = true;
0386 end
0387 end
0388 if sum(size(opt.target_size)) == 2
0389 opt.random_size = false;
0390 end
0391
0392
0393 Nelecs = length(fmdl.electrode);
0394 for i=1:Nelecs
0395 enodesi = fmdl.electrode(i).nodes;
0396 elec_loc(i,:) = mean( fmdl.nodes( enodesi,:),1 );
0397 end
0398 opt.elec_loc = elec_loc;
0399
0400 if ~isfield(opt, 'target_plane')
0401 opt.target_plane = mean(elec_loc(:,3));
0402 else
0403 t = opt.target_plane;
0404 minnode = min(fmdl.nodes);
0405 maxnode = max(fmdl.nodes);
0406 if t<minnode(3) || t>maxnode(3)
0407 warning('options.target_plane is outside the model!');
0408 eidors_msg('mk_GREIT_model: Resorting to default target_plane');
0409 opt.target_plane = mean(elec_loc(:,3));
0410 end
0411 end
0412 if ~isfield(opt, 'target_offset')
0413 opt.target_offset = 0;
0414 end
0415 if sum(size(opt.target_offset)) == 2
0416 if opt.target_offset < 0, opt.target_offset = 0; end
0417 opt.target_offset(2) = opt.target_offset(1);
0418 end
0419 if any(opt.target_offset > 0)
0420 opt.random_offset = true;
0421 else
0422 opt.random_offset = false;
0423 end
0424
0425 if ~isfield(opt,'noise_figure_targets');
0426 R = max(max(fmdl.nodes(:,1:2)) - min(fmdl.nodes(:,1:2)));
0427 xyzr = mean(fmdl.nodes);
0428 xyzr(3) = opt.target_plane;
0429 xyzr(4) = mean(opt.target_size)*0.5*R;
0430 opt.noise_figure_targets = xyzr;
0431 end
0432
0433
0434 if ~isfield(opt,'keep_intermediate_results');
0435 opt.keep_intermediate_results = false;
0436 end
0437
0438
0439 try, opt.normalize = fmdl.normalize_measurements;
0440 catch,
0441 opt.normalize = 0;
0442 eidors_msg('mk_GREIT_model: fmdl.normalize_measurements not specified, assuming 0');
0443 end
0444
0445
0446 slc = mdl_slice_mesher(fmdl,[inf inf opt.target_plane]);
0447 bnd = find_boundary(slc.fwd_model);
0448 opt.contour_boundary = order_loop(slc.fwd_model.nodes(unique(bnd),:));
0449
0450
0451 function do_unit_test
0452
0453 sidx= 1; subplot(4,4,sidx);
0454 do_performance_test;
0455
0456
0457 fmdl_1= ng_mk_ellip_models([1,1.2,0.8],[16,0.5],[0.1]);
0458
0459 extra={'ball','solid ball = sphere(0.5,0.5,0.5;0.1);'};
0460 [fmdl_2,mat_idx]= ng_mk_ellip_models([1,1.2,0.8],[16,0.5],[0.1],extra);
0461
0462 stim = mk_stim_patterns(16,1,[0,1],[0,1],{});
0463 fmdl_1.stimulation = stim;
0464 fmdl_2.stimulation = stim;
0465
0466 img = mk_image(fmdl_2, 0.5); vh = fwd_solve(img);
0467
0468 img.elem_data(mat_idx{2})= 1.0; vi = fwd_solve(img);
0469 sidx= sidx+1; subplot(4,4,sidx);
0470 show_fem(img);
0471
0472 imdl= mk_common_gridmdl('GREITc1');
0473 img= inv_solve(imdl,vh,vi);
0474 sidx= sidx+1; subplot(4,4,sidx);
0475 show_slices(img)
0476
0477
0478 opt.noise_figure = 0.5; opt.distr = 3;opt.square_pixels = 1;
0479 fmdl_2 = mdl_normalize(fmdl_2,0);
0480
0481 img_2 = mk_image(fmdl_2,0.5);
0482 imdl1 = mk_GREIT_model(img_2, 0.25, [], opt);
0483 img1= inv_solve(imdl1,vh,vi);
0484 sidx= sidx+1; subplot(4,4,sidx);
0485 show_slices(img1);
0486
0487
0488 opt = rmfield(opt,'noise_figure');
0489 opt.image_SNR = 1e-3;
0490 weight = 90;
0491 imdl1 = mk_GREIT_model(img_2, 0.25, weight, opt);
0492 img1= inv_solve(imdl1,vh,vi);
0493 sidx= sidx+1; subplot(4,4,sidx); show_slices(img1);
0494
0495 unit_test_cmp('Expect no PJT or X', ~isfield(imdl1.solve_use_matrix, 'PJt') & ...
0496 ~isfield(imdl1.solve_use_matrix, 'X'), true);
0497
0498 weight = [];
0499 opt.keep_intermediate_results = true;
0500 imdl1 = mk_GREIT_model(img_2, 0.25, weight, opt);
0501 img1= inv_solve(imdl1,vh,vi);
0502 sidx= sidx+1; subplot(4,4,sidx); show_slices(img1);
0503
0504 unit_test_cmp('Expect PJT and X', isfield(imdl1.solve_use_matrix, 'PJt') & ...
0505 isfield(imdl1.solve_use_matrix, 'X'), true);
0506
0507 opt = rmfield(opt,{'image_SNR', 'keep_intermediate_results'}); opt.noise_figure = 0.5;
0508
0509
0510 fmdl_1 = mdl_normalize(fmdl_1,0);
0511 imdl2 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0512 img2= inv_solve(imdl2,vh,vi);
0513 sidx= sidx+1; subplot(4,4,sidx); show_slices(img2);
0514
0515
0516
0517 opt.noise_figure_targets = [-.5 0 .5 .2;.5 0 .5 .2;];
0518 imdl3 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0519 img3= inv_solve(imdl3,vh,vi);
0520 sidx= sidx+1; subplot(4,4,sidx); show_slices(img3);
0521
0522 opt = rmfield(opt,'noise_figure_targets');
0523
0524
0525
0526 fmdl_2 = mdl_normalize(fmdl_2,1);
0527
0528 imdl3 = mk_GREIT_model(mk_image(fmdl_2,0.5), 0.25, [], opt);
0529 img3= inv_solve(imdl3,vh,vi);
0530
0531
0532 fmdl_1 = mdl_normalize(fmdl_1,1);
0533 imdl4 = mk_GREIT_model(mk_image(fmdl_1,0.5), 0.25, [], opt);
0534 img4= inv_solve(imdl4,vh,vi);
0535
0536 sidx= sidx+1; subplot(4,4,sidx);
0537 show_slices([img1 img2 img3 img4])
0538
0539
0540
0541 fmdl = mk_library_model('adult_male_16el_lungs');
0542 fmdl.stimulation = stim;
0543 fmdl = mdl_normalize(fmdl,1);
0544 img = mk_image(fmdl,1);
0545 img.elem_data([fmdl.mat_idx{2}; fmdl.mat_idx{3}],1) = 0.3;
0546 vh = fwd_solve(img);
0547 img.elem_data([fmdl.mat_idx{2}; fmdl.mat_idx{3}],1) = 0.4;
0548 vi = fwd_solve(img);
0549
0550
0551 fmdl2 = mk_library_model('adult_male_16el');
0552 fmdl2.stimulation = stim;
0553 fmdl2 = mdl_normalize(fmdl2,1);
0554
0555 opt.imgsz = [50 30];
0556 opt.square_pixels = 1;
0557 imdl = mk_GREIT_model(fmdl2,0.25,3,opt);
0558
0559 img = inv_solve(imdl,vh, vi);
0560 sidx= sidx+1; subplot(4,4,sidx);
0561 show_slices(img);
0562
0563
0564
0565 opt = rmfield(opt,'noise_figure');
0566 opt.image_SNR = 1e-4;
0567 weight = 0.5;
0568 imdl = mk_GREIT_model(fmdl2, 0.25, weight, opt);
0569 img = inv_solve(imdl,vh, vi);
0570 sidx= sidx+1; subplot(4,4,sidx); show_slices(img);
0571
0572 opt.image_SNR_targets = [0.3 0.3 0.5 0.05; 0.3 -0.3 0.5 0.05; ...
0573 0.3 -0.3 0.5 0.05; -0.3 -0.3 0.5 0.05; ...
0574 0.3 0 0.5 0.05; -0.3 0 0.5 0.05; ...
0575 0 0.3 0.5 0.05; 0 -0.3 0.5 0.05]';
0576 opt.image_SNR = 3e-4;
0577 weight = 1E-2;
0578 imdl = mk_GREIT_model(fmdl2, 0.25, weight, opt);
0579 img = inv_solve(imdl,vh, vi);
0580 sidx= sidx+1; subplot(4,4,sidx); show_slices(img);
0581
0582
0583 function do_performance_test
0584
0585 imdl_v1 = mk_common_gridmdl('GREITc1');
0586 imdl_v1.inv_solve.calc_solution_error = false;
0587
0588
0589 imdl_bp = mk_common_gridmdl('backproj');
0590
0591
0592
0593 fmdl = mk_library_model('cylinder_16x1el_fine');
0594 fmdl.nodes = fmdl.nodes/15;
0595 fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{'no_meas_current'}, 1);
0596 opt.noise_figure = 0.88;
0597 opt.target_size = 0.1;
0598 opt.distr = 0;
0599 imdl_gr = mk_GREIT_model(fmdl, 0.2, [], opt);
0600
0601 opt = struct();
0602 opt.noise_figure = 0.5;
0603 imdl_def = mk_GREIT_model(fmdl,0.2,[],opt);
0604
0605 opt.desired_solution_fn = 'GREIT_desired_img_original';
0606 imdl_org = mk_GREIT_model(fmdl,0.2,[],opt);
0607
0608 test_performance( { imdl_v1, imdl_gr, imdl_def, imdl_org},fmdl );