mk_tri_c2f

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

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 6152 2021-10-18 12:51:37Z aadler $
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             % move points to origin (helps for small elements at
0144             % large coordinates
0145             ctr = mean(pts);
0146             pts = bsxfun(@minus,pts,ctr);
0147             scale = max(abs(pts(:)));
0148             if scale == 0 %happens when there's only one point
0149                continue
0150             end
0151             % scale largest coordinate to 1 (helps with precision)
0152             pts = pts ./ scale;
0153             % force thorough search for initinal simplex and
0154             % supress precision warnings
0155 %           if any(any(dist2(pts))); continue; end
0156          pts= uniquetol(pts,1e-10,'ByRows',true,'DataScale',1);
0157          % This won't catch all cases
0158          if size(pts,1)<=size(pts,2); continue; end
0159          if any(std(pts)<1e-14); continue; end
0160          try
0161             [K, V(last_v)] = convhulln(pts,{'Qt Pp Qs'});
0162 
0163             V(last_v) = max(V(last_v),0); % numerical issues may produce tiny negative volume
0164             V(last_v) = V(last_v) * scale^2; % undo scaling
0165          catch err
0166 %    Save cases were errors called
0167 %           if ~exist('ok'); ptsi={}; end
0168 %           ptsi{end+1} = pts;
0169 %           save -mat CHP.mat ptsi;
0170             ok = false;
0171             if exist('OCTAVE_VERSION')
0172                if strcmp(err.message,'convhulln: qhull failed')
0173                   err.identifier =  'MATLAB:qhullmx:DegenerateData';
0174                end
0175                   
0176             end
0177             switch err.identifier
0178                case {'MATLAB:qhullmx:DegenerateData', 'MATLAB:qhullmx:UndefinedError'}
0179                   % border case point may be included multiple times.
0180                   % this is OK... though we may miss cases where more
0181                   % points should have been found but were not
0182                   u = uniquetol(pts*scale,1e-14,'ByRows',true,'DataScale', 1);
0183                   ok = ok | (size(u,1) < 3);
0184                   if ~ok
0185                      % test for colinearity
0186                      cp = bsxfun(@minus, u(2:end,:), u(1,:));
0187                      l = sqrt(sum(cp.^2,2));
0188                      cp = bsxfun(@rdivide, cp, l);
0189                      u = uniquetol(cp,1e-14,'ByRows',true,'DataScale',1);
0190                      ok = ok | size(u,1) == 1; % co-linear points
0191                   end
0192             end
0193             if ~ok
0194                if DEBUG || eidors_debug('query','mk_tet_c2f:convhulln');
0195                   tri.nodes = fmdl.nodes;
0196                   vox.nodes = rmdl.nodes;
0197                   tri.type = 'fwd_model';
0198                   vox.type = 'fwd_model';
0199                   vox.elems = rmdl.elems(v,:);
0200                   vox.boundary = vox.elems;
0201                   tri.elems = fmdl.elems(tri_todo(t),:);
0202                   clf
0203                   show_fem(vox)
0204                   hold on
0205                   h = show_fem(tri);
0206                   set(h,'EdgeColor','b')
0207                   pts = bsxfun(@plus,pts*scale,ctr);
0208                   plot(pts(:,1),pts(:,2),'o');
0209                   hold off
0210                   axis auto
0211                   keyboard
0212                else
0213                   fprintf('\n');
0214                   eidors_msg(['convhulln has thrown an error. (',e.message,')', ...
0215                      'Enable eidors_debug on mk_tri_c2f and re-run to see a debug plot'],0);
0216                   rethrow(err);
0217                end
0218             end
0219          end
0220          
0221       end
0222    end
0223    progress_msg(Inf);
0224    
0225    c2f = sparse(F,C,V,size(fmdl.elems,1),size(rmdl.elems,1));
0226    
0227    % add rtri contained in ftri
0228    try rmdl = rmfield(rmdl,'coarse2fine'); end % messes with volume
0229    c2f = c2f + bsxfun(@times, sparse(rtri_in_ftri), get_elem_volume(rmdl))';
0230    
0231    % normalize to fine volume
0232    vol = get_elem_volume(fmdl,-inf);
0233    c2f = bsxfun(@rdivide,c2f,vol);
0234    
0235    % count identical triangles just once
0236    ftri_in_rtri(rtri_in_ftri') = 0;
0237    
0238    % add fine elems contained in coarse elems
0239    c2f = c2f + ftri_in_rtri;
0240 
0241 end
0242 
0243 
0244 function [pts,FE2CE,FE2pts,CE2pts] = find_edge2edge_intersections(FE,FN,CE,CN)
0245    P1 = FN(FE(:,1),:);
0246    P2 = FN(FE(:,2),:);
0247    P3 = CN(CE(:,1),:);
0248    P4 = CN(CE(:,2),:);
0249    
0250    P21 = P2 - P1;
0251    P43 = P4 - P3;
0252    
0253    invden = ( bsxfun(@times, P21(:,1), P43(:,2)') - ...
0254               bsxfun(@times, P21(:,2), P43(:,1)')       ).^-1;
0255    P13_x = bsxfun(@minus,P1(:,1),P3(:,1)');
0256    P13_y = bsxfun(@minus,P1(:,2),P3(:,2)');
0257    s = ( bsxfun(@times,-P21(:,2), P13_x) + ...
0258          bsxfun(@times, P21(:,1), P13_y)) .* invden;
0259    t = ( bsxfun(@times,-P43(:,2)',P13_x) + ...
0260          bsxfun(@times, P43(:,1)',P13_y)) .* invden;
0261 
0262    FE2CE= s >= 0 & s <= 1 & t >= 0 & t <= 1;
0263    
0264    [fe, ce] = find(FE2CE);
0265    N_pts = size(fe,1);
0266    pts = zeros(N_pts,2);
0267    for i = 1:N_pts
0268       pts(i,:) = P1(fe(i),:) + t(fe(i),ce(i)) * P21(fe(i),:);
0269    end
0270    FE2CE = sparse(FE2CE);
0271    N_ce = size(CE,1);
0272    N_fe = size(FE,1);
0273    FE2pts = sparse(fe, 1:N_pts, ones(N_pts,1), N_fe, N_pts);
0274    CE2pts = sparse(ce, 1:N_pts, ones(N_pts,1), N_ce, N_pts);
0275 
0276 
0277 end
0278 
0279 %-------------------------------------------------------------------------%
0280 % Prepare model
0281 function fmdl = prepare_tri_mdl(fmdl)
0282    fmopt.elem2edge = true;
0283    fmopt.edge2elem = true;
0284 %    fmopt.face2elem = true;
0285    fmopt.node2elem = true;
0286 %    fmopt.normals   = true;
0287    fmopt.linear_reorder = false; % this is slow and not needed
0288    ll = eidors_msg('log_level',1);
0289    fmdl = fix_model(fmdl,fmopt);
0290    eidors_msg('log_level',ll);
0291    fmdl.node2elem = logical(fmdl.node2elem);
0292 
0293 end
0294 
0295 %-------------------------------------------------------------------------%
0296 % Remove obviously non-overlapping parts of the models
0297 function [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl)
0298 
0299 % instead of the usual approach, we could calculate the convex hull and
0300 % then use inpolygon... but that would require trusting inpolygon
0301 
0302    f_min = min(fmdl.nodes);
0303    f_max = max(fmdl.nodes);
0304    r_min = min(rmdl.nodes);
0305    r_max = max(rmdl.nodes);
0306    
0307    % nodes outside the bounding box of the other model
0308    f_gt  = bsxfun(@gt, fmdl.nodes, r_max);
0309    f_lt  = bsxfun(@lt, fmdl.nodes, r_min);
0310    r_gt  = bsxfun(@gt, rmdl.nodes, f_max);
0311    r_lt  = bsxfun(@lt, rmdl.nodes, f_min);
0312 
0313    % elems outside the bounding box of the other model
0314    re_gt = any(reshape(all(reshape(r_gt(rmdl.elems',:),3,[])),[],2),2);
0315    re_lt = any(reshape(all(reshape(r_lt(rmdl.elems',:),3,[])),[],2),2);
0316    fe_gt = any(reshape(all(reshape(f_gt(fmdl.elems',:),3,[])),[],2),2);
0317    fe_lt = any(reshape(all(reshape(f_lt(fmdl.elems',:),3,[])),[],2),2);
0318    
0319    % elems to keep
0320    rmdl_idx = ~(re_gt | re_lt);
0321    fmdl_idx = ~(fe_gt | fe_lt);
0322    
0323    % remove non-overlapping elems
0324    rmdl.elems = rmdl.elems(rmdl_idx,:);
0325    fmdl.elems = fmdl.elems(fmdl_idx,:);
0326    
0327    % remove unused nodes
0328    [r_used_nodes,jnk,r_n] = unique(rmdl.elems(:));
0329    [f_used_nodes,jnk,f_n] = unique(fmdl.elems(:));
0330    
0331    r_idx = 1:numel(r_used_nodes);
0332    f_idx = 1:numel(f_used_nodes);
0333    
0334    rmdl.elems = reshape(r_idx(r_n),[],3);
0335    fmdl.elems = reshape(f_idx(f_n),[],3);
0336    
0337    rmdl.nodes = rmdl.nodes(r_used_nodes,:);
0338    fmdl.nodes = fmdl.nodes(f_used_nodes,:);
0339    
0340    % for the benefit of any (debug) plots later on
0341    try, rmdl = rmfield(rmdl,'boundary'); end
0342    try, fmdl = rmfield(fmdl,'boundary'); end
0343     
0344 end
0345 
0346 %-------------------------------------------------------------------------%
0347 % Center scale models
0348 function[fmdl,rmdl] = center_scale_models(fmdl,rmdl, opt)
0349    ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0350    rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0351    fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0352    if isfield(opt,'do_not_scale') && opt.do_not_scale
0353       return
0354    end
0355    maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0356    scale = 1/maxnode;
0357    rmdl.nodes = scale*rmdl.nodes;
0358    fmdl.nodes = scale*fmdl.nodes;
0359    eidors_msg('@@ models scaled by %g', scale,2);
0360 end
0361 
0362 %-------------------------------------------------------------------------%
0363 % Parse option struct
0364  function opt = parse_opts(fmdl,rmdl, opt)
0365   
0366     if ~isfield(opt, 'tol_node2tri');
0367         opt.tol_node2tri = eps; % * max(rmdl_rng,fmdl_rng)^3;
0368     end
0369 
0370     eidors_msg('@@ node2tri  tolerance = %g', opt.tol_node2tri,2);
0371 
0372  end   
0373 
0374 function do_unit_test
0375    do_small_test
0376 end
0377 
0378 function do_small_test
0379    imdl = mk_common_model('a2c',16);
0380    rmdl = imdl.fwd_model;
0381 %    rmdl = ng_mk_cyl_models([0 .5],[],[]);
0382 %    rmdl.nodes(:,1) = rmdl.nodes(:,1) + .5;
0383 %    show_fem(rmdl,[0 0 1])
0384 %    hold all
0385    imdl = mk_common_model('d2c',16);
0386    fmdl = imdl.fwd_model;
0387 %    fmdl = ng_mk_cyl_models([0 .5 .1],[],[]);
0388 %    h = show_fem(fmdl);
0389 %    set(h,'edgecolor','b','facecolor','none');
0390 %    hold off
0391 %    axis auto
0392    c2f = mk_tri_c2f(fmdl,rmdl);
0393    
0394    
0395    clf
0396    img1 = mk_image(fmdl,0);
0397    img1.elem_data = sum(c2f,2);
0398    img1.calc_colous.ref_level = .5;
0399    img1.calc_colours.clim = .5;
0400    
0401    subplot(121)
0402    show_fem(img1);
0403    img2 = mk_image(rmdl,0);
0404    img2.elem_data = (c2f' * get_elem_volume(fmdl)) ./ get_elem_volume(rmdl);
0405    img2.calc_colous.ref_level = .5;
0406    img2.calc_colours.clim = .55;
0407    subplot(122)
0408    show_fem(img2);
0409    
0410    unit_test_cmp('Check C2F size', size(c2f),[length(fmdl.elems), length(rmdl.elems)]);
0411    unit_test_cmp('Check C2F max==1', max(c2f(:)), 1);
0412    unit_test_cmp('Check C2F min==0', min(c2f(:)), 0);
0413    
0414    f2c = bsxfun(@rdivide, bsxfun(@times, c2f, get_elem_volume(fmdl))', get_elem_volume(rmdl));
0415    unit_test_cmp('Check F2C max==1', max(sum(f2c,2)), 1, 1e-15);
0416    unit_test_cmp('Check F2C min==0', min(f2c(:)), 0);
0417 end