mk_tri2tet_c2f

PURPOSE ^

MK_TRI2TET_C2F - coarse2fine mapping between tri-based and tet-based models

SYNOPSIS ^

function c2f = mk_tri2tet_c2f(fmdl,rmdl, opt)

DESCRIPTION ^

MK_TRI2TET_C2F - coarse2fine mapping between tri-based and tet-based models
 C2F = MK_TRI2TET_C2F(FMDL,RMDL) returns in C2F the fraction of volume of
 each element of the fine (tet-based) model contained in each element of
 the coarse (tri-based) model. Each triangle is extruded along the z-axis
 to form a prism.
 Uses CONVHULLN to calculate the volume defined by a set of intersection
 points between individual tet and vox elements.

 C2F = MK_TRI2TET_C2F(FMDL,RMDL,OPT) allows specifying options.
 
 Inputs:
   FMDL - an EIDORS (tet-based) 3D forward model
   RMDL - an EIDORS (tri-based) 2D forward model
   OPT  - an option structure with the following fields and defaults:
      .z_depth
      .do_not_scale  - set to true to prevent scaling the models to unit
                       cube before any calculations, including thresholds.
                       Default: false
      .tol_node2tri  - minimum value of a barycentric coordinate to 
                       decide a tet node is lying inside a triangle and 
                       not on its edge. Default: eps
      .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 an edge-plane intersection point is lying 
                       inside a triangle. 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.
 
 See also MK_GRID_C2F, MK_TET_C2F, MK_TRI_C2F, MK_COARSE_FINE_MAPPING,
          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_tri2tet_c2f(fmdl,rmdl, opt)
