mk_tri2tet_c2f

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

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 6923 2024-05-31 10:07:13Z bgrychtol $
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    p=struct();
0081    p.debug = eidors_debug('query','mk_tri2tet_c2f:convhulln');
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          if p.debug
0327            p.fmdl = fmdl; p.rmdl = rmdl;
0328            p.v    = v;   p.t    = t;
0329            p.bot  = bot; p.top  = top; p.pts  = pts;
0330          end
0331          [~,V(last_v)] = convhulln_clean(pts,p);
0332       end
0333    end
0334    progress_msg(Inf);
0335    c2f = sparse(F,C,V,size(fmdl.elems,1),size(rmdl.elems,1));
0336    
0337    % add rtri contained in ftet
0338    if ~isinf(opt.z_depth)
0339       try rmdl = rmfield(rmdl,'coarse2fine'); end % messes with volume
0340       c2f = c2f + bsxfun(@times, sparse(tri_in_tet), opt.height * get_elem_volume(rmdl))';
0341    end
0342    % normalize to tet volume
0343    vol = get_elem_volume(fmdl);
0344    c2f = bsxfun(@rdivide,c2f,vol);
0345    
0346    % add tets contained in vox
0347    
0348    c2f = c2f + tet_in_tri;
0349    
0350 end
0351 
0352 function [intpts, tri2edge, tri2pts, edge2pts] = get_tet_intersection_pts(fmdl,rmdl,top,bot, epsilon)
0353    intpts = [];
0354    pt_idx = unique(rmdl.elems);
0355    pts = rmdl.nodes(pt_idx,1:2);
0356    
0357    bb_min = min(...
0358                 min(fmdl.nodes(fmdl.faces(:,1),1:2),...
0359                     fmdl.nodes(fmdl.faces(:,2),1:2)),...
0360                 fmdl.nodes(fmdl.faces(:,3),1:2));
0361    bb_max = max(...
0362                 max(fmdl.nodes(fmdl.faces(:,1),1:2),...
0363                     fmdl.nodes(fmdl.faces(:,2),1:2)),...
0364                 fmdl.nodes(fmdl.faces(:,3),1:2));
0365    todo = ~(  bsxfun(@gt,bb_min(:,1),pts(:,1)') ...
0366             | bsxfun(@gt,bb_min(:,2),pts(:,2)') ...
0367             | bsxfun(@lt,bb_max(:,1),pts(:,1)') ...
0368             | bsxfun(@lt,bb_max(:,2),pts(:,2)')); 
0369    [F,P] = find(todo);
0370    P = unique(P);
0371    in = false(size(fmdl.faces,1),size(pts,1));
0372    in(F,P) = point_in_triangle(pts(P,:),fmdl.faces(F,:),fmdl.nodes(:,1:2),epsilon)';
0373    
0374    [F,P] = find(in);
0375 
0376    % remove "vertical" faces
0377    vf = fmdl.normals(F,3) == 0;
0378    F(vf) = [];
0379    P(vf) = [];
0380    
0381    % project on faces
0382    % plane equation is ax+by+cz+d = 0, where d = -(ax0 + by0 + cz0)
0383    z = sum(fmdl.normals(F,:) .* fmdl.nodes(fmdl.faces(F,1),:),2);
0384 %    z = repmat(d,1,length(P));
0385    for j = 1:2
0386       z = z - fmdl.normals(F,j) .* pts(P,j);
0387    end
0388    z = z ./ fmdl.normals(F,3);
0389    out = z>top | z < bot;
0390    F(out) = [];
0391    P(out) = [];
0392    z(out) = [];
0393    
0394    intpts = [pts(P,:) z];
0395    I = (1:size(intpts,1))';
0396    F = cast(F,'like',pt_idx); % make sparse happy
0397    tri2edge = sparse(F,pt_idx(P),I,size(fmdl.faces,1),size(rmdl.nodes,1));
0398    I = cast(I,'like',pt_idx); % make sparse happy
0399    tri2pts = sparse(F,I,ones(size(I,1),1),size(fmdl.faces,1),size(intpts,1));
0400    edge2pts = sparse(pt_idx(P),I,ones(size(I,1),1),size(rmdl.nodes,1),size(intpts,1));
0401    
0402 end
0403 
0404 function [intpts, FE2CE, FE2pts, CE2pts] = ...
0405    find_edge2edge_intersections_2d(FE, FN, CE, CN, top, bot)
0406 
0407    P1 = FN(FE(:,1),:);
0408    P2 = FN(FE(:,2),:);
0409    P3 = CN(CE(:,1),:);
0410    P4 = CN(CE(:,2),:);
0411 
0412    C_bb = zeros(size(P3,1),4);
0413    C_bb(:,[1 3]) = min(P3(:,1:2),P4(:,1:2));
0414    C_bb(:,[2 4]) = max(P3(:,1:2),P4(:,1:2));
0415 
0416    F_bb = zeros(size(P1,1),4);
0417    F_bb(:,[1 3]) = min(P1(:,1:2),P2(:,1:2));
0418    F_bb(:,[2 4]) = max(P1(:,1:2),P2(:,1:2));
0419 
0420 
0421    todo =   bsxfun(@gt, F_bb(:,1), C_bb(:,2)') ...
0422           | bsxfun(@lt, F_bb(:,2), C_bb(:,1)') ...
0423           | bsxfun(@gt, F_bb(:,3), C_bb(:,4)') ...
0424           | bsxfun(@lt, F_bb(:,4), C_bb(:,3)');
0425    todo = ~todo;
0426 
0427    [T, V] = find(todo);
0428    
0429    S1 = P2(T,:) - P1(T,:);
0430    S2 = P4(V,:) - P3(V,:);
0431    
0432    denom =    S2(:,2) .* S1(:,1) - S2(:,1) .* S1(:,2);
0433 
0434    P13 = P1(T,1:2) - P3(V,1:2);
0435 
0436    num_a =    S2(:,1) .* P13(:,2) ...
0437             - S2(:,2) .* P13(:,1);
0438    num_b =    S1(:,1) .* P13(:,2) ...
0439             - S1(:,2) .* P13(:,1);
0440    
0441    mua = num_a ./ denom;
0442    mub = num_b ./ denom;
0443    
0444    IN = mua>0 & mua<1 & mub>0 & mub<1;
0445    T = T(IN);
0446    V = V(IN);
0447    intpts = P1(T,:) + bsxfun(@times, mua(IN), S1(IN,:));
0448    in = ~(intpts(:,3) > top | intpts(:,3) < bot);
0449    intpts = intpts(in,:);
0450    I = (1:size(intpts,1))';
0451    T = T(in);
0452    V = V(in);
0453    FE2CE = sparse(size(P1,1),size(P3,1));
0454    FE2CE(sub2ind(size(FE2CE),T,V)) = I;
0455    FE2pts = sparse(T,I,ones(size(I)),size(P1,1),size(I,1));
0456    CE2pts  = sparse(V,I,ones(size(I)),size(P3,1),size(I,1));
0457    
0458 end
0459 
0460 function [top_node2tet, top_nodes, nodes_above, nodes_below] = ...
0461    get_cap_nodes_in_tets(fmdl,rmdl,top,epsilon)
0462 
0463    top_nodes = rmdl.nodes;
0464    top_nodes(:,3) = top;
0465    nodes_above = fmdl.nodes(:,3) >= top;
0466    nodes_below = fmdl.nodes(:,3) <= top;
0467    
0468    % nodes in tets
0469    elidx = any(nodes_above(fmdl.elems),2) & any(nodes_below(fmdl.elems),2);
0470    top_node2tet = sparse(size(rmdl.nodes,1),size(fmdl.elems,1));
0471    if any(elidx)
0472       mdl = struct;
0473       mdl.nodes = fmdl.nodes;
0474       mdl.elems = fmdl.elems(elidx,:);
0475       top_node2tet(:,elidx) = point_in_tet(mdl,top_nodes,epsilon);
0476    end
0477 end
0478 
0479 %-------------------------------------------------------------------------%
0480 % Intersections between tet edges and tri faces
0481 function [intpts, edge2tri, edge2pts, tri2pts] = ...
0482          get_cap_tet_edge2tri_face_intersections(fmdl,rmdl,top,...
0483                                           nodes_above,nodes_below, epsilon)
0484    
0485    intpts = [];
0486    edge2tri = sparse(size(fmdl.edges,1),size(rmdl.elems,1));
0487    edge2pts = sparse(size(fmdl.edges,1),0);
0488    tri2pts  = sparse(size(rmdl.elems,1),0);
0489    
0490 
0491    % tet_edge2tri_face
0492    edgeidx = any(nodes_above(fmdl.edges),2) & any(nodes_below(fmdl.edges),2);
0493    % discard nodes on the plane
0494    edgeidx(edgeidx) = fmdl.nodes(fmdl.edges(edgeidx,1),3) ~= ...
0495                       fmdl.nodes(fmdl.edges(edgeidx,2),3);
0496    
0497    v = fmdl.nodes(fmdl.edges(edgeidx,2),:) - ...
0498        fmdl.nodes(fmdl.edges(edgeidx,1),:);
0499    u = (top - fmdl.nodes(fmdl.edges(edgeidx,1),3))  ./ v(:,3);
0500    pts = fmdl.nodes(fmdl.edges(edgeidx,1),:) + bsxfun(@times,u, v);
0501    t = point_in_triangle(pts(:,1:2),rmdl.elems,rmdl.nodes(:,1:2),-epsilon);
0502    % remove any edges that cross multiple faces -- they will be cought by
0503    % edge2edge intersections elsewhere
0504    t(sum(t,2)>1,:) = false;
0505    [E,T] = find(t);
0506    
0507    if any(t(:))
0508       intpts = pts(E,:);
0509       N = size(intpts,1);
0510       I = (1:N)';
0511       emap = find(edgeidx);
0512       E = emap(E);
0513       edge2tri = sparse(E,T,I,size(fmdl.edges,1), size(rmdl.elems,1));
0514       edge2pts = sparse(E,I,ones(size(I)),size(fmdl.edges,1), N);
0515       tri2pts  = sparse(T,I,ones(size(I)),size(rmdl.elems,1), N);
0516    end
0517    
0518 end
0519 %-------------------------------------------------------------------------%
0520 % Intersections between tet faces and tri edges
0521 function   [intpts, tri2edge,tri2intpt,edge2intpt] = ...
0522                            get_cap_tet_face2tri_edge_intersections( ...
0523                             fmdl,rmdl,top,nodes_above,nodes_below, epsilon)
0524    
0525    USE_POINT_IN_TRIANGLE = 0;
0526    
0527    intpts = [];
0528    tri2edge = sparse(size(fmdl.faces,1),size(rmdl.edges,1));
0529    tri2intpt = sparse(size(fmdl.faces,1),0);
0530    edge2intpt  = sparse(size(rmdl.edges,1),0);
0531 
0532                               
0533    %faces that cross the cap
0534    faceidx = any(nodes_above(fmdl.faces),2) & any(nodes_below(fmdl.faces),2);
0535  
0536    faces = fmdl.faces(faceidx,:);
0537    normals = fmdl.normals(faceidx,:);
0538    
0539    N_edges = size(rmdl.edges,1);
0540    N_faces = size(faces,1);
0541    
0542   
0543    face_bb = zeros(N_faces,6);
0544    face_bb(:,1) = min(reshape(fmdl.nodes(faces,1),N_faces,3),[],2);
0545    face_bb(:,2) = max(reshape(fmdl.nodes(faces,1),N_faces,3),[],2);
0546    face_bb(:,3) = min(reshape(fmdl.nodes(faces,2),N_faces,3),[],2);
0547    face_bb(:,4) = max(reshape(fmdl.nodes(faces,2),N_faces,3),[],2);
0548    
0549    edge_bb = zeros(N_edges,6);
0550    edge_bb(:,1) = min(reshape(rmdl.nodes(rmdl.edges,1),N_edges,2),[],2);
0551    edge_bb(:,2) = max(reshape(rmdl.nodes(rmdl.edges,1),N_edges,2),[],2);
0552    edge_bb(:,3) = min(reshape(rmdl.nodes(rmdl.edges,2),N_edges,2),[],2);
0553    edge_bb(:,4) = max(reshape(rmdl.nodes(rmdl.edges,2),N_edges,2),[],2);
0554    
0555    todo =   bsxfun(@gt, face_bb(:,1), edge_bb(:,2)') ...
0556           | bsxfun(@lt, face_bb(:,2), edge_bb(:,1)') ...
0557           | bsxfun(@gt, face_bb(:,3), edge_bb(:,4)') ...
0558           | bsxfun(@lt, face_bb(:,4), edge_bb(:,3)');
0559    todo = ~todo;
0560    e_todo = any(todo,1);
0561    f_todo = any(todo,2);
0562    faceidx(faceidx) = f_todo;
0563    faces = faces(f_todo,:);
0564    normals = normals(f_todo,:);
0565    [F,E] = find(todo);
0566    emap = zeros(size(e_todo,2),1);
0567    emap(e_todo) = 1:nnz(e_todo);
0568    E = emap(E);
0569    fmap = zeros(size(f_todo,1),1);
0570    fmap(f_todo) = 1:nnz(f_todo);
0571    F = fmap(F);
0572    
0573    P1 = rmdl.nodes(rmdl.edges(e_todo,1),:);
0574    P12 = P1 - rmdl.nodes(rmdl.edges(e_todo,2),:);
0575    P1(:,3) = top;
0576    P12(:,3) = 0;
0577 
0578    d = sum(fmdl.normals(faceidx,:) .* fmdl.nodes(fmdl.faces(faceidx,1),:),2);
0579    
0580    if ~USE_POINT_IN_TRIANGLE
0581       % for point_in_triangle
0582       nodes1 = fmdl.nodes(faces(:,1),:);
0583       v0 = fmdl.nodes(faces(:,3),:) - nodes1;
0584       v1 = fmdl.nodes(faces(:,2),:) - nodes1;
0585       dot00 = dot(v0, v0, 2);
0586       dot01 = dot(v0, v1, 2);
0587       % dot02 = dot(v0, v2, 2);
0588       dot11 = dot(v1, v1, 2);
0589       % dot12 = dot(v1, v2, 2);
0590       invDenom = 1 ./ (dot00 .* dot11 - dot01 .* dot01);
0591    end
0592    
0593    % find points of intersection between edges and face planes
0594    num = -d(F) + sum(normals(F,:).*P1(E,:),2);
0595    den = sum(normals(F,:) .* P12(E,:),2);
0596    u = num ./ den;
0597    % den == 0 => normal perpendicular to line
0598    idx = u >= 0 & u <= 1 & abs(den) >= eps;
0599    
0600    if any(idx)      
0601       F = F(idx);
0602       E = E(idx);
0603       ipts = P1(E,:) - bsxfun(@times, u(idx), P12(E,:));
0604       if USE_POINT_IN_TRIANGLE
0605          t = point_in_triangle(ipts,faces(F,:),fmdl.nodes,epsilon,'match');
0606       else
0607          v2 = ipts - fmdl.nodes(faces(F,1),:);
0608          dot02 = dot(v0(F,:),v2,2);
0609          dot12 = dot(v1(F,:),v2,2);
0610          % barycentric coordinates
0611          dot01 = dot01(F);
0612          invDenom = invDenom(F);
0613          u = (dot11(F) .* dot02 - dot01 .* dot12) .* invDenom;
0614          v = (dot00(F) .* dot12 - dot01 .* dot02) .* invDenom;
0615          t = u >= -epsilon & v >= -epsilon & (u+v-epsilon) <= 1;
0616       end
0617       
0618       if any(t)
0619          N = nnz(t);
0620          idv = (1:N)';
0621          intpts = ipts(t,:);
0622          I = idv;
0623          idx = find(faceidx); idx = idx(F);
0624          F = idx(t);
0625          eimap = find(emap); 
0626          E = eimap(E(t));
0627          
0628          tri2edge = sparse(F,E,I,size(fmdl.faces,1),size(rmdl.edges,1));
0629          tri2intpt = sparse(F,I,ones(size(I)),size(fmdl.faces,1),size(I,1));
0630          edge2intpt  = sparse(E,I,ones(size(I)),size(rmdl.edges,1),size(I,1));
0631       end
0632    end
0633 end
0634 
0635 %-------------------------------------------------------------------------%
0636 % Center scale models
0637 function [fmdl,rmdl, opt] = center_scale_models(fmdl,rmdl, opt)
0638    ctr = mean([min(rmdl.nodes);max(rmdl.nodes)]);
0639    rmdl.nodes = bsxfun(@minus,rmdl.nodes,ctr);
0640    fmdl.nodes = bsxfun(@minus,fmdl.nodes,ctr);
0641    opt.top = opt.top - ctr(3);
0642    opt.bot = opt.bot - ctr(3);
0643    if isfield(opt,'do_not_scale') && opt.do_not_scale
0644       return
0645    end
0646    maxnode = min( max(abs(rmdl.nodes(:))), max(abs(fmdl.nodes(:))));
0647    scale = 1/maxnode;
0648    rmdl.nodes = scale*rmdl.nodes;
0649    fmdl.nodes = scale*fmdl.nodes;
0650    opt.top    = scale*opt.top;
0651    opt.bot    = scale*opt.bot;
0652    opt.height = scale*opt.height;
0653    opt.z_depth= scale*opt.z_depth;
0654    eidors_msg('@@ models scaled by %g', scale,2);
0655 end
0656 
0657 %-------------------------------------------------------------------------%
0658 % Remove obviously non-overlapping parts of the models
0659 function [fmdl,rmdl,fmdl_idx,rmdl_idx] = crop_models(fmdl,rmdl, opt)
0660    f_min = min(fmdl.nodes);
0661    f_max = max(fmdl.nodes);
0662    r_min = min(rmdl.nodes);
0663    r_max = max(rmdl.nodes);
0664    
0665    if isfield(opt, 'z_depth') && ~isinf(opt.z_depth)
0666       lvl = mean(rmdl.nodes(:,3));
0667       r_max(3) = lvl + opt.z_depth;
0668       r_min(3) = lvl - opt.z_depth;
0669    else
0670       r_max(3) = f_max(3);
0671       r_min(3) = f_min(3);
0672    end
0673    
0674    % nodes outside the bounding box of the other model
0675    f_gt  = bsxfun(@gt, fmdl.nodes, r_max);
0676    f_lt  = bsxfun(@lt, fmdl.nodes, r_min);
0677    r_gt  = bsxfun(@gt, rmdl.nodes, f_max);
0678    r_lt  = bsxfun(@lt, rmdl.nodes, f_min);
0679    
0680    % elems outside the bounding box of the other model
0681    re_gt = any(reshape(all(reshape(r_gt(rmdl.elems',:),3,[])),[],3),2);
0682    re_lt = any(reshape(all(reshape(r_lt(rmdl.elems',:),3,[])),[],3),2);
0683    fe_gt = any(reshape(all(reshape(f_gt(fmdl.elems',:),4,[])),[],3),2);
0684    fe_lt = any(reshape(all(reshape(f_lt(fmdl.elems',:),4,[])),[],3),2);
0685    
0686    % elems to keep
0687    rmdl_idx = ~(re_gt | re_lt);
0688    fmdl_idx = ~(fe_gt | fe_lt);
0689    
0690    % remove non-overlapping elems
0691    rmdl.elems = rmdl.elems(rmdl_idx,:);
0692    fmdl.elems = fmdl.elems(fmdl_idx,:);
0693    
0694    % remove unused nodes
0695    [r_used_nodes,jnk,r_n] = unique(rmdl.elems(:));
0696    [f_used_nodes,jnk,f_n] = unique(fmdl.elems(:));
0697    
0698    r_idx = 1:numel(r_used_nodes);
0699    f_idx = 1:numel(f_used_nodes);
0700    
0701    rmdl.elems = reshape(r_idx(r_n),[],3);
0702    fmdl.elems = reshape(f_idx(f_n),[],4);
0703    
0704    rmdl.nodes = rmdl.nodes(r_used_nodes,:);
0705    fmdl.nodes = fmdl.nodes(f_used_nodes,:);
0706    
0707    % for the benefit of any (debug) plots later on
0708    try, rmdl = rmfield(rmdl,'boundary'); end
0709    try, fmdl = rmfield(fmdl,'boundary'); end
0710 end
0711 
0712 %-------------------------------------------------------------------------%
0713 % Prepare model
0714 function fmdl = prepare_tet_mdl(fmdl)
0715    fmopt.elem2edge = true;
0716    fmopt.edge2elem = true;
0717    fmopt.face2elem = true;
0718    fmopt.node2elem = true;
0719    fmopt.normals   = true;
0720    ll = eidors_msg('log_level',1);
0721    fmopt.linear_reorder = false; % this is slow and not needed
0722    fmdl = fix_model(fmdl,fmopt);
0723    eidors_msg('log_level',ll);
0724    fmdl.node2elem = logical(fmdl.node2elem);
0725    nElem = size(fmdl.elems,1);
0726    nFace = size(fmdl.faces,1);
0727    fmdl.elem2face = sparse(repmat((1:nElem)',1,4),double(fmdl.elem2face),true,nElem,nFace);
0728 end
0729 
0730 %-------------------------------------------------------------------------%
0731 % Prepare model
0732 function fmdl = prepare_tri_mdl(fmdl)
0733    fmopt.elem2edge = true;
0734    fmopt.edge2elem = true;
0735 %    fmopt.face2elem = true;
0736    fmopt.node2elem = true;
0737 %    fmopt.normals   = true;
0738    fmopt.linear_reorder = false; % this is slow and not needed
0739    ll = eidors_msg('log_level',1);
0740    fmdl = fix_model(fmdl,fmopt);
0741    eidors_msg('log_level',ll);
0742    fmdl.node2elem = logical(fmdl.node2elem);
0743 
0744 end
0745 
0746 %-------------------------------------------------------------------------%
0747 % Debug plot
0748 function debug_plot(fmdl,rmdl,v,t, bot, top, pts)
0749    tet.nodes = fmdl.nodes;
0750    tri.nodes = repmat(rmdl.nodes(rmdl.elems(v,:),:),2,1);
0751    tri.nodes(:,3) = [repmat(bot,3,1); repmat(top,3,1)];
0752    tri.elems = [ 1 2 5 4
0753                  2 3 6 5
0754                  3 1 4 6];
0755    tri.boundary = tri.elems;
0756    tet.type = 'fwd_model';
0757    tri.type = 'fwd_model';
0758    tet.elems = fmdl.elems(t,:);
0759    clf
0760    show_fem(tri)
0761    hold on
0762    h = show_fem(tet);
0763    set(h,'EdgeColor','b')
0764 %    pts = bsxfun(@plus,pts*scale,ctr);
0765    plot3(pts(:,1),pts(:,2),pts(:,3),'o');
0766    hold off
0767    axis auto
0768 end
0769 
0770 %-------------------------------------------------------------------------%
0771 % Parse option struct
0772 function opt = parse_opts(fmdl,rmdl, opt)
0773 
0774    lvl = mean(rmdl.nodes(:,3));
0775    
0776    if ~isfield(opt, 'z_depth')
0777       opt.z_depth = Inf;
0778    end
0779    if isinf(opt.z_depth)
0780       opt.top = max(fmdl.nodes(:,3));
0781       opt.bot = min(fmdl.nodes(:,3));
0782       opt.height = opt.top - opt.bot;
0783    else
0784       opt.top = lvl + opt.z_depth;
0785       opt.bot = lvl - opt.z_depth;
0786       opt.height = 2 * opt.z_depth;
0787    end
0788    if ~isfield(opt, 'tol_node2tri');
0789       opt.tol_node2tri = eps; % * max(rmdl_rng,fmdl_rng)^3;
0790    end
0791    if ~isfield(opt, 'tol_node2tet');
0792       opt.tol_node2tet = eps; % * max(rmdl_rng,fmdl_rng)^3;
0793    end
0794    if ~isfield(opt, 'tol_edge2edge')
0795       opt.tol_edge2edge = 6*eps(min(max(abs(fmdl.nodes(:))),max(abs(rmdl.nodes(:)))));
0796    end
0797    if ~isfield(opt, 'tol_edge2tri')
0798       opt.tol_edge2tri = eps; %1e-10
0799    end
0800    %     if ~isfield(opt, 'save_memory')
0801    %        opt.save_memory = 0;
0802    %     end
0803    eidors_msg('@@ node2tet  tolerance = %g', opt.tol_node2tet,2);
0804    eidors_msg('@@ edge2edge tolerance = %g', opt.tol_edge2edge,2);
0805    eidors_msg('@@ edge2tri  tolerance = %g', opt.tol_edge2tri,2);
0806  end
0807  
0808  
0809 function do_unit_test
0810    do_small_test
0811    do_inf_test
0812    do_realistic_test
0813    do_centre_slice; % failing code from tutorial
0814 end
0815 
0816 function do_inf_test
0817    jnk = mk_common_model('a2c',16);
0818    tri = jnk.fwd_model;
0819    tri.nodes = 1.1*tri.nodes;
0820    jnk = mk_common_model('c3cr',16);
0821    tet = jnk.fwd_model;
0822 
0823    c2f_a = mk_tri2tet_c2f(tet,tri);
0824    c2f_n = mk_approx_c2f(tet,tri);
0825    
0826    prob = c2f_n ~= 0 & c2f_a == 0;
0827    unit_test_cmp('Assert no missing intersections', nnz(prob),0);
0828    unit_test_cmp('mk_tri2tet_c2f v approximate', c2f_n,c2f_a, .5);
0829 end
0830 
0831 
0832 function do_small_test
0833    jnk = mk_common_model('a2c',16);
0834    tri = jnk.fwd_model;
0835    tri.nodes = 1.1*tri.nodes;
0836    jnk = mk_common_model('c3cr',16);
0837    tet = jnk.fwd_model;
0838 %    tet.elems = tet.elems(8081,:);
0839 %    tri.elems = tri.elems(1,:);
0840    opt.z_depth = 0.5;
0841    c2f = mk_tri2tet_c2f(tet,tri,opt);
0842  
0843    clf
0844    subplot(131)
0845    show_fem(tet);
0846    hold on
0847    show_fem(tri);
0848    view(2)
0849    subplot(132)
0850    img = mk_image(tet,1);
0851    img.elem_data = sum(c2f,2);
0852    img.calc_colours.clim = 1;
0853    img.calc_colours.ref_level = 0;
0854    show_fem(img);
0855    subplot(133)
0856    img = mk_image(tri,1);
0857    img.elem_data = sum(bsxfun(@times,c2f,get_elem_volume(tet)),1) ./ get_elem_volume(tri)';
0858    img.calc_colours.clim = 1;
0859    img.calc_colours.ref_level = 0;
0860    show_fem(img)
0861    
0862    unit_test_cmp('Check C2F size  ', size(c2f),[length(tet.elems), length(tri.elems)]);
0863    unit_test_cmp('Check C2F max==1', max(c2f(:)), 1);
0864    unit_test_cmp('Check C2F min==0', min(c2f(:)), 0);
0865    
0866    f2c = bsxfun(@rdivide, bsxfun(@times, c2f, get_elem_volume(tet))', get_elem_volume(tri));
0867    unit_test_cmp('Check F2C max==1', max(sum(f2c,2)), 1, 1e-14);
0868    unit_test_cmp('Check F2C min==0', min(f2c(:)), 0);
0869 end
0870  
0871 
0872 function do_realistic_test
0873       
0874    fmdl= ng_mk_cyl_models([2,2,.1],[16,1],[.1,0,.025]);
0875    xvec = -2:.5:2; %[-1.5 -.5:.2:.5 1.5];
0876    yvec = -2:.5:2; %[-1.6 -1:.2:1 1.6];
0877    zvec = [.5 1.5];
0878    grid2d = mk_grid_model([],xvec,yvec);
0879    grid3d = mk_grid_model([],xvec,yvec,zvec);
0880    eidors_cache clear
0881    tic
0882    opt.save_memory = 0;
0883    c2f_a = mk_grid_c2f(fmdl, grid3d,opt);
0884    t = toc;
0885    fprintf('Voxel:\tt=%f s\n',t);
0886 
0887    opt.z_depth = .5;
0888    grid2d.nodes(:,3) = 1;
0889    eidors_cache clear
0890    tic
0891 %    fmdl.elems = fmdl.elems(2764,:);
0892 %    grid2d.elems = grid2d.elems(4,:);
0893    c2f_b = mk_tri2tet_c2f(fmdl, grid2d, opt);
0894    t = toc;
0895    fprintf('tri2tet:\tt=%f \n',t);
0896    c2f_c = c2f_b * grid2d.coarse2fine;
0897    
0898    eidors_cache clear
0899 
0900    grid2d.mk_coarse_fine_mapping.z_depth = .5;
0901    grid2d.mk_coarse_fine_mapping.f2c_offset(3) = 1;
0902    grid2d = rmfield(grid2d,'coarse2fine');
0903    tic
0904    c2f_n = mk_approx_c2f(fmdl,grid2d);
0905    t = toc;
0906    fprintf('Approximate: t=%f s\n',t);
0907    
0908    unit_test_cmp('mk_tri2tet_c2f v mk_grid_c2f', c2f_c,c2f_a, 1e-5);
0909    unit_test_cmp('mk_tri2tet_c2f v approximate', c2f_n,c2f_b, .5);
0910    prob = c2f_n ~= 0 & c2f_b == 0;
0911    unit_test_cmp('Assert no missing intersections', nnz(prob),0);
0912 %    keyboard
0913 %    subplot(132)
0914 %    imagesc(c2f_c, [0 1]);
0915 %    subplot(133)
0916 %    imagesc(c2f_a, [0 1]);
0917 %    figure
0918 %    subplot(131)
0919 % show_fem(mk_image(fmdl,sum(c2f_a,2)));
0920 % img_a = mk_image(fmdl,sum(c2f_a,2));
0921 % img_a.calc_colours.ref_level = 0;
0922 % show_fem(img_a)
0923 % subplot(132)
0924 % img_c = mk_image(fmdl,sum(c2f_c,2));
0925 % img_c.calc_colours.ref_level = 0;
0926 % show_fem(img_c)
0927 % subplot(133)
0928 % img_n = mk_image(fmdl,sum(c2f_n,2));
0929 % img_n.calc_colours.ref_level = 0;
0930 % show_fem(img_n)
0931 %    keyboard
0932 end
0933  
0934  
0935 function do_centre_slice; % failing code from tutorial
0936    imdl = mk_common_model('b3cr',[16,2]);
0937    f_mdl= imdl.fwd_model;
0938    imdl2d= mk_common_model('b2c2',16);
0939    c_mdl= imdl2d.fwd_model;
0940    c2f= mk_coarse_fine_mapping( f_mdl, c_mdl);
0941 end