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)
                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,z_planes] 
     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.

 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 %                specify [] to use dafault
0011 %  maxh        = maximum edge length (if ng_opt_file is specified, maxh
0012 %                only applies to the volume and not the surface)
0013 % ELECTRODE POSITIONS:
0014 %  elec_pos = [n_elecs_per_plane,z_planes]
0015 %     OR
0016 %  elec_pos = [degrees,z] centres of each electrode (N_elecs x 2)
0017 %     OR
0018 %  elec_pos = [x y z] centres of each electrode (N_elecs x 3)
0019 %
0020 %  Note: N_elecs >= 2.
0021 %
0022 % ELECTRODE SHAPES::
0023 %  elec_shape = [width,height, maxsz]  % Rectangular elecs
0024 %     OR
0025 %  elec_shape = [radius, 0, maxsz ]    % Circular elecs
0026 %
0027 % NOTE that this function requires both Netgen and Gmsh.
0028 % It will completely re-mesh your model.
0029 % The code makes several assumptions about the output of Netgen, which it
0030 % attempts to control through the ng.opt file, but there will be meshes
0031 % for which this appraoch will fail. In this case, you can supply your own
0032 % file with options for Netgen (with a filename different than ng.opt), or
0033 % change your mesh and/or electrode locations. Most common problem is too
0034 % big electrode maxh value (must be significantly smaller than the smallest
0035 % element on which the electrode will fall).
0036 %
0037 % CITATION_REQUEST:
0038 % TITLE: FEM Electrode Refinement for Electrical Impedance Tomography
0039 % AUTHOR: B Grychtol and A Adler
0040 % JOURNAL: Engineering in Medicine and Biology Society (EMBC), 2013 Annual
0041 % International Conference of the IEEE
0042 % YEAR: 2013
0043 %
0044 % See also gmsh_stl2tet, ng_write_opt, merge_meshes
0045 
0046 % (C) Bartlomiej Grychtol and Andy Adler, 2012-2013. Licence: GPL v2 or v3
0047 % $Id: place_elec_on_surf.m 6352 2022-04-26 20:18:45Z bgrychtol $
0048 
0049 citeme(mfilename);
0050 
0051 if ischar(mdl) && strcmp(mdl, 'UNIT_TEST') do_unit_test; return; end;
0052 if nargin < 4
0053    ng_opt_file = '';
0054 end
0055 if nargin < 5
0056    maxh = [];
0057 end
0058 
0059 opt.fstr = 'place_elec_on_surf';
0060 mdl2 = eidors_cache(@do_place_elec_on_surf,{mdl,elec_pos, elec_spec,ng_opt_file, maxh},opt);
0061 
0062 
0063 
0064 function mdl2 = do_place_elec_on_surf(mdl,elec_pos, elec_spec,ng_opt_file, maxh)
0065 
0066 
0067 % filenames
0068 if do_debug; fnstem = 'tmp1';
0069 else;        fnstem = tempname;
0070 end
0071 
0072 stlfn = [fnstem,'.stl'];
0073 meshfn= [fnstem,'.vol'];
0074 
0075 if do_debug; fnstem = 'tmp2';
0076 else;        fnstem = tempname;
0077 end
0078 
0079 stlfn2 = [fnstem,'.stl'];
0080 
0081 % 1. Get a surface model
0082 mdl = prepare_surf_model(mdl);
0083 if isempty(maxh)
0084    maxh = max(mdl.edge_len);
0085 end
0086 
0087 elecs = parse_elecs(mdl,elec_pos,elec_spec);
0088 if isempty(elecs)
0089    error('EIDORS:WrongInput', 'Failed to parse electrode positions. Exiting');
0090 end
0091 
0092 
0093 % 2. Add extruded electrodes
0094 for i = 1:length(elecs)
0095    try
0096       [N, fc] = grow_neighbourhood(mdl,elecs(i));
0097       [mdl E1{i} E2{i} V{i}] = add_electrodes(mdl,N,fc, elecs(i));
0098    catch e
0099       eidors_msg('Failed to add electrode #%d',i,1);
0100       rethrow(e);
0101    end
0102 end
0103 
0104 % 3. Save as STL and mesh with NETGEN
0105 write_to_stl(mdl,stlfn);
0106 write_ng_opt_file(ng_opt_file, maxh)
0107 call_netgen(stlfn,meshfn);
0108 delete('ng.opt'); % clean up
0109 
0110 % 4. Extract surface
0111 fmdl=ng_mk_fwd_model(meshfn,[],[],[],[]);
0112 mdl = fix_model(fmdl);
0113 mdl = orient_boundary(mdl);
0114 mdl.elems = mdl.boundary;
0115 
0116 % 5. One by one, flatten the electrodes
0117 for i = 1:length(elecs)
0118    try 
0119       mdl = flatten_electrode(mdl,E1{i},E2{i}, V{i});
0120    catch e
0121       eidors_msg('Failed to flatten electrode #%d',i,1);
0122       rethrow(e);
0123    end
0124 end
0125 
0126 % 6. Keeping the surface intact, remesh the inside
0127 write_to_stl(mdl,stlfn2);
0128 mdl2 = gmsh_stl2tet(stlfn2, maxh);
0129 mdl2.electrode = mdl.electrode;
0130 
0131 % 7. Find all electrode nodes
0132 for i = 1:length(elecs)
0133    enodes = mdl.nodes(mdl.electrode(i).nodes,:);
0134    mdl2.electrode(i).nodes = find_matching_nodes(mdl2,enodes,1e-5);
0135 end
0136 
0137 function debugging = do_debug;
0138   debugging = eidors_debug('query','place_elec_on_surf');
0139 
0140 function write_to_stl(mdl,stlfn)
0141 STL.vertices = mdl.nodes;
0142 STL.faces    = mdl.elems;
0143 stl_write(STL,stlfn);
0144 
0145 function write_ng_opt_file(ng_opt_file, maxh)
0146 % these options are meant to insure that the electrode sides don't get
0147 % modified, but there's no guarantee
0148 if ~isempty(ng_opt_file)
0149    ng_write_opt(ng_opt_file);
0150 else
0151    opt.meshoptions.fineness = 6; % some options have no effect without this
0152    opt.options.curvaturesafety = 0.2;
0153    % small yangle preserves the original mesh, large encourages smoother
0154    % surface with nicer spreading of refinement
0155    opt.stloptions.yangle = 30; % was 10
0156  %    opt.stloptions.contyangle = 20;
0157    opt.stloptions.edgecornerangle = 0;
0158 %    opt.stloptions.chartangle = 0;
0159    opt.stloptions.outerchartangle = 120;
0160    opt.stloptions.resthchartdistenable = 1;
0161    opt.stloptions.resthchartdistfac = 2.0; % encourages slower increase of element size
0162    if ~isempty(maxh)
0163       opt.options.meshsize = maxh;
0164    end
0165    opt.meshoptions.laststep = 'mv'; % don't need volume optimization
0166    opt.options.optsteps2d =  5; % but we can up surface optimization
0167    opt.options.badellimit = 120; % decrease the maximum allowed angle
0168    ng_write_opt(opt);
0169 end
0170 
0171 
0172 % Extract a nice surface model from the one given
0173 function mdl = prepare_surf_model(mdl)
0174 try mdl = rmfield(mdl,'boundary');  end
0175 mdl = fix_model(mdl);
0176 mdl = orient_boundary(mdl);
0177 mdl.elems = mdl.boundary;
0178 mdl.faces = mdl.boundary;
0179 mdl.face_centre = mdl.face_centre(mdl.boundary_face,:);
0180 mdl.normals = mdl.normals(mdl.boundary_face,:);
0181 mdl = rmfield(mdl, 'inner_normal');
0182 mdl = rmfield(mdl, 'boundary_face');
0183 idx = nchoosek(1:3, 2);
0184 elem_sorted = sort(mdl.elems,2);
0185 [mdl.edges ib ia] = unique(reshape(elem_sorted(:,idx),[],2),'rows');
0186 D = mdl.nodes(mdl.edges(:,1),:) - mdl.nodes(mdl.edges(:,2),:);
0187 mdl.edge_len = sqrt(sum(D.^2,2)); 
0188 
0189 function mdl = orient_boundary(mdl)
0190 % consistently orient boundary elements
0191 flip = mdl.elem2face(logical(mdl.boundary_face(mdl.elem2face).*mdl.inner_normal));
0192 mdl.faces(flip,:) = mdl.faces(flip,[1 3 2]);
0193 mdl.normals(flip,:) = -mdl.normals(flip,:);
0194 mdl.boundary = mdl.faces(mdl.boundary_face,:);
0195 
0196 
0197 function mdl = flatten_electrode(mdl,inner,outer, V)
0198 n1 = find_matching_nodes(mdl,inner, 1e-2);
0199 n2 = find_matching_nodes(mdl,outer, 1e-5);
0200 % remove the side nodes of the electrode
0201 N1 = false(length(mdl.nodes),1);
0202 N1(n1) = true;
0203 N2 = false(length(mdl.nodes),1);
0204 N2(n2) = true;
0205 rm = sum(N1(mdl.elems),2)>0 & sum(N2(mdl.elems),2)>0;
0206 
0207 f = find(sum(N2(mdl.elems),2)>1 & ~rm,1,'first');
0208 B = find(mdl.boundary_face);
0209 p = mdl.face_centre(B(f),:);
0210 r = Inf;
0211 mdl.elems(rm,:) = [];
0212 mdl.boundary = mdl.elems;
0213 mdl.boundary_face(B(rm)) = [];
0214 mdl.face_centre(B(rm),:) = [];
0215 mdl.normals(B(rm),:)     = [];
0216 mdl.faces(B(rm),:)       = [];
0217 f = f - nnz(rm(1:f));
0218 N = grow_neighbourhood(mdl,f,p,r);
0219 
0220 % WARNING: Here we assume the sides of the electrode are one element high!
0221 
0222 %nodes to move
0223 ntm = unique(mdl.elems(N,:));
0224 mdl.nodes(ntm,:) = mdl.nodes(ntm,:) - repmat(V,length(ntm),1);
0225 e_nodes = ntm;
0226 
0227 %remap outer nodes to inner ones
0228 map = 1:length(mdl.nodes);
0229 map(n2) = n1;
0230 mdl.elems = map(mdl.elems);
0231 mdl.faces = map(mdl.faces);
0232 e_nodes = map(ntm);
0233 
0234 % remove the outer nodes
0235 m = true(length(mdl.nodes),1);
0236 m(n2) = false;
0237 map = zeros(size(m));
0238 map(m) = 1:nnz(m);
0239 
0240 mdl.nodes(n2,:) = [];
0241 mdl.elems = map(mdl.elems);
0242 mdl.faces = map(mdl.faces);
0243 e_nodes = map(e_nodes);
0244 
0245 mdl.boundary = mdl.elems;
0246 if ~isfield(mdl,'electrode')
0247    mdl.electrode = struct();
0248    l = 1;
0249 else
0250    l = length(mdl.electrode);
0251    % because we are changing the number of nodes, we need to correct the
0252    % electrodes that are there already
0253    for i = 1:l
0254       mdl.electrode(i).nodes = map(mdl.electrode(i).nodes);
0255    end
0256    l = l + 1;
0257 end
0258 mdl.electrode(l).nodes = double(e_nodes);
0259 mdl.electrode(l).z_contact = 0.01;
0260 
0261 if do_debug
0262    show_fem(mdl);
0263 end
0264 
0265 
0266 function match = find_matching_nodes(mdl, nodes,th)
0267 l0 = length(mdl.nodes);
0268 match = 0 * (1:length(nodes));
0269 for n = 1:length(nodes)
0270    D = mdl.nodes - repmat(nodes(n,:),l0,1);
0271    D = sqrt(sum(D.^2,2));
0272    [val p] = min(D);
0273    if val < th
0274       match(n) = p;
0275    end
0276 end
0277 
0278 % Returns a joint surface mesh and the list of nodes on the side of the
0279 % electrode
0280 function [joint EL1 EL2 V] = add_electrodes(mdl,N,fc,elecs)
0281 
0282 
0283 % fc = find_face_under_elec(mdl,elecs.pos);
0284 % N indexes the boundary, need index into faces
0285 % fcs = find(mdl.boundary_face);
0286 % fcs = fcs(N);
0287 fcs = N;
0288 used_nodes = unique(mdl.faces(fcs,:));
0289 node_map = zeros(1,size(mdl.nodes,1));
0290 node_map(used_nodes) = 1:numel(used_nodes);
0291 
0292 jnk.type = 'fwd_model';
0293 jnk.elems = node_map(mdl.boundary(N,:));
0294 jnk.nodes = mdl.nodes(used_nodes,:);
0295 jnk.boundary = jnk.elems;
0296 if do_debug
0297    show_fem(jnk);
0298    hold on
0299    plot3(elecs.points(:,1),elecs.points(:,2),elecs.points(:,3),'ro');
0300    % plot3(mdl.nodes(nn(outer),1), mdl.nodes(nn(outer),2), mdl.nodes(nn(outer),3),'bs')
0301    hold off
0302 end
0303 
0304 flat = level_model_slice(jnk,struct('centre',mdl.face_centre(fc,:),'normal',mdl.normals(fc,:)));
0305 elec_pts = level_model_slice(elecs.points,struct('centre',mdl.face_centre(fc,:),'normal',mdl.normals(fc,:)));
0306 elec_nodes = level_model_slice(elecs.nodes,struct('centre',mdl.face_centre(fc,:),'normal',mdl.normals(fc,:)));
0307 % now, the points are almost z = 0, so we can work in 2D
0308 warning off 'MATLAB:triangulation:PtsNotInTriWarnId'
0309 TR = triangulation(double(flat.elems), flat.nodes(:,1:2));
0310 warning on 'MATLAB:triangulation:PtsNotInTriWarnId'
0311 
0312 edges = TR.freeBoundary;
0313 
0314 % project all nodes of the faces in N onto the plane of the electrode
0315 PN = flat.nodes(:,1:2);
0316 
0317 % for every electrode point, find closest node
0318 neighbour = TR.nearestNeighbor(elec_pts(:,1:2));
0319 D = sqrt(sum((flat.nodes(neighbour,1:2) - elec_pts(:,1:2)).^2,2));
0320 rm = unique(neighbour(D < 2 * elecs.maxh));
0321 
0322 % we can only delete if it's not part of the boundary
0323 b = unique(edges(:));
0324 rm(ismember(rm,b)) = [];
0325 
0326 % remove and remap
0327 PN(rm,:) = [];
0328 used_nodes(rm) = [];
0329 
0330 n = size(flat.nodes,1);
0331 nodelist = 1:n;
0332 nodelist(rm) = [];
0333 map = zeros(n,1);
0334 map(nodelist) = 1:numel(nodelist);
0335 edges = map(edges); 
0336 
0337 points = [PN; elec_pts(:,1:2)];
0338 
0339 % constrained Delaunay triangulation in 2D
0340 f = length(PN) + (1:2);
0341 C = bsxfun(@plus, (0:length(elecs.points)-2)', f);
0342 [wtxt, wid] = lastwarn;
0343 lastwarn('','');
0344 warning off 'MATLAB:DelaunayTri:ConsConsSplitWarnId';
0345 D = DelaunayTri(points,[edges; C]);
0346 els = D.Triangulation(D.inOutStatus,:);
0347 [txt, id] = lastwarn;
0348 if strcmp(id,'MATLAB:DelaunayTri:ConsConsSplitWarnId')
0349     if do_debug
0350         keyboard
0351     else
0352         error(txt); % no point continuing
0353     end
0354 else
0355     lastwarn(wtxt,wid); % restore
0356 end
0357 warning on 'MATLAB:DelaunayTri:ConsConsSplitWarnId';
0358 % project all electrode points on the faces below them, using the normal of
0359 % the central face
0360 Ne = mdl.normals(fc,:);
0361 FN = TR.pointLocation(elec_nodes(:,1:2)); % face num under each electrode point
0362 FC = fcs(FN); % same, in original numbering
0363 for j = 1:length(elecs.nodes)
0364    Pe = elecs.nodes(j,:);
0365    Nf = mdl.normals(fcs(FN(j)),:);
0366    Cf = mdl.face_centre(fcs(FN(j)),:);
0367    Proj(j,:) = Pe + Ne * dot(Cf-Pe,Nf) / dot(Ne,Nf) ;
0368 end
0369 
0370 
0371 % this is just output
0372 EL1 = Proj(1:length(elecs.points),:);
0373 
0374 % remove any nodes inside the electrode
0375 ln = length(used_nodes);
0376 % IN = inpolygon(x(1:ln),y(1:ln),x(ln+1:end),y(ln+1:end));
0377 % nodes(IN) = [];
0378 
0379 add = elecs.maxh;
0380 
0381 nn = mdl.nodes(used_nodes,:);% + add * repmat(IN,1,3) .* repmat(Ne,ln,1);
0382 le = length(elecs.nodes);
0383 ne = Proj + add * repmat(Ne,le,1);
0384 
0385 %this is just output
0386 EL2 = ne(1:length(elecs.points),:);
0387 V = add*Ne;
0388 
0389 % the nodes of the electrode
0390 % IN = [IN; ones(le,1)];
0391 el_c = D.incenters;
0392 el_c(~D.inOutStatus,:) = [];
0393 e_el = inpolygon(el_c(:,1),el_c(:,2),points(ln+1:end,1),points(ln+1:end,2));
0394 els(e_el,:) = []; % els(e_el,:) + (els(e_el,:)>ln ) .* le;
0395 
0396 % add connecting elements
0397 E = [];
0398 le = length(elecs.points);
0399 f = ln + [ 1 le+1 le+2; le+2 2 1];
0400 for j = 0:(le-2)
0401    E = [E; j+f];
0402 end
0403 M = ln + [le+1 le 2*le; le le+1 1];
0404 E = [E; M];
0405 
0406 jnk.nodes = [nn ; Proj(1:le,:);  ne];
0407 jnk.elems = [ els; E; elecs.elems+ln+le];
0408 jnk.boundary = jnk.elems;
0409 if do_debug
0410    show_fem(jnk);
0411 end
0412 
0413 % remove the patch we're replacing
0414 big = mdl;
0415 big.boundary(N,:) = [];
0416 big.faces(N,:) = [];
0417 big.normals(N,:) = [];
0418 big.face_centre(N,:) = [];
0419 
0420 big.elems = big.boundary;
0421 
0422 joint = merge_meshes(big,jnk,0.001);
0423 joint.boundary = joint.elems;
0424 joint.faces = joint.boundary;
0425 opt.normals = true;
0426 opt.face_centre = true;
0427 joint = fix_model(joint,opt);
0428 
0429 
0430 function PN = project_nodes_on_elec(mdl,elecs,nodes)
0431 fc = find_face_under_elec(mdl,elecs.pos);
0432 Ne = mdl.normals(fc,:);
0433 Pe = elecs.pos;
0434 % for i = 1:length(nodes)
0435 %    P = mdl.nodes(nodes(i),:);
0436 %    PN(i,:) = P + dot(Pe - P, Ne) * Ne;
0437 % end
0438 N = mdl.nodes(nodes,:);
0439 % PN = N + sum((Pe-N) .* Ne,2) .* Ne;
0440 PN = N + bsxfun(@times,sum(bsxfun(@times,bsxfun(@minus,Pe,N), Ne),2), Ne);
0441 
0442 
0443 
0444 % OUTPUT:
0445 %  elecs(i).pos   = [x,y,z]
0446 %  elecs(i).shape = 'C' or 'R'
0447 %  elecs(i).dims  = [radius] or [width,height]
0448 %  elecs(i).maxh  = '-maxh=#' or '';
0449 %  elecs(i).points= list of points around the perimeter
0450 % Angles (th) are interpreted with the mean of boundary nodes as origin
0451 function [elecs] = parse_elecs(mdl, elec_pos, elec_shape )
0452 elecs = [];
0453 
0454 if size(elec_shape,2) < 3
0455    elec_shape(:,3) = elec_shape(:,1)/10;
0456 end
0457 
0458 have_xyz = 0;
0459 
0460 if size(elec_pos,1) == 1
0461    % Parse elec_pos = [n_elecs_per_plane,(0=equal angles,1=equal dist),z_planes]
0462    n_elecs= elec_pos(1); % per plane
0463    offset = elec_pos(2) - floor(elec_pos(2));
0464    switch floor(elec_pos(2))
0465       case 0
0466          th = linspace(0,2*pi, n_elecs+1)'; th(end)=[];
0467          th = th + offset*2*pi;
0468          ind = th >= 2*pi;
0469          th(ind) = th(ind) - 2*pi;
0470       case 1
0471          error('not implemented yet');
0472    end
0473    on_elecs = ones(n_elecs, 1);
0474    el_th = [];
0475    el_z  = [];
0476    for i=3:length(elec_pos)
0477       el_th = [el_th; th];
0478       el_z  = [el_z ; on_elecs*elec_pos(i)];
0479    end
0480 elseif size(elec_pos,2) == 2
0481    % elec_pos = [theta z];
0482    el_th = elec_pos(:,1)*2*pi/360;
0483    el_z  = elec_pos(:,2);
0484 elseif size(elec_pos,2) == 3
0485    % elec_pos = [x y z];
0486    have_xyz = 1;
0487    el_z  = elec_pos(:,3);
0488 end
0489 
0490 if ~have_xyz
0491    el_th(el_th>pi) =  el_th(el_th>pi) - 2*pi;
0492    el_th(el_th<-pi) = el_th(el_th<-pi) + 2*pi;
0493 end
0494 n_elecs= size(el_z,1);
0495 
0496 if size(elec_shape,1) == 1
0497    elec_shape = ones(n_elecs,1) * elec_shape;
0498 end
0499 
0500 for i = 1:n_elecs
0501    if ~have_xyz
0502       [fc elecs(i).pos] = find_elec_centre(mdl,el_th(i),el_z(i));
0503    else
0504       elecs(i).pos = elec_pos(i,:);
0505    end
0506 %    elecs(i).face = fc; % this changes too often to store!
0507    elecs(i).dims = elec_shape(i,1:2);
0508    elecs(i).dims(elecs(i).dims==0) = [];
0509    elecs(i).maxh = elec_shape(i,3);
0510    
0511    if elec_shape(i,2) == 0
0512       elecs(i).shape = 'C';
0513       r = elec_shape(i,1);
0514       n = ceil(2*pi*elec_shape(i,1) / elec_shape(i,3));
0515       t = linspace(0,2*pi,n+1); t(end) = [];
0516       x = r*sin(t); y = r*cos(t);
0517    else
0518       elecs(i).shape = 'R';
0519       height = elec_shape(i,1); width = elec_shape(i,2); d_org = elec_shape(i,3);
0520       % enforce a minimum of 5 nodes per side
0521       d = min( [ d_org , height/5, width/5]);
0522       if d < d_org
0523          elecs(i).maxh = d;
0524          eidors_msg('@@@ Decreased maxh of electrode %d from %f to %f',i,d_org, d,2);
0525       end
0526       nh = ceil(height/d)+1; nw = ceil(width/d)+1; 
0527       ph = linspace(-height/2,height/2,nh);
0528       pw = linspace(-width/2,width/2,nw);
0529       y = [ph, ph(end)*ones(1,nw-2), fliplr(ph), ph(1)*ones(1,nw-2)];
0530       x = [pw(1)*ones(1,nh-1), pw, pw(end)*ones(1,nh-2), fliplr(pw(2:end))];
0531       %    % we don't want real rectangles, because Netgen will merge coplanar
0532       %    % faces, so we create a nice superellipse instead
0533       n = 2*(nh+nw);
0534       %    t = linspace(2*pi,0,n); t(end) = [];
0535       %    N = 8;
0536       %    x = abs(cos(t)).^(2/N) * width/2  .* sign(cos(t));
0537       %    y = abs(sin(t)).^(2/N) * height/2 .* sign(sin(t));
0538       % superellipses are also bad, what about a wavy rectange?
0539 %       [pp] = fourier_fit([x; y]', min(size(x,2),18) );
0540 %       t = linspace(0,1,n+1); t(end) = [];
0541 %       xy = fourier_fit(pp,t);
0542 %       x = xy(:,1)'; y = xy(:,2)';
0543       % wavy rectangles are nice but don't guarantee absence of co-planar
0544       % faces
0545       % let's try a brute-force approach
0546       e = tand(0.5)*d;
0547       x = x + e* [0 power(-1,0:nh-3) zeros(1,nw)  power(-1,0:nh-3) zeros(1,nw-1)];
0548       y = y + e* [zeros(1,nh) power(-1,0:nw-3) zeros(1,nh) power(-1,0:nw-3)];
0549    end
0550    fc = find_face_under_elec(mdl,elecs(i).pos);
0551    [u v s] = get_face_basis(mdl, fc);
0552    
0553    np = length(x);
0554 %    elecs(i).points = flipud(ones(size(x))' * elecs(i).pos + x'*s + y'*v);
0555 
0556    ng_write_opt('meshoptions.fineness',1,'options.meshsize',1.2*elecs(i).maxh);
0557    emdl = ng_mk_2d_model(flipud([x', y']));
0558    x = emdl.nodes(:,1); y = emdl.nodes(:,2);
0559    elecs(i).nodes = ones(size(x)) * elecs(i).pos + x*s + y*v;
0560    elecs(i).elems = emdl.elems(:,[1 3 2]); % flip orientation to the outside
0561    elecs(i).points = elecs(i).nodes(1:np,:); % this must be the boundary
0562    % TODO: write code to check if this is true
0563    
0564 end
0565 delete('ng.opt');
0566 
0567 function [u v s] = get_face_basis(mdl, fc)
0568    u = mdl.normals(fc,:); % unit normal
0569 
0570    % project each coordinate axis on the plane
0571    I = eye(3);
0572    for i = 1:3
0573        proj(:,i) = I(:,i) - (dot(I(:,i),u')) * u';
0574    end
0575    norm_proj = vecnorm(proj);
0576    min_norm = min(norm_proj);
0577    
0578    % vertical vector on the plane of that surface triangle
0579    if norm_proj(3) ~= min_norm
0580       v = [0 0 1] - dot([0 0 1],u) *u;
0581    else
0582       % the element is essentially horizontal
0583 %       v = [0 1 0] - dot([0 1 0],u)*u;
0584 %TODO: need to expose an option to decide which it should be
0585       v = [1 0 0] - dot([1 0 0],u)*u;
0586    end
0587    v = v/norm(v);
0588    s = cross(u,v); s= s/norm(s);
0589 
0590 function [fc pos] = find_elec_centre(mdl, el_th,el_z)
0591 fc = [];
0592 pos = [];
0593 
0594 Ctr = mean(mdl.nodes(mdl.boundary,:));
0595 Ctr(3) = el_z;
0596 
0597 %1. Find edges that cross the z plane
0598 n_above = mdl.nodes(:,3) >= el_z;
0599 sum_above = sum(n_above(mdl.edges),2) ;
0600 edg = sum_above == 1;
0601 
0602 %2. Find an edge that crosses el_th
0603 n = unique(mdl.edges(edg,:));
0604 nn = mdl.nodes(n,1:2);
0605 nn = nn - repmat(Ctr(:,1:2),length(nn),1);
0606 th = cart2pol(nn(:,1),nn(:,2));
0607 th(:,2) = 1:length(th);
0608 th = sortrows(th);
0609 idx = find(th(:,1) > el_th,1,'first');
0610 if isempty(idx) || idx == 1
0611    n1 = n(th(1,2));
0612    n2 = n(th(end,2));
0613    % edges in edg that contain these nodes (they don't need to be on the
0614    % same element)
0615    ed = edg & sum( (mdl.edges == n1) + (mdl.edges == n2) ,2) > 0;
0616 else
0617 %    to_the_left = false(length(mdl.nodes),1);
0618 %    to_the_left(n(th(1:idx-1,2))) = true;
0619 %    sum_left = sum( to_the_left(mdl.boundary), 2);
0620 %    el = els & sum_left > 0 & sum_left < 3;
0621    n1 = n(th(idx-1,2));
0622    n2 = n(th(idx,  2));
0623    ed = edg & sum( (mdl.edges == n1) + (mdl.edges == n2) ,2) > 0;
0624 end
0625 
0626 el = false(length(mdl.boundary),1);
0627 for i = find(ed)'
0628    n1 = mdl.edges(i,1);
0629    n2 = mdl.edges(i,2);
0630    el = el | sum( (mdl.boundary == n1) + (mdl.boundary == n2), 2) == 2;
0631 end
0632 el = find(el);
0633 % fcs = find(mdl.boundary_face);
0634 % fcs = fcs(el);
0635 
0636 [De(1) De(2) De(3)]  = pol2cart(el_th,1, 0); 
0637 for i = 1:length(el)
0638    Nf = mdl.normals(el(i),:);
0639    Cf = mdl.face_centre(el(i),:);
0640    % the plane is (X - Cf).Nf = 0
0641    % the line is X = Ctr + tDe (through Ctr along De
0642    % We want X that satisfies both.
0643    % (Ctr +tDe -  Cf).Nf = 0
0644    % (Ctr - Cf).Nf + tDe.Nf = 0
0645    % t =
0646    % X = Ctr + De * (Cf-Ctr).Nf / (De.Nf)
0647    t = dot(Cf-Ctr,Nf) / dot(De,Nf);
0648    if t < 0, continue, end
0649    X = Ctr + De * t ;
0650    if point_in_triangle(X, mdl.faces(el(i),:), mdl.nodes)
0651       pos = X;
0652       fc = el(i);
0653       break;
0654    end
0655    
0656    % project the line on this element
0657    % check if it falls inside
0658 end
0659 if isempty(pos)
0660    keyboard
0661 end
0662 
0663 function [out, fc] = grow_neighbourhood(mdl, varargin)
0664 use_elec = false;
0665 if length(varargin) == 1
0666    use_elec = true;
0667    elecs = varargin{1};
0668    fc = find_face_under_elec(mdl,elecs.pos);
0669    p = elecs.pos;
0670    switch elecs.shape
0671       case 'R'
0672          r = sqrt(sum(elecs.dims.^2,2));
0673       case 'C'
0674          r = 2 * elecs.dims(1);
0675    end
0676 else
0677    fc = varargin{1};
0678    p = varargin{2};
0679    r = varargin{3};
0680 end
0681 
0682 done = false(length(mdl.boundary),1);
0683 todo = false(length(mdl.boundary),1);
0684 todo(fc) = true;
0685 bb = mdl.boundary;
0686 vv = mdl.nodes;
0687 % distance of each vertex to the line perpendicular to face fc passing
0688 % through p
0689 dv = vv - repmat(p,length(vv),1);
0690 nl = mdl.normals;
0691 nl = repmat(nl(fc,:),length(vv),1);
0692 dd = sqrt(sum( (dv - repmat(dot(dv,nl,2),1,3) .* nl).^2,2));
0693 dim = size(bb,2);
0694 first = true; % at first iteration, add all neighbours
0695 if use_elec
0696    PN = level_model_slice(mdl.nodes,struct('centre',mdl.face_centre(fc,:),'normal',mdl.normals(fc,:)));
0697    elec_pts = level_model_slice(elecs.points,struct('centre',mdl.face_centre(fc,:),'normal',mdl.normals(fc,:)));
0698    elec_pts = elec_pts(:,1:2);
0699    PN = PN(:,1:2);
0700    emin = min(elec_pts);
0701    emax = max(elec_pts);
0702    rng = emax-emin;
0703    emin = emin - 0.1*rng;
0704    emax = emax + 0.1*rng;
0705    toofar = false(size(mdl.boundary,1),1);
0706    
0707    for i = 1:2
0708       nodes = reshape(PN(mdl.boundary,i),[],3);
0709       toofar =  toofar |  sum(nodes > emax(i),2) == 3 | sum(nodes < emin(i),2) == 3;
0710    end
0711 end
0712 while any(todo)
0713    id = find(todo,1,'first');
0714    done(id) = 1;
0715    nn = find_neighbours(id,bb);
0716    if use_elec
0717       nn = nn & ~toofar;
0718    elseif first
0719       % include all neighbours
0720       first = false;
0721    else
0722       % at least one node must be close enough
0723       nn = nn & sum(dd(bb) <= r,2) > 0;
0724    end
0725    todo = todo | nn;
0726    todo(done) = 0;
0727 %    disp(sprintf('id: %d done: %d todo: %d',id, nnz(done),nnz(todo)));
0728 %    disp(find(todo)');
0729 %    disp(find(done)');
0730 end
0731 out = find(done);
0732 
0733 
0734 function nn =  find_neighbours(fc, bb);
0735 dim = size(bb,2);
0736 nn = false(length(bb),1);
0737 for i = 1:dim
0738    node = bb(fc,i);
0739    nn = nn | sum(bb == node,2) > 0;
0740 end
0741 nn(fc) = 0;
0742 
0743 function [e p] = find_face_under_elec(mdl, elec_pos)
0744 for i = 1:size(elec_pos,1)
0745    % 1. Project electrode on all faces
0746    ee = repmat(elec_pos(i,:),length(mdl.faces),1);
0747    fc = mdl.face_centre;
0748    n  = mdl.normals;
0749    proj1 = ee - repmat(dot(ee-fc, n,2),1,3) .* n;
0750    in1 = point_in_triangle(proj1,mdl.faces,mdl.nodes, 'match');
0751    dis1 = sqrt(sum((ee-proj1).^2,2));
0752    % 2. Project electrode on all edges
0753    edg = [mdl.faces(:,1:2);mdl.faces(:,2:3);mdl.faces(:,[3 1])];
0754    edg = sort(edg,2);
0755    [edg jnk e2f] = unique(edg,'rows');
0756    ee = repmat(elec_pos(i,:),length(edg),1);
0757    s = mdl.nodes(edg(:,2),:) - mdl.nodes(edg(:,1),:); %edge direction vector
0758    t = dot(ee-mdl.nodes(edg(:,1),:),s,2)./dot(s,s,2);
0759    in2 = t>=0 & t <=1;
0760    in2 = any(reshape(in2(e2f),[],3),2);
0761    proj2 = mdl.nodes(edg(:,1),:) + repmat(t,1,3).*s;
0762    dis = sqrt(sum((ee - proj2).^2,2));
0763    dis = repmat(dis,2,1);
0764    dis(t<0 | t > 1) = Inf;
0765    dis = reshape(dis(e2f),[],3);
0766    [jnk, pos] = min(dis,[],2);
0767    idx = sparse(1:length(pos),pos,1);
0768    dis = dis';
0769    dis2 = dis(logical(idx'));
0770 
0771    in = in1 | in2;
0772    if nnz(in) == 1
0773          e(i) = find(in1);  % this should be an index into mdl.boundary
0774          p(i,:) = proj1(in1,:);
0775    else
0776       % take the element that is closest to ee
0777       cand = find(in);
0778       dd(in1(cand)) = dis1(in1);
0779       dd(in2(cand)) = dis2(in2);
0780       [jnk pos] = min(dd);
0781       e(i) = cand(pos);
0782       p(i,:) = proj1(e(i),:);
0783    end
0784 
0785 end
0786 
0787 
0788 function do_unit_test
0789 xy= [ -0.89 -0.74 -0.21  0.31  0.79  0.96  0.67  0.05 -0.36 -0.97;
0790        0.14  0.51  0.35  0.50  0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0791 [fmdl] = ng_mk_extruded_model({2,xy,[4,80],},[],[]);
0792 elec_pos = [-0.5, -0.8, 1];
0793 % place_elec_on_surf(fmdl, elec_pos, [0.1 0 0.01]);
0794 % place_elec_on_surf(fmdl, elec_pos, [0.15 0.1 0.01]);
0795 mdl = place_elec_on_surf(fmdl, [16 0 1], [0.15 0.1 0.01]);
0796 % place_elec_on_surf(fmdl, [16 0 1], [0.1 0 0.01]);
0797 subplot(121)
0798 show_fem(mdl);
0799 
0800 mdl = place_elec_on_surf(fmdl, [16 0 1], [0.15 0.1 0.01],[],0.1);
0801 % place_elec_on_surf(fmdl, [16 0 1], [0.1 0 0.01]);
0802 subplot(122)
0803 show_fem(mdl);