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