point_in_tet

POINT_IN_TET test for points contained in elements
 point2tet = point_in_tet(fmdl,points, epsilon, exclude_nodes)

   fmdl: tet mesh to test against
   points: [x,y,z] of points to test
   epsilon: typically eps
   exclude_nodes: nodes to exclude

   fmdl.point_in_tet.algorithm: specify approach
     e.g. fmdl.point_in_tet.algorithm = 'do_point_in_tet_inequalities'
   
 point2tet: sparse matrix Npoints x Ntets 

 (C) 2015 Bartlomiej Grychtol
 License: GPL version 2 or 3
 $Id: point_in_tet.m 6943 2024-06-18 16:47:33Z aadler $

0001 function point2tet = point_in_tet(fmdl,points, epsilon, exclude_nodes)
0002 %POINT_IN_TET test for points contained in elements
0003 % point2tet = point_in_tet(fmdl,points, epsilon, exclude_nodes)
0004 %
0005 %   fmdl: tet mesh to test against
0006 %   points: [x,y,z] of points to test
0007 %   epsilon: typically eps
0008 %   exclude_nodes: nodes to exclude
0009 %
0010 %   fmdl.point_in_tet.algorithm: specify approach
0011 %     e.g. fmdl.point_in_tet.algorithm = 'do_point_in_tet_inequalities'
0012 %
0013 % point2tet: sparse matrix Npoints x Ntets
0014 %
0015 % (C) 2015 Bartlomiej Grychtol
0016 % License: GPL version 2 or 3
0017 % $Id: point_in_tet.m 6943 2024-06-18 16:47:33Z aadler $
0018 
0019 if nargin>=1 && ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test(); return; end
0020 
0021 if nargin < 4
0022     exclude_nodes = false;
0023 end
0024 
0025 ver = eidors_obj('interpreter_version');
0026 try 
0027    calc_fn= str2func(fmdl.point_in_tet.algorithm);
0028 catch % otherwise use triangulation unless octave
0029 if ver.isoctave 
0030    calc_fn= @do_point_in_tet_inequalities;
0031 else        
0032    calc_fn= @do_point_in_tet_triangulation;
0033 end
0034 end
0035 
0036 copt.fstr = 'point_in_tet';
0037 copt.cache_obj = {fmdl.nodes, fmdl.elems, points, epsilon, exclude_nodes, calc_fn};
0038 point2tet = eidors_cache(calc_fn,{fmdl, points, epsilon, exclude_nodes}, copt);
0039 
0040 end
0041 
0042 function point2tet = do_point_in_tet_triangulation(fmdl, points, epsilon, exclude_nodes)
0043    s = warning('off','MATLAB:triangulation:PtsNotInTriWarnId');
0044    % Triangularization requires doubles!
0045    % STUPID MATLAB puts isdouble() in a toolbox
0046    if ~isa(fmdl.elems,'double'); fmdl.elems=double(fmdl.elems); end
0047    TR = triangulation(fmdl.elems, fmdl.nodes);
0048    warning(s.state,'MATLAB:triangulation:PtsNotInTriWarnId')
0049    ID = pointLocation(TR, points);
0050    idx = ~isnan(ID);
0051    point2tet = sparse(find(idx),ID(idx),true,size(points,1),size(fmdl.elems,1));
0052 end
0053 
0054 function point2tet = do_point_in_tet_inequalities(fmdl, points, epsilon, exclude_nodes)
0055    % use the correct tolerance
0056    if isstruct(epsilon)
0057       epsilon = epsilon.tol_node2tet;
0058    end
0059    progress_msg('Tet inequalities');
0060    [A,b] = tet_to_inequal(fmdl.nodes,fmdl.elems);
0061    progress_msg(Inf);
0062    
0063    %split computation into managable chunks to fit in memory
0064    mem_req = size(points,1) * size(fmdl.elems,1) * 8; % in bytes
0065    mem_use = 2*(1024^3); % 2 GiB
0066    n_chunks = ceil(mem_req / mem_use);
0067    chunk_sz = ceil(size(points,1) / n_chunks);
0068    point2tet = logical(sparse(0, size(fmdl.elems,1)));
0069    chunk_end = 0;
0070    progress_msg('Point in tet',0,n_chunks);
0071    for c = 1:n_chunks
0072        progress_msg(c,n_chunks)
0073        chunk_start = chunk_end+1;
0074        chunk_end = min(chunk_start + chunk_sz, size(points,1));
0075        idx = chunk_start:chunk_end;
0076        if true
0077            A_= A(1:4:end,:)*points(idx,:)';
0078            b_= b(1:4:end);
0079            p2t = (bsxfun(@minus, A_,b_) <= epsilon)';
0080            %        progress_msg(.21);
0081            for i = 2:4
0082                % that was actually slower
0083                %good = find(any(p2t,2));
0084                %p2t(good,:) = p2t(good,:) & (bsxfun(@minus, A(i:4:end,:)*points(idx(good),:)',b(i:4:end)) <= epsilon)';
0085                A_= A(i:4:end,:)*points(idx,:)';
0086                b_= b(i:4:end);
0087                p2t = p2t & (bsxfun(@minus, A_, b_) <= epsilon)';
0088                %           progress_msg(.21 + (i-1)*.23);
0089            end
0090            point2tet = [point2tet; sparse(p2t)];
0091        else
0092            % slower...
0093            p2t = (bsxfun(@minus, A*points(idx,:)',b) <= epsilon)';
0094            point2tet = [point2tet; builtin('sparse',reshape(all(reshape(p2t',4,[])),[],length(idx))')];
0095        end
0096    end
0097    progress_msg(Inf);
0098    
0099    % exclude coinciding nodes
0100    %    ex= bsxfun(@eq,points(:,1),fmdl.nodes(:,1)') & ...
0101    %        bsxfun(@eq,points(:,2),fmdl.nodes(:,2)') & ...
0102    %        bsxfun(@eq,points(:,3),fmdl.nodes(:,3)');
0103    %    progress_msg(.94);
0104    %    point2tet(any(ex,2),:) = 0;
0105        
0106    if exclude_nodes    
0107        ex = ismember(points, fmdl.nodes, 'rows');
0108        point2tet(ex,:) = 0;
0109    end
0110 
0111 end
0112 
0113 function do_unit_test
0114    do_unit_test_3d
0115    do_unit_test_2d
0116 end
0117 
0118 function do_unit_test_2d
0119    fmdl = mk_common_model('a2c2',4);
0120    fmdl = fmdl.fwd_model;
0121    epsilon = 0.01;
0122    % NOTE: Triangularization gives different for this case
0123 %  [x,y] = meshgrid(-2:0.5:2,-2:0.5:2);
0124 %  vv=[41 42 50 51 49 40 33 32 31 43 59 39];
0125 %  hh=[ 4 26 29 30 31 32 34 35 36 51 55 59];
0126    [x,y] = meshgrid(-2:0.6:2,-2:0.6:2);
0127    vv=[26 32 25 33 31 19 24 18];
0128    hh=[ 9 11 16 30 43 46 59 63];
0129    in=point_in_tet(fmdl,[x(:),y(:)],epsilon);
0130    in_=sparse(vv,hh,1,length(x(:)),num_elems(fmdl));
0131    unit_test_cmp('2D test',in,in_);
0132 % in_ - in
0133 end
0134 
0135 function do_unit_test_3d
0136    fmdl = mk_common_model('n3r2',[16,2]);
0137    fmdl = fmdl.fwd_model;
0138    epsilon = 0;
0139    ll = linspace(-1,1,8);
0140    [x,y,z] = ndgrid(ll*1.8,ll*1.8,1.1:2.1);
0141    for ca = 1:2; switch ca
0142       case 1; fmdl.point_in_tet.algorithm = 'do_point_in_tet_inequalities'; 
0143       case 2; fmdl.point_in_tet.algorithm = 'do_point_in_tet_triangulation'; 
0144       end
0145       str = ['3D test:', fmdl.point_in_tet.algorithm];
0146       in=point_in_tet(fmdl,[x(:),y(:),z(:)],epsilon);
0147       vv=[  38  45  44  35  27  20  21];
0148       hh=[ 280 305 316 341 352 377 388];
0149       in_=sparse(vv,hh,1,length(x(:)),400);
0150       unit_test_cmp(str,all(sum(in,2)<=1), true);
0151       unit_test_cmp(str,in(:,1:400),in_);
0152 %     disp( in(:,1:400) )
0153    end
0154 end