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 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
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
0061 mem_req = size(points,1) * size(fmdl.elems,1) * 8;
0062 mem_use = 2*(1024^3);
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
0078 for i = 2:4
0079
0080
0081
0082 A_= A(i:4:end,:)*points(idx,:)';
0083 b_= b(i:4:end);
0084 p2t = p2t & (bsxfun(@minus, A_, b_) <= epsilon)';
0085
0086 end
0087 point2tet = [point2tet; builtin('sparse',p2t)];
0088 else
0089
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
0097
0098
0099
0100
0101
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
0120
0121
0122
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
0150 end
0151 end