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