mk_tet_c2f

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

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 5996 2019-06-27 17:39:56Z 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          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    chunk_size = 100;
0331    chunk_end = 0;
0332    N_pts = 0;
0333    for i = 1:N_edges
0334       if mod(i,mint)==0, progress_msg(i/N_edges); end
0335       
0336       if i > chunk_end
0337         chunk_start = chunk_end + 1;
0338         chunk_end = min(chunk_end+chunk_size,N_edges);
0339         chunk = chunk_start:chunk_end;
0340         excl =   face_bb(:,1) > edge_bb(chunk,2)' ...
0341             | face_bb(:,2) < edge_bb(chunk,1)' ...
0342             | face_bb(:,3) > edge_bb(chunk,4)' ...
0343             | face_bb(:,4) < edge_bb(chunk,3)' ...
0344             | face_bb(:,5) > edge_bb(chunk,6)' ...
0345             | face_bb(:,6) < edge_bb(chunk,5)';
0346         excl = ~excl;
0347         chunk_i=1;
0348       end
0349       fidx = excl(:,chunk_i);
0350       chunk_i = chunk_i + 1;
0351       
0352       if ~any(fidx), continue, end;
0353       
0354       num = -d(fidx) + sum(bsxfun(@times,fmdl.normals(fidx,:),P1(i,:)),2);
0355       
0356       den = sum(bsxfun(@times,fmdl.normals(fidx,:),P12(i,:)),2);
0357       
0358       u = num ./ den;
0359       % den == 0 => normal perpendicular to line
0360       idx = u >= 0 & u <= 1 & abs(den) > eps;
0361       
0362       % calculate the intersection points
0363       if any(idx)
0364          id = find(idx);
0365          ipts = bsxfun(@minus, P1(i,:), bsxfun(@times, u(id), P12(i,:)));
0366          
0367          if 1
0368             fcs = find(fidx);
0369             fid = fcs(id);
0370             % point in triangle test
0371             v2 = bsxfun(@minus,ipts,fmdl.nodes(fmdl.faces(fid,1),:));
0372             dot02 = dot(v0(fid,:),v2,2);
0373             dot12 = dot(v1(fid,:),v2,2);
0374             % barycentric coordinates
0375             u = (dot11(fid) .* dot02 - dot01(fid) .* dot12) .* invDenom(fid);
0376             v = (dot00(fid) .* dot12 - dot01(fid) .* dot02) .* invDenom(fid);
0377             t = u >= -epsilon & v >= -epsilon & (u+v-epsilon) <= 1; 
0378          else
0379             t = point_in_triangle(ipts,fmdl.faces(id,:),fmdl.nodes,epsilon,'match');
0380          end
0381          if any(t)
0382             N = nnz(t);
0383             if N_pts+N > N_alloc
0384                N_alloc = N_alloc + allocsz;
0385                intpts(N_alloc,3) = 0;
0386                I(N_alloc) = 0;
0387                T(N_alloc) = 0;
0388                E(N_alloc) = 0;
0389             end
0390             idv = N_pts + (1:N);
0391             intpts(idv,:) = ipts(t,:);
0392             I(idv) = idv;
0393             T(idv) = fid(t);
0394             E(idv) = i;
0395             N_pts = N_pts + N;
0396          end
0397       end
0398       
0399    end
0400    T = T(1:N_pts);
0401    E = E(1:N_pts);
0402    I = I(1:N_pts);
0403    intpts = intpts(1:N_pts,:);
0404    tri2edge = sparse(T,E,I,size(fmdl.faces,1),size(rmdl.edges,1));
0405    tri2intpt = sparse(T,I,ones(size(I)),size(fmdl.faces,1),size(I,1));
0406    edge2intpt  = sparse(E,I,ones(size(I)),size(rmdl.edges,1),size(I,1));    
0407    
0408 %-------------------------------------------------------------------------%
0409 % Assign each rmdl node to the tet it is in (nodes on tet faces are counted
0410 % mutltiple times)
0411 function rnode2tet = get_nodes_in_tets(fmdl,nodes, opt)
0412     
0413    [A,b] = tet_to_inequal(fmdl.nodes,fmdl.elems);
0414    progress_msg(.01);
0415    % This is split to decrease the memory footprint
0416    rnode2tet = (bsxfun(@minus, A(1:4:end,:)*nodes',b(1:4:end)) <= opt.tol_node2tet)';
0417    progress_msg(.21);
0418    for i = 2:4
0419       rnode2tet = rnode2tet & (bsxfun(@minus, A(i:4:end,:)*nodes',b(i:4:end)) <= opt.tol_node2tet)';
0420       progress_msg(.21 + (i-1)*.23);
0421    end
0422 
0423    % exclude coinciding nodes
0424    ex= bsxfun(@eq,nodes(:,1),fmdl.nodes(:,1)') & ...
0425        bsxfun(@eq,nodes(:,2),fmdl.nodes(:,2)') & ...
0426        bsxfun(@eq,nodes(:,3),fmdl.nodes(:,3)');
0427    progress_msg(.94);
0428    rnode2tet(any(ex,2),:) = 0;
0429    rnode2tet = sparse(rnode2tet);
0430    progress_msg(Inf);
0431 
0432 %-------------------------------------------------------------------------%
0433 % Prepare model
0434 function fmdl = prepare_tet_mdl(fmdl)
0435    fmopt.elem2edge = true;
0436    fmopt.edge2elem = true;
0437    fmopt.face2elem = true;
0438    fmopt.node2elem = true;
0439    fmopt.normals   = true;
0440    fmopt.linear_reorder = false; % this is slow and not needed
0441    ll = eidors_msg('log_level',1);
0442    fmdl = fix_model(fmdl,fmopt);
0443    eidors_msg('log_level',ll);
0444    fmdl.node2elem = logical(fmdl.node2elem);
0445    nElem = size(fmdl.elems,1);
0446    nFace = size(fmdl.faces,1);
0447    fmdl.elem2face = sparse(repmat((1:nElem)',1,4),double(fmdl.elem2face),true,nElem,nFace);
0448 
0449 %-------------------------------------------------------------------------%
0450 % Center scale models
0451 function[fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt)
0452    ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0453    rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0454    fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0455    if isfield(opt,'do_not_scale') && opt.do_not_scale
0456       return
0457    end
0458    maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0459    scale = 1/maxnode;
0460    rmdl.nodes = scale*rmdl.nodes;
0461    fmdl.nodes = scale*fmdl.nodes;
0462    eidors_msg('@@ models scaled by %g', scale,2);
0463 
0464 %-------------------------------------------------------------------------%
0465 % Remove obviously non-overlapping parts of the models
0466 function [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl)
0467    f_min = min(fmdl.nodes);
0468    f_max = max(fmdl.nodes);
0469    r_min = min(rmdl.nodes);
0470    r_max = max(rmdl.nodes);
0471    
0472    % nodes outside the bounding box of the other model
0473    f_gt  = bsxfun(@gt, fmdl.nodes, r_max);
0474    f_lt  = bsxfun(@lt, fmdl.nodes, r_min);
0475    r_gt  = bsxfun(@gt, rmdl.nodes, f_max);
0476    r_lt  = bsxfun(@lt, rmdl.nodes, f_min);
0477    
0478    % elems outside the bounding box of the other model
0479    re_gt = any(reshape(all(reshape(r_gt(rmdl.elems',:),4,[])),[],3),2);
0480    re_lt = any(reshape(all(reshape(r_lt(rmdl.elems',:),4,[])),[],3),2);
0481    fe_gt = any(reshape(all(reshape(f_gt(fmdl.elems',:),4,[])),[],3),2);
0482    fe_lt = any(reshape(all(reshape(f_lt(fmdl.elems',:),4,[])),[],3),2);
0483    
0484    % elems to keep
0485    rmdl_idx = ~(re_gt | re_lt);
0486    fmdl_idx = ~(fe_gt | fe_lt);
0487    
0488    % remove non-overlapping elems
0489    rmdl.elems = rmdl.elems(rmdl_idx,:);
0490    fmdl.elems = fmdl.elems(fmdl_idx,:);
0491    
0492    % remove unused nodes
0493    [r_used_nodes,jnk,r_n] = unique(rmdl.elems(:));
0494    [f_used_nodes,jnk,f_n] = unique(fmdl.elems(:));
0495    
0496    r_idx = 1:numel(r_used_nodes);
0497    f_idx = 1:numel(f_used_nodes);
0498    
0499    rmdl.elems = reshape(r_idx(r_n),[],4);
0500    fmdl.elems = reshape(f_idx(f_n),[],4);
0501    
0502    rmdl.nodes = rmdl.nodes(r_used_nodes,:);
0503    fmdl.nodes = fmdl.nodes(f_used_nodes,:);
0504    
0505    % for the benefit of any (debug) plots later on
0506    try, rmdl = rmfield(rmdl,'boundary'); end
0507    try, fmdl = rmfield(fmdl,'boundary'); end
0508     
0509 %-------------------------------------------------------------------------%
0510 % Parse option struct
0511  function opt = parse_opts(fmdl,rmdl, opt)
0512 
0513     
0514     if ~isfield(opt, 'tol_node2tet');
0515         opt.tol_node2tet = eps; % * max(rmdl_rng,fmdl_rng)^3;
0516     end
0517     if ~isfield(opt, 'tol_edge2edge')
0518         opt.tol_edge2edge = 6*eps(min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))));
0519     end
0520     if ~isfield(opt, 'tol_edge2tri')
0521         opt.tol_edge2tri = eps; %1e-10
0522     end
0523 %     if ~isfield(opt, 'save_memory')
0524 %        opt.save_memory = 0;
0525 %     end
0526     eidors_msg('@@ node2tet  tolerance = %g', opt.tol_node2tet,2);
0527     eidors_msg('@@ edge2edge tolerance = %g', opt.tol_edge2edge,2);
0528     eidors_msg('@@ edge2tri  tolerance = %g', opt.tol_edge2tri,2);
0529    
0530    
0531 %-------------------------------------------------------------------------%
0532 % Perfom unit tests
0533 function do_unit_test
0534    do_case_test;
0535    do_small_test;
0536    do_realistic_test;
0537 
0538    
0539 function do_case_test
0540    ll = eidors_msg('log_level');
0541    eidors_msg('log_level',1);
0542    t1.type = 'fwd_model';
0543    t1.elems = [1 2 3 4];
0544 
0545    X = 2; Y = 3; % subplot matrix
0546    for i = 1:30
0547         fprintf('%d\n',i);
0548         t1.nodes = [0 0 0; 0 1 0; 1 0 0; 0 0 1];
0549         t2 = t1;
0550         switch i
0551            case 1
0552               txt = 'identical';
0553               subplot(X,Y,i), show_test(t1,t2);
0554               c2f = mk_tet_c2f(t1,t2);
0555               unit_test_cmp(txt,c2f,1,eps);
0556            case 2
0557               txt = 'shared face';
0558               t2.nodes(end,end) = -1;
0559               subplot(X,Y,i), show_test(t1,t2);
0560               c2f = mk_tet_c2f(t1,t2);
0561               unit_test_cmp(txt,c2f,0,0);
0562            case 3
0563               txt = 'coplanar faces';
0564               t2.nodes(end,end) = -1;
0565               t1.nodes(:,1:2) = t1.nodes(:,1:2) + .2;
0566               subplot(X,Y,i), show_test(t1,t2);
0567               c2f = mk_tet_c2f(t1,t2);
0568               unit_test_cmp(txt,c2f,0,0);
0569            case 4
0570               txt = 'point on edge';
0571               t2.nodes(:,1) = t1.nodes(:,1) + 1;
0572               t2.nodes(:,2) = t1.nodes(:,2) - .3;
0573               subplot(X,Y,i), show_test(t1,t2);
0574               c2f = mk_tet_c2f(t1,t2);
0575               unit_test_cmp(txt,c2f,0,0);
0576            otherwise
0577              break;
0578         end
0579 
0580       
0581    end
0582    eidors_msg('log_level',ll);
0583 
0584 function do_small_test
0585    fmdl = ng_mk_cyl_models([1 .5],[],[]);
0586    show_fem(fmdl)
0587    v = -.5:.1:.5;
0588    rmdl = mk_grid_model([],v,v,0:.1:1);
0589    hold on
0590    h = show_fem(rmdl);
0591    set(h,'edgecolor','b');
0592    hold off
0593    c2f = mk_tet_c2f(fmdl,rmdl);
0594    tc2f = c2f * rmdl.coarse2fine;
0595    vc2f = mk_grid_c2f(fmdl,rmdl);
0596    unit_test_cmp('mk_tet_c2f v mk_grid_c2f', tc2f,vc2f, 1e-15);
0597 
0598 
0599 function do_realistic_test
0600    fmdl= ng_mk_cyl_models([2,2,.1],[16,1],[.1,0,.025]);
0601    xvec = [-1.5 -.5:.2:.5 1.5];
0602    yvec = [-1.6 -1:.2:1 1.6];
0603    zvec = 0:.25:2;
0604    rmdl = mk_grid_model([],xvec,yvec,zvec);
0605    tic
0606    opt.save_memory = 0;
0607    c2f_a = mk_grid_c2f(fmdl, rmdl,opt);
0608    t = toc;
0609    fprintf('Voxel: t=%f s\n',t);
0610 
0611    tic
0612    opt.save_memory = 0;
0613    c2f_b = mk_tet_c2f(fmdl, rmdl,opt);
0614    t = toc;
0615    fprintf('Tet: t=%f s\n',t);
0616 
0617    c2f_b = c2f_b * rmdl.coarse2fine;
0618    unit_test_cmp('mk_tet_c2f v mk_grid_c2f', c2f_b,c2f_a, 1e-5);
0619 
0620    tic
0621    c2f_n = mk_approx_c2f(fmdl,rmdl);
0622    t = toc;
0623    fprintf('Approximate: t=%f s\n',t);
0624 
0625  function show_test(vox,tet, pts)
0626     show_fem(vox);
0627     hold on
0628     h = show_fem(tet);
0629     set(h, 'EdgeColor','b');
0630     if nargin > 2
0631        plot3(pts(:,1),pts(:,2),pts(:,3),'o');
0632     end
0633     hold off
0634     axis auto