Description of mk_tet_c2f

0001 function [c2f] = mk_tet_c2f(fmdl, rmdl, opt)
0002 %MK_TET_C2F - calculate a coarse2fine mapping for two tet-based models.
0003 % C2F = MK_TET_C2F(FMDL,RMDL) returns in C2F the fraction of volume of
0004 % each element of the fine model FMDL contained in each element of
0005 % the coarse model RMDL.
0006 % Uses CONVHULLN to calculate the volume defined by a set of intersection
0007 % points between individual tet elements.
0008 %
0009 % C2F = MK_TET_C2F(FMDL,RMDL,OPT) allows specifying options.
0010 %
0011 % Inputs:
0012 %   FMDL - a (fine) EIDORS (tet-based) forward model
0013 %   RMDL - a (course) EIDORS (tet-based) forward model
0014 %   OPT  - an option structure with the following fields and defaults:
0015 %      .do_not_scale  - set to true to prevent scaling the models to unit
0016 %                       cube before any calculations, including thresholds.
0017 %                       Default: false
0018 %      .tol_node2tet  - tolerance for determinant <= 0 in testing for
0019 %                       points inside tets. Default: eps
0020 %      .tol_edge2edge - maximum distance between "intersecting" edges
0021 %                       Default: 6*sqrt(3)*eps(a), where a is
0022 %                       min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))
0023 %      .tol_edge2tri  - minimum value of a barycentric coordinate to
0024 %                       decide a point is lying inside a triangle and not
0025 %                       on its edge. Default: eps
0026 %
0027 % NOTE that for grid-based models, such as returned by MK_GRID_MODEL or
0028 % MK_VOXEL_VOLUME, MK_GRID_C2F is much faster.
0029 %
0030 % Set eidors_msg 'log level' < 2 to supress output to command line.
0031 %
0032 % Examples:
0033 %     fmdl = ng_mk_cyl_models([2,2,.2],[],[]);
0034 %     rmdl = ng_mk_cyl_models([2,2],[],[]);
0035 %     c2f  = mk_tet_c2f(fmdl,rmdl);
0036 %     h = show_fem(fmdl); set(h,'LineWidth',0.1)
0037 %     hold on
0038 %     h = show_fem(rmdl); set(h,'EdgeColor','b','LineWidth',2);
0039 %     hold off
0040 %
0041 % See also MK_GRID_C2F, FIND_EDGE2EDGE_INTERSECTIONS, CONVHULLN
0042 %     MK_COARSE_FINE_MAPPING, MK_APPROX_C2F, POINT_IN_TRIANGLE, EIDORS_MSG
0043 
0044 
0045 % (C) 2015 Bartlomiej Grychtol - all rights reserved by Swisstom AG
0046 % License: GPL version 2 or 3
0047 % $Id: mk_tet_c2f.m 5889 2018-12-23 19:57:54Z aadler $
0048 
0049 % >> SWISSTOM CONTRIBUTION <<
0050 
0051 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'), do_unit_test; return; end
0052 if nargin < 3
0053    opt = struct();
0054 end
0055 
0056 f_elems = size(fmdl.elems,1);
0057 r_elems = size(rmdl.elems,1);
0058 
0059 c2f = sparse(f_elems,r_elems);
0060 [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl);
0061 
0062 if ~(any(fmdl_idx) && any(rmdl_idx))
0063    eidors_msg('@@: models do not overlap, returning all-zeros');
0064    return
0065 end
0066 
0067 [fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt);
0068 
0069 opt = parse_opts(fmdl,rmdl, opt);
0070 
0071 
0072 copt.fstr = 'mk_tet_c2f';
0073 
0074 c2f(fmdl_idx,rmdl_idx) = eidors_cache(@do_mk_tet_c2f,{fmdl,rmdl,opt},copt);
0075 
0076 
0077 function c2f = do_mk_tet_c2f(fmdl,rmdl,opt)
0078    DEBUG = eidors_debug('query','mk_tet_c2f');
0079    
0080    c2f = sparse(0,0);
0081    progress_msg('Prepare fine model...');
0082    fmdl = prepare_tet_mdl(fmdl);
0083    progress_msg(Inf);
0084    
0085    progress_msg('Prepare course model...');
0086    rmdl = prepare_tet_mdl(rmdl);
0087    progress_msg(Inf);
0088    
0089    progress_msg('Find c_edge2f_face intersections...')
0090    [intpts1, fface2redge, fface2intpt1, redge2intpt1] = ...
0091       edge2face_intersections(fmdl,rmdl,opt);
0092    progress_msg(sprintf('Found %d', size(intpts1,1)), Inf);
0093 
0094    progress_msg('Find f_edge2c_face intersections...')
0095    [intpts2, rface2fedge, rface2intpt2, fedge2intpt2] = ...
0096       edge2face_intersections(rmdl,fmdl,opt);
0097    progress_msg(sprintf('Found %d', size(intpts2,1)), Inf);
0098 
0099    pmopt.final_msg = 'none';
0100    progress_msg('Find edge2edge intersections...',-1,pmopt)
0101    [intpts3, fedge2redge, fedge2intpt3, redge2intpt3] = ...
0102       find_edge2edge_intersections(fmdl.edges, fmdl.nodes, ...
0103                                    rmdl.edges, rmdl.nodes, ...
0104                                    opt.tol_edge2edge);
0105    progress_msg(sprintf('Found %d',size(intpts3,1)),Inf);
0106 
0107    progress_msg('Find c_nodes in f_tets...',pmopt);
0108    rnode2ftet = get_nodes_in_tets(fmdl,rmdl.nodes, opt);
0109    progress_msg(sprintf('Found %d', nnz(rnode2ftet)), Inf);
0110    
0111    
0112    progress_msg('Find c_elems in f_elems...',pmopt)
0113    rtet_in_ftet = (double(rmdl.node2elem') * rnode2ftet) == 4;
0114    progress_msg(sprintf('Found %d',nnz(rtet_in_ftet)), Inf);
0115    
0116    progress_msg('Find f_nodes in c_tets...',pmopt);
0117    fnode2rtet = get_nodes_in_tets(rmdl,fmdl.nodes, opt);
0118    progress_msg(sprintf('Found %d', nnz(fnode2rtet)), Inf);
0119 
0120    progress_msg('Find f_elems in c_elems...',pmopt)
0121    ftet_in_rtet = (double(fmdl.node2elem') * fnode2rtet) == 4;
0122    progress_msg(sprintf('Found %d',nnz(ftet_in_rtet)), Inf);
0123    
0124    progress_msg('Find total intersections...',pmopt);
0125    e2e = double(rmdl.edge2elem');
0126    rtet2ftet =  double(rmdl.elem2face) * (rface2fedge>0) * fmdl.edge2elem ...
0127                  | e2e * (fface2redge>0)' * fmdl.elem2face' ...
0128                  | e2e * fedge2redge' * fmdl.edge2elem;
0129    % exclude inclusion (dealt with separately)
0130    rtet2ftet = rtet2ftet & ~rtet_in_ftet & ~ftet_in_rtet'; 
0131    progress_msg(sprintf('Found %d',nnz(rtet2ftet)), Inf);
0132 
0133    
0134    progress_msg('Calculate intersection volumes...');
0135    % sparse logical multiplication doesn't exist
0136    rtet2intpt1 = logical(rmdl.edge2elem'*redge2intpt1)';
0137    ftet2intpt1 = logical(fmdl.elem2face *fface2intpt1)';
0138    
0139    rtet2intpt2 = logical(rmdl.elem2face * rface2intpt2)';
0140    ftet2intpt2 = logical(fmdl.edge2elem'* fedge2intpt2)';
0141    
0142    ftet2intpt3 = logical(fmdl.edge2elem'* fedge2intpt3)';
0143    rtet2intpt3 = logical(rmdl.edge2elem'* redge2intpt3)';
0144     
0145    rtet_todo = find(sum(rtet2ftet,2)>0);
0146    C = []; F = []; V = [];
0147    
0148    id = 0; N = length(rtet_todo);
0149    mint = ceil(N/100);
0150    
0151    problem = false;
0152    
0153    for v = rtet_todo'
0154       id = id+1;
0155       if mod(id,mint)==0, progress_msg(id/N); end
0156       tet_todo = find(rtet2ftet(v,:));
0157       common_intpts1 = bsxfun(@and,rtet2intpt1(:,v), ftet2intpt1(:,tet_todo));
0158       common_intpts2 = bsxfun(@and,rtet2intpt2(:,v), ftet2intpt2(:,tet_todo));
0159       common_intpts3 = bsxfun(@and,rtet2intpt3(:,v), ftet2intpt3(:,tet_todo));
0160       f_nodes     = bsxfun(@and,fnode2rtet(:,v), fmdl.node2elem(:,tet_todo));
0161       r_nodes     = bsxfun(@and,rnode2ftet(:,tet_todo), rmdl.node2elem(:,v));
0162       C = [C; v*ones(numel(tet_todo),1)];
0163       F = [F; tet_todo'];
0164       last_v = numel(V);
0165       V = [V; zeros(numel(tet_todo),1)]; % pre-allocate
0166       
0167       for t = 1:numel(tet_todo)
0168          pts = [ intpts1(common_intpts1(:,t),:);
0169             intpts2(common_intpts2(:,t),:);
0170             intpts3(common_intpts3(:,t),:);
0171             fmdl.nodes(f_nodes(:,t),:);
0172             rmdl.nodes(r_nodes(:,t),:)];
0173          last_v = last_v + 1;
0174          [~,V(last_v)] = convhulln_clean(pts,fmdl);
0175       end
0176    end
0177    progress_msg(Inf);
0178     
0179    
0180     c2f = sparse(F,C,V,size(fmdl.elems,1),size(rmdl.elems,1));
0181     
0182     % add rtet contained in ftet
0183     try rmdl = rmfield(rmdl,'coarse2fine'); end % messes with volume
0184     c2f = c2f + bsxfun(@times, sparse(rtet_in_ftet), get_elem_volume(rmdl))';
0185     
0186     % normalize to tet volume
0187     vol = get_elem_volume(fmdl);
0188     c2f = bsxfun(@rdivide,c2f,vol);
0189 
0190     % count identical tets just once
0191     ftri_in_rtri(rtet_in_ftet') = 0;
0192    
0193     
0194     % add tets contained in vox
0195     c2f = c2f + ftet_in_rtet;
0196     
0197     if problem
0198        warning('eidors:mk_tet_c2f:convhulln_issues', ...
0199           sprintf(['There were some problems with convhulln when running mk_tet_c2f. \n' ...
0200                    'Most often these are caused by numerical precision issues and can safely be ignored. \n' ...
0201                    'To save a plot of each problematic intersection, execute these commands:\n' ...
0202                    '  eidors_cache off\n' ...
0203                    '  eidors_debug on mk_tet_c2f\n' ...
0204                    'before re-running your code. Images will be saved to current directory\n' ...
0205                    'Alternatively, use the CHECK_C2F_QUALITY function.' ]));
0206     end                
0207           
0208        
0209        
0210     
0211     
0212     
0213 %-------------------------------------------------------------------------%
0214 % Calculate intersection points between faces and edges
0215 function [intpts, tri2edge, tri2intpt, edge2intpt] = edge2face_intersections(fmdl,rmdl,opt)
0216    N_edges = size(rmdl.edges,1);
0217    N_faces = size(fmdl.faces,1);
0218    
0219    face_bb = zeros(N_faces,6);
0220    face_bb(:,1) = min(reshape(fmdl.nodes(fmdl.faces,1),N_faces,3),[],2);
0221    face_bb(:,2) = max(reshape(fmdl.nodes(fmdl.faces,1),N_faces,3),[],2);
0222    face_bb(:,3) = min(reshape(fmdl.nodes(fmdl.faces,2),N_faces,3),[],2);
0223    face_bb(:,4) = max(reshape(fmdl.nodes(fmdl.faces,2),N_faces,3),[],2);
0224    face_bb(:,5) = min(reshape(fmdl.nodes(fmdl.faces,3),N_faces,3),[],2);
0225    face_bb(:,6) = max(reshape(fmdl.nodes(fmdl.faces,3),N_faces,3),[],2);
0226    
0227    edge_bb = zeros(N_edges,6);
0228    edge_bb(:,1) = min(reshape(rmdl.nodes(rmdl.edges,1),N_edges,2),[],2);
0229    edge_bb(:,2) = max(reshape(rmdl.nodes(rmdl.edges,1),N_edges,2),[],2);
0230    edge_bb(:,3) = min(reshape(rmdl.nodes(rmdl.edges,2),N_edges,2),[],2);
0231    edge_bb(:,4) = max(reshape(rmdl.nodes(rmdl.edges,2),N_edges,2),[],2);
0232    edge_bb(:,5) = min(reshape(rmdl.nodes(rmdl.edges,3),N_edges,2),[],2);
0233    edge_bb(:,6) = max(reshape(rmdl.nodes(rmdl.edges,3),N_edges,2),[],2);
0234    
0235    allocsz = max(N_edges,N_faces);
0236    N_alloc = allocsz;
0237    
0238    intpts = zeros(N_edges,3);
0239    T = zeros(N_edges,1);
0240    E = zeros(N_edges,1);
0241    I = zeros(N_edges,1);
0242   
0243    P1 = rmdl.nodes(rmdl.edges(:,1),:);
0244    P12 = P1 - rmdl.nodes(rmdl.edges(:,2),:);
0245 
0246    
0247    d = sum(fmdl.normals .* fmdl.nodes(fmdl.faces(:,1),:),2);
0248       
0249    mint = ceil(N_edges/100);
0250    
0251    % for point_in_triangle
0252    v0 = fmdl.nodes(fmdl.faces(:,3),:) - fmdl.nodes(fmdl.faces(:,1),:);
0253    v1 = fmdl.nodes(fmdl.faces(:,2),:) - fmdl.nodes(fmdl.faces(:,1),:);
0254    dot00 = dot(v0, v0, 2);
0255    dot01 = dot(v0, v1, 2);
0256    % dot02 = dot(v0, v2, 2);
0257    dot11 = dot(v1, v1, 2);
0258    % dot12 = dot(v1, v2, 2);
0259    invDenom = 1 ./ (dot00 .* dot11 - dot01 .* dot01);
0260    
0261    epsilon = opt.tol_edge2tri;
0262    
0263    excl =   bsxfun(@gt, face_bb(:,1), edge_bb(:,2)') ...
0264           | bsxfun(@lt, face_bb(:,2), edge_bb(:,1)') ...
0265           | bsxfun(@gt, face_bb(:,3), edge_bb(:,4)') ...
0266           | bsxfun(@lt, face_bb(:,4), edge_bb(:,3)') ...
0267           | bsxfun(@gt, face_bb(:,5), edge_bb(:,6)') ...
0268           | bsxfun(@lt, face_bb(:,6), edge_bb(:,5)');
0269    excl = ~excl;
0270    N_pts = 0;
0271    for i = 1:N_edges
0272       if mod(i,mint)==0, progress_msg(i/N_edges); end
0273      
0274       fidx = excl(:,i);
0275       if ~any(fidx), continue, end;
0276       
0277       num = -d(fidx) + sum(bsxfun(@times,fmdl.normals(fidx,:),P1(i,:)),2);
0278       
0279       den = sum(bsxfun(@times,fmdl.normals(fidx,:),P12(i,:)),2);
0280       
0281       u = num ./ den;
0282       % den == 0 => normal perpendicular to line
0283       idx = u >= 0 & u <= 1 & abs(den) > eps;
0284       
0285       % calculate the intersection points
0286       if any(idx)
0287          id = find(idx);
0288          ipts = bsxfun(@minus, P1(i,:), bsxfun(@times, u(id), P12(i,:)));
0289          
0290          if 1
0291             fcs = find(fidx);
0292             fid = fcs(id);
0293             % point in triangle test
0294             v2 = bsxfun(@minus,ipts,fmdl.nodes(fmdl.faces(fid,1),:));
0295             dot02 = dot(v0(fid,:),v2,2);
0296             dot12 = dot(v1(fid,:),v2,2);
0297             % barycentric coordinates
0298             u = (dot11(fid) .* dot02 - dot01(fid) .* dot12) .* invDenom(fid);
0299             v = (dot00(fid) .* dot12 - dot01(fid) .* dot02) .* invDenom(fid);
0300             t = u >= -epsilon & v >= -epsilon & (u+v-epsilon) <= 1; 
0301          else
0302             t = point_in_triangle(ipts,fmdl.faces(id,:),fmdl.nodes,epsilon,'match');
0303          end
0304          if any(t)
0305             N = nnz(t);
0306             if N_pts+N > N_alloc
0307                N_alloc = N_alloc + allocsz;
0308                intpts(N_alloc,3) = 0;
0309                I(N_alloc) = 0;
0310                T(N_alloc) = 0;
0311                E(N_alloc) = 0;
0312             end
0313             idv = N_pts + (1:N);
0314             intpts(idv,:) = ipts(t,:);
0315             I(idv) = idv;
0316             T(idv) = fid(t);
0317             E(idv) = i;
0318             N_pts = N_pts + N;
0319          end
0320       end
0321    end
0322    T = T(1:N_pts);
0323    E = E(1:N_pts);
0324    I = I(1:N_pts);
0325    intpts = intpts(1:N_pts,:);
0326    tri2edge = sparse(T,E,I,size(fmdl.faces,1),size(rmdl.edges,1));
0327    tri2intpt = sparse(T,I,ones(size(I)),size(fmdl.faces,1),size(I,1));
0328    edge2intpt  = sparse(E,I,ones(size(I)),size(rmdl.edges,1),size(I,1));    
0329    
0330 %-------------------------------------------------------------------------%
0331 % Assign each rmdl node to the tet it is in (nodes on tet faces are counted
0332 % mutltiple times)
0333 function rnode2tet = get_nodes_in_tets(fmdl,nodes, opt)
0334     
0335    [A,b] = tet_to_inequal(fmdl.nodes,fmdl.elems);
0336    progress_msg(.01);
0337    % This is split to decrease the memory footprint
0338    rnode2tet = (bsxfun(@minus, A(1:4:end,:)*nodes',b(1:4:end)) <= opt.tol_node2tet)';
0339    progress_msg(.21);
0340    for i = 2:4
0341       rnode2tet = rnode2tet & (bsxfun(@minus, A(i:4:end,:)*nodes',b(i:4:end)) <= opt.tol_node2tet)';
0342       progress_msg(.21 + (i-1)*.23);
0343    end
0344 
0345    % exclude coinciding nodes
0346    ex= bsxfun(@eq,nodes(:,1),fmdl.nodes(:,1)') & ...
0347        bsxfun(@eq,nodes(:,2),fmdl.nodes(:,2)') & ...
0348        bsxfun(@eq,nodes(:,3),fmdl.nodes(:,3)');
0349    progress_msg(.94);
0350    rnode2tet(any(ex,2),:) = 0;
0351    rnode2tet = sparse(rnode2tet);
0352    progress_msg(Inf);
0353 
0354 %-------------------------------------------------------------------------%
0355 % Prepare model
0356 function fmdl = prepare_tet_mdl(fmdl)
0357    fmopt.elem2edge = true;
0358    fmopt.edge2elem = true;
0359    fmopt.face2elem = true;
0360    fmopt.node2elem = true;
0361    fmopt.normals   = true;
0362    fmopt.linear_reorder = false; % this is slow and not needed
0363    ll = eidors_msg('log_level',1);
0364    fmdl = fix_model(fmdl,fmopt);
0365    eidors_msg('log_level',ll);
0366    fmdl.node2elem = logical(fmdl.node2elem);
0367    nElem = size(fmdl.elems,1);
0368    nFace = size(fmdl.faces,1);
0369    fmdl.elem2face = sparse(repmat((1:nElem)',1,4),double(fmdl.elem2face),true,nElem,nFace);
0370 
0371 %-------------------------------------------------------------------------%
0372 % Center scale models
0373 function[fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt)
0374    ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0375    rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0376    fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0377    if isfield(opt,'do_not_scale') && opt.do_not_scale
0378       return
0379    end
0380    maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0381    scale = 1/maxnode;
0382    rmdl.nodes = scale*rmdl.nodes;
0383    fmdl.nodes = scale*fmdl.nodes;
0384    eidors_msg('@@ models scaled by %g', scale,2);
0385 
0386 %-------------------------------------------------------------------------%
0387 % Remove obviously non-overlapping parts of the models
0388 function [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl)
0389    f_min = min(fmdl.nodes);
0390    f_max = max(fmdl.nodes);
0391    r_min = min(rmdl.nodes);
0392    r_max = max(rmdl.nodes);
0393    
0394    % nodes outside the bounding box of the other model
0395    f_gt  = bsxfun(@gt, fmdl.nodes, r_max);
0396    f_lt  = bsxfun(@lt, fmdl.nodes, r_min);
0397    r_gt  = bsxfun(@gt, rmdl.nodes, f_max);
0398    r_lt  = bsxfun(@lt, rmdl.nodes, f_min);
0399    
0400    % elems outside the bounding box of the other model
0401    re_gt = any(reshape(all(reshape(r_gt(rmdl.elems',:),4,[])),[],3),2);
0402    re_lt = any(reshape(all(reshape(r_lt(rmdl.elems',:),4,[])),[],3),2);
0403    fe_gt = any(reshape(all(reshape(f_gt(fmdl.elems',:),4,[])),[],3),2);
0404    fe_lt = any(reshape(all(reshape(f_lt(fmdl.elems',:),4,[])),[],3),2);
0405    
0406    % elems to keep
0407    rmdl_idx = ~(re_gt | re_lt);
0408    fmdl_idx = ~(fe_gt | fe_lt);
0409    
0410    % remove non-overlapping elems
0411    rmdl.elems = rmdl.elems(rmdl_idx,:);
0412    fmdl.elems = fmdl.elems(fmdl_idx,:);
0413    
0414    % remove unused nodes
0415    [r_used_nodes,jnk,r_n] = unique(rmdl.elems(:));
0416    [f_used_nodes,jnk,f_n] = unique(fmdl.elems(:));
0417    
0418    r_idx = 1:numel(r_used_nodes);
0419    f_idx = 1:numel(f_used_nodes);
0420    
0421    rmdl.elems = reshape(r_idx(r_n),[],4);
0422    fmdl.elems = reshape(f_idx(f_n),[],4);
0423    
0424    rmdl.nodes = rmdl.nodes(r_used_nodes,:);
0425    fmdl.nodes = fmdl.nodes(f_used_nodes,:);
0426    
0427    % for the benefit of any (debug) plots later on
0428    try, rmdl = rmfield(rmdl,'boundary'); end
0429    try, fmdl = rmfield(fmdl,'boundary'); end
0430     
0431 %-------------------------------------------------------------------------%
0432 % Parse option struct
0433  function opt = parse_opts(fmdl,rmdl, opt)
0434 
0435     
0436     if ~isfield(opt, 'tol_node2tet');
0437         opt.tol_node2tet = eps; % * max(rmdl_rng,fmdl_rng)^3;
0438     end
0439     if ~isfield(opt, 'tol_edge2edge')
0440         opt.tol_edge2edge = 6*eps(min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))));
0441     end
0442     if ~isfield(opt, 'tol_edge2tri')
0443         opt.tol_edge2tri = eps; %1e-10
0444     end
0445 %     if ~isfield(opt, 'save_memory')
0446 %        opt.save_memory = 0;
0447 %     end
0448     eidors_msg('@@ node2tet  tolerance = %g', opt.tol_node2tet,2);
0449     eidors_msg('@@ edge2edge tolerance = %g', opt.tol_edge2edge,2);
0450     eidors_msg('@@ edge2tri  tolerance = %g', opt.tol_edge2tri,2);
0451    
0452    
0453 %-------------------------------------------------------------------------%
0454 % Perfom unit tests
0455 function do_unit_test
0456    do_case_test;
0457    do_small_test;
0458    do_realistic_test;
0459 
0460    
0461 function do_case_test
0462    ll = eidors_msg('log_level');
0463    eidors_msg('log_level',1);
0464    t1.type = 'fwd_model';
0465    t1.elems = [1 2 3 4];
0466 
0467    X = 2; Y = 3; % subplot matrix
0468    for i = 1:30
0469         t1.nodes = [0 0 0; 0 1 0; 1 0 0; 0 0 1];
0470         t2 = t1;
0471         switch i
0472            case 1
0473               txt = 'identical';
0474               subplot(X,Y,i), show_test(t1,t2);
0475               c2f = mk_tet_c2f(t1,t2);
0476               tol = eps; cmptarg = 1;
0477            case 2
0478               txt = 'shared face';
0479               t2.nodes(end,end) = -1;
0480               subplot(X,Y,i), show_test(t1,t2);
0481               c2f = mk_tet_c2f(t1,t2);
0482               tol = 0; cmptarg = 0;
0483            case 3
0484               txt = 'coplanar faces';
0485               t2.nodes(end,end) = -1;
0486               t1.nodes(:,1:2) = t1.nodes(:,1:2) + .2;
0487               subplot(X,Y,i), show_test(t1,t2);
0488               c2f = mk_tet_c2f(t1,t2);
0489               tol = 0; cmptarg = 0;
0490            case 4
0491               txt = 'point on edge';
0492               t2.nodes(:,1) = t1.nodes(:,1) + 1;
0493               t2.nodes(:,2) = t1.nodes(:,2) - .3;
0494               subplot(X,Y,i), show_test(t1,t2);
0495               c2f = mk_tet_c2f(t1,t2);
0496               tol = 0; cmptarg = 0;
0497            otherwise
0498              break;
0499         end
0500         txt = sprintf('%02d: %s',i,txt);
0501         unit_test_cmp(txt,c2f,cmptarg,tol);
0502    end
0503    eidors_msg('log_level',ll);
0504 
0505 function do_small_test
0506    fmdl = ng_mk_cyl_models([1 .5],[],[]);
0507    show_fem(fmdl)
0508    v = -.5:.1:.5;
0509    rmdl = mk_grid_model([],v,v,0:.1:1);
0510    hold on
0511    h = show_fem(rmdl);
0512    set(h,'edgecolor','b');
0513    hold off
0514    c2f = mk_tet_c2f(fmdl,rmdl);
0515    tc2f = c2f * rmdl.coarse2fine;
0516    vc2f = mk_grid_c2f(fmdl,rmdl);
0517    unit_test_cmp('mk_tet_c2f v mk_grid_c2f', tc2f,vc2f, 1e-15);
0518 
0519 
0520 function do_realistic_test
0521    fmdl= ng_mk_cyl_models([2,2,.1],[16,1],[.1,0,.025]);
0522    xvec = [-1.5 -.5:.2:.5 1.5];
0523    yvec = [-1.6 -1:.2:1 1.6];
0524    zvec = 0:.25:2;
0525    rmdl = mk_grid_model([],xvec,yvec,zvec);
0526    tic
0527    opt.save_memory = 0;
0528    c2f_a = mk_grid_c2f(fmdl, rmdl,opt);
0529    t = toc;
0530    fprintf('Voxel: t=%f s\n',t);
0531 
0532    tic
0533    opt.save_memory = 0;
0534    c2f_b = mk_tet_c2f(fmdl, rmdl,opt);
0535    t = toc;
0536    fprintf('Tet: t=%f s\n',t);
0537 
0538    c2f_b = c2f_b * rmdl.coarse2fine;
0539    unit_test_cmp('mk_tet_c2f v mk_grid_c2f', c2f_b,c2f_a, 1e-5);
0540 
0541    tic
0542    c2f_n = mk_approx_c2f(fmdl,rmdl);
0543    t = toc;
0544    fprintf('Approximate: t=%f s\n',t);
0545 
0546  function show_test(vox,tet, pts)
0547     show_fem(vox);
0548     hold on
0549     h = show_fem(tet);
0550     set(h, 'EdgeColor','b');
0551     if nargin > 2
0552        plot3(pts(:,1),pts(:,2),pts(:,3),'o');
0553     end
0554     hold off
0555     axis auto
mk_tet_c2f_bak

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
