place_elec_on_surf

PLACE_ELEC_ON_SURF Place electrodes on the surface of a model
 mdl = place_elec_on_surf(mdl,elec_pos, elec_spec, ng_opt_file, maxh)
 INPUT:
  mdl         = an EIDORS fwd_model struct
  elec_pos    = an array specigying electrode positions
  elec_shape  = an array specifying electrode shape (can be different for
                each electrode)
  ng_opt_file = an alternative ng.opt file to use (OPTIONAL)
                filename or option struct (see ng_write_opt)
                specify [] to use dafault
  maxh        = maximum edge length (if ng_opt_file is specified, maxh 
                only applies to the volume and not the surface)
 ELECTRODE POSITIONS:
  elec_pos = [n_elecs_per_plane,(0=equal angles,1=equal dist),z1, z2, ...]
     OR
  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
     OR
  elec_pos = [x y z] centres of each electrode (N_elecs x 3)

  Note: N_elecs >= 2. Equal distance electrode spacing not implemented
  yet.

 ELECTRODE SHAPES::
  elec_shape = [width,height, maxsz]  % Rectangular elecs
     OR
  elec_shape = [radius, 0, maxsz ]    % Circular elecs

 NOTE that this function requires both Netgen and Gmsh.
 It will completely re-mesh your model.
 The code makes several assumptions about the output of Netgen, which it
 attempts to control through the ng.opt file, but there will be meshes 
 for which this appraoch will fail. In this case, you can supply your own 
 file with options for Netgen (with a filename different than ng.opt), or
 change your mesh and/or electrode locations. Most common problem is too 
 big electrode maxh value (must be significantly smaller than the smallest
 element on which the electrode will fall).

 CITATION_REQUEST:
 TITLE: FEM Electrode Refinement for Electrical Impedance Tomography 
 AUTHOR: B Grychtol and A Adler
 JOURNAL: Engineering in Medicine and Biology Society (EMBC), 2013 Annual 
 International Conference of the IEEE 
 YEAR: 2013

 See also gmsh_stl2tet, ng_write_opt, merge_meshes

