place_elec_on_surf

PURPOSE ^

PLACE_ELEC_ON_SURF Place electrodes on the surface of a model

SYNOPSIS ^

function mdl2 = place_elec_on_surf(mdl,elec_pos, elec_spec,ng_opt_file, maxh)

DESCRIPTION ^

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

Generated on Tue 12-May-2015 16:00:42 by m2html © 2005