mk_head_model_adult

PURPOSE ^

MK_HEAD_MODEL_ADULT Adult head model

SYNOPSIS ^

function out = mk_head_model_adult(varargin)

DESCRIPTION ^

MK_HEAD_MODEL_ADULT Adult head model 

 MK_HEAD_MODEL_ADULT provides predefined and custom head FEMs based on
 models contributed by Andrew Tizzard. You may be asked to download
 missing files.

 FMDL = MK_HEAD_MODEL_ADULT() 
 Returns a head fwd_model with no electrodes or tissues

 FMDL = MK_HEAD_MODEL_ADULT(elec_pos, elec_shape)
 FMDL = MK_HEAD_MODEL_ADULT(elec_pos, elec_shape, maxh)
 Allows specifying electrode configuration, see PLACE_ELEC_ON_SURF.
 Default maxh = 8.

 FMDL = MK_HEAD_MODEL_ADULT(elecstr)
 A predefined model with round electrodes (radius 2 mm). 
    '10-10'         - 73 labeled electrodes (10-10 placement system),
                      no stimulation.
    '2x16_planar'   - 2 rings of 16 electrodes, planar stimulation
    '2x16_odd-even' - 2 rings, each stimulation is cross-plane 
    '2x16_square'   - 2 rings, every second stimulation is cross-plane
    '2x16_adjacent' - like square, but with adjacent stimulation
    '1x32_ring'     - single ring of 32 electrodes
    'no_elecs'      - no electrodes

 IMG = MK_HEAD_MODEL_ADULT(... , tissuespec)
 Returns a head model containing separate regions defining scalp, skull, 
 CSF and brain as an image structure with conductivitity values. 
 TISSUESPEC must be one of:
    'coarse'        - (defualt) tissue surfaces with maxh = 5mm
    'fine'          - tissue surfaces with maxh = 2mm
    'void'          - no internal tissues (void), useful for testing
                      electrode positions

 MK_HEAD_MODEL('show_10-10') creates a figure showing the positions and
 labels of the 10-10 electrodes.

 
 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
 DOI: 10.1109/EMBC.2013.6611026

 CITATION_REQUEST:
 TITLE: Improving the Finite Element Forward Model of the Human Head by 
 Warping using Elastic Deformation
 AUTHOR: A Tizzard and RH Bayford
 JOURNAL: Physiol Meas
 YEAR: 2007
 VOL: 28
 PAGE: S163-S182
 DOI: 10.1088/0967-3334/28/7/S13

 See also PLAC_ELEC_ON_SURF

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function out = mk_head_model_adult(varargin)
0002 %MK_HEAD_MODEL_ADULT Adult head model
0003 %
0004 % MK_HEAD_MODEL_ADULT provides predefined and custom head FEMs based on
0005 % models contributed by Andrew Tizzard. You may be asked to download
0006 % missing files.
0007 %
0008 % FMDL = MK_HEAD_MODEL_ADULT()
0009 % Returns a head fwd_model with no electrodes or tissues
0010 %
0011 % FMDL = MK_HEAD_MODEL_ADULT(elec_pos, elec_shape)
0012 % FMDL = MK_HEAD_MODEL_ADULT(elec_pos, elec_shape, maxh)
0013 % Allows specifying electrode configuration, see PLACE_ELEC_ON_SURF.
0014 % Default maxh = 8.
0015 %
0016 % FMDL = MK_HEAD_MODEL_ADULT(elecstr)
0017 % A predefined model with round electrodes (radius 2 mm).
0018 %    '10-10'         - 73 labeled electrodes (10-10 placement system),
0019 %                      no stimulation.
0020 %    '2x16_planar'   - 2 rings of 16 electrodes, planar stimulation
0021 %    '2x16_odd-even' - 2 rings, each stimulation is cross-plane
0022 %    '2x16_square'   - 2 rings, every second stimulation is cross-plane
0023 %    '2x16_adjacent' - like square, but with adjacent stimulation
0024 %    '1x32_ring'     - single ring of 32 electrodes
0025 %    'no_elecs'      - no electrodes
0026 %
0027 % IMG = MK_HEAD_MODEL_ADULT(... , tissuespec)
0028 % Returns a head model containing separate regions defining scalp, skull,
0029 % CSF and brain as an image structure with conductivitity values.
0030 % TISSUESPEC must be one of:
0031 %    'coarse'        - (defualt) tissue surfaces with maxh = 5mm
0032 %    'fine'          - tissue surfaces with maxh = 2mm
0033 %    'void'          - no internal tissues (void), useful for testing
0034 %                      electrode positions
0035 %
0036 % MK_HEAD_MODEL('show_10-10') creates a figure showing the positions and
0037 % labels of the 10-10 electrodes.
0038 %
0039 %
0040 % CITATION_REQUEST:
0041 % TITLE: FEM Electrode Refinement for Electrical Impedance Tomography
0042 % AUTHOR: B Grychtol and A Adler
0043 % JOURNAL: Engineering in Medicine and Biology Society (EMBC), 2013 Annual
0044 % International Conference of the IEEE
0045 % YEAR: 2013
0046 % DOI: 10.1109/EMBC.2013.6611026
0047 %
0048 % CITATION_REQUEST:
0049 % TITLE: Improving the Finite Element Forward Model of the Human Head by
0050 % Warping using Elastic Deformation
0051 % AUTHOR: A Tizzard and RH Bayford
0052 % JOURNAL: Physiol Meas
0053 % YEAR: 2007
0054 % VOL: 28
0055 % PAGE: S163-S182
0056 % DOI: 10.1088/0967-3334/28/7/S13
0057 %
0058 % See also PLAC_ELEC_ON_SURF
0059 
0060 % (C) 2012-2024 Bartek Grychtol. License: GPL v2 or v3
0061 % $Id: mk_head_model_adult.m 7031 2024-11-28 21:12:36Z bgrychtol $
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; % Returns a head fwd_model with no electrodes or tissues
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       % nothing
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); % impose standard field order
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 % extract the outer skull surface from the scalp model with electrodes
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]); % flip normals
0221 
0222 % close surfaces from the bottom and create volume meshes
0223 p = bottom_z_from_xy ([bn1; bn2]); 
0224 skull = build_skull(skull1, skull2, p, maxh);
0225 
0226 % put it all together one by one
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 % show_fem(out)
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 % extract individual surfaces from the CSF model
0267 cond = @(xyz) xyz(:,3) < -122.8; % condition to apply to the *rotated* model
0268 [skull2, brain, bn2] = split_surface(CSF, cond);
0269 skull2.elems(:,2:3) = skull2.elems(:,[3 2]); % flip normals
0270 brain.elems(:,2:3) = brain.elems(:,[3 2]); % flip normals
0271 
0272 % close surface from the bottom and create volume meshe
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 % treat the bottom as an xy plane and mesh it
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 % merge surfaces and mesh the inside
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]); % reorient normals
0310 mrgd = merge_meshes(brain,bs,0.6);
0311 % stl_write(mrgd,'brain.stl');
0312 brain = gmsh_stl2tet(mrgd,[],[],2);
0313 
0314 function [upper, lower, bottom_nodes] = split_surface(mdl, fun)
0315 
0316 % find indecies of the bottom surface nodes
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 % remove the bottom surface from the boundary
0343 b = mdl.boundary;
0344 b = sum(bot(b),2)==3;
0345 mdl.boundary(b,:) = [];
0346 bottom_nodes = mdl.nodes(bot,:); % save bottom nodes for later
0347 
0348 % separate surfaces
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 % package
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 % find an equation relating y and z on the bottom surface
0371 function p = bottom_z_from_xy (bn)
0372 p = polyfit(bn(:,2),bn(:,3),6);
0373 % clf
0374 % plot(bn(:,2),bn(:,3),'.')
0375 % hold on
0376 % pp = polyval(p,-100:.1:100);
0377 % plot(-100:.1:100,pp,'r')
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 % boundary elements that contain that node
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'; % face2face connectivity matrix
0406 C = logical(spdiags(zeros(n_elem,1),0,C)); % remove diagonal
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 % fix orientation of the boundary
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); % re-meshes after merging duplicated surface nodes
0437 save(fpath, 'out');
0438 
0439 function mdl = fix_scalp(mdl)
0440 
0441 %skull has some funny holes in the eyes.... fix manually
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 % nodes 1065 and 1063 should be one
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 % mdl.boundary = mdl.elems(keep,:);
0476 % clf
0477 % show_fem(mdl);
0478 % hold on
0479 
0480 % we've only fixed the surface. Re-mesh the inside
0481 mdl = gmsh_stl2tet(mdl,[],[],2);
0482 mdl = fix_boundary(mdl);
0483 
0484 function opt = get_ng_opt(maxh)
0485    opt.meshoptions.fineness = 6; % some options have no effect without this
0486    opt.options.curvaturesafety = 0.2;
0487    % small yangle preserves the original mesh, large encourages smoother
0488    % surface with nicer spreading of refinement
0489    opt.stloptions.yangle = 10;
0490  %    opt.stloptions.contyangle = 20;
0491    opt.stloptions.edgecornerangle = 0;
0492 %    opt.stloptions.chartangle = 0;
0493    opt.stloptions.outerchartangle = 120;
0494    opt.stloptions.resthchartdistenable = 1;
0495    opt.stloptions.resthchartdistfac = 2.0; % encourages slower increase of element size
0496    opt.options.meshsize = maxh;
0497    opt.meshoptions.laststep = 'mv'; % don't need volume optimization
0498    opt.options.optsteps2d =  5; % but we can up surface optimization
0499    opt.options.badellimit = 120; % decrease the maximum allowed angle
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 %% ELECTRODE POSITIONS
0522 function [e_ar, names] = get_elec_pos(mdl, show)
0523 if nargin < 2, show = false; end
0524 % cut through the nose:
0525 slc = mdl_slice_mesher(mdl,[0 inf inf]);
0526 slc.fwd_model.boundary = find_boundary(slc.fwd_model);
0527 
0528 % exclude the nasal cavity
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 % define Nasion
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 % 0.457 seems to be a realistic place for inion, the other landmark
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    % show_fem(mdl);
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 % calculate a "horizontal plane between Nz and Iz
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 % load epos128;
0635 % phi = sign(epos(:,1))*90 - epos(:,2); phi = pi*phi/180;
0636 % the = epos(:,1); the(the>0) = -the(the>0); the = pi*the/180;
0637 % z = sin(the);
0638 % x = cos(the).*sin(phi);
0639 % y = cos(the).*cos(phi);
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 % function [x y z] = get_elec_cart(mdl,incl, azim)
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 % show_fem(slc);
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 %exclude nasal cavity
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 % find the position of the other electrode
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 % plot3(e(:,1),e(:,2),e(:,3),'o')
0699 % plot3(e(1,1),e(1,2),e(1,3),'o','MarkerFaceColor','b')
0700 
0701 
0702 %% UNIT TEST
0703 function do_unit_test
0704 
0705 warning('off','EIDORS:License');
0706 
0707 % FMDL = MK_HEAD_MODEL_ADULT()
0708 % Returns a head fwd_model with no electrodes or tissues
0709 fmdl = mk_head_model_adult();
0710 unit_test_cmp('fwd_model', fmdl.type, 'fwd_model');
0711 %
0712 % FMDL = MK_HEAD_MODEL_ADULT(elec_pos, elec_shape)
0713 % FMDL = MK_HEAD_MODEL_ADULT(elec_pos, elec_shape, maxh)
0714 % Allows specifying electrode configuration, see PLACE_ELEC_ON_SURF
0715 elec_pos = [8, 0, 0, 10]; % [N, equal angles, z1, z2];
0716 elec_shape = [3, 0, 1]; % [radius, 0, maxh]
0717 maxh = 10;
0718 img1 = mk_head_model_adult(elec_pos, elec_shape); % default max is 8
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 % FMDL = MK_HEAD_MODEL_ADULT(elecstr)
0724 % A predefined model with round electrodes (radius 2 mm).
0725 %    '10-10'         - 73 labeled electrodes (10-10 placement system),
0726 %                      no stimulation.
0727 %    '2x16_planar'   - 2 rings of 16 electrodes, planar stimulation
0728 %    '2x16_odd-even' - 2 rings, each stimulation is cross-plane
0729 %    '2x16_square'   - 2 rings, every second stimulation is cross-plane
0730 %    '2x16_adjacent' - like square, but with adjacent stimulation
0731 %    '1x32_ring'     - single ring of 32 electrodes
0732 %    'no_elecs'      - no electrodes
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 % IMG = MK_HEAD_MODEL_ADULT(... , tissuespec)
0744 % Returns a head model containing separate regions defining scalp, skull,
0745 % CSF and brain as an image structure with conductivitity values.
0746 % TISSUESPEC must be one of:
0747 %    'coarse'        - (defualt) tissue surfaces with maxh = 5mm
0748 %    'fine'          - tissue surfaces with maxh = 2mm
0749 %    'void'          - no internal tissues (void), useful for testing
0750 %                      electrode positions
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');

Generated on Sun 29-Dec-2024 11:41:59 by m2html © 2005