0001 function c2f = mk_tri2tet_c2f(fmdl,rmdl, opt)
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
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046 if nargin == 0 || (ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'))
0047 do_unit_test;
0048 return;
0049 end
0050 if nargin < 3
0051 opt = struct();
0052 end
0053
0054 f_elems = size(fmdl.elems,1);
0055 r_elems = size(rmdl.elems,1);
0056 c2f = sparse(f_elems,r_elems);
0057
0058 if size(rmdl.nodes,2) == 2
0059 rmdl.nodes(:,3) = max(fmdl.nodes(:,3))/2 + min(fmdl.nodes(:,3))/2;
0060 end
0061
0062 [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl, opt);
0063
0064 if ~(any(fmdl_idx) && any(rmdl_idx))
0065 eidors_msg('@@: models do not overlap, returning all-zeros');
0066 return
0067 end
0068
0069 opt = parse_opts(fmdl,rmdl, opt);
0070
0071 [fmdl,rmdl, opt] = center_scale_models(fmdl,rmdl, opt);
0072
0073 copt.fstr = 'mk_tri2tet_c2f';
0074 c2f(fmdl_idx,rmdl_idx) = eidors_cache(@do_mk_tri2tet_c2f,{fmdl,rmdl,opt},copt);
0075
0076
0077 end
0078
0079 function c2f = do_mk_tri2tet_c2f(fmdl,rmdl,opt)
0080
0081 DEBUG = eidors_debug('query','mk_tri2tet_c2f');
0082
0083 r_elems = size(rmdl.elems,1);
0084 f_elems = size(fmdl.elems,1);
0085
0086 c2f = sparse(f_elems, r_elems);
0087
0088 progress_msg('Prepare tet model...');
0089 fmdl = prepare_tet_mdl(fmdl);
0090 progress_msg(Inf);
0091
0092 progress_msg('Prepare tri model...');
0093 rmdl = prepare_tri_mdl(rmdl);
0094 progress_msg(Inf);
0095
0096 pmopt.final_msg = 'none';
0097 if ~isinf(opt.z_depth)
0098
0099 progress_msg('Find top cap nodes in tets...',-1,pmopt);
0100 [top_node2tet, top_nodes, top_nodes_above, top_nodes_below] = ...
0101 get_cap_nodes_in_tets(fmdl,rmdl,opt.top,opt.tol_node2tet);
0102 progress_msg(sprintf('Found %d', nnz(top_node2tet)), Inf);
0103
0104 progress_msg('Find bottom cap nodes in tets...',-1,pmopt);
0105 [bot_node2tet, bot_nodes, bot_nodes_above, bot_nodes_below] = ...
0106 get_cap_nodes_in_tets(fmdl,rmdl,opt.bot,opt.tol_node2tet);
0107 progress_msg(sprintf('Found %d', nnz(bot_node2tet)), Inf);
0108
0109
0110 progress_msg('Find tet_face2tri_edge intersections (top) ...');
0111 [intpts4, top_tet_face2tri_edge,top_tet_face2intpt4,top_tri_edge2intpt4] = ...
0112 get_cap_tet_face2tri_edge_intersections( fmdl,rmdl,opt.top, ...
0113 top_nodes_above,top_nodes_below, opt.tol_edge2tri);
0114 progress_msg(sprintf('Found %d', size(intpts4,1)),Inf);
0115
0116 progress_msg('Find tet_face2tri_edge intersections (bot) ...');
0117 [intpts5, bot_tet_face2tri_edge,bot_tet_face2intpt5,bot_tri_edge2intpt5] = ...
0118 get_cap_tet_face2tri_edge_intersections( fmdl,rmdl,opt.bot, ...
0119 bot_nodes_above,bot_nodes_below, opt.tol_edge2tri);
0120 progress_msg(sprintf('Found %d', size(intpts5,1)),Inf);
0121
0122
0123 progress_msg('Find tet_edge2tri_face intersections (top) ...');
0124 [intpts6, top_tet_edge2tri_face, top_tet_edge2intpt6, tri2intpt6] = ...
0125 get_cap_tet_edge2tri_face_intersections(fmdl,rmdl,opt.top, ...
0126 top_nodes_above, top_nodes_below, opt.tol_edge2tri);
0127 progress_msg(sprintf('Found %d', size(intpts6,1)),Inf);
0128
0129
0130 progress_msg('Find tet_edge2tri_face intersections (bot) ...');
0131 [intpts7, bot_tet_edge2tri_face, bot_tet_edge2intpt7, tri2intpt7] = ...
0132 get_cap_tet_edge2tri_face_intersections(fmdl,rmdl,opt.bot, ...
0133 bot_nodes_above, bot_nodes_below, opt.tol_edge2tri);
0134 progress_msg(sprintf('Found %d', size(intpts6,1)),Inf);
0135
0136
0137
0138 progress_msg('Find tet_edge2tri_edge intersections (top) ...',-1,pmopt);
0139 edgeidx = any(top_nodes_above(fmdl.edges),2) & any(top_nodes_below(fmdl.edges),2);
0140 nodes = rmdl.nodes;
0141 nodes(:,3) = opt.top;
0142 top_tet_edge2tri_edge = sparse(size(fmdl.edges,1),size(rmdl.edges,1));
0143 [intpts8, top_tet_edge2tri_edge(edgeidx,:), top_tet_edge2intpt8(edgeidx,:), top_tri_edge2intpt8] = ...
0144 find_edge2edge_intersections(fmdl.edges(edgeidx,:),fmdl.nodes, ...
0145 rmdl.edges,nodes, opt.tol_edge2edge);
0146 if size(top_tet_edge2intpt8,1)~=size(fmdl.edges,1)
0147 if ~isempty(intpts8)
0148 top_tet_edge2intpt8(size(fmdl.edges,1),1) = 0;
0149 else
0150 top_tet_edge2intpt8 = sparse(size(fmdl.edges,1),0);
0151 end
0152 end
0153 progress_msg(sprintf('Found %d', size(intpts8,1)),Inf);
0154
0155
0156 progress_msg('Find tet_edge2tri_edge intersections (bot) ...',-1,pmopt);
0157 edgeidx = any(bot_nodes_above(fmdl.edges),2) & any(bot_nodes_below(fmdl.edges),2);
0158 nodes = rmdl.nodes;
0159 nodes(:,3) = opt.bot;
0160 bot_tet_edge2tri_edge = sparse(size(fmdl.edges,1),size(rmdl.edges,1));
0161 [intpts9, bot_tet_edge2tri_edge(edgeidx,:), ...
0162 bot_tet_edge2intpt9(edgeidx,:), bot_tri_edge2intpt9] = ...
0163 find_edge2edge_intersections(fmdl.edges(edgeidx,:),fmdl.nodes, ...
0164 rmdl.edges,nodes, opt.tol_edge2edge);
0165 if size(bot_tet_edge2intpt9,1)~=size(fmdl.edges,1)
0166 if ~isempty(intpts9)
0167 bot_tet_edge2intpt9(size(fmdl.edges,1),1) = 0;
0168 else
0169 bot_tet_edge2intpt9 = sparse(size(fmdl.edges,1),0);
0170 end
0171 end
0172 progress_msg(sprintf('Found %d', size(intpts9,1)),Inf);
0173
0174
0175 progress_msg('Find tris contained in tet...')
0176 tri_in_tet = rmdl.node2elem' * (bot_node2tet + top_node2tet) == 6;
0177 progress_msg(sprintf('Found %d',nnz(tri_in_tet)), Inf);
0178 end
0179
0180
0181
0182 progress_msg('Find tet_nodes in tri_elems...');
0183 in = (fmdl.nodes(:,3) >= opt.bot) & (fmdl.nodes(:,3) <= opt.top);
0184 tet_node2tri = spalloc(size(fmdl.nodes,1),size(rmdl.elems,1),nnz(in));
0185 tet_node2tri(in,:) = point_in_triangle(fmdl.nodes(in,1:2), ...
0186 rmdl.elems, ...
0187 rmdl.nodes(:,1:2), ...
0188 opt.tol_node2tri);
0189 progress_msg(sprintf('Found %d', nnz(tet_node2tri)), Inf);
0190
0191 n_nodes = size(rmdl.nodes,1);
0192 vert_edges(:,1) = 1:n_nodes;
0193 vert_edges(:,2) = vert_edges(:,1) + n_nodes;
0194
0195 progress_msg('Find tet_edge2vert_edge intersections...',-1,pmopt);
0196 if isinf(opt.z_depth)
0197 bot_nodes = rmdl.nodes; bot_nodes(:,3) = opt.bot - 1;
0198 top_nodes = rmdl.nodes; top_nodes(:,3) = opt.top + 1;
0199 end
0200 [intpts1,tet_edge2vert_edge,tet_edge2intpt1,vert_edge2intpt1] = ...
0201 find_edge2edge_intersections( fmdl.edges, fmdl.nodes,...
0202 vert_edges, [bot_nodes; top_nodes], ...
0203 opt.tol_edge2edge);
0204 progress_msg(sprintf('Found %d',size(intpts1,1)),Inf);
0205
0206
0207
0208 progress_msg('Find tet_edge2vert_face intersections...');
0209 if isinf(opt.z_depth)
0210 top = opt.top + 1;
0211 bot = opt.bot - 1;
0212 else
0213 top = opt.top;
0214 bot = opt.bot;
0215 end
0216 [intpts2,tet_edge2tri_edge,tet_edge2intpt2,tri_edge2intpt2] = ...
0217 find_edge2edge_intersections_2d( fmdl.edges, fmdl.nodes(:,1:3), ...
0218 rmdl.edges, rmdl.nodes(:,1:2), ...
0219 top, bot);
0220 progress_msg(sprintf('Found %d',size(intpts2,1)),Inf);
0221
0222
0223 progress_msg('Find tet_face2vert_edge intersections...');
0224 [intpts3, tet_face2vert_edge, tet_face2intpt3, vert_edge2intpt3] = ...
0225 get_tet_intersection_pts(fmdl,rmdl,top,bot,opt.tol_edge2tri);
0226 progress_msg(sprintf('Found %d',size(intpts3,1)),Inf);
0227
0228
0229
0230 progress_msg('Find tets contained in tri...');
0231 tet_in_tri = (double(fmdl.node2elem') * tet_node2tri) == 4;
0232 progress_msg(sprintf('Found %d',nnz(tet_in_tri)), Inf);
0233
0234
0235 progress_msg('Find total tri v. tet intersections...');
0236 tri2intTet = double(rmdl.edge2elem') * tet_edge2tri_edge' * fmdl.edge2elem ...
0237 |double(rmdl.node2elem') * (tet_face2vert_edge>0)' * fmdl.elem2face';
0238 if ~isinf(opt.z_depth)
0239 tri2intTet = tri2intTet ...
0240 | double(rmdl.edge2elem') * (top_tet_face2tri_edge + bot_tet_face2tri_edge)' * fmdl.elem2face' ...
0241 | (top_tet_edge2tri_face + bot_tet_edge2tri_face)' * fmdl.edge2elem;
0242 end
0243
0244 tri2intTet = tri2intTet & ~tet_in_tri';
0245 if ~isinf(opt.z_depth)
0246 tri2intTet = tri2intTet & ~tri_in_tet;
0247 end
0248 progress_msg(sprintf('Found %d',nnz(tri2intTet)), Inf);
0249
0250
0251 progress_msg('Calculate intersection volumes...');
0252
0253 tri2intpt1 = logical(rmdl.node2elem'*vert_edge2intpt1)';
0254 tet2intpt1 = logical(fmdl.edge2elem'* tet_edge2intpt1)';
0255
0256 tet2intpt2 = logical(fmdl.edge2elem'*tet_edge2intpt2)';
0257 tri2intpt2 = logical(rmdl.edge2elem'*tri_edge2intpt2)';
0258
0259 tet2intpt3 = logical(double(fmdl.elem2face)*tet_face2intpt3)';
0260 tri2intpt3 = logical(rmdl.node2elem'*vert_edge2intpt3)';
0261 if ~isinf(opt.z_depth)
0262 tet2intpt4 = logical(double(fmdl.elem2face)*top_tet_face2intpt4)';
0263 tri2intpt4 = logical(rmdl.edge2elem'*top_tri_edge2intpt4)';
0264
0265 tet2intpt5 = logical(double(fmdl.elem2face)*bot_tet_face2intpt5)';
0266 tri2intpt5 = logical(rmdl.edge2elem'*bot_tri_edge2intpt5)';
0267
0268 tet2intpt6 = logical(fmdl.edge2elem'*top_tet_edge2intpt6)';
0269 tri2intpt6 = tri2intpt6';
0270
0271 tet2intpt7 = logical(fmdl.edge2elem'*bot_tet_edge2intpt7)';
0272 tri2intpt7 = tri2intpt7';
0273
0274 tet2intpt8 = logical(fmdl.edge2elem'*top_tet_edge2intpt8)';
0275 tri2intpt8 = logical(rmdl.edge2elem'*top_tri_edge2intpt8)';
0276
0277 tet2intpt9 = logical(fmdl.edge2elem'*bot_tet_edge2intpt9)';
0278 tri2intpt9 = logical(rmdl.edge2elem'*bot_tri_edge2intpt9)';
0279 end
0280 tri_todo = find(sum(tri2intTet,2)>0);
0281 C = []; F = []; V = [];
0282
0283 id = 0; lvox = length(tri_todo);
0284 mint = ceil(lvox/100);
0285
0286 for v = tri_todo'
0287 id = id+1;
0288 if mod(id,mint)==0, progress_msg(id/lvox); end
0289 tet_todo = find(tri2intTet(v,:));
0290 common_intpts1 = bsxfun(@and,tri2intpt1(:,v), tet2intpt1(:,tet_todo));
0291 common_intpts2 = bsxfun(@and,tri2intpt2(:,v), tet2intpt2(:,tet_todo));
0292 common_intpts3 = bsxfun(@and,tri2intpt3(:,v), tet2intpt3(:,tet_todo));
0293 if ~isinf(opt.z_depth)
0294 common_intpts4 = bsxfun(@and,tri2intpt4(:,v), tet2intpt4(:,tet_todo));
0295 common_intpts5 = bsxfun(@and,tri2intpt5(:,v), tet2intpt5(:,tet_todo));
0296 common_intpts6 = bsxfun(@and,tri2intpt6(:,v), tet2intpt6(:,tet_todo));
0297 common_intpts7 = bsxfun(@and,tri2intpt7(:,v), tet2intpt7(:,tet_todo));
0298 common_intpts8 = bsxfun(@and,tri2intpt8(:,v), tet2intpt8(:,tet_todo));
0299 common_intpts9 = bsxfun(@and,tri2intpt9(:,v), tet2intpt9(:,tet_todo));
0300 top_nodes_tet = bsxfun(@and,top_node2tet(:,tet_todo), rmdl.node2elem(:,v));
0301 bot_nodes_tet = bsxfun(@and,bot_node2tet(:,tet_todo), rmdl.node2elem(:,v));
0302 end
0303 tet_nodes = bsxfun(@and,tet_node2tri(:,v), fmdl.node2elem(:,tet_todo));
0304
0305 C = [C; v*ones(numel(tet_todo),1)];
0306 F = [F; tet_todo'];
0307 last_v = numel(V);
0308 V = [V; zeros(numel(tet_todo),1)];
0309 for t = 1:numel(tet_todo)
0310 pts = [ intpts1(common_intpts1(:,t),:);
0311 intpts2(common_intpts2(:,t),:);
0312 intpts3(common_intpts3(:,t),:);
0313 fmdl.nodes(tet_nodes(:,t),:);];
0314 if ~isinf(opt.z_depth)
0315 pts = [ pts;
0316 intpts4(common_intpts4(:,t),:);
0317 intpts5(common_intpts5(:,t),:);
0318 intpts6(common_intpts6(:,t),:);
0319 intpts7(common_intpts7(:,t),:);
0320 intpts8(common_intpts8(:,t),:);
0321 intpts9(common_intpts9(:,t),:);
0322 top_nodes(top_nodes_tet(:,t),:);
0323 bot_nodes(bot_nodes_tet(:,t),:)];
0324 end
0325 last_v = last_v + 1;
0326
0327 if size(pts,1) < 4
0328
0329
0330 continue
0331 end
0332 if any(isnan(pts(:))), keyboard, end
0333
0334 try
0335
0336
0337 ctr = mean(pts);
0338 if any(isnan(ctr)), keyboard,end
0339 pts = bsxfun(@minus,pts,ctr);
0340 scale = max(abs(pts(:)));
0341 if scale == 0
0342 continue
0343 end
0344
0345 pts = pts ./ scale;
0346
0347
0348 [K, V(last_v)] = convhulln(pts,{'Qt Pp Qs'});
0349 V(last_v) = max(V(last_v),0);
0350 V(last_v) = V(last_v) * scale^3;
0351 catch err
0352 ok = false;
0353 switch err.identifier
0354 case {'MATLAB:qhullmx:DegenerateData', 'MATLAB:qhullmx:UndefinedError'}
0355 pts = bsxfun(@plus, pts .* scale, ctr);
0356 u = uniquetol(pts,1e-14,'ByRows',true,'DataScale', 1);
0357 ok = ok | size(u,1) < 4;
0358 if ~ok
0359
0360 u12 = uniquetol(pts(:,1:2),1e-14,'ByRows',true,'DataScale',1);
0361 cp = bsxfun(@minus, u12(2:end,1:2), u12(1,1:2));
0362 l = sqrt(sum(cp.^2,2));
0363 cp = bsxfun(@rdivide, cp, l);
0364
0365 cp = abs(cp);
0366 un = uniquetol(cp,1e-12,'ByRows',true,'DataScale',1);
0367 ok = ok | size(un,1) == 1;
0368 end
0369 if ~ok
0370
0371 ok = ok | all(abs(pts(:,3) - top) < eps);
0372 ok = ok | all(abs(pts(:,3) - bot) < eps);
0373 end
0374 end
0375 if ~ok
0376 if DEBUG || eidors_debug('query','mk_tri2tet_c2f:convhulln')
0377 debug_plot(fmdl,rmdl,v,tet_todo(t), bot, top, pts)
0378 keyboard
0379 else
0380 fprintf('\n');
0381 eidors_msg(['convhulln has thrown an error. ' ...
0382 'Enable eidors_debug on mk_tri2tet_c2f and re-run to see a debug plot'],0);
0383 rethrow(err);
0384 end
0385 end
0386 end
0387 end
0388 end
0389 progress_msg(Inf);
0390 c2f = sparse(F,C,V,size(fmdl.elems,1),size(rmdl.elems,1));
0391
0392
0393 if ~isinf(opt.z_depth)
0394 try rmdl = rmfield(rmdl,'coarse2fine'); end
0395 c2f = c2f + bsxfun(@times, sparse(tri_in_tet), opt.height * get_elem_volume(rmdl))';
0396 end
0397
0398 vol = get_elem_volume(fmdl);
0399 c2f = bsxfun(@rdivide,c2f,vol);
0400
0401
0402
0403 c2f = c2f + tet_in_tri;
0404
0405 end
0406
0407 function [intpts, tri2edge, tri2pts, edge2pts] = get_tet_intersection_pts(fmdl,rmdl,top,bot, epsilon)
0408 intpts = [];
0409 pt_idx = unique(rmdl.elems);
0410 pts = rmdl.nodes(pt_idx,1:2);
0411
0412 bb_min = min(...
0413 min(fmdl.nodes(fmdl.faces(:,1),1:2),...
0414 fmdl.nodes(fmdl.faces(:,2),1:2)),...
0415 fmdl.nodes(fmdl.faces(:,3),1:2));
0416 bb_max = max(...
0417 max(fmdl.nodes(fmdl.faces(:,1),1:2),...
0418 fmdl.nodes(fmdl.faces(:,2),1:2)),...
0419 fmdl.nodes(fmdl.faces(:,3),1:2));
0420 todo = ~( bsxfun(@gt,bb_min(:,1),pts(:,1)') ...
0421 | bsxfun(@gt,bb_min(:,2),pts(:,2)') ...
0422 | bsxfun(@lt,bb_max(:,1),pts(:,1)') ...
0423 | bsxfun(@lt,bb_max(:,2),pts(:,2)'));
0424 [F,P] = find(todo);
0425 P = unique(P);
0426 in = false(size(fmdl.faces,1),size(pts,1));
0427 in(F,P) = point_in_triangle(pts(P,:),fmdl.faces(F,:),fmdl.nodes(:,1:2),epsilon)';
0428
0429 [F,P] = find(in);
0430
0431
0432 vf = fmdl.normals(F,3) == 0;
0433 F(vf) = [];
0434 P(vf) = [];
0435
0436
0437
0438 z = sum(fmdl.normals(F,:) .* fmdl.nodes(fmdl.faces(F,1),:),2);
0439
0440 for j = 1:2
0441 z = z - fmdl.normals(F,j) .* pts(P,j);
0442 end
0443 z = z ./ fmdl.normals(F,3);
0444 out = z>top | z < bot;
0445 F(out) = [];
0446 P(out) = [];
0447 z(out) = [];
0448
0449 intpts = [pts(P,:) z];
0450 I = (1:size(intpts,1))';
0451 tri2edge = sparse(F,pt_idx(P),I,size(fmdl.faces,1),size(rmdl.nodes,1));
0452 tri2pts = sparse(F,I,ones(size(I,1),1),size(fmdl.faces,1),size(intpts,1));
0453 edge2pts = sparse(pt_idx(P),I,ones(size(I,1),1),size(rmdl.nodes,1),size(intpts,1));
0454
0455 end
0456
0457 function [intpts, FE2CE, FE2pts, CE2pts] = ...
0458 find_edge2edge_intersections_2d(FE, FN, CE, CN, top, bot)
0459
0460 P1 = FN(FE(:,1),:);
0461 P2 = FN(FE(:,2),:);
0462 P3 = CN(CE(:,1),:);
0463 P4 = CN(CE(:,2),:);
0464
0465 C_bb = zeros(size(P3,1),4);
0466 C_bb(:,[1 3]) = min(P3(:,1:2),P4(:,1:2));
0467 C_bb(:,[2 4]) = max(P3(:,1:2),P4(:,1:2));
0468
0469 F_bb = zeros(size(P1,1),4);
0470 F_bb(:,[1 3]) = min(P1(:,1:2),P2(:,1:2));
0471 F_bb(:,[2 4]) = max(P1(:,1:2),P2(:,1:2));
0472
0473
0474 todo = bsxfun(@gt, F_bb(:,1), C_bb(:,2)') ...
0475 | bsxfun(@lt, F_bb(:,2), C_bb(:,1)') ...
0476 | bsxfun(@gt, F_bb(:,3), C_bb(:,4)') ...
0477 | bsxfun(@lt, F_bb(:,4), C_bb(:,3)');
0478 todo = ~todo;
0479
0480 [T, V] = find(todo);
0481
0482 S1 = P2(T,:) - P1(T,:);
0483 S2 = P4(V,:) - P3(V,:);
0484
0485 denom = S2(:,2) .* S1(:,1) - S2(:,1) .* S1(:,2);
0486
0487 P13 = P1(T,1:2) - P3(V,1:2);
0488
0489 num_a = S2(:,1) .* P13(:,2) ...
0490 - S2(:,2) .* P13(:,1);
0491 num_b = S1(:,1) .* P13(:,2) ...
0492 - S1(:,2) .* P13(:,1);
0493
0494 mua = num_a ./ denom;
0495 mub = num_b ./ denom;
0496
0497 IN = mua>0 & mua<1 & mub>0 & mub<1;
0498 T = T(IN);
0499 V = V(IN);
0500 intpts = P1(T,:) + bsxfun(@times, mua(IN), S1(IN,:));
0501 in = ~(intpts(:,3) > top | intpts(:,3) < bot);
0502 intpts = intpts(in,:);
0503 I = (1:size(intpts,1))';
0504 T = T(in);
0505 V = V(in);
0506 FE2CE = sparse(size(P1,1),size(P3,1));
0507 FE2CE(sub2ind(size(FE2CE),T,V)) = I;
0508 FE2pts = sparse(T,I,ones(size(I)),size(P1,1),size(I,1));
0509 CE2pts = sparse(V,I,ones(size(I)),size(P3,1),size(I,1));
0510
0511 end
0512
0513 function [top_node2tet, top_nodes, nodes_above, nodes_below] = ...
0514 get_cap_nodes_in_tets(fmdl,rmdl,top,epsilon)
0515
0516 top_nodes = rmdl.nodes;
0517 top_nodes(:,3) = top;
0518 nodes_above = fmdl.nodes(:,3) >= top;
0519 nodes_below = fmdl.nodes(:,3) <= top;
0520
0521
0522 elidx = any(nodes_above(fmdl.elems),2) & any(nodes_below(fmdl.elems),2);
0523 top_node2tet = sparse(size(rmdl.nodes,1),size(fmdl.elems,1));
0524 if any(elidx)
0525 mdl = struct;
0526 mdl.nodes = fmdl.nodes;
0527 mdl.elems = fmdl.elems(elidx,:);
0528 top_node2tet(:,elidx) = point_in_tet(mdl,top_nodes,epsilon);
0529 end
0530 end
0531
0532
0533
0534 function [intpts, edge2tri, edge2pts, tri2pts] = ...
0535 get_cap_tet_edge2tri_face_intersections(fmdl,rmdl,top,...
0536 nodes_above,nodes_below, epsilon)
0537
0538 intpts = [];
0539 edge2tri = sparse(size(fmdl.edges,1),size(rmdl.elems,1));
0540 edge2pts = sparse(size(fmdl.edges,1),0);
0541 tri2pts = sparse(size(rmdl.elems,1),0);
0542
0543
0544
0545 edgeidx = any(nodes_above(fmdl.edges),2) & any(nodes_below(fmdl.edges),2);
0546
0547 edgeidx(edgeidx) = fmdl.nodes(fmdl.edges(edgeidx,1),3) ~= ...
0548 fmdl.nodes(fmdl.edges(edgeidx,2),3);
0549
0550 v = fmdl.nodes(fmdl.edges(edgeidx,2),:) - ...
0551 fmdl.nodes(fmdl.edges(edgeidx,1),:);
0552 u = (top - fmdl.nodes(fmdl.edges(edgeidx,1),3)) ./ v(:,3);
0553 pts = fmdl.nodes(fmdl.edges(edgeidx,1),:) + bsxfun(@times,u, v);
0554 t = point_in_triangle(pts(:,1:2),rmdl.elems,rmdl.nodes(:,1:2),-epsilon);
0555
0556
0557 t(sum(t,2)>1,:) = false;
0558 [E,T] = find(t);
0559
0560 if any(t(:))
0561 intpts = pts(E,:);
0562 N = size(intpts,1);
0563 I = (1:N)';
0564 emap = find(edgeidx);
0565 E = emap(E);
0566 edge2tri = sparse(E,T,I,size(fmdl.edges,1), size(rmdl.elems,1));
0567 edge2pts = sparse(E,I,ones(size(I)),size(fmdl.edges,1), N);
0568 tri2pts = sparse(T,I,ones(size(I)),size(rmdl.elems,1), N);
0569 end
0570
0571 end
0572
0573
0574 function [intpts, tri2edge,tri2intpt,edge2intpt] = ...
0575 get_cap_tet_face2tri_edge_intersections( ...
0576 fmdl,rmdl,top,nodes_above,nodes_below, epsilon)
0577
0578 USE_POINT_IN_TRIANGLE = 0;
0579
0580 intpts = [];
0581 tri2edge = sparse(size(fmdl.faces,1),size(rmdl.edges,1));
0582 tri2intpt = sparse(size(fmdl.faces,1),0);
0583 edge2intpt = sparse(size(rmdl.edges,1),0);
0584
0585
0586
0587 faceidx = any(nodes_above(fmdl.faces),2) & any(nodes_below(fmdl.faces),2);
0588
0589 faces = fmdl.faces(faceidx,:);
0590 normals = fmdl.normals(faceidx,:);
0591
0592 N_edges = size(rmdl.edges,1);
0593 N_faces = size(faces,1);
0594
0595
0596 face_bb = zeros(N_faces,6);
0597 face_bb(:,1) = min(reshape(fmdl.nodes(faces,1),N_faces,3),[],2);
0598 face_bb(:,2) = max(reshape(fmdl.nodes(faces,1),N_faces,3),[],2);
0599 face_bb(:,3) = min(reshape(fmdl.nodes(faces,2),N_faces,3),[],2);
0600 face_bb(:,4) = max(reshape(fmdl.nodes(faces,2),N_faces,3),[],2);
0601
0602 edge_bb = zeros(N_edges,6);
0603 edge_bb(:,1) = min(reshape(rmdl.nodes(rmdl.edges,1),N_edges,2),[],2);
0604 edge_bb(:,2) = max(reshape(rmdl.nodes(rmdl.edges,1),N_edges,2),[],2);
0605 edge_bb(:,3) = min(reshape(rmdl.nodes(rmdl.edges,2),N_edges,2),[],2);
0606 edge_bb(:,4) = max(reshape(rmdl.nodes(rmdl.edges,2),N_edges,2),[],2);
0607
0608 todo = bsxfun(@gt, face_bb(:,1), edge_bb(:,2)') ...
0609 | bsxfun(@lt, face_bb(:,2), edge_bb(:,1)') ...
0610 | bsxfun(@gt, face_bb(:,3), edge_bb(:,4)') ...
0611 | bsxfun(@lt, face_bb(:,4), edge_bb(:,3)');
0612 todo = ~todo;
0613 e_todo = any(todo,1);
0614 f_todo = any(todo,2);
0615 faceidx(faceidx) = f_todo;
0616 faces = faces(f_todo,:);
0617 normals = normals(f_todo,:);
0618 [F,E] = find(todo);
0619 emap = zeros(size(e_todo,2),1);
0620 emap(e_todo) = 1:nnz(e_todo);
0621 E = emap(E);
0622 fmap = zeros(size(f_todo,1),1);
0623 fmap(f_todo) = 1:nnz(f_todo);
0624 F = fmap(F);
0625
0626 P1 = rmdl.nodes(rmdl.edges(e_todo,1),:);
0627 P12 = P1 - rmdl.nodes(rmdl.edges(e_todo,2),:);
0628 P1(:,3) = top;
0629 P12(:,3) = 0;
0630
0631 d = sum(fmdl.normals(faceidx,:) .* fmdl.nodes(fmdl.faces(faceidx,1),:),2);
0632
0633 if ~USE_POINT_IN_TRIANGLE
0634
0635 nodes1 = fmdl.nodes(faces(:,1),:);
0636 v0 = fmdl.nodes(faces(:,3),:) - nodes1;
0637 v1 = fmdl.nodes(faces(:,2),:) - nodes1;
0638 dot00 = dot(v0, v0, 2);
0639 dot01 = dot(v0, v1, 2);
0640
0641 dot11 = dot(v1, v1, 2);
0642
0643 invDenom = 1 ./ (dot00 .* dot11 - dot01 .* dot01);
0644 end
0645
0646
0647 num = -d(F) + sum(normals(F,:).*P1(E,:),2);
0648 den = sum(normals(F,:) .* P12(E,:),2);
0649 u = num ./ den;
0650
0651 idx = u >= 0 & u <= 1 & abs(den) >= eps;
0652
0653 if any(idx)
0654 F = F(idx);
0655 E = E(idx);
0656 ipts = P1(E,:) - bsxfun(@times, u(idx), P12(E,:));
0657 if USE_POINT_IN_TRIANGLE
0658 t = point_in_triangle(ipts,faces(F,:),fmdl.nodes,epsilon,'match');
0659 else
0660 v2 = ipts - fmdl.nodes(faces(F,1),:);
0661 dot02 = dot(v0(F,:),v2,2);
0662 dot12 = dot(v1(F,:),v2,2);
0663
0664 dot01 = dot01(F);
0665 invDenom = invDenom(F);
0666 u = (dot11(F) .* dot02 - dot01 .* dot12) .* invDenom;
0667 v = (dot00(F) .* dot12 - dot01 .* dot02) .* invDenom;
0668 t = u >= -epsilon & v >= -epsilon & (u+v-epsilon) <= 1;
0669 end
0670
0671 if any(t)
0672 N = nnz(t);
0673 idv = (1:N)';
0674 intpts = ipts(t,:);
0675 I = idv;
0676 idx = find(faceidx); idx = idx(F);
0677 F = idx(t);
0678 eimap = find(emap);
0679 E = eimap(E(t));
0680
0681 tri2edge = sparse(F,E,I,size(fmdl.faces,1),size(rmdl.edges,1));
0682 tri2intpt = sparse(F,I,ones(size(I)),size(fmdl.faces,1),size(I,1));
0683 edge2intpt = sparse(E,I,ones(size(I)),size(rmdl.edges,1),size(I,1));
0684 end
0685 end
0686 end
0687
0688
0689
0690 function [fmdl,rmdl, opt] = center_scale_models(fmdl,rmdl, opt)
0691 ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0692 rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0693 fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0694 opt.top = opt.top - ctr(3);
0695 opt.bot = opt.bot - ctr(3);
0696 if isfield(opt,'do_not_scale') && opt.do_not_scale
0697 return
0698 end
0699 maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0700 scale = 1/maxnode;
0701 rmdl.nodes = scale*rmdl.nodes;
0702 fmdl.nodes = scale*fmdl.nodes;
0703 opt.top = scale*opt.top;
0704 opt.bot = scale*opt.bot;
0705 opt.height = scale*opt.height;
0706 opt.z_depth= scale*opt.z_depth;
0707 eidors_msg('@@ models scaled by %g', scale,2);
0708 end
0709
0710
0711
0712 function [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl, opt)
0713 f_min = min(fmdl.nodes);
0714 f_max = max(fmdl.nodes);
0715 r_min = min(rmdl.nodes);
0716 r_max = max(rmdl.nodes);
0717
0718 if isfield(opt, 'z_depth') && ~isinf(opt.z_depth)
0719 lvl = mean(rmdl.nodes(:,3));
0720 r_max(3) = lvl + opt.z_depth;
0721 r_min(3) = lvl - opt.z_depth;
0722 else
0723 r_max(3) = f_max(3);
0724 r_min(3) = f_min(3);
0725 end
0726
0727
0728 f_gt = bsxfun(@gt, fmdl.nodes, r_max);
0729 f_lt = bsxfun(@lt, fmdl.nodes, r_min);
0730 r_gt = bsxfun(@gt, rmdl.nodes, f_max);
0731 r_lt = bsxfun(@lt, rmdl.nodes, f_min);
0732
0733
0734 re_gt = any(reshape(all(reshape(r_gt(rmdl.elems',:),3,[])),[],3),2);
0735 re_lt = any(reshape(all(reshape(r_lt(rmdl.elems',:),3,[])),[],3),2);
0736 fe_gt = any(reshape(all(reshape(f_gt(fmdl.elems',:),4,[])),[],3),2);
0737 fe_lt = any(reshape(all(reshape(f_lt(fmdl.elems',:),4,[])),[],3),2);
0738
0739
0740 rmdl_idx = ~(re_gt | re_lt);
0741 fmdl_idx = ~(fe_gt | fe_lt);
0742
0743
0744 rmdl.elems = rmdl.elems(rmdl_idx,:);
0745 fmdl.elems = fmdl.elems(fmdl_idx,:);
0746
0747
0748 [r_used_nodes,jnk,r_n] = unique(rmdl.elems(:));
0749 [f_used_nodes,jnk,f_n] = unique(fmdl.elems(:));
0750
0751 r_idx = 1:numel(r_used_nodes);
0752 f_idx = 1:numel(f_used_nodes);
0753
0754 rmdl.elems = reshape(r_idx(r_n),[],3);
0755 fmdl.elems = reshape(f_idx(f_n),[],4);
0756
0757 rmdl.nodes = rmdl.nodes(r_used_nodes,:);
0758 fmdl.nodes = fmdl.nodes(f_used_nodes,:);
0759
0760
0761 try, rmdl = rmfield(rmdl,'boundary'); end
0762 try, fmdl = rmfield(fmdl,'boundary'); end
0763 end
0764
0765
0766
0767 function fmdl = prepare_tet_mdl(fmdl)
0768 fmopt.elem2edge = true;
0769 fmopt.edge2elem = true;
0770 fmopt.face2elem = true;
0771 fmopt.node2elem = true;
0772 fmopt.normals = true;
0773 ll = eidors_msg('log_level',1);
0774 fmopt.linear_reorder = false;
0775 fmdl = fix_model(fmdl,fmopt);
0776 eidors_msg('log_level',ll);
0777 fmdl.node2elem = logical(fmdl.node2elem);
0778 nElem = size(fmdl.elems,1);
0779 nFace = size(fmdl.faces,1);
0780 fmdl.elem2face = sparse(repmat((1:nElem)',1,4),double(fmdl.elem2face),true,nElem,nFace);
0781 end
0782
0783
0784
0785 function fmdl = prepare_tri_mdl(fmdl)
0786 fmopt.elem2edge = true;
0787 fmopt.edge2elem = true;
0788
0789 fmopt.node2elem = true;
0790
0791 fmopt.linear_reorder = false;
0792 ll = eidors_msg('log_level',1);
0793 fmdl = fix_model(fmdl,fmopt);
0794 eidors_msg('log_level',ll);
0795 fmdl.node2elem = logical(fmdl.node2elem);
0796
0797 end
0798
0799
0800
0801 function debug_plot(fmdl,rmdl,v,t, bot, top, pts)
0802 tet.nodes = fmdl.nodes;
0803 tri.nodes = repmat(rmdl.nodes(rmdl.elems(v,:),:),2,1);
0804 tri.nodes(:,3) = [repmat(bot,3,1); repmat(top,3,1)];
0805 tri.elems = [ 1 2 5 4
0806 2 3 6 5
0807 3 1 4 6];
0808 tri.boundary = tri.elems;
0809 tet.type = 'fwd_model';
0810 tri.type = 'fwd_model';
0811 tet.elems = fmdl.elems(t,:);
0812 clf
0813 show_fem(tri)
0814 hold on
0815 h = show_fem(tet);
0816 set(h,'EdgeColor','b')
0817
0818 plot3(pts(:,1),pts(:,2),pts(:,3),'o');
0819 hold off
0820 axis auto
0821 end
0822
0823
0824
0825 function opt = parse_opts(fmdl,rmdl, opt)
0826
0827 lvl = mean(rmdl.nodes(:,3));
0828
0829 if ~isfield(opt, 'z_depth')
0830 opt.z_depth = Inf;
0831 end
0832 if isinf(opt.z_depth)
0833 opt.top = max(fmdl.nodes(:,3));
0834 opt.bot = min(fmdl.nodes(:,3));
0835 opt.height = opt.top - opt.bot;
0836 else
0837 opt.top = lvl + opt.z_depth;
0838 opt.bot = lvl - opt.z_depth;
0839 opt.height = 2 * opt.z_depth;
0840 end
0841 if ~isfield(opt, 'tol_node2tri');
0842 opt.tol_node2tri = eps;
0843 end
0844 if ~isfield(opt, 'tol_node2tet');
0845 opt.tol_node2tet = eps;
0846 end
0847 if ~isfield(opt, 'tol_edge2edge')
0848 opt.tol_edge2edge = 6*eps(min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))));
0849 end
0850 if ~isfield(opt, 'tol_edge2tri')
0851 opt.tol_edge2tri = eps;
0852 end
0853
0854
0855
0856 eidors_msg('@@ node2tet tolerance = %g', opt.tol_node2tet,2);
0857 eidors_msg('@@ edge2edge tolerance = %g', opt.tol_edge2edge,2);
0858 eidors_msg('@@ edge2tri tolerance = %g', opt.tol_edge2tri,2);
0859 end
0860
0861
0862 function do_unit_test
0863 do_small_test
0864 do_inf_test
0865 do_realistic_test
0866 do_centre_slice;
0867 end
0868
0869 function do_inf_test
0870 jnk = mk_common_model('a2c',16);
0871 tri = jnk.fwd_model;
0872 tri.nodes = 1.1*tri.nodes;
0873 jnk = mk_common_model('c3cr',16);
0874 tet = jnk.fwd_model;
0875
0876 c2f_a = mk_tri2tet_c2f(tet,tri);
0877 c2f_n = mk_approx_c2f(tet,tri);
0878
0879 prob = c2f_n ~= 0 & c2f_a == 0;
0880 unit_test_cmp('Assert no missing intersections', nnz(prob),0);
0881 unit_test_cmp('mk_tri2tet_c2f v approximate', c2f_n,c2f_a, .5);
0882 end
0883
0884
0885 function do_small_test
0886 jnk = mk_common_model('a2c',16);
0887 tri = jnk.fwd_model;
0888 tri.nodes = 1.1*tri.nodes;
0889 jnk = mk_common_model('c3cr',16);
0890 tet = jnk.fwd_model;
0891
0892
0893 opt.z_depth = 0.5;
0894 c2f = mk_tri2tet_c2f(tet,tri,opt);
0895
0896 clf
0897 subplot(131)
0898 show_fem(tet);
0899 hold on
0900 show_fem(tri);
0901 view(2)
0902 subplot(132)
0903 img = mk_image(tet,1);
0904 img.elem_data = sum(c2f,2);
0905 img.calc_colours.clim = 1;
0906 img.calc_colours.ref_level = 0;
0907 show_fem(img);
0908 subplot(133)
0909 img = mk_image(tri,1);
0910 img.elem_data = sum(bsxfun(@times,c2f,get_elem_volume(tet)),1) ./ get_elem_volume(tri)';
0911 img.calc_colours.clim = 1;
0912 img.calc_colours.ref_level = 0;
0913 show_fem(img)
0914
0915 unit_test_cmp('Check C2F size ', size(c2f),[length(tet.elems), length(tri.elems)]);
0916 unit_test_cmp('Check C2F max==1', max(c2f(:)), 1);
0917 unit_test_cmp('Check C2F min==0', min(c2f(:)), 0);
0918
0919 f2c = bsxfun(@rdivide, bsxfun(@times, c2f, get_elem_volume(tet))', get_elem_volume(tri));
0920 unit_test_cmp('Check F2C max==1', max(sum(f2c,2)), 1, 1e-14);
0921 unit_test_cmp('Check F2C min==0', min(f2c(:)), 0);
0922 end
0923
0924
0925 function do_realistic_test
0926
0927 fmdl= ng_mk_cyl_models([2,2,.1],[16,1],[.1,0,.025]);
0928 xvec = -2:.5:2;
0929 yvec = -2:.5:2;
0930 zvec = [.5 1.5];
0931 grid2d = mk_grid_model([],xvec,yvec);
0932 grid3d = mk_grid_model([],xvec,yvec,zvec);
0933 eidors_cache clear
0934 tic
0935 opt.save_memory = 0;
0936 c2f_a = mk_grid_c2f(fmdl, grid3d,opt);
0937 t = toc;
0938 fprintf('Voxel:\tt=%f s\n',t);
0939
0940 opt.z_depth = .5;
0941 grid2d.nodes(:,3) = 1;
0942 eidors_cache clear
0943 tic
0944
0945
0946 c2f_b = mk_tri2tet_c2f(fmdl, grid2d, opt);
0947 t = toc;
0948 fprintf('tri2tet:\tt=%f \n',t);
0949 c2f_c = c2f_b * grid2d.coarse2fine;
0950
0951 eidors_cache clear
0952
0953 grid2d.mk_coarse_fine_mapping.z_depth = .5;
0954 grid2d.mk_coarse_fine_mapping.f2c_offset(3) = 1;
0955 grid2d = rmfield(grid2d,'coarse2fine');
0956 tic
0957 c2f_n = mk_approx_c2f(fmdl,grid2d);
0958 t = toc;
0959 fprintf('Approximate: t=%f s\n',t);
0960
0961 unit_test_cmp('mk_tri2tet_c2f v mk_grid_c2f', c2f_c,c2f_a, 1e-5);
0962 unit_test_cmp('mk_tri2tet_c2f v approximate', c2f_n,c2f_b, .5);
0963 prob = c2f_n ~= 0 & c2f_b == 0;
0964 unit_test_cmp('Assert no missing intersections', nnz(prob),0);
0965
0966
0967
0968
0969
0970
0971
0972
0973
0974
0975
0976
0977
0978
0979
0980
0981
0982
0983
0984
0985 end
0986
0987
0988 function do_centre_slice;
0989 imdl = mk_common_model('b3cr',[16,2]);
0990 f_mdl= imdl.fwd_model;
0991 imdl2d= mk_common_model('b2c2',16);
0992 c_mdl= imdl2d.fwd_model;
0993 c2f= mk_coarse_fine_mapping( f_mdl, c_mdl);
0994 end