0001 function mdl2 = place_elec_on_surf(mdl,elec_pos, elec_spec,ng_opt_file, maxh)
0002 %PLACE_ELEC_ON_SURF Place electrodes on the surface of a model
0003 % mdl = place_elec_on_surf(mdl,elec_pos, elec_spec, ng_opt_file, maxh)
0004 % INPUT:
0005 %  mdl         = an EIDORS fwd_model struct
0006 %  elec_pos    = an array specigying electrode positions
0007 %  elec_shape  = an array specifying electrode shape (can be different for
0008 %                each electrode)
0009 %  ng_opt_file = an alternative ng.opt file to use (OPTIONAL)
0010 %                filename or option struct (see ng_write_opt)
0011 %                specify [] to use dafault
0012 %  maxh        = maximum edge length (if ng_opt_file is specified, maxh
0013 %                only applies to the volume and not the surface)
0014 % ELECTRODE POSITIONS:
0015 %  elec_pos = [n_elecs_per_plane,(0=equal angles,1=equal dist),z1, z2, ...]
0016 %     OR
0017 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0018 %     OR
0019 %  elec_pos = [x y z] centres of each electrode (N_elecs x 3)
0020 %
0021 %  Note: N_elecs >= 2. Equal distance electrode spacing not implemented
0022 %  yet.
0023 %
0024 % ELECTRODE SHAPES::
0025 %  elec_shape = [width,height, maxsz]  % Rectangular elecs
0026 %     OR
0027 %  elec_shape = [radius, 0, maxsz ]    % Circular elecs
0028 %
0029 % NOTE that this function requires both Netgen and Gmsh.
0030 % It will completely re-mesh your model.
0031 % The code makes several assumptions about the output of Netgen, which it
0032 % attempts to control through the ng.opt file, but there will be meshes
0033 % for which this appraoch will fail. In this case, you can supply your own
0034 % file with options for Netgen (with a filename different than ng.opt), or
0035 % change your mesh and/or electrode locations. Most common problem is too
0036 % big electrode maxh value (must be significantly smaller than the smallest
0037 % element on which the electrode will fall).
0038 %
0039 % CITATION_REQUEST:
0040 % TITLE: FEM Electrode Refinement for Electrical Impedance Tomography
0041 % AUTHOR: B Grychtol and A Adler
0042 % JOURNAL: Engineering in Medicine and Biology Society (EMBC), 2013 Annual
0043 % International Conference of the IEEE
0044 % YEAR: 2013
0045 %
0046 % See also gmsh_stl2tet, ng_write_opt, merge_meshes
0047 
0048 % (C) Bartlomiej Grychtol and Andy Adler, 2012-2013. Licence: GPL v2 or v3
0049 % $Id: place_elec_on_surf.m 6986 2024-11-14 20:26:39Z bgrychtol $
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 % filenames
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 % 1. Get a surface model
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 % 2. Add extruded electrodes
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 % 3. Save as STL and mesh with NETGEN
0110 stl_write(mdl,stlfn,'txt'); %netgen only works well with txt files
0111 write_ng_opt_file(ng_opt_file, maxh)
0112 call_netgen(stlfn,meshfn);
0113 delete('ng.opt'); % clean up
0114 
0115 % 4. Extract surface
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 % 5. One by one, flatten the electrodes
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 % 6. Keeping the surface intact, remesh the inside
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); % standard field order
0136 
0137 % 7. Find all electrode nodes
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 % 8. Clean up
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 % these options are meant to ensure that the electrode sides don't get
0155 % modified, but there's no guarantee
0156 if ~isempty(ng_opt_file)
0157    ng_write_opt(ng_opt_file);
0158 else
0159    opt.meshoptions.fineness = 6; % some options have no effect without this
0160    opt.options.curvaturesafety = 0.2;
0161    % small yangle preserves the original mesh, large encourages smoother
0162    % surface with nicer spreading of refinement
0163    opt.stloptions.yangle = 30; % was 10
0164  %    opt.stloptions.contyangle = 20;
0165    opt.stloptions.edgecornerangle = 0;
0166 %    opt.stloptions.chartangle = 0;
0167    opt.stloptions.outerchartangle = 120;
0168    opt.stloptions.resthchartdistenable = 1;
0169    opt.stloptions.resthchartdistfac = 2.0; % encourages slower increase of element size
0170    if ~isempty(maxh)
0171       opt.options.meshsize = maxh;
0172    end
0173    opt.meshoptions.laststep = 'mv'; % don't need volume optimization
0174    opt.options.optsteps2d =  5; % but we can up surface optimization
0175    opt.options.badellimit = 120; % decrease the maximum allowed angle
0176    ng_write_opt(opt);
0177 end
0178 
0179 
0180 % Extract a nice surface model from the one given
0181 function mdl = prepare_surf_model(mdl)
0182 % ensure mdl has a boundary, so fix_model doesn't recalculate
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 % consistently orient boundary elements
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); % how to figure out good values?
0211 % remove the side nodes of the electrode
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 % WARNING: Here we assume the sides of the electrode are one element high!
0232 
0233 %nodes to move
0234 ntm = unique(mdl.elems(N,:));
0235 mdl.nodes(ntm,:) = mdl.nodes(ntm,:) - repmat(V,length(ntm),1);
0236 
0237 %remap outer nodes to inner ones
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 % remove the outer nodes
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    % because we are changing the number of nodes, we need to correct the
0262    % electrodes that are there already
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 %    keyboard
0274 end
0275 
0276 %Find mdl.nodes less then th away from nodes
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 % Returns a joint surface mesh and the list of nodes on the side of the
0290 % electrode
0291 function [joint, EL1, EL2, V] = add_electrodes(mdl,N,fc,elecs)
0292 
0293 
0294 % fc = find_face_under_elec(mdl,elecs.pos);
0295 % N indexes the boundary, need index into faces
0296 % fcs = find(mdl.boundary_face);
0297 % fcs = fcs(N);
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    % plot3(mdl.nodes(nn(outer),1), mdl.nodes(nn(outer),2), mdl.nodes(nn(outer),3),'bs')
0312    hold off
0313 %    keyboard
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 % now, the points are almost z = 0, so we can work in 2D
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 % project all nodes of the faces in N onto the plane of the electrode
0327 PN = flat.nodes(:,1:2);
0328 
0329 % for every electrode point, find closest node
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 % we can only delete if it's not part of the boundary
0335 b = unique(edges(:));
0336 rm(ismember(rm,b)) = [];
0337 
0338 % remove and remap
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 % constrained Delaunay triangulation in 2D
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); % no point continuing
0365     end
0366 else
0367     lastwarn(wtxt,wid); % restore
0368 end
0369 warning on 'MATLAB:DelaunayTri:ConsConsSplitWarnId';
0370 % project all electrode points on the faces below them, using the normal of
0371 % the central face
0372 Ne = mdl.normals(fc,:);
0373 FN = TR.pointLocation(elec_nodes(:,1:2)); % face num under each electrode point
0374 %FC = fcs(FN); % same, in original numbering
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 % this is just output
0384 EL1 = Proj(1:length(elecs.points),:);
0385 
0386 % remove any nodes inside the electrode
0387 ln = length(used_nodes);
0388 % IN = inpolygon(x(1:ln),y(1:ln),x(ln+1:end),y(ln+1:end));
0389 % nodes(IN) = [];
0390 
0391 add = elecs.maxh;
0392 
0393 nn = mdl.nodes(used_nodes,:);% + add * repmat(IN,1,3) .* repmat(Ne,ln,1);
0394 le = length(elecs.nodes);
0395 ne = Proj + add * repmat(Ne,le,1);
0396 
0397 %this is just output
0398 EL2 = ne(1:length(elecs.points),:);
0399 V = add*Ne;
0400 
0401 % the nodes of the electrode
0402 % IN = [IN; ones(le,1)];
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,:) = []; % els(e_el,:) + (els(e_el,:)>ln ) .* le;
0407 
0408 % add connecting elements
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 %    keyboard
0424 end
0425 
0426 % remove the patch we're replacing
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); % don't want progress messages
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 % Not used now, but needed for some disabled options
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 % for i = 1:length(nodes)
0449 %    P = mdl.nodes(nodes(i),:);
0450 %    PN(i,:) = P + dot(Pe - P, Ne) * Ne;
0451 % end
0452 N = mdl.nodes(nodes,:);
0453 % PN = N + sum((Pe-N) .* Ne,2) .* Ne;
0454 PN = N + bsxfun(@times,sum(bsxfun(@times,bsxfun(@minus,Pe,N), Ne),2), Ne);
0455 
0456 
0457 
0458 % OUTPUT:
0459 %  elecs(i).pos   = [x,y,z]
0460 %  elecs(i).shape = 'C' or 'R'
0461 %  elecs(i).dims  = [radius] or [width,height]
0462 %  elecs(i).maxh  = '-maxh=#' or '';
0463 %  elecs(i).points= list of points around the perimeter
0464 % Angles (th) are interpreted with the mean of boundary nodes as origin
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    % Parse elec_pos = [n_elecs_per_plane,(0=equal angles,1=equal dist),z1, z2, ...]
0476    n_elecs= elec_pos(1); % per plane
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    % support for multiple z-values
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    % elec_pos = [theta z];
0496    el_th = elec_pos(:,1)*2*pi/360;
0497    el_z  = elec_pos(:,2);
0498 elseif size(elec_pos,2) == 3
0499    % elec_pos = [x y z];
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 %    elecs(i).face = fc; % this changes too often to store!
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       % enforce a minimum of 5 nodes per side
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       %    % we don't want real rectangles, because Netgen will merge coplanar
0546       %    % faces, so we create a nice superellipse instead
0547       %    n = 2*(nh+nw);
0548       %    t = linspace(2*pi,0,n); t(end) = [];
0549       %    N = 8;
0550       %    x = abs(cos(t)).^(2/N) * width/2  .* sign(cos(t));
0551       %    y = abs(sin(t)).^(2/N) * height/2 .* sign(sin(t));
0552       % superellipses are also bad, what about a wavy rectange?
0553       %       [pp] = fourier_fit([x; y]', min(size(x,2),18) );
0554       %       t = linspace(0,1,n+1); t(end) = [];
0555       %       xy = fourier_fit(pp,t);
0556       %       x = xy(:,1)'; y = xy(:,2)';
0557       % wavy rectangles are nice but don't guarantee absence of co-planar
0558       % faces
0559       % let's try a brute-force approach
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]); % flip orientation to the outside
0574    elecs(i).points = elecs(i).nodes(1:np,:); % this must be the boundary
0575    % TODO: write code to check if this is true
0576    
0577 end
0578 delete('ng.opt');
0579 
0580 % Basis vectors of a single face
0581 function [u, v, s] = get_face_basis(mdl, fc)
0582    u = mdl.normals(fc,:); % unit normal
0583 
0584    % project each coordinate axis on the plane
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    % vertical vector on the plane of that surface triangle
0593    if norm_proj(3) ~= min_norm
0594       v = [0 0 1] - dot([0 0 1],u) *u;
0595    else
0596       % the element is essentially horizontal
0597 %       v = [0 1 0] - dot([0 1 0],u)*u;
0598 %TODO: need to expose an option to decide which it should be
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 % Find the center of the electrode on the model and the face it falls on
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 %1. Find edges that cross the z plane
0613 n_above = mdl.nodes(:,3) >= el_z;
0614 sum_above = sum(n_above(mdl.edges),2) ;
0615 edg = sum_above == 1;
0616 
0617 %2. Find an edge that crosses el_th
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    % edges in edg that contain these nodes (they don't need to be on the
0629    % same element)
0630    ed = edg & sum( (mdl.edges == n1) + (mdl.edges == n2) ,2) > 0;
0631 else
0632 %    to_the_left = false(length(mdl.nodes),1);
0633 %    to_the_left(n(th(1:idx-1,2))) = true;
0634 %    sum_left = sum( to_the_left(mdl.boundary), 2);
0635 %    el = els & sum_left > 0 & sum_left < 3;
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 % Find elements that contain that edge
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 % Find face under elec center
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    % the plane is (X - Cf).Nf = 0
0656    % the line is X = Ctr + tDe (through Ctr along De
0657    % We want X that satisfies both.
0658    % (Ctr +tDe -  Cf).Nf = 0
0659    % (Ctr - Cf).Nf + tDe.Nf = 0
0660    % t =
0661    % X = Ctr + De * (Cf-Ctr).Nf / (De.Nf)
0662    t = dot(Cf-Ctr,Nf) / dot(De,Nf);
0663    if t < 0, continue, end
0664    % check if it falls inside
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 %grow_neighbourhood(mdl, elecs_struct)
0679 %grow_neighbourhood(mdl, fc, pos, radius)
0680 % Find faces with nodes within radius of the electrode center without
0681 % computing the distance to all nodes (expensive)
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 % distance of each vertex to the line perpendicular to face fc passing
0707 % through p
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; % at first iteration, add all neighbours
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       % include all neighbours
0740       first = false;
0741    else
0742       % at least one node must be close enough
0743       nn = nn & near_nodes;
0744    end
0745    todo = todo | nn;
0746    todo(done) = 0;
0747 %    disp(sprintf('id: %d done: %d todo: %d',id, nnz(done),nnz(todo)));
0748 %    disp(find(todo)');
0749 %    disp(find(done)');
0750 end
0751 out = find(done);
0752 
0753 % Find neighbouring faces, return logical vector
0754 function nn =  find_neighbours(fc, bb)
0755 % nn = any(ismember(bb,bb(fc,:)),2); % <- slower
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; % sum is faster than any
0761 end
0762 nn(fc) = 0;
0763 
0764 % Find the face under each elec_pos coordinate
0765 function [e, p] = find_face_under_elec(mdl, elec_pos)
0766 
0767 for i = 1:size(elec_pos,1)
0768    % 1. Project electrode on all faces
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    % 2. Project electrode on all edges
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),:); %edge direction vector
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);  % this should be an index into mdl.boundary
0797          p(i,:) = proj1(in1,:);
0798    else
0799       % take the element that is closest to ee
0800       cand = find(in);
0801       % pre-allocation is not worth it
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 % elec_pos = [-0.5, -0.8, 1; 0.5, -0.8, 1];
0820 % place_elec_on_surf(fmdl, elec_pos, [0.1 0 0.01]);
0821 % mdl = place_elec_on_surf(fmdl, elec_pos, [0.15 0.1 0.01]);
0822 mdl = place_elec_on_surf(fmdl, [16 0 1],[0.15 0.1 0.01]);
0823 % place_elec_on_surf(fmdl, [16 0 1], [0.1 0 0.01]);
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 % place_elec_on_surf(fmdl, [16 0 1], [0.1 0 0.01]);
0829 subplot(122)
0830 show_fem(mdl);