0001 function [pts,FE2CE,FE2pts,CE2pts] = find_edge2edge_intersections(FE,FN,CE,CN, epsilon)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031 if size(CE,1) > size(FE,1)
0032
0033 [pts,FE2CE,CE2pts,FE2pts] = find_edge2edge_intersections(CE,CN,FE,FN, epsilon);
0034 FE2CE = FE2CE';
0035 return
0036 end
0037
0038 epsilon = epsilon.^2;
0039
0040
0041 try
0042
0043
0044 [pts,FE2CE,FE2pts,CE2pts] = edge2edge_intersections_serial(FE,FN,CE,CN,epsilon);
0045 catch err
0046 if strcmp(err.identifier, 'MATLAB:nomem')
0047
0048 [pts,FE2CE,FE2pts,CE2pts] = edge2edge_intersections_wrapper(FE,FN,CE,CN,epsilon);
0049 else
0050 rethrow(err);
0051 end
0052 end
0053
0054
0055
0056
0057 function [intpts, FE2CE, FE2pts, CE2pts] = edge2edge_intersections_wrapper(FE,FN,CE,CN, epsilon)
0058
0059
0060 sz = size(FE,1) * size(CE,1) * 8;
0061 desired_mem = 2*(1024^3 );
0062
0063 n_chunks = ceil(10*sz / desired_mem);
0064
0065 len_chnk = ceil(size(CE,1) / n_chunks);
0066 intpts = [];
0067 FE2CE = sparse(0,0);
0068 FE2pts = sparse(0,0);
0069 CE2pts = sparse(0,0);
0070 for c = 1:n_chunks
0071 eidors_msg('@@@: chunk %d of %d',c,n_chunks,2);
0072 start = 1 + (c-1)*len_chnk;
0073 stop = min(1 + c*len_chnk, size(CE,1));
0074 rng = start:stop;
0075 [ip, te2ve, te2ip, ve2ip] = edge2edge_intersections(FE,FN,CE(rng,:),CN, epsilon);
0076 len = size(intpts,1);
0077 intpts = [intpts; ip];
0078 idx = te2ve>0;
0079 te2ve(idx) = len + te2ve(idx);
0080 FE2CE = [FE2CE te2ve];
0081 FE2pts = [FE2pts te2ip];
0082 if size(ve2ip,1) < size(CE2pts,1)
0083 ve2ip(size(CE2pts,1),1) = 0;
0084 end
0085 CE2pts = [CE2pts ve2ip];
0086 end
0087
0088
0089
0090
0091 function [intpts, FE2CE, FE2pts, CE2pts] = edge2edge_intersections_serial(FE,FN,CE,CN, epsilon)
0092
0093
0094
0095
0096
0097 progress_msg(0);
0098
0099 P1 = FN(FE(:,1),:);
0100 P2 = FN(FE(:,2),:);
0101 P3 = CN(CE(:,1),:);
0102 P4 = CN(CE(:,2),:);
0103
0104 C_bb = zeros(size(P3,1),6);
0105 C_bb(:,[1 3 5]) = min(P3,P4);
0106 C_bb(:,[2 4 6]) = max(P3,P4);
0107
0108 F_bb = zeros(size(P1,1),6);
0109 F_bb(:,[1 3 5]) = min(P1,P2);
0110 F_bb(:,[2 4 6]) = max(P1,P2);
0111
0112
0113 excl = bsxfun(@gt, F_bb(:,1), C_bb(:,2)') ...
0114 | bsxfun(@lt, F_bb(:,2), C_bb(:,1)') ...
0115 | bsxfun(@gt, F_bb(:,3), C_bb(:,4)') ...
0116 | bsxfun(@lt, F_bb(:,4), C_bb(:,3)') ...
0117 | bsxfun(@gt, F_bb(:,5), C_bb(:,6)') ...
0118 | bsxfun(@lt, F_bb(:,6), C_bb(:,5)');
0119 excl = ~excl;
0120
0121
0122
0123
0124
0125
0126 p43 = P4 - P3;
0127 p21 = P2 - P1;
0128
0129 d4343 = sum(p43.^2,2);
0130 d2121 = sum(p21.^2,2);
0131
0132
0133
0134 IN = false(size(P1,1),size(P3,1));
0135
0136 intpts = [];
0137
0138 todo = size(CE,1);
0139 id = 0;
0140 step = ceil(todo/100);
0141
0142 pa = zeros(size(P1));
0143
0144 for v = 1:size(CE,1)
0145 id = id+1;
0146 if mod(id,step)==0, progress_msg(id/todo); end
0147
0148 fid = excl(:,v);
0149 if ~any(fid), continue, end;
0150
0151 d1343 = 0; d4321 = 0; d1321 = 0;
0152 for i = 1:3
0153 p13 = P1(fid,i) - P3(v,i);
0154 d1343 = d1343 + p13*p43(v,i);
0155 d4321 = d4321 + p21(fid,i)*p43(v,i);
0156 d1321 = d1321 + p13 .* p21(fid,i);
0157 end
0158
0159
0160
0161 denom = d2121(fid) * d4343(v) - d4321.^2;
0162
0163
0164 numer = d1343 .* d4321 - d1321 * d4343(v);
0165
0166 mua = numer./denom;
0167 mub = (d1343 + d4321.*mua) / d4343(v);
0168
0169
0170 D = 0;
0171 for i = 1:3
0172 pa(fid,i) = P1(fid,i) + p21(fid,i).*mua;
0173 pb = P3(v,i) + p43(v,i) * mub;
0174 D = D + (pa(fid,i) - pb).^2;
0175 end
0176
0177
0178 IN(fid,v) = ((mua>0) + (mua<1) + (mub>0) + (mub<1) + (D<=epsilon)) == 5;
0179 if any(IN(fid,v))
0180 intpts = [intpts; pa(IN(:,v),:)];
0181 end
0182 end
0183 [T, V] = find(IN);
0184 I = (1:size(intpts,1))';
0185 FE2CE = sparse(size(P1,1),size(P3,1));
0186 FE2CE(sub2ind(size(FE2CE),T,V)) = I;
0187 FE2pts = sparse(T,I,ones(size(I)),size(P1,1),size(I,1));
0188 CE2pts = sparse(V,I,ones(size(I)),size(P3,1),size(I,1));
0189
0190 if eidors_msg('log_level')>1
0191 progress_msg(Inf);
0192 end
0193
0194
0195
0196
0197 function [intpts, FE2CE, FE2pts, CE2pts] = edge2edge_intersections(FE,FN,CE,CN, epsilon)
0198
0199
0200
0201
0202
0203
0204
0205 P1 = FN(FE(:,1),:);
0206 P2 = FN(FE(:,2),:);
0207 P3 = CN(CE(:,1),:);
0208 P4 = CN(CE(:,2),:);
0209
0210
0211
0212
0213
0214
0215 p43 = P4 - P3;
0216 p21 = P2 - P1;
0217
0218
0219 JNK = bsxfun(@or,sparse(all(abs(p21)<eps,2)),sparse(all(abs(p43)<eps,2)'));
0220
0221 f = {'x','y','z'};
0222 d1343 = 0; d4321 = 0; d1321 = 0;
0223
0224 for i = 1:3;
0225 p13 = bsxfun(@minus, P1(:,i), P3(:,i)');
0226 d1343 = d1343 + bsxfun(@times,p13,p43(:,i)');
0227 d4321 = d4321 + bsxfun(@times,p21(:,i),p43(:,i)');
0228 d1321 = d1321 + bsxfun(@times,p13,p21(:,i));
0229 end
0230
0231 clear p13
0232 d4343 = sum(p43.^2,2);
0233 d2121 = sum(p21.^2,2);
0234
0235 denom = bsxfun(@times,d2121,d4343') - d4321.^2;
0236 JNK = JNK | abs(denom)<eps;
0237 numer = d1343 .* d4321 - bsxfun(@times,d1321,d4343');
0238
0239 mua = numer./denom;
0240 mub = bsxfun(@rdivide, d1343 + d4321.*mua , d4343');
0241
0242 clear d1343 d4321 d1321 d4343 d2121
0243
0244 pa = struct; pb = struct;
0245 D = 0;
0246 for i = 1:3
0247 pa.(f{i}) = bsxfun(@plus,P1(:,i),bsxfun(@times,p21(:,i),mua));
0248 pb = bsxfun(@plus,P3(:,i)',bsxfun(@times,p43(:,i)',mub));
0249 D = D + (pa.(f{i}) - pb).^2;
0250 end
0251
0252 D(JNK) = Inf;
0253
0254
0255
0256 IN = mua>0 & mua <1 & mub>0 & mub<1 & D<=epsilon;
0257 nPts = nnz(IN);
0258 intpts = zeros(nPts,3);
0259 for i = 1:3
0260 intpts(:,i) = pa.(f{i})(IN);
0261 end
0262
0263 [T, V] = find(IN);
0264 I = (1:nPts)';
0265 FE2CE = sparse(size(P1,1),size(P3,1));
0266 FE2CE(IN) = I;
0267
0268 FE2pts = sparse(T,I,ones(size(I)),size(P1,1),size(I,1));
0269 CE2pts = sparse(V,I,ones(size(I)),size(P3,1),size(I,1));
0270