mk_tet_c2f

PURPOSE ^

MK_TET_C2F - calculate a coarse2fine mapping for two tet-based models.

SYNOPSIS ^

function [c2f] = mk_tet_c2f(fmdl, rmdl, opt)

DESCRIPTION ^

MK_TET_C2F - calculate a coarse2fine mapping for two tet-based models.
 C2F = MK_TET_C2F(FMDL,RMDL) returns in C2F the fraction of volume of
 each element of the fine model FMDL contained in each element of
 the coarse model RMDL.
 Uses CONVHULLN to calculate the volume defined by a set of intersection
 points between individual tet elements.

 C2F = MK_TET_C2F(FMDL,RMDL,OPT) allows specifying options.
 
 Inputs:
   FMDL - a (fine) EIDORS (tet-based) forward model
   RMDL - a (course) EIDORS (tet-based) forward model
   OPT  - an option structure with the following fields and defaults:
      .do_not_scale  - set to true to prevent scaling the models to unit
                       cube before any calculations, including thresholds.
                       Default: false
      .tol_node2tet  - tolerance for determinant <= 0 in testing for
                       points inside tets. Default: eps
      .tol_edge2edge - maximum distance between "intersecting" edges
                       Default: 6*sqrt(3)*eps(a), where a is
                       min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))
      .tol_edge2tri  - minimum value of a barycentric coordinate to 
                       decide a point is lying inside a triangle and not
                       on its edge. Default: eps

 NOTE that for grid-based models, such as returned by MK_GRID_MODEL or
 MK_VOXEL_VOLUME, MK_GRID_C2F is much faster.

 Set eidors_msg 'log level' < 2 to supress output to command line.

 Examples:
     fmdl = ng_mk_cyl_models([2,2,.2],[],[]);
     rmdl = ng_mk_cyl_models([2,2],[],[]);
     c2f  = mk_tet_c2f(fmdl,rmdl);
     h = show_fem(fmdl); set(h,'LineWidth',0.1)
     hold on
     h = show_fem(rmdl); set(h,'EdgeColor','b','LineWidth',2);
     hold off

 See also MK_GRID_C2F, FIND_EDGE2EDGE_INTERSECTIONS, CONVHULLN
     MK_COARSE_FINE_MAPPING, MK_APPROX_C2F, POINT_IN_TRIANGLE, EIDORS_MSG

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 5336 2017-02-16 18:38:35Z bgrychtol $
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          if size(pts,1) < 4 
0175             % there are some degenerate cases, sometimes caused by
0176             % numerical issues alone
0177             continue
0178          end
0179          try
0180             % move points to origin (helps for small elements at
0181             % large coordinates
0182             ctr = mean(pts);
0183             pts = bsxfun(@minus,pts,ctr);
0184             scale = max(abs(pts(:)));
0185             if scale == 0 %happens when there's only one point
0186                continue
0187             end
0188             % scale largest coordinate to 1 (helps with precision)
0189             pts = pts ./ scale;
0190             % force thorough search for initinal simplex and
0191             % supress precision warnings
0192             [K, V(last_v)] = convhulln(pts,{'Qt Pp Qs'});
0193             V(last_v) = V(last_v) * scale^3; % undo scaling
0194          catch err
0195             ok = false;
0196             switch err.identifier
0197                case {'MATLAB:qhullmx:DegenerateData', 'MATLAB:qhullmx:UndefinedError'}
0198                   if size(pts,1) > 3
0199                      u = uniquetol(pts*scale,6*eps,'ByRows',true,'DataScale', 1);
0200                      ok = ok | size(u,1) < 4;
0201                   end
0202             end
0203             if ~ok
0204                if DEBUG || eidors_debug('query','mk_tet_c2f:convhulln')
0205                   tet.nodes = fmdl.nodes;
0206                   vox.nodes = rmdl.nodes;
0207                   tet.type = 'fwd_model';
0208                   vox.type = 'fwd_model';
0209                   vox.elems = rmdl.faces(logical(rmdl.elem2face(v,:)),:);
0210                   vox.boundary = vox.elems;
0211                   tet.elems = fmdl.elems(tet_todo(t),:);
0212                   clf
0213                   pts = bsxfun(@plus,pts*scale,ctr);
0214                   subplot(221)
0215                   show_test(vox,tet,pts);
0216                   
0217                   subplot(222)
0218                   show_test(vox,tet,pts);
0219                   view(90,0)
0220                   
0221                   subplot(223)
0222                   show_test(vox,tet,pts);
0223                   view(0,90)
0224                   
0225                   subplot(224)
0226                   show_test(vox,tet,pts);
0227                   view(0,0)
0228                   
0229                   
0230                   str = sprintf('mk_tet_c2f problem fe %d ce %d', v, tet_todo(t));
0231                   print(gcf,'-dpng',str);
0232 %                   keyboard
0233                else
0234                   problem = true;
0235 %                   fprintf('\n');
0236 %                   eidors_msg(['convhulln has thrown an error. ' ...
0237 %                      'Enable eidors_debug on mk_tet_c2f and re-run to see a debug plot'],0);
0238 %                   rethrow(err);
0239                end
0240             end
0241          end
0242       end
0243    end
0244    progress_msg(Inf);
0245     
0246    
0247     c2f = sparse(F,C,V,size(fmdl.elems,1),size(rmdl.elems,1));
0248     
0249     % add rtet contained in ftet
0250     try rmdl = rmfield(rmdl,'coarse2fine'); end % messes with volume
0251     c2f = c2f + bsxfun(@times, sparse(rtet_in_ftet), get_elem_volume(rmdl))';
0252     
0253     % normalize to tet volume
0254     vol = get_elem_volume(fmdl);
0255     c2f = bsxfun(@rdivide,c2f,vol);
0256 
0257     % count identical tets just once
0258     ftri_in_rtri(rtet_in_ftet') = 0;
0259    
0260     
0261     % add tets contained in vox
0262     c2f = c2f + ftet_in_rtet;
0263     
0264     if problem
0265        warning('eidors:mk_tet_c2f:convhulln_issues', ...
0266           sprintf(['There were some problems with convhulln when running mk_tet_c2f. \n' ...
0267                    'Most often these are caused by numerical precision issues and can safely be ignored. \n' ...
0268                    'To save a plot of each problematic intersection, execute these commands:\n' ...
0269                    '  eidors_cache off\n' ...
0270                    '  eidors_debug on mk_tet_c2f\n' ...
0271                    'before re-running your code. Images will be saved to current directory\n' ...
0272                    'Alternatively, use the CHECK_C2F_QUALITY function.' ]));
0273     end                
0274           
0275        
0276        
0277     
0278     
0279     
0280 %-------------------------------------------------------------------------%
0281 % Calculate intersection points between faces and edges
0282 function [intpts, tri2edge, tri2intpt, edge2intpt] = edge2face_intersections(fmdl,rmdl,opt)
0283    N_edges = size(rmdl.edges,1);
0284    N_faces = size(fmdl.faces,1);
0285    
0286    face_bb = zeros(N_faces,6);
0287    face_bb(:,1) = min(reshape(fmdl.nodes(fmdl.faces,1),N_faces,3),[],2);
0288    face_bb(:,2) = max(reshape(fmdl.nodes(fmdl.faces,1),N_faces,3),[],2);
0289    face_bb(:,3) = min(reshape(fmdl.nodes(fmdl.faces,2),N_faces,3),[],2);
0290    face_bb(:,4) = max(reshape(fmdl.nodes(fmdl.faces,2),N_faces,3),[],2);
0291    face_bb(:,5) = min(reshape(fmdl.nodes(fmdl.faces,3),N_faces,3),[],2);
0292    face_bb(:,6) = max(reshape(fmdl.nodes(fmdl.faces,3),N_faces,3),[],2);
0293    
0294    edge_bb = zeros(N_edges,6);
0295    edge_bb(:,1) = min(reshape(rmdl.nodes(rmdl.edges,1),N_edges,2),[],2);
0296    edge_bb(:,2) = max(reshape(rmdl.nodes(rmdl.edges,1),N_edges,2),[],2);
0297    edge_bb(:,3) = min(reshape(rmdl.nodes(rmdl.edges,2),N_edges,2),[],2);
0298    edge_bb(:,4) = max(reshape(rmdl.nodes(rmdl.edges,2),N_edges,2),[],2);
0299    edge_bb(:,5) = min(reshape(rmdl.nodes(rmdl.edges,3),N_edges,2),[],2);
0300    edge_bb(:,6) = max(reshape(rmdl.nodes(rmdl.edges,3),N_edges,2),[],2);
0301    
0302    allocsz = max(N_edges,N_faces);
0303    N_alloc = allocsz;
0304    
0305    intpts = zeros(N_edges,3);
0306    T = zeros(N_edges,1);
0307    E = zeros(N_edges,1);
0308    I = zeros(N_edges,1);
0309   
0310    P1 = rmdl.nodes(rmdl.edges(:,1),:);
0311    P12 = P1 - rmdl.nodes(rmdl.edges(:,2),:);
0312 
0313    
0314    d = sum(fmdl.normals .* fmdl.nodes(fmdl.faces(:,1),:),2);
0315       
0316    mint = ceil(N_edges/100);
0317    
0318    % for point_in_triangle
0319    v0 = fmdl.nodes(fmdl.faces(:,3),:) - fmdl.nodes(fmdl.faces(:,1),:);
0320    v1 = fmdl.nodes(fmdl.faces(:,2),:) - fmdl.nodes(fmdl.faces(:,1),:);
0321    dot00 = dot(v0, v0, 2);
0322    dot01 = dot(v0, v1, 2);
0323    % dot02 = dot(v0, v2, 2);
0324    dot11 = dot(v1, v1, 2);
0325    % dot12 = dot(v1, v2, 2);
0326    invDenom = 1 ./ (dot00 .* dot11 - dot01 .* dot01);
0327    
0328    epsilon = opt.tol_edge2tri;
0329    
0330    excl =   bsxfun(@gt, face_bb(:,1), edge_bb(:,2)') ...
0331           | bsxfun(@lt, face_bb(:,2), edge_bb(:,1)') ...
0332           | bsxfun(@gt, face_bb(:,3), edge_bb(:,4)') ...
0333           | bsxfun(@lt, face_bb(:,4), edge_bb(:,3)') ...
0334           | bsxfun(@gt, face_bb(:,5), edge_bb(:,6)') ...
0335           | bsxfun(@lt, face_bb(:,6), edge_bb(:,5)');
0336    excl = ~excl;
0337    N_pts = 0;
0338    for i = 1:N_edges
0339       if mod(i,mint)==0, progress_msg(i/N_edges); end
0340      
0341       fidx = excl(:,i);
0342       if ~any(fidx), continue, end;
0343       
0344       num = -d(fidx) + sum(bsxfun(@times,fmdl.normals(fidx,:),P1(i,:)),2);
0345       
0346       den = sum(bsxfun(@times,fmdl.normals(fidx,:),P12(i,:)),2);
0347       
0348       u = num ./ den;
0349       % den == 0 => normal perpendicular to line
0350       idx = u >= 0 & u <= 1 & abs(den) > eps;
0351       
0352       % calculate the intersection points
0353       if any(idx)
0354          id = find(idx);
0355          ipts = bsxfun(@minus, P1(i,:), bsxfun(@times, u(id), P12(i,:)));
0356          
0357          if 1
0358             fcs = find(fidx);
0359             fid = fcs(id);
0360             % point in triangle test
0361             v2 = bsxfun(@minus,ipts,fmdl.nodes(fmdl.faces(fid,1),:));
0362             dot02 = dot(v0(fid,:),v2,2);
0363             dot12 = dot(v1(fid,:),v2,2);
0364             % barycentric coordinates
0365             u = (dot11(fid) .* dot02 - dot01(fid) .* dot12) .* invDenom(fid);
0366             v = (dot00(fid) .* dot12 - dot01(fid) .* dot02) .* invDenom(fid);
0367             t = u >= -epsilon & v >= -epsilon & (u+v-epsilon) <= 1; 
0368          else
0369             t = point_in_triangle(ipts,fmdl.faces(id,:),fmdl.nodes,epsilon,'match');
0370          end
0371          if any(t)
0372             N = nnz(t);
0373             if N_pts+N > N_alloc
0374                N_alloc = N_alloc + allocsz;
0375                intpts(N_alloc,3) = 0;
0376                I(N_alloc) = 0;
0377                T(N_alloc) = 0;
0378                E(N_alloc) = 0;
0379             end
0380             idv = N_pts + (1:N);
0381             intpts(idv,:) = ipts(t,:);
0382             I(idv) = idv;
0383             T(idv) = fid(t);
0384             E(idv) = i;
0385             N_pts = N_pts + N;
0386          end
0387       end
0388    end
0389    T = T(1:N_pts);
0390    E = E(1:N_pts);
0391    I = I(1:N_pts);
0392    intpts = intpts(1:N_pts,:);
0393    tri2edge = sparse(T,E,I,size(fmdl.faces,1),size(rmdl.edges,1));
0394    tri2intpt = sparse(T,I,ones(size(I)),size(fmdl.faces,1),size(I,1));
0395    edge2intpt  = sparse(E,I,ones(size(I)),size(rmdl.edges,1),size(I,1));    
0396    
0397 %-------------------------------------------------------------------------%
0398 % Assign each rmdl node to the tet it is in (nodes on tet faces are counted
0399 % mutltiple times)
0400 function rnode2tet = get_nodes_in_tets(fmdl,nodes, opt)
0401     
0402    [A,b] = tet_to_inequal(fmdl.nodes,fmdl.elems);
0403    progress_msg(.01);
0404    % This is split to decrease the memory footprint
0405    rnode2tet = (bsxfun(@minus, A(1:4:end,:)*nodes',b(1:4:end)) <= opt.tol_node2tet)';
0406    progress_msg(.21);
0407    for i = 2:4
0408       rnode2tet = rnode2tet & (bsxfun(@minus, A(i:4:end,:)*nodes',b(i:4:end)) <= opt.tol_node2tet)';
0409       progress_msg(.21 + (i-1)*.23);
0410    end
0411 
0412    % exclude coinciding nodes
0413    ex= bsxfun(@eq,nodes(:,1),fmdl.nodes(:,1)') & ...
0414        bsxfun(@eq,nodes(:,2),fmdl.nodes(:,2)') & ...
0415        bsxfun(@eq,nodes(:,3),fmdl.nodes(:,3)');
0416    progress_msg(.94);
0417    rnode2tet(any(ex,2),:) = 0;
0418    rnode2tet = sparse(rnode2tet);
0419    progress_msg(Inf);
0420 
0421 %-------------------------------------------------------------------------%
0422 % Prepare model
0423 function fmdl = prepare_tet_mdl(fmdl)
0424    fmopt.elem2edge = true;
0425    fmopt.edge2elem = true;
0426    fmopt.face2elem = true;
0427    fmopt.node2elem = true;
0428    fmopt.normals   = true;
0429    fmopt.linear_reorder = false; % this is slow and not needed
0430    ll = eidors_msg('log_level',1);
0431    fmdl = fix_model(fmdl,fmopt);
0432    eidors_msg('log_level',ll);
0433    fmdl.node2elem = logical(fmdl.node2elem);
0434    nElem = size(fmdl.elems,1);
0435    nFace = size(fmdl.faces,1);
0436    fmdl.elem2face = sparse(repmat((1:nElem)',1,4),double(fmdl.elem2face),true,nElem,nFace);
0437 
0438 %-------------------------------------------------------------------------%
0439 % Center scale models
0440 function[fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt)
0441    ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0442    rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0443    fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0444    if isfield(opt,'do_not_scale') && opt.do_not_scale
0445       return
0446    end
0447    maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0448    scale = 1/maxnode;
0449    rmdl.nodes = scale*rmdl.nodes;
0450    fmdl.nodes = scale*fmdl.nodes;
0451    eidors_msg('@@ models scaled by %g', scale,2);
0452 
0453 %-------------------------------------------------------------------------%
0454 % Remove obviously non-overlapping parts of the models
0455 function [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl)
0456    f_min = min(fmdl.nodes);
0457    f_max = max(fmdl.nodes);
0458    r_min = min(rmdl.nodes);
0459    r_max = max(rmdl.nodes);
0460    
0461    % nodes outside the bounding box of the other model
0462    f_gt  = bsxfun(@gt, fmdl.nodes, r_max);
0463    f_lt  = bsxfun(@lt, fmdl.nodes, r_min);
0464    r_gt  = bsxfun(@gt, rmdl.nodes, f_max);
0465    r_lt  = bsxfun(@lt, rmdl.nodes, f_min);
0466    
0467    % elems outside the bounding box of the other model
0468    re_gt = any(reshape(all(reshape(r_gt(rmdl.elems',:),4,[])),[],3),2);
0469    re_lt = any(reshape(all(reshape(r_lt(rmdl.elems',:),4,[])),[],3),2);
0470    fe_gt = any(reshape(all(reshape(f_gt(fmdl.elems',:),4,[])),[],3),2);
0471    fe_lt = any(reshape(all(reshape(f_lt(fmdl.elems',:),4,[])),[],3),2);
0472    
0473    % elems to keep
0474    rmdl_idx = ~(re_gt | re_lt);
0475    fmdl_idx = ~(fe_gt | fe_lt);
0476    
0477    % remove non-overlapping elems
0478    rmdl.elems = rmdl.elems(rmdl_idx,:);
0479    fmdl.elems = fmdl.elems(fmdl_idx,:);
0480    
0481    % remove unused nodes
0482    [r_used_nodes,jnk,r_n] = unique(rmdl.elems(:));
0483    [f_used_nodes,jnk,f_n] = unique(fmdl.elems(:));
0484    
0485    r_idx = 1:numel(r_used_nodes);
0486    f_idx = 1:numel(f_used_nodes);
0487    
0488    rmdl.elems = reshape(r_idx(r_n),[],4);
0489    fmdl.elems = reshape(f_idx(f_n),[],4);
0490    
0491    rmdl.nodes = rmdl.nodes(r_used_nodes,:);
0492    fmdl.nodes = fmdl.nodes(f_used_nodes,:);
0493    
0494    % for the benefit of any (debug) plots later on
0495    try, rmdl = rmfield(rmdl,'boundary'); end
0496    try, fmdl = rmfield(fmdl,'boundary'); end
0497     
0498 %-------------------------------------------------------------------------%
0499 % Parse option struct
0500  function opt = parse_opts(fmdl,rmdl, opt)
0501 
0502     
0503     if ~isfield(opt, 'tol_node2tet');
0504         opt.tol_node2tet = eps; % * max(rmdl_rng,fmdl_rng)^3;
0505     end
0506     if ~isfield(opt, 'tol_edge2edge')
0507         opt.tol_edge2edge = 6*eps(min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))));
0508     end
0509     if ~isfield(opt, 'tol_edge2tri')
0510         opt.tol_edge2tri = eps; %1e-10
0511     end
0512 %     if ~isfield(opt, 'save_memory')
0513 %        opt.save_memory = 0;
0514 %     end
0515     eidors_msg('@@ node2tet  tolerance = %g', opt.tol_node2tet,2);
0516     eidors_msg('@@ edge2edge tolerance = %g', opt.tol_edge2edge,2);
0517     eidors_msg('@@ edge2tri  tolerance = %g', opt.tol_edge2tri,2);
0518    
0519    
0520 %-------------------------------------------------------------------------%
0521 % Perfom unit tests
0522 function do_unit_test
0523    do_case_test;
0524    do_small_test;
0525    do_realistic_test;
0526 
0527    
0528 function do_case_test
0529    ll = eidors_msg('log_level');
0530    eidors_msg('log_level',1);
0531    t1.type = 'fwd_model';
0532    t1.elems = [1 2 3 4];
0533 
0534    X = 2; Y = 3; % subplot matrix
0535    for i = 1:30
0536         fprintf('%d\n',i);
0537         t1.nodes = [0 0 0; 0 1 0; 1 0 0; 0 0 1];
0538         t2 = t1;
0539         switch i
0540            case 1
0541               txt = 'identical';
0542               subplot(X,Y,i), show_test(t1,t2);
0543               c2f = mk_tet_c2f(t1,t2);
0544               unit_test_cmp(txt,c2f,1,eps);
0545            case 2
0546               txt = 'shared face';
0547               t2.nodes(end,end) = -1;
0548               subplot(X,Y,i), show_test(t1,t2);
0549               c2f = mk_tet_c2f(t1,t2);
0550               unit_test_cmp(txt,c2f,0,0);
0551            case 3
0552               txt = 'coplanar faces';
0553               t2.nodes(end,end) = -1;
0554               t1.nodes(:,1:2) = t1.nodes(:,1:2) + .2;
0555               subplot(X,Y,i), show_test(t1,t2);
0556               c2f = mk_tet_c2f(t1,t2);
0557               unit_test_cmp(txt,c2f,0,0);
0558            case 4
0559               txt = 'point on edge';
0560               t2.nodes(:,1) = t1.nodes(:,1) + 1;
0561               t2.nodes(:,2) = t1.nodes(:,2) - .3;
0562               subplot(X,Y,i), show_test(t1,t2);
0563               c2f = mk_tet_c2f(t1,t2);
0564               unit_test_cmp(txt,c2f,0,0);
0565            otherwise
0566              break;
0567         end
0568 
0569       
0570    end
0571    eidors_msg('log_level',ll);
0572 
0573 function do_small_test
0574    fmdl = ng_mk_cyl_models([1 .5],[],[]);
0575    show_fem(fmdl)
0576    v = -.5:.1:.5;
0577    rmdl = mk_grid_model([],v,v,0:.1:1);
0578    hold on
0579    h = show_fem(rmdl);
0580    set(h,'edgecolor','b');
0581    hold off
0582    c2f = mk_tet_c2f(fmdl,rmdl);
0583    tc2f = c2f * rmdl.coarse2fine;
0584    vc2f = mk_grid_c2f(fmdl,rmdl);
0585    unit_test_cmp('mk_tet_c2f v mk_grid_c2f', tc2f,vc2f, 1e-15);
0586 
0587 
0588 function do_realistic_test
0589    fmdl= ng_mk_cyl_models([2,2,.1],[16,1],[.1,0,.025]);
0590    xvec = [-1.5 -.5:.2:.5 1.5];
0591    yvec = [-1.6 -1:.2:1 1.6];
0592    zvec = 0:.25:2;
0593    rmdl = mk_grid_model([],xvec,yvec,zvec);
0594    tic
0595    opt.save_memory = 0;
0596    c2f_a = mk_grid_c2f(fmdl, rmdl,opt);
0597    t = toc;
0598    fprintf('Voxel: t=%f s\n',t);
0599 
0600    tic
0601    opt.save_memory = 0;
0602    c2f_b = mk_tet_c2f(fmdl, rmdl,opt);
0603    t = toc;
0604    fprintf('Tet: t=%f s\n',t);
0605 
0606    c2f_b = c2f_b * rmdl.coarse2fine;
0607    unit_test_cmp('mk_tet_c2f v mk_grid_c2f', c2f_b,c2f_a, 1e-5);
0608 
0609    tic
0610    c2f_n = mk_approx_c2f(fmdl,rmdl);
0611    t = toc;
0612    fprintf('Approximate: t=%f s\n',t);
0613 
0614  function show_test(vox,tet, pts)
0615     show_fem(vox);
0616     hold on
0617     h = show_fem(tet);
0618     set(h, 'EdgeColor','b');
0619     if nargin > 2
0620        plot3(pts(:,1),pts(:,2),pts(:,3),'o');
0621     end
0622     hold off
0623     axis auto

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005