0001 function point2tet = point_in_tet(fmdl,points, epsilon, exclude_nodes)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
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
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
0045
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
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
0064 mem_req = size(points,1) * size(fmdl.elems,1) * 8;
0065 mem_use = 2*(1024^3);
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
0081 for i = 2:4
0082
0083
0084
0085 A_= A(i:4:end,:)*points(idx,:)';
0086 b_= b(i:4:end);
0087 p2t = p2t & (bsxfun(@minus, A_, b_) <= epsilon)';
0088
0089 end
0090 point2tet = [point2tet; sparse(p2t)];
0091 else
0092
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
0100
0101
0102
0103
0104
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
0123
0124
0125
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
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
0153 end
0154 end