0001 function out = mk_head_model_adult(varargin)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063 citeme(mfilename)
0064
0065 if nargin==1 && ischar(varargin{1}) && strcmp(varargin{1},'UNIT_TEST')
0066 do_unit_test; return;
0067 end
0068
0069 if nargin == 0
0070 out = get_base_model;
0071 end
0072
0073 copyright = '';
0074 organspec = 'coarse';
0075 if nargin > 1 && ischar(varargin{end})
0076 organspec = varargin{end};
0077 varargin(end) = [];
0078 end
0079
0080 switch length(varargin)
0081 case 0
0082
0083 case 1
0084 if ischar(varargin{1})
0085 if strcmp(varargin{1},'show_10-10')
0086 [out, img, ~] = get_base_model;
0087 get_elec_pos(img, true);
0088 return
0089 end
0090 copyright = '(C) EIDORS Project';
0091 out = build_predef_elecs(varargin{1}, organspec);
0092 else
0093 error('EIDORS:WrongInput','Electrode specification not understood');
0094 end
0095 case {2, 3}
0096 [out, img, materials] = get_base_model;
0097 out = build_head_elecs(out, varargin{:});
0098 out = mk_shell_image(img, materials, out, organspec);
0099
0100 otherwise
0101 error('EIDORS:WrongInput', 'Too many inputs');
0102 end
0103
0104 out = eidors_obj('set', out);
0105
0106 if ~isempty(copyright)
0107 out.copyright = copyright;
0108 end
0109
0110 out.attribution = ...
0111 ['Modified from "SAH262" by Andrew Tizzard '...
0112 '(Tizzard, A. and Bayford, R.H. (2007). Improving the Finite Element ' ...
0113 'Forward Model of the Human Head by Warping using Elastic Deformation. '...
0114 'Physiol.Meas. 28:S163-S182)'];
0115
0116 if strcmp(out.type, 'image')
0117 try out.fwd_model.copyright = out.copyright; end
0118 out.fwd_model.attribution = out.attribution;
0119 end
0120
0121 ws = warning('off', 'backtrace');
0122 warning('EIDORS:License', ...
0123 ['The head model is a derivative of worked licensed under CC BY.\n' ...
0124 'Please make sure to attribute correctly. ' ...
0125 'See the attribution field for details.']);
0126 warning(ws.state, 'backtrace');
0127
0128 function out = build_predef_elecs(stimpat, organspec)
0129 names = [];
0130 switch stimpat
0131 case {'10-10', 'no_elecs'}
0132 mdl_name = stimpat;
0133 case {'2x16_planar'; '2x16_odd-even'; '2x16_square';
0134 '2x16_adjacent';'1x32_ring'}
0135 mdl_name = '3_rings';
0136 otherwise
0137 error('EIDORS:WrongInput', 'Don''t understand model %s.', stimpat )
0138 end
0139
0140 cache_path = get_cache_path;
0141 fname = sprintf('head_model_%s',mdl_name);
0142 suffix = '';
0143 if ~isempty(organspec)
0144 fname = [fname '_' organspec];
0145 suffix = [' - ' organspec];
0146 end
0147 name = ['Adult Head Model' suffix ' - ' stimpat];
0148
0149 fpath = [cache_path filesep fname '.mat'];
0150 if exist(fpath,'file')
0151 eidors_msg('@@ Using stored model.', 3);
0152 load(fpath, 'out');
0153 else
0154 [out, img, materials] = get_base_model;
0155 switch stimpat
0156 case 'no_elecs'
0157 elec_pos = [];
0158 case '10-10'
0159 [elec_pos, names] = get_elec_pos(img);
0160 elec_spec = [2 0 0.2];
0161 maxh = 8;
0162 case {'2x16_planar'; '2x16_odd-even'; '2x16_square';
0163 '2x16_adjacent';'1x32_ring'}
0164 Nel = 32;
0165 eth32 = -pi/2 : 2*pi/Nel : 2*pi-pi/2; eth32(1)=[];
0166 eth16 = eth32(1:2:end);
0167 eth32 = eth32 - pi/Nel;
0168 t = [eth16, eth32, eth16];
0169 xm=0;
0170 ym=12;
0171 a=80; b=95;
0172 elec_pos(:,1) = xm+a*cos(t);
0173 elec_pos(:,2) = ym+b*sin(t);
0174 elec_pos(1:16,3) = 0;
0175 elec_pos(17:48,3) = 5;
0176 elec_pos(49:64,3) = 10;
0177 elec_spec = [2 0 0.2];
0178 maxh = 8;
0179 otherwise
0180 error('EIDORS:WrongInput','Parameter electsr not understood');
0181 end
0182 if ~isempty(elec_pos)
0183 out = build_head_elecs(out, elec_pos, elec_spec, maxh);
0184 end
0185 out.name = name;
0186 if ~isempty(names)
0187 for i = 1:length(out.electrode)
0188 out.electrode(i).label = names{i};
0189 end
0190 end
0191 out = mk_shell_image(img, materials, out, organspec);
0192 save(fpath,'out');
0193
0194 end
0195
0196 out.name = name;
0197
0198 if strcmp(mdl_name,'3_rings')
0199 out = set_predef_stim(out, stimpat);
0200 end
0201
0202 function out = mk_shell_image(img, materials, out, organspec)
0203 if isempty(organspec), return, end
0204 switch organspec
0205 case 'void'
0206 return
0207 case 'coarse'
0208 maxh = 5;
0209 case 'fine'
0210 maxh = 2;
0211 otherwise
0212 error('EIDORS:WrongInput', 'Parameter tissuespec not understood.');
0213 end
0214 [CSF, brain, skull2, bn2] = build_brain_CSF(img,materials,maxh);
0215
0216
0217 cond = @(xyz) xyz(:,3) < -121.1 & ...
0218 xyz(:,3) < (abs(xyz(:,2).^1/11) -123);
0219 [skull1, ~, bn1] = split_surface(out, cond);
0220 skull1.elems(:,2:3) = skull1.elems(:,[3 2]);
0221
0222
0223 p = bottom_z_from_xy ([bn1; bn2]);
0224 skull = build_skull(skull1, skull2, p, maxh);
0225
0226
0227 out.mat_idx = {1:length(out.elems)};
0228 l1 = length(out.boundary);
0229 out.boundary_numbers = uint32(ones(l1,1));
0230 out = merge_meshes(out,skull,0.05);
0231 l1 = length(out.boundary);
0232 out.boundary_numbers(end+1:l1) = 2;
0233 out = merge_meshes(out,CSF,0.025);
0234 l1 = length(out.boundary);
0235 out.boundary_numbers(end+1:l1) = 3;
0236 out = merge_meshes(out,brain);
0237 l1 = length(out.boundary);
0238 out.boundary_numbers(end+1:l1) = 4;
0239 d = sqrt(sum(out.nodes.^2,2));
0240 [~, out.gnd_node] = min(d);
0241 mats = [3 1 4 2];
0242 out.mat_names = materials.names(mats);
0243 out = mk_image(out,.5);
0244
0245 for i = 1:4
0246 out.elem_data(out.fwd_model.mat_idx{i}) = materials.sigma(mats(i));
0247 end
0248
0249 out.name = 'Head image with conductivities';
0250
0251 function [CSF, brain, skull2, bn2] = build_brain_CSF(img,materials,maxh)
0252
0253 cache_path = get_cache_path;
0254 fname = sprintf('brain_C2F_maxh=%.1f.mat', maxh);
0255 fpath = [cache_path filesep fname];
0256 if exist(fpath, 'file')
0257 eidors_msg('@@ Using stored internal meshes.', 3);
0258 load(fpath, 'CSF', 'brain', 'skull2', 'bn2');
0259 return
0260 end
0261
0262 idx = find(strcmp(materials.names, 'CSF'));
0263 CSF = img.fwd_model; CSF.elems = CSF.elems(materials.ref==idx,:);
0264 CSF = remesh_CSF(CSF, maxh);
0265
0266
0267 cond = @(xyz) xyz(:,3) < -122.8;
0268 [skull2, brain, bn2] = split_surface(CSF, cond);
0269 skull2.elems(:,2:3) = skull2.elems(:,[3 2]);
0270 brain.elems(:,2:3) = brain.elems(:,[3 2]);
0271
0272
0273 p = bottom_z_from_xy (bn2);
0274 brain = build_brain(brain,p);
0275
0276 save(fpath, 'CSF', 'brain', 'skull2', 'bn2');
0277
0278 function skull = build_skull(skull1, skull2, p, maxh)
0279
0280 rz = rotation_matrix;
0281 rot = rz*skull1.nodes'; rot = rot';
0282 b1 = rot(:,3) < -121.1 & ...
0283 rot(:,3) < (abs(rot(:,2).^1/11) -123);
0284 rot = rz*skull2.nodes'; rot = rot';
0285 b2 = rot(:,3) < -122.7;
0286 bn1 = skull1.nodes(b1,:);order_loop(bn1,1);
0287 bn2 = skull2.nodes(b2,:);order_loop(bn2,2);
0288 l1 = bn1(:,1:2); l1 = order_loop(l1,-1);
0289 l2 = bn2(:,1:2); l2 = order_loop(l2,-1);
0290 ng_write_opt('meshoptions.fineness',1);
0291 bs = ng_mk_2d_model({l1, l2});
0292 delete('ng.opt');
0293 bs.nodes(:,3) = polyval(p,bs.nodes(:,2));
0294 bs.elems(:,2:3) = bs.elems(:,[3 2]);
0295
0296 skullsrf = merge_meshes(bs,skull1,skull2,min(maxh/5,1));
0297 skull = gmsh_stl2tet(skullsrf,[],[],2);
0298
0299 function brain = build_brain(brain,p)
0300
0301 rz = rotation_matrix;
0302 rot= rz*brain.nodes'; rot = rot';
0303 bot = rot(:,3) < -123;
0304 loop = order_loop(brain.nodes(bot,1:2),-1);
0305 ng_write_opt('meshoptions.fineness',1);
0306 bs = ng_mk_2d_model(loop);
0307 delete('ng.opt');
0308 bs.nodes(:,3) = polyval(p,bs.nodes(:,2));
0309 bs.elems(:,2:3) = bs.elems(:,[3 2]);
0310 mrgd = merge_meshes(brain,bs,0.6);
0311
0312 brain = gmsh_stl2tet(mrgd,[],[],2);
0313
0314 function [upper, lower, bottom_nodes] = split_surface(mdl, fun)
0315
0316
0317 rz = rotation_matrix;
0318 rot= rz*mdl.nodes';
0319 rot = rot';
0320 bot = false(size(mdl.nodes,1),1);
0321 bot(unique(mdl.boundary)) = true;
0322 bot = bot & fun(rot);
0323
0324
0325 bot_idx = find(bot);
0326 on_boundary = ismember(bot_idx, unique(mdl.boundary));
0327 bot = false(size(bot));
0328 bot(bot_idx(on_boundary)) = true;
0329
0330 if eidors_debug('query','mk_head_model_adult:split_surface')
0331 clf
0332 jnk = mdl;
0333 jnk.nodes = rot;
0334 show_fem(jnk)
0335 crop_model([], @(x,y,z) z > -120)
0336 axis tight
0337 hold on
0338 plot3(rot(bot,1),rot(bot,2),rot(bot,3),'o');
0339 hold off
0340 end
0341
0342
0343 b = mdl.boundary;
0344 b = sum(bot(b),2)==3;
0345 mdl.boundary(b,:) = [];
0346 bottom_nodes = mdl.nodes(bot,:);
0347
0348
0349 bnd = unique(mdl.boundary);
0350 [jnk, pos] = min(abs(mdl.nodes(bnd,1))+abs(mdl.nodes(bnd,2))+mdl.nodes(bnd,3));
0351 idx = find_connected_bnd_elems(mdl,bnd(pos));
0352
0353
0354 upper = build_shell(mdl, idx);
0355 log_idx = true(size(mdl.boundary,1),1);
0356 log_idx(idx) = false;
0357 lower = build_shell(mdl, log_idx);
0358
0359 function shell = build_shell(mdl, idx)
0360 shell.nodes = mdl.nodes;
0361 shell.elems = mdl.boundary(idx,:);
0362 shell.type = 'fwd_model';
0363 shell = remove_unused_nodes(shell);
0364 shell.boundary = shell.elems;
0365
0366 function rz = rotation_matrix
0367 tz = -5.8/180 * pi;
0368 rz = [1 0 0; 0 cos(tz) -sin(tz); 0 sin(tz) cos(tz) ;];
0369
0370
0371 function p = bottom_z_from_xy (bn)
0372 p = polyfit(bn(:,2),bn(:,3),6);
0373
0374
0375
0376
0377
0378
0379 function mdl = remesh_CSF(mdl, maxh)
0380
0381 opt.options.meshsize = maxh;
0382 opt.options.curvaturesafety = 3.0;
0383 opt.stloptions.yangle = 20.0;
0384 opt.stloptions.contyangle = 50.0;
0385 opt.stloptions.edgecornerangle = 100;
0386 opt.stloptions.chartangle = 30;
0387
0388 mdl = fix_boundary(mdl);
0389 mdl = ng_stl2tet(mdl, 'moderate', opt);
0390
0391 function idx = find_connected_bnd_elems(mdl,nd)
0392 n_elem = length(mdl.boundary);
0393 todo = false(1,n_elem);
0394 done = false(1,n_elem);
0395 tmp = false(1,length(mdl.nodes));
0396 tmp(nd) = true;
0397
0398 id = any(tmp(mdl.boundary),2);
0399 todo(id) = true;
0400 mdl.elems = mdl.boundary;
0401 mdl = fix_model(mdl, struct('face2edge',true));
0402 f2e = mdl.face2edge;
0403 i = repmat((1:length(mdl.boundary))',1,3);
0404 S = sparse(i(:),f2e(:),true);
0405 C = S*S';
0406 C = logical(spdiags(zeros(n_elem,1),0,C));
0407 while any(todo)
0408 [~, id] = find(C(todo,:));
0409 done(todo) = true;
0410 todo(id) = true;
0411 todo(done) = false;
0412 end
0413 idx = find(done);
0414
0415 function mdle = build_head_elecs(scalp, elec_pos, elec_spec, maxh)
0416 if nargin < 4, maxh = 8; end
0417 ng_opt = get_ng_opt(maxh);
0418 mdle = place_elec_on_surf(scalp, elec_pos, elec_spec,ng_opt);
0419
0420 mdle = fix_boundary(mdle); mdl = mdle;
0421
0422 function [out, img, materials] = get_base_model()
0423 [img, materials] = get_at_model();
0424
0425 cache_path = get_cache_path;
0426 fname = 'base_model.mat';
0427 fpath = [cache_path filesep fname];
0428 if exist(fpath, 'file')
0429 eidors_msg('@@ Using stored base model.', 3);
0430 load(fpath, 'out');
0431 return
0432 end
0433
0434 idx = find(strcmp(materials.names, 'scalp'));
0435 scalp = img.fwd_model; scalp.elems = scalp.elems(materials.ref==idx,:);
0436 out = fix_scalp(scalp);
0437 save(fpath, 'out');
0438
0439 function mdl = fix_scalp(mdl)
0440
0441
0442 mdl = fix_model(mdl);
0443 flip = mdl.elem2face(logical(mdl.boundary_face(mdl.elem2face).*mdl.inner_normal));
0444 mdl.faces(flip,:) = mdl.faces(flip,[1 3 2]);
0445 mdl.normals(flip,:) = -mdl.normals(flip,:);
0446 mdl.boundary = mdl.faces(mdl.boundary_face,:);
0447 mdl.elems = mdl.boundary;
0448
0449 mdl.nodes(1063,:) = mean(mdl.nodes([1063 1065],:));
0450 mdl.nodes(1065,:) = mean(mdl.nodes([109 111], :));
0451 mdl.nodes(1045,:) = mean(mdl.nodes([1043 1045],:));
0452 mdl.nodes(1043,:) = mean(mdl.nodes([101 103], :));
0453
0454 keep = any(reshape(mdl.nodes(mdl.elems,2) < -65 & ...
0455 mdl.nodes(mdl.elems,2) > -75 & ...
0456 mdl.nodes(mdl.elems,3) < -15 & ...
0457 mdl.nodes(mdl.elems,3) > -30 & ...
0458 mdl.nodes(mdl.elems,1) < -10 & ...
0459 mdl.nodes(mdl.elems,1) > -20 ,[],3),2);
0460 h63 = find(any(mdl.elems == 1063,2));
0461 h65 = find(any(mdl.elems == 1065,2));
0462 mdl.elems(h63(2),mdl.elems(h63(2),:)== 1063) = 1065;
0463 mdl.elems(h63(4),mdl.elems(h63(4),:)== 1063) = 1065;
0464 mdl.elems(h65(1),mdl.elems(h65(1),:)== 1065) = 1063;
0465 mdl.elems(h65(3),mdl.elems(h65(3),:)== 1065) = 1063;
0466 mdl.elems(h65(5),mdl.elems(h65(5),:)== 1065) = 1063;
0467 h43 = find(any(mdl.elems == 1043,2));
0468 h45 = find(any(mdl.elems == 1045,2));
0469 mdl.elems(h43(1),mdl.elems(h43(1),:)== 1043) = 1045;
0470 mdl.elems(h43(3),mdl.elems(h43(3),:)== 1043) = 1045;
0471 mdl.elems(h43(5),mdl.elems(h43(5),:)== 1043) = 1045;
0472 mdl.elems(h45(2),mdl.elems(h45(2),:)== 1045) = 1043;
0473 mdl.elems(h45(4),mdl.elems(h45(4),:)== 1045) = 1043;
0474
0475
0476
0477
0478
0479
0480
0481 mdl = gmsh_stl2tet(mdl,[],[],2);
0482 mdl = fix_boundary(mdl);
0483
0484 function opt = get_ng_opt(maxh)
0485 opt.meshoptions.fineness = 6;
0486 opt.options.curvaturesafety = 0.2;
0487
0488
0489 opt.stloptions.yangle = 10;
0490
0491 opt.stloptions.edgecornerangle = 0;
0492
0493 opt.stloptions.outerchartangle = 120;
0494 opt.stloptions.resthchartdistenable = 1;
0495 opt.stloptions.resthchartdistfac = 2.0;
0496 opt.options.meshsize = maxh;
0497 opt.meshoptions.laststep = 'mv';
0498 opt.options.optsteps2d = 5;
0499 opt.options.badellimit = 120;
0500
0501 function [img, materials] = get_at_model()
0502
0503 contrib = 'at-head-mesh'; file = 'SAH262.mat';
0504 load(get_contrib_data(contrib,file),'tri', 'vtx', 'materials');
0505
0506 mdl = eidors_obj('fwd_model','sah262','nodes',vtx,'elems',tri);
0507 img = mk_image(mdl,1);
0508 for i = 1:4
0509 img.elem_data(materials.ref==i) = materials.sigma(i);
0510 end
0511
0512 materials.names = cellstr(materials.names);
0513
0514 function out = get_cache_path
0515 global eidors_objects
0516 out = [eidors_objects.model_cache filesep 'head_model_adult'];
0517 if ~exist(out, 'dir')
0518 mkdir(out);
0519 end
0520
0521
0522 function [e_ar, names] = get_elec_pos(mdl, show)
0523 if nargin < 2, show = false; end
0524
0525 slc = mdl_slice_mesher(mdl,[0 inf inf]);
0526 slc.fwd_model.boundary = find_boundary(slc.fwd_model);
0527
0528
0529 b = unique(slc.fwd_model.boundary);
0530 box = [-87.6 -63.5;-88.2 -28;9.6 -26.8;11.6 -67.5];
0531 b = b(~inpolygon(slc.fwd_model.nodes(b,2),slc.fwd_model.nodes(b,3),box(:,1),box(:,2)));
0532
0533 outln = order_loop(slc.fwd_model.nodes(b,:),-1);
0534 nasion= [-81.754 -8.157];
0535 pp = fourier_fit(outln(:,2:3),51);
0536
0537 fr = linspace(0,0.456,1001); fr(end) = [];
0538 cap = fourier_fit(pp,fr,nasion);
0539 e_ar = [];
0540 names = {};
0541 [e_ar,names] = add_elec(e_ar,names, [0 cap(1,:) ],'Nz');
0542 [e_ar,names] = add_elec(e_ar,names, [0 cap(100,:)],'Fpz');
0543 [e_ar,names] = add_elec(e_ar,names, [0 cap(200,:)],'AFz');
0544 [e_ar,names] = add_elec(e_ar,names, [0 cap(300,:)],'Fz');
0545 [e_ar,names] = add_elec(e_ar,names, [0 cap(400,:)],'FCz');
0546 [e_ar,names] = add_elec(e_ar,names, [0 cap(500,:)],'Cz');
0547 [e_ar,names] = add_elec(e_ar,names, [0 cap(600,:)],'CPz');
0548 [e_ar,names] = add_elec(e_ar,names, [0 cap(700,:)],'Pz');
0549 [e_ar,names] = add_elec(e_ar,names, [0 cap(800,:)],'POz');
0550 [e_ar,names] = add_elec(e_ar,names, [0 cap(900,:)],'Oz');
0551 [e_ar,names] = add_elec(e_ar,names, [0 cap(end,:)],'Iz');
0552
0553 if show
0554 clf
0555
0556 axis equal
0557 hold on
0558 slc.calc_colours.transparency_thresh = -1;
0559 slc.calc_colours.ref_level = 0.25;
0560 show_fem(slc);
0561 end
0562
0563
0564 e1 = get_elec(e_ar,names,'Fpz');
0565 e2 = get_elec(e_ar,names,'Oz');
0566 d = e2 - e1;
0567 t2 = -e1(2)/d(2);
0568 t3 = -e1(3)/d(3);
0569 z = e1(3) + t2*d(3);
0570 y = e1(2) + t3*d(2);
0571 slc = mdl_slice_mesher(mdl,[inf y z]);
0572 slc.calc_colours.transparency_thresh = -1;
0573 slc.calc_colours.ref_level = 0.25;
0574 if show
0575 show_fem(slc);
0576 view(3)
0577 end
0578 slc.fwd_model.boundary = find_boundary(slc.fwd_model);
0579 b = unique(slc.fwd_model.boundary);
0580 outln = order_loop(slc.fwd_model.nodes(b,:),-1);
0581 pp = fourier_fit(outln(:,1:2),51);
0582 fr = linspace(0,1,21); fr(end) = [];
0583 pl = fourier_fit(pp,fr,e1(:,1:2));
0584 pl(:,3) = e1(3) + d(3)*(pl(:,2) - e1(2))/d(2);
0585 pl([1 11],:) = [];
0586 nm = {'Fp1','AF7','F7','FT7','T7','TP7','P7','PO7','O1',...
0587 'O2','PO8','P8','TP8','T8','FT8','F8','AF8','Fp2'};
0588 [e_ar names] = add_elecs(e_ar,names,pl,nm);
0589
0590
0591 nm = {'AF3','AF4'};
0592 [slc el] = define_plane(mdl,e_ar, names, 'AFz','AF7');
0593 el = el([3 6],:);
0594 [e_ar names] = add_elecs(e_ar,names,el,nm);
0595 if show, show_fem(slc); end
0596
0597 nm = {'F9','F5','F3','F1','F2','F4','F6','F10'};
0598 [slc el] = define_plane(mdl,e_ar, names, 'Fz','F7');
0599 [e_ar names] = add_elecs(e_ar,names,el,nm);
0600 if show, show_fem(slc); end
0601
0602 nm = {'FT9','FC5','FC3','FC1','FC2','FC4','FC6','FT10'};
0603 [slc el] = define_plane(mdl,e_ar, names, 'FCz','FT7');
0604 [e_ar names] = add_elecs(e_ar,names,el,nm);
0605 if show, show_fem(slc); end
0606
0607 nm = {'T9','C5','C3','C1','C2','C4','C6','T10'};
0608 [slc el] = define_plane(mdl,e_ar, names, 'Cz','T7');
0609 [e_ar names] = add_elecs(e_ar,names,el,nm);
0610 if show, show_fem(slc); end
0611
0612 nm = {'TP9','CP5','CP3','CP1','CP2','CP4','CP6','TP10'};
0613 [slc el] = define_plane(mdl,e_ar, names, 'CPz','TP7');
0614 [e_ar names] = add_elecs(e_ar,names,el,nm);
0615 if show, show_fem(slc); end
0616
0617 nm = {'P9','P5','P3','P1','P2','P4','P6','P10'};
0618 [slc el] = define_plane(mdl,e_ar, names, 'Pz','P7');
0619 [e_ar names] = add_elecs(e_ar,names,el,nm);
0620 if show, show_fem(slc); end
0621
0622 nm = {'PO3','PO4'};
0623 [slc el] = define_plane(mdl,e_ar, names, 'POz','PO7');
0624 el = el([3 6],:);
0625 [e_ar names] = add_elecs(e_ar,names,el,nm);
0626 if show, show_fem(slc); end
0627 if show
0628 for i = 1:numel(names)
0629 plot_elec(e_ar,names,names{i});
0630 end
0631 hold off
0632 end
0633
0634
0635
0636
0637
0638
0639
0640
0641 function [e_ar names] = add_elecs(e_ar,names, e, name)
0642 for i = 1:size(e,1)
0643 [e_ar names] = add_elec(e_ar,names,e(i,:),name{i});
0644 end
0645
0646 function [e_ar names] = add_elec(e_ar,names,e,name);
0647 e_ar(end+1,:) = e;
0648 names{end+1} = name;
0649
0650
0651 function [e1] = get_elec(e_ar,names,name)
0652 idx = find(ismember(names,name));
0653 e1 = e_ar(idx,:);
0654
0655 function plot_elec(e_ar,names,name)
0656 e = get_elec(e_ar,names,name);
0657 if size(e,2) == 3
0658 plot3(e(:,1),e(:,2),e(:,3),'bo','MarkerFaceColor','b','MarkerSize',10);
0659 else
0660 plot(e(:,1),e(:,2),'bo','MarkerFaceColor','b','MarkerSize',10);
0661 end
0662
0663 function [slc, e] = define_plane(mdl, e_ar,names,n1, n2);
0664 e1 = get_elec(e_ar,names,n1);
0665 e2 = get_elec(e_ar,names,n2);
0666 d = e2 - e1;
0667 t2 = -e1(2)/d(2);
0668 t3 = -e1(3)/d(3);
0669 z = e1(3) + t2*d(3);
0670 y = e1(2) + t3*d(2);
0671 slc = mdl_slice_mesher(mdl,[inf y z]);
0672
0673 slc.calc_colours.transparency_thresh = -1;
0674 slc.calc_colours.ref_level = 0.25;
0675 slc.fwd_model.boundary = find_boundary(slc.fwd_model);
0676
0677
0678 b = unique(slc.fwd_model.boundary);
0679 box = [ -16.0484 -19.0245
0680 22.6613 -15.7946
0681 23.7903 -92.9910
0682 -20.0806 -91.6990];
0683 b = b(~inpolygon(slc.fwd_model.nodes(b,1),slc.fwd_model.nodes(b,3),box(:,1),box(:,2)));
0684 outln = order_loop(slc.fwd_model.nodes(b,:),-1);
0685 pp = fourier_fit(outln(:,[1 3]),51);
0686 fr = linspace(0,1,10001); fr(end) = [];
0687 pl = fourier_fit(pp,fr,e1(:,[1 3]));
0688
0689 di = sqrt(sum((pl - repmat(e2([1 3]),length(pl),1)).^2,2));
0690 [jnk, p] = min(di);
0691 p = p/10000;
0692 if p > 0.5, p = 1-p; end
0693 fr = [-5 -3 -2 -1 1 2 3 5] ./4 .* p;
0694 el = fourier_fit(pp, fr, e1(:,[1 3]));
0695 e(:,[1 3]) = el;
0696 e(:,2) = e1(2) + d(2).*(e(:,3) - e1(3))./d(3);
0697
0698
0699
0700
0701
0702
0703 function do_unit_test
0704
0705 warning('off','EIDORS:License');
0706
0707
0708
0709 fmdl = mk_head_model_adult();
0710 unit_test_cmp('fwd_model', fmdl.type, 'fwd_model');
0711
0712
0713
0714
0715 elec_pos = [8, 0, 0, 10];
0716 elec_shape = [3, 0, 1];
0717 maxh = 10;
0718 img1 = mk_head_model_adult(elec_pos, elec_shape);
0719 img2 = mk_head_model_adult(elec_pos, elec_shape, maxh);
0720 unit_test_cmp('custom + maxh', num_elems(img1) > num_elems(img2), true);
0721
0722
0723
0724
0725
0726
0727
0728
0729
0730
0731
0732
0733
0734 img10 = mk_head_model_adult('10-10');
0735 unit_test_cmp('10-10 # elecs', numel(img10.fwd_model.electrode),73);
0736 img1 = mk_head_model_adult('1x32_ring');
0737 unit_test_cmp('1x32 # elecs', numel(img1.fwd_model.electrode),32);
0738 img1 = mk_head_model_adult('no_elecs');
0739 unit_test_cmp('no elecs', isfield(img1.fwd_model, 'electrode'),false);
0740 unit_test_cmp('materials', numel(img1.fwd_model.mat_idx), 4);
0741
0742
0743
0744
0745
0746
0747
0748
0749
0750
0751
0752 fmdl = mk_head_model_adult(elec_pos, elec_shape, maxh, 'void');
0753 unit_test_cmp('void type', fmdl.type, 'fwd_model');
0754 img1 = mk_head_model_adult('10-10', 'fine');
0755 unit_test_cmp('fine', num_elems(img10) < num_elems(img1), true);
0756
0757 warning('on','EIDORS:License');