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 6096 2021-03-11 21:45:34Z aadler $
0024 
0025 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0026 
0027 % Calculate electrodes on faces or nodes
0028 elec_faces = true;
0029 try 
0030    elec_faces = ~fmdl.mat_idx_to_electrode.nodes_electrode;
0031 end
0032 
0033 rm_elems = [];
0034 for i = 1:length(mat_idxes)
0035    if elec_faces 
0036       [fmdl,rm_elemi] = create_electrode_faces_from_mat_idx(fmdl, mat_idxes{i});
0037    else
0038       [fmdl,rm_elemi] = create_electrode_nodes_from_mat_idx(fmdl, mat_idxes{i});
0039    end
0040    rm_elems = union(rm_elems,rm_elemi); 
0041 end
0042 
0043 fmdl = remove_unused_nodes(fmdl);
0044 fmdl = remove_unused_boundary(fmdl);
0045 if any(fmdl.boundary(:) == 0)
0046    eidors_msg('WARNING: PROBLEM WITH BOUNDARY',1)
0047 end
0048 for i=1:num_elecs(fmdl)
0049     zeroface = false;
0050     try  
0051        zeroface = any(fmdl.electrode(i).faces(:) == 0);
0052     end
0053     if zeroface || any(fmdl.electrode(i).nodes(:) == 0)
0054        eidors_msg('WARNING: PROBLEM WITH BOUNDARY %d',i,1)
0055        keyboard
0056     end
0057 end
0058 
0059 % This code adds the new electrode as a elec(i).nodes
0060 % adds electrode to the end
0061 % femobj is removed elems
0062 function [fmdl,femobj] = create_electrode_nodes_from_mat_idx(fmdl,nmat_idx);
0063    zc = 1e-5;
0064    try zc = fmdl.mat_idx_to_electrode.z_contact;
0065    end
0066    femobj = vertcat(fmdl.mat_idx{nmat_idx});
0067    faces = calc_elec_faces(fmdl.elems,femobj);
0068 
0069    femobjnodes = fmdl.elems(femobj,:);
0070    [fmdl,faces] = rm_elems( fmdl, femobj, faces);
0071 
0072    % Add to the boundary with the boundary of the new electrode
0073    fmdl.boundary = unique([fmdl.boundary;faces],'rows');
0074 
0075    % vtx= intersect(fmdl.boundary,femobjnodes);
0076    % This doesn't seem needed, just unique faces
0077    % vt = find_bdynodes(fmdl,femobjnodes);
0078    vt = unique(faces);
0079 
0080    elstr =  struct('nodes',vt(:)','z_contact',zc);
0081    fmdl = add_elec(fmdl, elstr);
0082 
0083 function fmdl = remove_unused_boundary(fmdl);
0084    elems= fmdl.elems;
0085    fmdl.boundary = unique(sort(fmdl.boundary,2),'rows');
0086    switch mdl_dim(fmdl)
0087       case 2; selem = [sort(elems(:,[1,2]),2);
0088                        sort(elems(:,[1,3]),2);
0089                        sort(elems(:,[2,3]),2)];
0090       case 3; selem = [sort(elems(:,[1,2,3]),2);
0091                        sort(elems(:,[1,2,4]),2);
0092                        sort(elems(:,[1,3,4]),2);
0093                        sort(elems(:,[2,3,4]),2)];
0094       otherwise; error('Dims must be 2 or 3');
0095    end
0096    [~,idx] = setdiff(fmdl.boundary, selem,'rows');
0097 %  disp( fmdl.boundary(idx,:) )
0098    fmdl.boundary(idx,:) = [];
0099 
0100 function [fmdl,faces] = rm_elems( fmdl, femobj, faces);
0101    % fix the mat_idx object, since we remove femobj
0102    for i=1:length(fmdl.mat_idx) 
0103       els = false(num_elems(fmdl),1);
0104       els(fmdl.mat_idx{i}) = true;
0105       els(femobj) = [];
0106       fmdl.mat_idx{i} = find(els);
0107    end
0108    fmdl.elems(femobj,:) = [];
0109 
0110 % Taken care of by remove_unused_boundary
0111 %  % remove faces that are no longer on the boundary
0112 %  if nargin>2
0113 %     inels = reshape(...
0114 %          ismember(faces(:),fmdl.elems(:)), size(faces));
0115 %     faces(any(~inels,2),:) = [];
0116 %     sfaces = sort(faces,2);
0117 %     selem1 = sort(fmdl.elems(:,[1,2,3]),2);
0118 %     selem2 = sort(fmdl.elems(:,[1,2,4]),2);
0119 %     selem3 = sort(fmdl.elems(:,[1,3,4]),2);
0120 %     selem4 = sort(fmdl.elems(:,[2,3,4]),2);
0121 %     [~,idx] = setdiff(sfaces,[selem1;selem2;selem3;selem4],'rows');
0122 %     faces(idx,:) = [];
0123 %  end
0124 
0125 % This code adds the new electrode as a elec(i).faces
0126 % adds electrode to the end
0127 % femobj is removed elems
0128 function [fmdl,femobj] = create_electrode_faces_from_mat_idx(fmdl,nmat_idx);
0129    zc = 1e-5;
0130    try zc = fmdl.mat_idx_to_electrode.z_contact;
0131    end
0132    femobj = vertcat(fmdl.mat_idx{nmat_idx});
0133    faces = calc_elec_faces(fmdl.elems,femobj);
0134    [fmdl,faces] = rm_elems( fmdl, femobj, faces);
0135    fmdl.boundary = unique([fmdl.boundary;faces],'rows');
0136 
0137    elstr =  struct('nodes',[],'z_contact',zc,'faces',faces);
0138 
0139    fmdl = add_elec(fmdl, elstr);
0140 
0141 function fmdl = add_elec(fmdl, elstr)
0142    if isfield(fmdl,'electrode')
0143 %    Stupid matlab forces you to add this
0144      if isfield(elstr,'faces') && ~isfield(fmdl.electrode,'faces');
0145         [fmdl.electrode(:).faces] = deal([]);
0146      end
0147      fmdl.electrode(end+1) = elstr;
0148    else
0149      fmdl.electrode(    1) = elstr;
0150    end
0151 
0152 % Detect outer face of electrode by
0153 %  removing all tet faces which exist more
0154 %  than once
0155 function faces = calc_elec_faces(elems,femobj);
0156    thisels = elems(femobj,:);
0157    switch size(thisels,2)  % dimentions
0158      case 2; % 1D
0159        allfaces = thisels;
0160      case 3; % 2D
0161        allfaces = [thisels(:,[1,2]);
0162                    thisels(:,[2,3]);
0163                    thisels(:,[1,3])];
0164      case 4; % 3D
0165        allfaces = [thisels(:,[1,2,3]);
0166                    thisels(:,[1,2,4]);
0167                    thisels(:,[1,3,4]);
0168                    thisels(:,[2,3,4])];
0169      otherwise;
0170        error 'Dimensions of elems';
0171    end
0172    allfaces = sort(allfaces,2);
0173    % remove all faces which exist more than once
0174    [faces,ia,ib] = unique(allfaces,'rows');
0175    exist_faces = sparse(1,ib,1);
0176    faces(exist_faces>1,:) = [];
0177 
0178 % This functions doesn't seem necessary ... remove
0179 % AA: mar 2021
0180 %function vt = find_bdynodes(fmdl,femobjnodes)
0181 %% Slow way
0182 %   vt = intersect(femobjnodes,find_boundary(fmdl));
0183 %% Fast: pre-process
0184 %   usenodes = reshape( ismember( ...
0185 %      fmdl.elems, femobjnodes), size(fmdl.elems));
0186 %   fmdl.elems(~any(usenodes,2),:) = [];
0187 %   vt = intersect(femobjnodes,find_boundary(fmdl));
0188       
0189 
0190 function do_unit_test
0191    clf; do_unit_test_contained_shape(true)
0192    clf; do_unit_test_contained_shape(false)
0193    return
0194    clf; subplot(221);
0195    do_unit_test_2d(true)
0196    do_unit_test_2d(false)
0197    subplot(223);
0198    do_unit_test_3d_netgen(true)
0199    do_unit_test_3d_netgen(false)
0200    do_unit_test_3d_netgen2(true)
0201    do_unit_test_3d_netgen2(false)
0202 
0203 function do_unit_test_2d(nodes_electrode)
0204    fmdl = getfield(mk_common_model('a2c2',1),'fwd_model');
0205    fmdl.mat_idx_to_electrode.nodes_electrode = nodes_electrode;
0206    fmdl.electrode(1).nodes = 26:41;
0207    fmdl.mat_idx{1} = [1:4];
0208    fmdl= mat_idx_to_electrode(fmdl, {1});
0209    unit_test_cmp('a2c2-01',num_elems(fmdl),64-4);
0210    if nodes_electrode;
0211       unit_test_cmp('a2c2-02',length(fmdl.electrode(2).nodes),4);
0212    else
0213       unit_test_cmp('a2c2-02a',length(fmdl.electrode(2).nodes),0);
0214       unit_test_cmp('a2c2-02b',size(fmdl.electrode(2).faces),[4,2]);
0215    end
0216 
0217    fmdl.stimulation = stim_meas_list([1,2,1,2]);
0218    img = mk_image(fmdl,1);
0219    img.fwd_solve.get_all_meas = true;
0220    vh = fwd_solve(img);
0221 
0222    img = rmfield(img,'elem_data');
0223    img.node_data = vh.volt;
0224    
0225    img.fwd_model = rmfield(img.fwd_model,'electrode');
0226    show_fem(img,[0,0,2]);
0227 %  plot(sqrt(sum(fmdl.nodes.^2,2)),vh.volt,'*');
0228    unit_test_cmp('a2c2-03', std(vh.volt(13:24)),0,0.002);
0229    vd = mean(vh.volt(5:12)) - mean(vh.volt(13:24));
0230    unit_test_cmp('a2c2-03', vd,0.065251,1e-6);
0231 
0232 
0233 
0234 function do_unit_test_3d_netgen(node_elecs_flag)
0235    extra={'ball_inside','ball_surface', [ ...
0236           'solid ball_inside  = sphere(-0.4,0,0.5;0.05);' ...
0237           'solid ball_surface = sphere( 0.4,0,1.0;0.05) -maxh=.005;' ...
0238           ]};
0239    fmdl= ng_mk_cyl_models(1,[8,.5],[.1],extra); 
0240    fmdl.mat_idx_to_electrode.nodes_electrode = node_elecs_flag;
0241 
0242    fmd_= mat_idx_to_electrode(fmdl, {2,3});
0243    unit_test_cmp('cyl_models 01',num_elecs(fmd_),10);
0244 
0245    fmd_= mat_idx_to_electrode(fmdl, {2:3});
0246    unit_test_cmp('cyl_models 02',num_elecs(fmd_),9);
0247 
0248    fmd_= mat_idx_to_electrode(fmdl, {2});
0249    unit_test_cmp('cyl_models 03',num_elecs(fmd_),9);
0250 
0251    img = mk_image(fmd_,1);
0252    img.elem_data(fmd_.mat_idx{3})= 1.1;
0253    img.calc_colours.ref_level = 1;
0254    unit_test_cmp('cyl_models 04',num_elecs(img),9);
0255 
0256    img = mk_image(fmdl,1);
0257    img.elem_data(vertcat(fmdl.mat_idx{2:3}))= 1.1;
0258    img.fwd_model.mat_idx_to_electrode.nodes_electrode=node_elecs_flag;
0259    [img.fwd_model,rm_elems]= mat_idx_to_electrode( ...
0260         img.fwd_model, {2});
0261    img.elem_data(rm_elems) = [];
0262    unit_test_cmp('cyl_models 05',num_elecs(img),9);
0263    vol = get_elem_volume(img);
0264    vvs = sum(vol(find(img.elem_data == 1.1)));
0265    unit_test_cmp('cyl_models 05a',vvs,4/3*pi*.05^3,1e-5);
0266   
0267 
0268 
0269    img.calc_colours.ref_level = 1;
0270    show_fem(img); view(3,12);
0271 
0272 
0273 function do_unit_test_3d_netgen2(node_elecs_flag)
0274    shape_str = ['solid sqelec = orthobrick(-1,-1,-1;1,1,0); tlo sqelec;' ...
0275                 'solid mainobj= orthobrick(-5,-5,-5;5,5,0) and not sqelec;'];
0276    fmdl = ng_mk_gen_models(shape_str, [],[],'');
0277    fmdl.mat_idx_to_electrode.nodes_electrode = node_elecs_flag;
0278    fmdl = mat_idx_to_electrode(fmdl,{1});
0279    elimnodes = [2 35 36; 4 34 36; 6 31 35; 8 31 34; 31 34 35; 34 35 36];
0280    sd = intersect(elimnodes,sort(fmdl.boundary,2),'rows');
0281    unit_test_cmp('cut boundary', size(sd,1),0);
0282 
0283 function do_unit_test_contained_shape(node_elecs_flag)
0284     splus = [ ' -maxh=10.02; tlo shape;' ...
0285        'solid el1 = sphere(1,0,0;0.11); tlo el1;' ...
0286        'solid el2 = sphere(0,1,0;0.11); tlo el2;' ...
0287        'solid nots = not( shape or el1 or el2);' ...
0288        'solid ground_ = sphere(0,0,0;2);' ...
0289        'solid ground  = ground_ and nots; tlo ground;' ...
0290        'solid mainobj = sphere(0,0,0;3) and not ground_;'];
0291     shape = 'solid shape = sphere(0.5,0.0,.5;0.2)';
0292     shape = 'solid shape = orthobrick(0.4,-0.1,.4;0.6,0.1,.6)';
0293     shape_str = [shape, splus ];
0294     fmdl = ng_mk_gen_models(shape_str,[],[],'');
0295     fmdl.mat_idx_to_electrode.nodes_electrode = node_elecs_flag;
0296     fmdl = mat_idx_to_electrode(fmdl,{1,2,3,5});
0297     if node_elecs_flag
0298         sizes = [8,61,61,134];
0299         for i=1:4; 
0300            unit_test_cmp('Contained ', ...
0301               size(fmdl.electrode(i).nodes),[1,sizes(i)])
0302         end
0303     else
0304         sizes = [12,118,118,264];
0305         for i=1:4; 
0306            unit_test_cmp('Contained ', ...
0307               size(fmdl.electrode(i).faces),[sizes(i),3])
0308         end
0309     end