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 6493 2022-12-29 12:51:51Z 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 = point_in_tet(fmdl,rmdl.nodes, opt, true);
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 = point_in_tet(rmdl,fmdl.nodes, opt, true);
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     ftet_in_rtet(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 %-------------------------------------------------------------------------%
0410 % Prepare model
0411 function fmdl = prepare_tet_mdl(fmdl)
0412    fmopt.elem2edge = true;
0413    fmopt.edge2elem = true;
0414    fmopt.face2elem = true;
0415    fmopt.node2elem = true;
0416    fmopt.normals   = true;
0417    fmopt.linear_reorder = false; % this is slow and not needed
0418    ll = eidors_msg('log_level',1);
0419    fmdl = fix_model(fmdl,fmopt);
0420    eidors_msg('log_level',ll);
0421    fmdl.node2elem = logical(fmdl.node2elem);
0422    nElem = size(fmdl.elems,1);
0423    nFace = size(fmdl.faces,1);
0424    fmdl.elem2face = sparse(repmat((1:nElem)',1,4),double(fmdl.elem2face),true,nElem,nFace);
0425 
0426 %-------------------------------------------------------------------------%
0427 % Center scale models
0428 function[fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt)
0429    ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0430    rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0431    fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0432    if isfield(opt,'do_not_scale') && opt.do_not_scale
0433       return
0434    end
0435    maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0436    scale = 1/maxnode;
0437    rmdl.nodes = scale*rmdl.nodes;
0438    fmdl.nodes = scale*fmdl.nodes;
0439    eidors_msg('@@ models scaled by %g', scale,2);
0440 
0441 %-------------------------------------------------------------------------%
0442 % Remove obviously non-overlapping parts of the models
0443 function [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl)
0444    f_min = min(fmdl.nodes);
0445    f_max = max(fmdl.nodes);
0446    r_min = min(rmdl.nodes);
0447    r_max = max(rmdl.nodes);
0448    
0449    % nodes outside the bounding box of the other model
0450    f_gt  = bsxfun(@gt, fmdl.nodes, r_max);
0451    f_lt  = bsxfun(@lt, fmdl.nodes, r_min);
0452    r_gt  = bsxfun(@gt, rmdl.nodes, f_max);
0453    r_lt  = bsxfun(@lt, rmdl.nodes, f_min);
0454    
0455    % elems outside the bounding box of the other model
0456    re_gt = any(reshape(all(reshape(r_gt(rmdl.elems',:),4,[])),[],3),2);
0457    re_lt = any(reshape(all(reshape(r_lt(rmdl.elems',:),4,[])),[],3),2);
0458    fe_gt = any(reshape(all(reshape(f_gt(fmdl.elems',:),4,[])),[],3),2);
0459    fe_lt = any(reshape(all(reshape(f_lt(fmdl.elems',:),4,[])),[],3),2);
0460    
0461    % elems to keep
0462    rmdl_idx = ~(re_gt | re_lt);
0463    fmdl_idx = ~(fe_gt | fe_lt);
0464    
0465    % remove non-overlapping elems
0466    rmdl.elems = rmdl.elems(rmdl_idx,:);
0467    fmdl.elems = fmdl.elems(fmdl_idx,:);
0468    
0469    % remove unused nodes
0470    [r_used_nodes,jnk,r_n] = unique(rmdl.elems(:));
0471    [f_used_nodes,jnk,f_n] = unique(fmdl.elems(:));
0472    
0473    r_idx = 1:numel(r_used_nodes);
0474    f_idx = 1:numel(f_used_nodes);
0475    
0476    rmdl.elems = reshape(r_idx(r_n),[],4);
0477    fmdl.elems = reshape(f_idx(f_n),[],4);
0478    
0479    rmdl.nodes = rmdl.nodes(r_used_nodes,:);
0480    fmdl.nodes = fmdl.nodes(f_used_nodes,:);
0481    
0482    % for the benefit of any (debug) plots later on
0483    try, rmdl = rmfield(rmdl,'boundary'); end
0484    try, fmdl = rmfield(fmdl,'boundary'); end
0485     
0486 %-------------------------------------------------------------------------%
0487 % Parse option struct
0488  function opt = parse_opts(fmdl,rmdl, opt)
0489 
0490     
0491     if ~isfield(opt, 'tol_node2tet');
0492         opt.tol_node2tet = eps; % * max(rmdl_rng,fmdl_rng)^3;
0493     end
0494     if ~isfield(opt, 'tol_edge2edge')
0495         opt.tol_edge2edge = 6*eps(min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))));
0496     end
0497     if ~isfield(opt, 'tol_edge2tri')
0498         opt.tol_edge2tri = eps; %1e-10
0499     end
0500 %     if ~isfield(opt, 'save_memory')
0501 %        opt.save_memory = 0;
0502 %     end
0503     eidors_msg('@@ node2tet  tolerance = %g', opt.tol_node2tet,2);
0504     eidors_msg('@@ edge2edge tolerance = %g', opt.tol_edge2edge,2);
0505     eidors_msg('@@ edge2tri  tolerance = %g', opt.tol_edge2tri,2);
0506    
0507    
0508 %-------------------------------------------------------------------------%
0509 % Perfom unit tests
0510 function do_unit_test
0511    do_case_test;
0512 %  do_small_test;
0513 %  do_realistic_test;
0514 
0515    
0516 function do_case_test
0517    ll = eidors_msg('log_level');
0518    eidors_msg('log_level',1);
0519    t1.type = 'fwd_model';
0520    t1.elems = [1 2 3 4];
0521 
0522    X = 2; Y = 3; % subplot matrix
0523    for i = 1:30
0524         fprintf('%d\n',i);
0525         t1.nodes = [0 0 0; 0 1 0; 1 0 0; 0 0 1];
0526         t2 = t1;
0527         switch i
0528            case 1
0529               txt = 'identical';
0530               subplot(X,Y,i), show_test(t1,t2);
0531               c2f = mk_tet_c2f(t1,t2);
0532               unit_test_cmp(txt,c2f,1,eps);
0533            case 2
0534               txt = 'shared face';
0535               t2.nodes(end,end) = -1;
0536               subplot(X,Y,i), show_test(t1,t2);
0537               c2f = mk_tet_c2f(t1,t2);
0538               unit_test_cmp(txt,c2f,0,0);
0539            case 3
0540               txt = 'coplanar faces';
0541               t2.nodes(end,end) = -1;
0542               t1.nodes(:,1:2) = t1.nodes(:,1:2) + .2;
0543               subplot(X,Y,i), show_test(t1,t2);
0544               c2f = mk_tet_c2f(t1,t2);
0545               unit_test_cmp(txt,c2f,0,0);
0546            case 4
0547               txt = 'point on edge';
0548               t2.nodes(:,1) = t1.nodes(:,1) + 1;
0549               t2.nodes(:,2) = t1.nodes(:,2) - .3;
0550               subplot(X,Y,i), show_test(t1,t2);
0551               c2f = mk_tet_c2f(t1,t2);
0552               unit_test_cmp(txt,c2f,0,0);
0553            case 5
0554               txt = 'tet split into four';
0555               t2.elems = [1 2 3 5;1 2 4 5;1 3 4 5;2 3 4 5];
0556               t2.nodes = [0 0 0;0 5 0;5 0 0;0 0 5;1 1 1]/5;
0557               subplot(X,Y,i), show_test(t1,t2);
0558               c2f = mk_tet_c2f(t1,t2);
0559               unit_test_cmp(txt,c2f*10,[2,2,2,4],10*eps);
0560            otherwise
0561              break;
0562         end
0563 
0564       
0565    end
0566    eidors_msg('log_level',ll);
0567 
0568 function do_small_test
0569    ll = eidors_msg('log_level');
0570    eidors_msg('log_level',1);
0571    fmdl = ng_mk_cyl_models([1 .5],[],[]);
0572    show_fem(fmdl)
0573    v = -.5:.1:.5;
0574    rmdl = mk_grid_model([],v,v,0:.1:1);
0575    hold on
0576    h = show_fem(rmdl);
0577    set(h,'edgecolor','b');
0578    hold off
0579    c2f = mk_tet_c2f(fmdl,rmdl);
0580    tc2f = c2f * rmdl.coarse2fine;
0581    vc2f = mk_grid_c2f(fmdl,rmdl);
0582    unit_test_cmp('mk_tet_c2f v mk_grid_c2f', tc2f,vc2f, 1e-15);
0583    eidors_msg('log_level',ll);
0584 
0585 
0586 function do_realistic_test
0587    ll = eidors_msg('log_level');
0588    eidors_msg('log_level',1);
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    eidors_msg('log_level',ll);
0614 
0615  function show_test(vox,tet, pts)
0616     show_fem(vox);
0617     hold on
0618     h = show_fem(tet);
0619     set(h, 'EdgeColor','b');
0620     if nargin > 2
0621        plot3(pts(:,1),pts(:,2),pts(:,3),'o');
0622     end
0623     hold off
0624     axis auto
Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005