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 5207 2016-03-06 21:30:16Z 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 = grow_neighbourhood(mdl,elecs(i));
0097       [mdl E1{i} E2{i} V{i}] = add_electrodes(mdl,N,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,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 
0289 jnk.type = 'fwd_model';
0290 jnk.elems = mdl.boundary(N,:);
0291 jnk.nodes = mdl.nodes;
0292 jnk.boundary = jnk.elems;
0293 img = mk_image(jnk,1);
0294 if do_debug
0295    show_fem(jnk);
0296    hold on
0297    plot3(elecs.points(:,1),elecs.points(:,2),elecs.points(:,3),'ro');
0298    % plot3(mdl.nodes(nn(outer),1), mdl.nodes(nn(outer),2), mdl.nodes(nn(outer),3),'bs')
0299    hold off
0300 end
0301 
0302 
0303 %nodes used
0304 [nn,I, J] = unique(mdl.faces(fcs,:));
0305 outer = true(size(nn));
0306 for i = 1:length(nn)
0307    if sum(J==i) == sum(mdl.boundary(:) == nn(i))
0308       outer(i) = false;
0309    end
0310 end
0311 % we want to keep the ones on the outside
0312 keep = false(size(mdl.nodes,1),1);
0313 keep(nn) = outer;
0314 keep = keep(jnk.elems);
0315 
0316 % this will not catch the situation where the element reaches from boundary
0317 % to boundary and the electrode is in the middle (small electrode, big
0318 % element). Fortunately, in these cases the edge will be there twice
0319 edges = reshape(jnk.elems(:,[1 2 2 3 3 1])',2,[])';
0320 if size(keep,2) == 1; keep = shiftdim(keep,1); end
0321 keep = reshape(keep(:,[1 2 2 3 3 1])',2,[])';
0322 rm = sum(keep,2)<2;
0323 edges(rm,:) = [];
0324 
0325 % detect and remove double entries
0326 rm = ismember(edges,edges(:,[2 1]),'rows');
0327 edges(rm,:) = [];
0328 
0329 % project all nodes of the faces in N onto the plane of the electrode
0330 nodes = unique(mdl.faces(fcs,:));
0331 PN = project_nodes_on_elec(mdl,elecs,nodes);
0332 
0333 % electrode coordinate system
0334 [u v s] = get_face_basis(mdl,fc);
0335 % u = mdl.normals(fc,:); % unit normal
0336 % % vertical vector on the plane of that surface triangle
0337 % v = [0 0 1] - dot([0 0 1],u) *u; v = v/norm(v);
0338 % s = cross(u,v); s= s/norm(s);
0339 
0340 % mark nodes that are too close to elecs.points for removal
0341 rm = false(length(PN),1);
0342 for i = 1:length(PN)
0343    D = repmat(PN(i,:),length(elecs.points),1) - elecs.points;
0344    D = sqrt(sum(D.^2,2));
0345    if any(D < 2*elecs.maxh)
0346       rm(i) = true;
0347    end
0348 end
0349 
0350 % we can only delete if it's not part of the boundary
0351 b = unique(edges(:));
0352 rm = find(rm);
0353 rm(ismember(nodes(rm),b)) = [];
0354 
0355 % remove and remap
0356 PN(rm,:) = [];
0357 nodes(rm) = [];
0358 
0359 points = [PN; elecs.points];
0360 np = size(points,1);
0361 x = dot(points,repmat(v,np,1),2);
0362 y = dot(points,repmat(s,np,1),2);
0363 
0364 map(nodes) = 1:length(nodes);
0365 edges = map(edges); %
0366 
0367 % constrained Delaunay triangulation in 2D
0368 f = length(PN) +(1:2);
0369 C = [];
0370 for i= 0:length(elecs.points)-2
0371    C = [C; i+f];
0372 end
0373 D = DelaunayTri([x y],[edges; C]);
0374 els = D.Triangulation(D.inOutStatus,:);
0375 
0376 
0377 % project all electrode points on all triangles, using the normal of the central elem
0378 Ne = mdl.normals(fc,:);
0379 for j = 1:length(elecs.nodes)
0380    Pe = elecs.nodes(j,:);
0381    for i = 1:length(fcs)
0382       Nf = mdl.normals(fcs(i),:);
0383       Cf = mdl.face_centre(fcs(i),:);
0384       % the plane is (X - Cf).Nf = 0
0385       % the line is X = Pe + tNe (through Pe perpendicular to the main elec
0386       % face
0387       % We want X that satisfies both.
0388       % (Pe +tNe -  Cf).Nf = 0
0389       % (Pe - Cf).Nf + tNe.Nf = 0
0390       % t = (Cf-Pe).Nf / (Ne.Nf)
0391       % X = Pe + Ne * (Cf-Pe).Nf / (Ne.Nf)
0392       X = Pe + Ne * dot(Cf-Pe,Nf) / dot(Ne,Nf) ;
0393       if point_in_triangle(X, mdl.faces(fcs(i),:), mdl.nodes)
0394          Proj(j,:) = X;
0395          FC(j) = fcs(i);
0396          break;
0397       end
0398    end
0399 end
0400 
0401 % this is just output
0402 EL1 = Proj(1:length(elecs.points),:);
0403 
0404 % remove any nodes inside the electrode
0405 ln = length(nodes);
0406 % IN = inpolygon(x(1:ln),y(1:ln),x(ln+1:end),y(ln+1:end));
0407 % nodes(IN) = [];
0408 
0409 add = elecs.maxh;
0410 
0411 nn = mdl.nodes(nodes,:);% + add * repmat(IN,1,3) .* repmat(Ne,ln,1);
0412 le = length(elecs.nodes);
0413 ne = Proj + add * repmat(Ne,le,1);
0414 
0415 %this is just output
0416 EL2 = ne(1:length(elecs.points),:);
0417 V = add*Ne;
0418 
0419 % the nodes of the electrode
0420 % IN = [IN; ones(le,1)];
0421 el_c = D.incenters;
0422 el_c(~D.inOutStatus,:) = [];
0423 e_el = inpolygon(el_c(:,1),el_c(:,2),x(ln+1:end),y(ln+1:end));
0424 els(e_el,:) = []; % els(e_el,:) + (els(e_el,:)>ln ) .* le;
0425 
0426 % add connecting elements
0427 E = [];
0428 le = length(elecs.points);
0429 f = ln + [ 1 le+1 le+2; le+2 2 1];
0430 for j = 0:(le-2)
0431    E = [E; j+f];
0432 end
0433 M = ln + [le+1 le 2*le; le le+1 1];
0434 E = [E; M];
0435 
0436 jnk.nodes = [nn ; Proj(1:le,:);  ne];
0437 jnk.elems = [ els; E; elecs.elems+ln+le];
0438 jnk.boundary = jnk.elems;
0439 if do_debug
0440    show_fem(jnk);
0441 end
0442 
0443 % remove the patch we're replacing
0444 big = mdl;
0445 big.boundary(N,:) = [];
0446 big.faces(N,:) = [];
0447 big.normals(N,:) = [];
0448 big.face_centre(N,:) = [];
0449 
0450 big.elems = big.boundary;
0451 
0452 joint = merge_meshes(big,jnk,0.001);
0453 joint.boundary = joint.elems;
0454 joint.faces = joint.boundary;
0455 opt.normals = true;
0456 opt.face_centre = true;
0457 joint = fix_model(joint,opt);
0458 
0459 function PN = project_nodes_on_elec(mdl,elecs,nodes)
0460 fc = find_face_under_elec(mdl,elecs.pos);
0461 Ne = mdl.normals(fc,:);
0462 Pe = elecs.pos;
0463 for i = 1:length(nodes)
0464    P = mdl.nodes(nodes(i),:);
0465    PN(i,:) = P + dot(Pe - P, Ne) * Ne;
0466 end
0467 
0468 % OUTPUT:
0469 %  elecs(i).pos   = [x,y,z]
0470 %  elecs(i).shape = 'C' or 'R'
0471 %  elecs(i).dims  = [radius] or [width,height]
0472 %  elecs(i).maxh  = '-maxh=#' or '';
0473 %  elecs(i).points= list of points around the perimeter
0474 % Angles (th) are interpreted with the mean of boundary nodes as origin
0475 function [elecs] = parse_elecs(mdl, elec_pos, elec_shape )
0476 elecs = [];
0477 
0478 if size(elec_shape,2) < 3
0479    elec_shape(:,3) = elec_shape(:,1)/10;
0480 end
0481 
0482 have_xyz = 0;
0483 
0484 if size(elec_pos,1) == 1
0485    % Parse elec_pos = [n_elecs_per_plane,(0=equal angles,1=equal dist),z_planes]
0486    n_elecs= elec_pos(1); % per plane
0487    offset = elec_pos(2) - floor(elec_pos(2));
0488    switch floor(elec_pos(2))
0489       case 0
0490          th = linspace(0,2*pi, n_elecs+1)'; th(end)=[];
0491          th = th + offset*2*pi;
0492          ind = th >= 2*pi;
0493          th(ind) = th(ind) - 2*pi;
0494       case 1
0495          error('not implemented yet');
0496    end
0497    on_elecs = ones(n_elecs, 1);
0498    el_th = [];
0499    el_z  = [];
0500    for i=3:length(elec_pos)
0501       el_th = [el_th; th];
0502       el_z  = [el_z ; on_elecs*elec_pos(i)];
0503    end
0504 elseif size(elec_pos,2) == 2
0505    % elec_pos = [theta z];
0506    el_th = elec_pos(:,1)*2*pi/360;
0507    el_z  = elec_pos(:,2);
0508 elseif size(elec_pos,2) == 3
0509    % elec_pos = [x y z];
0510    have_xyz = 1;
0511    el_z  = elec_pos(:,3);
0512 end
0513 
0514 if ~have_xyz
0515    el_th(el_th>pi) =  el_th(el_th>pi) - 2*pi;
0516    el_th(el_th<-pi) = el_th(el_th<-pi) + 2*pi;
0517 end
0518 n_elecs= size(el_z,1);
0519 
0520 if size(elec_shape,1) == 1
0521    elec_shape = ones(n_elecs,1) * elec_shape;
0522 end
0523 
0524 for i = 1:n_elecs
0525    if ~have_xyz
0526       [fc elecs(i).pos] = find_elec_centre(mdl,el_th(i),el_z(i));
0527    else
0528       elecs(i).pos = elec_pos(i,:);
0529    end
0530 %    elecs(i).face = fc; % this changes too often to store!
0531    elecs(i).dims = elec_shape(i,1:2);
0532    elecs(i).dims(elecs(i).dims==0) = [];
0533    elecs(i).maxh = elec_shape(i,3);
0534    
0535    if elec_shape(i,2) == 0
0536       elecs(i).shape = 'C';
0537       r = elec_shape(i,1);
0538       n = ceil(2*pi*elec_shape(i,1) / elec_shape(i,3));
0539       t = linspace(0,2*pi,n+1); t(end) = [];
0540       x = r*sin(t); y = r*cos(t);
0541    else
0542       elecs(i).shape = 'R';
0543       height = elec_shape(i,1); width = elec_shape(i,2); d_org = elec_shape(i,3);
0544       % enforce a minimum of 5 nodes per side
0545       d = min( [ d_org , height/5, width/5]);
0546       if d < d_org
0547          elecs(i).maxh = d;
0548          eidors_msg('@@@ Decreased maxh of electrode %d from %f to %f',i,d_org, d,2);
0549       end
0550       nh = ceil(height/d)+1; nw = ceil(width/d)+1; 
0551       ph = linspace(-height/2,height/2,nh);
0552       pw = linspace(-width/2,width/2,nw);
0553       y = [ph, ph(end)*ones(1,nw-2), fliplr(ph), ph(1)*ones(1,nw-2)];
0554       x = [pw(1)*ones(1,nh-1), pw, pw(end)*ones(1,nh-2), fliplr(pw(2:end))];
0555       %    % we don't want real rectangles, because Netgen will merge coplanar
0556       %    % faces, so we create a nice superellipse instead
0557       n = 2*(nh+nw);
0558       %    t = linspace(2*pi,0,n); t(end) = [];
0559       %    N = 8;
0560       %    x = abs(cos(t)).^(2/N) * width/2  .* sign(cos(t));
0561       %    y = abs(sin(t)).^(2/N) * height/2 .* sign(sin(t));
0562       % superellipses are also bad, what about a wavy rectange?
0563 %       [pp] = fourier_fit([x; y]', min(size(x,2),18) );
0564 %       t = linspace(0,1,n+1); t(end) = [];
0565 %       xy = fourier_fit(pp,t);
0566 %       x = xy(:,1)'; y = xy(:,2)';
0567       % wavy rectangles are nice but don't guarantee absence of co-planar
0568       % faces
0569       % let's try a brute-force approach
0570       e = tand(0.5)*d;
0571       x = x + e* [0 power(-1,0:nh-3) zeros(1,nw)  power(-1,0:nh-3) zeros(1,nw-1)];
0572       y = y + e* [zeros(1,nh) power(-1,0:nw-3) zeros(1,nh) power(-1,0:nw-3)];
0573    end
0574    fc = find_face_under_elec(mdl,elecs(i).pos);
0575    [u v s] = get_face_basis(mdl, fc);
0576    
0577    np = length(x);
0578 %    elecs(i).points = flipud(ones(size(x))' * elecs(i).pos + x'*s + y'*v);
0579 
0580    ng_write_opt('meshoptions.fineness',1,'options.meshsize',1.2*elecs(i).maxh);
0581    emdl = ng_mk_2d_model(flipud([x', y']));
0582    x = emdl.nodes(:,1); y = emdl.nodes(:,2);
0583    elecs(i).nodes = ones(size(x)) * elecs(i).pos + x*s + y*v;
0584    elecs(i).elems = emdl.elems(:,[1 3 2]); % flip orientation to the outside
0585    elecs(i).points = elecs(i).nodes(1:np,:); % this must be the boundary
0586    % TODO: write code to check if this is true
0587    
0588 end
0589 delete('ng.opt');
0590 
0591 function [u v s] = get_face_basis(mdl, fc)
0592    u = mdl.normals(fc,:); % unit normal
0593    % vertical vector on the plane of that surface triangle
0594    v = [0 0 1] - dot([0 0 1],u) *u;
0595    if norm(v) == 0
0596       % the element is horizontal
0597       v = [0 1 0] - dot([0 1 0],u)*u;
0598    end
0599    v = v/norm(v);
0600    s = cross(u,v); s= s/norm(s);
0601 
0602 function [fc pos] = find_elec_centre(mdl, el_th,el_z)
0603 fc = [];
0604 pos = [];
0605 
0606 Ctr = mean(mdl.nodes(mdl.boundary,:));
0607 Ctr(3) = el_z;
0608 
0609 %1. Find edges that cross the z plane
0610 n_above = mdl.nodes(:,3) >= el_z;
0611 sum_above = sum(n_above(mdl.edges),2) ;
0612 edg = sum_above == 1;
0613 
0614 %2. Find an edge that crosses el_th
0615 n = unique(mdl.edges(edg,:));
0616 nn = mdl.nodes(n,1:2);
0617 nn = nn - repmat(Ctr(:,1:2),length(nn),1);
0618 th = cart2pol(nn(:,1),nn(:,2));
0619 th(:,2) = 1:length(th);
0620 th = sortrows(th);
0621 idx = find(th(:,1) > el_th,1,'first');
0622 if isempty(idx) || idx == 1
0623    n1 = n(th(1,2));
0624    n2 = n(th(end,2));
0625    % edges in edg that contain these nodes (they don't need to be on the
0626    % same element)
0627    ed = edg & sum( (mdl.edges == n1) + (mdl.edges == n2) ,2) > 0;
0628 else
0629 %    to_the_left = false(length(mdl.nodes),1);
0630 %    to_the_left(n(th(1:idx-1,2))) = true;
0631 %    sum_left = sum( to_the_left(mdl.boundary), 2);
0632 %    el = els & sum_left > 0 & sum_left < 3;
0633    n1 = n(th(idx-1,2));
0634    n2 = n(th(idx,  2));
0635    ed = edg & sum( (mdl.edges == n1) + (mdl.edges == n2) ,2) > 0;
0636 end
0637 
0638 el = false(length(mdl.boundary),1);
0639 for i = find(ed)'
0640    n1 = mdl.edges(i,1);
0641    n2 = mdl.edges(i,2);
0642    el = el | sum( (mdl.boundary == n1) + (mdl.boundary == n2), 2) == 2;
0643 end
0644 el = find(el);
0645 % fcs = find(mdl.boundary_face);
0646 % fcs = fcs(el);
0647 
0648 [De(1) De(2) De(3)]  = pol2cart(el_th,1, 0); 
0649 for i = 1:length(el)
0650    Nf = mdl.normals(el(i),:);
0651    Cf = mdl.face_centre(el(i),:);
0652    % the plane is (X - Cf).Nf = 0
0653    % the line is X = Ctr + tDe (through Ctr along De
0654    % We want X that satisfies both.
0655    % (Ctr +tDe -  Cf).Nf = 0
0656    % (Ctr - Cf).Nf + tDe.Nf = 0
0657    % t =
0658    % X = Ctr + De * (Cf-Ctr).Nf / (De.Nf)
0659    t = dot(Cf-Ctr,Nf) / dot(De,Nf);
0660    if t < 0, continue, end
0661    X = Ctr + De * t ;
0662    if point_in_triangle(X, mdl.faces(el(i),:), mdl.nodes)
0663       pos = X;
0664       fc = el(i);
0665       break;
0666    end
0667    
0668    % project the line on this element
0669    % check if it falls inside
0670 end
0671 if isempty(pos)
0672    keyboard
0673 end
0674 
0675 function out = grow_neighbourhood(mdl, varargin)
0676 use_elec = false;
0677 if length(varargin) == 1
0678    use_elec = true;
0679    elecs = varargin{1};
0680    fc = find_face_under_elec(mdl,elecs.pos);
0681    p = elecs.pos;
0682    switch elecs.shape
0683       case 'R'
0684          r = sqrt(sum(elecs.dims.^2,2));
0685       case 'C'
0686          r = 2 * elecs.dims(1);
0687    end
0688 else
0689    fc = varargin{1};
0690    p = varargin{2};
0691    r = varargin{3};
0692 end
0693 
0694 done = false(length(mdl.boundary),1);
0695 todo = false(length(mdl.boundary),1);
0696 todo(fc) = true;
0697 bb = mdl.boundary;
0698 vv = mdl.nodes;
0699 % distance of each vertex to the line perpendicular to face fc passing
0700 % through p
0701 dv = vv - repmat(p,length(vv),1);
0702 nl = mdl.normals;
0703 nl = repmat(nl(fc,:),length(vv),1);
0704 dd = sqrt(sum( (dv - repmat(dot(dv,nl,2),1,3) .* nl).^2,2));
0705 dim = size(bb,2);
0706 first = true; % at first iteration, add all neighbours
0707 if use_elec
0708    PN = project_nodes_on_elec(mdl,elecs,1:length(mdl.nodes));
0709    emin = min(elecs.points);
0710    emax = max(elecs.points);
0711    rng = emax-emin;
0712    emin = emin - 0.1*rng;
0713    emax = emax + 0.1*rng;
0714    toofar = false(size(mdl.boundary,1),1);
0715    
0716    for i = 1:3
0717       nodes = reshape(PN(mdl.boundary,i),[],3);
0718       toofar =  toofar |  sum(nodes > emax(i),2) == 3 | sum(nodes < emin(i),2) == 3;
0719    end
0720 end
0721 while any(todo)
0722    id = find(todo,1,'first');
0723    done(id) = 1;
0724    nn = find_neighbours(id,bb);
0725    if use_elec
0726       nn = nn & ~toofar;
0727    elseif first
0728       % include all neighbours
0729       first = false;
0730    else
0731       % at least one node must be close enough
0732       nn = nn & sum(dd(bb) <= r,2) > 0;
0733    end
0734    todo = todo | nn;
0735    todo(done) = 0;
0736 %    disp(sprintf('id: %d done: %d todo: %d',id, nnz(done),nnz(todo)));
0737 %    disp(find(todo)');
0738 %    disp(find(done)');
0739 end
0740 out = find(done);
0741 
0742 
0743 function nn =  find_neighbours(fc, bb);
0744 dim = size(bb,2);
0745 nn = false(length(bb),1);
0746 for i = 1:dim
0747    node = bb(fc,i);
0748    nn = nn | sum(bb == node,2) > 0;
0749 end
0750 nn(fc) = 0;
0751 
0752 function [e p] = find_face_under_elec(mdl, elec_pos)
0753 for i = 1:size(elec_pos,1)
0754    % 1. Project electrode on all faces
0755    ee = repmat(elec_pos(i,:),length(mdl.faces),1);
0756    fc = mdl.face_centre;
0757    n  = mdl.normals;
0758    proj1 = ee - repmat(dot(ee-fc, n,2),1,3) .* n;
0759    in1 = point_in_triangle(proj1,mdl.faces,mdl.nodes, 'match');
0760    dis1 = sqrt(sum((ee-proj1).^2,2));
0761    % 2. Project electrode on all edges
0762    edg = [mdl.faces(:,1:2);mdl.faces(:,2:3);mdl.faces(:,[3 1])];
0763    edg = sort(edg,2);
0764    [edg jnk e2f] = unique(edg,'rows');
0765    ee = repmat(elec_pos(i,:),length(edg),1);
0766    s = mdl.nodes(edg(:,2),:) - mdl.nodes(edg(:,1),:); %edge direction vector
0767    t = dot(ee-mdl.nodes(edg(:,1),:),s,2)./dot(s,s,2);
0768    in2 = t>=0 & t <=1;
0769    in2 = any(reshape(in2(e2f),[],3),2);
0770    proj2 = mdl.nodes(edg(:,1),:) + repmat(t,1,3).*s;
0771    dis = sqrt(sum((ee - proj2).^2,2));
0772    dis = repmat(dis,2,1);
0773    dis(t<0 | t > 1) = Inf;
0774    dis = reshape(dis(e2f),[],3);
0775    [jnk, pos] = min(dis,[],2);
0776    idx = sparse(1:length(pos),pos,1);
0777    dis = dis';
0778    dis2 = dis(logical(idx'));
0779 
0780    in = in1 | in2;
0781    if nnz(in) == 1
0782          e(i) = find(in1);  % this should be an index into mdl.boundary
0783          p(i,:) = proj1(in1,:);
0784    else
0785       % take the element that is closest to ee
0786       cand = find(in);
0787       dd(in1(cand)) = dis1(in1);
0788       dd(in2(cand)) = dis2(in2);
0789       [jnk pos] = min(dd);
0790       e(i) = cand(pos);
0791       p(i,:) = proj1(e(i),:);
0792    end
0793 
0794 end
0795 
0796 
0797 function do_unit_test
0798 xy= [ -0.89 -0.74 -0.21  0.31  0.79  0.96  0.67  0.05 -0.36 -0.97;
0799        0.14  0.51  0.35  0.50  0.27 -0.23 -0.86 -0.69 -0.85 -0.46]';
0800 [fmdl] = ng_mk_extruded_model({2,xy,[4,80],},[],[]);
0801 elec_pos = [-0.5, -0.8, 1];
0802 % place_elec_on_surf(fmdl, elec_pos, [0.1 0 0.01]);
0803 % place_elec_on_surf(fmdl, elec_pos, [0.15 0.1 0.01]);
0804 mdl = place_elec_on_surf(fmdl, [16 0 1], [0.15 0.1 0.01]);
0805 % place_elec_on_surf(fmdl, [16 0 1], [0.1 0 0.01]);
0806 subplot(121)
0807 show_fem(mdl);
0808 
0809 mdl = place_elec_on_surf(fmdl, [16 0 1], [0.15 0.1 0.01],[],0.1);
0810 % place_elec_on_surf(fmdl, [16 0 1], [0.1 0 0.01]);
0811 subplot(122)
0812 show_fem(mdl);