mat_idx_to_electrode

 MAT_IDX_TO_ELECTRODE: create electrodes from mat_idx values
 fmdl = mat_idx_to_electrode(fmdl, mat_idxes)
 fmdl: input and output fmdl
 Options: fmdl.mat_idx_to_electrode.z_contact  (z_contact value to use ... or default)
 mat_idxes = cell array of mat_idxes => convert electrode
    example mat_idxes = {1:2, 5, [12,14]}

 by default, adds a 'faces'-type electrode (i.e. an
    electrode defined in terms of its faces)
 Parameter
   fmdl.mat_idx_to_electrode.nodes_electrode = true
      Add a 'nodes'-type electrode

 To work with an image object, use:
  [img.fwd_model,rm_elems] = mat_idx_to_electrode(img.fwd_model,{mat_idxes});
  img.elem_data(rm_elems) = [];

 Notes:
   mat_idx regions cannot be directly touching each other.

0001 function [fmdl,rm_elems] = mat_idx_to_electrode(fmdl, mat_idxes)
0002 % MAT_IDX_TO_ELECTRODE: create electrodes from mat_idx values
0003 % fmdl = mat_idx_to_electrode(fmdl, mat_idxes)
0004 % fmdl: input and output fmdl
0005 % Options: fmdl.mat_idx_to_electrode.z_contact  (z_contact value to use ... or default)
0006 % mat_idxes = cell array of mat_idxes => convert electrode
0007 %    example mat_idxes = {1:2, 5, [12,14]}
0008 %
0009 % by default, adds a 'faces'-type electrode (i.e. an
0010 %    electrode defined in terms of its faces)
0011 % Parameter
0012 %   fmdl.mat_idx_to_electrode.nodes_electrode = true
0013 %      Add a 'nodes'-type electrode
0014 %
0015 % To work with an image object, use:
0016 %  [img.fwd_model,rm_elems] = mat_idx_to_electrode(img.fwd_model,{mat_idxes});
0017 %  img.elem_data(rm_elems) = [];
0018 %
0019 % Notes:
0020 %   mat_idx regions cannot be directly touching each other.
0021 
0022 % (C) 2019 Andy Adler. License: GPL v2 or v3.
0023 % $Id: mat_idx_to_electrode.m 5968 2019-06-05 16:58:26Z aadler $
0024 
0025 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0026 
0027 elec_faces = true;
0028 try 
0029    elec_faces = ~fmdl.mat_idx_to_electrode.nodes_electrode;
0030 end
0031 
0032 rm_elems = [];
0033 for i = 1:length(mat_idxes)
0034    if elec_faces 
0035       [fmdl,rm_elemi] = create_electrode_faces_from_mat_idx(fmdl, mat_idxes{i});
0036    else
0037       [fmdl,rm_elemi] = create_electrode_nodes_from_mat_idx(fmdl, mat_idxes{i});
0038    end
0039    rm_elems = union(rm_elems,rm_elemi); 
0040 end
0041 
0042 fmdl = remove_unused_nodes(fmdl);
0043 if any(fmdl.boundary(:) == 0)
0044    keyboard
0045 end
0046 
0047 % This code adds the new electrode as a elec(i).nodes
0048 % adds electrode to the end
0049 % femobj is removed elems
0050 function [fmdl,femobj] = create_electrode_nodes_from_mat_idx(fmdl,nmat_idx);
0051    zc = 1e-5;
0052    try zc = fmdl.mat_idx_to_electrode.z_contact;
0053    end
0054    femobj = vertcat(fmdl.mat_idx{nmat_idx});
0055    faces = calc_elec_faces(fmdl.elems,femobj);
0056 
0057    femobjnodes = fmdl.elems(femobj,:);
0058    [fmdl,faces] = rm_elems( fmdl, femobj, faces);
0059 
0060    fmdl.boundary = unique([fmdl.boundary;faces],'rows');
0061 
0062    vt = find_bdynodes(fmdl,femobjnodes);
0063    elstr =  struct('nodes',vt(:)','z_contact',zc);
0064    fmdl = add_elec(fmdl, elstr);
0065 
0066 function [fmdl,faces] = rm_elems( fmdl, femobj, faces);
0067    % fix the mat_idx object, since we remove femobj
0068    for i=1:length(fmdl.mat_idx) 
0069       els = false(num_elems(fmdl),1);
0070       els(fmdl.mat_idx{i}) = true;
0071       els(femobj) = [];
0072       fmdl.mat_idx{i} = find(els);
0073    end
0074    fmdl.elems(femobj,:) = [];
0075 
0076    if nargin>2
0077       inels = reshape(... 
0078            ismember(faces(:),fmdl.elems(:)), size(faces));
0079       faces(any(~inels,2),:) = [];
0080    end
0081 
0082 % This code adds the new electrode as a elec(i).faces
0083 % adds electrode to the end
0084 % femobj is removed elems
0085 function [fmdl,femobj] = create_electrode_faces_from_mat_idx(fmdl,nmat_idx);
0086    zc = 1e-5;
0087    try zc = fmdl.mat_idx_to_electrode.z_contact;
0088    end
0089    femobj = vertcat(fmdl.mat_idx{nmat_idx});
0090    faces = calc_elec_faces(fmdl.elems,femobj);
0091    [fmdl,faces] = rm_elems( fmdl, femobj, faces);
0092    fmdl.boundary = unique([fmdl.boundary;faces],'rows');
0093 
0094    elstr =  struct('nodes',[],'z_contact',zc,'faces',faces);
0095 
0096    fmdl = add_elec(fmdl, elstr);
0097 
0098 function fmdl = add_elec(fmdl, elstr)
0099    if isfield(fmdl,'electrode')
0100 %    Stupid matlab forces you to add this
0101      if isfield(elstr,'faces') && ~isfield(fmdl.electrode,'faces');
0102         [fmdl.electrode(:).faces] = deal([]);
0103      end
0104      fmdl.electrode(end+1) = elstr;
0105    else
0106      fmdl.electrode(    1) = elstr;
0107    end
0108 
0109 function faces = calc_elec_faces(elems,femobj);
0110    thisels = elems(femobj,:);
0111    switch size(thisels,2)  % dimentions
0112      case 2; % 1D
0113        allfaces = thisels;
0114      case 3; % 2D
0115        allfaces = [thisels(:,[1,2]);
0116                    thisels(:,[2,3]);
0117                    thisels(:,[1,3])];
0118      case 4; % 3D
0119        allfaces = [thisels(:,[1,2,3]);
0120                    thisels(:,[1,2,4]);
0121                    thisels(:,[1,3,4]);
0122                    thisels(:,[2,3,4])];
0123      otherwise;
0124        error 'Dimensions of elems';
0125    end
0126    allfaces = sort(allfaces,2);
0127    % remove all faces which exist more than once
0128    [faces,ia,ib] = unique(allfaces,'rows');
0129    exist_faces = sparse(1,ib,1);
0130    faces(exist_faces>1,:) = [];
0131 
0132 function vt = find_bdynodes(fmdl,femobjnodes)
0133 % Slow way
0134 %  vt = intersect(femobjnodes,find_boundary(fmdl));
0135 % Fast: pre-process
0136    usenodes = reshape( ismember( ...
0137       fmdl.elems, femobjnodes), size(fmdl.elems));
0138    fmdl.elems(~any(usenodes,2),:) = [];
0139    vt = intersect(femobjnodes,find_boundary(fmdl));
0140       
0141 
0142 function do_unit_test
0143    clf; subplot(221);
0144    do_unit_test_2d(true)
0145    do_unit_test_2d(false)
0146    subplot(223);
0147 %  do_unit_test_3d_netgen
0148 
0149 function do_unit_test_2d(nodes_electrode)
0150    fmdl = getfield(mk_common_model('a2c2',1),'fwd_model');
0151    fmdl.mat_idx_to_electrode.nodes_electrode = nodes_electrode;
0152    fmdl.electrode(1).nodes = 26:41;
0153    fmdl.mat_idx{1} = [1:4];
0154    fmdl= mat_idx_to_electrode(fmdl, {1});
0155    unit_test_cmp('a2c2-01',num_elems(fmdl),64-4);
0156    if nodes_electrode;
0157       unit_test_cmp('a2c2-02',length(fmdl.electrode(2).nodes),4);
0158    else
0159       unit_test_cmp('a2c2-02a',length(fmdl.electrode(2).nodes),0);
0160       unit_test_cmp('a2c2-02b',size(fmdl.electrode(2).faces),[4,2]);
0161    end
0162 
0163    fmdl.stimulation = stim_meas_list([1,2,1,2]);
0164    img = mk_image(fmdl,1);
0165    img.fwd_solve.get_all_meas = true;
0166    vh = fwd_solve(img);
0167 
0168    img = rmfield(img,'elem_data');
0169    img.node_data = vh.volt;
0170    
0171    img.fwd_model = rmfield(img.fwd_model,'electrode');
0172    show_fem(img,[0,0,2]);
0173 %  plot(sqrt(sum(fmdl.nodes.^2,2)),vh.volt,'*');
0174    unit_test_cmp('a2c2-03', std(vh.volt(13:24)),0,0.002);
0175    vd = mean(vh.volt(5:12)) - mean(vh.volt(13:24));
0176    unit_test_cmp('a2c2-03', vd,0.065251,1e-6);
0177 
0178 
0179 
0180 function do_unit_test_3d_netgen
0181    extra={'ball_inside','ball_surface', [ ...
0182           'solid ball_inside  = sphere(-0.4,0,0.5;0.05);' ...
0183           'solid ball_surface = sphere( 0.4,0,1.0;0.05) -maxh=.005;' ...
0184           ]};
0185    fmdl= ng_mk_cyl_models(1,[8,.5],[.1],extra); 
0186    fmdl.mat_idx_to_electrode.nodes_electrode = true;
0187 
0188    fmd_= mat_idx_to_electrode(fmdl, {2,3});
0189    unit_test_cmp('cyl_models 01',num_elecs(fmd_),10);
0190 
0191    fmd_= mat_idx_to_electrode(fmdl, {2:3});
0192    unit_test_cmp('cyl_models 02',num_elecs(fmd_),9);
0193 
0194    fmd_= mat_idx_to_electrode(fmdl, {2});
0195    unit_test_cmp('cyl_models 03',num_elecs(fmd_),9);
0196 
0197    img = mk_image(fmd_,1);
0198    img.elem_data(fmd_.mat_idx{3})= 1.1;
0199    img.calc_colours.ref_level = 1;
0200    unit_test_cmp('cyl_models 04',num_elecs(img),9);
0201 
0202    img = mk_image(fmdl,1);
0203    img.elem_data(vertcat(fmdl.mat_idx{2:3}))= 1.1;
0204    [img.fwd_model,rm_elems]= mat_idx_to_electrode( ...
0205         img.fwd_model, {2});
0206    img.elem_data(rm_elems) = [];
0207    unit_test_cmp('cyl_models 05',num_elecs(img),9);
0208    vol = get_elem_volume(img);
0209    vvs = sum(vol(find(img.elem_data == 1.1)));
0210    unit_test_cmp('cyl_models 05a',vvs,4/3*pi*.05^3,1e-5);
0211   
0212 
0213 
0214    img.calc_colours.ref_level = 1;
0215    show_fem(img); view(3,12);
0216