0001 function [c2f, m] = mk_grid_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
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061 persistent save_memory;
0062
0063 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'), do_unit_test; return; end
0064 if ischar(fmdl) && strcmp(fmdl,'PROFILE'), do_small_test; return; end
0065 if ischar(fmdl) && strcmp(fmdl,'save_memory')
0066 if nargin ~= 2, error('Expected a value for save_memory option'); end
0067 if ischar(rmdl), rmdl = str2double(rmdl); end
0068 save_memory = rmdl;
0069 return
0070 end
0071 if nargin < 3
0072 opt = struct;
0073 end
0074 if ~isempty(save_memory)
0075 try
0076 opt.save_memory;
0077 catch
0078 opt.save_memory = save_memory;
0079 end
0080 end
0081
0082 [fmdl,rmdl] = eidors_cache(@center_scale_models,{fmdl,rmdl, opt});
0083
0084 opt = eidors_cache(@parse_opts,{fmdl,rmdl, opt});
0085
0086 copt.cache_obj = {fmdl.nodes,
0087 fmdl.elems,
0088 rmdl.nodes,
0089 rmfield(opt,'save_memory')};
0090
0091 copt.fstr = 'mk_grid_c2f';
0092
0093 [c2f, m] = eidors_cache(@do_mk_grid_c2f,{fmdl,rmdl,opt},copt);
0094
0095
0096
0097
0098
0099 function [c2f, m]= do_mk_grid_c2f(fmdl0,rmdl0,opt0)
0100 eidors_msg('@@@',2);
0101
0102 m = [];
0103 c2f = sparse(opt0.nTet,opt0.nVox);
0104
0105 progress_msg('Prepare tet model...');
0106 fmdl0 = prepare_fmdl(fmdl0);
0107 progress_msg(Inf);
0108
0109 if opt0.save_memory == 0
0110 relem_idx = ones(opt0.nVox,1);
0111 [fmdl, rmdl, opt, felem_idx] = crop_models(fmdl0,rmdl0,opt0,relem_idx);
0112 if any(felem_idx)
0113 [tmp, m] = separable_calculations(fmdl,rmdl0,opt);
0114 c2f = combine_c2f(c2f, tmp,felem_idx,relem_idx);
0115 end
0116 elseif opt0.save_memory == 10
0117
0118
0119
0120 logmsg(' Saving memory mode level %d\n',opt0.save_memory);
0121 n_elems = num_elems(fmdl0); step = 5e4;
0122 max_iter = floor(n_elems/step);
0123 pmopt.final_msg = '';
0124 ctr = interp_mesh(fmdl0);
0125 [~,xidx] = sort(ctr(:,1));
0126 [xl] = ndgrid(opt0.xvec(1:end-1),opt0.yvec(1:end-1),opt0.zvec(1:end-1)); xl=xl(:);
0127 [xu] = ndgrid(opt0.xvec(2:end-0),opt0.yvec(2:end-0),opt0.zvec(2:end-0)); xu=xu(:);
0128 progress_msg('Progress:',0,max_iter+1,pmopt);
0129 progress = 0;
0130 for k = 0:max_iter;
0131 eidx = true(n_elems,1);
0132 eidx( xidx((k*step+1):min((k+1)*step,n_elems)) ) = false;
0133 fmdl = fmdl0; fmdl.elems(eidx,:) = [];
0134 fmdl.edge2elem(:,eidx)= [];
0135 fmdl.node2elem(:,eidx)= [];
0136 fmdl.elem2face(eidx,:)= [];
0137 opt = opt0; opt.nTet = sum(eidx==0);
0138 if 1
0139 xvals = fmdl.nodes(fmdl.elems(:),1);
0140 ridx = true(size(xu));
0141 ridx(xl>max(xvals)) = false;
0142 ridx(xu<min(xvals)) = false;
0143 rmdl = rmdl0;
0144 idx = rmdl0.coarse2fine * ridx > 0;
0145 rmdl.elems = rmdl0.elems(idx,:);
0146 [node_idx, ~,Nn] = unique(rmdl.elems(:));
0147 rmdl.elems = reshape(Nn,size(rmdl.elems));
0148 rmdl.nodes = rmdl0.nodes(node_idx,:);
0149 opt.xvec = unique(rmdl.nodes( rmdl.elems(:), 1));
0150 opt.Xsz = length(opt.xvec)-1;
0151 opt.ystep = opt.xstep*opt.Xsz+1;
0152 opt.zstep = opt.ystep*(opt.Ysz+1);
0153 opt.nVox = opt.Xsz*opt.Ysz*opt.Zsz;
0154
0155 if 0
0156 fmdl.boundary = find_boundary(fmdl);
0157 rmdl.boundary = find_boundary(rmdl);
0158 hh=show_fem(rmdl); set(hh,'EdgeColor',[0,0,1]);
0159 hold on; show_fem(fmdl); hold off; keyboard
0160 end
0161 else
0162 rmdl = rmdl0; ridx= true(opt.nVox,1);
0163 end
0164 if isempty(rmdl.elems); continue; end
0165
0166 eidors_msg('log_level',eidors_msg('log_level')-2);
0167 c2f(~eidx,ridx) = separable_calculations(fmdl,rmdl,opt);
0168 eidors_msg('log_level',eidors_msg('log_level')+2);
0169 progress_msg(k+1, max_iter+1);
0170 end
0171 progress_msg(Inf);
0172 else
0173 logmsg(' Saving memory mode level %d\n',opt0.save_memory);
0174 max_iter = opt0.Zsz;
0175 if opt0.save_memory >= 2, max_iter = max_iter * opt0.Ysz; end
0176 if opt0.save_memory == 3, max_iter = max_iter * opt0.Xsz; end
0177 pmopt.final_msg = '';
0178 progress_msg('Progress:',0,max_iter,pmopt);
0179 progress = 0;
0180 opt = opt0;
0181 m = prepare_vox_mdl(rmdl0,opt0);
0182 for z = 1:opt0.Zsz
0183 opt.Zsz = 1;
0184 opt.zvec = opt0.zvec(z:z+1);
0185 opt.xplane = opt0.xplane(z,:);
0186 opt.yplane = opt0.yplane(z,:);
0187 relem_idx = false(opt0.Xsz*opt0.Ysz*opt0.Zsz,1);
0188 relem_idx((z-1)*opt0.Xsz*opt0.Ysz + (1:opt0.Xsz*opt0.Ysz)) = true;
0189 if opt0.save_memory > 1
0190 for y = 1:opt0.Ysz
0191 opt.Ysz = 1;
0192 opt.yvec = opt0.yvec(y:y+1);
0193 opt.zplane = opt0.zplane(y,:);
0194 opt.xplane = opt0.xplane(z,y);
0195 opt.zstep = opt.ystep*(opt.Ysz+1);
0196 relem_idx_y = false(opt0.Xsz*opt0.Ysz,1);
0197 relem_idx_y((y-1)*opt0.Xsz + (1:opt0.Xsz)) = true;
0198 relem_idx_y = relem_idx & repmat(relem_idx_y,opt0.Zsz,1);
0199 if opt0.save_memory > 2
0200 for x = 1:opt0.Xsz
0201 opt.Xsz = 1;
0202 opt.xvec = opt0.xvec(x:x+1);
0203 opt.yplane = opt0.yplane(z,x);
0204 opt.zplane = opt0.zplane(y,x);
0205 opt.ystep = opt.xstep*opt.Xsz+1;
0206 opt.zstep = opt.ystep*(opt.Ysz+1);
0207 relem_idx_x = false(opt0.Xsz,1);
0208 relem_idx_x(x) = true;
0209 relem_idx_x = relem_idx_y & repmat(relem_idx_x,opt0.Ysz*opt0.Zsz,1);
0210 [fmdl, rmdl, opt, felem_idx] = crop_models(fmdl0,rmdl0,opt,relem_idx_x);
0211 if any(felem_idx)
0212 eidors_msg('log_level',eidors_msg('log_level')-2);
0213 tmp = separable_calculations(fmdl,rmdl,opt);
0214 eidors_msg('log_level',eidors_msg('log_level')+2);
0215 c2f = combine_c2f(c2f, tmp,felem_idx,relem_idx_x);
0216 end
0217 progress = progress + 1;
0218 progress_msg(progress,max_iter);
0219 end
0220 else
0221 [fmdl, rmdl, opt, felem_idx] = crop_models(fmdl0,rmdl0,opt,relem_idx_y);
0222 if any(felem_idx)
0223 eidors_msg('log_level',eidors_msg('log_level')-2);
0224 tmp = separable_calculations(fmdl,rmdl,opt);
0225 eidors_msg('log_level',eidors_msg('log_level')+2);
0226 c2f = combine_c2f(c2f, tmp,felem_idx,relem_idx_y);
0227 end
0228 progress = progress + 1;
0229 progress_msg(progress, max_iter);
0230 end
0231 end
0232
0233 else
0234 [fmdl, rmdl, opt, felem_idx] = crop_models(fmdl0,rmdl0,opt,relem_idx);
0235 if any(felem_idx)
0236 eidors_msg('log_level',eidors_msg('log_level')-2);
0237 tmp = separable_calculations(fmdl,rmdl,opt);
0238 eidors_msg('log_level',eidors_msg('log_level')+2);
0239 c2f = combine_c2f(c2f, tmp,felem_idx,relem_idx);
0240 end
0241 progress = progress + 1;
0242 progress_msg(progress, max_iter);
0243 end
0244
0245
0246 end
0247 progress_msg(Inf);
0248 end
0249
0250
0251
0252 function c2f = combine_c2f(c2f, tmp,felem_idx,relem_idx)
0253 fidx = find(felem_idx);
0254 ridx = find(relem_idx);
0255 c2f(fidx,ridx) = c2f(fidx,ridx) + tmp;
0256
0257
0258
0259 function [c2f, m] = separable_calculations(fmdl,rmdl,opt)
0260 DEBUG = eidors_debug('query','mk_grid_c2f');
0261
0262 progress_msg('Prepare vox model...');
0263 m = prepare_vox_mdl(rmdl,opt);
0264 progress_msg(Inf);
0265
0266 try
0267
0268
0269 progress_msg('Find tet_edge2vox_face intersections...')
0270 [intpts1, rec2tedge, rec2intpt1, tedge2intpt1] = get_voxel_intersection_points(fmdl,m.faces,opt);
0271 progress_msg(sprintf('Found %d', size(intpts1,1)), Inf);
0272
0273
0274 progress_msg('Find vox_edge2tet_face intersections...',0,3)
0275 [intpts2, tri2vedge, tri2intpt2, vedge2intpt2] = get_tet_intersection_points(fmdl,m,opt);
0276 progress_msg(sprintf('Found %d', size(intpts2,1)), Inf);
0277
0278
0279
0280
0281
0282
0283
0284 pmopt.final_msg = 'none';
0285 progress_msg('Find tet_edge2vox_edge intersections...',-1,pmopt)
0286 [intpts3, tedge2vedge, tedge2intpt3, vedge2intpt3] =find_edge2edge_intersections(fmdl.edges,fmdl.nodes,m.edges,rmdl.nodes, opt.tol_edge2edge);
0287 progress_msg(sprintf('Found %d',size(intpts3,1)),Inf);
0288
0289
0290 progress_msg('Find tet_nodes in voxels...')
0291 [fnode2vox] = get_nodes_in_voxels(fmdl,rmdl);
0292 progress_msg(sprintf('Found %d', nnz(fnode2vox)), Inf);
0293
0294
0295 progress_msg('Find vox_nodes in tets...');
0296 vnode2tet = get_nodes_in_tets(fmdl,rmdl, opt);
0297 progress_msg(sprintf('Found %d', nnz(vnode2tet)), Inf);
0298
0299
0300 progress_msg('Find vox contained in tet...')
0301 vox_in_tet = (m.node2vox' * vnode2tet) == 8;
0302 progress_msg(sprintf('Found %d',nnz(vox_in_tet)), Inf);
0303
0304
0305 progress_msg('Find tets contained in vox...');
0306 tet_in_vox = (double(fmdl.node2elem') * fnode2vox) == 4;
0307 progress_msg(sprintf('Found %d',nnz(tet_in_vox)), Inf);
0308
0309
0310
0311 progress_msg('Find total vox v. tet intersections...');
0312 vox2intTet = m.vox2face * (rec2tedge>0) * fmdl.edge2elem ...
0313 | m.vox2edge * (tri2vedge>0)' * fmdl.elem2face' ...
0314 | m.vox2edge * tedge2vedge' * fmdl.edge2elem;
0315
0316 vox2intTet = vox2intTet & ~vox_in_tet & ~tet_in_vox';
0317 progress_msg(sprintf('Found %d',nnz(vox2intTet)), Inf);
0318
0319 catch err
0320 if (strcmp(err.identifier,'MATLAB:nomem'))
0321 msg = sprintf('%s', ...
0322 'Matlab ran out of memory. Consider setting save_memory > 0.\n', ...
0323 'See HELP MK_GRID_C2F for details.');
0324 error('EIDORS:mk_grid_c2f:nomem',msg);
0325 else
0326 rethrow(err);
0327 end
0328 end
0329
0330 progress_msg('Calculate intersection volumes...');
0331
0332 vox2intpt1 = logical(m.vox2face*rec2intpt1)';
0333 tet2intpt1 = logical(fmdl.edge2elem'*tedge2intpt1)';
0334
0335 tet2intpt2 = logical(fmdl.elem2face*tri2intpt2)';
0336 vox2intpt2 = logical(m.vox2edge*vedge2intpt2)';
0337
0338 tet2intpt3 = logical(fmdl.edge2elem'*tedge2intpt3)';
0339 vox2intpt3 = logical(m.vox2edge*vedge2intpt3)';
0340
0341 vox_todo = find(sum(vox2intTet,2)>0);
0342 C = []; F = []; V = [];
0343
0344 id = 0; lvox = length(vox_todo);
0345 mint = ceil(lvox/100);
0346 for v = vox_todo'
0347 id = id+1;
0348 if mod(id,mint)==0, progress_msg(id/lvox); end
0349 tet_todo = find(vox2intTet(v,:));
0350 common_intpts1 = bsxfun(@and,vox2intpt1(:,v), tet2intpt1(:,tet_todo));
0351 common_intpts2 = bsxfun(@and,vox2intpt2(:,v), tet2intpt2(:,tet_todo));
0352 common_intpts3 = bsxfun(@and,vox2intpt3(:,v), tet2intpt3(:,tet_todo));
0353 tet_nodes = bsxfun(@and,fnode2vox(:,v), fmdl.node2elem(:,tet_todo));
0354 vox_nodes = bsxfun(@and,vnode2tet(:,tet_todo), m.node2vox(:,v));
0355
0356 C = [C; v*ones(numel(tet_todo),1)];
0357 F = [F; tet_todo'];
0358 last_v = numel(V);
0359 V = [V; zeros(numel(tet_todo),1)];
0360
0361 for t = 1:numel(tet_todo)
0362 pts = [ intpts1(common_intpts1(:,t),:);
0363 intpts2(common_intpts2(:,t),:);
0364 intpts3(common_intpts3(:,t),:);
0365 fmdl.nodes(tet_nodes(:,t),:);
0366 rmdl.nodes(vox_nodes(:,t),:)];
0367 last_v = last_v + 1;
0368 ok = false;
0369 if size(pts,1) < 4
0370
0371
0372 continue
0373 end
0374 if size(pts,1) < 4
0375
0376
0377 E = fmdl.edges(fmdl.elem2edge(tet_todo(t),:),:);
0378 P1 = fmdl.nodes(E(:,1),:);
0379 P2 = fmdl.nodes(E(:,2),:);
0380
0381
0382
0383 D = P1-P2;
0384 ok = any(abs(D(:)) <= eps);
0385 end
0386 if ok; continue, end
0387 try
0388
0389
0390 ctr = mean(pts);
0391 pts = bsxfun(@minus,pts,ctr);
0392 scale = max(abs(pts(:)));
0393 if scale == 0
0394 continue
0395 end
0396
0397 pts = pts ./ scale;
0398
0399
0400 [~, V(last_v)] = convhulln(pts,{'Qt Pp Qs'});
0401 V(last_v) = V(last_v) * scale^3;
0402 catch err
0403 ok = false;
0404 switch err.identifier
0405 case {'MATLAB:qhullmx:DegenerateData'}
0406
0407 ok = 1;
0408
0409 case {'MATLAB:qhullmx:DegenerateData', 'MATLAB:qhullmx:UndefinedError', ...
0410 'MATLAB:cgprechecks:NotEnoughPts'}
0411
0412
0413 E = fmdl.edges(fmdl.elem2edge(tet_todo(t),:),:);
0414 P1 = fmdl.nodes(E(:,1),:);
0415 P2 = fmdl.nodes(E(:,2),:);
0416
0417
0418
0419 D = P1-P2;
0420 ok = any(abs(D) < 10*eps);
0421
0422
0423
0424 u = uniquetol(pts*scale,10*eps,'ByRows',true,'DataScale', 1);
0425 ok = ok | size(u,1) < 4;
0426 end
0427 if ~ok & size(u,1) == 4
0428 ok = ok | det([ones(4,1),u])<eps;
0429 end
0430 if ~ok
0431 if DEBUG || eidors_debug('query','mk_grid_c2f:convhulln')
0432 tet.nodes = fmdl.nodes;
0433 vox.nodes = rmdl.nodes;
0434 tet.type = 'fwd_model';
0435 vox.type = 'fwd_model';
0436 vox.elems = m.faces(logical(m.vox2face(v,:)),:);
0437 vox.boundary = vox.elems;
0438 tet.elems = fmdl.elems(tet_todo(t),:);
0439
0440 clf
0441 show_fem(vox)
0442 hold on
0443 h = show_fem(tet);
0444 set(h,'EdgeColor','b')
0445 pts = bsxfun(@plus,pts*scale,ctr);
0446 plot3(pts(:,1),pts(:,2),pts(:,3),'o');
0447
0448
0449 hold off
0450 axis auto
0451 keyboard
0452 else
0453 fprintf('\n');
0454 eidors_msg(['convhulln has thrown an error. ' ...
0455 'Enable eidors_debug on mk_grid_c2f:convhulln and re-run to see a debug plot'],0);
0456 rethrow(err);
0457 end
0458 end
0459 end
0460 end
0461 end
0462 progress_msg(Inf);
0463 c2f = sparse(F,C,V,opt.nTet,opt.nVox);
0464
0465
0466 c2f = c2f + bsxfun(@times, sparse(vox_in_tet), m.volume)';
0467
0468
0469 vol = get_elem_volume(fmdl);
0470 c2f = bsxfun(@rdivide,c2f,vol);
0471
0472
0473
0474 c2f = c2f + tet_in_vox;
0475
0476
0477
0478 function [fmdl, rmdl, opt, felem_idx] = crop_models(fmdl0,rmdl0,opt, relem_idx)
0479
0480 fmdl = fmdl0;
0481 rmdl = rmdl0;
0482
0483 opt.nVox = opt.Xsz*opt.Ysz*opt.Zsz;
0484 idx = rmdl0.coarse2fine * relem_idx > 0;
0485 rmdl.elems = rmdl0.elems(idx,:);
0486 [node_idx, m,n] = unique(rmdl.elems(:));
0487 rmdl.elems = reshape(n,size(rmdl.elems));
0488 rmdl.nodes = rmdl0.nodes(node_idx,:);
0489
0490 fnode_above = fmdl0.nodes(:,3) > opt.zvec(end);
0491
0492 felem_above = all(reshape(fnode_above(fmdl0.elems), size(fmdl0.elems)),2);
0493
0494 fnode_below = fmdl0.nodes(:,3) < opt.zvec(1);
0495 felem_below = all(reshape(fnode_below(fmdl0.elems), size(fmdl0.elems)),2);
0496
0497 felem_idx = ~(felem_below | felem_above);
0498 fmdl.elems = fmdl0.elems(felem_idx,:);
0499 fnode_idx = unique(fmdl.elems(:));
0500 fnode_idx_map = zeros(size(fmdl0.nodes,1),1);
0501 fnode_idx_map(fnode_idx) = 1:length(fnode_idx);
0502 fmdl.elems = reshape(fnode_idx_map(fmdl.elems),size(fmdl.elems));
0503 fmdl.edges = reshape(fnode_idx_map(fmdl0.edges),size(fmdl0.edges));
0504 fmdl.nodes = fmdl0.nodes(fnode_idx,:);
0505 fmdl.node2elem = fmdl0.node2elem(fnode_idx,felem_idx);
0506
0507
0508 fface_idx = sum(fmdl0.elem2face(felem_idx,:),1)>0;
0509 fmdl.elem2face = fmdl0.elem2face(felem_idx,fface_idx);
0510 felem_idx_map = felem_idx;
0511 felem_idx_map(felem_idx) = 1:nnz(felem_idx);
0512 felem_idx_map = [0; felem_idx_map];
0513 fmdl.face2elem = fmdl0.face2elem(fface_idx,:);
0514 fmdl.face2elem = reshape(felem_idx_map(fmdl.face2elem + 1),size(fmdl.face2elem));
0515 fmdl.normals = fmdl0.normals(fface_idx,:);
0516 fmdl.faces = fmdl0.faces(fface_idx,:);
0517 fmdl.faces = reshape(fnode_idx_map(fmdl.faces),size(fmdl.faces));
0518
0519 fedge_idx = unique(fmdl0.elem2edge(felem_idx,:));
0520 fedge_idx_map = zeros(size(fmdl0.edges,1),1);
0521 fedge_idx_map(fedge_idx) = 1:length(fedge_idx);
0522 fmdl.elem2edge = fedge_idx_map(fmdl0.elem2edge(felem_idx,:));
0523 fmdl.edge2elem = fmdl0.edge2elem(fedge_idx,felem_idx);
0524 fmdl.edges = fmdl.edges(fedge_idx,:);
0525
0526 opt.nTet = size(fmdl.elems,1);
0527
0528
0529
0530 function fmdl = prepare_fmdl(fmdl)
0531 fmopt.elem2edge = true;
0532 fmopt.edge2elem = true;
0533 fmopt.face2elem = true;
0534 fmopt.node2elem = true;
0535 fmopt.normals = true;
0536 fmopt.linear_reorder = false;
0537 ll = eidors_msg('log_level',1);
0538 fmdl = fix_model(fmdl,fmopt);
0539 eidors_msg('log_level',ll);
0540 fmdl.node2elem = logical(fmdl.node2elem);
0541 nElem = size(fmdl.elems,1);
0542 nFace = size(fmdl.faces,1);
0543 fmdl.elem2face = sparse(repmat((1:nElem)',1,4),double(fmdl.elem2face),true,nElem,nFace);
0544
0545
0546
0547
0548 function m = prepare_vox_mdl(rmdl,opt)
0549
0550 DEBUG = eidors_debug('query','mk_grid_c2f');
0551
0552 [voxels, m.node2vox] = mk_voxels(opt);
0553 if DEBUG || eidors_debug('query','mk_grid_c2f:mk_voxels')
0554 show_voxels(rmdl,voxels); title('mk\_voxels');
0555 end
0556
0557 m.faces = mk_faces(voxels,opt);
0558 if DEBUG || eidors_debug('query','mk_grid_c2f:mk_faces')
0559 test_faces(rmdl,m.faces,opt); title('mk\_faces');
0560 end
0561
0562 m.vox2face = mk_vox2face(opt);
0563 if DEBUG || eidors_debug('query','mk_grid_c2f:mk_vox2face')
0564
0565 show_voxels(rmdl,m.faces(any(m.vox2face),:));title('mk\_vox2face');
0566 end
0567 m.edges = mk_edges(voxels,opt);
0568
0569 m.edge_length = rmdl.nodes(m.edges(:,1),:) - rmdl.nodes(m.edges(:,2),:);
0570 m.edge_length = sqrt(sum(m.edge_length.^2,2));
0571
0572 [m.vox2edge, m.volume]= mk_vox2edge(m,opt);
0573
0574
0575
0576
0577
0578 function rnode2tet = get_nodes_in_tets(fmdl,rmdl, opt)
0579
0580 [A,b] = tet_to_inequal(fmdl.nodes,fmdl.elems);
0581 progress_msg(.01);
0582
0583 rnode2tet = (bsxfun(@minus, A(1:4:end,:)*rmdl.nodes',b(1:4:end)) <= opt.tol_node2tet)';
0584 progress_msg(.21);
0585 for i = 2:4
0586 rnode2tet = rnode2tet & (bsxfun(@minus, A(i:4:end,:)*rmdl.nodes',b(i:4:end)) <= opt.tol_node2tet)';
0587 progress_msg(.21 + (i-1)*.23);
0588 end
0589
0590
0591 ex= bsxfun(@eq,rmdl.nodes(:,1),fmdl.nodes(:,1)') & ...
0592 bsxfun(@eq,rmdl.nodes(:,2),fmdl.nodes(:,2)') & ...
0593 bsxfun(@eq,rmdl.nodes(:,3),fmdl.nodes(:,3)');
0594 progress_msg(.94);
0595 rnode2tet(any(ex,2),:) = 0;
0596 progress_msg(1);
0597
0598
0599
0600
0601
0602 function [insnode] = get_nodes_in_voxels(fmdl,rmdl)
0603
0604 E = reshape(rmdl.elems',4*6,[])';
0605 E = E(:,[1 2 3 4 8 12 16 23]);
0606
0607 NE = size(E,1);
0608 xnodes = reshape(rmdl.nodes(E,1),NE,[]);
0609 ynodes = reshape(rmdl.nodes(E,2),NE,[]);
0610 znodes = reshape(rmdl.nodes(E,3),NE,[]);
0611 minx = min(xnodes,[],2);
0612 maxx = max(xnodes,[],2);
0613 miny = min(ynodes,[],2);
0614 maxy = max(ynodes,[],2);
0615 minz = min(znodes,[],2);
0616 maxz = max(znodes,[],2);
0617
0618 leftof = bsxfun(@lt, fmdl.nodes(:,1)+eps, minx');
0619 rightof = bsxfun(@gt, fmdl.nodes(:,1)-eps, maxx');
0620 infront = bsxfun(@lt, fmdl.nodes(:,2)+eps, miny');
0621 behind = bsxfun(@gt, fmdl.nodes(:,2)-eps, maxy');
0622 below = bsxfun(@lt, fmdl.nodes(:,3)+eps, minz');
0623 above = bsxfun(@gt, fmdl.nodes(:,3)-eps, maxz');
0624
0625 outnode = leftof | rightof | behind | infront | below | above;
0626 insnode = sparse(~outnode);
0627
0628
0629
0630 function [intpts, face2edge, face2intpt, edge2intpt] = get_voxel_intersection_points(fmdl,faces,opt)
0631 edges = fmdl.edges;
0632 nodes = fmdl.nodes;
0633 dir = nodes(edges(:,2),:) - nodes(edges(:,1),:);
0634 intpts = [];
0635 F = []; E = []; I = [];
0636
0637 SZ = [opt.Xsz, opt.Ysz, opt.Zsz];
0638 VEC = {opt.xvec, opt.yvec, opt.zvec};
0639
0640
0641
0642
0643
0644 todo = sum(SZ)+3;
0645 id = 0;
0646 step = ceil(todo/100);
0647
0648 for d = 1:3
0649 for x = 0:SZ(d)
0650 id = id+1;
0651 if mod(id,step)==0, progress_msg(id/todo); end
0652
0653 t = (VEC{d}(x+1) - nodes(edges(:,1),d)) ./ dir(:,d);
0654 crossing = t>0 & t<1;
0655 crossed = find(crossing);
0656 if isempty(crossed), continue, end;
0657
0658 tmp = nodes(edges(crossing,1),:) + bsxfun(@times,t(crossing),dir(crossing,:));
0659 face_coord = zeros(length(crossed),3);
0660 rmv = false(length(crossed),1);
0661 for i = 1:3
0662 if i == d
0663 face_coord(:,i) = x+1;
0664 else
0665 face_coord(:,i) = sum(bsxfun(@times,diff(bsxfun(@lt, tmp(:,i), VEC{i}'),1,2),(1:SZ(i))),2);
0666 rmv = rmv | face_coord(:,i) == 0;
0667 end
0668 end
0669
0670 same_x = any(bsxfun(@eq,tmp(:,1),opt.xvec'),2);
0671 same_y = any(bsxfun(@eq,tmp(:,2),opt.yvec'),2);
0672 same_z = any(bsxfun(@eq,tmp(:,3),opt.zvec'),2);
0673 rmv = rmv | (same_x + same_y + same_z) > 1;
0674 if nnz(rmv) == numel(rmv), continue, end
0675 tmp(rmv,:) = [];
0676 face_coord(rmv,:) = [];
0677 crossed(rmv) = [];
0678 face_coord = face_coord - 1;
0679 face_idx = d + face_coord(:,1)*opt.xstep + face_coord(:,2)*opt.ystep + face_coord(:,3)*opt.zstep;
0680 I = [I; (1:size(tmp,1))' + size(I,1)];
0681 F = [F; face_idx];
0682 E = [E; crossed];
0683 intpts = [intpts; tmp];
0684
0685
0686
0687
0688
0689
0690
0691
0692
0693
0694
0695
0696
0697
0698
0699
0700
0701 end
0702 end
0703 face2edge = sparse(F,E,I,size(faces,1),size(edges,1));
0704 face2intpt = sparse(F,I,ones(size(I)),size(faces,1),size(I,1));
0705 edge2intpt = sparse(E,I,ones(size(I)),size(edges,1),size(I,1));
0706
0707
0708
0709 function [intpts, tri2edge, tri2intpt, edge2intpt] = get_tet_intersection_points(fmdl,m,opt)
0710
0711 intpts = [];
0712 T = []; E = []; I = [];
0713 Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0714 SZ = [opt.Xsz, opt.Ysz, opt.Zsz];
0715 VEC = {opt.xvec, opt.yvec, opt.zvec};
0716 STEP(1) = opt.xstep;
0717 STEP(2) = STEP(1)*(Xsz+1);
0718 STEP(3) = STEP(2)*(Ysz+1);
0719
0720 d = sum(fmdl.normals .* fmdl.nodes(fmdl.faces(:,1),:),2);
0721
0722 line_axis = [3 1 2];
0723 for i = 1:3
0724 progress_msg(i,3);
0725
0726 idx = 1:3;
0727 op = line_axis(i);
0728 idx(op) = [];
0729
0730 pts = [repmat(VEC{idx(1)},SZ(idx(2))+1,1) kron(VEC{idx(2)},ones(SZ(idx(1))+1,1)) ];
0731 pt_idx = uint32(repmat((0:SZ(idx(1)))',SZ(idx(2))+1,1)*STEP(idx(1)) ...
0732 + kron((0:SZ(idx(2)))', ones(SZ(idx(1))+1,1))*STEP(idx(2)));
0733
0734
0735
0736 z = repmat(d,1,size(pts,1));
0737 for j = 1:2
0738 z = z - bsxfun(@times,fmdl.normals(:,idx(j)),pts(:,j)');
0739 end
0740 z = z ./ repmat(fmdl.normals(:,op),1,size(pts,1));
0741 in = point_in_triangle(pts,fmdl.faces,fmdl.nodes(:,idx),2*opt.tol_edge2tri)';
0742
0743
0744
0745 v = VEC{op}';
0746 in = in & ~reshape(any(bsxfun(@eq,z(:),v),2),size(in));
0747 M = bsxfun(@lt, z(:),v);
0748 M = xor(M(:,1:end-1), M(:,2:end));
0749
0750 edge_num = reshape(uint32(M*(1:SZ(op))'), size(z));
0751
0752 in = in & edge_num > 0;
0753 if nnz(in) == 0, continue, end
0754 edge_num(~in) = 0;
0755
0756 edge_num(in) = (edge_num(in)-1) * STEP(op);
0757 edge_idx = edge_num + bsxfun(@times,uint32(full(in)), uint32(i) + pt_idx');
0758
0759 [t, p] = find(in);
0760 tmp = zeros(length(p),3);
0761 tmp(:,idx) = pts(p,1:2);
0762 tmp(:,op) = z(in);
0763
0764 I = [I; (1:size(tmp,1))' + size(I,1)];
0765 T = [T; t];
0766 E = [E; edge_idx(in)];
0767 intpts = [intpts; tmp];
0768 end
0769
0770 tri2edge = sparse(T,E,I,size(fmdl.faces,1),size(m.edges,1));
0771 tri2intpt = sparse(T,I,ones(size(I)),size(fmdl.faces,1),size(I,1));
0772 edge2intpt = sparse(E,I,ones(size(I)),size(m.edges,1),size(I,1));
0773
0774
0775
0776
0777 function [voxels, node2vox] = mk_voxels(opt)
0778 Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0779 Xp = Xsz+1; Yp = Ysz+1; Zp = Zsz+1;
0780
0781 up = Xp*Yp;
0782 vox = [ 1 1+up 1+up+Xp 1+Xp;
0783 1 2 2+up 1+up;
0784 1 1+Xp 2+Xp 2;
0785 2 2+up 2+up+Xp 2+Xp;
0786 1+Xp 2+Xp 2+up+Xp 1+up+Xp;
0787 1+up 1+Xp+up 2+Xp+up 2+up];
0788
0789 voxrow = bsxfun(@plus, repmat(vox, Xsz,1), kron(0:Xsz-1 ,ones(1,6))');
0790 voxplane = bsxfun(@plus, repmat(voxrow, Ysz,1), kron(Xp*(0:Ysz-1) , ones(1, 6*Xsz))');
0791 voxels = bsxfun(@plus, repmat(voxplane,Zsz,1), kron(up*(0:Zsz-1) , ones(1, 6*Xsz*Ysz))');
0792 nVox = Xsz * Ysz * Zsz;
0793 node2vox = sparse(voxels(1:3:end,:),kron((1:nVox)',ones(2,4)),1);
0794
0795
0796
0797 function faces = mk_faces(voxels,opt)
0798 Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0799
0800 facerow = [nonzeros(bsxfun(@plus,(1:3)',0:6:(Xsz-1)*6)); (Xsz-1)*6+4];
0801 faceplane = bsxfun(@plus,facerow,0:6*Xsz:6*Xsz*(Ysz-1));
0802
0803 cap = kron(ones(3,1),faceplane(2:3:end,end)')+3;
0804 faceplane = [faceplane,[ cap(:); cap(end)]];
0805 faces = bsxfun(@plus, faceplane(:), 0:6*Xsz*Ysz:6*Xsz*Ysz*(Zsz-1));
0806
0807 cap = reshape(kron(ones(3,1), 3:3:3*Xsz*Ysz'),3*Xsz,[]);
0808 cap = bsxfun(@plus, cap, 0:Ysz-1);
0809 cap(end+1,:)= cap(end,:);
0810 faces = [faces(:); faces(cap(:),end)+3];
0811 faces = voxels(faces,:);
0812
0813
0814
0815 function edges = mk_edges(voxels, opt)
0816 Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0817
0818 edgerow = voxels([nonzeros(bsxfun(@plus,(1:3)',0:6:(Xsz-1)*6)); (Xsz-1)*6+4],1:2);
0819 edgerow(end+1,:) = edgerow(end,:);
0820 edgerow(end+1,:) = voxels((Xsz-1)*6+4, [1 4]);
0821 edgeplane = repmat(edgerow,Ysz+1,1) + kron((Xsz+1)*(0:Ysz)', ones(size(edgerow)));
0822
0823 edgeplane((size(edgerow,1)*Ysz +3):3:end,:) = edgeplane((size(edgerow,1)*Ysz +2):3:end,:);
0824
0825 up = (Xsz+1)*(Ysz+1);
0826 edges = repmat(edgeplane,Zsz+1,1) + kron((0:Zsz)'*up, ones(size(edgeplane)));
0827
0828
0829 start = (size(edgeplane,1)*Zsz);
0830 edges( (start +1):3:end,:) = edges( (start +2):3:end,:);
0831 start = (size(edgeplane,1)*Zsz) + 3*Xsz;
0832 step = size(edgerow,1);
0833 edges( (start +1):step:end,:) = edges( (start +3):step:end,:);
0834 edges( (start +2):step:end,:) = edges( (start +3):step:end,:);
0835 edges( end-2:end,:) = edges(end-5:end-3,:);
0836
0837
0838
0839 function vox2face = mk_vox2face(opt)
0840 Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0841 xstep = opt.xstep; ystep = opt.ystep; zstep = opt.zstep;
0842
0843 vox = [1 2 3 1+xstep 2+ystep 3+zstep];
0844 voxrow = bsxfun(@plus,(0:Xsz-1)'*xstep,vox)';
0845 voxplane = bsxfun(@plus,(0:Ysz-1)*ystep,voxrow(:));
0846 voxvol = bsxfun(@plus,(0:Zsz-1)*zstep,voxplane(:));
0847 voxvol = reshape(voxvol, 6, [])';
0848 I = kron((1:size(voxvol,1))',ones(1,6));
0849 nFaces = Zsz*(Ysz+1)*(3*Xsz+1) + Ysz*(3*Xsz+1);
0850 nVox = Xsz*Ysz*Zsz;
0851 vox2face = sparse(I,voxvol,ones(size(I)),nVox,nFaces);
0852
0853
0854
0855 function [vox2edge, vol] = mk_vox2edge(m,opt)
0856 Xsz = opt.Xsz; Ysz = opt.Ysz; Zsz = opt.Zsz;
0857 xstep = opt.xstep; ystep = xstep*(Xsz+1); zstep = ystep*(Ysz+1);
0858
0859 vox = [1 2 3 1+xstep 3+xstep 1+ystep 2+ystep 1+ystep+xstep ...
0860 2+zstep 3+zstep 3+zstep+xstep 2+zstep+ystep];
0861 voxrow = bsxfun(@plus,(0:Xsz-1)'*xstep,vox)';
0862 voxplane = bsxfun(@plus,(0:Ysz-1)*ystep,voxrow(:));
0863 voxvol = bsxfun(@plus,(0:Zsz-1)*zstep,voxplane(:));
0864 voxvol = reshape(voxvol, 12, [])';
0865 I = kron((1:size(voxvol,1))',ones(1,12));
0866 nEdges = 3*(Zsz+1)*(Ysz+1)*(Xsz+1);
0867 nVox = Xsz*Ysz*Zsz;
0868 vox2edge = sparse(I,voxvol,ones(size(I)),nVox,nEdges);
0869 vol = m.edge_length(voxvol(:,1)) ...
0870 .* m.edge_length(voxvol(:,2)) ...
0871 .* m.edge_length(voxvol(:,3));
0872
0873
0874
0875
0876 function show_voxels(rmdl,voxels)
0877 mdl=rmdl;
0878 mdl.elems = voxels;
0879 mdl.boundary = mdl.elems;
0880 clf
0881 show_fem(mdl,[0 0 1]);
0882
0883
0884
0885 function test_faces(rmdl, faces, opt)
0886 mdl = rmdl;
0887 mdl.elems = faces;
0888 mdl.boundary = mdl.elems;
0889 img = mk_image(mdl,0);
0890 img.elem_data(opt.xplane + (randi(opt.Xsz+1)-1)*opt.xstep) = 1;
0891 img.elem_data(opt.yplane + (randi(opt.Ysz+1)-1)*opt.ystep) = 2;
0892 img.elem_data(opt.zplane + (randi(opt.Zsz+1)-1)*opt.zstep) = 3;
0893 show_fem(img,[0 0 1]);
0894
0895
0896
0897 function[fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt)
0898 ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0899 rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0900 fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0901 if isfield(opt,'do_not_scale') && opt.do_not_scale
0902 return
0903 end
0904 maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0905 scale = 1/maxnode;
0906 rmdl.nodes = scale*rmdl.nodes;
0907 fmdl.nodes = scale*fmdl.nodes;
0908 eidors_msg('@@ models scaled by %g', scale,2);
0909
0910
0911
0912 function opt = parse_opts(fmdl,rmdl, opt)
0913
0914 opt.xvec = unique(rmdl.nodes(:,1));
0915 opt.yvec = unique(rmdl.nodes(:,2));
0916 opt.zvec = unique(rmdl.nodes(:,3));
0917 opt.Xsz = numel(opt.xvec)-1;
0918 opt.Ysz = numel(opt.yvec)-1;
0919 opt.Zsz = numel(opt.zvec)-1;
0920 opt.xstep = 3;
0921 opt.ystep = opt.xstep*opt.Xsz+1;
0922 opt.zstep = opt.ystep*(opt.Ysz+1);
0923 xrow = 1 + (0:opt.Ysz-1)*opt.ystep;
0924 opt.xplane = bsxfun(@plus, (0:opt.Zsz-1)'*opt.zstep,xrow);
0925 yrow = 2 + (0:opt.Xsz-1)*opt.xstep;
0926 opt.yplane = bsxfun(@plus, (0:opt.Zsz-1)'*opt.zstep,yrow);
0927 zrow = 3 + (0:opt.Xsz-1)*opt.xstep;
0928 opt.zplane = bsxfun(@plus, (0:opt.Ysz-1)'*opt.ystep,zrow);
0929 opt.nVox = opt.Xsz*opt.Ysz*opt.Zsz;
0930 opt.nTet = size(fmdl.elems,1);
0931
0932 if ~isfield(opt, 'tol_node2tet');
0933 opt.tol_node2tet = eps;
0934 end
0935 if ~isfield(opt, 'tol_edge2edge')
0936 opt.tol_edge2edge = 6*sqrt(3)*eps(min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))));
0937 end
0938 if ~isfield(opt, 'tol_edge2tri')
0939 opt.tol_edge2tri = eps;
0940 end
0941 if ~isfield(opt, 'save_memory')
0942 opt.save_memory = 0;
0943 end
0944 eidors_msg('@@ node2tet tolerance = %g', opt.tol_node2tet,2);
0945 eidors_msg('@@ edge2edge tolerance = %g', opt.tol_edge2edge,2);
0946 eidors_msg('@@ edge2tri tolerance = %g', opt.tol_edge2tri,2);
0947
0948
0949
0950 function logmsg(varargin)
0951 if eidors_msg('log_level') >= 2
0952 fprintf(varargin{:});
0953 end
0954
0955
0956
0957
0958
0959 function do_unit_test
0960 mk_grid_c2f('save_memory',0); do_small_test;
0961 mk_grid_c2f('save_memory',1); do_small_test;
0962 mk_grid_c2f('save_memory',10); do_small_test;
0963 return
0964 do_realistic_test;
0965 figure
0966 do_case_tests;
0967 do_edge2edge_timing_test;
0968
0969
0970
0971 function do_case_tests
0972 ll = eidors_msg('log_level');
0973 eidors_msg('log_level',1);
0974 vox = mk_grid_model([],0:1,0:1,0:1);
0975 tet.type = 'fwd_model';
0976 tet.elems = [1 2 3 4];
0977
0978 X = 4; Y = 6;
0979 for i = 1:30
0980 tet.nodes = [0 0 0; 0 1 0; 1 0 0; 0 0 1];
0981 fprintf('%d\n',i);
0982 switch i
0983 case 1
0984 txt = 'Nothing in common';
0985 tet.nodes = tet.nodes + 2;
0986 subplot(X,Y,i), show_test(vox,tet);
0987 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
0988 unit_test_cmp(txt,size(intpts,1),0,0);
0989 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
0990 unit_test_cmp(txt,size(intpts,1),0,0);
0991 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
0992 unit_test_cmp(txt,size(intpts,1),0,0);
0993 mk_grid_c2f(tet,vox);
0994
0995 case 2
0996 tet.nodes = tet.nodes + 1;
0997 subplot(X,Y,i), show_test(vox,tet);
0998 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
0999 unit_test_cmp('Comm node tedge2rec',size(intpts,1),0,0);
1000 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1001 unit_test_cmp('Comm node vedge2tri',size(intpts,1),0,0);
1002 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1003 unit_test_cmp('Comm node edge2edge',size(intpts,1),0,0);
1004 [fnode2vox] = test_tnode_in_vox(vox,tet);
1005 unit_test_cmp('Comm node tnode2vox',fnode2vox,[1;0;0;0],0);
1006 rnode2tet = test_vnode_in_tet(vox,tet);
1007 unit_test_cmp('Comm node vnode2tet',nnz(rnode2tet),0,0);
1008 mk_grid_c2f(tet,vox);
1009
1010 case 3
1011 txt = 'tet_edge v vox_node';
1012 tet.nodes(:,1:2) = tet.nodes(:,1:2) - 0.5;
1013 subplot(X,Y,i), show_test(vox,tet);
1014 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1015 unit_test_cmp(txt,size(intpts,1),0,0);
1016 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1017 unit_test_cmp(txt,size(intpts,1),0,0);
1018 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1019 unit_test_cmp(txt,size(intpts,1),0,0);
1020 [fnode2vox] = test_tnode_in_vox(vox,tet);
1021 unit_test_cmp(txt,nnz(fnode2vox),0,0);
1022 rnode2tet = test_vnode_in_tet(vox,tet);
1023 res = zeros(8,1); res(1) = 1;
1024 unit_test_cmp(txt,rnode2tet,res,0);
1025 mk_grid_c2f(tet,vox);
1026
1027 case 4
1028 txt = 'tet_edge v vox_edge';
1029 tet.nodes(:,1:2) = tet.nodes(:,1:2) - 0.5;
1030 tet.nodes(:,3) = tet.nodes(:,3) + 0.5;
1031 subplot(X,Y,i), show_test(vox,tet);
1032 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1033 unit_test_cmp(txt,size(intpts,1),0,0);
1034 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1035 unit_test_cmp(txt,size(intpts,1),0,0);
1036 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1037 unit_test_cmp(txt,size(intpts,1),1,0);
1038 [fnode2vox] = test_tnode_in_vox(vox,tet);
1039 unit_test_cmp(txt,nnz(fnode2vox),0,0);
1040 rnode2tet = test_vnode_in_tet(vox,tet);
1041 unit_test_cmp(txt,nnz(rnode2tet),0,0);
1042 mk_grid_c2f(tet,vox);
1043
1044 case 5
1045 txt = 'tet_edge on vox_face';
1046 tet.nodes(:,1:2) = tet.nodes(:,1:2) - 0.3;
1047 subplot(X,Y,i), show_test(vox,tet);
1048 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1049 unit_test_cmp(txt,size(intpts,1),0,0);
1050 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1051 unit_test_cmp(txt,size(intpts,1),1,0);
1052 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1053 unit_test_cmp(txt,size(intpts,1),2,0);
1054 [fnode2vox] = test_tnode_in_vox(vox,tet);
1055 unit_test_cmp(txt,nnz(fnode2vox),0,0);
1056 rnode2tet = test_vnode_in_tet(vox,tet);
1057 unit_test_cmp(txt,nnz(rnode2tet),1,0);
1058 mk_grid_c2f(tet,vox);
1059
1060 case 6
1061 txt = 'vox_node on tet_face';
1062 tet.nodes(:,1:2) = tet.nodes(:,1:2) - 0.3;
1063 tet.nodes(:,3) = tet.nodes(:,3) -0.4;
1064 subplot(X,Y,i), show_test(vox,tet);
1065 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1066 unit_test_cmp(txt,size(intpts,1),0,0);
1067 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1068 unit_test_cmp(txt,size(intpts,1),0,0);
1069 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1070 unit_test_cmp(txt,size(intpts,1),0,0);
1071 [fnode2vox] = test_tnode_in_vox(vox,tet);
1072 unit_test_cmp(txt,nnz(fnode2vox),0,0);
1073 rnode2tet = test_vnode_in_tet(vox,tet);
1074 unit_test_cmp(txt,nnz(rnode2tet),1,0);
1075 mk_grid_c2f(tet,vox);
1076
1077 case 7
1078 txt = 'tet_node on vox_face';
1079 tet.nodes(:,1) = tet.nodes(:,1) - 1;
1080 tet.nodes(:,2:3) = tet.nodes(:,2:3) + .5;
1081 subplot(X,Y,i), show_test(vox,tet);
1082 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1083 unit_test_cmp(txt,size(intpts,1),0,0);
1084 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1085 unit_test_cmp(txt,size(intpts,1),0,0);
1086 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1087 unit_test_cmp(txt,size(intpts,1),0,0);
1088 [fnode2vox] = test_tnode_in_vox(vox,tet);
1089 unit_test_cmp(txt,nnz(fnode2vox),1,0);
1090 rnode2tet = test_vnode_in_tet(vox,tet);
1091 unit_test_cmp(txt,nnz(rnode2tet),0,0);
1092 mk_grid_c2f(tet,vox);
1093
1094
1095 case 8
1096 txt = 'tet_node on vox_edge';
1097 tet.nodes(:,1) = tet.nodes(:,1) - 1;
1098 tet.nodes(:,2) = tet.nodes(:,2) + .5;
1099 subplot(X,Y,i), show_test(vox,tet);
1100 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1101 unit_test_cmp(txt,size(intpts,1),0,0);
1102 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1103 unit_test_cmp(txt,size(intpts,1),0,0);
1104 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1105 unit_test_cmp(txt,size(intpts,1),0,0);
1106 [fnode2vox] = test_tnode_in_vox(vox,tet);
1107 unit_test_cmp(txt,nnz(fnode2vox),1,0);
1108 rnode2tet = test_vnode_in_tet(vox,tet);
1109 unit_test_cmp(txt,nnz(rnode2tet),0,0);
1110 mk_grid_c2f(tet,vox);
1111
1112 case 9
1113 txt = 'all nodes common';
1114 subplot(X,Y,i), show_test(vox,tet);
1115 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1116 unit_test_cmp(txt,size(intpts,1),0,0);
1117 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1118 unit_test_cmp(txt,size(intpts,1),0,0);
1119 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1120 unit_test_cmp(txt,size(intpts,1),0,0);
1121 [fnode2vox] = test_tnode_in_vox(vox,tet);
1122 unit_test_cmp(txt,nnz(fnode2vox),4,0);
1123 rnode2tet = test_vnode_in_tet(vox,tet);
1124 unit_test_cmp(txt,nnz(rnode2tet),0,0);
1125 mk_grid_c2f(tet,vox);
1126
1127 case 10
1128 txt = 'tet in vox';
1129 tet.nodes = .25 + .5*tet.nodes;
1130 subplot(X,Y,i), show_test(vox,tet);
1131 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1132 unit_test_cmp(txt,size(intpts,1),0,0);
1133 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1134 unit_test_cmp(txt,size(intpts,1),0,0);
1135 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1136 unit_test_cmp(txt,size(intpts,1),0,0);
1137 [fnode2vox] = test_tnode_in_vox(vox,tet);
1138 unit_test_cmp(txt,nnz(fnode2vox),4,0);
1139 rnode2tet = test_vnode_in_tet(vox,tet);
1140 unit_test_cmp(txt,nnz(rnode2tet),0,0);
1141 mk_grid_c2f(tet,vox);
1142 c2f = mk_grid_c2f(tet,vox);
1143 max_vox = max(vox.nodes);
1144 min_vox = min(vox.nodes);
1145 edges = max_vox-min_vox;
1146 vox_vol = edges(:,1) .* edges(:,2) .* edges(:,3);
1147 tet_vol = get_elem_volume(tet);
1148 unit_test_cmp(txt,c2f, 1, 0);
1149
1150 case 11
1151 txt = 'vox in tet';
1152 tet.nodes = 4*tet.nodes - .25;
1153 subplot(X,Y,i), show_test(vox,tet);
1154 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1155 unit_test_cmp(txt,size(intpts,1),0,0);
1156 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1157 unit_test_cmp(txt,size(intpts,1),0,0);
1158 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1159 unit_test_cmp(txt,size(intpts,1),0,0);
1160 [fnode2vox] = test_tnode_in_vox(vox,tet);
1161 unit_test_cmp(txt,nnz(fnode2vox),0,0);
1162 rnode2tet = test_vnode_in_tet(vox,tet);
1163 unit_test_cmp(txt,nnz(rnode2tet),8,0);
1164 c2f = mk_grid_c2f(tet,vox);
1165 max_vox = max(vox.nodes);
1166 min_vox = min(vox.nodes);
1167 edges = max_vox-min_vox;
1168 vox_vol = edges(:,1) .* edges(:,2) .* edges(:,3);
1169 tet_vol = get_elem_volume(tet);
1170 unit_test_cmp(txt,c2f, vox_vol / tet_vol, 0);
1171
1172
1173 case 12
1174 txt = 'tet_edge v. vox_face';
1175 tet.nodes(:,1:2) = tet.nodes(:,1:2) - 0.3;
1176 tet.nodes(:,3) = tet.nodes(:,3) + .4;
1177 subplot(X,Y,i), show_test(vox,tet);
1178 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1179 unit_test_cmp(txt,size(intpts,1),2,0);
1180 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1181 unit_test_cmp(txt,size(intpts,1),2,0);
1182 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1183 unit_test_cmp(txt,size(intpts,1),0,0);
1184 [fnode2vox] = test_tnode_in_vox(vox,tet);
1185 unit_test_cmp(txt,nnz(fnode2vox),0,0);
1186 rnode2tet = test_vnode_in_tet(vox,tet);
1187 unit_test_cmp(txt,nnz(rnode2tet),0,0);
1188 mk_grid_c2f(tet,vox);
1189
1190
1191 case 13
1192 txt = 'everything';
1193 tet.nodes = tet.nodes + 0.7;
1194 subplot(X,Y,i), show_test(vox,tet);
1195 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1196 unit_test_cmp(txt,size(intpts,1),3,0);
1197 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1198 unit_test_cmp(txt,size(intpts,1),3,0);
1199 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1200 unit_test_cmp(txt,size(intpts,1),0,0);
1201 [fnode2vox] = test_tnode_in_vox(vox,tet);
1202 unit_test_cmp(txt,nnz(fnode2vox),1,0);
1203 rnode2tet = test_vnode_in_tet(vox,tet);
1204 unit_test_cmp(txt,nnz(rnode2tet),1,0);
1205 mk_grid_c2f(tet,vox);
1206
1207
1208 case 14
1209 txt = 'DG1: Common edge';
1210 tet.nodes(:,1:2) = tet.nodes(:,1:2) + 1;
1211 subplot(X,Y,i), show_test(vox,tet);
1212 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1213 unit_test_cmp(txt,size(intpts,1),0,0);
1214 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1215 unit_test_cmp(txt,size(intpts,1),0,0);
1216 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1217 unit_test_cmp(txt,size(intpts,1),0,0);
1218 [fnode2vox] = test_tnode_in_vox(vox,tet);
1219 unit_test_cmp(txt,fnode2vox,[1;0;0;1],0);
1220 rnode2tet = test_vnode_in_tet(vox,tet);
1221 unit_test_cmp(txt,nnz(rnode2tet),0,0);
1222 mk_grid_c2f(tet,vox);
1223
1224 case 15
1225 txt = 'DG2: Common edge';
1226 tet.nodes(:,1:2) = tet.nodes(:,1:2) + 1;
1227 tet.nodes(:,3) = tet.nodes(:,3) + 0.5;
1228 subplot(X,Y,i), show_test(vox,tet);
1229 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1230 unit_test_cmp(txt,size(intpts,1),0,0);
1231 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1232 unit_test_cmp(txt,size(intpts,1),0,0);
1233 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1234 unit_test_cmp(txt,size(intpts,1),0,0);
1235 [fnode2vox] = test_tnode_in_vox(vox,tet);
1236 unit_test_cmp(txt,nnz(fnode2vox),1,0);
1237 rnode2tet = test_vnode_in_tet(vox,tet);
1238 unit_test_cmp(txt,nnz(rnode2tet),1,0);
1239 mk_grid_c2f(tet,vox);
1240
1241 case 16
1242 txt = 'DG3: Common edge';
1243 tet.nodes(:,1:2) = tet.nodes(:,1:2) + 1;
1244 z = tet.nodes(:,3) == 0;
1245 tet.nodes(z,3) = -.5;
1246 tet.nodes(~z,3) = 1.5;
1247 subplot(X,Y,i), show_test(vox,tet);
1248 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1249 unit_test_cmp(txt,size(intpts,1),0,0);
1250 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1251 unit_test_cmp(txt,size(intpts,1),0,0);
1252 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1253 unit_test_cmp(txt,size(intpts,1),0,0);
1254 [fnode2vox] = test_tnode_in_vox(vox,tet);
1255 unit_test_cmp(txt,nnz(fnode2vox),0,0);
1256 rnode2tet = test_vnode_in_tet(vox,tet);
1257 unit_test_cmp(txt,nnz(rnode2tet),2,0);
1258 mk_grid_c2f(tet,vox);
1259
1260 case 17
1261 txt = 'DG4: edge on face';
1262 tet.nodes = [0 0 0; 1 .5 0; 0 1 0; 1 .5 1];
1263 tet.nodes(:,1) = tet.nodes(:,1) - 1;
1264 subplot(X,Y,i), show_test(vox,tet);
1265 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1266 unit_test_cmp(txt,size(intpts,1),0,0);
1267 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1268 unit_test_cmp(txt,size(intpts,1),0,0);
1269 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1270 unit_test_cmp(txt,size(intpts,1),0,0);
1271 [fnode2vox] = test_tnode_in_vox(vox,tet);
1272 unit_test_cmp(txt,nnz(fnode2vox),2,0);
1273 rnode2tet = test_vnode_in_tet(vox,tet);
1274 unit_test_cmp(txt,nnz(rnode2tet),0,0);
1275 mk_grid_c2f(tet,vox);
1276
1277 case 18
1278 txt = 'DG5: edge2edge only';
1279 tet.nodes = [0 0 0; 1 .5 0; 0 1 0; 1 .5 1];
1280 tet.nodes(:,1) = tet.nodes(:,1) - 1;
1281 z = tet.nodes(:,3) == 0;
1282 tet.nodes(z,3) = -.5;
1283 tet.nodes(~z,3) = 1.5;
1284 subplot(X,Y,i), show_test(vox,tet);
1285 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1286 unit_test_cmp(txt,size(intpts,1),0,0);
1287 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1288 unit_test_cmp(txt,size(intpts,1),0,0);
1289 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1290 unit_test_cmp(txt,size(intpts,1),2,0);
1291 [fnode2vox] = test_tnode_in_vox(vox,tet);
1292 unit_test_cmp(txt,nnz(fnode2vox),0,0);
1293 rnode2tet = test_vnode_in_tet(vox,tet);
1294 unit_test_cmp(txt,nnz(rnode2tet),0,0);
1295 mk_grid_c2f(tet,vox);
1296
1297 case 19
1298 txt = 'DG6: edge on face';
1299 tet.nodes = [0 0 0; 1 .5 0; 0 1 0; 1 .5 1];
1300 tet.nodes(:,1) = tet.nodes(:,1) - 1;
1301 tet.nodes(:,3) = tet.nodes(:,3) - .5;
1302 subplot(X,Y,i), show_test(vox,tet);
1303 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1304 unit_test_cmp(txt,size(intpts,1),0,0);
1305 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1306 unit_test_cmp(txt,size(intpts,1),0,0);
1307 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1308 unit_test_cmp(txt,size(intpts,1),1,0);
1309 [fnode2vox] = test_tnode_in_vox(vox,tet);
1310 unit_test_cmp(txt,nnz(fnode2vox),1,0);
1311 rnode2tet = test_vnode_in_tet(vox,tet);
1312 unit_test_cmp(txt,nnz(rnode2tet),0,0);
1313 mk_grid_c2f(tet,vox);
1314
1315 case 20
1316 txt = 'DG7: face on face';
1317 tet.nodes(:,1) = tet.nodes(:,1) + 1;
1318 tet.nodes(:,2) = tet.nodes(:,2) + 0.5;
1319 tet.nodes(:,3) = tet.nodes(:,3) + 0.5;
1320 subplot(X,Y,i), show_test(vox,tet);
1321 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1322 unit_test_cmp(txt,size(intpts,1),0,0);
1323 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1324 unit_test_cmp(txt,size(intpts,1),0,0);
1325 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1326 unit_test_cmp(txt,size(unique(intpts,'rows'),1),2,0);
1327 [fnode2vox] = test_tnode_in_vox(vox,tet);
1328 unit_test_cmp(txt,nnz(fnode2vox),1,0);
1329 rnode2tet = test_vnode_in_tet(vox,tet);
1330 unit_test_cmp(txt,nnz(rnode2tet),1,0);
1331 mk_grid_c2f(tet,vox);
1332
1333 case 21
1334 txt = 'DG8: face on face';
1335 tet.nodes(:,1) = tet.nodes(:,1) + 1;
1336 tet.nodes(:,2) = tet.nodes(:,2) + 0.4;
1337 tet.nodes(:,3) = tet.nodes(:,3) - 0.6;
1338 subplot(X,Y,i), show_test(vox,tet);
1339 [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(vox,tet);
1340 unit_test_cmp(txt,size(intpts,1),0,0);
1341 [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(vox,tet);
1342 unit_test_cmp(txt,size(intpts,1),0,0);
1343 [intpts, face2edge, face2intpt, edge2intpt] = test_tedge2vedge_intersections(vox,tet);
1344 unit_test_cmp(txt,size(unique(intpts,'rows'),1),2,0);
1345 [fnode2vox] = test_tnode_in_vox(vox,tet);
1346 unit_test_cmp(txt,nnz(fnode2vox),1,0);
1347 rnode2tet = test_vnode_in_tet(vox,tet);
1348 unit_test_cmp(txt,nnz(rnode2tet),0,0);
1349 mk_grid_c2f(tet,vox);
1350
1351 otherwise
1352 break;
1353 end
1354 end
1355 eidors_msg('log_level',ll);
1356
1357 function do_small_test
1358 fmdl= ng_mk_cyl_models([2,2,.1],[16,1],[.1,0,.025]);
1359 xvec = [-1.5:1:1.5];
1360 yvec = [-1.6:.8:1.6];
1361 zvec = -1:1:2;
1362 rmdl = mk_grid_model([],xvec,yvec,zvec);
1363
1364 eidors_cache clear
1365 tic
1366 c2f_a = mk_grid_c2f(fmdl, rmdl);
1367 toc
1368 eidors_cache clear
1369 tic
1370 opt.save_memory = 2;
1371 c2f_b = mk_grid_c2f(fmdl, rmdl, opt);
1372 toc
1373 eidors_cache clear
1374 tic
1375 c2f_c = mk_approx_c2f(fmdl, rmdl);
1376 toc
1377 assert(any(c2f_a(:) ~= c2f_b(:)) == 0);
1378
1379
1380
1381
1382 function do_realistic_test
1383 fmdl= ng_mk_cyl_models([2,2,.1],[16,1],[.1,0,.025]);
1384 xvec = [-1.5 -.5:.2:.5 1.5];
1385 yvec = [-1.6 -1:.2:1 1.6];
1386 zvec = 0:.25:2;
1387 rmdl = mk_grid_model([],xvec,yvec,zvec);
1388
1389 tic
1390 opt.save_memory = 0;
1391 c2f_a = mk_grid_c2f(fmdl, rmdl,opt);
1392 t = toc;
1393 fprintf('Analytic: t=%f s\n',t);
1394
1395 tic
1396 c2f_n = mk_approx_c2f(fmdl,rmdl);
1397 t = toc;
1398 fprintf('Approximate: t=%f s\n',t);
1399
1400
1401 tetvol = get_elem_volume(fmdl);
1402 opt = parse_opts(fmdl,rmdl);
1403 m = prepare_vox_mdl(rmdl,opt);
1404
1405
1406 f2c_a = bsxfun(@times, c2f_a, tetvol);
1407 f2c_a = bsxfun(@rdivide,f2c_a', m.volume);
1408 img = mk_image(rmdl,0);
1409 img.elem_data = f2c_a*ones(size(fmdl.elems,1),1);
1410 subplot(132)
1411 show_fem(img);
1412
1413 f2c_n = bsxfun(@times, c2f_n, tetvol);
1414 f2c_n = bsxfun(@rdivide,f2c_n', m.volume);
1415 img = mk_image(rmdl,0);
1416 img.elem_data = f2c_n*ones(size(fmdl.elems,1),1);
1417 subplot(133)
1418 show_fem(img);
1419 subplot(131);
1420 h = show_fem(fmdl);
1421 set(h,'LineWidth',0.1)
1422 hold on
1423 h = show_fem(rmdl);
1424 set(h,'EdgeColor','b','LineWidth',2);
1425 hold off
1426
1427 function do_edge2edge_timing_test
1428 fmdl= ng_mk_cyl_models([2,2],[8,1],[.15,0,.05]);
1429 xvec = linspace(-2,2,33);
1430 yvec = linspace(-2,2,33);
1431 zvec = 0:.5:2;
1432 rmdl = mk_grid_model([],xvec,yvec,zvec);
1433 tic
1434 test_tedge2vedge_intersections(rmdl,fmdl);
1435 toc
1436
1437 function rnode2tet = test_vnode_in_tet(rmdl,fmdl)
1438 opt = parse_opts(fmdl,rmdl);
1439 fmdl = prepare_fmdl(fmdl);
1440 rnode2tet = get_nodes_in_tets(fmdl,rmdl, opt);
1441
1442 function [fnode2vox] = test_tnode_in_vox(rmdl,fmdl)
1443 fmdl = prepare_fmdl(fmdl);
1444 [fnode2vox] = get_nodes_in_voxels(fmdl,rmdl);
1445
1446 function [intpts, face2edge, face2intpt, edge2intpt] = test_rec_tedge_intersections(rmdl,fmdl)
1447 opt = parse_opts(fmdl,rmdl);
1448 m = prepare_vox_mdl(rmdl,opt);
1449 fmdl = prepare_fmdl(fmdl);
1450 [intpts, face2edge, face2intpt, edge2intpt] = get_voxel_intersection_points(fmdl,m.faces,opt);
1451
1452 function [intpts, face2edge, face2intpt, edge2intpt] = test_tet_vedge_intersections(rmdl,fmdl)
1453 opt = parse_opts(fmdl,rmdl);
1454 m = prepare_vox_mdl(rmdl,opt);
1455 fmdl = prepare_fmdl(fmdl);
1456 [intpts, face2edge, face2intpt, edge2intpt] = get_tet_intersection_points(fmdl,m,opt);
1457
1458 function [intpts, tedge2vedge, tedge2intpt, vedge2intpt] = test_tedge2vedge_intersections(rmdl,fmdl)
1459 opt = parse_opts(fmdl,rmdl);
1460 m = prepare_vox_mdl(rmdl,opt);
1461 fmdl = prepare_fmdl(fmdl);
1462 [intpts, tedge2vedge, tedge2intpt, vedge2intpt] = find_edge2edge_intersections(fmdl.edges,fmdl.nodes,m.edges,rmdl.nodes, opt.tol_node2tet);
1463
1464 function show_test(vox,tet)
1465 show_fem(vox);
1466 hold on
1467 h = show_fem(tet);
1468 set(h, 'EdgeColor','b');
1469 hold off
1470 axis auto