0001 function mdl2 = place_elec_on_surf(mdl,elec_pos, elec_spec,ng_opt_file, maxh)
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 citeme(mfilename);
0052
0053 if ischar(mdl) && strcmp(mdl, 'UNIT_TEST'), do_unit_test; return; end
0054 if nargin < 4
0055 ng_opt_file = '';
0056 end
0057 if nargin < 5
0058 maxh = [];
0059 end
0060
0061 opt.fstr = 'place_elec_on_surf';
0062 mdl2 = eidors_cache(@do_place_elec_on_surf,{mdl,elec_pos, elec_spec,ng_opt_file, maxh},opt);
0063
0064 mdl2 = mdl_normalize(mdl2, 'default');
0065 try name = mdl.name; catch, name = 'Unnamed model'; end
0066 mdl2 = eidors_obj('fwd_model', mdl2, 'name', [name ' + place_elec_on_surf']);
0067
0068
0069 function mdl2 = do_place_elec_on_surf(mdl,elec_pos, elec_spec,ng_opt_file, maxh)
0070
0071
0072
0073 if do_debug; fnstem = 'tmp1';
0074 else; fnstem = tempname;
0075 end
0076
0077 stlfn = [fnstem,'.stl'];
0078 meshfn= [fnstem,'.vol'];
0079
0080 if do_debug; fnstem = 'tmp2';
0081 else; fnstem = tempname;
0082 end
0083
0084 stlfn2 = [fnstem,'.stl'];
0085
0086
0087 mdl = prepare_surf_model(mdl);
0088 if isempty(maxh)
0089 maxh = max(mdl.edge_len);
0090 end
0091 elecs = [];
0092 elecs = parse_elecs(mdl,elec_pos,elec_spec);
0093 if isempty(elecs)
0094 error('EIDORS:WrongInput', 'Failed to parse electrode positions. Exiting');
0095 end
0096
0097
0098
0099 for i = 1:length(elecs)
0100 try
0101 [N, fc] = grow_neighbourhood(mdl,elecs(i));
0102 [mdl, E1{i}, E2{i}, V{i}] = add_electrodes(mdl,N,fc, elecs(i));
0103 catch e
0104 eidors_msg('Failed to add electrode #%d',i,1);
0105 rethrow(e);
0106 end
0107 end
0108
0109
0110 stl_write(mdl,stlfn,'txt');
0111 write_ng_opt_file(ng_opt_file, maxh)
0112 call_netgen(stlfn,meshfn);
0113 delete('ng.opt');
0114
0115
0116 fmdl=ng_mk_fwd_model(meshfn,[],[],[],[]);
0117 mdl = fix_model(fmdl, struct('inner_normal', true, 'boundary_face', true));
0118 mdl = orient_boundary(mdl);
0119 mdl.elems = mdl.boundary;
0120
0121
0122 for i = 1:length(elecs)
0123 try
0124 mdl = flatten_electrode(mdl,E1{i},E2{i}, V{i});
0125 catch e
0126 eidors_msg('Failed to flatten electrode #%d',i,1);
0127 rethrow(e);
0128 end
0129 end
0130
0131
0132 stl_write(mdl,stlfn2);
0133 mdl2 = gmsh_stl2tet(stlfn2, maxh);
0134 mdl2 = fix_boundary(mdl2);
0135 mdl2 = eidors_obj('set',mdl2,'electrode',mdl.electrode);
0136
0137
0138 for i = 1:length(elecs)
0139 enodes = mdl.nodes(mdl.electrode(i).nodes,:);
0140 mdl2.electrode(i).nodes = find_matching_nodes(mdl2,enodes,1e-4);
0141 end
0142
0143
0144 if ~do_debug
0145 delete(stlfn);
0146 delete(meshfn);
0147 delete(stlfn2);
0148 end
0149
0150 function debugging = do_debug
0151 debugging = eidors_debug('query','place_elec_on_surf');
0152
0153 function write_ng_opt_file(ng_opt_file, maxh)
0154
0155
0156 if ~isempty(ng_opt_file)
0157 ng_write_opt(ng_opt_file);
0158 else
0159 opt.meshoptions.fineness = 6;
0160 opt.options.curvaturesafety = 0.2;
0161
0162
0163 opt.stloptions.yangle = 30;
0164
0165 opt.stloptions.edgecornerangle = 0;
0166
0167 opt.stloptions.outerchartangle = 120;
0168 opt.stloptions.resthchartdistenable = 1;
0169 opt.stloptions.resthchartdistfac = 2.0;
0170 if ~isempty(maxh)
0171 opt.options.meshsize = maxh;
0172 end
0173 opt.meshoptions.laststep = 'mv';
0174 opt.options.optsteps2d = 5;
0175 opt.options.badellimit = 120;
0176 ng_write_opt(opt);
0177 end
0178
0179
0180
0181 function mdl = prepare_surf_model(mdl)
0182
0183 mdl.boundary = 1;
0184 opt = struct('inner_normal', true, 'boundary_face', true);
0185 mdl = fix_model(mdl,opt);
0186 mdl.boundary = mdl.faces(mdl.boundary_face,:);
0187 mdl = orient_boundary(mdl);
0188 mdl.elems = mdl.boundary;
0189 mdl.faces = mdl.boundary;
0190 mdl.face_centre = mdl.face_centre(mdl.boundary_face,:);
0191 mdl.normals = mdl.normals(mdl.boundary_face,:);
0192 mdl = rmfield(mdl, {'inner_normal','boundary_face','elem_centre',...
0193 'face2elem','elem2face'});
0194 idx = nchoosek(1:3, 2);
0195 elem_sorted = sort(mdl.elems,2);
0196 [mdl.edges, ~, ~] = unique(reshape(elem_sorted(:,idx),[],2),'rows');
0197 D = mdl.nodes(mdl.edges(:,1),:) - mdl.nodes(mdl.edges(:,2),:);
0198 mdl.edge_len = sqrt(sum(D.^2,2));
0199
0200 function mdl = orient_boundary(mdl)
0201
0202 flip = mdl.elem2face(logical(mdl.boundary_face(mdl.elem2face).*mdl.inner_normal));
0203 mdl.faces(flip,:) = mdl.faces(flip,[1 3 2]);
0204 mdl.normals(flip,:) = -mdl.normals(flip,:);
0205 mdl.boundary = mdl.faces(mdl.boundary_face,:);
0206
0207
0208 function mdl = flatten_electrode(mdl,inner,outer, V)
0209 n1 = find_matching_nodes(mdl,inner, 1e-2);
0210 n2 = find_matching_nodes(mdl,outer, 1e-4);
0211
0212 N1 = false(length(mdl.nodes),1);
0213 N1(n1) = true;
0214 N2 = false(length(mdl.nodes),1);
0215 N2(n2) = true;
0216 rm = sum(N1(mdl.elems),2)>0 & sum(N2(mdl.elems),2)>0;
0217
0218 f = find(sum(N2(mdl.elems),2)>1 & ~rm,1,'first');
0219 B = find(mdl.boundary_face);
0220 p = mdl.face_centre(B(f),:);
0221 r = Inf;
0222 mdl.elems(rm,:) = [];
0223 mdl.boundary = mdl.elems;
0224 mdl.boundary_face(B(rm)) = [];
0225 mdl.face_centre(B(rm),:) = [];
0226 mdl.normals(B(rm),:) = [];
0227 mdl.faces(B(rm),:) = [];
0228 f = f - nnz(rm(1:f));
0229 N = grow_neighbourhood(mdl,f,p,r);
0230
0231
0232
0233
0234 ntm = unique(mdl.elems(N,:));
0235 mdl.nodes(ntm,:) = mdl.nodes(ntm,:) - repmat(V,length(ntm),1);
0236
0237
0238 map = 1:length(mdl.nodes);
0239 map(n2) = n1;
0240 mdl.elems = map(mdl.elems);
0241 mdl.faces = map(mdl.faces);
0242 e_nodes = map(ntm);
0243
0244
0245 m = true(length(mdl.nodes),1);
0246 m(n2) = false;
0247 map = zeros(size(m));
0248 map(m) = 1:nnz(m);
0249
0250 mdl.nodes(n2,:) = [];
0251 mdl.elems = map(mdl.elems);
0252 mdl.faces = map(mdl.faces);
0253 e_nodes = map(e_nodes);
0254
0255 mdl.boundary = mdl.elems;
0256 if ~isfield(mdl,'electrode')
0257 mdl.electrode = struct();
0258 l = 1;
0259 else
0260 l = length(mdl.electrode);
0261
0262
0263 for i = 1:l
0264 mdl.electrode(i).nodes = map(mdl.electrode(i).nodes);
0265 end
0266 l = l + 1;
0267 end
0268 mdl.electrode(l).nodes = double(e_nodes);
0269 mdl.electrode(l).z_contact = 0.01;
0270
0271 if do_debug
0272 show_fem(mdl);
0273
0274 end
0275
0276
0277 function match = find_matching_nodes(mdl, nodes,th)
0278 len = length(mdl.nodes);
0279 match = 0 * (1:length(nodes));
0280 for n = 1:length(nodes)
0281 D = mdl.nodes - repmat(nodes(n,:),len,1);
0282 D = sqrt(sum(D.^2,2));
0283 [val, p] = min(D);
0284 if val < th
0285 match(n) = p;
0286 end
0287 end
0288
0289
0290
0291 function [joint, EL1, EL2, V] = add_electrodes(mdl,N,fc,elecs)
0292
0293
0294
0295
0296
0297
0298 fcs = N;
0299 used_nodes = unique(mdl.faces(fcs,:));
0300 node_map = zeros(1,size(mdl.nodes,1));
0301 node_map(used_nodes) = 1:numel(used_nodes);
0302
0303 jnk.type = 'fwd_model';
0304 jnk.elems = node_map(mdl.boundary(N,:));
0305 jnk.nodes = mdl.nodes(used_nodes,:);
0306 jnk.boundary = jnk.elems;
0307 if do_debug
0308 show_fem(jnk);
0309 hold on
0310 plot3(elecs.points(:,1),elecs.points(:,2),elecs.points(:,3),'ro');
0311
0312 hold off
0313
0314 end
0315
0316 flat = level_model_slice(jnk,struct('centre',mdl.face_centre(fc,:),'normal',mdl.normals(fc,:)));
0317 elec_pts = level_model_slice(elecs.points,struct('centre',mdl.face_centre(fc,:),'normal',mdl.normals(fc,:)));
0318 elec_nodes = level_model_slice(elecs.nodes,struct('centre',mdl.face_centre(fc,:),'normal',mdl.normals(fc,:)));
0319
0320 warning off 'MATLAB:triangulation:PtsNotInTriWarnId'
0321 TR = triangulation(double(flat.elems), flat.nodes(:,1:2));
0322 warning on 'MATLAB:triangulation:PtsNotInTriWarnId'
0323
0324 edges = TR.freeBoundary;
0325
0326
0327 PN = flat.nodes(:,1:2);
0328
0329
0330 neighbour = TR.nearestNeighbor(elec_pts(:,1:2));
0331 D = sqrt(sum((flat.nodes(neighbour,1:2) - elec_pts(:,1:2)).^2,2));
0332 rm = unique(neighbour(D < 2 * elecs.maxh));
0333
0334
0335 b = unique(edges(:));
0336 rm(ismember(rm,b)) = [];
0337
0338
0339 PN(rm,:) = [];
0340 used_nodes(rm) = [];
0341
0342 n = size(flat.nodes,1);
0343 nodelist = 1:n;
0344 nodelist(rm) = [];
0345 map = zeros(n,1);
0346 map(nodelist) = 1:numel(nodelist);
0347 edges = map(edges);
0348
0349 points = [PN; elec_pts(:,1:2)];
0350
0351
0352 f = length(PN) + (1:2);
0353 C = bsxfun(@plus, (0:length(elecs.points)-2)', f);
0354 [wtxt, wid] = lastwarn;
0355 lastwarn('','');
0356 warning off 'MATLAB:DelaunayTri:ConsConsSplitWarnId';
0357 D = DelaunayTri(points,[edges; C]);
0358 els = D.Triangulation(D.inOutStatus,:);
0359 [txt, id] = lastwarn;
0360 if strcmp(id,'MATLAB:DelaunayTri:ConsConsSplitWarnId')
0361 if do_debug
0362 keyboard
0363 else
0364 error(txt);
0365 end
0366 else
0367 lastwarn(wtxt,wid);
0368 end
0369 warning on 'MATLAB:DelaunayTri:ConsConsSplitWarnId';
0370
0371
0372 Ne = mdl.normals(fc,:);
0373 FN = TR.pointLocation(elec_nodes(:,1:2));
0374
0375 for j = 1:length(elecs.nodes)
0376 Pe = elecs.nodes(j,:);
0377 Nf = mdl.normals(fcs(FN(j)),:);
0378 Cf = mdl.face_centre(fcs(FN(j)),:);
0379 Proj(j,:) = Pe + Ne * dot(Cf-Pe,Nf) / dot(Ne,Nf) ;
0380 end
0381
0382
0383
0384 EL1 = Proj(1:length(elecs.points),:);
0385
0386
0387 ln = length(used_nodes);
0388
0389
0390
0391 add = elecs.maxh;
0392
0393 nn = mdl.nodes(used_nodes,:);
0394 le = length(elecs.nodes);
0395 ne = Proj + add * repmat(Ne,le,1);
0396
0397
0398 EL2 = ne(1:length(elecs.points),:);
0399 V = add*Ne;
0400
0401
0402
0403 el_c = D.incenters;
0404 el_c(~D.inOutStatus,:) = [];
0405 e_el = inpolygon(el_c(:,1),el_c(:,2),points(ln+1:end,1),points(ln+1:end,2));
0406 els(e_el,:) = [];
0407
0408
0409 E = [];
0410 le = length(elecs.points);
0411 f = ln + [ 1 le+1 le+2; le+2 2 1];
0412 for j = 0:(le-2)
0413 E = [E; j+f];
0414 end
0415 M = ln + [le+1 le 2*le; le le+1 1];
0416 E = [E; M];
0417
0418 jnk.nodes = [nn ; Proj(1:le,:); ne];
0419 jnk.elems = [ els; E; elecs.elems+ln+le];
0420 jnk.boundary = jnk.elems;
0421 if do_debug
0422 show_fem(jnk);
0423
0424 end
0425
0426
0427 big = mdl;
0428 big.boundary(N,:) = [];
0429 big.faces(N,:) = [];
0430 big.normals(N,:) = [];
0431 big.face_centre(N,:) = [];
0432
0433 big.elems = big.boundary;
0434 log_level = eidors_msg('log_level',1);
0435 joint = merge_meshes(big,jnk,0.001);
0436 eidors_msg('log_level', log_level);
0437 joint.boundary = joint.elems;
0438 joint.faces = joint.boundary;
0439 opt.normals = true;
0440 opt.face_centre = true;
0441 joint = fix_model(joint,opt);
0442
0443
0444 function PN = project_nodes_on_elec(mdl,elecs,nodes)
0445 fc = find_face_under_elec(mdl,elecs.pos);
0446 Ne = mdl.normals(fc,:);
0447 Pe = elecs.pos;
0448
0449
0450
0451
0452 N = mdl.nodes(nodes,:);
0453
0454 PN = N + bsxfun(@times,sum(bsxfun(@times,bsxfun(@minus,Pe,N), Ne),2), Ne);
0455
0456
0457
0458
0459
0460
0461
0462
0463
0464
0465 function [elecs] = parse_elecs(mdl, elec_pos, elec_shape )
0466 elecs = [];
0467
0468 if size(elec_shape,2) < 3
0469 elec_shape(:,3) = elec_shape(:,1)/10;
0470 end
0471
0472 have_xyz = 0;
0473
0474 if size(elec_pos,1) == 1
0475
0476 n_elecs= elec_pos(1);
0477 offset = elec_pos(2) - floor(elec_pos(2));
0478 switch floor(elec_pos(2))
0479 case 0
0480 th = linspace(0,2*pi, n_elecs+1)'; th(end)=[];
0481 th = th + offset*2*pi;
0482 ind = th >= 2*pi;
0483 th(ind) = th(ind) - 2*pi;
0484 case 1
0485 error('not implemented yet');
0486 end
0487 on_elecs = ones(n_elecs, 1);
0488
0489 n_rings = length(elec_pos) - 2;
0490 el_th = repmat(th,n_rings, 1);
0491 el_z = on_elecs * elec_pos(3:end);
0492 el_z = el_z(:);
0493
0494 elseif size(elec_pos,2) == 2
0495
0496 el_th = elec_pos(:,1)*2*pi/360;
0497 el_z = elec_pos(:,2);
0498 elseif size(elec_pos,2) == 3
0499
0500 have_xyz = 1;
0501 el_z = elec_pos(:,3);
0502 end
0503
0504 if ~have_xyz
0505 el_th(el_th>pi) = el_th(el_th>pi) - 2*pi;
0506 el_th(el_th<-pi) = el_th(el_th<-pi) + 2*pi;
0507 end
0508 n_elecs= size(el_z,1);
0509
0510 if size(elec_shape,1) == 1
0511 elec_shape = ones(n_elecs,1) * elec_shape;
0512 end
0513
0514 for i = 1:n_elecs
0515 if ~have_xyz
0516 [fc, elecs(i).pos] = find_elec_centre(mdl,el_th(i),el_z(i));
0517 else
0518 elecs(i).pos = elec_pos(i,:);
0519 end
0520
0521 elecs(i).dims = elec_shape(i,1:2);
0522 elecs(i).dims(elecs(i).dims==0) = [];
0523 elecs(i).maxh = elec_shape(i,3);
0524
0525 if elec_shape(i,2) == 0
0526 elecs(i).shape = 'C';
0527 r = elec_shape(i,1);
0528 n = ceil(2*pi*elec_shape(i,1) / elec_shape(i,3));
0529 t = linspace(0,2*pi,n+1); t(end) = [];
0530 x = r*sin(t); y = r*cos(t);
0531 else
0532 elecs(i).shape = 'R';
0533 height = elec_shape(i,1); width = elec_shape(i,2); d_org = elec_shape(i,3);
0534
0535 d = min( [ d_org , height/5, width/5]);
0536 if d < d_org
0537 elecs(i).maxh = d;
0538 eidors_msg('@@@ Decreased maxh of electrode %d from %f to %f',i,d_org, d,2);
0539 end
0540 nh = ceil(height/d)+1; nw = ceil(width/d)+1;
0541 ph = linspace(-height/2,height/2,nh);
0542 pw = linspace(-width/2,width/2,nw);
0543 y = [ph, ph(end)*ones(1,nw-2), fliplr(ph), ph(1)*ones(1,nw-2)];
0544 x = [pw(1)*ones(1,nh-1), pw, pw(end)*ones(1,nh-2), fliplr(pw(2:end))];
0545
0546
0547
0548
0549
0550
0551
0552
0553
0554
0555
0556
0557
0558
0559
0560 e = tand(0.5)*d;
0561 x = x + e* [0 power(-1,0:nh-3) zeros(1,nw) power(-1,0:nh-3) zeros(1,nw-1)];
0562 y = y + e* [zeros(1,nh) power(-1,0:nw-3) zeros(1,nh) power(-1,0:nw-3)];
0563 end
0564 fc = find_face_under_elec(mdl,elecs(i).pos);
0565 [u, v, s] = get_face_basis(mdl, fc);
0566
0567 np = length(x);
0568
0569 ng_write_opt('meshoptions.fineness',1,'options.meshsize',1.2*elecs(i).maxh);
0570 emdl = ng_mk_2d_model(flipud([x', y']));
0571 x = emdl.nodes(:,1); y = emdl.nodes(:,2);
0572 elecs(i).nodes = ones(size(x)) * elecs(i).pos + x*s + y*v;
0573 elecs(i).elems = emdl.elems(:,[1 3 2]);
0574 elecs(i).points = elecs(i).nodes(1:np,:);
0575
0576
0577 end
0578 delete('ng.opt');
0579
0580
0581 function [u, v, s] = get_face_basis(mdl, fc)
0582 u = mdl.normals(fc,:);
0583
0584
0585 I = eye(3);
0586 for i = 1:3
0587 proj(:,i) = I(:,i) - (dot(I(:,i),u')) * u';
0588 end
0589 norm_proj = vecnorm(proj);
0590 min_norm = min(norm_proj);
0591
0592
0593 if norm_proj(3) ~= min_norm
0594 v = [0 0 1] - dot([0 0 1],u) *u;
0595 else
0596
0597
0598
0599 v = [1 0 0] - dot([1 0 0],u)*u;
0600 end
0601 v = v/norm(v);
0602 s = cross(u,v); s= s/norm(s);
0603
0604
0605 function [fc, pos] = find_elec_centre(mdl, el_th,el_z)
0606 fc = [];
0607 pos = [];
0608
0609 Ctr = mean(mdl.nodes(mdl.boundary,:));
0610 Ctr(3) = el_z;
0611
0612
0613 n_above = mdl.nodes(:,3) >= el_z;
0614 sum_above = sum(n_above(mdl.edges),2) ;
0615 edg = sum_above == 1;
0616
0617
0618 n = unique(mdl.edges(edg,:));
0619 nn = mdl.nodes(n,1:2);
0620 nn = nn - repmat(Ctr(:,1:2),length(nn),1);
0621 th = cart2pol(nn(:,1),nn(:,2));
0622 th(:,2) = 1:length(th);
0623 th = sortrows(th);
0624 idx = find(th(:,1) > el_th,1,'first');
0625 if isempty(idx) || idx == 1
0626 n1 = n(th(1,2));
0627 n2 = n(th(end,2));
0628
0629
0630 ed = edg & sum( (mdl.edges == n1) + (mdl.edges == n2) ,2) > 0;
0631 else
0632
0633
0634
0635
0636 n1 = n(th(idx-1,2));
0637 n2 = n(th(idx, 2));
0638 ed = edg & sum( (mdl.edges == n1) + (mdl.edges == n2) ,2) > 0;
0639 end
0640
0641
0642 el = false(length(mdl.boundary),1);
0643 for i = find(ed)'
0644 n1 = mdl.edges(i,1);
0645 n2 = mdl.edges(i,2);
0646 el = el | sum( (mdl.boundary == n1) + (mdl.boundary == n2), 2) == 2;
0647 end
0648 el = find(el);
0649
0650
0651 [De(1), De(2), De(3)] = pol2cart(el_th,1, 0);
0652 for i = 1:length(el)
0653 Nf = mdl.normals(el(i),:);
0654 Cf = mdl.face_centre(el(i),:);
0655
0656
0657
0658
0659
0660
0661
0662 t = dot(Cf-Ctr,Nf) / dot(De,Nf);
0663 if t < 0, continue, end
0664
0665 X = Ctr + De * t ;
0666 if point_in_triangle(X, mdl.faces(el(i),:), mdl.nodes)
0667 pos = X;
0668 fc = el(i);
0669 break;
0670 end
0671
0672 end
0673 if isempty(pos)
0674 error('Error identifying central face under electrode');
0675 end
0676
0677
0678
0679
0680
0681
0682 function [out, fc] = grow_neighbourhood(mdl, varargin)
0683 use_elec = false;
0684 if length(varargin) == 1
0685 use_elec = true;
0686 elecs = varargin{1};
0687 fc = find_face_under_elec(mdl,elecs.pos);
0688 p = elecs.pos;
0689 switch elecs.shape
0690 case 'R'
0691 r = sqrt(sum(elecs.dims.^2,2));
0692 case 'C'
0693 r = 2 * elecs.dims(1);
0694 end
0695 else
0696 fc = varargin{1};
0697 p = varargin{2};
0698 r = varargin{3};
0699 end
0700
0701 done = false(length(mdl.boundary),1);
0702 todo = false(length(mdl.boundary),1);
0703 todo(fc) = true;
0704 bb = mdl.boundary;
0705 vv = mdl.nodes;
0706
0707
0708 dv = vv - repmat(p,length(vv),1);
0709 nl = mdl.normals;
0710 nl = repmat(nl(fc,:),length(vv),1);
0711 dd = sqrt(sum( (dv - repmat(dot(dv,nl,2),1,3) .* nl).^2,2));
0712 dim = size(bb,2);
0713 first = true;
0714 if use_elec
0715 PN = level_model_slice(mdl.nodes,struct('centre',mdl.face_centre(fc,:),'normal',mdl.normals(fc,:)));
0716 elec_pts = level_model_slice(elecs.points,struct('centre',mdl.face_centre(fc,:),'normal',mdl.normals(fc,:)));
0717 elec_pts = elec_pts(:,1:2);
0718 PN = PN(:,1:2);
0719 emin = min(elec_pts);
0720 emax = max(elec_pts);
0721 rng = emax-emin;
0722 emin = emin - 0.1*rng;
0723 emax = emax + 0.1*rng;
0724 toofar = false(size(mdl.boundary,1),1);
0725
0726 for i = 1:2
0727 nodes = reshape(PN(mdl.boundary,i),[],3);
0728 toofar = toofar | sum(nodes > emax(i),2) == 3 | sum(nodes < emin(i),2) == 3;
0729 end
0730 end
0731 near_nodes = any(dd(bb) <= r,2);
0732 while any(todo)
0733 id = find(todo,1,'first');
0734 done(id) = 1;
0735 nn = find_neighbours(id,bb);
0736 if use_elec
0737 nn = nn & ~toofar;
0738 elseif first
0739
0740 first = false;
0741 else
0742
0743 nn = nn & near_nodes;
0744 end
0745 todo = todo | nn;
0746 todo(done) = 0;
0747
0748
0749
0750 end
0751 out = find(done);
0752
0753
0754 function nn = find_neighbours(fc, bb)
0755
0756 dim = size(bb,2);
0757 nn = false(length(bb),1);
0758 for i = 1:dim
0759 node = bb(fc,i);
0760 nn = nn | sum(bb == node,2) > 0;
0761 end
0762 nn(fc) = 0;
0763
0764
0765 function [e, p] = find_face_under_elec(mdl, elec_pos)
0766
0767 for i = 1:size(elec_pos,1)
0768
0769 ee = repmat(elec_pos(i,:),length(mdl.faces),1);
0770 fc = mdl.face_centre;
0771 n = mdl.normals;
0772 proj1 = ee - repmat(dot(ee-fc, n,2),1,3) .* n;
0773 in1 = point_in_triangle(proj1,mdl.faces,mdl.nodes, 'match');
0774 dis1 = sqrt(sum((ee-proj1).^2,2));
0775
0776 edg = [mdl.faces(:,1:2);mdl.faces(:,2:3);mdl.faces(:,[3 1])];
0777 edg = sort(edg,2);
0778 [edg, ~, e2f] = unique(edg,'rows');
0779 ee = repmat(elec_pos(i,:),length(edg),1);
0780 s = mdl.nodes(edg(:,2),:) - mdl.nodes(edg(:,1),:);
0781 t = dot(ee-mdl.nodes(edg(:,1),:),s,2)./dot(s,s,2);
0782 in2 = t>=0 & t <=1;
0783 in2 = any(reshape(in2(e2f),[],3),2);
0784 proj2 = mdl.nodes(edg(:,1),:) + repmat(t,1,3).*s;
0785 dis = sqrt(sum((ee - proj2).^2,2));
0786 dis = repmat(dis,2,1);
0787 dis(t<0 | t > 1) = Inf;
0788 dis = reshape(dis(e2f),[],3);
0789 [~, pos] = min(dis,[],2);
0790 idx = sparse(1:length(pos),pos,1);
0791 dis = dis';
0792 dis2 = dis(logical(idx'));
0793
0794 in = in1 | in2;
0795 if nnz(in) == 1
0796 e(i) = find(in1);
0797 p(i,:) = proj1(in1,:);
0798 else
0799
0800 cand = find(in);
0801
0802 dd(in1(cand)) = dis1(in1);
0803 dd(in2(cand)) = dis2(in2);
0804 [~, pos] = min(dd);
0805 e(i) = cand(pos);
0806 p(i,:) = proj1(e(i),:);
0807 end
0808
0809 end
0810
0811
0812 function do_unit_test
0813 xy= [ -0.89 -0.74 -0.21 0.31 0.79 0.96 0.67 0.05 -0.36 -0.97;
0814 0.14 0.51 0.35 0.50 0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0815 if exist('ng.opt','file')
0816 delete ng.opt
0817 end
0818 [fmdl] = ng_mk_extruded_model({2,xy,[4,80],},[],[]);
0819
0820
0821
0822 mdl = place_elec_on_surf(fmdl, [16 0 1],[0.15 0.1 0.01]);
0823
0824 subplot(121)
0825 show_fem(mdl);
0826
0827 mdl = place_elec_on_surf(fmdl, [16 0 1], [0.15 0.1 0.01],[],0.1);
0828
0829 subplot(122)
0830 show_fem(mdl);