Description of point_in

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 6516 2022-12-30 20:53:41Z 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    TR = triangulation(fmdl.elems, fmdl.nodes);
0045    warning(s.state,'MATLAB:triangulation:PtsNotInTriWarnId')
0046    ID = pointLocation(TR, points);
0047    idx = ~isnan(ID);
0048    point2tet = builtin('sparse',find(idx),ID(idx),1,size(points,1),size(fmdl.elems,1));
0049 end
0050 
0051 function point2tet = do_point_in_tet_inequalities(fmdl, points, epsilon, exclude_nodes)
0052    % use the correct tolerance
0053    if isstruct(epsilon)
0054       epsilon = epsilon.tol_node2tet;
0055    end
0056    progress_msg('Tet inequalities');
0057    [A,b] = tet_to_inequal(fmdl.nodes,fmdl.elems);
0058    progress_msg(Inf);
0059    
0060    %split computation into managable chunks to fit in memory
0061    mem_req = size(points,1) * size(fmdl.elems,1) * 8; % in bytes
0062    mem_use = 2*(1024^3); % 2 GiB
0063    n_chunks = ceil(mem_req / mem_use);
0064    chunk_sz = ceil(size(points,1) / n_chunks);
0065    point2tet = logical(builtin('sparse',0, size(fmdl.elems,1)));
0066    chunk_end = 0;
0067    progress_msg('Point in tet',0,n_chunks);
0068    for c = 1:n_chunks
0069        progress_msg(c,n_chunks)
0070        chunk_start = chunk_end+1;
0071        chunk_end = min(chunk_start + chunk_sz, size(points,1));
0072        idx = chunk_start:chunk_end;
0073        if true
0074            A_= A(1:4:end,:)*points(idx,:)';
0075            b_= b(1:4:end);
0076            p2t = (bsxfun(@minus, A_,b_) <= epsilon)';
0077            %        progress_msg(.21);
0078            for i = 2:4
0079                % that was actually slower
0080                %good = find(any(p2t,2));
0081                %p2t(good,:) = p2t(good,:) & (bsxfun(@minus, A(i:4:end,:)*points(idx(good),:)',b(i:4:end)) <= epsilon)';
0082                A_= A(i:4:end,:)*points(idx,:)';
0083                b_= b(i:4:end);
0084                p2t = p2t & (bsxfun(@minus, A_, b_) <= epsilon)';
0085                %           progress_msg(.21 + (i-1)*.23);
0086            end
0087            point2tet = [point2tet; builtin('sparse',p2t)];
0088        else
0089            % slower...
0090            p2t = (bsxfun(@minus, A*points(idx,:)',b) <= epsilon)';
0091            point2tet = [point2tet; builtin('sparse',reshape(all(reshape(p2t',4,[])),[],length(idx))')];
0092        end
0093    end
0094    progress_msg(Inf);
0095    
0096    % exclude coinciding nodes
0097    %    ex= bsxfun(@eq,points(:,1),fmdl.nodes(:,1)') & ...
0098    %        bsxfun(@eq,points(:,2),fmdl.nodes(:,2)') & ...
0099    %        bsxfun(@eq,points(:,3),fmdl.nodes(:,3)');
0100    %    progress_msg(.94);
0101    %    point2tet(any(ex,2),:) = 0;
0102        
0103    if exclude_nodes    
0104        ex = ismember(points, fmdl.nodes, 'rows');
0105        point2tet(ex,:) = 0;
0106    end
0107 
0108 end
0109 
0110 function do_unit_test
0111    do_unit_test_3d
0112    do_unit_test_2d
0113 end
0114 
0115 function do_unit_test_2d
0116    fmdl = mk_common_model('a2c2',4);
0117    fmdl = fmdl.fwd_model;
0118    epsilon = 0.01;
0119    % NOTE: Triangularization gives different for this case
0120 %  [x,y] = meshgrid(-2:0.5:2,-2:0.5:2);
0121 %  vv=[41 42 50 51 49 40 33 32 31 43 59 39];
0122 %  hh=[ 4 26 29 30 31 32 34 35 36 51 55 59];
0123    [x,y] = meshgrid(-2:0.6:2,-2:0.6:2);
0124    vv=[26 32 25 33 31 19 24 18];
0125    hh=[ 9 11 16 30 43 46 59 63];
0126    in=point_in_tet(fmdl,[x(:),y(:)],epsilon);
0127    in_=sparse(vv,hh,1,length(x(:)),num_elems(fmdl));
0128    unit_test_cmp('2D test',in,in_);
0129 in_ - in
0130 end
0131 
0132 function do_unit_test_3d
0133    fmdl = mk_common_model('n3r2',[16,2]);
0134    fmdl = fmdl.fwd_model;
0135    epsilon = 0;
0136    ll = linspace(-1,1,8);
0137    [x,y,z] = ndgrid(ll*1.8,ll*1.8,1.1:2.1);
0138    for ca = 1:2; switch ca
0139       case 1; fmdl.point_in_tet.algorithm = 'do_point_in_tet_inequalities'; 
0140       case 2; fmdl.point_in_tet.algorithm = 'do_point_in_tet_triangulation'; 
0141       end
0142       str = ['3D test:', fmdl.point_in_tet.algorithm];
0143       in=point_in_tet(fmdl,[x(:),y(:),z(:)],epsilon);
0144       vv=[  38  45  44  35  27  20  21];
0145       hh=[ 280 305 316 341 352 377 388];
0146       in_=sparse(vv,hh,1,length(x(:)),400);
0147       unit_test_cmp(str,all(sum(in,2)<=1), true);
0148       unit_test_cmp(str,in(:,1:400),in_);
0149 %     disp( in(:,1:400) )
0150    end
0151 end
point_in_tet

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
