mk_tri_c2f

PURPOSE ^

MK_TRI_C2F - calculate a coarse2fine mapping for triangle-based models.

SYNOPSIS ^

function c2f = mk_tri_c2f(fmdl,rmdl,opt)

DESCRIPTION ^

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

 C2F = MK_TRI_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
                       square before any calculations, including 
                       thresholds. Default: false
      .tol_node2tri  - 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([0,2,.2],[],[]);
     rmdl = ng_mk_cyl_models([0,2],[],[]);
     c2f  = mk_tri_c2f(fmdl,rmdl);
     h = show_fem(fmdl); set(h,'LineWidth',0.1,'FaceColor','none')
     hold on
     h = show_fem(rmdl); set(h,'EdgeColor','b','LineWidth',2,...
        'FaceColor','none');
     hold off

 See also MK_GRID_C2F, FIND_EDGE2EDGE_INTERSECTIONS, CONVHULLN
          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_tri_c2f(fmdl,rmdl,opt)
0002 %MK_TRI_C2F - calculate a coarse2fine mapping for triangle-based models.
0003 % C2F = MK_TRI_C2F(FMDL,RMDL) returns in C2F the fraction of area 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 area defined by a set of intersection
0007 % points between individual triangles.
0008 %
0009 % C2F = MK_TRI_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 %                       square before any calculations, including
0017 %                       thresholds. Default: false
0018 %      .tol_node2tri  - minimum value of a barycentric coordinate to
0019 %                       decide a point is lying inside a triangle and not
0020 %                       on its edge. Default: eps
0021 %
0022 % NOTE that for grid-based models, such as returned by MK_GRID_MODEL or
0023 % MK_VOXEL_VOLUME, MK_GRID_C2F is much faster.
0024 %
0025 % Set eidors_msg 'log level' < 2 to supress output to command line.
0026 %
0027 % Examples:
0028 %     fmdl = ng_mk_cyl_models([0,2,.2],[],[]);
0029 %     rmdl = ng_mk_cyl_models([0,2],[],[]);
0030 %     c2f  = mk_tri_c2f(fmdl,rmdl);
0031 %     h = show_fem(fmdl); set(h,'LineWidth',0.1,'FaceColor','none')
0032 %     hold on
0033 %     h = show_fem(rmdl); set(h,'EdgeColor','b','LineWidth',2,...
0034 %        'FaceColor','none');
0035 %     hold off
0036 %
0037 % See also MK_GRID_C2F, FIND_EDGE2EDGE_INTERSECTIONS, CONVHULLN
0038 %          MK_APPROX_C2F, POINT_IN_TRIANGLE, EIDORS_MSG
0039 
0040 % (C) 2015 Bartlomiej Grychtol
0041 % License: GPL version 2 or 3
0042 % $Id: mk_tri_c2f.m 5068 2015-05-30 12:07:52Z bgrychtol $
0043 
0044 
0045 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'), do_unit_test; return; end
0046 if nargin < 3
0047    opt = struct();
0048 end
0049 
0050 f_elems = size(fmdl.elems,1);
0051 r_elems = size(rmdl.elems,1);
0052 
0053 c2f = sparse(f_elems,r_elems);
0054 [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl);
0055 
0056 if ~(any(fmdl_idx) && any(rmdl_idx))
0057    eidors_msg('@@: models do not overlap, returning all-zeros');
0058    return
0059 end
0060 
0061 [fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt);
0062 
0063 opt = parse_opts(fmdl,rmdl, opt);
0064 
0065 
0066 copt.fstr = 'mk_tri_c2f';
0067 
0068 c2f(fmdl_idx,rmdl_idx) = eidors_cache(@do_mk_tri_c2f,{fmdl,rmdl,opt},copt);
0069 end
0070 
0071 function c2f = do_mk_tri_c2f(fmdl,rmdl,opt)
0072    DEBUG = eidors_debug('query','mk_tri_c2f');
0073    
0074    c2f = sparse(0,0);
0075    progress_msg('Prepare fine model...');
0076    fmdl = prepare_tri_mdl(fmdl);
0077    progress_msg(Inf);
0078    
0079    progress_msg('Prepare course model...');
0080    rmdl = prepare_tri_mdl(rmdl);
0081    progress_msg(Inf);
0082    
0083    progress_msg('Find edge2edge intersections...');
0084    [intpts,fedge2redge,fedge2intpts,redge2intpts] = ...
0085       find_edge2edge_intersections( fmdl.edges,fmdl.nodes,...
0086                                     rmdl.edges,rmdl.nodes);
0087    progress_msg(sprintf('Found %d',size(intpts,1)),Inf);
0088 
0089    
0090    progress_msg('Find c_nodes in f_elems...');
0091    rnode2ftri = point_in_triangle(rmdl.nodes, fmdl.elems, fmdl.nodes, opt.tol_node2tri);
0092    progress_msg(sprintf('Found %d', nnz(rnode2ftri)), Inf);
0093    
0094    progress_msg('Find c_elems in f_elems...')
0095    rtri_in_ftri = (double(rmdl.node2elem') * rnode2ftri) == 3;
0096    progress_msg(sprintf('Found %d',nnz(rtri_in_ftri)), Inf);
0097    
0098    progress_msg('Find f_nodes in c_elems...');
0099    fnode2rtri = point_in_triangle(fmdl.nodes, rmdl.elems, rmdl.nodes, opt.tol_node2tri);
0100    progress_msg(sprintf('Found %d', nnz(fnode2rtri)), Inf);
0101    
0102    progress_msg('Find f_elems in c_elems...')
0103    ftri_in_rtri = (double(fmdl.node2elem') * fnode2rtri) == 3;
0104    progress_msg(sprintf('Found %d',nnz(ftri_in_rtri)), Inf);
0105 
0106    progress_msg('Find total intersections...');
0107    rtri2ftri = double(rmdl.edge2elem') * fedge2redge' * fmdl.edge2elem;
0108    % exclude inclusion (dealt with separately)
0109    rtri2ftri = rtri2ftri & ~rtri_in_ftri & ~ftri_in_rtri'; 
0110    progress_msg(sprintf('Found %d',nnz(rtri2ftri)), Inf);
0111    
0112    progress_msg('Calculate intersection volumes...');
0113    % sparse logical multiplication doesn't exist
0114    rtri2intpt = logical(rmdl.edge2elem'*redge2intpts)';
0115    ftri2intpt = logical(fmdl.edge2elem'*fedge2intpts)';
0116    
0117    rtri_todo = find(sum(rtri2ftri,2)>0);
0118    C = []; F = []; V = [];
0119 
0120    id = 0; N = length(rtri_todo);
0121    mint = ceil(N/100);
0122    for v = rtri_todo'
0123       id = id+1;
0124       if mod(id,mint)==0, progress_msg(id/N); end
0125       tri_todo = find(rtri2ftri(v,:));
0126 
0127       common_intpts = bsxfun(@and,rtri2intpt(:,v), ftri2intpt(:,tri_todo));
0128       
0129       f_nodes     = bsxfun(@and,fnode2rtri(:,v), fmdl.node2elem(:,tri_todo));
0130       r_nodes     = bsxfun(@and,rnode2ftri(:,tri_todo), rmdl.node2elem(:,v));
0131       
0132       C = [C; v*ones(numel(tri_todo),1)];
0133       F = [F; tri_todo'];
0134       last_v = numel(V);
0135       V = [V; zeros(numel(tri_todo),1)]; % pre-allocate
0136       
0137       for t = 1:numel(tri_todo)
0138          pts = [ intpts(common_intpts(:,t),:);
0139                   fmdl.nodes(f_nodes(:,t),:);
0140                   rmdl.nodes(r_nodes(:,t),:)];
0141          last_v = last_v + 1;
0142          if size(pts,1) < 3, continue, end 
0143          try
0144             % move points to origin (helps for small elements at
0145             % large coordinates
0146             ctr = mean(pts);
0147             pts = bsxfun(@minus,pts,ctr);
0148             scale = max(abs(pts(:)));
0149             if scale == 0 %happens when there's only one point
0150                continue
0151             end
0152             % scale largest coordinate to 1 (helps with precision)
0153             pts = pts ./ scale;
0154             % force thorough search for initinal simplex and
0155             % supress precision warnings
0156             [K, V(last_v)] = convhulln(pts,{'Qt Pp Qs'});
0157             V(last_v) = max(V(last_v),0); % numerical issues may produce tiny negative volume
0158             V(last_v) = V(last_v) * scale^2; % undo scaling
0159          catch err
0160             ok = false;
0161             switch err.identifier
0162                case {'MATLAB:qhullmx:DegenerateData', 'MATLAB:qhullmx:UndefinedError'}
0163                   % border case point may be included multiple times.
0164                   % this is OK... though we may miss cases where more
0165                   % points should have been found but were not
0166                   u = uniquetol(pts*scale,1e-14,'ByRows',true,'DataScale', 1);
0167                   ok = ok | (size(u,1) < 3);
0168                   if ~ok
0169                      % test for colinearity
0170                      cp = bsxfun(@minus, u(2:end,:), u(1,:));
0171                      l = sqrt(sum(cp.^2,2));
0172                      cp = bsxfun(@rdivide, cp, l);
0173                      u = uniquetol(cp,1e-14,'ByRows',true,'DataScale',1);
0174                      ok = ok | size(u,1) == 1; % co-linear points
0175                   end
0176             end
0177             if ~ok
0178                if DEBUG || eidors_debug('query','mk_tet_c2f:convhulln');
0179                   tri.nodes = fmdl.nodes;
0180                   vox.nodes = rmdl.nodes;
0181                   tri.type = 'fwd_model';
0182                   vox.type = 'fwd_model';
0183                   vox.elems = rmdl.elems(v,:);
0184                   vox.boundary = vox.elems;
0185                   tri.elems = fmdl.elems(tri_todo(t),:);
0186                   clf
0187                   show_fem(vox)
0188                   hold on
0189                   h = show_fem(tri);
0190                   set(h,'EdgeColor','b')
0191                   pts = bsxfun(@plus,pts*scale,ctr);
0192                   plot(pts(:,1),pts(:,2),'o');
0193                   hold off
0194                   axis auto
0195                   keyboard
0196                else
0197                   fprintf('\n');
0198                   eidors_msg(['convhulln has thrown an error. ' ...
0199                      'Enable eidors_debug on mk_tri_c2f and re-run to see a debug plot'],0);
0200                   rethrow(err);
0201                end
0202             end
0203          end
0204          
0205       end
0206    end
0207    progress_msg(Inf);
0208    
0209    c2f = sparse(F,C,V,size(fmdl.elems,1),size(rmdl.elems,1));
0210    
0211    % add rtri contained in ftri
0212    try rmdl = rmfield(rmdl,'coarse2fine'); end % messes with volume
0213    c2f = c2f + bsxfun(@times, sparse(rtri_in_ftri), get_elem_volume(rmdl))';
0214    
0215    % normalize to fine volume
0216    vol = get_elem_volume(fmdl);
0217    c2f = bsxfun(@rdivide,c2f,vol);
0218    
0219    % count identical triangles just once
0220    ftri_in_rtri(rtri_in_ftri') = 0;
0221    
0222    % add fine elems contained in coarse elems
0223    c2f = c2f + ftri_in_rtri;
0224 
0225 end
0226 
0227 
0228 function [pts,FE2CE,FE2pts,CE2pts] = find_edge2edge_intersections(FE,FN,CE,CN)
0229    P1 = FN(FE(:,1),:);
0230    P2 = FN(FE(:,2),:);
0231    P3 = CN(CE(:,1),:);
0232    P4 = CN(CE(:,2),:);
0233    
0234    P21 = P2 - P1;
0235    P43 = P4 - P3;
0236    
0237    invden = ( bsxfun(@times, P21(:,1), P43(:,2)') - ...
0238               bsxfun(@times, P21(:,2), P43(:,1)')       ).^-1;
0239    P13_x = bsxfun(@minus,P1(:,1),P3(:,1)');
0240    P13_y = bsxfun(@minus,P1(:,2),P3(:,2)');
0241    s = ( bsxfun(@times,-P21(:,2), P13_x) + ...
0242          bsxfun(@times, P21(:,1), P13_y)) .* invden;
0243    t = ( bsxfun(@times,-P43(:,2)',P13_x) + ...
0244          bsxfun(@times, P43(:,1)',P13_y)) .* invden;
0245 
0246    FE2CE= s >= 0 & s <= 1 & t >= 0 & t <= 1;
0247    
0248    [fe, ce] = find(FE2CE);
0249    N_pts = size(fe,1);
0250    pts = zeros(N_pts,2);
0251    for i = 1:N_pts
0252       pts(i,:) = P1(fe(i),:) + t(fe(i),ce(i)) * P21(fe(i),:);
0253    end
0254    FE2CE = sparse(FE2CE);
0255    N_ce = size(CE,1);
0256    N_fe = size(FE,1);
0257    FE2pts = sparse(fe, 1:N_pts, ones(N_pts,1), N_fe, N_pts);
0258    CE2pts = sparse(ce, 1:N_pts, ones(N_pts,1), N_ce, N_pts);
0259 
0260 
0261 end
0262 
0263 %-------------------------------------------------------------------------%
0264 % Prepare model
0265 function fmdl = prepare_tri_mdl(fmdl)
0266    fmopt.elem2edge = true;
0267    fmopt.edge2elem = true;
0268 %    fmopt.face2elem = true;
0269    fmopt.node2elem = true;
0270 %    fmopt.normals   = true;
0271    fmopt.linear_reorder = false; % this is slow and not needed
0272    ll = eidors_msg('log_level',1);
0273    fmdl = fix_model(fmdl,fmopt);
0274    eidors_msg('log_level',ll);
0275    fmdl.node2elem = logical(fmdl.node2elem);
0276 
0277 end
0278 
0279 %-------------------------------------------------------------------------%
0280 % Remove obviously non-overlapping parts of the models
0281 function [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl)
0282 
0283 % instead of the usual approach, we could calculate the convex hull and
0284 % then use inpolygon... but that would require trusting inpolygon
0285 
0286    f_min = min(fmdl.nodes);
0287    f_max = max(fmdl.nodes);
0288    r_min = min(rmdl.nodes);
0289    r_max = max(rmdl.nodes);
0290    
0291    % nodes outside the bounding box of the other model
0292    f_gt  = bsxfun(@gt, fmdl.nodes, r_max);
0293    f_lt  = bsxfun(@lt, fmdl.nodes, r_min);
0294    r_gt  = bsxfun(@gt, rmdl.nodes, f_max);
0295    r_lt  = bsxfun(@lt, rmdl.nodes, f_min);
0296 
0297    % elems outside the bounding box of the other model
0298    re_gt = any(reshape(all(reshape(r_gt(rmdl.elems',:),3,[])),[],2),2);
0299    re_lt = any(reshape(all(reshape(r_lt(rmdl.elems',:),3,[])),[],2),2);
0300    fe_gt = any(reshape(all(reshape(f_gt(fmdl.elems',:),3,[])),[],2),2);
0301    fe_lt = any(reshape(all(reshape(f_lt(fmdl.elems',:),3,[])),[],2),2);
0302    
0303    % elems to keep
0304    rmdl_idx = ~(re_gt | re_lt);
0305    fmdl_idx = ~(fe_gt | fe_lt);
0306    
0307    % remove non-overlapping elems
0308    rmdl.elems = rmdl.elems(rmdl_idx,:);
0309    fmdl.elems = fmdl.elems(fmdl_idx,:);
0310    
0311    % remove unused nodes
0312    [r_used_nodes,jnk,r_n] = unique(rmdl.elems(:));
0313    [f_used_nodes,jnk,f_n] = unique(fmdl.elems(:));
0314    
0315    r_idx = 1:numel(r_used_nodes);
0316    f_idx = 1:numel(f_used_nodes);
0317    
0318    rmdl.elems = reshape(r_idx(r_n),[],3);
0319    fmdl.elems = reshape(f_idx(f_n),[],3);
0320    
0321    rmdl.nodes = rmdl.nodes(r_used_nodes,:);
0322    fmdl.nodes = fmdl.nodes(f_used_nodes,:);
0323    
0324    % for the benefit of any (debug) plots later on
0325    try, rmdl = rmfield(rmdl,'boundary'); end
0326    try, fmdl = rmfield(fmdl,'boundary'); end
0327     
0328 end
0329 
0330 %-------------------------------------------------------------------------%
0331 % Center scale models
0332 function[fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt)
0333    ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0334    rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0335    fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0336    if isfield(opt,'do_not_scale') && opt.do_not_scale
0337       return
0338    end
0339    maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0340    scale = 1/maxnode;
0341    rmdl.nodes = scale*rmdl.nodes;
0342    fmdl.nodes = scale*fmdl.nodes;
0343    eidors_msg('@@ models scaled by %g', scale,2);
0344 end
0345 
0346 %-------------------------------------------------------------------------%
0347 % Parse option struct
0348  function opt = parse_opts(fmdl,rmdl, opt)
0349   
0350     if ~isfield(opt, 'tol_node2tri');
0351         opt.tol_node2tri = eps; % * max(rmdl_rng,fmdl_rng)^3;
0352     end
0353 
0354     eidors_msg('@@ node2tri  tolerance = %g', opt.tol_node2tri,2);
0355 
0356  end   
0357 
0358 function do_unit_test
0359    do_small_test
0360 end
0361 
0362 function do_small_test
0363    imdl = mk_common_model('a2c',16);
0364    rmdl = imdl.fwd_model;
0365 %    rmdl = ng_mk_cyl_models([0 .5],[],[]);
0366 %    rmdl.nodes(:,1) = rmdl.nodes(:,1) + .5;
0367 %    show_fem(rmdl,[0 0 1])
0368 %    hold all
0369    imdl = mk_common_model('d2c',16);
0370    fmdl = imdl.fwd_model;
0371 %    fmdl = ng_mk_cyl_models([0 .5 .1],[],[]);
0372 %    h = show_fem(fmdl);
0373 %    set(h,'edgecolor','b','facecolor','none');
0374 %    hold off
0375 %    axis auto
0376    c2f = mk_tri_c2f(fmdl,rmdl);
0377    
0378    
0379    clf
0380    img1 = mk_image(fmdl,0);
0381    img1.elem_data = sum(c2f,2);
0382    img1.calc_colous.ref_level = .5;
0383    img1.calc_colours.clim = .5;
0384    
0385    subplot(121)
0386    show_fem(img1);
0387    img2 = mk_image(rmdl,0);
0388    img2.elem_data = (c2f' * get_elem_volume(fmdl)) ./ get_elem_volume(rmdl);
0389    img2.calc_colous.ref_level = .5;
0390    img2.calc_colours.clim = .55;
0391    subplot(122)
0392    show_fem(img2);
0393    
0394    unit_test_cmp('Check C2F size', size(c2f),[length(fmdl.elems), length(rmdl.elems)]);
0395    unit_test_cmp('Check C2F max==1', max(c2f(:)), 1);
0396    unit_test_cmp('Check C2F min==0', min(c2f(:)), 0);
0397    
0398    f2c = bsxfun(@rdivide, bsxfun(@times, c2f, get_elem_volume(fmdl))', get_elem_volume(rmdl));
0399    unit_test_cmp('Check F2C max==1', max(sum(f2c,2)), 1, 1e-15);
0400    unit_test_cmp('Check F2C min==0', min(f2c(:)), 0);
0401 end

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