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.html 2819 2011-09-07 16:43:11Z aadler $
0013 
0014 if isstr(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 else
0022    eidors_msg('find_boundary: WARNING: not 2D or 3D simplices',1);
0023    srf=[]; return;
0024 end
0025 
0026 % sort surfaces. If there is more than one, its not on the boundary
0027 function [srf,idx]= find_2or3d_boundary(simp,dim);
0028    if size(simp,1) < 4e9 % max of uint32
0029       % convert to integer to make sort faster
0030       simp = uint32( simp );
0031    end
0032    localface = nchoosek(1:dim+1,dim);
0033    srf_local= simp(:,localface');
0034    srf_local= reshape( srf_local', dim, []); % D x 3E
0035    srf_local= sort(srf_local)'; % Sort each row
0036    [sort_srl,sort_idx] = sortrows( srf_local );
0037 
0038    % Fine the ones that are the same
0039    first_ones =  sort_srl(1:end-1,:);
0040    next_ones  =  sort_srl(2:end,:);
0041    same_srl = find( all( first_ones == next_ones, 2) );
0042 
0043    % Assume they're all different. then find the same ones
0044    diff_srl = logical(ones(size(srf_local,1),1));
0045    diff_srl(same_srl) = 0;
0046    diff_srl(same_srl+1) = 0;
0047 
0048    srf= sort_srl( diff_srl,: );
0049    idx= sort_idx( diff_srl);
0050    idx= ceil(idx/(dim+1));
0051 
0052 function do_unit_test
0053 ok=1;
0054 
0055 %2D Test:
0056 mdl = mk_common_model('c2c',16);
0057 bdy = find_boundary(mdl.fwd_model.elems);
0058 bdy = sort_boundary(bdy);
0059 bdyc= sort_boundary(mdl.fwd_model.boundary);
0060 
0061 ok= match(bdy,bdyc,ok,'2D test');
0062 
0063 %3D Test:
0064 mdl = mk_common_model('n3r2',16);
0065 bdy = find_boundary(mdl.fwd_model.elems);
0066 bdy = sort_boundary(bdy);
0067 bdyc= sort_boundary(mdl.fwd_model.boundary);
0068 
0069 ok= match(bdy,bdyc,ok,'3D test n3r2');
0070 
0071 %3D Test:
0072 mdl = mk_common_model('a3cr',16);
0073 bdy = find_boundary(mdl.fwd_model.elems);
0074 bdy = sort_boundary(bdy);
0075 bdyc= sort_boundary(mdl.fwd_model.boundary);
0076 
0077 ok= match(bdy,bdyc,ok,'3D test a3c2');
0078 
0079 %3D Test:
0080 mdl = mk_common_model('b3cr',16);
0081 bdy = find_boundary(mdl.fwd_model.elems);
0082 bdy = sort_boundary(bdy);
0083 bdyc= sort_boundary(mdl.fwd_model.boundary);
0084 
0085 ok= match(bdy,bdyc,ok,'3D test b3c2');
0086 
0087 function bdy= sort_boundary(bdy)
0088    bdy = sort(bdy,2);
0089    bdy = sortrows(bdy);
0090 
0091 function ok= match( pat1, pat2, ok, descr)
0092     ok =  all(pat1(:) == pat2(:));
0093     fprintf('find_bounday_test: ok=%d\n',ok);
0094 
0095 % function ok= match( pat1, pat2, ok, descr)
0096 %    if ~all(pat1(:) == pat2(:))
0097 %       ok=0;
0098 %       eidors_msg('find_bounday_test: fail %s',descr,1);
0099 %    end
0100 
0101