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

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