find_boundary

 [srf, idx] = find_boundary(simp);

Caclulates the boundary faces of a given 3D volume.
Usefull in electrode assignment.

srf  =  array of elements on each boundary simplex
        boundary simplices are of 1 lower dimention than simp
idx  =  index of simplex to which each boundary belongs
simp = The simplices matrix

 Also:
  fmdl.boundary = find_boundary(fmdl_old)

0001 function [srf, idx] = find_boundary(simp);
0002 % [srf, idx] = find_boundary(simp);
0003 %
0004 %Caclulates the boundary faces of a given 3D volume.
0005 %Usefull in electrode assignment.
0006 %
0007 %srf  =  array of elements on each boundary simplex
0008 %        boundary simplices are of 1 lower dimention than simp
0009 %idx  =  index of simplex to which each boundary belongs
0010 %simp = The simplices matrix
0011 %
0012 % Also:
0013 %  fmdl.boundary = find_boundary(fmdl_old)
0014 
0015 % $Id: find_boundary.m 6354 2022-05-01 16:15:26Z aadler $
0016 
0017 if ischar(simp) && strcmp(simp,'UNIT_TEST'); do_unit_test; return; end
0018 if isstruct(simp) && strcmp(simp.type,'fwd_model'); simp= simp.elems; end
0019 
0020 wew = size(simp,2) - 1;
0021 
0022 if wew==3 || wew==2
0023    [srf,idx]= find_2or3d_boundary(simp,wew);
0024 elseif wew == 1
0025    [srf,idx]= find_1d_boundary(simp);
0026 else
0027    eidors_msg('find_boundary: WARNING: not 1D, 2D or 3D simplices',1);
0028    srf=[]; return;
0029 end
0030 
0031 % sort surfaces. If there is more than one, its not on the boundary
0032 function [srf,idx]= find_2or3d_boundary(simp,dim);
0033    if size(simp,1) < 4e9 % max of uint32
0034       % convert to integer to make sort faster
0035       simp = uint32( simp );
0036    end
0037    localface = nchoosek(1:dim+1,dim);
0038    srf_local= simp(:,localface');
0039    srf_local= reshape( srf_local', dim, []); % D x 3E
0040    srf_local= sort(srf_local)'; % Sort each row
0041    [sort_srl,sort_idx] = sortrows( srf_local );
0042 
0043    % Fine the ones that are the same
0044    first_ones =  sort_srl(1:end-1,:);
0045    next_ones  =  sort_srl(2:end,:);
0046    same_srl = find( all( first_ones == next_ones, 2) );
0047 
0048    % Assume they're all different. then find the same ones
0049    diff_srl = logical(ones(size(srf_local,1),1));
0050    diff_srl(same_srl) = 0;
0051    diff_srl(same_srl+1) = 0;
0052 
0053    srf= sort_srl( diff_srl,: );
0054    idx= sort_idx( diff_srl);
0055    idx= ceil(idx/(dim+1));
0056 
0057 function [srf,idx]= find_1d_boundary(simp);
0058    if size(simp,1) < 4e9 % max of uint32
0059       % convert to integer to make sort faster
0060       simp = uint32( simp );
0061    end
0062    % we expect two nodes as a result
0063    idx = find(isunique(simp(:)) == 1);
0064    srf = simp(idx);
0065    idx = rem(idx-1,size(simp,1))+1;
0066 
0067 function x = isunique(a);
0068    u=unique(a);
0069    n=histc(a,u);
0070    x=ismember(a,u(n==1));
0071 
0072 function do_unit_test
0073 
0074 %2D Test:
0075 mdl = mk_common_model('c2c',16);
0076 bdy = find_boundary(mdl.fwd_model.elems);
0077 bdy = sort_boundary(bdy);
0078 bdyc= sort_boundary(mdl.fwd_model.boundary);
0079 
0080 unit_test_cmp('2D test', bdy, bdyc);
0081 
0082 %3D Test:
0083 mdl = mk_common_model('n3r2',[16,2]);
0084 bdy = find_boundary(mdl.fwd_model.elems);
0085 bdy = sort_boundary(bdy);
0086 bdyc= sort_boundary(mdl.fwd_model.boundary);
0087 
0088 unit_test_cmp( '3D test n3r2', bdy,bdyc);
0089 
0090 %3D Test:
0091 mdl = mk_common_model('a3cr',16);
0092 bdy = find_boundary(mdl.fwd_model.elems);
0093 bdy = sort_boundary(bdy);
0094 bdyc= sort_boundary(mdl.fwd_model.boundary);
0095 
0096 unit_test_cmp('3D test a3c2', bdy, bdyc);
0097 
0098 %3D Test:
0099 mdl = mk_common_model('b3cr',16);
0100 bdy = find_boundary(mdl.fwd_model.elems);
0101 bdy = sort_boundary(bdy);
0102 bdyc= sort_boundary(mdl.fwd_model.boundary);
0103 
0104 unit_test_cmp('3D test b3c2', bdy, bdyc);
0105 
0106 simp = [  10 190; ...
0107          182 183; ...
0108          183 184; ...
0109          184 185; ...
0110           11 182; ...
0111          185 186; ...
0112          187 186; ...
0113          187 188; ...
0114          188 189; ...
0115          189 190];
0116 [bdy, idx] = find_boundary(simp);
0117 unit_test_cmp('1D bdy', bdy,[10;11]);
0118 unit_test_cmp('1D bdy', idx,[1;5]);
0119 
0120 function bdy= sort_boundary(bdy)
0121    bdy = sort(bdy,2);
0122    bdy = sortrows(bdy);
0123