0002 %MK_TRI2TET_C2F - coarse2fine mapping between tri-based and tet-based models
0003 % C2F = MK_TRI2TET_C2F(FMDL,RMDL) returns in C2F the fraction of volume of
0004 % each element of the fine (tet-based) model contained in each element of
0005 % the coarse (tri-based) model. Each triangle is extruded along the z-axis
0006 % to form a prism.
0007 % Uses CONVHULLN to calculate the volume defined by a set of intersection
0008 % points between individual tet and vox elements.
0009 %
0010 % C2F = MK_TRI2TET_C2F(FMDL,RMDL,OPT) allows specifying options.
0011 %
0012 % Inputs:
0013 %   FMDL - an EIDORS (tet-based) 3D forward model
0014 %   RMDL - an EIDORS (tri-based) 2D forward model
0015 %   OPT  - an option structure with the following fields and defaults:
0016 %      .z_depth
0017 %      .do_not_scale  - set to true to prevent scaling the models to unit
0018 %                       cube before any calculations, including thresholds.
0019 %                       Default: false
0020 %      .tol_node2tri  - minimum value of a barycentric coordinate to
0021 %                       decide a tet node is lying inside a triangle and
0022 %                       not on its edge. Default: eps
0023 %      .tol_node2tet  - tolerance for determinant <= 0 in testing for
0024 %                       points inside tets. Default: eps
0025 %      .tol_edge2edge - maximum distance between "intersecting" edges
0026 %                       Default: 6*sqrt(3)*eps(a), where a is
0027 %                       min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))
0028 %      .tol_edge2tri  - minimum value of a barycentric coordinate to
0029 %                       decide an edge-plane intersection point is lying
0030 %                       inside a triangle. Default: eps
0031 %
0032 % NOTE that for grid-based models, such as returned by MK_GRID_MODEL or
0033 % MK_VOXEL_VOLUME, MK_GRID_C2F is much faster.
0034 %
0035 % Set eidors_msg 'log level' < 2 to supress output to command line.
0036 %
0037 % See also MK_GRID_C2F, MK_TET_C2F, MK_TRI_C2F, MK_COARSE_FINE_MAPPING,
0038 %          FIND_EDGE2EDGE_INTERSECTIONS, CONVHULLN, MK_APPROX_C2F,
0039 %          POINT_IN_TRIANGLE, EIDORS_MSG
0040 
0041 % (C) 2015 Bartlomiej Grychtol
0042 % License: GPL version 2 or 3
0043 % $Id: mk_tri2tet_c2f.m 5548 2017-06-18 14:52:36Z aadler $
0044 
0045 
0046 if nargin == 0 || (ischar(fmdl) && strcmp(fmdl,'UNIT_TEST')) 
0047    do_unit_test; 
0048    return; 
0049 end
0050 if nargin < 3
0051    opt = struct();
0052 end
0053 
0054 f_elems = size(fmdl.elems,1);
0055 r_elems = size(rmdl.elems,1);
0056 c2f = sparse(f_elems,r_elems);
0057 
0058 if size(rmdl.nodes,2) == 2
0059    rmdl.nodes(:,3) = max(fmdl.nodes(:,3))/2 + min(fmdl.nodes(:,3))/2;
0060 end
0061 
0062 [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl, opt);
0063 
0064 if ~(any(fmdl_idx) && any(rmdl_idx))
0065    eidors_msg('@@: models do not overlap, returning all-zeros');
0066    return
0067 end
0068 
0069 opt = parse_opts(fmdl,rmdl, opt);
0070 
0071 [fmdl,rmdl, opt] = center_scale_models(fmdl,rmdl, opt);
0072 
0073 copt.fstr = 'mk_tri2tet_c2f';
0074 c2f(fmdl_idx,rmdl_idx) = eidors_cache(@do_mk_tri2tet_c2f,{fmdl,rmdl,opt},copt);
0075 % c2f = eidors_cache(@do_mk_tri2tet_c2f,{fmdl,rmdl,opt},copt);
0076 
0077 end
0078 
0079 function c2f = do_mk_tri2tet_c2f(fmdl,rmdl,opt)
0080 
0081    DEBUG =  eidors_debug('query','mk_tri2tet_c2f');
0082 
0083    r_elems = size(rmdl.elems,1);
0084    f_elems = size(fmdl.elems,1);
0085  
0086    c2f = sparse(f_elems, r_elems);
0087    
0088    progress_msg('Prepare tet model...');
0089    fmdl = prepare_tet_mdl(fmdl);
0090    progress_msg(Inf);
0091    
0092    progress_msg('Prepare tri model...');
0093    rmdl = prepare_tri_mdl(rmdl);
0094    progress_msg(Inf);
0095 
0096    pmopt.final_msg = 'none';
0097    if ~isinf(opt.z_depth)
0098       % tri nodes on top and bot planes inside tets
0099       progress_msg('Find top cap nodes in tets...',-1,pmopt);
0100       [top_node2tet, top_nodes, top_nodes_above, top_nodes_below] = ...
0101          get_cap_nodes_in_tets(fmdl,rmdl,opt.top,opt.tol_node2tet);
0102       progress_msg(sprintf('Found %d', nnz(top_node2tet)), Inf);
0103       
0104       progress_msg('Find bottom cap nodes in tets...',-1,pmopt);
0105       [bot_node2tet, bot_nodes, bot_nodes_above, bot_nodes_below] = ...
0106          get_cap_nodes_in_tets(fmdl,rmdl,opt.bot,opt.tol_node2tet);
0107       progress_msg(sprintf('Found %d', nnz(bot_node2tet)), Inf);
0108       
0109       
0110       progress_msg('Find tet_face2tri_edge intersections (top) ...');
0111       [intpts4, top_tet_face2tri_edge,top_tet_face2intpt4,top_tri_edge2intpt4] = ...
0112          get_cap_tet_face2tri_edge_intersections( fmdl,rmdl,opt.top, ...
0113          top_nodes_above,top_nodes_below, opt.tol_edge2tri);
0114       progress_msg(sprintf('Found %d', size(intpts4,1)),Inf);
0115       
0116       progress_msg('Find tet_face2tri_edge intersections (bot) ...');
0117       [intpts5, bot_tet_face2tri_edge,bot_tet_face2intpt5,bot_tri_edge2intpt5] = ...
0118          get_cap_tet_face2tri_edge_intersections( fmdl,rmdl,opt.bot, ...
0119          bot_nodes_above,bot_nodes_below, opt.tol_edge2tri);
0120       progress_msg(sprintf('Found %d', size(intpts5,1)),Inf);
0121       
0122       % find tet_edge2tri_face intersections
0123       progress_msg('Find tet_edge2tri_face intersections (top) ...');
0124       [intpts6, top_tet_edge2tri_face, top_tet_edge2intpt6, tri2intpt6] = ...
0125          get_cap_tet_edge2tri_face_intersections(fmdl,rmdl,opt.top, ...
0126          top_nodes_above, top_nodes_below, opt.tol_edge2tri);
0127       progress_msg(sprintf('Found %d', size(intpts6,1)),Inf);
0128       
0129       
0130       progress_msg('Find tet_edge2tri_face intersections (bot) ...');
0131       [intpts7, bot_tet_edge2tri_face, bot_tet_edge2intpt7, tri2intpt7] = ...
0132          get_cap_tet_edge2tri_face_intersections(fmdl,rmdl,opt.bot, ...
0133          bot_nodes_above, bot_nodes_below, opt.tol_edge2tri);
0134       progress_msg(sprintf('Found %d', size(intpts6,1)),Inf);
0135       
0136       
0137       
0138       progress_msg('Find tet_edge2tri_edge intersections (top) ...',-1,pmopt);
0139       edgeidx = any(top_nodes_above(fmdl.edges),2) & any(top_nodes_below(fmdl.edges),2);
0140       nodes = rmdl.nodes;
0141       nodes(:,3) = opt.top;
0142       top_tet_edge2tri_edge = sparse(size(fmdl.edges,1),size(rmdl.edges,1));
0143       [intpts8, top_tet_edge2tri_edge(edgeidx,:), top_tet_edge2intpt8(edgeidx,:), top_tri_edge2intpt8] = ...
0144          find_edge2edge_intersections(fmdl.edges(edgeidx,:),fmdl.nodes, ...
0145          rmdl.edges,nodes, opt.tol_edge2edge);
0146       if size(top_tet_edge2intpt8,1)~=size(fmdl.edges,1)
0147          if ~isempty(intpts8)
0148             top_tet_edge2intpt8(size(fmdl.edges,1),1) = 0;
0149          else
0150             top_tet_edge2intpt8 = sparse(size(fmdl.edges,1),0);
0151          end
0152       end
0153       progress_msg(sprintf('Found %d', size(intpts8,1)),Inf);
0154       
0155       
0156       progress_msg('Find tet_edge2tri_edge intersections (bot) ...',-1,pmopt);
0157       edgeidx = any(bot_nodes_above(fmdl.edges),2) & any(bot_nodes_below(fmdl.edges),2);
0158       nodes = rmdl.nodes;
0159       nodes(:,3) = opt.bot;
0160       bot_tet_edge2tri_edge = sparse(size(fmdl.edges,1),size(rmdl.edges,1));
0161       [intpts9, bot_tet_edge2tri_edge(edgeidx,:), ...
0162          bot_tet_edge2intpt9(edgeidx,:), bot_tri_edge2intpt9] = ...
0163          find_edge2edge_intersections(fmdl.edges(edgeidx,:),fmdl.nodes, ...
0164          rmdl.edges,nodes, opt.tol_edge2edge);
0165       if size(bot_tet_edge2intpt9,1)~=size(fmdl.edges,1)
0166          if ~isempty(intpts9)
0167             bot_tet_edge2intpt9(size(fmdl.edges,1),1) = 0;
0168          else
0169             bot_tet_edge2intpt9 = sparse(size(fmdl.edges,1),0);
0170          end
0171       end
0172       progress_msg(sprintf('Found %d', size(intpts9,1)),Inf);
0173       
0174       % tri contained in tet
0175       progress_msg('Find tris contained in tet...')
0176       tri_in_tet = rmdl.node2elem' * (bot_node2tet + top_node2tet) == 6;
0177       progress_msg(sprintf('Found %d',nnz(tri_in_tet)), Inf);
0178    end
0179 
0180    % tet nodes inside triangles
0181    % positive tolerance to catch also points on the edges.
0182    progress_msg('Find tet_nodes in tri_elems...');
0183    in = (fmdl.nodes(:,3) >= opt.bot) & (fmdl.nodes(:,3) <= opt.top);
0184    tet_node2tri = spalloc(size(fmdl.nodes,1),size(rmdl.elems,1),nnz(in));
0185    tet_node2tri(in,:) = point_in_triangle(fmdl.nodes(in,1:2), ...
0186                                     rmdl.elems, ...
0187                                     rmdl.nodes(:,1:2), ...
0188                                     opt.tol_node2tri);      
0189    progress_msg(sprintf('Found %d', nnz(tet_node2tri)), Inf);
0190   
0191    n_nodes = size(rmdl.nodes,1);
0192    vert_edges(:,1) = 1:n_nodes;
0193    vert_edges(:,2) = vert_edges(:,1) + n_nodes;
0194    
0195    progress_msg('Find tet_edge2vert_edge intersections...',-1,pmopt);
0196    if isinf(opt.z_depth)
0197       bot_nodes = rmdl.nodes; bot_nodes(:,3) = opt.bot - 1;
0198       top_nodes = rmdl.nodes; top_nodes(:,3) = opt.top + 1;
0199    end
0200    [intpts1,tet_edge2vert_edge,tet_edge2intpt1,vert_edge2intpt1] = ...
0201       find_edge2edge_intersections( fmdl.edges, fmdl.nodes,...
0202                                     vert_edges, [bot_nodes; top_nodes], ...
0203                                     opt.tol_edge2edge);
0204    progress_msg(sprintf('Found %d',size(intpts1,1)),Inf);
0205    
0206    % find tet_edges that cross the ractangular faces between the top and
0207    % bottom caps
0208    progress_msg('Find tet_edge2vert_face intersections...');
0209    if isinf(opt.z_depth)
0210       top = opt.top + 1;
0211       bot = opt.bot - 1;
0212    else
0213       top = opt.top;
0214       bot = opt.bot;
0215    end
0216    [intpts2,tet_edge2tri_edge,tet_edge2intpt2,tri_edge2intpt2] = ...
0217       find_edge2edge_intersections_2d( fmdl.edges, fmdl.nodes(:,1:3), ...
0218                                        rmdl.edges, rmdl.nodes(:,1:2), ...
0219                                        top, bot);
0220    progress_msg(sprintf('Found %d',size(intpts2,1)),Inf);
0221 
0222    
0223    progress_msg('Find tet_face2vert_edge intersections...');
0224    [intpts3, tet_face2vert_edge, tet_face2intpt3, vert_edge2intpt3] = ...
0225       get_tet_intersection_pts(fmdl,rmdl,top,bot,opt.tol_edge2tri);
0226    progress_msg(sprintf('Found %d',size(intpts3,1)),Inf);
0227    
0228       
0229    % tet contained in tri
0230    progress_msg('Find tets contained in tri...');
0231    tet_in_tri = (double(fmdl.node2elem') * tet_node2tri) == 4;
0232    progress_msg(sprintf('Found %d',nnz(tet_in_tri)), Inf);
0233    
0234    % tets and tris that intersect
0235    progress_msg('Find total tri v. tet intersections...');
0236    tri2intTet = double(rmdl.edge2elem') * tet_edge2tri_edge' * fmdl.edge2elem ...
0237       |double(rmdl.node2elem') * (tet_face2vert_edge>0)' * fmdl.elem2face';
0238    if ~isinf(opt.z_depth)
0239       tri2intTet = tri2intTet ...
0240       | double(rmdl.edge2elem') * (top_tet_face2tri_edge + bot_tet_face2tri_edge)' * fmdl.elem2face' ...
0241       | (top_tet_edge2tri_face + bot_tet_edge2tri_face)' * fmdl.edge2elem;
0242    end
0243    % exclude complete inclusion
0244    tri2intTet = tri2intTet & ~tet_in_tri';
0245    if ~isinf(opt.z_depth)
0246       tri2intTet = tri2intTet & ~tri_in_tet;
0247    end
0248    progress_msg(sprintf('Found %d',nnz(tri2intTet)), Inf);
0249    
0250    
0251    progress_msg('Calculate intersection volumes...');
0252    % sparse logical multiplication doesn't exist
0253    tri2intpt1 = logical(rmdl.node2elem'*vert_edge2intpt1)';
0254    tet2intpt1 = logical(fmdl.edge2elem'* tet_edge2intpt1)';
0255    
0256    tet2intpt2 = logical(fmdl.edge2elem'*tet_edge2intpt2)';
0257    tri2intpt2 = logical(rmdl.edge2elem'*tri_edge2intpt2)';
0258    
0259    tet2intpt3 = logical(double(fmdl.elem2face)*tet_face2intpt3)';
0260    tri2intpt3 = logical(rmdl.node2elem'*vert_edge2intpt3)';
0261    if ~isinf(opt.z_depth)
0262       tet2intpt4  = logical(double(fmdl.elem2face)*top_tet_face2intpt4)';
0263       tri2intpt4  = logical(rmdl.edge2elem'*top_tri_edge2intpt4)';
0264       
0265       tet2intpt5  = logical(double(fmdl.elem2face)*bot_tet_face2intpt5)';
0266       tri2intpt5  = logical(rmdl.edge2elem'*bot_tri_edge2intpt5)';
0267       
0268       tet2intpt6  = logical(fmdl.edge2elem'*top_tet_edge2intpt6)';
0269       tri2intpt6  = tri2intpt6';
0270       
0271       tet2intpt7  = logical(fmdl.edge2elem'*bot_tet_edge2intpt7)';
0272       tri2intpt7  = tri2intpt7';
0273    
0274       tet2intpt8  = logical(fmdl.edge2elem'*top_tet_edge2intpt8)';
0275       tri2intpt8  = logical(rmdl.edge2elem'*top_tri_edge2intpt8)';
0276       
0277       tet2intpt9  = logical(fmdl.edge2elem'*bot_tet_edge2intpt9)';
0278       tri2intpt9  = logical(rmdl.edge2elem'*bot_tri_edge2intpt9)';
0279    end
0280    tri_todo = find(sum(tri2intTet,2)>0);
0281    C = []; F = []; V = [];
0282    
0283    id = 0; lvox = length(tri_todo);
0284    mint = ceil(lvox/100);
0285    
0286    for v = tri_todo'
0287       id = id+1;
0288       if mod(id,mint)==0, progress_msg(id/lvox); end
0289       tet_todo = find(tri2intTet(v,:));
0290       common_intpts1 = bsxfun(@and,tri2intpt1(:,v), tet2intpt1(:,tet_todo));
0291       common_intpts2 = bsxfun(@and,tri2intpt2(:,v), tet2intpt2(:,tet_todo));
0292       common_intpts3 = bsxfun(@and,tri2intpt3(:,v), tet2intpt3(:,tet_todo));
0293       if ~isinf(opt.z_depth)
0294          common_intpts4 = bsxfun(@and,tri2intpt4(:,v), tet2intpt4(:,tet_todo));
0295          common_intpts5 = bsxfun(@and,tri2intpt5(:,v), tet2intpt5(:,tet_todo));
0296          common_intpts6 = bsxfun(@and,tri2intpt6(:,v), tet2intpt6(:,tet_todo));
0297          common_intpts7 = bsxfun(@and,tri2intpt7(:,v), tet2intpt7(:,tet_todo));
0298          common_intpts8 = bsxfun(@and,tri2intpt8(:,v), tet2intpt8(:,tet_todo));
0299          common_intpts9 = bsxfun(@and,tri2intpt9(:,v), tet2intpt9(:,tet_todo));
0300          top_nodes_tet = bsxfun(@and,top_node2tet(:,tet_todo), rmdl.node2elem(:,v));
0301          bot_nodes_tet = bsxfun(@and,bot_node2tet(:,tet_todo), rmdl.node2elem(:,v));
0302       end
0303       tet_nodes     = bsxfun(@and,tet_node2tri(:,v), fmdl.node2elem(:,tet_todo));
0304 
0305       C = [C; v*ones(numel(tet_todo),1)];
0306       F = [F; tet_todo'];
0307       last_v = numel(V);
0308       V = [V; zeros(numel(tet_todo),1)]; % pre-allocate
0309       for t = 1:numel(tet_todo)
0310          pts = [  intpts1(common_intpts1(:,t),:);
0311                   intpts2(common_intpts2(:,t),:);
0312                   intpts3(common_intpts3(:,t),:);
0313                   fmdl.nodes(tet_nodes(:,t),:);];
0314          if ~isinf(opt.z_depth) 
0315             pts = [ pts;
0316                   intpts4(common_intpts4(:,t),:);
0317                   intpts5(common_intpts5(:,t),:);
0318                   intpts6(common_intpts6(:,t),:);
0319                   intpts7(common_intpts7(:,t),:);
0320                   intpts8(common_intpts8(:,t),:);
0321                   intpts9(common_intpts9(:,t),:);
0322                   top_nodes(top_nodes_tet(:,t),:);
0323                   bot_nodes(bot_nodes_tet(:,t),:)];
0324          end
0325          last_v = last_v + 1;
0326          
0327          if size(pts,1) < 4
0328 %             debug_plot(fmdl,rmdl,v,tet_todo(t), bot, top, pts)
0329 %             keyboard
0330             continue
0331          end
0332          if any(isnan(pts(:))), keyboard, end
0333 %          debug_plot(fmdl,rmdl,v,tet_todo(t), bot, top, pts);
0334          try
0335             % move points to origin (helps for small elements at
0336             % large coordinates
0337             ctr = mean(pts);
0338             if any(isnan(ctr)), keyboard,end
0339             pts = bsxfun(@minus,pts,ctr);
0340             scale = max(abs(pts(:)));
0341             if scale == 0 %happens when there's only one point
0342                continue
0343             end
0344             % scale largest coordinate to 1 (helps with precision)
0345             pts = pts ./ scale;
0346             % force thorough search for initinal simplex and
0347             % supress precision warnings
0348             [K, V(last_v)] = convhulln(pts,{'Qt Pp Qs'});
0349             V(last_v) = max(V(last_v),0); % numerical issues may produce tiny negative volume
0350             V(last_v) = V(last_v) * scale^3; % undo scaling
0351          catch err
0352             ok = false;
0353             switch err.identifier
0354                case {'MATLAB:qhullmx:DegenerateData', 'MATLAB:qhullmx:UndefinedError'}
0355                   pts = bsxfun(@plus, pts .* scale, ctr);
0356                   u = uniquetol(pts,1e-14,'ByRows',true,'DataScale', 1);
0357                   ok = ok | size(u,1) < 4;
0358                   if ~ok
0359                      % test for colinearity in the xy plane
0360                      u12 = uniquetol(pts(:,1:2),1e-14,'ByRows',true,'DataScale',1);
0361                      cp = bsxfun(@minus, u12(2:end,1:2), u12(1,1:2));
0362                      l = sqrt(sum(cp.^2,2));
0363                      cp = bsxfun(@rdivide, cp, l);
0364                      % counteract colinear vectors in different directions
0365                      cp = abs(cp); 
0366                      un = uniquetol(cp,1e-12,'ByRows',true,'DataScale',1);
0367                      ok = ok | size(un,1) == 1; % co-linear points
0368                   end
0369                   if ~ok
0370                      % test if all points lie on the top or bottom caps
0371                      ok = ok | all(abs(pts(:,3) - top) < eps);
0372                      ok = ok | all(abs(pts(:,3) - bot) < eps);
0373                   end
0374             end
0375             if ~ok
0376                if DEBUG || eidors_debug('query','mk_tri2tet_c2f:convhulln')
0377                   debug_plot(fmdl,rmdl,v,tet_todo(t), bot, top, pts)
0378                   keyboard
0379                else
0380                   fprintf('\n');
0381                   eidors_msg(['convhulln has thrown an error. ' ...
0382                      'Enable eidors_debug on mk_tri2tet_c2f and re-run to see a debug plot'],0);
0383                   rethrow(err);
0384                end
0385             end
0386          end
0387       end
0388    end
0389    progress_msg(Inf);
0390    c2f = sparse(F,C,V,size(fmdl.elems,1),size(rmdl.elems,1));
0391    
0392    % add rtri contained in ftet
0393    if ~isinf(opt.z_depth)
0394       try rmdl = rmfield(rmdl,'coarse2fine'); end % messes with volume
0395       c2f = c2f + bsxfun(@times, sparse(tri_in_tet), opt.height * get_elem_volume(rmdl))';
0396    end
0397    % normalize to tet volume
0398    vol = get_elem_volume(fmdl);
0399    c2f = bsxfun(@rdivide,c2f,vol);
0400    
0401    % add tets contained in vox
0402    
0403    c2f = c2f + tet_in_tri;
0404    
0405 end
0406 
0407 function [intpts, tri2edge, tri2pts, edge2pts] = get_tet_intersection_pts(fmdl,rmdl,top,bot, epsilon)
0408    intpts = [];
0409    pt_idx = unique(rmdl.elems);
0410    pts = rmdl.nodes(pt_idx,1:2);
0411    
0412    bb_min = min(...
0413                 min(fmdl.nodes(fmdl.faces(:,1),1:2),...
0414                     fmdl.nodes(fmdl.faces(:,2),1:2)),...
0415                 fmdl.nodes(fmdl.faces(:,3),1:2));
0416    bb_max = max(...
0417                 max(fmdl.nodes(fmdl.faces(:,1),1:2),...
0418                     fmdl.nodes(fmdl.faces(:,2),1:2)),...
0419                 fmdl.nodes(fmdl.faces(:,3),1:2));
0420    todo = ~(  bsxfun(@gt,bb_min(:,1),pts(:,1)') ...
0421             | bsxfun(@gt,bb_min(:,2),pts(:,2)') ...
0422             | bsxfun(@lt,bb_max(:,1),pts(:,1)') ...
0423             | bsxfun(@lt,bb_max(:,2),pts(:,2)')); 
0424    [F,P] = find(todo);
0425    P = unique(P);
0426    in = false(size(fmdl.faces,1),size(pts,1));
0427    in(F,P) = point_in_triangle(pts(P,:),fmdl.faces(F,:),fmdl.nodes(:,1:2),epsilon)';
0428    
0429    [F,P] = find(in);
0430 
0431    % remove "vertical" faces
0432    vf = fmdl.normals(F,3) == 0;
0433    F(vf) = [];
0434    P(vf) = [];
0435    
0436    % project on faces
0437    % plane equation is ax+by+cz+d = 0, where d = -(ax0 + by0 + cz0)
0438    z = sum(fmdl.normals(F,:) .* fmdl.nodes(fmdl.faces(F,1),:),2);
0439 %    z = repmat(d,1,length(P));
0440    for j = 1:2
0441       z = z - fmdl.normals(F,j) .* pts(P,j);
0442    end
0443    z = z ./ fmdl.normals(F,3);
0444    out = z>top | z < bot;
0445    F(out) = [];
0446    P(out) = [];
0447    z(out) = [];
0448    
0449    intpts = [pts(P,:) z];
0450    I = (1:size(intpts,1))';
0451    tri2edge = sparse(F,pt_idx(P),I,size(fmdl.faces,1),size(rmdl.nodes,1));
0452    tri2pts = sparse(F,I,ones(size(I,1),1),size(fmdl.faces,1),size(intpts,1));
0453    edge2pts = sparse(pt_idx(P),I,ones(size(I,1),1),size(rmdl.nodes,1),size(intpts,1));
0454    
0455 end
0456 
0457 function [intpts, FE2CE, FE2pts, CE2pts] = ...
0458    find_edge2edge_intersections_2d(FE, FN, CE, CN, top, bot)
0459 
0460    P1 = FN(FE(:,1),:);
0461    P2 = FN(FE(:,2),:);
0462    P3 = CN(CE(:,1),:);
0463    P4 = CN(CE(:,2),:);
0464 
0465    C_bb = zeros(size(P3,1),4);
0466    C_bb(:,[1 3]) = min(P3(:,1:2),P4(:,1:2));
0467    C_bb(:,[2 4]) = max(P3(:,1:2),P4(:,1:2));
0468 
0469    F_bb = zeros(size(P1,1),4);
0470    F_bb(:,[1 3]) = min(P1(:,1:2),P2(:,1:2));
0471    F_bb(:,[2 4]) = max(P1(:,1:2),P2(:,1:2));
0472 
0473 
0474    todo =   bsxfun(@gt, F_bb(:,1), C_bb(:,2)') ...
0475           | bsxfun(@lt, F_bb(:,2), C_bb(:,1)') ...
0476           | bsxfun(@gt, F_bb(:,3), C_bb(:,4)') ...
0477           | bsxfun(@lt, F_bb(:,4), C_bb(:,3)');
0478    todo = ~todo;
0479 
0480    [T, V] = find(todo);
0481    
0482    S1 = P2(T,:) - P1(T,:);
0483    S2 = P4(V,:) - P3(V,:);
0484    
0485    denom =    S2(:,2) .* S1(:,1) - S2(:,1) .* S1(:,2);
0486 
0487    P13 = P1(T,1:2) - P3(V,1:2);
0488 
0489    num_a =    S2(:,1) .* P13(:,2) ...
0490             - S2(:,2) .* P13(:,1);
0491    num_b =    S1(:,1) .* P13(:,2) ...
0492             - S1(:,2) .* P13(:,1);
0493    
0494    mua = num_a ./ denom;
0495    mub = num_b ./ denom;
0496    
0497    IN = mua>0 & mua<1 & mub>0 & mub<1;
0498    T = T(IN);
0499    V = V(IN);
0500    intpts = P1(T,:) + bsxfun(@times, mua(IN), S1(IN,:));
0501    in = ~(intpts(:,3) > top | intpts(:,3) < bot);
0502    intpts = intpts(in,:);
0503    I = (1:size(intpts,1))';
0504    T = T(in);
0505    V = V(in);
0506    FE2CE = sparse(size(P1,1),size(P3,1));
0507    FE2CE(sub2ind(size(FE2CE),T,V)) = I;
0508    FE2pts = sparse(T,I,ones(size(I)),size(P1,1),size(I,1));
0509    CE2pts  = sparse(V,I,ones(size(I)),size(P3,1),size(I,1));
0510    
0511 end
0512 
0513 function [top_node2tet, top_nodes, nodes_above, nodes_below] = ...
0514    get_cap_nodes_in_tets(fmdl,rmdl,top,epsilon)
0515 
0516    top_nodes = rmdl.nodes;
0517    top_nodes(:,3) = top;
0518    nodes_above = fmdl.nodes(:,3) >= top;
0519    nodes_below = fmdl.nodes(:,3) <= top;
0520    
0521    % nodes in tets
0522    elidx = any(nodes_above(fmdl.elems),2) & any(nodes_below(fmdl.elems),2);
0523    top_node2tet = sparse(size(rmdl.nodes,1),size(fmdl.elems,1));
0524    if any(elidx)
0525       mdl = struct;
0526       mdl.nodes = fmdl.nodes;
0527       mdl.elems = fmdl.elems(elidx,:);
0528       top_node2tet(:,elidx) = point_in_tet(mdl,top_nodes,epsilon);
0529    end
0530 end
0531 
0532 %-------------------------------------------------------------------------%
0533 % Intersections between tet edges and tri faces
0534 function [intpts, edge2tri, edge2pts, tri2pts] = ...
0535          get_cap_tet_edge2tri_face_intersections(fmdl,rmdl,top,...
0536                                           nodes_above,nodes_below, epsilon)
0537    
0538    intpts = [];
0539    edge2tri = sparse(size(fmdl.edges,1),size(rmdl.elems,1));
0540    edge2pts = sparse(size(fmdl.edges,1),0);
0541    tri2pts  = sparse(size(rmdl.elems,1),0);
0542    
0543 
0544    % tet_edge2tri_face
0545    edgeidx = any(nodes_above(fmdl.edges),2) & any(nodes_below(fmdl.edges),2);
0546    % discard nodes on the plane
0547    edgeidx(edgeidx) = fmdl.nodes(fmdl.edges(edgeidx,1),3) ~= ...
0548                       fmdl.nodes(fmdl.edges(edgeidx,2),3);
0549    
0550    v = fmdl.nodes(fmdl.edges(edgeidx,2),:) - ...
0551        fmdl.nodes(fmdl.edges(edgeidx,1),:);
0552    u = (top - fmdl.nodes(fmdl.edges(edgeidx,1),3))  ./ v(:,3);
0553    pts = fmdl.nodes(fmdl.edges(edgeidx,1),:) + bsxfun(@times,u, v);
0554    t = point_in_triangle(pts(:,1:2),rmdl.elems,rmdl.nodes(:,1:2),-epsilon);
0555    % remove any edges that cross multiple faces -- they will be cought by
0556    % edge2edge intersections elsewhere
0557    t(sum(t,2)>1,:) = false;
0558    [E,T] = find(t);
0559    
0560    if any(t(:))
0561       intpts = pts(E,:);
0562       N = size(intpts,1);
0563       I = (1:N)';
0564       emap = find(edgeidx);
0565       E = emap(E);
0566       edge2tri = sparse(E,T,I,size(fmdl.edges,1), size(rmdl.elems,1));
0567       edge2pts = sparse(E,I,ones(size(I)),size(fmdl.edges,1), N);
0568       tri2pts  = sparse(T,I,ones(size(I)),size(rmdl.elems,1), N);
0569    end
0570    
0571 end
0572 %-------------------------------------------------------------------------%
0573 % Intersections between tet faces and tri edges
0574 function   [intpts, tri2edge,tri2intpt,edge2intpt] = ...
0575                            get_cap_tet_face2tri_edge_intersections( ...
0576                             fmdl,rmdl,top,nodes_above,nodes_below, epsilon)
0577    
0578    USE_POINT_IN_TRIANGLE = 0;
0579    
0580    intpts = [];
0581    tri2edge = sparse(size(fmdl.faces,1),size(rmdl.edges,1));
0582    tri2intpt = sparse(size(fmdl.faces,1),0);
0583    edge2intpt  = sparse(size(rmdl.edges,1),0);
0584 
0585                               
0586    %faces that cross the cap
0587    faceidx = any(nodes_above(fmdl.faces),2) & any(nodes_below(fmdl.faces),2);
0588  
0589    faces = fmdl.faces(faceidx,:);
0590    normals = fmdl.normals(faceidx,:);
0591    
0592    N_edges = size(rmdl.edges,1);
0593    N_faces = size(faces,1);
0594    
0595   
0596    face_bb = zeros(N_faces,6);
0597    face_bb(:,1) = min(reshape(fmdl.nodes(faces,1),N_faces,3),[],2);
0598    face_bb(:,2) = max(reshape(fmdl.nodes(faces,1),N_faces,3),[],2);
0599    face_bb(:,3) = min(reshape(fmdl.nodes(faces,2),N_faces,3),[],2);
0600    face_bb(:,4) = max(reshape(fmdl.nodes(faces,2),N_faces,3),[],2);
0601    
0602    edge_bb = zeros(N_edges,6);
0603    edge_bb(:,1) = min(reshape(rmdl.nodes(rmdl.edges,1),N_edges,2),[],2);
0604    edge_bb(:,2) = max(reshape(rmdl.nodes(rmdl.edges,1),N_edges,2),[],2);
0605    edge_bb(:,3) = min(reshape(rmdl.nodes(rmdl.edges,2),N_edges,2),[],2);
0606    edge_bb(:,4) = max(reshape(rmdl.nodes(rmdl.edges,2),N_edges,2),[],2);
0607    
0608    todo =   bsxfun(@gt, face_bb(:,1), edge_bb(:,2)') ...
0609           | bsxfun(@lt, face_bb(:,2), edge_bb(:,1)') ...
0610           | bsxfun(@gt, face_bb(:,3), edge_bb(:,4)') ...
0611           | bsxfun(@lt, face_bb(:,4), edge_bb(:,3)');
0612    todo = ~todo;
0613    e_todo = any(todo,1);
0614    f_todo = any(todo,2);
0615    faceidx(faceidx) = f_todo;
0616    faces = faces(f_todo,:);
0617    normals = normals(f_todo,:);
0618    [F,E] = find(todo);
0619    emap = zeros(size(e_todo,2),1);
0620    emap(e_todo) = 1:nnz(e_todo);
0621    E = emap(E);
0622    fmap = zeros(size(f_todo,1),1);
0623    fmap(f_todo) = 1:nnz(f_todo);
0624    F = fmap(F);
0625    
0626    P1 = rmdl.nodes(rmdl.edges(e_todo,1),:);
0627    P12 = P1 - rmdl.nodes(rmdl.edges(e_todo,2),:);
0628    P1(:,3) = top;
0629    P12(:,3) = 0;
0630 
0631    d = sum(fmdl.normals(faceidx,:) .* fmdl.nodes(fmdl.faces(faceidx,1),:),2);
0632    
0633    if ~USE_POINT_IN_TRIANGLE
0634       % for point_in_triangle
0635       nodes1 = fmdl.nodes(faces(:,1),:);
0636       v0 = fmdl.nodes(faces(:,3),:) - nodes1;
0637       v1 = fmdl.nodes(faces(:,2),:) - nodes1;
0638       dot00 = dot(v0, v0, 2);
0639       dot01 = dot(v0, v1, 2);
0640       % dot02 = dot(v0, v2, 2);
0641       dot11 = dot(v1, v1, 2);
0642       % dot12 = dot(v1, v2, 2);
0643       invDenom = 1 ./ (dot00 .* dot11 - dot01 .* dot01);
0644    end
0645    
0646    % find points of intersection between edges and face planes
0647    num = -d(F) + sum(normals(F,:).*P1(E,:),2);
0648    den = sum(normals(F,:) .* P12(E,:),2);
0649    u = num ./ den;
0650    % den == 0 => normal perpendicular to line
0651    idx = u >= 0 & u <= 1 & abs(den) >= eps;
0652    
0653    if any(idx)      
0654       F = F(idx);
0655       E = E(idx);
0656       ipts = P1(E,:) - bsxfun(@times, u(idx), P12(E,:));
0657       if USE_POINT_IN_TRIANGLE
0658          t = point_in_triangle(ipts,faces(F,:),fmdl.nodes,epsilon,'match');
0659       else
0660          v2 = ipts - fmdl.nodes(faces(F,1),:);
0661          dot02 = dot(v0(F,:),v2,2);
0662          dot12 = dot(v1(F,:),v2,2);
0663          % barycentric coordinates
0664          dot01 = dot01(F);
0665          invDenom = invDenom(F);
0666          u = (dot11(F) .* dot02 - dot01 .* dot12) .* invDenom;
0667          v = (dot00(F) .* dot12 - dot01 .* dot02) .* invDenom;
0668          t = u >= -epsilon & v >= -epsilon & (u+v-epsilon) <= 1;
0669       end
0670       
0671       if any(t)
0672          N = nnz(t);
0673          idv = (1:N)';
0674          intpts = ipts(t,:);
0675          I = idv;
0676          idx = find(faceidx); idx = idx(F);
0677          F = idx(t);
0678          eimap = find(emap); 
0679          E = eimap(E(t));
0680          
0681          tri2edge = sparse(F,E,I,size(fmdl.faces,1),size(rmdl.edges,1));
0682          tri2intpt = sparse(F,I,ones(size(I)),size(fmdl.faces,1),size(I,1));
0683          edge2intpt  = sparse(E,I,ones(size(I)),size(rmdl.edges,1),size(I,1));
0684       end
0685    end
0686 end
0687 
0688 %-------------------------------------------------------------------------%
0689 % Center scale models
0690 function [fmdl,rmdl, opt] = center_scale_models(fmdl,rmdl, opt)
0691    ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0692    rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0693    fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0694    opt.top = opt.top - ctr(3);
0695    opt.bot = opt.bot - ctr(3);
0696    if isfield(opt,'do_not_scale') && opt.do_not_scale
0697       return
0698    end
0699    maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0700    scale = 1/maxnode;
0701    rmdl.nodes = scale*rmdl.nodes;
0702    fmdl.nodes = scale*fmdl.nodes;
0703    opt.top    = scale*opt.top;
0704    opt.bot    = scale*opt.bot;
0705    opt.height = scale*opt.height;
0706    opt.z_depth= scale*opt.z_depth;
0707    eidors_msg('@@ models scaled by %g', scale,2);
0708 end
0709 
0710 %-------------------------------------------------------------------------%
0711 % Remove obviously non-overlapping parts of the models
0712 function [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl, opt)
0713    f_min = min(fmdl.nodes);
0714    f_max = max(fmdl.nodes);
0715    r_min = min(rmdl.nodes);
0716    r_max = max(rmdl.nodes);
0717    
0718    if isfield(opt, 'z_depth') && ~isinf(opt.z_depth)
0719       lvl = mean(rmdl.nodes(:,3));
0720       r_max(3) = lvl + opt.z_depth;
0721       r_min(3) = lvl - opt.z_depth;
0722    else
0723       r_max(3) = f_max(3);
0724       r_min(3) = f_min(3);
0725    end
0726    
0727    % nodes outside the bounding box of the other model
0728    f_gt  = bsxfun(@gt, fmdl.nodes, r_max);
0729    f_lt  = bsxfun(@lt, fmdl.nodes, r_min);
0730    r_gt  = bsxfun(@gt, rmdl.nodes, f_max);
0731    r_lt  = bsxfun(@lt, rmdl.nodes, f_min);
0732    
0733    % elems outside the bounding box of the other model
0734    re_gt = any(reshape(all(reshape(r_gt(rmdl.elems',:),3,[])),[],3),2);
0735    re_lt = any(reshape(all(reshape(r_lt(rmdl.elems',:),3,[])),[],3),2);
0736    fe_gt = any(reshape(all(reshape(f_gt(fmdl.elems',:),4,[])),[],3),2);
0737    fe_lt = any(reshape(all(reshape(f_lt(fmdl.elems',:),4,[])),[],3),2);
0738    
0739    % elems to keep
0740    rmdl_idx = ~(re_gt | re_lt);
0741    fmdl_idx = ~(fe_gt | fe_lt);
0742    
0743    % remove non-overlapping elems
0744    rmdl.elems = rmdl.elems(rmdl_idx,:);
0745    fmdl.elems = fmdl.elems(fmdl_idx,:);
0746    
0747    % remove unused nodes
0748    [r_used_nodes,jnk,r_n] = unique(rmdl.elems(:));
0749    [f_used_nodes,jnk,f_n] = unique(fmdl.elems(:));
0750    
0751    r_idx = 1:numel(r_used_nodes);
0752    f_idx = 1:numel(f_used_nodes);
0753    
0754    rmdl.elems = reshape(r_idx(r_n),[],3);
0755    fmdl.elems = reshape(f_idx(f_n),[],4);
0756    
0757    rmdl.nodes = rmdl.nodes(r_used_nodes,:);
0758    fmdl.nodes = fmdl.nodes(f_used_nodes,:);
0759    
0760    % for the benefit of any (debug) plots later on
0761    try, rmdl = rmfield(rmdl,'boundary'); end
0762    try, fmdl = rmfield(fmdl,'boundary'); end
0763 end
0764 
0765 %-------------------------------------------------------------------------%
0766 % Prepare model
0767 function fmdl = prepare_tet_mdl(fmdl)
0768    fmopt.elem2edge = true;
0769    fmopt.edge2elem = true;
0770    fmopt.face2elem = true;
0771    fmopt.node2elem = true;
0772    fmopt.normals   = true;
0773    ll = eidors_msg('log_level',1);
0774    fmopt.linear_reorder = false; % this is slow and not needed
0775    fmdl = fix_model(fmdl,fmopt);
0776    eidors_msg('log_level',ll);
0777    fmdl.node2elem = logical(fmdl.node2elem);
0778    nElem = size(fmdl.elems,1);
0779    nFace = size(fmdl.faces,1);
0780    fmdl.elem2face = sparse(repmat((1:nElem)',1,4),double(fmdl.elem2face),true,nElem,nFace);
0781 end
0782 
0783 %-------------------------------------------------------------------------%
0784 % Prepare model
0785 function fmdl = prepare_tri_mdl(fmdl)
0786    fmopt.elem2edge = true;
0787    fmopt.edge2elem = true;
0788 %    fmopt.face2elem = true;
0789    fmopt.node2elem = true;
0790 %    fmopt.normals   = true;
0791    fmopt.linear_reorder = false; % this is slow and not needed
0792    ll = eidors_msg('log_level',1);
0793    fmdl = fix_model(fmdl,fmopt);
0794    eidors_msg('log_level',ll);
0795    fmdl.node2elem = logical(fmdl.node2elem);
0796 
0797 end
0798 
0799 %-------------------------------------------------------------------------%
0800 % Debug plot
0801 function debug_plot(fmdl,rmdl,v,t, bot, top, pts)
0802    tet.nodes = fmdl.nodes;
0803    tri.nodes = repmat(rmdl.nodes(rmdl.elems(v,:),:),2,1);
0804    tri.nodes(:,3) = [repmat(bot,3,1); repmat(top,3,1)];
0805    tri.elems = [ 1 2 5 4
0806                  2 3 6 5
0807                  3 1 4 6];
0808    tri.boundary = tri.elems;
0809    tet.type = 'fwd_model';
0810    tri.type = 'fwd_model';
0811    tet.elems = fmdl.elems(t,:);
0812    clf
0813    show_fem(tri)
0814    hold on
0815    h = show_fem(tet);
0816    set(h,'EdgeColor','b')
0817 %    pts = bsxfun(@plus,pts*scale,ctr);
0818    plot3(pts(:,1),pts(:,2),pts(:,3),'o');
0819    hold off
0820    axis auto
0821 end
0822 
0823 %-------------------------------------------------------------------------%
0824 % Parse option struct
0825 function opt = parse_opts(fmdl,rmdl, opt)
0826 
0827    lvl = mean(rmdl.nodes(:,3));
0828    
0829    if ~isfield(opt, 'z_depth')
0830       opt.z_depth = Inf;
0831    end
0832    if isinf(opt.z_depth)
0833       opt.top = max(fmdl.nodes(:,3));
0834       opt.bot = min(fmdl.nodes(:,3));
0835       opt.height = opt.top - opt.bot;
0836    else
0837       opt.top = lvl + opt.z_depth;
0838       opt.bot = lvl - opt.z_depth;
0839       opt.height = 2 * opt.z_depth;
0840    end
0841    if ~isfield(opt, 'tol_node2tri');
0842       opt.tol_node2tri = eps; % * max(rmdl_rng,fmdl_rng)^3;
0843    end
0844    if ~isfield(opt, 'tol_node2tet');
0845       opt.tol_node2tet = eps; % * max(rmdl_rng,fmdl_rng)^3;
0846    end
0847    if ~isfield(opt, 'tol_edge2edge')
0848       opt.tol_edge2edge = 6*eps(min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))));
0849    end
0850    if ~isfield(opt, 'tol_edge2tri')
0851       opt.tol_edge2tri = eps; %1e-10
0852    end
0853    %     if ~isfield(opt, 'save_memory')
0854    %        opt.save_memory = 0;
0855    %     end
0856    eidors_msg('@@ node2tet  tolerance = %g', opt.tol_node2tet,2);
0857    eidors_msg('@@ edge2edge tolerance = %g', opt.tol_edge2edge,2);
0858    eidors_msg('@@ edge2tri  tolerance = %g', opt.tol_edge2tri,2);
0859  end
0860  
0861  
0862 function do_unit_test
0863    do_small_test
0864    do_inf_test
0865    do_realistic_test
0866    do_centre_slice; % failing code from tutorial
0867 end
0868 
0869 function do_inf_test
0870    jnk = mk_common_model('a2c',16);
0871    tri = jnk.fwd_model;
0872    tri.nodes = 1.1*tri.nodes;
0873    jnk = mk_common_model('c3cr',16);
0874    tet = jnk.fwd_model;
0875 
0876    c2f_a = mk_tri2tet_c2f(tet,tri);
0877    c2f_n = mk_approx_c2f(tet,tri);
0878    
0879    prob = c2f_n ~= 0 & c2f_a == 0;
0880    unit_test_cmp('Assert no missing intersections', nnz(prob),0);
0881    unit_test_cmp('mk_tri2tet_c2f v approximate', c2f_n,c2f_a, .5);
0882 end
0883 
0884 
0885 function do_small_test
0886    jnk = mk_common_model('a2c',16);
0887    tri = jnk.fwd_model;
0888    tri.nodes = 1.1*tri.nodes;
0889    jnk = mk_common_model('c3cr',16);
0890    tet = jnk.fwd_model;
0891 %    tet.elems = tet.elems(8081,:);
0892 %    tri.elems = tri.elems(1,:);
0893    opt.z_depth = 0.5;
0894    c2f = mk_tri2tet_c2f(tet,tri,opt);
0895  
0896    clf
0897    subplot(131)
0898    show_fem(tet);
0899    hold on
0900    show_fem(tri);
0901    view(2)
0902    subplot(132)
0903    img = mk_image(tet,1);
0904    img.elem_data = sum(c2f,2);
0905    img.calc_colours.clim = 1;
0906    img.calc_colours.ref_level = 0;
0907    show_fem(img);
0908    subplot(133)
0909    img = mk_image(tri,1);
0910    img.elem_data = sum(bsxfun(@times,c2f,get_elem_volume(tet)),1) ./ get_elem_volume(tri)';
0911    img.calc_colours.clim = 1;
0912    img.calc_colours.ref_level = 0;
0913    show_fem(img)
0914    
0915    unit_test_cmp('Check C2F size  ', size(c2f),[length(tet.elems), length(tri.elems)]);
0916    unit_test_cmp('Check C2F max==1', max(c2f(:)), 1);
0917    unit_test_cmp('Check C2F min==0', min(c2f(:)), 0);
0918    
0919    f2c = bsxfun(@rdivide, bsxfun(@times, c2f, get_elem_volume(tet))', get_elem_volume(tri));
0920    unit_test_cmp('Check F2C max==1', max(sum(f2c,2)), 1, 1e-14);
0921    unit_test_cmp('Check F2C min==0', min(f2c(:)), 0);
0922 end
0923  
0924 
0925 function do_realistic_test
0926       
0927    fmdl= ng_mk_cyl_models([2,2,.1],[16,1],[.1,0,.025]);
0928    xvec = -2:.5:2; %[-1.5 -.5:.2:.5 1.5];
0929    yvec = -2:.5:2; %[-1.6 -1:.2:1 1.6];
0930    zvec = [.5 1.5];
0931    grid2d = mk_grid_model([],xvec,yvec);
0932    grid3d = mk_grid_model([],xvec,yvec,zvec);
0933    eidors_cache clear
0934    tic
0935    opt.save_memory = 0;
0936    c2f_a = mk_grid_c2f(fmdl, grid3d,opt);
0937    t = toc;
0938    fprintf('Voxel:\tt=%f s\n',t);
0939 
0940    opt.z_depth = .5;
0941    grid2d.nodes(:,3) = 1;
0942    eidors_cache clear
0943    tic
0944 %    fmdl.elems = fmdl.elems(2764,:);
0945 %    grid2d.elems = grid2d.elems(4,:);
0946    c2f_b = mk_tri2tet_c2f(fmdl, grid2d, opt);
0947    t = toc;
0948    fprintf('tri2tet:\tt=%f \n',t);
0949    c2f_c = c2f_b * grid2d.coarse2fine;
0950    
0951    eidors_cache clear
0952 
0953    grid2d.mk_coarse_fine_mapping.z_depth = .5;
0954    grid2d.mk_coarse_fine_mapping.f2c_offset(3) = 1;
0955    grid2d = rmfield(grid2d,'coarse2fine');
0956    tic
0957    c2f_n = mk_approx_c2f(fmdl,grid2d);
0958    t = toc;
0959    fprintf('Approximate: t=%f s\n',t);
0960    
0961    unit_test_cmp('mk_tri2tet_c2f v mk_grid_c2f', c2f_c,c2f_a, 1e-5);
0962    unit_test_cmp('mk_tri2tet_c2f v approximate', c2f_n,c2f_b, .5);
0963    prob = c2f_n ~= 0 & c2f_b == 0;
0964    unit_test_cmp('Assert no missing intersections', nnz(prob),0);
0965 %    keyboard
0966 %    subplot(132)
0967 %    imagesc(c2f_c, [0 1]);
0968 %    subplot(133)
0969 %    imagesc(c2f_a, [0 1]);
0970 %    figure
0971 %    subplot(131)
0972 % show_fem(mk_image(fmdl,sum(c2f_a,2)));
0973 % img_a = mk_image(fmdl,sum(c2f_a,2));
0974 % img_a.calc_colours.ref_level = 0;
0975 % show_fem(img_a)
0976 % subplot(132)
0977 % img_c = mk_image(fmdl,sum(c2f_c,2));
0978 % img_c.calc_colours.ref_level = 0;
0979 % show_fem(img_c)
0980 % subplot(133)
0981 % img_n = mk_image(fmdl,sum(c2f_n,2));
0982 % img_n.calc_colours.ref_level = 0;
0983 % show_fem(img_n)
0984 %    keyboard
0985 end
0986  
0987  
0988 function do_centre_slice; % failing code from tutorial
0989    imdl = mk_common_model('b3cr',[16,2]);
0990    f_mdl= imdl.fwd_model;
0991    imdl2d= mk_common_model('b2c2',16);
0992    c_mdl= imdl2d.fwd_model;
0993    c2f= mk_coarse_fine_mapping( f_mdl, c_mdl);
0994 end

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